说个我经常遇到的场景一个刚接触电力系统动态分析的学弟拿到“多机系统静态稳定仿真——MATLAB编程实践”这个题目第一反应是打开MATLAB然后发现不知道从哪下手。查资料看到一堆特征值、状态空间、功角振荡的说法脑子里全是概念就是落不了地。如果你也卡在这一步那这篇笔记应该正好适合你。我会用一台普通笔记本和MATLAB R2021b之后的任意版本从电力系统最经典的发电机二阶模型出发把一个三机九节点系统的静态稳定仿真完整跑通从潮流计算得到初始运行点到网络化简、同步功率系数矩阵、状态矩阵组装最后用特征值判断系统在小扰动下稳不稳定。整个过程不依赖任何现成的电力系统分析工具箱代码逻辑完全掌握在自己手里。多机系统静态稳定仿真这个题目核心不是“能不能算出稳定/不稳定”而是你能不能把每一步背后的物理和数学说出来。我见过不少同学用现成工具包一键出图问他特征值里实部为正代表什么、阻尼比怎么算答不上来。这不是他的问题是现代工具给了太多黑盒反而不利于建立直觉。下面我按自己写这套程序的顺序把关键环节拆开讲。1. 从“静态稳定”到“矩阵特征值”多机系统仿真的核心逻辑1.1 为什么“静态稳定”这件事值得自己写代码先回顾单机无穷大系统的经典结论。一台发电机通过线路接到无穷大母线输出功率 $P_e$ 和功角 $\delta$ 的关系是正弦曲线。静态稳定的经典判据是 $dP_e/d\delta 0$即运行点落在功角特性的上升段。小扰动下电磁功率增量与功角增量方向相反形成恢复力系统就能维持稳定。这个判据简单直接很多教材都拿它作为静态稳定的入门图景。但多机系统没有这么舒服。发电机之间有电气耦合一台机子的功角变化会通过电网影响其他机子不存在一条干净的 $P_e\delta$ 曲线。你可以把系统看成一组非线性微分方程在某个运行点附近做线性化然后用特征值判断稳定性。这就是李雅普诺夫间接法的思路只要状态矩阵的所有特征值实部都小于零系统在这个运行点就是渐近稳定的出现正实部则静态失稳。这样说下来多机系统的静态稳定分析框架其实就是“建模型、找平衡点、线性化、算特征值”四件事。MATLAB 恰好是完成这四件事非常趁手的工具。1.2 多机系统稳定分析的统一框架特征值判据状态方程写成标准形式 $\dot{x} A x$ 后特征值决定稳定性特征向量决定状态变量之间的参与关系。与单机判据相比特征值分析有两个天然优势第一它能统一处理“非周期失稳”和“振荡失稳”。如果出现正实特征值对应单调增长的失稳模式如果出现实部为正的共轭复特征值对应增幅振荡。这两种失稳在物理表象上完全不同但数学上都能被特征值的位置刻画。第二它能帮助定位薄弱环节。通过计算参与因子、特征值对参数的灵敏度我们能够知道是哪个发电机的哪个状态变量主导了这个振荡模式。这在设计电力系统稳定器、选址配置无功补偿等工程场景中非常有用。1.3 整体流水线与各模块的职责边界我把一套可复现的MATLAB多机静态稳定仿真拆成下面六步后文将严格按这个顺序展开构建系统数据母线、支路、变压器、发电机参数牛顿-拉夫逊潮流计算得到各母线电压幅值和相角这是线性化的初始运行点发电机内电势计算由潮流解算出每台发电机的 $E_i\angle\delta_i$化简网络导纳矩阵消去所有非发电机节点得到仅保留发电机内电势节点的等效导纳矩阵 $Y_{red}$组装同步功率系数矩阵和状态矩阵 $A$对 $A$ 求特征值判稳并分析振荡模式。这套流程不依赖任何电力系统专用工具箱只要掌握矩阵运算和基本数值方法就能跑通。我建议初学者不要一上来加励磁、调速器、动态负荷这些细节先把经典二阶模型玩透再逐步增加复杂度。2. 数学模型落地从摇摆方程到状态方程2.1 发电机经典二阶模型假设与适用范围多机静态稳定仿真最常用的发电机模型是经典二阶模型也叫摇摆方程模型。它的本质是把转子运动方程和功率平衡方程结合忽略励磁绕组动态、阻尼绕组效应和凸极效应认为发电机暂态电动势 $E_i$ 在扰动过程中保持恒定。这个假设在分析第一摆稳定、研究机电振荡的固有特征时是够用的尤其适合理解“多机之间的功角摇摆”这个核心问题。第 $i$ 台发电机的转子运动方程写为$$ \frac{d\delta_i}{dt} \Delta\omega_i $$$$ M_i\frac{d\Delta\omega_i}{dt} P_{mi} - P_{ei} - D_i\Delta\omega_i $$其中$\delta_i$ 是转子角$\Delta\omega_i$ 是转速偏差$P_{mi}$ 是机械功率$P_{ei}$ 是电磁功率$D_i$ 是阻尼系数$M_i$ 是惯性时间常数标幺制下通常取 $M_i 2H_i / \omega_s$$\omega_s$ 是同步角速度。注意不同教材书写习惯有差异有的把 $M_i$ 写成 $2H_i/\omega_s$有的用 $2H_i$ 直接参与运算这会导致数值上差一个 $\omega_s$ 倍是编程时最常见的坑之一后面我会专门说。2.2 网络方程与导纳矩阵化简怎么从“电网”视角看发电机把多台发电机接在同一个电网上彼此之间的电气联系需要通过网络方程表达。若保留全部母线系统节点导纳矩阵 $Y_{bus}$ 的规模很大但我们只关心发电机内电势节点之间的等效关系。处理思路是把每台发电机的内电势节点看作一个独立节点先通过一个纯电抗支路$jxd$连到机端母线再把所有非发电机节点通过高斯消元消去最终得到只含 $n$ 个发电机内电势节点的等效导纳矩阵 $Y{red}$。写成公式就是子矩阵运算$$ Y_{red} Y_{EE} - Y_{EN},Y_{NN}^{-1},Y_{NE} $$其中下标 $E$ 表示内部节点集合$N$ 表示所有需要消去的网络节点集合。这个步骤在编程里看起来只是在做分块矩阵消元物理上却很有意义它把“多机通过传输网络相互影响”这件事浓缩成了发电机内节点之间的一组等效阻抗实部 $G_{ij}$ 对应有功传递和损耗虚部 $B_{ij}$ 对应无功和相位关系。得到 $Y_{red}$ 后第 $i$ 台发电机输出的电磁功率为$$ P_{ei} E_i^2 G_{ii} \sum_{j \neq i} E_iE_j\left[ G_{ij}\cos(\delta_i - \delta_j) B_{ij}\sin(\delta_i - \delta_j) \right] $$这个式子是多机系统里最常用的功率表达式。很多延伸分析包括潮流可行性、极限传输功率、暂态稳定中的等面积定则都建立在这个表达式之上。2.3 线性化与同步功率系数矩阵静态稳定分析关心小扰动下的行为因此要在平衡点附近做线性化。设平衡点处各发电机的内电势幅值为 $E_{i0}$功角为 $\delta_{i0}$。所谓静态稳定就是看 $\delta$ 和 $\omega$ 离开平衡点后能不能回来。先对功率方程求偏导定义第 $i$ 台机对第 $j$ 台机的同步功率系数$$ K_{ij} E_iE_j\left[ B_{ij}\cos(\delta_{i0} - \delta_{j0}) - G_{ij}\sin(\delta_{i0} - \delta_{j0}) \right], \quad i \neq j $$对角元满足 $K_{ii} -\sum_{j \neq i} K_{ij}$。这样处理之后电磁功率的线性化增量可以写成$$ \Delta P_{ei} \sum_{j \neq i} K_{ij}(\Delta\delta_i - \Delta\delta_j) $$可以看到同步功率系数 $K_{ij}$ 像一组“弹性系数”把功角偏差映射成电磁功率偏差。单机无穷大系统的 $dP_e/d\delta$本质上是这台发电机与参考机之间的一个同步功率系数多机系统则有一组这样的系数构成一个矩阵。这也是为什么多机静态稳定问题的数学模型最终会落到一个矩阵的特征值分析上。2.4 状态矩阵组装与参考机选择得到同步功率系数矩阵 $K$ 以后就可以组装线性化状态方程。需要注意的是功角是相对量必须选一台参考机令其 $\Delta\delta_{ref}0$。转速偏差则是绝对的每台机的 $\Delta\omega_i$ 都可以作为状态变量。以 $n$ 台发电机为例取第 $n$ 台机为参考机状态向量为$$ x [\Delta\delta_1, \dots, \Delta\delta_{n-1}, \Delta\omega_1, \dots, \Delta\omega_n]^T $$线性化后的状态矩阵为$$ A \begin{bmatrix} 0 I_{n-1} 0_{n-1} \ -M^{-1}K_{col} -M^{-1}D \end{bmatrix} $$其中 $K_{col}$ 是 $K$ 去掉参考机那一列得到的子矩阵维数为 $n\times(n-1)$。因为参考机 $\Delta\delta_n 0$所以那一列对功率增量没有贡献。$M$ 和 $D$ 都是 $n\times n$ 对角矩阵。这个矩阵的维数是 $(2n-1)\times(2n-1)$对三机系统就是 $5\times5$手算也能验一部分结果。如果不理解参考机的处理后面算出来的特征值很可能是错的。最典型的表现是算出来一个零特征值或者数值很小的特征值看着像临界稳定其实只是坐标冗余没处理干净。3. MATLAB代码实现完整可复现的编程主线3.1 算例数据三机九节点系统的接线与参数我用的算例是电力系统课程里非常经典的三机九节点系统也就是常说的WSCC 3-machine 9-bus基准算例。基准容量100 MVA基准频率60 Hz母线1为平衡节点母线2、3为发电机节点其余为负荷或联络节点。母线数据见下表。母线类型P_G(p.u.)Q_G(p.u.)P_L(p.u.)Q_L(p.u.)V(p.u.)相角(°)1slack--001.04002PV1.63-001.025-3PV0.85-001.025-4PQ00001.0-5PQ001.250.51.0-6PQ000.90.31.0-7PQ00001.0-8PQ001.00.351.0-9PQ00001.0-支路和变压器参数按标幺值整理如下。变压器支路的电阻取0电抗分别为0.0576、0.0625、0.0586。首端母线末端母线R(p.u.)X(p.u.)B/2(p.u.)140.00000.05760450.01000.08500.088560.01700.09200.079360.00000.05860670.03200.16100.153780.00850.07200.0745820.00000.06250890.01190.10080.1045940.02200.17600.179发电机参数取经典值母线1的 $x_d0.0608$$H23.64$母线2的 $x_d0.1198$$H6.4$母线3的 $x_d0.1813$$H3.01$阻尼系数 $D$ 统一取0.1。注意 $H$ 的单位是秒不是无量纲数。3.2 潮流计算求初始运行点静态稳定分析依托于潮流解给出的运行点所以第一步是解潮流。三机九节点系统用牛顿-拉夫逊法求解收敛后得到所有母线电压幅值和相角。我记得第一次写这个程序时最耗时间的不是牛顿迭代本身而是雅可比矩阵的索引关系母线编号和数组下标对不上就会报维度错误。给出核心框架function V nr_powerflow(busData, lineData) N size(busData, 1); V busData(:, 8) .* exp(1j * deg2rad(busData(:, 9))); Y buildYbus(N, lineData); tol 1e-8; for iter 1:30 [Pcal, Qcal] calcPQ(V, Y); dP (busData(:,3) - busData(:,5)) - Pcal; dQ (busData(:,4) - busData(:,6)) - Qcal; % 平衡节点不参与PV节点Q方程不参与自行通过索引屏蔽 [J] calcJacobian(V, Y, busData); dx J \ [dP; dQ]; [V, flag] updateV(V, dx, busData); if max(abs([dP; dQ])) tol, break; end end if iter 30, warning(潮流未完全收敛请检查初值); end end这里 $P_{cal}$、$Q_{cal}$ 的计算公式是标准的节点功率注入方程。雅可比矩阵按四个分块组合$\partial P/\partial\theta$、$\partial P/\partial V$、$\partial Q/\partial\theta$、$\partial Q/\partial V$。初值方面所有PQ节点电压幅值给1.0、相角给0PV节点给指定电压幅值和0初相角通常迭代五六次就收敛。3.3 内电势与化简导纳矩阵潮流解得到机端电压 $V_i$ 后要换算出发电机内电势。先根据潮流结果计算注入电流$$ I_i \frac{P_i - jQ_i}{\bar{V}_i} $$其中 $P_i$、$Q_i$ 是发电机向网络注入的有功和无功。再计算内电势$$ E_i V_i j x_d I_i $$相角 $\delta_{i0} \angle E_i$幅值 $E_i |E_i|$。这段代码非常短但却是很多马虎之处有人忘取共轭有人把 $V_i$ 写成标量导致内电势全错。导纳矩阵化简的MATLAB函数我习惯写成这样function Yred reduceYbus(Y, genBus, xd) N size(Y, 1); n length(genBus); Yext [Y, zeros(N, n); zeros(n, N), zeros(n, n)]; for k 1:n b genBus(k); yg 1 / (1j * xd(k)); Yext(b, b) Yext(b, b) yg; Yext(Nk, Nk) yg; Yext(b, Nk) -yg; Yext(Nk, b) -yg; end Ynn Yext(1:N, 1:N); Yne Yext(1:N, N1:end); Yen Yext(N1:end, 1:N); Yee Yext(N1:end, N1:end); Yred Yee - Yen * (Ynn \ Yne); end注意 $Y_{nn}$ 是包含发电机内阻抗支路贡献后的网络节点子矩阵一般不会奇异。如果遇到奇异多半是网络中存在孤岛节点需要先修正接线数据而不是硬解。3.4 状态矩阵与特征值分析主程序有了 $E_i$、$\delta_{i0}$ 和化简后的 $Y_{red}$就可计算同步功率系数矩阵并组装状态矩阵。我给出一段浓缩了完整链路的主程序% 原始数据 genBus [1 2 3]; xd [0.0608 0.1198 0.1813]; H [23.64 6.4 3.01]; D [0.1 0.1 0.1]; omega_s 2*pi*60; n length(genBus); % 潮流与内电势 Ybus buildYbus(N, lineData); V nr_powerflow(busData, lineData); I_gen conj((P_gen - 1j*Q_gen) ./ V(genBus)); E V(genBus) 1j .* xd .* I_gen; delta0 angle(E); Emag abs(E); % 化简导纳矩阵并拆分 Yred reduceYbus(Ybus, genBus, xd); G real(Yred); B imag(Yred); % 同步功率系数矩阵 K zeros(n, n); for i 1:n for j 1:n if i ~ j d delta0(i) - delta0(j); K(i,j) Emag(i)*Emag(j)*(B(i,j)*cos(d) - G(i,j)*sin(d)); end end end for i 1:n K(i,i) -sum(K(i,:)); end % 状态矩阵参考机取第n台 Kcol K(:, 1:n-1); Mmat diag(2*H./omega_s); Dmat diag(D); A [zeros(n-1) eye(n-1) zeros(n-1,1); -Mmat\Kcol, -Mmat\Dmat]; % 特征值、阻尼比、振荡频率 [VEC, DIA] eig(A); lambda diag(DIA); sigma real(lambda); omega imag(lambda); zeta -sigma ./ sqrt(sigma.^2 omega.^2); freq_hz abs(omega) / (2*pi);这段程序跑完你能直接得到一个特征值列表。如果全部实部为负说明在当前运行点下系统静态稳定发现正实部就要警惕失稳。4. 仿真结果的解读与工况对比不止是“稳定/不稳定”四个字4.1 基准运行点的特征值结果与振荡模式我在三机九节点基准参数下跑出来的结果特征值通常由一对复共轭特征值和几个负实特征值组成。复特征值对应机电振荡模式频率一般落在0.5 Hz到2 Hz之间这是电力系统低频振荡的典型频段。阻尼比是评估小扰动稳定性的重要指标工程上一般认为阻尼比低于2%就需要关注低于0就是负阻尼系统一旦受到扰动振荡会持续增长。以我的算例为例基准工况下两个主要振荡模式的频率大约在1.1 Hz和1.5 Hz阻尼比约5%~10%。第一个模式通常与发电机1、2之间的功角摇摆相关第二个模式与发电机3相对其他机的摇摆相关。可以借用特征值分布图直观判断所有特征值都在左半平面系统稳定某对特征值离虚轴越近对应的模式越危险。4.2 负荷加重/单线断开时的稳定性演变静态稳定仿真最有价值的用途之一是观察系统在什么运行条件下会失稳。我做了两组典型实验第一组把所有负荷等比例增加到1.6倍。潮流解中功角差明显变大同步功率系数矩阵的数值变小某些 $K_{ij}$ 甚至变为负数。特征值随之向虚轴移动先是一对复特征值的实部从负变正系统呈现增幅振荡失稳。这是典型的“静态功角失稳”直观对应单机曲线里的运行点越过功角特性峰值。第二组把母线5-6间的线路断开模拟N-1事故。由于网络拓扑变化$Y_{red}$ 发生变化发电机之间的电气联系减弱同步功率系数下降特征值实部往往更靠近虚轴。这说明单一输电通道断开后系统静态稳定裕度显著下降。这两组实验做下来你对“稳定裕度”的感受会比只报一个稳定/不稳定的结论深刻得多。4.3 阻尼、惯性时间常数对特征值的影响把 $D$ 从0.1改成0特征值实部会明显变小甚至接近零把 $D$ 改成0.3阻尼比明显提升。这告诉大家一个工程常识阻尼是维持小扰动稳定的重要因素。经典模型中 $D$ 是个聚合参数并不真的对应某个物理阻尼器而是反映了励磁系统、调速器、负荷特性等综合效应。所以在后续做详细建模时励磁系统的负阻尼效应有时会让系统失稳这就是电力系统稳定器PSS存在的意义。惯性时间常数 $H$ 的影响也很有趣。小电机通常 $H$ 小转速对功率不平衡更敏感振荡频率偏高大电机 $H$ 大振荡频率偏低。改变 $H$同一振荡模式的频率会明显变化但特征值实部未必成比例变化。这些实验在MATLAB里做都很简单改一个参数再运行就行适合用来建立物理直觉。4.4 特征向量与参与因子快速找出“危险区域”特征值只看“稳不稳”特征向量可以回答“谁在晃”。计算状态矩阵特征值对应的特征向量观察各台发电机 $\Delta\omega_i$ 分量的幅值和相位就能判断某个振荡模式里哪些发电机是同调摇摆、哪些是反向摇摆。正式做法是计算参与因子即特征向量左、右矩阵对应元素的乘积。举例来说如果第一个振荡模式的参与因子在发电机2上最大说明这个模式主要是发电机2相对其余机组的局部振荡。配置稳定器时优先在参与因子大的机组上安装效果最明显。这一层分析把静态稳定仿真从“稳定分析”延伸到了“稳定控制设计”是后续学习的小入口也是面试或答辩时非常容易出彩的扩展点。5. 调试与改进这套仿真最容易翻车的几个地方5.1 潮流初值和收敛性静态稳定仿真全链路里潮流不收敛是最常见的第一个障碍。原因多半出在初值上PV节点电压幅值给的太离谱或者负荷水平太高导致潮流无解。调试技巧很简单先把所有PQ节点电压幅值设为1.0相角设为0逐步增加负荷水平看哪个工况开始不收敛。那通常就是静态稳定极限附近恰好是失稳前兆。另外牛顿-拉夫逊法的雅可比矩阵如果不对会表现为迭代振荡或不收敛。建议在每一步迭代后打印最大功率不平衡量如果从$10^{-2}$量级降到$10^{-8}$量级说明雅可比是对的如果卡在某个值不动多半是有功、无功方程的索引屏蔽写错尤其是PV节点和无功方程的关系。5.2 导纳化简的奇异问题化简导纳矩阵用到 $Y_{nn}^{-1}$。如果网络数据存在孤岛或者某台发电机的内阻抗支路没有把机端母线连进网络$Y_{nn}$ 可能是奇异的。我在处理某个自定义算例时遇到过两台发电机之间只有电容器支路没有实际电抗回路潮流本身才勉强有解化简时矩阵接近奇异。处理方式是回到母线数据检查连通性给孤岛节点加合理的并联导纳或者换用带微小的正则化项但这只能作为临时方案不能掩盖数据问题。5.3 单位与量纲H、M、ω_s、角度这是我反复强调的坑。用 $M 2H/\omega_s$ 时$\omega_s$ 必须以 rad/s 为单位60 Hz系统的 $\omega_s 120\pi \approx 376.99$。如果你直接用 $M 2H$特征值会整体放大 $\omega_s$ 倍振荡频率直接从0.1 Hz变成10 Hz量级一看就知道不对。角度的弧度/度混用也很常见尤其是潮流结果里电压相角习惯用度而内电势计算和 $K$ 矩阵里必须用弧度。我在程序里统一用deg2rad转换减少出错。整理一个自查表量常用单位说明$\delta$rad状态方程里必须用弧度$H$s发电机惯性时间常数$M$s$M2H/\omega_s$标幺系统下$D$p.u.阻尼系数$\omega_s$rad/s同步角速度60Hz时为$2\pi\times60$特征值实部1/s反映衰减速度阻尼比无量纲$-\sigma/\sqrt{\sigma^2\omega^2}$5.4 特征值结果的工程校验和物理直觉、商业软件互证程序写完之后不要急着下结论。我习惯做三组校验第一基准工况必须稳定。三机九节点系统在额定数据下不应该出现正实部特征值如果你算出不稳定先检查内电势计算和 $K$ 矩阵而不是怀疑系统。第二单机无穷大极限可以用手算对比。把多机降到单机模型$K$ 直接对应 $dP_e/d\delta$特征值的振荡频率可以用 $\sqrt{\omega_s K / (2H)}$ 估算相差应该在百分之几以内。第三条件允许的话用商业电力系统分析软件或者经典教材算例结果做一次交叉验证。特征值差个零点几可以接受但如果符号都弄反了肯定有bug。5.5 下一步扩展方向经典二阶模型只是入门。等你把这套代码跑通最自然的扩展方向有三个一是加入励磁系统动态用三阶或四阶发电机模型多几个状态变量特征值分析会从纯机电模式扩展到励磁模式静态稳定的物理图景更完整。二是把负荷改成电压相关模型甚至动态负荷模型这在重负荷工况下会对特征值有明显影响。三是基于参与因子设计PSS参数把正阻尼注入到失稳模式上把静态稳定分析和控制设计串起来。这样一套下来你就从“会跑仿真”进化到“能分析、能设计、能解释”了。我个人的习惯是每写一套仿真程序都保留一个“基准工况一个失稳工况一个临界工况”的最小复现集后续改模型、改参数时非常方便。这篇文章里的步骤和代码骨架基本就是我每次搭静态稳定仿真平台时的起点。照着这个思路走一遍多机系统静态稳定仿真的原理、实现和工程对应关系应该能建立一个比较扎实的整体认知。
多机系统静态稳定仿真：MATLAB编程与特征值分析实战
说个我经常遇到的场景一个刚接触电力系统动态分析的学弟拿到“多机系统静态稳定仿真——MATLAB编程实践”这个题目第一反应是打开MATLAB然后发现不知道从哪下手。查资料看到一堆特征值、状态空间、功角振荡的说法脑子里全是概念就是落不了地。如果你也卡在这一步那这篇笔记应该正好适合你。我会用一台普通笔记本和MATLAB R2021b之后的任意版本从电力系统最经典的发电机二阶模型出发把一个三机九节点系统的静态稳定仿真完整跑通从潮流计算得到初始运行点到网络化简、同步功率系数矩阵、状态矩阵组装最后用特征值判断系统在小扰动下稳不稳定。整个过程不依赖任何现成的电力系统分析工具箱代码逻辑完全掌握在自己手里。多机系统静态稳定仿真这个题目核心不是“能不能算出稳定/不稳定”而是你能不能把每一步背后的物理和数学说出来。我见过不少同学用现成工具包一键出图问他特征值里实部为正代表什么、阻尼比怎么算答不上来。这不是他的问题是现代工具给了太多黑盒反而不利于建立直觉。下面我按自己写这套程序的顺序把关键环节拆开讲。1. 从“静态稳定”到“矩阵特征值”多机系统仿真的核心逻辑1.1 为什么“静态稳定”这件事值得自己写代码先回顾单机无穷大系统的经典结论。一台发电机通过线路接到无穷大母线输出功率 $P_e$ 和功角 $\delta$ 的关系是正弦曲线。静态稳定的经典判据是 $dP_e/d\delta 0$即运行点落在功角特性的上升段。小扰动下电磁功率增量与功角增量方向相反形成恢复力系统就能维持稳定。这个判据简单直接很多教材都拿它作为静态稳定的入门图景。但多机系统没有这么舒服。发电机之间有电气耦合一台机子的功角变化会通过电网影响其他机子不存在一条干净的 $P_e\delta$ 曲线。你可以把系统看成一组非线性微分方程在某个运行点附近做线性化然后用特征值判断稳定性。这就是李雅普诺夫间接法的思路只要状态矩阵的所有特征值实部都小于零系统在这个运行点就是渐近稳定的出现正实部则静态失稳。这样说下来多机系统的静态稳定分析框架其实就是“建模型、找平衡点、线性化、算特征值”四件事。MATLAB 恰好是完成这四件事非常趁手的工具。1.2 多机系统稳定分析的统一框架特征值判据状态方程写成标准形式 $\dot{x} A x$ 后特征值决定稳定性特征向量决定状态变量之间的参与关系。与单机判据相比特征值分析有两个天然优势第一它能统一处理“非周期失稳”和“振荡失稳”。如果出现正实特征值对应单调增长的失稳模式如果出现实部为正的共轭复特征值对应增幅振荡。这两种失稳在物理表象上完全不同但数学上都能被特征值的位置刻画。第二它能帮助定位薄弱环节。通过计算参与因子、特征值对参数的灵敏度我们能够知道是哪个发电机的哪个状态变量主导了这个振荡模式。这在设计电力系统稳定器、选址配置无功补偿等工程场景中非常有用。1.3 整体流水线与各模块的职责边界我把一套可复现的MATLAB多机静态稳定仿真拆成下面六步后文将严格按这个顺序展开构建系统数据母线、支路、变压器、发电机参数牛顿-拉夫逊潮流计算得到各母线电压幅值和相角这是线性化的初始运行点发电机内电势计算由潮流解算出每台发电机的 $E_i\angle\delta_i$化简网络导纳矩阵消去所有非发电机节点得到仅保留发电机内电势节点的等效导纳矩阵 $Y_{red}$组装同步功率系数矩阵和状态矩阵 $A$对 $A$ 求特征值判稳并分析振荡模式。这套流程不依赖任何电力系统专用工具箱只要掌握矩阵运算和基本数值方法就能跑通。我建议初学者不要一上来加励磁、调速器、动态负荷这些细节先把经典二阶模型玩透再逐步增加复杂度。2. 数学模型落地从摇摆方程到状态方程2.1 发电机经典二阶模型假设与适用范围多机静态稳定仿真最常用的发电机模型是经典二阶模型也叫摇摆方程模型。它的本质是把转子运动方程和功率平衡方程结合忽略励磁绕组动态、阻尼绕组效应和凸极效应认为发电机暂态电动势 $E_i$ 在扰动过程中保持恒定。这个假设在分析第一摆稳定、研究机电振荡的固有特征时是够用的尤其适合理解“多机之间的功角摇摆”这个核心问题。第 $i$ 台发电机的转子运动方程写为$$ \frac{d\delta_i}{dt} \Delta\omega_i $$$$ M_i\frac{d\Delta\omega_i}{dt} P_{mi} - P_{ei} - D_i\Delta\omega_i $$其中$\delta_i$ 是转子角$\Delta\omega_i$ 是转速偏差$P_{mi}$ 是机械功率$P_{ei}$ 是电磁功率$D_i$ 是阻尼系数$M_i$ 是惯性时间常数标幺制下通常取 $M_i 2H_i / \omega_s$$\omega_s$ 是同步角速度。注意不同教材书写习惯有差异有的把 $M_i$ 写成 $2H_i/\omega_s$有的用 $2H_i$ 直接参与运算这会导致数值上差一个 $\omega_s$ 倍是编程时最常见的坑之一后面我会专门说。2.2 网络方程与导纳矩阵化简怎么从“电网”视角看发电机把多台发电机接在同一个电网上彼此之间的电气联系需要通过网络方程表达。若保留全部母线系统节点导纳矩阵 $Y_{bus}$ 的规模很大但我们只关心发电机内电势节点之间的等效关系。处理思路是把每台发电机的内电势节点看作一个独立节点先通过一个纯电抗支路$jxd$连到机端母线再把所有非发电机节点通过高斯消元消去最终得到只含 $n$ 个发电机内电势节点的等效导纳矩阵 $Y{red}$。写成公式就是子矩阵运算$$ Y_{red} Y_{EE} - Y_{EN},Y_{NN}^{-1},Y_{NE} $$其中下标 $E$ 表示内部节点集合$N$ 表示所有需要消去的网络节点集合。这个步骤在编程里看起来只是在做分块矩阵消元物理上却很有意义它把“多机通过传输网络相互影响”这件事浓缩成了发电机内节点之间的一组等效阻抗实部 $G_{ij}$ 对应有功传递和损耗虚部 $B_{ij}$ 对应无功和相位关系。得到 $Y_{red}$ 后第 $i$ 台发电机输出的电磁功率为$$ P_{ei} E_i^2 G_{ii} \sum_{j \neq i} E_iE_j\left[ G_{ij}\cos(\delta_i - \delta_j) B_{ij}\sin(\delta_i - \delta_j) \right] $$这个式子是多机系统里最常用的功率表达式。很多延伸分析包括潮流可行性、极限传输功率、暂态稳定中的等面积定则都建立在这个表达式之上。2.3 线性化与同步功率系数矩阵静态稳定分析关心小扰动下的行为因此要在平衡点附近做线性化。设平衡点处各发电机的内电势幅值为 $E_{i0}$功角为 $\delta_{i0}$。所谓静态稳定就是看 $\delta$ 和 $\omega$ 离开平衡点后能不能回来。先对功率方程求偏导定义第 $i$ 台机对第 $j$ 台机的同步功率系数$$ K_{ij} E_iE_j\left[ B_{ij}\cos(\delta_{i0} - \delta_{j0}) - G_{ij}\sin(\delta_{i0} - \delta_{j0}) \right], \quad i \neq j $$对角元满足 $K_{ii} -\sum_{j \neq i} K_{ij}$。这样处理之后电磁功率的线性化增量可以写成$$ \Delta P_{ei} \sum_{j \neq i} K_{ij}(\Delta\delta_i - \Delta\delta_j) $$可以看到同步功率系数 $K_{ij}$ 像一组“弹性系数”把功角偏差映射成电磁功率偏差。单机无穷大系统的 $dP_e/d\delta$本质上是这台发电机与参考机之间的一个同步功率系数多机系统则有一组这样的系数构成一个矩阵。这也是为什么多机静态稳定问题的数学模型最终会落到一个矩阵的特征值分析上。2.4 状态矩阵组装与参考机选择得到同步功率系数矩阵 $K$ 以后就可以组装线性化状态方程。需要注意的是功角是相对量必须选一台参考机令其 $\Delta\delta_{ref}0$。转速偏差则是绝对的每台机的 $\Delta\omega_i$ 都可以作为状态变量。以 $n$ 台发电机为例取第 $n$ 台机为参考机状态向量为$$ x [\Delta\delta_1, \dots, \Delta\delta_{n-1}, \Delta\omega_1, \dots, \Delta\omega_n]^T $$线性化后的状态矩阵为$$ A \begin{bmatrix} 0 I_{n-1} 0_{n-1} \ -M^{-1}K_{col} -M^{-1}D \end{bmatrix} $$其中 $K_{col}$ 是 $K$ 去掉参考机那一列得到的子矩阵维数为 $n\times(n-1)$。因为参考机 $\Delta\delta_n 0$所以那一列对功率增量没有贡献。$M$ 和 $D$ 都是 $n\times n$ 对角矩阵。这个矩阵的维数是 $(2n-1)\times(2n-1)$对三机系统就是 $5\times5$手算也能验一部分结果。如果不理解参考机的处理后面算出来的特征值很可能是错的。最典型的表现是算出来一个零特征值或者数值很小的特征值看着像临界稳定其实只是坐标冗余没处理干净。3. MATLAB代码实现完整可复现的编程主线3.1 算例数据三机九节点系统的接线与参数我用的算例是电力系统课程里非常经典的三机九节点系统也就是常说的WSCC 3-machine 9-bus基准算例。基准容量100 MVA基准频率60 Hz母线1为平衡节点母线2、3为发电机节点其余为负荷或联络节点。母线数据见下表。母线类型P_G(p.u.)Q_G(p.u.)P_L(p.u.)Q_L(p.u.)V(p.u.)相角(°)1slack--001.04002PV1.63-001.025-3PV0.85-001.025-4PQ00001.0-5PQ001.250.51.0-6PQ000.90.31.0-7PQ00001.0-8PQ001.00.351.0-9PQ00001.0-支路和变压器参数按标幺值整理如下。变压器支路的电阻取0电抗分别为0.0576、0.0625、0.0586。首端母线末端母线R(p.u.)X(p.u.)B/2(p.u.)140.00000.05760450.01000.08500.088560.01700.09200.079360.00000.05860670.03200.16100.153780.00850.07200.0745820.00000.06250890.01190.10080.1045940.02200.17600.179发电机参数取经典值母线1的 $x_d0.0608$$H23.64$母线2的 $x_d0.1198$$H6.4$母线3的 $x_d0.1813$$H3.01$阻尼系数 $D$ 统一取0.1。注意 $H$ 的单位是秒不是无量纲数。3.2 潮流计算求初始运行点静态稳定分析依托于潮流解给出的运行点所以第一步是解潮流。三机九节点系统用牛顿-拉夫逊法求解收敛后得到所有母线电压幅值和相角。我记得第一次写这个程序时最耗时间的不是牛顿迭代本身而是雅可比矩阵的索引关系母线编号和数组下标对不上就会报维度错误。给出核心框架function V nr_powerflow(busData, lineData) N size(busData, 1); V busData(:, 8) .* exp(1j * deg2rad(busData(:, 9))); Y buildYbus(N, lineData); tol 1e-8; for iter 1:30 [Pcal, Qcal] calcPQ(V, Y); dP (busData(:,3) - busData(:,5)) - Pcal; dQ (busData(:,4) - busData(:,6)) - Qcal; % 平衡节点不参与PV节点Q方程不参与自行通过索引屏蔽 [J] calcJacobian(V, Y, busData); dx J \ [dP; dQ]; [V, flag] updateV(V, dx, busData); if max(abs([dP; dQ])) tol, break; end end if iter 30, warning(潮流未完全收敛请检查初值); end end这里 $P_{cal}$、$Q_{cal}$ 的计算公式是标准的节点功率注入方程。雅可比矩阵按四个分块组合$\partial P/\partial\theta$、$\partial P/\partial V$、$\partial Q/\partial\theta$、$\partial Q/\partial V$。初值方面所有PQ节点电压幅值给1.0、相角给0PV节点给指定电压幅值和0初相角通常迭代五六次就收敛。3.3 内电势与化简导纳矩阵潮流解得到机端电压 $V_i$ 后要换算出发电机内电势。先根据潮流结果计算注入电流$$ I_i \frac{P_i - jQ_i}{\bar{V}_i} $$其中 $P_i$、$Q_i$ 是发电机向网络注入的有功和无功。再计算内电势$$ E_i V_i j x_d I_i $$相角 $\delta_{i0} \angle E_i$幅值 $E_i |E_i|$。这段代码非常短但却是很多马虎之处有人忘取共轭有人把 $V_i$ 写成标量导致内电势全错。导纳矩阵化简的MATLAB函数我习惯写成这样function Yred reduceYbus(Y, genBus, xd) N size(Y, 1); n length(genBus); Yext [Y, zeros(N, n); zeros(n, N), zeros(n, n)]; for k 1:n b genBus(k); yg 1 / (1j * xd(k)); Yext(b, b) Yext(b, b) yg; Yext(Nk, Nk) yg; Yext(b, Nk) -yg; Yext(Nk, b) -yg; end Ynn Yext(1:N, 1:N); Yne Yext(1:N, N1:end); Yen Yext(N1:end, 1:N); Yee Yext(N1:end, N1:end); Yred Yee - Yen * (Ynn \ Yne); end注意 $Y_{nn}$ 是包含发电机内阻抗支路贡献后的网络节点子矩阵一般不会奇异。如果遇到奇异多半是网络中存在孤岛节点需要先修正接线数据而不是硬解。3.4 状态矩阵与特征值分析主程序有了 $E_i$、$\delta_{i0}$ 和化简后的 $Y_{red}$就可计算同步功率系数矩阵并组装状态矩阵。我给出一段浓缩了完整链路的主程序% 原始数据 genBus [1 2 3]; xd [0.0608 0.1198 0.1813]; H [23.64 6.4 3.01]; D [0.1 0.1 0.1]; omega_s 2*pi*60; n length(genBus); % 潮流与内电势 Ybus buildYbus(N, lineData); V nr_powerflow(busData, lineData); I_gen conj((P_gen - 1j*Q_gen) ./ V(genBus)); E V(genBus) 1j .* xd .* I_gen; delta0 angle(E); Emag abs(E); % 化简导纳矩阵并拆分 Yred reduceYbus(Ybus, genBus, xd); G real(Yred); B imag(Yred); % 同步功率系数矩阵 K zeros(n, n); for i 1:n for j 1:n if i ~ j d delta0(i) - delta0(j); K(i,j) Emag(i)*Emag(j)*(B(i,j)*cos(d) - G(i,j)*sin(d)); end end end for i 1:n K(i,i) -sum(K(i,:)); end % 状态矩阵参考机取第n台 Kcol K(:, 1:n-1); Mmat diag(2*H./omega_s); Dmat diag(D); A [zeros(n-1) eye(n-1) zeros(n-1,1); -Mmat\Kcol, -Mmat\Dmat]; % 特征值、阻尼比、振荡频率 [VEC, DIA] eig(A); lambda diag(DIA); sigma real(lambda); omega imag(lambda); zeta -sigma ./ sqrt(sigma.^2 omega.^2); freq_hz abs(omega) / (2*pi);这段程序跑完你能直接得到一个特征值列表。如果全部实部为负说明在当前运行点下系统静态稳定发现正实部就要警惕失稳。4. 仿真结果的解读与工况对比不止是“稳定/不稳定”四个字4.1 基准运行点的特征值结果与振荡模式我在三机九节点基准参数下跑出来的结果特征值通常由一对复共轭特征值和几个负实特征值组成。复特征值对应机电振荡模式频率一般落在0.5 Hz到2 Hz之间这是电力系统低频振荡的典型频段。阻尼比是评估小扰动稳定性的重要指标工程上一般认为阻尼比低于2%就需要关注低于0就是负阻尼系统一旦受到扰动振荡会持续增长。以我的算例为例基准工况下两个主要振荡模式的频率大约在1.1 Hz和1.5 Hz阻尼比约5%~10%。第一个模式通常与发电机1、2之间的功角摇摆相关第二个模式与发电机3相对其他机的摇摆相关。可以借用特征值分布图直观判断所有特征值都在左半平面系统稳定某对特征值离虚轴越近对应的模式越危险。4.2 负荷加重/单线断开时的稳定性演变静态稳定仿真最有价值的用途之一是观察系统在什么运行条件下会失稳。我做了两组典型实验第一组把所有负荷等比例增加到1.6倍。潮流解中功角差明显变大同步功率系数矩阵的数值变小某些 $K_{ij}$ 甚至变为负数。特征值随之向虚轴移动先是一对复特征值的实部从负变正系统呈现增幅振荡失稳。这是典型的“静态功角失稳”直观对应单机曲线里的运行点越过功角特性峰值。第二组把母线5-6间的线路断开模拟N-1事故。由于网络拓扑变化$Y_{red}$ 发生变化发电机之间的电气联系减弱同步功率系数下降特征值实部往往更靠近虚轴。这说明单一输电通道断开后系统静态稳定裕度显著下降。这两组实验做下来你对“稳定裕度”的感受会比只报一个稳定/不稳定的结论深刻得多。4.3 阻尼、惯性时间常数对特征值的影响把 $D$ 从0.1改成0特征值实部会明显变小甚至接近零把 $D$ 改成0.3阻尼比明显提升。这告诉大家一个工程常识阻尼是维持小扰动稳定的重要因素。经典模型中 $D$ 是个聚合参数并不真的对应某个物理阻尼器而是反映了励磁系统、调速器、负荷特性等综合效应。所以在后续做详细建模时励磁系统的负阻尼效应有时会让系统失稳这就是电力系统稳定器PSS存在的意义。惯性时间常数 $H$ 的影响也很有趣。小电机通常 $H$ 小转速对功率不平衡更敏感振荡频率偏高大电机 $H$ 大振荡频率偏低。改变 $H$同一振荡模式的频率会明显变化但特征值实部未必成比例变化。这些实验在MATLAB里做都很简单改一个参数再运行就行适合用来建立物理直觉。4.4 特征向量与参与因子快速找出“危险区域”特征值只看“稳不稳”特征向量可以回答“谁在晃”。计算状态矩阵特征值对应的特征向量观察各台发电机 $\Delta\omega_i$ 分量的幅值和相位就能判断某个振荡模式里哪些发电机是同调摇摆、哪些是反向摇摆。正式做法是计算参与因子即特征向量左、右矩阵对应元素的乘积。举例来说如果第一个振荡模式的参与因子在发电机2上最大说明这个模式主要是发电机2相对其余机组的局部振荡。配置稳定器时优先在参与因子大的机组上安装效果最明显。这一层分析把静态稳定仿真从“稳定分析”延伸到了“稳定控制设计”是后续学习的小入口也是面试或答辩时非常容易出彩的扩展点。5. 调试与改进这套仿真最容易翻车的几个地方5.1 潮流初值和收敛性静态稳定仿真全链路里潮流不收敛是最常见的第一个障碍。原因多半出在初值上PV节点电压幅值给的太离谱或者负荷水平太高导致潮流无解。调试技巧很简单先把所有PQ节点电压幅值设为1.0相角设为0逐步增加负荷水平看哪个工况开始不收敛。那通常就是静态稳定极限附近恰好是失稳前兆。另外牛顿-拉夫逊法的雅可比矩阵如果不对会表现为迭代振荡或不收敛。建议在每一步迭代后打印最大功率不平衡量如果从$10^{-2}$量级降到$10^{-8}$量级说明雅可比是对的如果卡在某个值不动多半是有功、无功方程的索引屏蔽写错尤其是PV节点和无功方程的关系。5.2 导纳化简的奇异问题化简导纳矩阵用到 $Y_{nn}^{-1}$。如果网络数据存在孤岛或者某台发电机的内阻抗支路没有把机端母线连进网络$Y_{nn}$ 可能是奇异的。我在处理某个自定义算例时遇到过两台发电机之间只有电容器支路没有实际电抗回路潮流本身才勉强有解化简时矩阵接近奇异。处理方式是回到母线数据检查连通性给孤岛节点加合理的并联导纳或者换用带微小的正则化项但这只能作为临时方案不能掩盖数据问题。5.3 单位与量纲H、M、ω_s、角度这是我反复强调的坑。用 $M 2H/\omega_s$ 时$\omega_s$ 必须以 rad/s 为单位60 Hz系统的 $\omega_s 120\pi \approx 376.99$。如果你直接用 $M 2H$特征值会整体放大 $\omega_s$ 倍振荡频率直接从0.1 Hz变成10 Hz量级一看就知道不对。角度的弧度/度混用也很常见尤其是潮流结果里电压相角习惯用度而内电势计算和 $K$ 矩阵里必须用弧度。我在程序里统一用deg2rad转换减少出错。整理一个自查表量常用单位说明$\delta$rad状态方程里必须用弧度$H$s发电机惯性时间常数$M$s$M2H/\omega_s$标幺系统下$D$p.u.阻尼系数$\omega_s$rad/s同步角速度60Hz时为$2\pi\times60$特征值实部1/s反映衰减速度阻尼比无量纲$-\sigma/\sqrt{\sigma^2\omega^2}$5.4 特征值结果的工程校验和物理直觉、商业软件互证程序写完之后不要急着下结论。我习惯做三组校验第一基准工况必须稳定。三机九节点系统在额定数据下不应该出现正实部特征值如果你算出不稳定先检查内电势计算和 $K$ 矩阵而不是怀疑系统。第二单机无穷大极限可以用手算对比。把多机降到单机模型$K$ 直接对应 $dP_e/d\delta$特征值的振荡频率可以用 $\sqrt{\omega_s K / (2H)}$ 估算相差应该在百分之几以内。第三条件允许的话用商业电力系统分析软件或者经典教材算例结果做一次交叉验证。特征值差个零点几可以接受但如果符号都弄反了肯定有bug。5.5 下一步扩展方向经典二阶模型只是入门。等你把这套代码跑通最自然的扩展方向有三个一是加入励磁系统动态用三阶或四阶发电机模型多几个状态变量特征值分析会从纯机电模式扩展到励磁模式静态稳定的物理图景更完整。二是把负荷改成电压相关模型甚至动态负荷模型这在重负荷工况下会对特征值有明显影响。三是基于参与因子设计PSS参数把正阻尼注入到失稳模式上把静态稳定分析和控制设计串起来。这样一套下来你就从“会跑仿真”进化到“能分析、能设计、能解释”了。我个人的习惯是每写一套仿真程序都保留一个“基准工况一个失稳工况一个临界工况”的最小复现集后续改模型、改参数时非常方便。这篇文章里的步骤和代码骨架基本就是我每次搭静态稳定仿真平台时的起点。照着这个思路走一遍多机系统静态稳定仿真的原理、实现和工程对应关系应该能建立一个比较扎实的整体认知。