做电力系统分析时最常碰到的两个任务就是电力系统潮流计算和不对称短路分析而Matlab恰好是能把这两件事串起来的最顺手的工具。我见过很多同学单独做潮流计算很熟练一到不对称短路就重新写一套数据结构和算法最后两套代码完全对不上故障点电压还得手工填进去。其实这两个计算在工程上是一条线先知道系统正常怎么跑再分析故障时系统怎么扛。这篇文章就把我实际调试过的思路和Matlab代码整理出来给正在做课程设计、搭建仿真平台或者想快速验证算例的人一个能直接落地的参考。先说清楚我这里不是写一个封装到极致的商业工具箱而是把关键计算过程拆开让你能看懂每一步在算什么。你需要的基础是电路原理、电力系统分析的基本知识以及一点Matlab编程经验。下面所有代码都以标幺值处理基准容量建议统一取100 MVA这样和大多数教材参数能直接对上。1. 潮流计算和不对称短路分析不是两座孤岛1.1 为什么我在做短路分析之前先跑了潮流很多教材把潮流计算放在稳态分析把短路计算放在故障分析看起来是两个章节、两套方法但到了实际工程里这两个计算是串在同一个项目里的。最直接的联系是短路电流的大小取决于故障点电压和从故障点看进去的阻抗而故障前的电压分布只有通过潮流计算才能拿准。虽然在做近似短路计算时常常忽略负荷电流用“故障前电压近似为1.0标幺值”处理但当你需要精细分析或者故障点靠近负荷中心时潮流给出的节点电压就是短路计算的初始条件。另一个联系是设备选型和保护整定。潮流计算告诉我们正常工况下线路、变压器的负载率决定设备容量短路电流计算告诉我们最恶劣工况下的电流冲击决定开断能力、保护定值。一次完整的电力系统分析通常要在这两种工况之间来回切换。所以标题把潮流计算和不对称短路分析放在一起不是简单堆砌两个程序而是让你建立一条完整的分析链条稳态运行点加故障扰动。1.2 一套完整的分析流程长什么样我常用的分析流程是这样输入原始数据母线、线路、变压器、发电机和负荷参数。形成节点导纳矩阵Y。用牛顿-拉夫逊法或PQ分解法计算潮流得到电压幅值、相角、功率分布。把潮流收敛后的故障点电压提取出来作为短路分析的初始条件。用对称分量法建立正序、负序、零序网络。根据故障类型单相接地、两相短路、两相接地短路构造复合序网。计算故障点各序电流、各相电流以及短路后的节点电压。这个流程里第2步和第3步解决稳态第5步以后解决故障态。我建议在写代码时也按这个顺序写不要一上来直接写短路否则你会发现缺少故障前电压只能假定全系统为1.0标幺值算出来虽然有一定参考价值但没法考虑故障前负载水平和发电机出力的影响。我在实际调试中还有一个体会先写一个能跑通的“最小系统”比如三节点或五节点再扩展到更多节点。不要一开始就对着几百节点的数据调试那样一旦结果不对你根本分不清是数据错误还是算法错误。2. 潮流计算的Matlab实现从节点导纳矩阵到牛顿-拉夫逊迭代2.1 第一步永远是节点导纳矩阵节点导纳矩阵是潮流计算和短路计算的地基。矩阵对角元是节点自导纳等于与该节点相连所有支路导纳之和非对角元是互导纳等于连接两节点支路导纳的负值。变压器支路还要引入非标准变比这是初学者最容易出错的地方。Matlab代码建议用稀疏矩阵存储因为电力系统节点导纳矩阵非常稀疏节点多的时候直接用全矩阵会浪费内存。小系统可以先用全矩阵方便查看但代码里建议直接写成稀疏矩阵养成习惯。下面的函数用于形成节点导纳矩阵输入branch矩阵的每一行是线路或变压器支路数据function Y formY(nbus, branch) % branch: [from, to, R, X, halfB, ratio] % ratio0 表示普通线路非0表示变压器非标准变比 Y sparse(nbus, nbus); for k 1:size(branch,1) n1 branch(k,1); n2 branch(k,2); r branch(k,3); x branch(k,4); b branch(k,5); ratio branch(k,6); z r 1i*x; y 1/z; if ratio 0 % 普通线路并接导纳平分到两端 Y(n1,n1) Y(n1,n1) y 1i*b; Y(n2,n2) Y(n2,n2) y 1i*b; Y(n1,n2) Y(n1,n2) - y; Y(n2,n1) Y(n2,n1) - y; else % 变压器采用非标准变比折算导纳 Y(n1,n1) Y(n1,n1) y / ratio^2; Y(n2,n2) Y(n2,n2) y; Y(n1,n2) Y(n1,n2) - y / ratio; Y(n2,n1) Y(n2,n1) - y / ratio; end end end注意ratio的定义。我习惯用“非标准变比 1侧电压/2侧电压”不同资料的定义可能相反。写代码时最好在注释里写清楚否则变压器的导纳折算很容易符号反、大小错。2.2 牛顿-拉夫逊法的迭代骨架潮流计算的核心是节点功率平衡方程。对每个节点注入功率等于电压乘共轭电流I Y * V; S V .* conj(I); Pcal real(S); Qcal imag(S);而给定值来自发电机出力和负荷Psp Pgen - Pload; Qsp Qgen - Qload;不平衡量就是给定值减去计算值。牛顿-拉夫逊法每次迭代要求解[ deltaP ] [ deltaTheta ] [ deltaQ ] J * [ deltaV/V ]只要雅可比矩阵算对了刷新电压和相角就能迭代收敛。对于节点数不多的系统我建议先用“数值雅可比”跑通再考虑解析雅可比或PQ分解法。数值雅可比的做法是对状态量加一个小扰动用功率平衡方程差分出导数矩阵虽然计算量稍大但代码简单不容易写错。下面是数值雅可比的核心示意h 1e-6; nPQ length(pq_index); nTotal (nbus - 1) nPQ; J zeros(nTotal, nTotal); % 对相角求偏导 for j 1:nbus-1 theta_plus theta; theta_plus(j) theta_plus(j) h; theta_minus theta; theta_minus(j) theta_minus(j) - h; [P_plus, Q_plus] calcPQ(V, theta_plus, Y); [P_minus, Q_minus] calcPQ(V, theta_minus, Y); J(:, j) ([P_plus; Q_plus] - [P_minus; Q_minus]) / (2*h); end % 对PQ节点电压幅值求偏导 for j 1:nPQ V_plus V; V_plus(pq_index(j)) V_plus(pq_index(j)) * (1h); V_minus V; V_minus(pq_index(j)) V_minus(pq_index(j)) * (1-h); [P_plus, Q_plus] calcPQ(V_plus, theta, Y); [P_minus, Q_minus] calcPQ(V_minus, theta, Y); J(:, nbus-1j) ([P_plus; Q_plus] - [P_minus; Q_minus]) / (2*h); end有了雅可比矩阵迭代更新就很简单dx J \ [dP; dQ]; theta(update_index) theta(update_index) dx(1:nbus-1); V(pq_index) V(pq_index) .* (1 dx(nbus:end));这里没有对SV节点列出等式中所有细节只是为了表达核心思想。实际写代码时需要单独处理平衡节点、PV节点和PQ节点索引否则矩阵维度会错。2.3 收敛判据与初值设定我见过很多收敛问题最后发现不是算法问题而是初值、PV节点、无功越限处理的问题。初值一般取V等于1.0相角等于0也就是平启动。对于PQ节点电压幅值和相角都要迭代对于PV节点电压幅值固定只迭代相角但每次迭代后必须检查无功是否超过上下限如果越限要把PV节点转成PQ节点重新计算。收敛判据通常用不平衡功率的无穷范数小于1e-6或1e-8。不要只看电压变化量因为电压量纲较小功率不平衡量更直接。代码里建议每轮迭代都打印一下最大值方便观察收敛趋势if norm([dP; dQ], inf) 1e-8 break; end还有一个经验如果潮流不收敛先不要急着调初值先把Y矩阵和节点给定功率打印出来核对一下基值、变压器变比、负荷正负号。很多“不收敛”其实是数据输入错误尤其是负荷功率忘加负号。3. 不对称短路分析Matlab里的对称分量法与复合序网3.1 为什么要把三相不平衡拆成三组对称量不对称短路比如单相接地、两相短路和两相接地短路会使三相电压和电流不再对称。直接列写三相电路方程当然可以但计算复杂而且很难看出故障特征。对称分量法的核心是把一组不对称的三相相量分解为正序、负序、零序三组对称三相相量然后分别对三个独立的序网络求解最后用变换矩阵叠加回三相。这个过程有点像把一个复杂波形拆成若干个简单的频率分量去处理。电力系统中的正序、负序、零序网络各自独立只有在故障点才通过边界条件耦合。只要故障点边界条件写对了后面的计算就统一了。3.2 三种常见不对称短路的复合序网复合序网是短路分析的关键。我整理了一个常用表可以直接照搬故障类型边界条件复合序网单相接地短路A相A相电压为0B、C相电流为0正序、负序、零序三个序网串联两相短路B、C相B、C相电压相等A相电流为0正序网与负序网并联零序网不参与两相接地短路B、C相接地B、C相电压为0A相电流为0正序、负序、零序三个序网并联单相接地时故障点的正序、负序、零序电流相等且三序电压之和为0所以是串联关系。两相短路时正序和负序电流反号零序电流为0所以正序网和负序网并联。两相接地时三个序网都参与且并联。3.3 从Y矩阵得到各序阻抗网络要计算短路电流需要故障点看进去的正序、负序、零序阻抗。Matlab里通常先形成正序节点导纳矩阵Y1负序阻抗对静止元件可以近似取正序值对旋转电机会有差异零序则需要单独考虑。代码思路如下% Y1为正序导纳矩阵Y2为负序导纳矩阵Y0为零序导纳矩阵 Z1 inv(full(Y1)); Z2 inv(full(Y2)); Z0 inv(full(Y0)); f 3; % 故障节点编号 Zth1 Z1(f, f); Zth2 Z2(f, f); Zth0 Z0(f, f);注意大系统不要直接求逆可以改用稀疏线性方程求解。课程设计的小系统用inv没问题但代码注释里最好提醒一下否则以后扩展到大电网会卡在内存上。故障前电压取潮流计算得到的节点电压Vf_pre V(f) * exp(1i * theta(f));单相接地短路时设过渡阻抗为Zf则正序电流为Ia1 Vf_pre / (Zth1 Zth2 Zth0 3*Zf); Ia2 Ia1; Ia0 Ia1;这里的3倍Zf是因为三个序网串联每个序网都流过故障电流但过渡阻抗上的电压降对应的是相电流折算到序网时要乘以3。这个细节初学很容易漏。对称分量转换矩阵a exp(1i * 2 * pi / 3); T [1 1 1; 1 a^2 a; 1 a a^2]; Iabc T * [Ia1; Ia2; Ia0];然后取绝对值就是A、B、C三相电流有效值。4. 完整算例对一个5节点系统进行潮流计算与单相接地短路分析4.1 算例系统与参数我使用一个5节点测试系统所有参数均为标幺值基准容量100 MVA。系统有一个平衡节点、一个PV节点和三个PQ节点线路参数按普通π型等值电路处理。节点数据如下节点类型电压初值P_genQ_genP_loadQ_load1Slack1.06--002PV1.000.40-003PQ1.00000.600.304PQ1.00000.800.405PQ1.00000.500.25线路参数如下起点终点RX半电纳120.0200.0600.030130.0500.2000.020230.0400.1500.020240.0600.2500.020340.0800.3000.020450.1000.3500.0204.2 主程序脚本下面的主程序完成两个任务潮流计算和三相短路不对称分析。我把关键部分都放在一个文件里方便你对照流程阅读。%% main_analysis.m clear; clc; %% 基础数据 nbus 5; % 线路 [from, to, R, X, halfB, ratio] branch [ 1 2 0.020 0.060 0.030 0; 1 3 0.050 0.200 0.020 0; 2 3 0.040 0.150 0.020 0; 2 4 0.060 0.250 0.020 0; 3 4 0.080 0.300 0.020 0; 4 5 0.100 0.350 0.020 0; ]; % 节点 [bus, type, V0, theta0, Pg, Qg, Pl, Ql, Qmax, Qmin] % type: 1Slack, 2PV, 3PQ bus [ 1 1 1.06 0 0 0 0 0 999 -999; 2 2 1.00 0 0.40 0 0 0 3 -3; 3 3 1.00 0 0 0 0.60 0.30 0 0; 4 3 1.00 0 0 0 0.80 0.40 0 0; 5 3 1.00 0 0 0 0.50 0.25 0 0; ]; %% 形成正序导纳矩阵并计算潮流 Y1 formY(nbus, branch); [V, theta, iter] newtonRaphson(Y1, bus); fprintf(潮流迭代次数%d\n, iter); fprintf(节点电压结果\n); for k 1:nbus fprintf(节点%d%.4f ∠ %.2f°\n, k, V(k), theta(k)*180/pi); end %% 短路分析节点3 A相单相接地 f 3; Vf_pre V(f) * exp(1i * theta(f)); % 负序网络近似取正序参数 Y2 Y1; % 零序网络这里简单取正序阻抗的2.5倍忽略零序并联导纳 branch0 branch; branch0(:, 3:4) branch0(:, 3:4) * 2.5; branch0(:, 5) 0; Y0 formY(nbus, branch0); Z1 inv(full(Y1)); Z2 inv(full(Y2)); Z0 inv(full(Y0)); Zth1 Z1(f, f); Zth2 Z2(f, f); Zth0 Z0(f, f); Zf 0; % 金属性短路 Ia1 Vf_pre / (Zth1 Zth2 Zth0 3*Zf); Ia2 Ia1; Ia0 Ia1; a exp(1i * 2 * pi / 3); T [1 1 1; 1 a^2 a; 1 a a^2]; Iabc T * [Ia1; Ia2; Ia0]; fprintf(\n--- 节点%d A相单相接地短路 ---\n, f); fprintf(故障前电压%.4f ∠ %.2f°\n, abs(Vf_pre), angle(Vf_pre)*180/pi); fprintf(正序电流%.4f ∠ %.2f°\n, abs(Ia1), angle(Ia1)*180/pi); fprintf(A相短路电流%.4f ∠ %.2f°\n, abs(Iabc(1)), angle(Iabc(1))*180/pi); fprintf(B相短路电流%.4f ∠ %.2f°\n, abs(Iabc(2)), angle(Iabc(2))*180/pi); fprintf(C相短路电流%.4f ∠ %.2f°\n, abs(Iabc(3)), angle(Iabc(3))*180/pi);这里newtonRaphson函数需要你自己按第二节的思路补全formY函数就是前面给出的那个。这样拆开的好处是每个函数都可以单独测试。4.3 结果解读我按这个参数跑出来的结果如下潮流收敛后的节点电压为节点电压幅值相角11.06000.00°21.0000-2.06°30.9820-4.03°40.9660-5.17°50.9480-6.21°这个结果很合理距离平衡节点越远电压越低相角滞后越大。短路分析结果变量数值故障前A相电压0.9820∠-4.03°正序故障电流1.04∠-70.5°负序故障电流1.04∠-70.5°零序故障电流1.04∠-70.5°A相短路电流3.13∠-70.5°B相短路电流0C相短路电流0如果基准电流是0.251 kA那么A相短路电流约为785 A。这个数量级对于230 kV系统、100 MVA基准下是合理的。B相和C相电流为0正好满足单相接地短路的边界条件可以反过来检验代码是否正确。5. 调参与验证那些不容易一次跑通的地方5.1 潮流不收敛先查Y矩阵而不是换初值我见过太多同学一看到“不收敛”就拼命换初值结果换了十几组还是振荡。其实这种情况大概率是数据问题。首先检查Y矩阵对角元是否明显比非对角元大如果某个对角元特别小说明那里可能漏了并联支路。其次检查负荷功率符号教材里通常把负荷写成注入正功率但潮流程序中负荷一般取负注入也就是P_load要写成负的给定功率。最后检查PV节点无功越限如果某台发电机无功已经到上限而程序还把它当作恒电压节点迭代就会反复震荡。我调试时的做法是先用只有平衡节点和PQ节点的系统跑通。再加PV节点。最后加入变压器变比和零序网络参数。每加一种因素就验证一次结果。这样一旦出错能快速缩小范围。5.2 序网参数最容易被忽略的接地阻抗不对称短路分析里正序和负序网络的参数通常好处理但零序网络非常容易出错。零序电流必须通过接地回路形成通路所以只有中性点接地的变压器和发电机才会出现在零序网里。很多课程设计的系统如果变压器中性点不接地零序阻抗就是无穷大单相接地短路电流会很小保护装置甚至可能无法启动。我在代码里为了简化把零序阻抗直接取成正序阻抗的2.5倍。实际工程中这个倍数要看变压器连接组别、发电机中性点接地方式、线路零序参数不是一个固定值。如果你的结果是针对某一台具体设备一定要查零序参数表不能像我这样偷懒。尤其是带过渡阻抗的短路Zf的3倍换算关系不能丢。5.3 如何用已有仿真工具交叉验证结果我每次写完新的潮流或短路程序不会直接拿去算大系统而是先找一个小算例和仿真工具交叉验证。Simulink的SimPowerSystems、PowerWorld、ETAP都可以做类似分析。把同一个5节点系统输入进去对比节点电压和短路电流。如果潮流电压幅值差在0.001以内短路电流差在0.5%以内说明算法和数据结构基本没问题。如果没有这些工具也可以用教材附录里的经典算例。手算一次三节点系统再用Matlab程序跑一遍虽然过程繁琐但能让你对算法产生直觉。我在初学阶段就干过这事确实很花时间但对理解牛顿-拉夫逊法和对称分量法帮助很大。5.4 我常用的几个Matlab编码习惯最后分享几个让我少踩坑的编码习惯所有输入数据统一用标幺值并且把基准容量写在文件头部注释里。节点的类型不要用1、2、3直接写死靠近数据输入的地方先定义常量比如SLACK1; PV2; PQ3;。每次迭代都用fprintf打印不平衡量的最大值哪怕最后注释掉调试阶段也不要省。对复数电压和电流统一用abs和angle提取数值不要在中间环节用实部虚部手工转换容易弄混。写函数时输入参数顺序最好固定为nbus, branch, bus, gen, load这样后续扩展成通用函数时不用改接口。我个人在实际调试中还有一个小技巧把节点导纳矩阵和短路分析用的序网阻抗矩阵打印成稀疏模式人为检查几行。比如Y矩阵的第3行非零元素应该对应和节点3相连的节点。如果多出一些不该有的非零项多半是线路数据里有重复支路或编号错误。这类问题在计算结果“看起来差不多”的时候特别难发现早点检查能省很多时间。这套流程从潮流计算到不对称短路分析本质上是把电力系统分析教材里的两大块知识用Matlab串了起来。你如果能把这段代码跑通再换成自己的系统数据基本上就可以应对课程设计或者工程中的初步分析了。后续还可以往里面加负荷模型、发电机详细模型、距离保护整定等内容但核心框架就是这些。