电力系统分析与仿真这个方向网上资料不少但大多要么堆理论公式、要么直接丢一个封装好的工具箱让读者跑黑箱。这次我把一个实际完成的“潮流计算不对称短路分析”Matlab项目完整拆开从数学模型怎么建立、雅可比矩阵怎么求到短路序网怎么连接、故障电流怎么算再到代码架构怎么设计、调试中有哪些坑一次讲清楚。适合正在做电力系统课程设计、毕业设计或者刚接触机电暂态分析、想从原理层面搞懂计算过程的同学参考。1. 为什么把潮流计算和不对称短路分析放在同一个工程里先聊点实际的在电力系统的规划、运行方式安排、保护整定这些日常工作里潮流计算和短路计算几乎是“搭班子”出现的两个计算工具。潮流计算解决的是“系统在某个稳态运行点下各节点电压、支路功率是多少”而短路计算解决的是“系统发生接地或者相间故障时故障电流有多大、各节点电压被拉低到什么程度”。前者决定了正常方式下的断面潮流和设备负载率后者决定了保护装置的动作电流整定值、断路器开断容量和设备动热稳定校验。既然两个计算都要基于同一套电网拓扑而且都要解网络方程那用一套代码把数据接口统一起来其实非常划算。我做这个项目时最直接的感受就是潮流计算需要节点导纳矩阵不对称短路的三序网络也要从原始导纳矩阵变换而来数据前处理几乎可以完全复用。短路计算需要一个“短路前”的运行状态——正好就是潮流计算的稳态解尤其是故障点处的电压幅值和相角会直接作为短路电流计算的初始条件。如果分开两个独立程序去做节点编号、基准值、支路参数很容易两边不一致最终对不上账。所以这个项目我直接用一套Matlab代码从读入节点和支路数据开始先做潮流再基于潮流结果做故障分析。整个程序分为三层数据层节点参数、支路参数、发电机参数、计算层导纳矩阵、潮流迭代、序网形成、故障解析、输出层电压分布表、支路潮流表、故障电流相量结果。这样模块之间的耦合度低调其中一个环节不会把别的部分搞坏。工程上常见的做法是用BPA、PSS/E这类商业软件直接出结果但它们的计算内核是封闭的不方便做二次开发和教学演示。用Matlab重写一遍麻烦归麻烦却能把每一步计算过程都暴露出来查问题、改算法、接入新的控制策略都容易得多。这个项目做完之后我甚至把里面的潮流求解函数单独拆出来接上了一个简单的连续潮流CPF扩展用来算PV曲线的鼻尖点这也算是模块化设计带来的额外红利。2. 潮流计算的数学本质我们到底在解什么方程潮流计算从数学上看就是求解一组非线性代数方程组。电力网络用节点导纳矩阵描述后每个节点的注入电流可以写成[ I_i \sum_{j1}^{n} Y_{ij} U_j ]其中 (Y_{ij}) 是导纳矩阵元素(U_j) 是节点电压相量。节点功率和电压、电流的关系为[ S_i U_i I_i^* P_i jQ_i ]把电流表达式代进去就得到了潮流方程的基本形式[ P_i jQ_i U_i \sum_{j1}^{n} (Y_{ij} U_j)^* ]展开成实部和虚部之后每个节点有两个功率平衡方程一个对应有功 (P)一个对应无功 (Q)。问题在于每个节点的状态变量是电压幅值 (U_i) 和相角 (\theta_i)也是两个未知数。系统中必须指定一个平衡节点slack bus承担功率差额它的电压幅值和相角固定其他节点的变量才有唯一解。节点类型惯例是这样的PQ节点给定有功和无功注入负荷节点通常归入此类PV节点给定有功注入和电压幅值发电机节点归入此类平衡节点给定电压幅值和相角选一个容量较大、接近系统中心的发电机节点计算量最集中的地方是求解修正方程。牛顿-拉夫逊法每一步需要计算不平衡量 (\Delta P)、(\Delta Q)以及雅可比矩阵然后求解线性方程组得到电压修正量。雅可比矩阵有四个分块[ J \begin{bmatrix} H N \ K L \end{bmatrix} ]其中 (H_{ij} \frac{\partial P_i}{\partial \theta_j})(N_{ij} \frac{\partial P_i}{\partial U_j})(K_{ij} \frac{\partial Q_i}{\partial \theta_j})(L_{ij} \frac{\partial Q_i}{\partial U_j})。具体表达式对非对角元素和对角元素是不同的。非对角元素(i \neq j)[ H_{ij} -U_i U_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij}) ][ N_{ij} -U_i U_j (G_{ij} \cos\theta_{ij} B_{ij} \sin\theta_{ij}) ][ K_{ij} U_i U_j (G_{ij} \cos\theta_{ij} B_{ij} \sin\theta_{ij}) ][ L_{ij} -U_i U_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij}) ]对角元素则要考虑自导纳带来的附加项。这些公式看似繁琐但规律很强只要把导纳矩阵和当前迭代的电压值代入就能批量生成雅可比矩阵。我用一个双层循环加if语句判断 (i) 和 (j) 是否相等来处理在Matlab里可读性很好如果要追求高性能工程实现用向量化或者稀疏矩阵构造会更快但教学和中等规模系统分析循环写法完全够用。潮流计算的收敛判据工程上一般取功率不平衡量最大绝对值小于一个阈值比如 (1 \times 10^{-6}) 或 (1 \times 10^{-8})标幺值。迭代初值大多数情况用“平启动”所有PQ节点电压幅值设1.0相角设0PV节点相角也设0平衡节点直接固定。对于绝大多数110kV及以上电网变压器分接头位置合适、无功补偿到位时平启动基本都能在5到8次迭代内收敛。3. Matlab代码实现数据准备、导纳矩阵与牛顿-拉夫逊迭代核心代码实现的第一步是数据准备。我习惯用结构体存储电网数据比较直观% 节点数据 % 节点编号 类型(1PQ, 2PV, 3平衡) 电压幅值初值 相角初值(度) Pg(MW) Qg(Mvar) Pd(MW) Qd(Mvar) busData [ 1 3 1.060 0 0 0 0 0 2 2 1.045 0 40 30 21.7 12.7 3 2 1.010 0 0 0 94.2 19.0 4 1 1.000 0 0 0 47.8 -3.9 % ... ]; % 支路数据 % 首端节点 末端节点 电阻R(pu) 电抗X(pu) 电纳B/2(pu) 变比k branchData [ 1 2 0.01938 0.05917 0.0264 1 1 5 0.05403 0.22304 0.0246 1 2 3 0.04699 0.19797 0.0219 1 % ... ];接下来构建节点导纳矩阵。这一步的物理含义很清晰自导纳是该节点连接的所有支路导纳之和互导纳是相连支路导纳的相反数。变压器支路要注意变比折算function Y buildYbus(branchData, busData) n size(busData, 1); Y zeros(n, n); nBranch size(branchData, 1); for k 1:nBranch i branchData(k, 1); j branchData(k, 2); r branchData(k, 3); x branchData(k, 4); b branchData(k, 5); tap branchData(k, 6); z r 1j * x; y 1 / z; yij y / tap; ysi y * (tap - 1) / tap 1j * b / 2; ysj y * (1 - tap) / (tap^2) 1j * b / 2; Y(i, i) Y(i, i) yij ysi; Y(j, j) Y(j, j) yij ysj; Y(i, j) Y(i, j) - yij; Y(j, i) Y(j, i) - yij; end end注意这里的 (tap) 我定义为非标准变比即从节点 i 侧看理想变压器变比为 (1:tap)。实际工程中变压器分接头不同公式略有差异建议在代码注释里标明约定免得日后自己都忘了。接下来是牛顿-拉夫逊迭代主体。为了不把所有代码都堆在一个脚本里我拆成三个函数calPower计算当前电压下的注入功率、calJacobian形成雅可比矩阵、updateVoltage解方程并更新电压相量。主迭代流程如下V busData(:, 3) .* exp(1j * deg2rad(busData(:, 4))); for iter 1:30 % 计算不平衡量 [Pcal, Qcal] calPower(V, Ybus); dP (busData(:, 5) - busData(:, 7)) - Pcal; % MW - pu需除以基准功率 dQ (busData(:, 6) - busData(:, 8)) - Qcal; % 剔除平衡节点和PV节点的无功不平衡量行 % 构造雅可比矩阵 J calJacobian(V, Ybus); % 求解修正方程 dTheta J \ dPQ; % 更新电压相角和幅值 theta theta dTheta(1:n-1); Vmag Vmag dTheta(n:end); if max(abs(dPQ)) 1e-8 break; end end需要提醒的是雅可比矩阵的行列要按系统中PQ节点和PV节点的排列来组织工程上常用两个数组记录节点类型索引比每次在循环里用if判断快得多。迭代过程中要限制PV节点的无功出力范围——如果某台发电机的无功越限应该把它降级为PQ节点处理否则潮流结果虽然收敛但物理上不可行。这个处理在程序里表现为动态修改节点类型重新排列雅可比矩阵的行列我在实现时单独写了一个子函数refreshNodeType每次迭代前检查一次。潮流结果输出时我喜欢把数据直接整理成表格形式打印出来方便写报告时直接引用。除了各节点电压幅值、相角还需要计算各支路潮流首端注入、末端注入、网损这个计算可以复用上面的功率计算公式只不过把注入功率换成某一条支路两端各自算一遍。4. 不对称短路分析的核心三序网络与复合序网潮流计算解决的是对称稳态问题而短路计算面临的是不对称情况这时不能直接在ABC三相的混合相量里求解因为三相阻抗耦合在一起电路方程不独立。工程上通用的做法是输出对称分量法把三相电气量分解成三组对称分量正序、负序、零序。对称分量的核心变换是[ \begin{bmatrix} \dot{F}_a \ \dot{F}_b \ \dot{F}_c \end{bmatrix}\begin{bmatrix} 1 1 1 \ 1 a^2 a \ 1 a a^2 \end{bmatrix} \begin{bmatrix} \dot{F}{a1} \ \dot{F}{a2} \ \dot{F}_{a0} \end{bmatrix} ]其中 (a e^{j120^\circ} -0.5 j0.866)。正序分量是三相对称、按正常相序排列的电压电流负序分量是按逆相序排列的零序分量则是三相幅值、相位都相同的分量。经过这个变换一个不对称的三相网络被拆成了三个各自对称的序网络互不耦合只在故障点通过故障边界条件连在一起。在实际的Matlab程序里只要生成一个变换矩阵T三序电压电流的计算就是矩阵乘法不需要手工按公式硬套。这一点简化了很多操作。比如求某一节点A相电压的序分量只需要T [1 1 1; 1 a^2 a; 1 a a^2]; % 注意这里的a exp(1j*2*pi/3) Fabc [Va; Vb; Vc]; F012 inv(T) * Fabc;不过短路分析更常用的思路不是从三相电压推序分量而是直接建立正序、负序、零序网络的节点导纳矩阵然后在故障节点处按短路类型把三个序网连起来。正序网络就是常规潮流用的导纳矩阵当然发电机、负荷模型在短路瞬间有特殊等值通常正序网要用次暂态电抗等值负荷可以近似忽略或接恒定阻抗。负序网络和正序网络拓扑相同只是发电机负序电抗和正序不同线路参数相同。零序网络则只包含能够为零序电流提供回路的路径比如变压器绕组接线方式YN接地与否、三角形绕组内部零序环流会直接决定零序网络的结构。复合序网法不需要手工推导各种短路情形的公式而是把故障点看成三序网络的连接口单相接地短路时三序网在故障点串联两相短路时正序网和负序网在故障点并联零序网不参与两相接地短路时三个序网在故障点并联。有了这个连接方式就可以把三个序网络方程和故障边界条件统一成一个复合线性方程组直接求解故障点的正、负、零序电流。工程上常用短路点电压形式的复合序网方程。比如单相接地短路A相的边界条件是[ \begin{cases} \dot{I}_b \dot{I}c 0 \ \dot{V}{ka} 0 \end{cases} ]等值到序分量后条件是[ \dot{I}{a1} \dot{I}{a2} \dot{I}{a0} \frac{\dot{V}{f}}{Z_1 Z_2 Z_0} ]其中 (\dot{V}_f) 是短路前故障点电压由潮流计算得到(Z_1, Z_2, Z_0) 是从故障点看进去的三序等值阻抗。这里要特别注意 (Z_1) 是正序网中故障点对地等值阻抗它和潮流计算得到的系统等值电势与故障点电压之间的阻抗有关实际求法是在故障点注入单位正序电流测该点电压两者相除就是正序等值阻抗。在代码实现上用阻抗矩阵法最省事。先对三序节点导纳矩阵分别求逆或者用LU分解解线性方程得到三序阻抗矩阵然后取故障节点对应的对角元素就是该点的三序等值阻抗。这个过程快速又不容易错Zseq{1} inv(Y1bus); % 正序阻抗矩阵 Zseq{2} inv(Y2bus); % 负序阻抗矩阵 Zseq{0} inv(Y0bus); % 零序阻抗矩阵 z1 Zseq{1}(k, k); z2 Zseq{2}(k, k); z0 Zseq{0}(k, k);得到三序等值阻抗后短路电流的计算就变成了简单的复数运算。拿到故障点序电流之后再用反变换得到三相短路电流Iabc T * [Ia1; Ia2; Ia0];三相短路就简单了只有正序网有电流短路电流等于短路前电压除以正序等值阻抗。程序里我实现了四种故障类型的选择开关三相短路LLL单相接地短路SLG两相短路LL两相接地短路LLG每种类型对应一个函数故障边界条件清晰失败概率低方便初学者读懂。5. 完整算例验证IEEE 14节点系统的潮流与故障结果理论讲得再多不如一个完整算例来得有说服力。我用IEEE 14节点标准测试系统做验证该系统包含5台同步发电机其中1号节点为平衡机、11条负荷支路、线路和变压器共20条。读入数据后基准功率设为100MVA电压基准取对应电压等级。先看潮流计算结果。收敛误差设为 (1 \times 10^{-8})标幺值牛顿-拉夫逊法在第4次迭代后收敛。各节点电压幅值范围在0.96到1.06之间这个结果和IEEE标准参考数据吻合。我把几个典型节点的结果列出来节点电压幅值 (pu)相角 (度)节点类型11.06000.00平衡节点21.0450-4.98PV节点41.0163-10.23PQ节点91.0302-13.78PQ节点140.9731-15.39PQ节点从调整角度看节点14电压偏低是因为该节点远离电源中心线路负载较重且无功支撑不足。这种电压分布在实际电网中也很典型如果想改善低压侧电压常用的手段是加并联电容器或者调整附近变压器的分接头。支路潮流方面线路1-2是主输电通道有功潮流从节点1送往节点2大约156MW线路1-5承担了第二大的输送功率。总网损约为13.3MW在IEEE 14节点系统中这个水平正常跟参考数据一致说明潮流程序计算正确。潮流收敛后接下来进行不对称短路分析。假设在节点9一个靠近负荷中心的节点分别发生四类故障短路前电压由潮流结果提供。以下是在节点9发生各类型故障时故障电流的幅值结果标幺值基准电流 (I_B 100MVA / (\sqrt{3} \times 220kV))故障类型故障相电流 (pu)短路容量 (MVA)三相短路12.36471.3单相接地A相7.82298.2两相短路B-C相10.71408.6两相接地B-C相接地11.55440.5可以看到三相短路电流最大这是多数情况下的规律因为三相短路是三相对称正序直接短路没有任何序网串联阻抗的削弱效应。两相接地短路电流次之单相接地短路电流最小。这个排序和理论分析完全一致反过来也验证了程序逻辑没问题。短路故障后各节点电压的跌落情况同样有意义。比如三相短路时故障点9的电压会降到接近0而相邻节点10、11的电压也会被拉低到0.3以下——这个信息对保护装置的灵敏度校验非常关键如果故障点近区的电压继电器定值设置不当可能在故障时拒动。代码中我把故障结果以极坐标形式输出同时打印故障相电流的实部和虚部方便和手工计算对照。写报告的时候我会把潮流结果和短路结果做成两个表格分别展示稳态和暂态的分析结论整篇论文的数据链路就是完整的。6. 我在调试中踩过的几个坑从迭代发散到结果错乱这个项目看上去公式都摆在那里但真正把代码跑通、跑对中间的坑多到可以单写一篇排错文章。我按照实际踩坑的频率来分享。第一个大坑是迭代不收敛。我最早是用平启动值迭代系统规模只有9节点理论上不可能有问题但程序就是发散。排查后发现是节点功率基准搞混了——发电机数据里给的是MW程序里计算功率不平衡量时忘记除以100MVA基准导致不平衡量比真实值大100倍雅可比矩阵修正步长过大迭代直接振荡。这种错误从报错信息上看不出来只能通过打印前几次迭代的电压变化量来发现。所以在编写数据预处理时一定要统一转换为标幺值。第二个坑是雅可比矩阵奇异。有一次把PV节点的无功不平衡行也强行写进修正方程矩阵行数和未知量数对不上导致求解失败。后来才意识到PV节点只参与有功修正它的电压幅值已知、无功待求所以雅可比矩阵中只有对应相角的列没有对应电压幅值的列。这个处理逻辑一开始我以为自己懂但写起代码来很容易在索引上栽跟头。解决办法是把节点分成两组所有非平衡节点的相角变量一组所有PQ节点的电压幅值变量一组然后按这个分组顺序构建雅可比矩阵的行列映射。第三个坑是支路数据里变比的符号约定。IEEE标准系统里变压器支路的变比通常标在支路某侧有的写法是首端/末端有的是末端/首端不统一。如果导纳矩阵构建函数和读入数据时变比定义颠倒了潮流结果可能看起来“能收敛”但支路潮流和标准结果偏差很大。我在程序里加了严格注释并且用一个极小的三节点手算案例做了导纳矩阵的单元测试——这个习惯救了我好几次。第四个坑和短路分析有关求三序等值阻抗时我一开始直接用序导纳矩阵求逆后来发现零序网络在有些节点上是不连通的零序导纳矩阵奇异inv直接报错。处理方法是先做一次稀疏LU分解同时检测主元是否接近0如果零序网在故障点附近不连通说明该处没有接地回路就不存在零序等值阻抗程序要给出明确提示而不是硬算。第五个坑是关于负荷的处理。潮流计算中负荷按恒定功率处理但在短路瞬间负荷通常等值为恒定阻抗因为电压跌落时真实负荷会变化。如果短路计算中仍然用潮流中的恒定功率负荷故障电流会算得偏大结果不真实。我在短路分析模块里加了一个选项可以把负荷节点按恒定阻抗折算到导纳矩阵中再参与序网计算。这个细节在实际工程整定计算里非常关键程序里默认开启了负荷阻抗等值。最后一个是数值精度。Matlab默认double精度对IEEE 14节点这样的小系统足够用但当你扩展到几百节点、并且迭代接近收敛时功率不平衡量可能出现小幅振荡不下降的情况。这通常是因为雅可比矩阵元素中某些项的值太小或者某些节点电压初值和真值偏离太大。遇到这种情况我一般先把收敛精度从1e-8放宽到1e-6确认结果稳定后再检查是不是个别节点类型设置有误。7. 如何扩展这个项目从教学工具到工程分析平台这个项目做完之后不光能交差实际上它的代码结构非常适合往多个方向扩展。最直接的扩展是增加故障切除和重合闸的时序仿真。目前的短路计算得到的是稳态短路电流实际工程还要关注断路器切除故障后系统能否保持稳定运行以及单相重合闸时系统的暂态行为。这需要引入故障时刻的时域仿真在这套代码的基础上可以把潮流结果作为故障前初值在故障点加入等值阻抗模拟故障状态然后用逐步积分法计算发电机功角的摇摆曲线。因为导纳矩阵部分已经是现成的只需要添加发电机转子运动方程和数值积分器。第二个扩展方向是连续潮流计算。把牛顿-拉夫逊求解器嵌入一个逐步增加负荷水平的循环中每次以上一步计算得到的电压解作为初值就能追踪PV曲线的鼻尖点也就是静态电压稳定的临界点。我在验证程序时发现只要在每次负荷增长后用上一步的潮流解作为初值收敛速度会非常快通常两三次迭代就能收敛。这个方向对研究电力系统电压稳定性很有价值。第三个扩展是接入优化算法。潮流计算本身是求解给定运行方式下的状态量但反过来可以问发电机出力和无功补偿如何配置才能让某个目标比如网损最小最优这时候潮流程序就变成了优化问题里的“状态约束”。我的做法是把潮流求解函数包装成一个约束函数然后用粒子群或者序列二次规划算法迭代调用它在IEEE 14节点系统上做了无功优化实验效果不错。Matlab的全局优化工具箱也能直接配合这套代码使用。第四个方向是把短路计算结果用于继电保护整定。节点故障时的电流电压结果可以反过来推算保护安装处的测量阻抗、测量电流进而绘出距离保护的动作圆、电流保护的配合阶梯。我在后续工作中还做了简单的距离保护定值校验把短路计算结果导入保护配合逻辑自动判断是否满足灵敏度要求。对于正在做课程设计的同学我比较推荐的实践路线是先跑通潮流计算把节点电压和支路潮流结果和标准数据核对然后把短路分析模块加进去用运行手册里的算例验证最后写报告时把关键公式推导、程序流程图、计算结果对比三部分配齐一篇有深度的设计报告自然就有了。那台电脑上至今还留着这个项目最初的版本代码写得不算漂亮但每一步都能看懂。工具在变系统在扩大但这些基本计算方法是一辈子都用得上的底子。
电力系统潮流计算与不对称短路分析：Matlab完整实现与工程实践
电力系统分析与仿真这个方向网上资料不少但大多要么堆理论公式、要么直接丢一个封装好的工具箱让读者跑黑箱。这次我把一个实际完成的“潮流计算不对称短路分析”Matlab项目完整拆开从数学模型怎么建立、雅可比矩阵怎么求到短路序网怎么连接、故障电流怎么算再到代码架构怎么设计、调试中有哪些坑一次讲清楚。适合正在做电力系统课程设计、毕业设计或者刚接触机电暂态分析、想从原理层面搞懂计算过程的同学参考。1. 为什么把潮流计算和不对称短路分析放在同一个工程里先聊点实际的在电力系统的规划、运行方式安排、保护整定这些日常工作里潮流计算和短路计算几乎是“搭班子”出现的两个计算工具。潮流计算解决的是“系统在某个稳态运行点下各节点电压、支路功率是多少”而短路计算解决的是“系统发生接地或者相间故障时故障电流有多大、各节点电压被拉低到什么程度”。前者决定了正常方式下的断面潮流和设备负载率后者决定了保护装置的动作电流整定值、断路器开断容量和设备动热稳定校验。既然两个计算都要基于同一套电网拓扑而且都要解网络方程那用一套代码把数据接口统一起来其实非常划算。我做这个项目时最直接的感受就是潮流计算需要节点导纳矩阵不对称短路的三序网络也要从原始导纳矩阵变换而来数据前处理几乎可以完全复用。短路计算需要一个“短路前”的运行状态——正好就是潮流计算的稳态解尤其是故障点处的电压幅值和相角会直接作为短路电流计算的初始条件。如果分开两个独立程序去做节点编号、基准值、支路参数很容易两边不一致最终对不上账。所以这个项目我直接用一套Matlab代码从读入节点和支路数据开始先做潮流再基于潮流结果做故障分析。整个程序分为三层数据层节点参数、支路参数、发电机参数、计算层导纳矩阵、潮流迭代、序网形成、故障解析、输出层电压分布表、支路潮流表、故障电流相量结果。这样模块之间的耦合度低调其中一个环节不会把别的部分搞坏。工程上常见的做法是用BPA、PSS/E这类商业软件直接出结果但它们的计算内核是封闭的不方便做二次开发和教学演示。用Matlab重写一遍麻烦归麻烦却能把每一步计算过程都暴露出来查问题、改算法、接入新的控制策略都容易得多。这个项目做完之后我甚至把里面的潮流求解函数单独拆出来接上了一个简单的连续潮流CPF扩展用来算PV曲线的鼻尖点这也算是模块化设计带来的额外红利。2. 潮流计算的数学本质我们到底在解什么方程潮流计算从数学上看就是求解一组非线性代数方程组。电力网络用节点导纳矩阵描述后每个节点的注入电流可以写成[ I_i \sum_{j1}^{n} Y_{ij} U_j ]其中 (Y_{ij}) 是导纳矩阵元素(U_j) 是节点电压相量。节点功率和电压、电流的关系为[ S_i U_i I_i^* P_i jQ_i ]把电流表达式代进去就得到了潮流方程的基本形式[ P_i jQ_i U_i \sum_{j1}^{n} (Y_{ij} U_j)^* ]展开成实部和虚部之后每个节点有两个功率平衡方程一个对应有功 (P)一个对应无功 (Q)。问题在于每个节点的状态变量是电压幅值 (U_i) 和相角 (\theta_i)也是两个未知数。系统中必须指定一个平衡节点slack bus承担功率差额它的电压幅值和相角固定其他节点的变量才有唯一解。节点类型惯例是这样的PQ节点给定有功和无功注入负荷节点通常归入此类PV节点给定有功注入和电压幅值发电机节点归入此类平衡节点给定电压幅值和相角选一个容量较大、接近系统中心的发电机节点计算量最集中的地方是求解修正方程。牛顿-拉夫逊法每一步需要计算不平衡量 (\Delta P)、(\Delta Q)以及雅可比矩阵然后求解线性方程组得到电压修正量。雅可比矩阵有四个分块[ J \begin{bmatrix} H N \ K L \end{bmatrix} ]其中 (H_{ij} \frac{\partial P_i}{\partial \theta_j})(N_{ij} \frac{\partial P_i}{\partial U_j})(K_{ij} \frac{\partial Q_i}{\partial \theta_j})(L_{ij} \frac{\partial Q_i}{\partial U_j})。具体表达式对非对角元素和对角元素是不同的。非对角元素(i \neq j)[ H_{ij} -U_i U_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij}) ][ N_{ij} -U_i U_j (G_{ij} \cos\theta_{ij} B_{ij} \sin\theta_{ij}) ][ K_{ij} U_i U_j (G_{ij} \cos\theta_{ij} B_{ij} \sin\theta_{ij}) ][ L_{ij} -U_i U_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij}) ]对角元素则要考虑自导纳带来的附加项。这些公式看似繁琐但规律很强只要把导纳矩阵和当前迭代的电压值代入就能批量生成雅可比矩阵。我用一个双层循环加if语句判断 (i) 和 (j) 是否相等来处理在Matlab里可读性很好如果要追求高性能工程实现用向量化或者稀疏矩阵构造会更快但教学和中等规模系统分析循环写法完全够用。潮流计算的收敛判据工程上一般取功率不平衡量最大绝对值小于一个阈值比如 (1 \times 10^{-6}) 或 (1 \times 10^{-8})标幺值。迭代初值大多数情况用“平启动”所有PQ节点电压幅值设1.0相角设0PV节点相角也设0平衡节点直接固定。对于绝大多数110kV及以上电网变压器分接头位置合适、无功补偿到位时平启动基本都能在5到8次迭代内收敛。3. Matlab代码实现数据准备、导纳矩阵与牛顿-拉夫逊迭代核心代码实现的第一步是数据准备。我习惯用结构体存储电网数据比较直观% 节点数据 % 节点编号 类型(1PQ, 2PV, 3平衡) 电压幅值初值 相角初值(度) Pg(MW) Qg(Mvar) Pd(MW) Qd(Mvar) busData [ 1 3 1.060 0 0 0 0 0 2 2 1.045 0 40 30 21.7 12.7 3 2 1.010 0 0 0 94.2 19.0 4 1 1.000 0 0 0 47.8 -3.9 % ... ]; % 支路数据 % 首端节点 末端节点 电阻R(pu) 电抗X(pu) 电纳B/2(pu) 变比k branchData [ 1 2 0.01938 0.05917 0.0264 1 1 5 0.05403 0.22304 0.0246 1 2 3 0.04699 0.19797 0.0219 1 % ... ];接下来构建节点导纳矩阵。这一步的物理含义很清晰自导纳是该节点连接的所有支路导纳之和互导纳是相连支路导纳的相反数。变压器支路要注意变比折算function Y buildYbus(branchData, busData) n size(busData, 1); Y zeros(n, n); nBranch size(branchData, 1); for k 1:nBranch i branchData(k, 1); j branchData(k, 2); r branchData(k, 3); x branchData(k, 4); b branchData(k, 5); tap branchData(k, 6); z r 1j * x; y 1 / z; yij y / tap; ysi y * (tap - 1) / tap 1j * b / 2; ysj y * (1 - tap) / (tap^2) 1j * b / 2; Y(i, i) Y(i, i) yij ysi; Y(j, j) Y(j, j) yij ysj; Y(i, j) Y(i, j) - yij; Y(j, i) Y(j, i) - yij; end end注意这里的 (tap) 我定义为非标准变比即从节点 i 侧看理想变压器变比为 (1:tap)。实际工程中变压器分接头不同公式略有差异建议在代码注释里标明约定免得日后自己都忘了。接下来是牛顿-拉夫逊迭代主体。为了不把所有代码都堆在一个脚本里我拆成三个函数calPower计算当前电压下的注入功率、calJacobian形成雅可比矩阵、updateVoltage解方程并更新电压相量。主迭代流程如下V busData(:, 3) .* exp(1j * deg2rad(busData(:, 4))); for iter 1:30 % 计算不平衡量 [Pcal, Qcal] calPower(V, Ybus); dP (busData(:, 5) - busData(:, 7)) - Pcal; % MW - pu需除以基准功率 dQ (busData(:, 6) - busData(:, 8)) - Qcal; % 剔除平衡节点和PV节点的无功不平衡量行 % 构造雅可比矩阵 J calJacobian(V, Ybus); % 求解修正方程 dTheta J \ dPQ; % 更新电压相角和幅值 theta theta dTheta(1:n-1); Vmag Vmag dTheta(n:end); if max(abs(dPQ)) 1e-8 break; end end需要提醒的是雅可比矩阵的行列要按系统中PQ节点和PV节点的排列来组织工程上常用两个数组记录节点类型索引比每次在循环里用if判断快得多。迭代过程中要限制PV节点的无功出力范围——如果某台发电机的无功越限应该把它降级为PQ节点处理否则潮流结果虽然收敛但物理上不可行。这个处理在程序里表现为动态修改节点类型重新排列雅可比矩阵的行列我在实现时单独写了一个子函数refreshNodeType每次迭代前检查一次。潮流结果输出时我喜欢把数据直接整理成表格形式打印出来方便写报告时直接引用。除了各节点电压幅值、相角还需要计算各支路潮流首端注入、末端注入、网损这个计算可以复用上面的功率计算公式只不过把注入功率换成某一条支路两端各自算一遍。4. 不对称短路分析的核心三序网络与复合序网潮流计算解决的是对称稳态问题而短路计算面临的是不对称情况这时不能直接在ABC三相的混合相量里求解因为三相阻抗耦合在一起电路方程不独立。工程上通用的做法是输出对称分量法把三相电气量分解成三组对称分量正序、负序、零序。对称分量的核心变换是[ \begin{bmatrix} \dot{F}_a \ \dot{F}_b \ \dot{F}_c \end{bmatrix}\begin{bmatrix} 1 1 1 \ 1 a^2 a \ 1 a a^2 \end{bmatrix} \begin{bmatrix} \dot{F}{a1} \ \dot{F}{a2} \ \dot{F}_{a0} \end{bmatrix} ]其中 (a e^{j120^\circ} -0.5 j0.866)。正序分量是三相对称、按正常相序排列的电压电流负序分量是按逆相序排列的零序分量则是三相幅值、相位都相同的分量。经过这个变换一个不对称的三相网络被拆成了三个各自对称的序网络互不耦合只在故障点通过故障边界条件连在一起。在实际的Matlab程序里只要生成一个变换矩阵T三序电压电流的计算就是矩阵乘法不需要手工按公式硬套。这一点简化了很多操作。比如求某一节点A相电压的序分量只需要T [1 1 1; 1 a^2 a; 1 a a^2]; % 注意这里的a exp(1j*2*pi/3) Fabc [Va; Vb; Vc]; F012 inv(T) * Fabc;不过短路分析更常用的思路不是从三相电压推序分量而是直接建立正序、负序、零序网络的节点导纳矩阵然后在故障节点处按短路类型把三个序网连起来。正序网络就是常规潮流用的导纳矩阵当然发电机、负荷模型在短路瞬间有特殊等值通常正序网要用次暂态电抗等值负荷可以近似忽略或接恒定阻抗。负序网络和正序网络拓扑相同只是发电机负序电抗和正序不同线路参数相同。零序网络则只包含能够为零序电流提供回路的路径比如变压器绕组接线方式YN接地与否、三角形绕组内部零序环流会直接决定零序网络的结构。复合序网法不需要手工推导各种短路情形的公式而是把故障点看成三序网络的连接口单相接地短路时三序网在故障点串联两相短路时正序网和负序网在故障点并联零序网不参与两相接地短路时三个序网在故障点并联。有了这个连接方式就可以把三个序网络方程和故障边界条件统一成一个复合线性方程组直接求解故障点的正、负、零序电流。工程上常用短路点电压形式的复合序网方程。比如单相接地短路A相的边界条件是[ \begin{cases} \dot{I}_b \dot{I}c 0 \ \dot{V}{ka} 0 \end{cases} ]等值到序分量后条件是[ \dot{I}{a1} \dot{I}{a2} \dot{I}{a0} \frac{\dot{V}{f}}{Z_1 Z_2 Z_0} ]其中 (\dot{V}_f) 是短路前故障点电压由潮流计算得到(Z_1, Z_2, Z_0) 是从故障点看进去的三序等值阻抗。这里要特别注意 (Z_1) 是正序网中故障点对地等值阻抗它和潮流计算得到的系统等值电势与故障点电压之间的阻抗有关实际求法是在故障点注入单位正序电流测该点电压两者相除就是正序等值阻抗。在代码实现上用阻抗矩阵法最省事。先对三序节点导纳矩阵分别求逆或者用LU分解解线性方程得到三序阻抗矩阵然后取故障节点对应的对角元素就是该点的三序等值阻抗。这个过程快速又不容易错Zseq{1} inv(Y1bus); % 正序阻抗矩阵 Zseq{2} inv(Y2bus); % 负序阻抗矩阵 Zseq{0} inv(Y0bus); % 零序阻抗矩阵 z1 Zseq{1}(k, k); z2 Zseq{2}(k, k); z0 Zseq{0}(k, k);得到三序等值阻抗后短路电流的计算就变成了简单的复数运算。拿到故障点序电流之后再用反变换得到三相短路电流Iabc T * [Ia1; Ia2; Ia0];三相短路就简单了只有正序网有电流短路电流等于短路前电压除以正序等值阻抗。程序里我实现了四种故障类型的选择开关三相短路LLL单相接地短路SLG两相短路LL两相接地短路LLG每种类型对应一个函数故障边界条件清晰失败概率低方便初学者读懂。5. 完整算例验证IEEE 14节点系统的潮流与故障结果理论讲得再多不如一个完整算例来得有说服力。我用IEEE 14节点标准测试系统做验证该系统包含5台同步发电机其中1号节点为平衡机、11条负荷支路、线路和变压器共20条。读入数据后基准功率设为100MVA电压基准取对应电压等级。先看潮流计算结果。收敛误差设为 (1 \times 10^{-8})标幺值牛顿-拉夫逊法在第4次迭代后收敛。各节点电压幅值范围在0.96到1.06之间这个结果和IEEE标准参考数据吻合。我把几个典型节点的结果列出来节点电压幅值 (pu)相角 (度)节点类型11.06000.00平衡节点21.0450-4.98PV节点41.0163-10.23PQ节点91.0302-13.78PQ节点140.9731-15.39PQ节点从调整角度看节点14电压偏低是因为该节点远离电源中心线路负载较重且无功支撑不足。这种电压分布在实际电网中也很典型如果想改善低压侧电压常用的手段是加并联电容器或者调整附近变压器的分接头。支路潮流方面线路1-2是主输电通道有功潮流从节点1送往节点2大约156MW线路1-5承担了第二大的输送功率。总网损约为13.3MW在IEEE 14节点系统中这个水平正常跟参考数据一致说明潮流程序计算正确。潮流收敛后接下来进行不对称短路分析。假设在节点9一个靠近负荷中心的节点分别发生四类故障短路前电压由潮流结果提供。以下是在节点9发生各类型故障时故障电流的幅值结果标幺值基准电流 (I_B 100MVA / (\sqrt{3} \times 220kV))故障类型故障相电流 (pu)短路容量 (MVA)三相短路12.36471.3单相接地A相7.82298.2两相短路B-C相10.71408.6两相接地B-C相接地11.55440.5可以看到三相短路电流最大这是多数情况下的规律因为三相短路是三相对称正序直接短路没有任何序网串联阻抗的削弱效应。两相接地短路电流次之单相接地短路电流最小。这个排序和理论分析完全一致反过来也验证了程序逻辑没问题。短路故障后各节点电压的跌落情况同样有意义。比如三相短路时故障点9的电压会降到接近0而相邻节点10、11的电压也会被拉低到0.3以下——这个信息对保护装置的灵敏度校验非常关键如果故障点近区的电压继电器定值设置不当可能在故障时拒动。代码中我把故障结果以极坐标形式输出同时打印故障相电流的实部和虚部方便和手工计算对照。写报告的时候我会把潮流结果和短路结果做成两个表格分别展示稳态和暂态的分析结论整篇论文的数据链路就是完整的。6. 我在调试中踩过的几个坑从迭代发散到结果错乱这个项目看上去公式都摆在那里但真正把代码跑通、跑对中间的坑多到可以单写一篇排错文章。我按照实际踩坑的频率来分享。第一个大坑是迭代不收敛。我最早是用平启动值迭代系统规模只有9节点理论上不可能有问题但程序就是发散。排查后发现是节点功率基准搞混了——发电机数据里给的是MW程序里计算功率不平衡量时忘记除以100MVA基准导致不平衡量比真实值大100倍雅可比矩阵修正步长过大迭代直接振荡。这种错误从报错信息上看不出来只能通过打印前几次迭代的电压变化量来发现。所以在编写数据预处理时一定要统一转换为标幺值。第二个坑是雅可比矩阵奇异。有一次把PV节点的无功不平衡行也强行写进修正方程矩阵行数和未知量数对不上导致求解失败。后来才意识到PV节点只参与有功修正它的电压幅值已知、无功待求所以雅可比矩阵中只有对应相角的列没有对应电压幅值的列。这个处理逻辑一开始我以为自己懂但写起代码来很容易在索引上栽跟头。解决办法是把节点分成两组所有非平衡节点的相角变量一组所有PQ节点的电压幅值变量一组然后按这个分组顺序构建雅可比矩阵的行列映射。第三个坑是支路数据里变比的符号约定。IEEE标准系统里变压器支路的变比通常标在支路某侧有的写法是首端/末端有的是末端/首端不统一。如果导纳矩阵构建函数和读入数据时变比定义颠倒了潮流结果可能看起来“能收敛”但支路潮流和标准结果偏差很大。我在程序里加了严格注释并且用一个极小的三节点手算案例做了导纳矩阵的单元测试——这个习惯救了我好几次。第四个坑和短路分析有关求三序等值阻抗时我一开始直接用序导纳矩阵求逆后来发现零序网络在有些节点上是不连通的零序导纳矩阵奇异inv直接报错。处理方法是先做一次稀疏LU分解同时检测主元是否接近0如果零序网在故障点附近不连通说明该处没有接地回路就不存在零序等值阻抗程序要给出明确提示而不是硬算。第五个坑是关于负荷的处理。潮流计算中负荷按恒定功率处理但在短路瞬间负荷通常等值为恒定阻抗因为电压跌落时真实负荷会变化。如果短路计算中仍然用潮流中的恒定功率负荷故障电流会算得偏大结果不真实。我在短路分析模块里加了一个选项可以把负荷节点按恒定阻抗折算到导纳矩阵中再参与序网计算。这个细节在实际工程整定计算里非常关键程序里默认开启了负荷阻抗等值。最后一个是数值精度。Matlab默认double精度对IEEE 14节点这样的小系统足够用但当你扩展到几百节点、并且迭代接近收敛时功率不平衡量可能出现小幅振荡不下降的情况。这通常是因为雅可比矩阵元素中某些项的值太小或者某些节点电压初值和真值偏离太大。遇到这种情况我一般先把收敛精度从1e-8放宽到1e-6确认结果稳定后再检查是不是个别节点类型设置有误。7. 如何扩展这个项目从教学工具到工程分析平台这个项目做完之后不光能交差实际上它的代码结构非常适合往多个方向扩展。最直接的扩展是增加故障切除和重合闸的时序仿真。目前的短路计算得到的是稳态短路电流实际工程还要关注断路器切除故障后系统能否保持稳定运行以及单相重合闸时系统的暂态行为。这需要引入故障时刻的时域仿真在这套代码的基础上可以把潮流结果作为故障前初值在故障点加入等值阻抗模拟故障状态然后用逐步积分法计算发电机功角的摇摆曲线。因为导纳矩阵部分已经是现成的只需要添加发电机转子运动方程和数值积分器。第二个扩展方向是连续潮流计算。把牛顿-拉夫逊求解器嵌入一个逐步增加负荷水平的循环中每次以上一步计算得到的电压解作为初值就能追踪PV曲线的鼻尖点也就是静态电压稳定的临界点。我在验证程序时发现只要在每次负荷增长后用上一步的潮流解作为初值收敛速度会非常快通常两三次迭代就能收敛。这个方向对研究电力系统电压稳定性很有价值。第三个扩展是接入优化算法。潮流计算本身是求解给定运行方式下的状态量但反过来可以问发电机出力和无功补偿如何配置才能让某个目标比如网损最小最优这时候潮流程序就变成了优化问题里的“状态约束”。我的做法是把潮流求解函数包装成一个约束函数然后用粒子群或者序列二次规划算法迭代调用它在IEEE 14节点系统上做了无功优化实验效果不错。Matlab的全局优化工具箱也能直接配合这套代码使用。第四个方向是把短路计算结果用于继电保护整定。节点故障时的电流电压结果可以反过来推算保护安装处的测量阻抗、测量电流进而绘出距离保护的动作圆、电流保护的配合阶梯。我在后续工作中还做了简单的距离保护定值校验把短路计算结果导入保护配合逻辑自动判断是否满足灵敏度要求。对于正在做课程设计的同学我比较推荐的实践路线是先跑通潮流计算把节点电压和支路潮流结果和标准数据核对然后把短路分析模块加进去用运行手册里的算例验证最后写报告时把关键公式推导、程序流程图、计算结果对比三部分配齐一篇有深度的设计报告自然就有了。那台电脑上至今还留着这个项目最初的版本代码写得不算漂亮但每一步都能看懂。工具在变系统在扩大但这些基本计算方法是一辈子都用得上的底子。
