MATLAB下的3机9节点电力系统暂态稳定分析程序实现与验证

简介：针对电力系统暂态稳定分析需求，基于MATLAB的3机9节点系统暂态稳定计算程序完整实现了暂态稳定计算流程，适合电力专业学生、研究人员及工程师用于教学自学与工程验证。压缩包共30个文件，以18个m源文件为主，涵盖数据导入、初始值计算、潮流求解、雅可比矩阵、故障模拟及结果绘图等模块，另有9个asv自动备份文件、2个doc配套文档和1个txt数据文件，整体仅219KB。其中《暂态稳定分析程序报告》详解了发电机动态建模、网络方程和ode45求解思路，《数据格式说明》则帮助使用者按规范整理输入数据。程序采用模块化设计，通过主函数串联各子模块，便于修改参数和观察不同扰动条件下的机组功角、电压等动态响应；主函数与子函数、数据文件分离，也方便二次开发。资源已有302人学习下载，对理解暂态稳定机理和提升MATLAB电力系统编程能力都具实用价值。

1. 从3机9节点到暂态稳定：先理解这个算例在算什么

电力系统暂态稳定分析里，3机9节点系统出现频率可能比任何IEEE标准算例都高。规模不大——3台发电机、9条母线、3个负荷，但发电机动态、网络代数约束、故障与切除操作全部包含，且参数有公开参考值，结果可互相校验。很多研究生的第一个暂态稳定程序，就是在这个系统上跑通的。用MATLAB实现的好处很直接：矩阵运算是原生能力，数值积分可以少量代码手写，不需要编译，也不依赖第三方仿真平台。标题里的计算程序，本质上要做的事就是：给定系统参数和运行工况，人为设置一个扰动（通常是三相短路），然后逐时步求解描述发电机转子运动的微分方程，观察功角能否恢复同步。下面这套流程，按我平时的工程习惯，从模型建立、初始化、仿真求解到故障操作一次讲完，并把单位、初值、判据这几个最容易出错的地方单独点出来。

2. 暂态稳定模型：经典二阶模型与网络代数方程

2.1 发电机用经典二阶模型，为什么够用

暂态稳定关心的是扰动后1到10秒内发电机的转子运动，电磁暂态细节可以忽略，因此最常用的是经典二阶模型，也叫摇摆方程。它把每台发电机等效为一个暂态电抗后的恒定电压源，转子运动用一对微分方程描述：

• 功角对时间的变化率等于转子转速偏差
• 转速变化率取决于机械功率与电磁功率之差，除以惯性时间常数

机械功率一般在几秒内认为不变，电磁功率则由网络方程求出。这样整个系统被拆成一个微分方程块加一个代数方程块，合称微分代数方程组。求解的关键是，每个积分步里先解代数方程求出电磁功率，再推进微分方程。

function dx = swing_dynamics(t, x, Ybus_red, Pm, M) % 状态向量 x = [delta1 delta2 delta3; omega1 omega2 omega3] delta = x(1:3); omega = x(3:4); E = abs(Ep) * exp(1j * delta); % 暂态电动势相量 I = Ybus_red * E(:); % 消去网络后的注入电流 Pe = real(E(:) .* conj(I)); % 电磁功率 dx = [omega; (Pm - Pe) ./ M]; end

代码里的Ybus_red是消去负荷节点后的发电机内电势节点等值导纳矩阵。M是惯性时间常数换算后的转动惯量标幺值。实际工程中这一步通常还要加一个阻尼项，系数取0到2之间，用来模拟转子阻尼绕组和机械阻尼的作用，不加阻尼的系统在扰动后功角会持续摆动不衰减，这不影响稳定性判断，但影响曲线观感。

2.2 网络方程如何处理：把负荷变成恒阻抗

暂态稳定计算里网络方程是线性的，前提是把负荷处理成恒阻抗。这样负荷可以折算成接地导纳并入导纳矩阵，网络节点只剩下发电机内电势节点，最终得到一个节点数为发电机台数的低阶等值导纳矩阵。这一步做对了，后面每个积分步的代数求解就退化成一次矩阵乘法，计算量非常小。

常见做法是：先形成完整节点导纳矩阵，把负荷按初始电压和初始功率折算成阻抗，加到对应母线对角元上，再通过高斯消去法消去无源母线，只保留发电机内电势节点。消去后的矩阵随系统拓扑固定不变，故障期间和故障切除后只是局部修改然后重新消去。

2.3 九节点系统的典型参数与基值选择

3机9节点（常称为WSCC 9-bus）的标准参数组合，在IEEE和大多数教材中基本一致：基准容量100 MVA，基准电压230 kV，三台发电机分别接在母线1、2、3上，变压器变比和线路阻抗都有明确的标幺值。我一般直接把参数写成MATLAB脚本里的数组，集中管理，方便修改。

注意H的单位是秒，要和基值功率配合换算。转动惯量M = 2H/ωs，在50 Hz系统里ωs = 2π×50，标幺化后M的数量级通常在几十到几百之间，这个数量级对积分步长选择很敏感，后面第4章专门说。

3. MATLAB实现：导纳矩阵形成、潮流初值与程序骨架

3.1 数据组织方式与基值换算

程序第一步是组织数据。我的习惯是把线路、变压器、发电机参数分别写成嵌套结构体或表格，便于对照原始算例。线路参数统一用标幺值存储，全部按100 MVA基值折算。如果原始数据是实际值，折算公式是：

Z_pu = Z_actual × S_base / V_base²

这一步最容易出错的是线路充电电容和变压器变比。九节点系统里变压器接在发电机升压侧，变比不是1:1，处理导纳矩阵时必须先把变压器导纳折算到同一侧，再参与节点导纳组装。忽视变比直接代数值，初始化相位就会错。

3.2 形成节点导纳矩阵的核心函数

function Y = build_ybus(bus, branch) % 输入: bus节点表, branch支路表[首端 末端 电阻 电抗 半电纳/2] n = size(bus, 1); Y = zeros(n, n); for k = 1:size(branch, 1) f = branch(k, 1); t = branch(k, 2); z = branch(k, 3) + 1j * branch(k, 4); y = 1 / z; b = branch(k, 5); 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; end % 负荷折算为恒阻抗并加到对角元 for i = 1:size(bus, 1) if bus(i, 4) ~= 0 % 有功负荷 Y(i, i) = Y(i, i) + conj(bus(i,4) + 1j * bus(i,5)) / (bus(i,6)^2); end end end

这段代码思路是：先按支路串联阻抗形成导纳，把半电纳加到两端节点对角元，然后把已知运行电压下的负荷功率折算成阻抗。这里功率和电压必须用复数共轭相除，因为负荷是“从节点吸收功率”，电流方向与注入相反。很多计算对不上，问题都出在这个共轭上。

3.3 初值计算：必须从潮流结果出发

暂态稳定不是从零开始算，而是从稳态工况出发。初值包括每台发电机的暂态电动势幅值和初相角。求法是在潮流计算结果基础上，对每台发电机做一步戴维南等效：

E' = Vt + jXd' × I_gen

其中Vt是机端电压，I_gen是发电机注入电流。机端电压和注入功率是潮流给出的结果，所以整个程序的前置步骤必须有一个潮流计算，哪怕是简化版的高斯-赛德尔法也行。我一般在做暂态稳定之前先跑一遍牛顿-拉夫逊潮流，把各母线电压的幅值和相角存下来，作为初始化输入。

% 由潮流结果求发电机内电势 for g = 1:ng Vt = Vbus(gen_bus(g)); S = Pgen(g) + 1j * Qgen(g); % 发电机注入复功率 Igen = conj(S / Vt); % 注意电流方向 Ep(g) = Vt + 1j * Xd(g) * Igen; % 暂态电动势 end

这里要特别注意功率方向约定：潮流计算里发电机节点通常是PQ或PV节点，正方向是注入网络，所以电流等于复功率共轭除以电压。如果这一行写反，初始功角会整体偏移，接下来仿真很难收敛。

3.4 主程序框架按拓扑变化分段

故障仿真的特点是系统拓扑随时间变化：故障前、故障中、故障后三个时间段的导纳矩阵不同。常见的做法是提前把三个矩阵都算好，仿真循环里只做一个切换判断，而不是每个积分步重新组装矩阵。

% 主时间推进循环 t = 0; x = x0; h = 0.005; T_end = 5; while t < T_end if t < t_fault Yr = Y_pre; elseif t < t_clear Yr = Y_fault; else Yr = Y_post; end x = euler_mod(t, x, h, Yr, Pm, M); t = t + h; % 记录功角, 判断是否失稳 end

故障期间导纳矩阵的形成方式取决于故障类型。三相短路通常模拟为故障点对地接入一个极小阻抗，这样故障点电压几乎为零，发电机输出的电磁功率骤降，导致转子加速。切除故障则相当于把故障点和相关支路同时从网络中移除，系统的电气距离变大，传输能力下降。

4. 暂态仿真求解：改进欧拉法与暂稳判据

4.1 为什么不用ode45而用改进欧拉

MATLAB自带的ode45在普通微分方程上表现很好，但暂态稳定仿真并不推荐直接用它。原因有两点：微分代数方程组在切换时刻存在非光滑跳变，变步长求解器容易因为步长振荡降低效率，甚至导致事件检测误差；另外，暂态稳定程序经常要和后续的优化、批量扫描配合，固定步长便于代码结构统一。改进欧拉法，也就是二阶龙格-库塔法，精度足够，实现简单，是电力系统暂态稳定计算程序里最常见的选择。

function x_new = euler_mod(t, x, h, Yr, Pm, M) % 预测步 kx1 = swing_dynamics(t, x, Yr, Pm, M); x_pred = x + h * kx1; % 校正步 kx2 = swing_dynamics(t + h, x_pred, Yr, Pm, M); x_new = x + h * (kx1 + kx2) / 2; end

两个关键参数需要关注：步长h和总仿真时长。步长通常取0.005秒到0.01秒。步长太大会使功角曲线出现虚假发散或阻尼，太小则计算量倍增。总时长一般取5秒左右，既能覆盖暂态过程又不至于算太久。

4.2 电磁功率计算细节

每个时步里电磁功率的计算，是整个程序的计算瓶颈，也是和网络方程交互的部分。计算过程是先由当前功角合成发电机内电势相量，再乘以等值导纳矩阵得到电流，最后取实部得到功率。这里有一个容易忽略的问题：导纳矩阵的维度和排序。发电机节点顺序必须与状态向量里功角的顺序完全一致，我建议在初始化阶段显示输出一次矩阵维度，避免后续索引错位。

电磁功率表达式也可以写成显式形式：Pe_i = Ei² × Gii + ΣEiEj × (Gij cosδij + Bij sinδij)。很多教材直接给这个公式，实际编程时用复数相量计算更简洁，两种结果完全等价，但复数计算要注意MATLAB底层是列向量还是行向量，E(:)这一步是为了强制把电动势排成列向量，保证矩阵乘法维度正确。

4.3 稳定判据：不只看绝对功角

暂态稳定最经典的判据是观察发电机间相对功角随时间的变化。在九节点系统里，通常把G2或G3作为参考机，观察其他发电机相对它的功角差。如果功角差随时间单调增大，穿越180度后不回落，基本可以认定为失稳。

具体阈值建议如下，注意这是工程经验值，不是解析边界：

实践中更常用的做法是自动搜索临界切除时间CCT。固定故障位置，逐步增加切除时间，观察系统首次失稳的临界点。这个搜索过程可以用二分法，大幅减少仿真次数。

% CCT搜索, 使用二分法 t_low = 0.1; t_high = 0.5; tolerance = 0.005; while (t_high - t_low) > tolerance t_clear = (t_low + t_high) / 2; stable = run_simulation(t_clear); % 返回是否稳定 if stable t_low = t_clear; else t_high = t_clear; end end CCT = (t_low + t_high) / 2;

注意run_simulation内部需要重新初始化状态，因为每一次故障切除时间不同，初始条件完全相同但扰动过程不同。我用这个方式做批量稳定评估时，一个CCT扫描通常要跑几十次完整仿真，在MATLAB里耗时大概几十秒到几分钟，比手动试快得多。

4.4 结果解读与曲线绘制

仿真完成后，画功角曲线的标准做法是以母线1的相角作为参考，把三台发电机的功角画在同一张图上。MATLAB里plot(t, delta*180/pi)即可，注意把弧度转成度方便阅读。稳定情况下，三台发电机的功角最终会以相同斜率增长或收敛到新的平衡点；失稳情况下，至少有一台发电机的功角相对其他机单调上升，曲线呈喇叭口展开。

5. 排错与扩展：从离线计算到更实用工程

5.1 三个高频错误与定位方法

程序跑不通，最常出问题的位置基本集中在初始化、矩阵形成和数值发散三个环节。

第一是初始化阶段的内电势幅值异常。检查方法很简单：用计算出的E'回代潮流公式，看发电机输出功率是否与设定值一致，偏差超过0.1%就说明初值算错，最常见原因是功率单位没统一或负荷恒阻抗折算方向写反。

第二是导纳矩阵奇异。发生这种情况多半是消去节点时保留节点集合为空，或者负荷阻抗过大造成病态。建议在消去前后输出矩阵条件数，条件数超过1e10就要警惕数值稳定性。九节点系统规模小，理论上不会出现这个问题，但参数改动过大时概率很高。

第三是数值发散，典型现象是仿真刚开始几步功角就跳到几千度。排查顺序：先减小步长到0.001秒看是否改善，再把阻尼项调大看趋势，最后检查故障期间导纳矩阵是否出现了对地阻抗为零的节点。多数情况下是故障矩阵形成错误，而不是积分器的问题。

5.2 验证程序的正确性：先跑无故障再跑扰动

进阶验证方法建议分三步走，保证程序可靠。

第一步，无故障仿真。在没有任何扰动的情况下运行5秒，功角应该完全不变。如果曲线出现漂移，说明初始条件本身就不是平衡点，问题出在潮流或内电势计算，不用往下查。

第二步，设一个很小的扰动并快速切除，比如0.01秒三相短路后0.02秒切除，系统应该能够恢复稳定，功角最终收敛到一个新的平衡点或绕原平衡点附近小幅摆动。这一步验证的是故障矩阵和切除逻辑是否正确。

第三步，与已知算例结果对照。九节点系统在IEEE标准算例中有公开的临界切除时间参考值，通常在三相短路位于某条主要输电线路首端时，CCT在0.3秒到0.4秒之间。如果搜索结果偏差超过0.05秒，基本可以确认模型参数或潮流初值有问题。

5.3 从离线仿真走向批量评估和在线应用

程序跑通后，真正的工程价值在于批量化和扩展性。常见的扩展方向有三个。

第一，枚举故障扫描。把九节点里的所有线路和母线分别设为故障位置，形成故障场景矩阵，循环调用仿真函数，输出每个场景的CCT或裕度指标。这个扫描在MATLAB里用parfor并行化，可以提速数倍。注意parfor循环内不要修改全局变量，所有中间数据通过函数参数传递。

第二，模型升级。把发电机从经典二阶模型升级为详细模型，在微分方程中加入励磁系统动态和调速器动态。MATLAB里实现方式是扩展状态向量，把励磁电压作为额外状态加入swing_dynamics函数。九节点系统升级到详细模型后，仿真结果更贴近实际，但调试复杂度明显上升。

第三，直接调用Matpower等工具做潮流初始化，再把自己的暂态稳定函数挂接上去。这样可以把验证过的暂态稳定计算能力复用到任意规模系统，不需要重新为每个算例写潮流程序。需要注意不同工具之间的接口命名和单位约定，我一般统一在入口处做一次数据清洗，后续代码不再处理单位问题。

最后留一个实践技巧：在仿真主循环中加入每100步输出一次当前最大功角差的逻辑，既能监控进度，又能在批量扫描中快速定位失稳时刻。对运行时间敏感的场景，用tic/toc统计每个时步的耗时，把瓶颈定位在矩阵乘法还是代数方程求解上，再做针对性优化。这套验证和优化流程，比单纯对着结果看曲线要有效得多。

本文还有配套的精品资源，点击获取