PQ分解法潮流计算原理与C++实现：从B‘矩阵到IEEE 33节点配电网验证

简介一份面向电力系统学习与开发者的基础资料包围绕PQ分解法潮流计算原理及C实现展开适合具备一定电力系统分析基础、希望结合代码理解潮流算法的读者。包内文件总共32个主要以txt输入输出数据、C源码与工程调试文件为构成另有exe可执行程序与doc程序说明便于对照验证。压缩包约519KB已有360人学习浏览。内容覆盖14、33、57、115、300节点等多个规模系统的计算配有原始数据与结果文件可直观对比PQ分解法在不同网络规模下的求解效果其中33节点潮流案例完整帮助读者掌握PQ节点与PV节点的处理方式、迭代收敛判断、C编程中节点和线路的数据结构设计以及求解流程。通过运行exe或研究源码读者可以复现计算结果并进一步扩展为教学演示或工程分析模块是理解潮流计算与数值方法结合的可操作参考。1. 33节点PQ分解法潮流程序核心不在“33”而在“分解”拿到一个叫 PQflow.rar 的工程文件名已经把要素写全33节点、PQ分解法、潮流、C。这种组合在课程设计和毕业设计里非常常见小系统用PQ分解法既避开牛顿法每次迭代重新形成雅可比矩阵的繁琐又能把稀疏矩阵和LU分解在C里完整练一遍。对5年以上经验的工程师这套程序的意义不在算法新鲜而在怎么把教材里的修正方程翻译成能编译、能收敛、能换数据文件的工程代码。本文按一线工程师的做法从原理推到代码再落到调参和验证。新手按步骤能跑通熟手可以重点看收敛判据和配电网的坑。2. PQ分解法原理从牛顿法修正方程里“剪”出来的两次近似2.1 极坐标修正方程与P-Q解耦先复习极坐标下牛顿-拉夫逊法的修正方程。对节点i有功和无功偏差写成ΔP H Δθ N ΔV / V ΔQ M Δθ L ΔV / VH、N、M、L是雅可比矩阵的四个分块。PQ分解法要做两次近似第一认为输电线路电阻R远小于电抗X有功主要受电压相角影响无功主要受电压幅值影响于是把N和M两块直接置零实现P-Q解耦第二在解耦后的方程里把H和L中的电压幅值按标幺值近似为1再把每轮迭代都变化的元素固定成只跟网络参数有关的常数矩阵得到B和B。两步之后修正方程变成B Δθ ΔP / V B ΔV ΔQ / VB和B在迭代开始前形成一次、做一次LU分解之后每一轮只做前代回代和功率更新这就是“量大又快”的由来。对比之下牛顿法每轮都要重新计算雅可比矩阵并重新分解单轮成本高出两个量级。2.2 B和B的形成规则不是简单取虚部B和B都来自节点导纳矩阵的虚部但取值细节不同。常见实现是先构建Y矩阵记YGjB再取B_imagimag(Y)。B取-B_imag同时忽略对地电容支路和非标准变比变压器的等效导纳然后删掉平衡节点所在的行和列。B同样取-b_imag但保留对地支路删掉平衡节点以及所有PV节点对应的行和列。对于纯PQ节点的33节点系统B和B的阶数相同都是32阶。但代码里要分开构造不要共用一个矩阵否则后续扩展成含PV节点的系统时会漏改。我一般会在构造时额外传一个布尔参数分别控制“是否包含对地支路”和“是否保留PV节点”这样两个矩阵共用同一段代码扩展成本很低。B和B的物理含义也值得说一句它们本质上是网络电纳矩阵决定了ΔP与Δθ、ΔQ与ΔV之间的线性化关系。迭代过程中这两个矩阵不更新所以收敛速度达不到牛顿法的二阶但胜在单轮计算量小、内存占用低。对这个规模LU分解做一次只要数毫秒。2.3 33节点配电网的R/X问题PQ分解还成立吗这里有个常被忽略的坑。课本推导PQ分解法时默认RX那是输电网的假设。而IEEE 33节点是经典配电网算例不少支路的R/X大于1比如标准参数里2-3节点之间R0.3660、X0.1864电阻接近电抗的两倍。直接套用PQ分解法严格说不满足推导前提。实际跑起来会怎样33节点规模小、PQ节点多、又没有PV节点用平启动初值PQ分解法通常仍能收敛只是迭代次数会比输电网算例多不少。真正危险的是重载工况当节点5之后的馈线负荷全部加满时电压可能低于0.95此时ΔP/V和ΔQ/V方程里“V近似等于1”的假设偏离加大可能出现收敛变慢甚至振荡。这时候不要急着改算法先试下面两个手段都比换牛顿法省事用BX法替代XB法即B用各支路X串联形成的网络B仍用Y矩阵虚部。配电网场景下BX法稳定性显著更好。给更新量加松弛因子比如Δθ只取原值的0.8倍。用牺牲迭代次数换稳定性在这个规模下完全值得。2.4 三种潮流算法选型对比方法单轮迭代成本典型迭代次数初值敏感度实现复杂度高斯-赛德尔低数百低低牛顿-拉夫逊高3~8高中PQ分解法低10~50低中表格说明一点PQ分解法和牛顿法在数学上是“同一族”解算器差别在于前者把雅可比矩阵“冻结”了。所以当牛顿法能收敛时PQ分解法通常也能收敛只是慢反过来PQ分解法对初值更宽容在配电网重载场景下反而更容易从平启动直接拉回收敛点。这就是标题这种“小系统配电网”组合选它的真正理由。3. C实现PQflow数据结构、Y矩阵与迭代主循环3.1 工程结构与节点/支路建模标题里的 PQflow.rar解压后常见结构是一个数据文件夹存放节点和支路参数一个主程序文件负责读数和迭代一个结果文件输出节点电压与网损。这种结构我一般会拆成三个文件bus.h、branch.h和main.cpp别嫌多余后面换算例文件的时候能省大量时间。节点和支路的数据结构如下struct Bus { int id; // 节点编号1基方便对照IEEE 33节点标准参数 int type; // 0PQ节点, 1PV节点, 2平衡节点 double P, Q; // 注入功率标幺值负荷为正 double V, theta; // 电压幅值和相角初值由平启动设定 }; struct Branch { int from, to; // 两端节点编号 double r, x; // 电阻、电抗标幺值 double g, b; // 对地导纳标幺值 };注意三点一是节点编号统一用1基与公开算例表对应避免把IEEE 33节点的“节点0”和C数组下标0混在一起二是有功和无功都按“电网流入为正”的约定读取负荷文件时要在前面取负三是Branch里保留g、b字段即使配电网没有对地支路也留着处理输电网算例时不用改结构。3.2 构建节点导纳矩阵Y33节点用稠密矩阵完全够用vectorvectorcomplexdouble即可没必要上稀疏矩阵。但代码要写得能直接切到CSR格式关键就是“先建支路、再组装Y”的两段式写法vectorvectorcomplexdouble buildY(const vectorBranch br, int n) { vectorvectorcomplexdouble Y(n 1, vectorcomplexdouble(n 1, complexdouble(0, 0))); for (const auto b : br) { complexdouble y(0, 0); if (b.r ! 0 || b.x ! 0) { y 1.0 / complexdouble(b.r, b.x); // 线路导纳 } Y[b.from][b.from] y complexdouble(b.g, b.b); Y[b.to][b.to] y complexdouble(b.g, b.b); Y[b.from][b.to] - y; Y[b.to][b.from] - y; } return Y; }这段代码做了三件小事忽略r和x同时为0的支路避免除零对角元累加线路导纳和对地导纳非对角元取负导纳。IEEE 33节点所有支路都是串联阻抗、无对地导纳因此g、b为0但这个函数直接拿去做IEEE 118节点输电网算例也能跑不用改代码。3.3 从Y矩阵提取B、B并做LU分解B和B构造是纯数字操作先提取Y的虚部取负再删行删列。下面用函数生成缩小后的矩阵并当场做LU分解vectorvectordouble buildB(const vectorvectorcomplexdouble Y, int n, bool keepShunt, int slack, const vectorint pvList) { vectorvectordouble B getNodalSusceptance(n, Y, keepShunt); removeRowsCols(B, slack); // 去掉平衡节点 for (int pv : pvList) removeRowsCols(B, pv); // 去掉PV节点 return B; } bool luDecompose(vectorvectordouble A) { int n A.size(); for (int k 0; k n - 1; k) { if (fabs(A[k][k]) 1e-12) { // 简易主元选取 int p k; while (p n fabs(A[p][k]) 1e-12) p; if (p n) return false; swap(A[k], A[p]); } for (int i k 1; i n; i) { A[i][k] / A[k][k]; for (int j k 1; j n; j) { A[i][j] - A[i][k] * A[k][j]; } } } return true; }getNodalSusceptance的实现在这里不多贴逻辑就是遍历Y的对角和非对角元取imag(Y[i][j])再加负号。removeRowsCols负责把指定编号的行列整体删除这里要注意编号映射删除后剩余节点序号代替不了原编号迭代时计算功率必须用原始编号索引Y矩阵。LU分解在程序启动阶段对B和B各执行一次即可之后不再进入迭代循环。这是PQ分解法高频面试题“为什么PQ分解比牛顿法快”的真正答案因子表固定。3.4 潮流主循环先算P更新θ再算Q更新V主循环按标准快速分解法的P-Q交替顺序写。先算不平衡量、更新相角再算无功、更新电压幅值。核心代码如下for (int iter 0; iter maxIter; iter) { calcPower(Y, V, theta, Pcalc, Qcalc, n); // 计算各节点注入功率 double dPmax 0; for (int i : pqList) { // 去掉平衡节点和PV节点的PQ节点 dP[i] (Pspec[i] - Pcalc[i]) / V[i]; dPmax max(dPmax, fabs(dP[i])); } auto dTheta luSolve(LU_Bp, pack(dP, idxMap)); for (size_t k 0; k dTheta.size(); k) { theta[originId[k]] dTheta[k]; // 把解写回完整相角数组 } calcPower(Y, V, theta, Pcalc, Qcalc, n); double dQmax 0; for (int i : pqList) { dQ[i] (Qspec[i] - Qcalc[i]) / V[i]; dQmax max(dQmax, fabs(dQ[i])); } auto dV luSolve(LU_Bpp, pack(dQ, idxMap)); for (size_t k 0; k dV.size(); k) { V[originId[k]] dV[k]; } if (dPmax tol dQmax tol) break; // 同时满足才认为收敛 }luSolve是前代回代函数代码很短但必须和luDecompose配对着看分解时原地修改了矩阵A回代时要按同一矩阵的上下三角结构计算。这里有一个容易踩的坑B和B如果阶数都是32那么dTheta、dV向量也是32维写回完整数组时靠originId做映射如果直接按下标顺序写回节点编号和数组下标错位结果会莫名其妙发散。3.5 初始化和迭代参数设置平启动是这类程序稳定收敛的关键所有PQ节点V1.0、theta0平衡节点保持V1.0、theta0。PS不要以为指定了负荷就能直接算初值质量决定PQ分解法前10轮迭代是稳步下降还是原地振荡。迭代参数经验值给一组maxIter设100次tol设1e-6标幺值。IEEE 33节点配电网从平启动出发一般20~40次收敛具体取决于负荷倍率。电压最小值出现在节点17附近约0.913这个值可以作为程序是否写对的“指纹”。4. 用IEEE 33节点算例验证参数、数据读取与结果对照4.1 标准算例参数与输入文件写法IEEE 33节点配电系统的公开参数很稳定基准电压12.66kV基准容量10MVA总负荷3715kW和2300kvar32条支路全部为串联阻抗。网上流传的参数表略有出入不同版本电阻值在小数点后第三位会不同不影响算法验证。关键是支路1、2、3的参数必须准确因为它们是馈线首段电压降落主要由这三条支路决定。支路号首端末端RΩXΩ1010.09220.04702120.49300.25113230.36600.18644340.38110.19415450.81900.70706560.18720.6188支路5的R/X为1.16支路6的R/X为0.30配电网就是这样电阻占比偏高。输入文件常见的格式是逗号分隔首行注释或直接跳过每行表示一条支路例from,to,r(ohm),x(ohm) 0,1,0.0922,0.0470 1,2,0.4930,0.2511读取时先用基准阻抗把欧姆值换成标幺值Zbase12.66²/1016.025Ω。不做这一步直接代入B矩阵会差好几个数量级收敛判断形同虚设。读取代码用ifstream加getline即可不推荐用scanf因为配电网参数文件常有多余空格或空行ifstream f(branch33.csv); string line; getline(f, line); // 跳过表头 while (getline(f, line)) { if (line.empty()) continue; replace(line.begin(), line.end(), ,, ); stringstream ss(line); Branch b; ss b.from b.to b.r b.x; b.r / Zbase; b.x / Zbase; branches.push_back(b); }注意从“节点编号0”到内部“编号1”的转换这里最简单有效的做法读取时直接b.from、b.to后续所有节点编号就都变成1基避免在迭代循环里反复判断是否为平衡节点。4.2 计算结果对照网损、最低电压与牛顿法对比程序输出的两个关键指标是全网总网损和最低节点电压。网损可以逐支路算double lossKW 0; for (const auto b : branches) { double I abs( (V[b.from] * exp(i * theta[b.from]) - V[b.to] * exp(i * theta[b.to])) / complexdouble(b.r, b.x) ); lossKW I * I * b.r * Sbase; // Sbase10MVA }IEEE 33节点标准算例下网损约202kW最低电压约0.913pu。下表给出一组典型的对比例子前提是同一套参数、同样的平启动初值解算方法网损kW最低电压pu迭代次数单轮耗时牛顿法202.680.91303~4高PQ分解法202.680.913015~40低高斯-赛德尔202.680.9130200低注意同一算例三种方法的最终解应当几乎一致差别只在迭代路径。如果你的PQ分解法结果和牛顿法差超过0.001pu不要怀疑算法理论先查B矩阵构造最常见的错误是把PV节点的行删错了或者删平衡节点时连带把电压幅值初始值也删了。还有一种隐蔽错误是Q计算时没有排除平衡节点导致平衡节点无功一直参与迭代结果不收敛也不报错。4.3 输出节点电压并画“电压走廊”算完保存结果通常写成node_id, V, theta_deg三列。33节点配电网最有辨识度的输出是电压曲线横轴节点编号纵轴电压幅值从0号节点的1.0一路下降到17号节点附近的0.913末段轻微回升。这条曲线是判断程序正确与否最直观的参考。写文件时建议保留至少6位有效数字并且把网损也附加在末尾。很多同学在验收时只贴一张控制台截图没有把节点电压输出成文件遇到结果异常无法复盘。我习惯把每轮迭代的最大不平衡量也追加到日志文件这样能清楚看到是P不收敛还是Q不收敛是线性下降还是振荡。5. 收敛慢时先改这四个地方再考虑换算法5.1 松弛因子最简单有效的“救火开关”PQ分解法更新量有时过冲尤其在重载工况下。给Δθ和ΔV乘一个介于0.6到0.9之间的系数能显著改善稳定性。实现只需一行theta[k] omega * dTheta[k]; V[k] omega * dV[k];代价是收敛速度下降但对33节点这种规模完全可接受。我在IEEE 33节点满载工况下用0.85的松弛因子迭代次数会从25次增加到35次左右但中途不再出现dPmax反弹到10⁻³以上的现象。如果系统负荷倍率到1.5倍不乘松弛因子直接发散的情况很常见。5.2 判据不要只盯ΔV很多初学者只在P-Q交替循环里判断最大电压偏差电压一稳定就提前跳出导致P不平衡量还挂着10⁻²量级的误差。正确做法是ΔPmax和ΔQmax同时小于阈值才算收敛电压变化只能作为参考不能作为唯一判据。原因是PQ分解法本质上是分别对P方程和Q方程迭代P方程收敛但Q方程不收敛的情况在配电网高无功负荷下确实会发生。我一般把tol拆成两个值有功tolP1e-6无功tolQ1e-6但在代码里同时判断且把两者的历史曲线都打印出来。如果dQmax一直降不到1e-6以下多半是B矩阵构造有误或无功初值偏离太多。5.3 平启动初值与“热启动”连续潮流连续潮流法的核心是逐步增加负荷系数λ每一步以上一步的电压解作为初值。PQ分解法在这里有个天然优势B和B不随负荷变化外循环把λ从0.2逐步升到1.5时内部每次潮流迭代都不用重新分解矩阵只改Pspec[i]乘λ再做前代回代。用这个思路可以从单一潮流程序平滑过渡到连续潮流法。5.4 验证程序的“电量”方法从单案例泛化拿到一个陌生算例文件先确认它的基准电压和基准功率是否和程序默认值一致再确认节点编号是否从0开始。最常见的翻车现场是支路表用欧姆负荷表用有名值程序里却按标幺值直接参赛。我的习惯是数据文件里显式写明基准值程序解析时一并读入并检查避免在代码里硬编码。最后分享一个验证技巧拿你最关心的末端节点电压33节点就是17号节点画一条“迭代次数-电压值”曲线。正确程序的表现是单调逼近0.913振荡递减或早早停在0.95以上都说明代码有逻辑错误。这种图形化验证比单纯看网损数值要灵敏得多因为在某些错误状态下网损也能算出一个看似合理的值。本文还有配套的精品资源点击获取