简介面向电力系统稳态与暂态分析学习者的Matlab代码资料包聚焦潮流计算与不对称短路分析的编程实现适合具备电力系统基础与Matlab能力的高校研究生、科研人员及能源系统工程技术人员。压缩包共1个docx文件约79KB以文档形式组织分析思路、代码实现方案与案例索引便于按主题查阅和对照练习。目前已有109人学习下载可作为课程设计、论文复现与科研建模的入门参考。内容在潮流与短路计算之外还延伸至碳排放流计算、无功补偿、配电网重构、微电网优化、储能配置、经济调度及多目标优化等场景并引入遗传算法、粒子群优化、支持向量机与神经网络等建模方法。读者可据此理解算法流程与模型构建逻辑借助Simulink仿真实例完成搭建、验证与排错并将优化算法和仿真模型作为科研项目的二次开发基础模块。1. 潮流算得对不对短路电流会替你打脸做电力系统潮流计算和不对称短路分析很多人习惯分成两个文件夹来做一个跑稳态电压一个套对称分量法公式。真到把 Matlab 代码实现跑通、拿结果去写报告的时候才发现这两件事根本拆不开——不对称短路里那个故障点开路电压就是潮流算出来的节点电压潮流给出的相角还间接决定了短路瞬间非周期分量往哪个方向偏。这份资源把两套东西放在一起一套是基于节点导纳矩阵和牛顿-拉夫逊法的潮流求解脚本一套是用对称分量法加复合序网做单相接地、两相短路、两相接地短路的故障计算再配一个 Simulink 模型把稳态初值灌进去、在指定时刻触发故障看波形。适合已经学过电路和电力系统分析、但自己从零写代码总是卡在收敛不了短路电流大了一个数量级的人。下面按我拆包时的顺序讲先立住原理再落到每一行参数怎么填。2. 潮流计算从支路表到 Ybus再到牛顿迭代收敛2.1 支路表怎么变成节点导纳矩阵Matlab 里没有天然的电网对象所有潮流程序的起点都是两张表。节点表我一般按[编号 类型 Pd Qd Vm0 角度0]排类型用 1 表示平衡节点、2 表示 PV 节点、3 表示 PQ 节点支路表按[i j R X B tap]排全部是标幺值。这一步看起来枯燥但九成的不收敛都出在这里R、X 忘了除以基准阻抗或者变压器变比的方向搞反了。标幺化的基准要先定死。取Sb 100 MVA各电压等级的Ub取平均额定电压那么阻抗基准Zb Ub^2/Sb支路表里的 R、X 就得先除以 Zb 再填进去。变压器如果按非标准变比在 i 侧建模tap 直接填Ui/Uj的实际比值不要填成标幺后的相对值否则最后算出来的电压会整体偏一截。function Y build_ybus(branch, nb) % branch 每行: [i j R X B tap]全部为标幺值 % B 为线路总充电电纳两端各挂一半tap 为理想变比无变压器填 1 Y zeros(nb, nb); for k 1:size(branch,1) i branch(k,1); j branch(k,2); z branch(k,3) 1j*branch(k,4); y 1/z; % 串联导纳 bc 1j*branch(k,5)/2; % 对地充电电纳的一半 t branch(k,6); if t 0, t 1; end % 容错空值按无变压器处理 Y(i,i) Y(i,i) (y bc)/(t^2); % 变比侧对角要除以 t 的平方 Y(j,j) Y(j,j) (y bc); Y(i,j) Y(i,j) - y/conj(t); % 非对角元素与变比相位有关 Y(j,i) Y(j,i) - y/t; end这段代码的逻辑是标准的 π 型等值串联支路进对角两次、进非对角两次取负并联充电电纳只进对角。参数上要注意三处。第一bc用的是总电纳的一半因为线路两端各有一半对地电容如果你手上的数据本身就是单端电纳那就要去掉这个/2。第二变比t只对 i 侧做修正Y(i,i)除以t^2是因为理想变压器把电压折算后电流也要折算这一项漏掉含变压器的算例必然发散。第三conj(t)只有在变比带相位移相变压器时才和t不同普通变压器两者相等写conj是为了保留通用性。拼完 Ybus先做一次自检sum(Y,2)应该接近 0不含对地支路时Ybus 必须对称对角元的虚部应该是负的感性网络。这三条对不上后面迭代写得再漂亮也是白搭。2.2 极坐标牛顿-拉夫逊四个分块矩阵加一次索引切片潮流方程本质是求解一组非线性代数方程牛顿法的核心是每次迭代用雅可比矩阵把非线性问题线性化。手写的时候最容易出错的是雅可比矩阵的对角元——用通式代入 ij 会漏掉一项导致收敛变慢甚至振荡。我的做法是先用矩阵运算一次性算出所有元素的通式值再单独覆盖对角元代码短且不容易错。function [Vm, th, P, Q, it] nr_pf(bus, Y, tol, maxit) nb size(bus,1); typ bus(:,2); pq find(typ 3); % PQ 节点下标 pvpq find(typ ~ 1); % 除平衡节点外全部参与有功不平衡 Vm bus(:,5); th deg2rad(bus(:,6)); Psp bus(:,3); Qsp bus(:,4); G real(Y); B imag(Y); for it 1:maxit V Vm .* exp(1j*th); S V .* conj(Y*V); % 各节点注入复功率 mis [Psp(pvpq) - real(S(pvpq)); Qsp(pq) - imag(S(pq))]; if max(abs(mis)) tol, break; end dth th - th.; % θij 矩阵 Vij Vm * Vm.; % Vi*Vj 矩阵 H Vij .* (G.*sin(dth) - B.*cos(dth)); N Vm .* (G.*cos(dth) B.*sin(dth)); M -Vij .* (G.*cos(dth) B.*sin(dth)); L Vm .* (G.*sin(dth) - B.*cos(dth)); % 对角元单独覆盖通式在 ij 处不成立 H(1:nb1:end) -imag(S) - B(1:nb1:end).*Vm.^2; N(1:nb1:end) real(S)./Vm G(1:nb1:end).*Vm; M(1:nb1:end) real(S) - G(1:nb1:end).*Vm.^2; L(1:nb1:end) imag(S)./Vm - B(1:nb1:end).*Vm; J [H(pvpq,pvpq) N(pvpq,pq); M(pq,pvpq) L(pq,pq)]; dx J \ mis; th(pvpq) th(pvpq) dx(1:numel(pvpq)); Vm(pq) Vm(pq) dx(numel(pvpq)1:end); end P real(S); Q imag(S);代码的关键点在J的拼装行按有功不平衡在前、无功不平衡在后排列按相角修正在前、幅值修正在后排四个分块的大小必须是[npvpq×npvpq, npvpq×npq; npq×npvpq, npq×npq]。如果这里行列顺序搞反J \ mis依然能算出一个数来但物理意义全错表现为迭代十几次还在一团乱跳。参数上tol取1e-8到1e-10标幺比较合适再小就进入浮点噪声区maxit给 20 足够正常算例 4 到 6 次就收敛。初值用平启动所有 PQ 节点Vm1.0、th0在大多数算例上够用只有重载或者长距离辐射状网络才需要先用直流潮流估一个相角。另外Vm .* (...)这一行依赖 Matlab 的隐式扩展R2016b 之前的版本要改成diag(Vm) * (...)否则会报维度不匹配——这类报错在老代码迁移时很常见。2.3 PV 节点无功越限与三个校验指标发电机节点在迭代中要先按 PV 处理每轮算完检查它的无功出力有没有超出上下限。超了就把这个节点临时转成 PQ把无功钉在限值上继续迭代收敛后再转回来核对一遍。这个逻辑不加含无功限值的算例会出现电压算出 1.15 甚至更高的假收敛结果看起来收敛了实际不可行。Qg imag(S) bus(:,4); % 发电机无功 注入无功 负荷无功 over typ2 (Qg bus(:,7) | Qg bus(:,8)); if any(over) bus(over,2) 3; % 转 PQ Qsp(over) bus(over,7); % 钉在上限或下限 typ bus(:,2); pq find(typ3); pvpq find(typ~1); end收敛之后不要急着出图先过三个校验。第一个是残差max(abs(mis))要小于设定容差这个程序自己会保证。第二个是功率闭合把所有节点的注入复功率加起来实部应该等于总网损虚部接近 0因为平衡节点注入的功率要正好覆盖负荷加损耗。第三个是迭代次数牛顿法是二次收敛残差应该呈每轮平方下降的形态5 节点算例 4 次左右、14 节点算例 5 次左右。如果迭代次数到了十几次还在慢慢降说明雅可比有错或者数据里有异常小阻抗。V Vm .* exp(1j*th); Sb 100; % 基准容量 MVA Sloss sum(V .* conj(Y*V)) * Sb; % 全网损耗含对地充电 Pgen real(V(bus(:,2)1) .* conj(Y(bus(:,2)1,:)*V)) * Sb; fprintf(网损 %.3f j%.3f MVA平衡节点注入 %.3f MW\n, ... real(Sloss), imag(Sloss), Pgen);网损这一项特别值得看一眼虚部。纯有功负荷、没有充电电纳的网络网损虚部应该是很小的负数无功损耗如果虚部大得离谱八成是某条支路的 B 填成了正数以外的量级或者标幺基准在某个电压等级上没统一。这套校验习惯我从第一次做配网重构算例时就养成了后面接手的每个潮流脚本都强制跑一遍省下来的返工时间比写代码本身多得多。3. 不对称短路三序网拼装与复合序网的 Matlab 落地3.1 正序、负序、零序三套阻抗不是同一个东西对称分量法的前提是把一组不对称三相量分解成正序、负序、零序三组对称分量每组分量各自服从独立的三相电路规律。落到计算上就是画三张序网图每张图的阻抗参数都不一样。正序网就是潮流用的那张网电源用故障点的开路电压等效各元件用正序阻抗。负序网没有电源电动势所有电源的负序电势为零只有负序阻抗通路同步机的负序阻抗大约取(XdXq)/2工程上常简化成正序的 0.9 到 1.0 倍线路的负序阻抗等于正序阻抗。零序网最麻烦它取决于有没有零序电流通路变压器绕组接法决定了零序能不能穿过线路的零序阻抗通常是有架空地线时正序的 2 到 3 倍零序电容也比正序大不少。变压器接法零序能否从该侧流入说明YNyn两侧都能中性点直接接地零序可穿越YNd仅 YN 侧能Δ 侧把零序短接在内部不出口YNy仅 YN 侧能另一侧无中性点引出Yy / Dd都不能无接地中性点零序网断开D 侧不能角接绕组不对外提供零序通路这张表是零序网搭建的全部依据。我见过最常见的翻车是把 YNd 变压器当成两侧都通处理结果零序阻抗直接少了一大块单相接地电流算出来偏大三四成还找不到原因。3.2 三类故障的边界条件与复合序网连接方式不对称短路的求解套路是固定的写出故障点的边界条件把它翻译成序分量之间的关系再把三张序网按这个关系连起来得到一个可以手算的等效电路。三类常见故障的结论如下。故障类型边界条件序网连接正序电流单相接地A 相IbIc0Va0三序网串联I1 Vf/(Z1Z2Z0)两相短路B、C 相Ia0Ib-Ic正序与负序并联零序断开I1 Vf/(Z1Z2)两相接地B、C 相接地Ia0VbVc0三序网并联I1 Vf/(Z1Z2·Z0/(Z2Z0))单相接地时三序电流相等故障相电流是正序电流的三倍这也是为什么中性点直接接地系统的单相接地电流往往比三相短路还大。两相短路没有零序通路故障相电流是正序电流的 √3 倍。两相接地则是三者中最复杂的零序和负序阻抗形成并联分流接地电流的大小取决于Z0/Z2的比值——这个比值在电缆网络和架空线网络里差别很大同一个算例换个线路模型结果能差一倍。这里有个经常被忽略的点Vf是故障发生前故障点的电压不是额定电压 1.0。正确做法是从潮流结果里取该节点的实际电压相量。如果偷懒直接用 1.0∠0°在重载或者长线路末端短路电流会被低估 5% 到 15%而短路电流计算本来就是按偏保守取值低估方向出错比算不准更麻烦。3.3 从序分量到三相电流代码与相序变换把上面的公式写成函数输入是故障点电压和三个序阻抗输出是三相电流相量。注意正序阻抗要用戴维南等效阻抗也就是把故障点看进去的等效阻抗做法是在潮流 Ybus 的基础上补上发电机次暂态电抗作为对地支路再求逆取对角元。function [Iabc, I012] asym_fault(Vf, Z1, Z2, Z0, ftype) % Vf 故障前故障点电压相量Z1/Z2/Z0 为故障点看进去的三序等效阻抗 a exp(1j*2*pi/3); switch upper(ftype) case SLG % 单相接地三序串联 I1 Vf/(Z1 Z2 Z0); I2 I1; I0 I1; case LL % 两相短路正负序并联 I1 Vf/(Z1 Z2); I2 -I1; I0 0; case LLG % 两相接地三序并联 I1 Vf/(Z1 Z2*Z0/(Z2Z0)); I2 -I1*Z0/(Z2Z0); I0 -I1*Z2/(Z2Z0); otherwise error(ftype 只支持 SLG / LL / LLG); end I012 [I0; I1; I2]; A [1 1 1; 1 a^2 a; 1 a a^2]; % 对称分量反变换矩阵 Iabc A * I012; % 得到 A/B/C 三相电流相量 end这段代码里A矩阵是标准的对称分量反变换第一行对应 A 相三序分量同相相加第二行对应 B 相正序乘 a²、负序乘 a第三行对应 C 相。之所以写成矩阵乘法而不是三行显式相加是为了让相序关系和变量名对得上改故障类型时不用动这块。参数上有两个坑要交代。第一Vf是相量带角度不能只传幅值否则算出来的三相电流相位关系是错的接到 Simulink 里看波形会发现三相不对称的方向反了。第二Z1必须是包含发电机次暂态电抗的戴维南阻抗如果你直接把潮流 Ybus 求逆取对角元得到的是不含电源内阻抗的网络阻抗数值会明显偏小。补发电机支路的写法很短Y1 Y; for k 1:numel(gen_bus) Y1(gen_bus(k),gen_bus(k)) Y1(gen_bus(k),gen_bus(k)) 1/(1j*xdpp(k)); end Z1 inv(Y1); Z1f Z1(fault_bus, fault_bus); % 故障点正序戴维南阻抗负序网把xdpp换成负序电抗再走一遍零序网则要把没有零序通路的支路整条剔掉只保留按 3.1 节那张表判断为可通的路径线路参数也换成零序参数。三张网各自求逆取故障点自阻抗就是公式里的Z1、Z2、Z0。最后提醒一点量纲标幺短路电流要乘基准电流才是有名值Ib Sb/(√3·Ub)其中 Ub 用故障点所在电压等级的基准电压。计算报告里写的短路电流 12.5 kA和程序里输出的 3.2标幺之间差的就是这个系数忘了乘就会得出一个看起来很安全、实际很危险的结论。4. Simulink 联调把潮流初值灌进模型再触发故障4.1 初值传递的三条路径纯 m 文件算出来的是相量稳态解要看暂态波形必须上 Simulink。麻烦在于电机和负荷模块如果从零状态启动前几个周期全是启动暂态根本分不清哪段是故障响应。解决办法是把潮流结果作为初始条件送进模型。第一条路径是assignin(base, ...)把Vm、th写进基础工作区模型里的电压源模块直接引用变量名。第二条路径是在模型的InitFcn回调里调用潮流脚本每次仿真自动重算一遍改参数不用手动同步。第三条是利用 powergui 自带的潮流工具在图形界面上对电机做初始化适合模型里有同步机、需要设置转子初始角的场景。% run_case.m先算潮流再把初值送进 Simulink 模型 [Vm, th, ~, ~] nr_pf(bus, Y, 1e-8, 20); assignin(base, Vm0, Vm); assignin(base, th0, th); assignin(base, Vf_fault, Vm(fault_bus) * exp(1j*th(fault_bus))); set_param(grid_fault_demo, StopTime, 0.25); set_param(grid_fault_demo, SimulationCommand, start);这段脚本的顺序不能反必须先算完潮流再assignin最后才start否则模型读到的是上一轮的旧变量。Vf_fault单独存一份是为了和 3.3 节的解析解对齐仿真跑完直接拿它做对照。仿真时长给 0.25 s 对应 50 Hz 系统的 12.5 个周期足够覆盖故障前稳态、故障暂态和故障后恢复三段。4.2 步长、故障时刻与故障电阻的取值Simulink 里参数填错最典型的后果是波形全是毛刺或者干脆发散。下面这张表是我自己调模型时总结的一组稳定组合算例规模不一样可以微调但数量级不要乱动。参数建议取值影响求解器ode23tb 或 ode15s刚性系统用变步长纯电磁暂态也可用离散化最大步长1e-5 ~ 5e-5 s太大丢高频分量太小仿真慢到跑不动故障起始时刻0.10 s对应 5 个工频周期暂态已衰减完故障持续时间0.05 ~ 0.10 s至少覆盖 2 到 3 个周期才能读到稳态峰值故障电阻 Rf0.001 ~ 0.01 Ω直接填 0 容易造成数值奇异powergui 类型Discrete 或 Continuous与步长设置必须匹配故障时刻的选取有个细节值得说如果只关心短路电流幅值故障在电压峰值时刻触发非周期分量最小波形干净好读数如果要做设备校验反而要让故障发生在电压过零附近这时直流分量最大冲击电流最恶劣。两种时刻各跑一遍取最大值才是保守结果。4.3 波形与解析解的对照口径波形出来之后别只看形状要落到三个数字上比对。第一个是故障相电流的稳态有效值和 3.3 节abs(Iabc)的结果比偏差应该在 5% 以内超了先查模型里的阻抗参数是不是和 m 文件一致。第二个是冲击系数也就是峰值除以稳态有效值的 √2 倍高压网络一般落在 1.8 附近低压配电网络因为阻抗比不同会低一些到 1.5 至 1.6。第三个是三相不对称度单相接地时非故障相的电压会升高中性点出现位移升高幅度和零序阻抗与正序阻抗的比值直接相关。t out.Ia.time; ia out.Ia.signals.values; win t 0.15 t 0.20; % 故障稳态段 Irms rms(ia(win)); Ipeak max(abs(ia(win))); kappa Ipeak / (sqrt(2)*Irms); % 反推冲击系数 fprintf(稳态 %.2f A峰值 %.2f A冲击系数 %.2f\n, Irms, Ipeak, kappa);如果kappa算出来是 1.41说明窗口里根本没有非周期分量多半是故障时刻正好卡在电压峰值如果算出 2.5 以上通常是窗口选到了故障刚发生的第一个半波把暂态峰值也统计进去了。窗口的选取规则是故障起始时刻加上两个周期之后开始到故障结束前为止。5. 避坑与排查五个把结果带偏一个数量级的细节5.1 数据和建模层的问题现象一潮流迭代十几次还在小幅振荡残差降不下去。原因多半是网络里有异常小的支路阻抗比如某条支路的 X 填了 0.0001 标幺导致 Ybus 条件数极差或者是 PV 节点的电压设定值互相冲突两个相邻的 PV 节点设了不同的目标电压但中间只有一条短线。解决办法是先把所有支路阻抗打印出来排序看有没有比同电压等级正常值小两个数量级的再把 PV 节点逐个放开成 PQ 跑一遍看是哪个节点引起的冲突。现象二潮流收敛了但平衡节点注入功率和总负荷加网损对不上。这是符号约定翻车的典型。负荷在节点表里填的是正数程序里S V.*conj(Y*V)算的是注入功率那么不平衡量应该写成Psp - real(S)其中Psp对负荷节点是负值。如果节点表里负荷直接填了正数又没取负潮流会往反方向收最终算出的电压偏高、网损偏小看起来结果很漂亮实际全错。我的习惯是在读数据的那一步就把负荷取负后面所有公式统一按注入理解。现象三单相接地电流算出来比三相短路还小一个数量级。原因通常是零序网建错了把 YNd 变压器的 Δ 侧也当成了零序通路或者线路的零序阻抗直接沿用了正序值。解决方式是单独把Z0打印出来看中性点直接接地的高压网络Z0/Z1一般落在 1 到 3 之间如果算出 0.1 以下几乎可以肯定是零序通路多接了一条如果算出 10 以上反而是漏了接地支路。这个比值是最快的体检指标比逐条支路查快得多。5.2 求解和仿真层的问题现象四m 文件算出的短路电流和 Simulink 波形差了三成以上。先别怀疑算法八成是两边的阻抗参数不是同一套。m 文件里用的是标幺值加戴维南等效Simulink 模型里填的是有名值中间隔着基准电压和基准容量两道换算。还有一种情况是模型里变压器的接线组别设成了 Yy 而 m 文件里按 YNd 处理零序通路直接对不上。排查办法是把两边的Z1、Z2、Z0都换算到同一基准后打印出来对比三个数都对上了再看波形。现象五仿真波形出现密集毛刺或者跑到一半报步长过小。毛刺通常来自步长和开关器件的配合变步长求解器遇上理想开关会产生数值跳变把最大步长收到1e-5以下往往能压住。报步长过小则是模型里有代数环或者某个模块的采样时间和求解器不匹配检查故障模块的接触器电阻是不是设成了 0以及 powergui 的离散化步长和模型步长是否一致。这类问题属于玄学范畴最快的办法是先把模型简化成单个电源加一条线路跑通了再往上加元件。6. 进阶把算例做成可回归的校验脚本代码写到能跑通只是第一步真正省时间的是把它变成改完参数一按就跑、跑完自动告诉你哪里不对的回归脚本。我的做法是给每个算例目录配一个check_case.m固定做三件事跑潮流、跑三类短路、把结果和一份基准结果文件比对。function ok check_case(case_name, tol) % 读入算例 → 潮流 → 三类短路 → 与基准比对 [bus, branch] load_case(case_name); Y build_ybus(branch, size(bus,1)); [Vm, th] nr_pf(bus, Y, 1e-8, 20); res struct(Vm, Vm, th, th); for ft {SLG,LL,LLG} [Iabc, ~] asym_fault(Vm(f)*exp(1j*th(f)), Z1, Z2, Z0, ft{1}); res.(ft{1}) abs(Iabc); end ref load([ref_ case_name .mat]); ok all(abs([res.SLG; res.LL; res.LLG] - ... [ref.SLG; ref.LL; ref.LLG]) tol); end这里的tol我一般取1e-4标幺因为不同 Matlab 版本的求解器在浮点末位上会有微小差异卡到1e-10会频繁误报。基准文件ref_*.mat只在确认结果正确时生成一次之后每次改代码都拿它当后悔药哪一步改坏了立刻能定位到是潮流部分还是短路部分。更实用的一层是加个交叉验证同一个算例用三种口径各算一遍看它们的偏差是否在预期范围内。解析解复合序网公式、仿真波形峰值、短路容量法粗估三者的典型偏差如下表。如果某次改动之后三者突然从彼此吻合变成各自为政不用怀疑一定是哪个环节的单位或者基准被动了。验证口径相对解析解的典型偏差适用场景复合序网解析解基准校验程序的正确性Simulink 波形峰值2% ~ 8%校验模型参数与暂态过程短路容量法粗估10% ~ 20%初设阶段快速估算我踩过最深的一次坑是改了一版负序阻抗的默认系数之后忘了同步基准文件结果后面连着三个算例的报告数据全是错的一直到答辩前一天晚上对着波形本文还有配套的精品资源点击获取