企业官网建设流程全解析

说个我经常遇到的场景：一个刚接触电力系统动态分析的学弟，拿到“多机系统静态稳定仿真——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多机静态稳定仿真拆成下面六步，后文将严格按这个顺序展开：

1. 构建系统数据：母线、支路、变压器、发电机参数；
2. 牛顿-拉夫逊潮流计算：得到各母线电压幅值和相角，这是线性化的初始运行点；
3. 发电机内电势计算：由潮流解算出每台发电机的 $E'_i\angle\delta_i$；
4. 化简网络导纳矩阵：消去所有非发电机节点，得到仅保留发电机内电势节点的等效导纳矩阵 $Y_{red}$；
5. 组装同步功率系数矩阵和状态矩阵 $A$；
6. 对 $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}$ 的规模很大，但我们只关心发电机内电势节点之间的等效关系。处理思路是：把每台发电机的内电势节点看作一个独立节点，先通过一个纯电抗支路（$jx'd$）连到机端母线，再把所有非发电机节点通过高斯消元消去，最终得到只含 $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为发电机节点，其余为负荷或联络节点。母线数据见下表。

支路和变压器参数按标幺值整理如下。变压器支路的电阻取0，电抗分别为0.0576、0.0625、0.0586。

发电机参数取经典值：母线1的 $x'_d=0.0608$，$H=23.64$；母线2的 $x'_d=0.1198$，$H=6.4$；母线3的 $x'_d=0.1813$，$H=3.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、相角给0，PV节点给指定电压幅值和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(N+k, N+k) = yg; Yext(b, N+k) = -yg; Yext(N+k, b) = -yg; end Ynn = Yext(1:N, 1:N); Yne = Yext(1:N, N+1:end); Yen = Yext(N+1:end, 1:N); Yee = Yext(N+1:end, N+1: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转换，减少出错。

整理一个自查表：

5.4 特征值结果的工程校验：和物理直觉、商业软件互证

程序写完之后不要急着下结论。我习惯做三组校验：

第一，基准工况必须稳定。三机九节点系统在额定数据下不应该出现正实部特征值，如果你算出不稳定，先检查内电势计算和 $K$ 矩阵，而不是怀疑系统。

第二，单机无穷大极限可以用手算对比。把多机降到单机模型，$K$ 直接对应 $dP_e/d\delta$，特征值的振荡频率可以用 $\sqrt{\omega_s K / (2H)}$ 估算，相差应该在百分之几以内。

第三，条件允许的话，用商业电力系统分析软件或者经典教材算例结果做一次交叉验证。特征值差个零点几可以接受，但如果符号都弄反了，肯定有bug。

5.5 下一步扩展方向

经典二阶模型只是入门。等你把这套代码跑通，最自然的扩展方向有三个：

一是加入励磁系统动态，用三阶或四阶发电机模型，多几个状态变量，特征值分析会从纯机电模式扩展到励磁模式，静态稳定的物理图景更完整。

二是把负荷改成电压相关模型甚至动态负荷模型，这在重负荷工况下会对特征值有明显影响。

三是基于参与因子设计PSS参数，把正阻尼注入到失稳模式上，把静态稳定分析和控制设计串起来。这样一套下来，你就从“会跑仿真”进化到“能分析、能设计、能解释”了。

我个人的习惯是，每写一套仿真程序，都保留一个“基准工况+一个失稳工况+一个临界工况”的最小复现集，后续改模型、改参数时非常方便。这篇文章里的步骤和代码骨架，基本就是我每次搭静态稳定仿真平台时的起点。照着这个思路走一遍，多机系统静态稳定仿真的原理、实现和工程对应关系，应该能建立一个比较扎实的整体认知。