前几天一个做变电站设计的师弟拿着短路电流计算结果来找我说Matlab跑出来的单相接地短路电流跟手册差了一大截。我看了半天代码问题出在他短路计算用的电源电势直接取了1.0压根没走潮流计算——故障点的实际电压根本不是标幺值1.0。这其实是把“电力系统潮流计算”和“不对称短路分析”两个环节割裂开来的典型错误。潮流计算和不对称短路分析是电力系统稳态分析里最核心的两块内容前者解决“系统正常运行时各节点电压、支路潮流是多少”后者解决“系统发生不对称故障时短路电流、故障相电压怎么算”。很多教材把这两个问题分开讲但工程实际里它们必须串联起来短路点的故障前电压恰恰来自潮流计算结果。用Matlab把这两块写到一套代码里既能讲清楚算法原理又能直接用于课程设计、毕业设计甚至小型工程校核。这篇文章我会从一个5节点测试系统出发带你从节点导纳矩阵建起写完牛顿-拉夫逊法潮流再叠加正序、负序、零序网络做不对称短路分析最后把两者串成完整的分析流程。1. 潮流计算到底在算什么先搞清分析对象与工程场景1.1 为什么要自己写潮流程序做电网规划、运行方式安排、保护定值整定、新能源并网分析第一步几乎都是算潮流。所谓潮流计算就是在给定网络拓扑、发电机有功出力和机端电压、负荷有功无功的条件下求解各个节点的电压幅值、相角以及各条支路的功率分布。算完之后你才知道哪条线路重载、哪个节点电压偏低、系统网损是多少。商业软件比如BPA、PSS/E、PSASP当然能直接出结果但如果你是学生或者正在做算法研究自己用Matlab写一遍的价值完全不同。第一你能真正理解牛顿-拉夫逊法的迭代过程、雅可比矩阵的物理含义而不是把软件当黑盒第二潮流程序写完后可以跟短路计算联合使用做故障前电压的精确计算第三后续加分布式电源、调压策略、优化算法所有逻辑都在自己手里自由度非常高。1.2 节点类型设计与一个可复现的5节点算例潮流计算一般把节点分成三类这是必须刻在脑子里的基本概念节点类型已知量待求量典型元件PQ节点P、QV、δ负荷节点、无调压能力的发电节点PV节点P、VQ、δ有自动电压调节器的发电机节点平衡节点SlackV、δP、Q承担系统功率差额的等值大系统平衡节点在数学上必不可少因为系统总有网损所有发电机的有功不能提前精确给定必须有一个节点来吸收或者补足不平衡功率。一般选择主网架中容量较大、调频能力强的节点。为了让你能直接复现我搭了一个5节点测试系统所有参数用标幺值给出基准容量取100MVA节点1平衡节点V1.05∠0°节点2PV节点P0.5puV1.02pu节点3PQ节点负荷0.8j0.25pu节点4PQ节点负荷0.4j0.15pu节点5PQ节点负荷0.3j0.1pu线路参数如下R、X、B/2均为标幺值支路首端末端RXB/2L1120.020.060.030L2130.050.200.030L3230.040.150.020L4240.060.250.030L5250.050.180.030L6340.030.120.020L7450.040.160.025这个系统的规模不大不小既能手工验算中间结果又能体现算法特征后面所有代码和结果都基于它展开。2. 节点导纳矩阵的Matlab建模后续一切计算的地基2.1 自导纳与互导纳的定义节点导纳矩阵Y矩阵是所有潮流和短路计算的基础。它的形成规则非常简单对角线元素Yii是连接到节点i的所有支路导纳之和包括接地支路非对角线元素Yij是节点i和节点j之间支路导纳的负值。用公式表示就是Yii Σj yij yi0其中yi0是节点i的对地导纳线路充电电容、变压器励磁导纳等Yij -yij这里的yij是支路导纳等于1/(RjX)如果是变压器支路还涉及变比折算。理解Y矩阵的关键在于它把整个网络的电气连接关系压缩成了一个矩阵矩阵的稀疏程度跟网络的联络紧密程度直接相关。2.2 带非标准变比变压器的处理如果你的系统里有变压器Y矩阵不能简单按线路处理。三相变压器用标幺值表示时如果变比不等于1:1即非标准变比工程上通常用π型等效电路来处理。设变压器支路标准变比为kk≠1短路阻抗为zt那么等值电路中的三个导纳参数为串联臂yt/k左侧并联臂(1-k)yt/k²右侧并联臂(k-1)yt/k实际程序里为了简化处理很多教材和代码直接令k1把变压器当成普通阻抗支路处理。但如果你的课题涉及多电压等级网络非标准变比必须实现否则算出来的电压分布会跟实际相差很远。变压器变比档位调节在潮流计算中本身就是一种调压手段这一点在做毕业设计时尤其容易忽略。2.3 用稀疏矩阵高效装配Y矩阵的代码模板Matlab中直接用一个双重循环组装Y矩阵在节点数少的时候看不出问题但到了几百上千节点的大电网效率就很低。工程上推荐用稀疏矩阵的思路先记录非零元素的位置和值最后用sparse一次性构建。下面是我常用的装配模板% 输入数据格式bus节点表、branch支路表 % branch列[from, to, R, X, 0.5*B, k] % k为变压器非标准变比若为线路则k1 nb size(bus, 1); nl size(branch, 1); Y zeros(nb, nb); % 先用普通矩阵装配小系统方便调试 for i 1:nl f branch(i, 1); t branch(i, 2); R branch(i, 3); X branch(i, 4); B branch(i, 5); k branch(i, 6); z R 1j * X; y 1 / z; if k 1 % 线路 Y(f, f) Y(f, f) y 1j * B; Y(t, t) Y(t, t) y 1j * B; Y(f, t) Y(f, t) - y; Y(t, f) Y(t, f) - y; else % 变压器非标准变比 Y(f, f) Y(f, f) y / k^2; Y(t, t) Y(t, t) y; Y(f, t) Y(f, t) - y / k; Y(t, f) Y(t, f) - y / k; % 变压器π型等效注意内部节点处理这里是把变比折算到f侧 end end % 大系统可以改用稀疏模式 % I []; J []; S []; % 收集 (i, j, value) 三元组后用 sparse(I, J, S, nb, nb)这里有一个非常关键的细节并联电容B的单位。很多教材给的是线路总的充电电纳B而组装时每端只加B/2所以输入数据里我习惯直接给B/2代码里不再除以2。如果你从PSASP或BPA导出数据要注意原始数据到底是总电纳还是半电纳这个搞错了潮流结果会直接出现莫名其妙的无功分布。Y矩阵建好之后建议做一个自检把Y矩阵打印出来检查每行非对角线元素之和再加对角线元素如果网络没有接地支路理论应等于0有接地支路时这个值等于接地导纳之和。3. 牛顿-拉夫逊法潮流计算迭代到底在迭代什么3.1 极坐标形式的功率不平衡方程牛顿-拉夫逊法是目前潮流计算的主流算法其核心思想是先给一组电压初始值然后根据节点功率平衡条件求解修正量不断迭代直到功率不平衡量小到可以接受。用极坐标表示节点电压Vi Vi∠δi那么节点注入功率方程为节点i的有功和无功可以写成Pi Vi Σj Vj (Gij cosδij Bij sinδij)Qi Vi Σj Vj (Gij sinδij - Bij cosδij)其中δij δi - δj。注意我用了Vi、δiGij、Bij是Y矩阵的实部虚部。所谓功率不平衡量就是ΔPi Pi_spec - Pi_calcΔQi Qi_spec - Qi_calcPi_spec是给定的注入功率发电出力减负荷Pi_calc是用当前电压初值代入上式计算出的注入功率。我们的目标就是通过迭代把ΔP、ΔQ压到零。3.2 雅可比矩阵的解析与修正方程把功率不平衡方程对电压幅值和相角求偏导就得到雅可比矩阵J。极坐标下J可以分成四个块子块表达式物理含义H∂P/∂δ有功对相角的敏感度NV·∂P/∂V有功对电压幅值的敏感度K∂Q/∂δ无功对相角的敏感度LV·∂Q/∂V无功对电压幅值的敏感度修正方程的形式是[ ΔP ] [ H N ] [ Δδ ] [ ] -[ ] * [ ] [ ΔQ ] [ K L ] [ ΔV/V ]注意右边修正量用的是ΔV/V而不是直接的ΔV这是极坐标形式的特点。V·∂P/∂V这种形式的偏导数在推导时更简洁迭代中解出ΔV/V之后实际电压修正为V_new V ΔV。这算是一个容易让人困惑的细节很多新手在这里对着公式莫名其妙。雅可比矩阵中各元素的偏导公式比较长但Matlab实现时可以直接用数值方法验证。工程上为了稳妥我一般先用解析表达式计算再用简单的有限差分核对一遍子矩阵元素确认无误后才继续往下写。这属于排错成本最低的手段。3.3 完整迭代主循环与算例结果迭代主循环的逻辑如下初始化所有PQ节点电压幅值取1.0相角取0平启动PV节点电压取给定值计算ΔP、ΔQPV节点只算ΔP平衡节点两者都不算判断max(|ΔP|, |ΔQ|)是否小于阈值比如1e-6满足则结束组装雅可比矩阵求解修正方程得到Δδ和ΔV/V更新电压δ δ ΔδV V ΔV返回第2步核心代码片段V bus(:, 3); % 电压幅值初值 theta bus(:, 4); % 相角初值 tol 1e-6; iter 0; maxiter 20; PQ find(bus(:, 2) 1); % PQ节点集合 PV find(bus(:, 2) 2); % PV节点集合 ref find(bus(:, 2) 3); % 平衡节点 while iter maxiter % 计算P、Q Pcal zeros(nb, 1); Qcal zeros(nb, 1); for i 1:nb for k 1:nb Pcal(i) Pcal(i) V(i)*V(k)*(G(i,k)*cos(theta(i)-theta(k)) B(i,k)*sin(theta(i)-theta(k))); Qcal(i) Qcal(i) V(i)*V(k)*(G(i,k)*sin(theta(i)-theta(k)) - B(i,k)*cos(theta(i)-theta(k))); end end dP Pspec - Pcal; dQ Qspec - Qcal; % 划去平衡节点PV节点只留dP dP(ref) []; dQ(PV) []; if max(abs([dP; dQ])) tol break; end % 组装雅可比矩阵并求解 J buildJacobian(V, theta, G, B, PQ, PV, ref); dX J \ [dP; dQ]; nP length(PV) length(PQ); dTheta dX(1:nP); dVnorm dX(nP1:end); % 更新 theta(PQ) theta(PQ) dTheta(1:length(PQ)); theta(PV) theta(PV) dTheta(length(PQ)1:end); V(PQ) V(PQ) .* (1 dVnorm); iter iter 1; end这个代码里buildJacobian是需要专门编写的函数篇幅有限不展开全部元素但有一点必须提醒雅可比矩阵每次迭代都要重新计算因为它依赖当前电压幅值和相角。有些初学者为了省事把它当成常数矩阵迭代几步就开始发散这是最常见的错误。用上面5节点系统跑完收敛后得到的结果大致如下节点电压幅值/pu相角/°11.0500.0021.020-1.2330.982-3.4540.975-4.0250.988-3.28这个结果符合物理直觉离平衡节点越近的节点电压越高负荷重的节点3、4电压偏低系统整体没有越限。潮流计算这部分到这里就闭环了。4. 不对称短路的序网络法三序分量怎么拆怎么合4.1 对称分量法与序阻抗三相短路是对称故障用单相电路就能解决。但实际电网中单相接地、两相短路、两相接地是更常见的不对称故障这时候三相电路互有耦合不能直接拆成单相分析。对称分量法的核心思想是把一组不对称的三相电气量分解成三组对称的分量——正序、负序、零序然后分别在三序网络中计算最后叠加回三相量。以电流为例Ia Ia1 Ia2 Ia0Ib a²·Ia1 a·Ia2 Ia0Ic a·Ia1 a²·Ia2 Ia0其中a e^(j120°) -1/2 j√3/2a² e^(j240°)。正序分量是幅值相等、相位依次超前120°的对称量负序分量是相位依次滞后120°的对称量零序分量是三相同相位的量。分解用的变换矩阵为[ Ia0 ] [ 1 1 1 ] [ Ia ] [ Ia1 ] [ 1 a a² ] [ Ib ] [ Ia2 ] [ 1 a² a ] [ Ic ]三序网络本质上是三个独立的单相网络正序网络与正常运行时的网络完全相同包含所有电源负序网络不含电源元件的负序阻抗与正序略有差异零序网络最难处理因为零序电流的中性点通路决定了它的拓扑结构。4.2 三类序网络的特性和零序通路注意事项三序网络的区别用一张表说明项目正序网络负序网络零序网络电源有发电机电动势无无与中性点接地关系无关无关强相关元件序阻抗正序阻抗z1负序阻抗z2≈z1零序阻抗z02~4倍z1变压器接线影响经过变压器正常传递正常传递只有YN侧能提供零序通路零序网络是初学者最容易失分的地方。零序电流只能在有接地中性点的绕组中流通比如YNd接线的变压器YN侧发生接地故障时零序电流可以流过变压器而三角形绕组侧由于没有中性点接地零序电流在该侧形不成通路。这会导致零序网络在某些节点处直接断开了所以在构建零序网络时必须根据变压器的接线方式决定是否把变压器支路纳入并且在短路计算之前先检查整个网络中是否存在完整的零序通路。如果故障点所在的系统零序阻抗无穷大那单相接地短路电流反而可能比三相短路电流还小这是完全反直觉的但物理上完全合理。4.3 边界条件与复合序网不同类型的短路对应不同的故障点边界条件而边界条件决定了三个序网络怎么连接。把故障端电压、电流的序分量关系整理出来然后用序网络的串并联表达就得到复合序网。故障类型边界条件特征复合序网连接方式短路电流公式三相短路三相对称只有正序网络If E / Z1单相接地短路Ia≠0IbIc0三序网络串联If 3E / (Z1Z2Z0)两相短路Ib-IcIa0正序、负序并联If √3·E / (Z1Z2)两相接地短路Ia0IbIcIg三序网络先并后串If 3E·Z2 / [Z1(Z2Z0)Z2Z0]这些公式看起来是结论但推导过程才体现理解程度。比如单相接地边界条件意味着三序电流相等Ia1Ia2Ia0E/(Z1Z2Z0)所以故障相电流是三者之和的3倍。两相短路时正序电流等于负序电流的相反数故障相的短路电流幅值是正序电流的√3倍。4.4 Matlab中按序阻抗矩阵求解短路电流在计算机程序中一般不用手工推导的戴维南等值阻抗而是先生成三序网络的节点导纳矩阵Y1、Y2、Y0再求逆得到节点阻抗矩阵Z1、Z2、Z0直接取故障点k的对角元素Zkk就得到了从故障点看进去的戴维南等值阻抗。这个方法的优势在于不需要做网络化简程序通用性极强。% 假设Y1、Y2、Y0已经构造完成 Z1f full(inv(Y1)); % 故障点k为第k个对角线元素 Z2f full(inv(Y2)); Z0f full(inv(Y0)); % 故障点编号 kf z1kk Z1f(kf, kf); z2kk Z2f(kf, kf); z0kk Z0f(kf, kf); Ef 1.0; % 故障前电压后续会替换为潮流结果 switch faultType case 3ph If Ef / z1kk; I1 If; I2 0; I0 0; case singleLG I1 Ef / (z1kk z2kk z0kk); I0 I1; I2 I1; If 3 * I1; case twoPhase I1 Ef / (z1kk z2kk); I2 -I1; I0 0; If sqrt(3) * abs(I1); case twoPhaseGnd I1 Ef / (z1kk z2kk * z0kk / (z2kk z0kk)); I2 -I1 * z0kk / (z2kk z0kk); I0 -I1 * z2kk / (z2kk z0kk); If 3 * abs(I0); end这个框架可以覆盖所有常规短路类型。注意短路电流求出来后是标幺值要换算成有名值必须乘以基准电流Ib Sbase / (√3·Ubase)。我在实际项目中见过有人忘了这一步把0.05pu的电流直接当成500A上报结果整定计算差了一个数量级。5. 从潮流到短路的联合分析故障前电压的正确用法5.1 为什么不能把电源电势直接取1.0教材里讲不对称短路计算时为了方便推导通常默认故障前系统空载即故障点电压为额定电压并且所有电源电动势都取1.0∠0°。但实际系统带负荷运行时由于线路压降和变压器阻抗的影响故障点电压往往偏离1.0。前面5节点潮流的结果已经显示节点4的电压只有0.975pu节点3才0.982pu。短路前电压取1.0还是取0.975对短路电流的计算结果影响可达3%~8%对保护整定来说这个误差绝对不能忽略。从物理本质上讲故障前系统内各节点的电压是全网功率分布的结果发电机节点电压可以靠自动电压调节器维持在1.0以上但负荷节点电压就是会跌。短路计算本质上是在稳态运行点的基础上叠加一个故障扰动故障点的戴维南等值电动势就应该取故障前该点的实际电压。把潮流和短路分成两套独立程序、互不传递数据是我见过最多的工程计算错误。5.2 联合计算的完整程序流程把两部分集成到一起后完整的流程是输入网络拓扑、线路/变压器参数、发电机出力、负荷数据构造节点导纳矩阵Y1正序网络的基础运行牛顿-拉夫逊法潮流计算得到所有节点电压V和相角θ记录故障点kf的故障前电压Uf V(kf)·e^(jθ(kf))由正序网络参数构造Y1由负序网络参数构造Y2通常与Y1结构相同但阻抗取z2由零序网络参数构造Y0根据变压器接线方式调整拓扑对Y1、Y2、Y0求逆取故障点对角元得到z1kk、z2kk、z0kk代入复合序网公式求短路电流序分量通过对称分量反变换得到三相短路电流和各节点故障电压这个流程里第4步是潮流与短路的接口也是最容易被忽略的一步。5.3 算例对比三相短路与单相接地的差异在5节点系统里假设节点4发生三相短路和单相接地短路分别用故障前电压Uf潮流结果替代1.0观察差异。先看故障前节点4电压Uf 0.975∠-4.02°。三相短路时用Uf0.975计算If 0.975 / |z1kk|假设z1kk ≈ j0.15pu则If ≈ 6.5pu用Uf1.0计算If ≈ 6.67pu误差约2.6%单相接地时短路电流公式为If 3Uf / |z1z2z0|故障电流会汇集到零序通路数值上可能达到三相短路的1.5~2.5倍跟系统接地方式关系很大更重要的是短路后的节点电压计算需要用到故障前电压作为初值。用潮流结果计算故障后节点电压分布更符合实际节点电压可能不是单纯跌落或者抬升还会伴随相角偏移。如果你的短路计算结果用来校验母线电压、保护动作行为这些细节直接决定结果是否可信。我在这里强烈建议你把潮流结果、短路电流、短路后节点电压这三个输出画在同一张图上对比。Matlab里用plot或者bar画节点电压幅值曲线能看到故障前后电压分布的明显差异这也是毕设和论文里很有说服力的图形素材。6. 工程实操中的常见坑与收敛性调优经验6.1 初值与平启动牛顿-拉夫逊法对初值敏感。默认的平启动所有PQ节点V1.0δ0在大多数输电网中都能收敛但在配电网这类高R/X比的网络里平启动经常发散或者收敛到奇怪的解。解决思路有两个第一是改用保留非线性项的潮流算法或者PQ分解法第二是先用直流潮流算一遍得到相角初值再带入牛顿-拉夫逊法迭代。在写程序的时候建议把初值作为可配置参数暴露出来方便调试不同网络。6.2 无功越限和PV/PQ节点动态转换PV节点之所以能维持电压幅值是因为发电机具备无功调节能力。但发电机的无功出力有上下限当迭代过程中某台发电机的无功Q越出上限或下限时实际运行中这台机组就丧失了电压调节能力节点应该从PV节点转换为PQ节点即不再固定V而是固定Q在越限边界值。这个转换逻辑在潮流主循环里要用一个标志位记录每次迭代都检查一旦触发就一直保持到计算结束。很多教材代码没有处理这一条算出来的PV节点无功超限也不报错这在工程上是不能接受的。6.3 高阻网络与雅可比矩阵病态配电网、微电网、部分新能源汇集系统的R/X比很高雅可比矩阵条件数变大导致修正量振荡。实操中我常用的手段是加一个阻尼因子α把修正方程改成ΔX_new α · J \ ΔFα取值范围0.5~0.9牺牲一点收敛速度换来稳定性。另一个做法是换用直角坐标形式或者改用带最优乘子的牛顿法。这些方法不复杂但能在关键时刻救回一个发散的case。6.4 标幺值、基准值和零序通路的细节最后汇总一下我踩过的坑基准值不一致潮流计算中所有数据必须统一到同一套基准容量和基准电压下。线路参数如果是从有名值换算来的换算系数要先核对否则短路电流结果错得离谱还很难发现。零序网络构建不要直接拿正序Y矩阵的拓扑当零序网络用。必须逐个检查变压器支路的中性点接地方式、变压器的接线组别以及负荷的中性点接地情况。漏了一条零序通路计算结果可能差出几倍。短路电流的数值计算得到的短路电流是稳态分量工程上还要考虑非周期分量和冲击系数来求冲击电流用于开关设备开断能力校验。把稳态短路电流直接当冲击电流用会出大问题。矩阵求逆方式节点多的时候不要用inv(Y)求Z矩阵用Y \ eye(n)或者先做LU分解再对指定列求解效率高一个数量级。坦白说我自己第一次把潮流和不对称短路写成联合程序时也经历过从“算出数就高兴”到“发现数值看起来对但物理上不对”的过程。根源就是没有串起故障前电压、没有仔细检查零序通路。这些细节光看书本很难意识到只有在跑实际数据、对着参考手册反复比对的时候才能体会。如果现在的你正在做类似的Matlab实现我的建议是先把潮流计算跑稳把节点电压和支路潮流跟教材或者商业软件对上再去做短路分析短路分析先做三相短路再做单相接地和两相短路一步一步扩展千万不要一上来就追求“一键出所有结果”。算法这东西一旦理解透了底层逻辑后面换什么网络、加什么故障类型都是水到渠成的事。