简介一份基于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 矩阵的特征值决定了扰动后振荡的形态和衰减速度。输出量计算方法在小信号稳定中的作用特征值 λσ±jωeig(A)σ0 系统稳定σ0 临界σ0 失稳振荡频率 ffω/2π区分局部模式与区间模式阻尼比 ζ-σ/√(σ²ω²)工程判据低于 3% 视为弱阻尼参与因子 Pv_ki·w_ik判断振荡由哪些状态变量主导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 系数组装状态矩阵 AK 常数可以直接由潮流工作点和网络导纳阵算出但为了先把流程跑通先给一组有代表性的单机模块数值。状态变量取四阶Δω、Δδ、Δ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/s60 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 模型包先不急着信结果第一步永远是把特征值重新算一遍、参与因子重新排一遍确认前几步没被脚本的黑匣子糊弄过去。只有阻尼比、频率、参与因子三者能互相自洽才会再去调控制参数。这个习惯救过我很多次。希望这些整理能帮到你。本文还有配套的精品资源点击获取