简介：一份基于Kundur两区域四机三模式经典模型的小信号稳定性分析MATLAB资源，专为电力系统专业学生、研究生及工程师设计，可用于理解《电力系统稳定与控制》中的核心算例与线性化分析方法。资源包共2个文件，含1个MATLAB脚本和1个配套说明文档，脚本覆盖系统参数设置、动态模型定义、扰动施加与仿真输出全过程，整体仅6KB，轻量易用。目前已有211人学习下载。通过运行脚本，可直观观察小扰动下系统频率、电压的动态响应曲线，并结合特征值实部正负快速判断系统稳定性，从而掌握线性化建模与状态空间分析的核心流程。资源完整复现了三模式运行工况，既能帮助学习者打通从理论模型到代码实践的关键链条，也能作为课堂教学演示或科研验证的实用工具，在负荷波动、局部故障等场景下揭示小信号稳定的核心机理。

1. kundur_12o3_v4.zip 到底是什么：小信号稳定分析的通用底座

拿到kundur_12o3_v4.zip这类文件，大多数人第一反应是解压、找 main 脚本、直接跑出几条曲线。但真正决定这个包值不值得留的，是里面那套“小信号稳定”状态矩阵组装逻辑能不能改、能不能信。Kundur 体系几乎是电力系统动态稳定方向的默认基准：四机两区域模型、区间低频振荡、励磁与 PSS 的交互，全都在这套框架里讲透。标题里的kundur锁定的是这个体系，而不是某个不可替代的软件版本。

我准备把这类模型包背后通用的一条落地方案讲清楚：从潮流平衡点出发，把发电机与励磁系统线性化为状态矩阵 A，用 MATLAB 的特征值分解拿到振荡频率与阻尼比，再用参与因子定位“振荡到底从哪台机来”。后文会具体拆解12o3_v4这个版本号里的常见改动方向，并给出可以直接跑的脚本骨架。适合谁看：研究生做稳定分析、工程师审 PSS 参数、写论文需要复现经典算例的从业者。

2. 先看懂 Kundur 两区域四机系统：小信号稳定为什么拿它当基准

2.1 两区域四机系统的结构与振荡模式划分

Kundur 体系里最常出现的作战场地就是两区域四机系统：区域 1 两台发电机，区域 2 两台发电机，中间通过双回联络线传递功率，典型工况是区域 1 向区域 2 送电几百兆瓦。这个拓扑看着简单，但它把小信号稳定的两类核心问题都装进去了——局部振荡模式，单台机对相邻机组来回摆；区间振荡模式，区域 1 整体对区域 2 整体摆动，频率一般落在 0.5~0.7 Hz 之间。

工程上判断小信号稳定是否达标，看的不是特征值实部到底负多少，而是模式阻尼比。阻尼比低于 3% 就属于弱阻尼，低于 2% 基本要安排 PSS 或调整励磁参数。频率决定这是局部模式还是区间模式，阻尼比决定要不要加控制措施，两者缺一不可。拿到任何 kundur 类模型包，第一步就该把这两张表算出来，而不是急着看扰动仿真曲线。

2.2 小信号稳定的数学骨架：从非线性微分方程到状态矩阵 A

小信号稳定的本质是把电力系统在稳态平衡点附近做泰勒展开，丢掉二次项，得到一个线性时不变系统：Δẋ = AΔx。状态变量通常包括发电机功角 δ、转速偏差 Δω、暂态电势 Eq′，励磁系统加一个 Efd 状态，PSS 投入的话再挂两到三个控制状态。这个 A 矩阵的特征值，决定了扰动后振荡的形态和衰减速度。

A 矩阵的组装配料表在业界有个更常听到的名字：Heffron-Phillips 模型。它把电磁转矩拆成同步转矩与阻尼转矩两部分，用 K1~K6 六个偏导系数描述功角、暂态电势、励磁电压之间的耦合关系。K1 是电磁转矩对功角的偏导，K2 是电磁转矩对 Eq′ 的偏导，K3、K4 描述励磁绕组动态，K5、K6 描述机端电压反馈。后面写代码时，这几个系数就是状态矩阵里的原料。

2.3 12o3_v4 里最需要先确认的三个模型开关

12o3_v4这个编号我不做过度解读，按这类仿真包的命名习惯，多半是状态维数、区域数或迭代版本的简写。但版本迭代最常见的改动点，集中在这三个开关上。

第一是励磁机模型。用直流励磁 DC1A 还是静止励磁 ST1A，直接决定励磁时间常数与增益范围。静止励磁响应快、增益高，更容易在区间模式上贡献负阻尼。第二是 PSS 默认投不投入。同一套 Kundur 参数，有 PSS 和没 PSS 的阻尼比可以差出三倍。第三是负荷模型用恒功率还是恒阻抗。恒功率负荷对振荡的阻尼效应与恒阻抗负荷差异明显。拿到包内脚本前，先把这三个点勾出来，不然后面调参数翻车了都不知道在调什么。

3. 用 MATLAB 把模型跑起来：平衡点、状态矩阵与特征值提取

3.1 平衡点初值准备：潮流输出怎么变成状态初值

小信号分析最容易被忽略的前置步骤是平衡点初始化。A 矩阵必须在稳态工作点处线性化，平衡点给错了，后面算出来的特征值全是空中楼阁。常见做法是先用潮流工具算出母线电压幅值、相角和发电机出力，再把这些量换算成状态变量的初值。

% 以两区域四机系统、公共基准 100 MVA 为例的一组典型平衡点初值 % 电压单位：pu，功角单位：rad V = [1.030; 1.010; 1.030; 1.010]; % 四台机端电压幅值 theta = [0.12; 0.32; -0.14; -0.06]; % 四台机端电压相角，注意差异要符合潮流方向 Pe = [700; 700; 700; 700] / 900; % 发电机有功出力，折到公共基准后的 pu 值 Pd = [967; 967] / 900; % 两个区域负荷，pu % 状态初值顺序：delta、omega、Eqp delta0 = theta - theta(1); % 以 1 号机功角为参考点 omega0 = zeros(4, 1); % 稳态下转速偏差为 0 % Eqp 理论上是潮流隐式解，实际应由 dq 变换给出 Eqp0 = V + [0.35; 0.32; 0.31; 0.33]; % 暂态电势近似初值

这段代码的关键逻辑在于：功角初值用的是相对值而不是绝对值，因为摇摆方程只对功角差敏感；转速偏差在稳态下必须清零，这是线性化点被严格定义的前提；Eqp0 在简化脚本里是估算值，真实项目中要从 dq 分解得到，否则励磁回路的状态方程会带着一个虚假的初始误差进入线性化。我一般会先用潮流工具核一遍联络线功率，偏差小于 0.5% 再往下走。

3.2 用 Heffron-Phillips 系数组装状态矩阵 A

K 常数可以直接由潮流工作点和网络导纳阵算出，但为了先把流程跑通，先给一组有代表性的单机模块数值。状态变量取四阶：Δω、Δδ、ΔEq′、ΔEfd。这个模块组装出来以后，多机系统就是把每台机的局部块摆到对角位置，再把跨机耦合项填到对应位置。

% 单机无穷大母线系统的 Heffron-Phillips 状态矩阵组装 % 状态顺序: [Delta_w; Delta_delta; Delta_Eqp; Delta_Efd] H = 6.5; % 惯性时间常数，s，基准容量 100 MVA D = 0; % 阻尼系数，pu，多数小信号算例取 0 omega0 = 2*pi*60; % 同步转速，rad/s，60 Hz 系统 K1 = 1.0; K2 = 0.8; K3 = 0.3; K4 = 1.4; K5 = -0.06; K6 = 0.5; Td0 = 7.5; % 励磁绕组开路时间常数，s TA = 0.05; % 励磁机等效时间常数，s KA = 50; % 励磁电压增益，典型范围 20~200 Efd0 = 2.5; % 稳态励磁电压，pu，平衡点处由潮流决定 A_local = zeros(4, 4); A_local(1,1) = -D/(2*H); % 阻尼项 A_local(1,2) = -K1/(2*H); % 同步转矩对功角 A_local(1,3) = -K2/(2*H); % 暂态电势产生的电磁转矩 A_local(2,1) = omega0; % 功角与转速的积分关系 A_local(3,2) = -K4/Td0; % 功角变化引起的励磁绕组磁链变化 A_local(3,3) = -1/(K3*Td0); % 励磁绕组自身衰减 A_local(4,3) = -KA*K6/TA; % 机端电压偏差经 AVR 反馈到励磁电压 A_local(4,4) = -1/TA; % 励磁机惯性

这段代码的物理含义要拆开看：第 1 行是摇摆方程，机械功率增量与电磁功率增量之差除以 2H 得到转速变化率；第 2 行是功角-转速的运动学关系；第 3 行描述励磁绕组暂态；第 4 行是自动电压调节器。需要强调的是，K5 在单机无穷大母线系统中通常为负，这正是高增益励磁带来负阻尼的根源，代码里它通过第 4 行的反馈路径影响整个闭环。

多机版本的扩展不用重写公式，只需要在 A 矩阵的跨块位置填入“第 i 台机电磁功率对第 j 台机功角的偏导”，数值由网络雅可比矩阵给出。规模变大后建议用稀疏矩阵组装，四机十二阶系统还好，扩展到几十台机时稀疏结构能省掉大量内存和计算时间。

3.3 特征值、振荡频率与阻尼比的提取脚本

拿到 A 矩阵后的第一件事就是全量特征值分解，然后从中筛出振荡模式。振荡模式的标志是特征值虚部显著不为零，工程上用一个很小的阈值过滤数值噪声。

% 特征值分解与振荡模式筛选 [V, D] = eig(A); lambda = diag(D); tol = 1e-6; % 虚部过滤阈值，低于此值视为非振荡模式 osc_idx = find(abs(imag(lambda)) > tol); lambda_osc = lambda(osc_idx); sigma = real(lambda_osc); freq_hz = abs(imag(lambda_osc)) / (2*pi); zeta = -sigma ./ sqrt(sigma.^2 + (2*pi*freq_hz).^2); % 按阻尼比从低到高排序，输出最危险的三个模式 [zeta_sorted, I] = sort(zeta); for k = 1:min(3, length(zeta_sorted)) fprintf('模式 %d: sigma=%.4f 1/s, f=%.3f Hz, zeta=%.4f\n', ... k, sigma(I(k)), freq_hz(I(k)), zeta_sorted(k)); end

这段脚本的关键参数是 tol。设得太大会漏掉低频区间模式，设得太小会把数值噪声当成振荡模式。工程上我习惯先用 1e-6 过滤，如果没有筛出任何模式再把阈值降到 1e-8，同时检查 A 矩阵是否存在病态。阻尼比排序是后续所有工作的起点：排在最前面的那几个模式，就是 PSS、储能、FACTS 设备应该去压制的对象。

大型系统里eig(A)会算全部特征值，算到几百阶矩阵时效率明显下降。此时改用eigs(A, 10, 'smallestabs')求模最小的十个特征值，往往更高效。需要提醒的是，不同 MATLAB 版本对位移参数的关键字有差异，老版本用'sm'，新版本推荐'smallestabs'，换版本跑脚本时容易在这里踩坑。

4. 参与因子：从一堆特征值里定位振荡源头

4.1 参与因子的计算：左右特征向量一起用

特征值只能告诉你振荡存不存在，不能告诉你在哪里。定位任务是参与因子的职责。参与因子 P_ki 由左特征向量和右特征向量对应分量相乘得到，它衡量第 i 个状态变量与第 k 个模式之间的关联强度。MATLAB 的eig返回的是右特征向量，左特征向量需要对 A 转置再求，或者直接用右特征向量矩阵的逆的共轭来得到。

% 从特征值分解结果计算参与因子矩阵 % V: 右特征向量矩阵，D: 特征值对角阵 [V, D] = eig(A); W = conj(inv(V)); % 左特征向量矩阵的便捷求法 Pmat = abs(W' * V); % 参与因子矩阵，行为模式，列为状态变量 % 查看摸个特定模式 k 下参与程度最高的状态量 k = 1; % 这里用阻尼比最差的模式为例 [P_sorted, idx] = sort(Pmat(k, :), 'descend'); for j = 1:3 fprintf('模式 %d 主导状态 #%d, 参与因子 %.4f\n', k, idx(j), P_sorted(j)); end

这里有个细节：W = conj(inv(V))得到的左特征向量方向可能与eig(A.')的结果相差一个常数因子，但参与因子对左右特征向量做了整体归一化，只要左右特征向量配对正确，最终数值是稳定的。脚本里的k是一个模式索引，对应第 3 章里排序后的危险模式编号。

4.2 用参与因子排序识别区间振荡模式与局部振荡模式

拿到参与因子矩阵后，识别逻辑就变得机械了。对某个振荡模式，如果参与度最高的状态来自多台机的功角和转速，且区域 1 的机组与区域 2 的机组相位正好相反，这就是区间振荡模式。如果参与度几乎集中在单台机上，就是局部振荡模式。频率也能辅助判断：低频 0.3~0.8 Hz 偏区间，1 Hz 以上偏局部。

实际项目中我习惯把参与因子结果与特征向量相位一起看。参与因子告诉你在哪个物理量上作用，特征向量相位告诉你参与的形式是同相还是反相。两者合成，才能决定 PSS 应该把信号加在哪里、加到什么方向。只盯着参与因子的数值大小，不看特征向量的符号，是最常见的误用方式。

4.3 PSS 接入点选择：看参与因子，不靠猜

PSS 安装位置的标准答案来自参与因子排序：哪个机组的功角与转速状态在目标模式中参与因子最大，PSS 就装在那台机组上。增益初值从 2~5 起步，然后通过根轨迹调整。所谓根轨迹在这里就是固定系统其他参数，只扫描 PSS 增益，观察目标特征值向左还是向右移动。

经验是：PSS 增益不是越大越好。超过某一点后，特征值反而向右拐，出现新的高频模式。这个拐点就是当前控制器结构的极限。实际调参时先把增益从 0 慢慢往上加，每步都重算区间模式阻尼比，看到阻尼比开始回落就停住。比直接拍一个大增益再试错要省一半时间。

5. 小信号稳定分析避坑与常见问题排障：翻车比较多的四类情况

5.1 平衡点初始化不对：正实部特征值未必代表系统失稳

现象：算出来的特征值里出现一个很大的正实部，比如 +0.8，但时域仿真曲线并没有立刻发散成那种夸张的速度，甚至某些变量还能稳住一小段时间。原因大概率是线性化点不在真正的潮流平衡点上，励磁初值、功角初值或负荷模型与潮流解对不上。解决：先把联络线功率与潮流设定值核对，误差控制在 0.5% 以内；再把 Eqp0 从估算值换成由 dq 分解得到的精确初值。多数翻车都是这一步偷懒造成的，别跳过。

5.2 单位基准没统一：振荡频率偏一倍的血泪经验

现象：区间振荡频率算出来约 1.2 Hz，而教科书和公开文献里 Kundur 两区域系统的区间模式大约是 0.6 Hz 左右，正好差一倍。原因几乎总是单位基准问题：发电机额定容量与系统公共基准容量不一致，惯性常数 H 没有按容量比折算；或者把 60 Hz 系统的 ω0 写成了 314。解决：所有发电机常数统一折算到公共基准，公式为 H_基准 = H_额定 × S_额定 / S_基准；ω0 直接用 2π×60，不要默认为 314。

5.3 特征值筛选漏掉目标模式：阈值太小看不见振荡

现象：设置了 tol = 1e-6，但筛出来的振荡模式全是局部模式，低频区间模式找不到。原因：区间模式对应的特征值虚部很小，数值上只有 0.05 左右，如果 A 矩阵的数值尺度偏大，这个模式很容易被当成实特征值。解决：先看一眼全部特征值的虚部分布，把 tol 从 1e-6 降到 1e-8 再筛；如果还是找不到，用eigs(A, 20, 'smallestabs')配合位移点近似求取，不要死磕全量eig。

5.4 励磁增益加大，阻尼反而变差：这是机理，不是代码 bug

现象：把励磁增益 KA 从 50 加到 400，区间振荡模式阻尼比从 0.06 掉到 0.02，甚至变成负的。初看像代码错误，其实是经典的小信号机理：高增益快速励磁通过 K5 路径引入负阻尼，抵消了同步电机的固有阻尼。解决：不要试图靠降低励磁增益来提升阻尼，那会牺牲电压调节动态性能。正确做法是加 PSS，或者对励磁系统设置合理的暂态增益降低环节。判断到底是机理还是 bug 的方法很简单：把 KA 调回 50，看阻尼比是否恢复，恢复就说明代码没问题，是控制参数交互在起作用。

6. 用根轨迹与时域校核给 12o3_v4 收尾

一套 kundur 模型值不值得投入时间，最终要看它的结论能不能经受双重校核。第一重是模态校核：把特征值分析得到的频率和阻尼比，与非线性时域仿真里振荡曲线的包络衰减对照。具体做法是给定一个小的功角阶跃扰动，记录功角差曲线，用对数包络拟合衰减斜率，看它是否与特征值实部一致。这个环节能戳破状态矩阵组装里隐藏的错误。

第二重校核是 PSS 参数扫描。以阻尼比为纵轴、PSS 增益为横轴画一条曲线，这条曲线的形状直接决定调参空间。扫参脚本很短：

Kpss_list = 0:0.5:15; zeta_result = zeros(size(Kpss_list)); for ii = 1:length(Kpss_list) A_pss = assemble_system_A(Kpss_list(ii)); % 重建含 PSS 的 A 矩阵 lambda = eig(A_pss); % 选择区间模式，计算阻尼比后存入 zeta_result(ii) end plot(Kpss_list, zeta_result, 'o-'); xlabel('PSS 增益'); ylabel('区间模式阻尼比');

我自己的习惯是，拿到这类 kundur 模型包先不急着信结果，第一步永远是把特征值重新算一遍、参与因子重新排一遍，确认前几步没被脚本的黑匣子糊弄过去。只有阻尼比、频率、参与因子三者能互相自洽，才会再去调控制参数。这个习惯救过我很多次。希望这些整理能帮到你。

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