做电力系统分析时，最常碰到的两个任务就是电力系统潮流计算和不对称短路分析，而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] % ratio=0 表示普通线路；非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)) * (1+h); 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-1+j) = ([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相电压为0，B、C相电流为0 | 正序、负序、零序三个序网串联 |
| 两相短路（B、C相） | B、C相电压相等，A相电流为0 | 正序网与负序网并联，零序网不参与 |
| 两相接地短路（B、C相接地） | B、C相电压为0，A相电流为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_gen | Q_gen | P_load | Q_load |
|---|---|---|---|---|---|---|
| 1 | Slack | 1.06 | - | - | 0 | 0 |
| 2 | PV | 1.00 | 0.40 | - | 0 | 0 |
| 3 | PQ | 1.00 | 0 | 0 | 0.60 | 0.30 |
| 4 | PQ | 1.00 | 0 | 0 | 0.80 | 0.40 |
| 5 | PQ | 1.00 | 0 | 0 | 0.50 | 0.25 |
线路参数如下：
| 起点 | 终点 | R | X | 半电纳 |
|---|---|---|---|---|
| 1 | 2 | 0.020 | 0.060 | 0.030 |
| 1 | 3 | 0.050 | 0.200 | 0.020 |
| 2 | 3 | 0.040 | 0.150 | 0.020 |
| 2 | 4 | 0.060 | 0.250 | 0.020 |
| 3 | 4 | 0.080 | 0.300 | 0.020 |
| 4 | 5 | 0.100 | 0.350 | 0.020 |
4.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: 1=Slack, 2=PV, 3=PQ 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 结果解读
我按这个参数跑出来的结果如下，潮流收敛后的节点电压为：
| 节点 | 电压幅值 | 相角 |
|---|---|---|
| 1 | 1.0600 | 0.00° |
| 2 | 1.0000 | -2.06° |
| 3 | 0.9820 | -4.03° |
| 4 | 0.9660 | -5.17° |
| 5 | 0.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相短路电流 | 0 |
| C相短路电流 | 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直接写死，靠近数据输入的地方先定义常量，比如
SLACK=1; PV=2; PQ=3;。 - 每次迭代都用
fprintf打印不平衡量的最大值，哪怕最后注释掉，调试阶段也不要省。 - 对复数电压和电流统一用
abs和angle提取数值，不要在中间环节用实部虚部手工转换，容易弄混。 - 写函数时，输入参数顺序最好固定为
nbus, branch, bus, gen, load，这样后续扩展成通用函数时不用改接口。
我个人在实际调试中还有一个小技巧：把节点导纳矩阵和短路分析用的序网阻抗矩阵打印成稀疏模式，人为检查几行。比如Y矩阵的第3行非零元素应该对应和节点3相连的节点。如果多出一些不该有的非零项，多半是线路数据里有重复支路或编号错误。这类问题在计算结果“看起来差不多”的时候特别难发现，早点检查能省很多时间。
这套流程从潮流计算到不对称短路分析，本质上是把电力系统分析教材里的两大块知识用Matlab串了起来。你如果能把这段代码跑通，再换成自己的系统数据，基本上就可以应对课程设计或者工程中的初步分析了。后续还可以往里面加负荷模型、发电机详细模型、距离保护整定等内容，但核心框架就是这些。
