简介面向电力系统小信号稳定性研究的MATLAB算例包基于Kundur经典两区域四机三模式模型适合电力专业高年级本科生、研究生及工程技术人员用于掌握小扰动下的系统动态分析与特征值判别方法。压缩包内共2个文件包含一个MATLAB脚本和一个说明文档脚本实现系统参数设定、状态矩阵构建、扰动施加及仿真结果输出说明文档则对模型结构与使用流程进行补充整体约6KB轻量便于快速运行与修改。目前已有211人浏览学习可作为课堂配套练习或自学仿真素材。通过运行脚本可复现不同运行模式下频率、功角与电压的动态响应曲线并结合特征值实部正负判断系统稳定裕度同时学习将非线性模型线性化处理的状态空间建模思路为后续电力系统小信号稳定分析提供完整的方法示范与可拓展的代码基础。1. 这串字符背后是电力系统低频振荡分析的通用试验台字符串里的 kundur 指电力系统小信号分析里绕不开的四机两区域测试系统简单说就是两台发电机模拟一个区域两个区域之间用一条联络线连起来用来复现真实的区间低频振荡场景。这个问题从几十年前研究到今天结论一直是靠这个系统验证的。后面的 matlab 和小信号稳定明示了这是一整套能直接跑的源码包从潮流计算、状态矩阵构建到特征值分析和参与因子计算全部围绕“振荡怎么来、怎么定位、怎么抑制”展开。如果你是刚接触小信号稳定分析的工程师或研究生想搞懂低频振荡和 PSS 参数整定但没有头绪拿这套系统的代码包入手是成本最低的路径。我在实际项目中多次用它来验证控制算法至少能帮你节省大量搭仿真环境的时间下面把拆包、跑通和踩坑的思路一次讲清楚。2. 认识 Kundur 系统小信号稳定分析绕不开的标准测试床2.1 系统拓扑与经典参数四台机、两个区域、一条联络线小信号稳定研究最头疼的问题不是算法而是没有一个“大家都认”的测试系统。你用自己的电网模型结果没法跟别人对比你用一个随手改的算例审稿人又不认。Kundur 系统的价值就在这里拓扑简单、参数公开、行为典型既有本地振荡模式又有区间振荡模式恰好覆盖小信号分析需要验证的两类核心现象。我拿到压缩包后做的第一件事不是急着跑脚本而是先看数据文件里母线、线路、发电机的定义是否与我手上的版本一致。这个系统存在 11 节点和 12 节点等多个版本你手上的 kundur_12o3_v4 按命名推断是某个 12 节点子版本。不同版本差别主要在负荷分配和联络线长度上不影响下面讲的建模思路。系统拓扑可以概括为区域 1 两台发电机带一个负荷节点群区域 2 同样两台发电机带一个负荷节点群两区域之间通过一条双回联络线连接。每个区域内部还留有本地负荷节点用来制造本地振荡模式和区间振荡模式的频率差。% kundur_sys.m —— 两区域四机系统数据模板以包内数据为准 % 母线数据: [母线号 区域号 电压幅值(pu) 电压相角(度) 有功负荷(pu) 无功负荷(pu)] bus_data [ 1 1 1.03 0.0 0.0 0.0; % 区域1机端母线G1 2 1 1.01 -5.0 9.4 1.5; % 区域1负荷母线L1 3 1 1.01 -7.0 9.4 1.5; % 区域1负荷母线L2 4 2 1.01 -8.0 0.0 0.0; % 区域2机端母线G4 5 2 1.01 -10.0 9.4 1.5; % 区域2负荷母线L3 6 2 1.01 -12.0 9.4 1.5; % 区域2负荷母线L4 7 1 1.02 -3.0 0.0 0.0; % 区域1内第二台发电机G2 8 2 1.02 -3.5 0.0 0.0; % 区域2内第二台发电机G3 ]; % 线路数据: [起点 终点 电阻(pu) 电抗(pu) 对地电纳一半(pu)] line_data [ 1 2 0.0025 0.025 0.043; % 区域1内部线 1 3 0.0010 0.010 0.017; 2 3 0.0025 0.025 0.043; 3 4 0.0010 0.010 0.017; % 联络线起点 4 5 0.0025 0.025 0.043; 4 6 0.0010 0.010 0.017; 5 6 0.0025 0.025 0.043; 3 5 0.0110 0.110 0.191; % 两区域间双回联络线 ];这段代码把系统抽象成两类数据母线决定节点类型和区域归属线路决定网络拓扑和电气距离。注意联络线母线 3 到 5的电抗 0.11 pu 明显大于区域内部线路的 0.025 pu这个电抗差是区间振荡模式产生的根源——电抗大意味着两个区域的功角耦合刚度弱同步转矩不足以把两群发电机拉在一起于是出现慢速的区间摇摆。参数说明坐标都是标幺值基准容量取全系统统一的 100 MVA所以任何发电机参数在写入状态矩阵前都要先换算到这个基准下。母线 1 设置为平衡节点承揽系统的有功差额。电压相角初始值不是潮流计算的结果而是潮流迭代的初值真正求解时会覆盖这里示意用。2.2 状态方程里藏着的机电振荡本质小信号稳定的核心是判断平衡点的局部渐近稳定性。这句话翻译成人话就是系统受到一个小扰动后是逐渐回到原来的运行点还是越摆越大最终失步。回到 Kundur 系统四台发电机就是四个摆动体它们通过电网电磁功率耦合在一起构成一个多自由度摆系统。经典二阶模型下每台发电机的转子运动方程线性化后写成Δδ̇ ω0 Δω2H · Δω̇ ΔTm − ΔTe − D · Δω其中 δ 是功角ω 是角速度偏差H 是惯性时间常数D 是机械阻尼系数ΔTe 是电磁功率变化量。电磁功率本身是 δ 的非线性函数在平衡点附近泰勒展开后ΔTe 可以写成一堆 Δδ 的线性组合。把这些方程和励磁系统动态方程拼在一起就得到标准的线性状态方程Δẋ A · ΔxA 矩阵的特征值决定了系统的动态行为实部为负的所有特征值系统小扰动稳定实部为正的特征值系统存在增幅振荡或单调失稳模式。特征值的虚部对应振荡频率实部和虚部一起算出阻尼比。Kundur 系统最吸引人的地方在于它的 A 矩阵特征值里天生就有两组典型模式——区域之间约 0.5 Hz 的区间振荡模式和区域内部约 1.1 Hz 的本地振荡模式跟真实电网中观测到的现象高度一致。2.3 用 MATLAB 组织系统数据从散乱数组到结构体数据文件在压缩包里可能是 .m、.mat 或者 .txt 格式拿到后别急着改数值先读一遍结构。我最开始的习惯是把数据一股脑塞进工作区变量后来发现改一个参数要重新跑整个脚本这种工作方式在小信号分析里极容易翻车还是老老实实换成结构体承载全部系统信息。% load_kundur.m —— 将原始数据整理为统一结构体 sys struct(); sys.bus bus_data; % 母线数据 sys.line line_data; % 线路参数 sys.baseMVA 100; % 基准容量 sys.gen [ 7 900 13.8; % G1: 母线号额定容量MVA机端电压kV 8 900 13.8; % G2: 同上 1 900 13.8; % G3: 同上 4 900 13.8; % G4: 同上 ]; sys.freq 60; % 系统额定频率 Hz sys.omega0 2*pi*60; % 同步角速度 rad/s sys.genData.H [6.5 6.5 6.5 6.5]; % 惯性时间常数 s sys.genData.Ra [0 0 0 0 ]; % 电枢电阻 pu sys.genData.xd [1.8 1.8 1.8 1.8]; % 直轴同步电抗 sys.genData.xq [1.7 1.7 1.7 1.7]; % 交轴同步电抗 sys.genData.xdp [0.3 0.3 0.3 0.3]; % 直轴暂态电抗 sys.genData.Td0p [8.0 8.0 8.0 8.0]; % 直轴开路时间常数 s结构体方案的好处是参数传递不会错位。四台机组参数一致是 Kundur 系统的惯例这也让振荡模式的分析更干净因为对称性模式天然解耦参与因子不会因为参数不对称而糊成一片。后面构建状态矩阵和 PSS 设计时只需要把这个 sys 对象传进函数改参数只需要动一处。参数说明H 取 6.5 s 是原系统参考值表示大机组惯量偏大的特征。Ra 是定子电阻在暂态模型中通常忽略所以留 0。xd 和 xq 是稳态同步电抗决定机电功率传输的最大能力xdp 是暂态电抗决定暂态过程中发电机端电压与内部电势的耦合强度。Td0p 是励磁绕组时间常数和励磁系统动态一起影响振荡模式的阻尼。这些值在实际包内可能有差异以你手上的 m 文件为准。3. 从潮流到特征值把小信号分析全流程跑通3.1 先算潮流状态矩阵的初始化前提小信号分析的最大前提是“平衡点”。线性化方程里的 A 矩阵每一项都在某个稳态运行点取值而这个运行点由潮流解确定。很多人直接从状态方程出发写程序不先算潮流结果 A 矩阵里的网络电压和功角初值全是猜的算出来的特征值自然不对。我一般把潮流计算看成两个任务一是算各母线的电压幅值和相角二是从电压算出发电机端电流、内部电势和励磁电压初值。后者在状态矩阵构建里特别关键做不好后面全是错。% powerflow_kundur.m —— 用牛顿-拉夫逊法求平衡点 % 输入: sys结构体; 输出: x0为潮流解向量, V/theta为母线电压 function [V, theta, x0] powerflow_kundur(sys) nbus size(sys.bus, 1); V sys.bus(:, 4); theta deg2rad(sys.bus(:, 5)); % 节点类型: 1平衡 2PQ 3PV这里按母线数据第三列约定 type sys.bus(:, 3); Psp sys.bus(:, 6); % 有功注入给定值 Qsp sys.bus(:, 7); % 无功注入给定值 % 生成导纳矩阵 Y Y ybus_from_lines(sys.line, nbus); % 牛顿法迭代误差收敛到1e-9 maxIter 30; tol 1e-9; for iter 1:maxIter [P, Q] power_injection(V, theta, Y); mismatch [Psp - P; Qsp - Q]; % 剔除平衡节点的无功方程和电压幅值方程得到指标准则 mismatch(type1) 0; % 平衡节点不参与误差计算 mismatch(nbusfind(type1)) 0; if max(abs(mismatch)) tol, break; end J jacobian_power(V, theta, Y); dX J \ mismatch; theta theta dX(1:nbus); V V dX(nbus1:end); end x0 [theta; V]; % 平衡点向量 end这里我没把雅可比矩阵函数贴全它的结构是标准牛顿法对角元是各母线功率对本身电压相角/幅值的偏导数非对角元是互偏导。潮流收敛后V 和 theta 就是初始平衡点。从工程角度讲如果潮流不收敛首先检查线路数据和负荷单位是不是统一其次检查 PV 节点的无功越限处理。参数说明tol 取 1e-9 是为了后续特征值计算的数值稳定性因为对平衡点的微小偏差会被状态矩阵放大maxIter 取 30 对四机系统绰绰有余。J 矩阵在 MATLAB 里直接用稀疏矩阵存储会快很多四机系统节点少全矩阵也无所谓。潮流求解完成后紧接着要做发电机初值计算反算每台发电机的内部暂态电势 Eq、Ed 和励磁电压 Efd。这里的物理意义是发电机端电压和电流从潮流结果里提取后通过相量方程求出内部状态励磁电压初始值则决定了励磁系统稳态输出。没有这一步后面构建状态方程时 A 矩阵的初始化全部失真。3.2 状态矩阵构建的核心逻辑与代码状态矩阵构建有两种路线。一是手推线性化表达式把转子运动方程、励磁动态方程和网络方程联立消去代数变量后得到显式 A 矩阵二是构建微分代数方程组 DAE然后用数值方法消去代数变量。手推路线推导量特别大六阶模型下很容易出错。我常用的是数值雅可比路线先写出微分方程右端函数 f_sys(x)其中代数变量在每一步通过潮流方程即时求解然后对状态向量做小扰动用差分拟合出 A 矩阵。% build_A_matrix.m —— 数值雅可比法构建状态矩阵 % 状态向量 x: 每台发电机5个状态20状态系统 % 顺序: [delta_i; omega_i; Eq_i; Ed_i; Efd_i], i1..4 function A build_A_matrix(x0, sys) nx length(x0); eps_diff 1e-6; A zeros(nx, nx); for j 1:nx xp x0; xm x0; xp(j) xp(j) eps_diff; xm(j) xm(j) - eps_diff; % f_sys: 计算各发电机状态变量的导数 fxp f_sys(xp, sys); fxm f_sys(xm, sys); A(:, j) (fxp - fxm) / (2 * eps_diff); end % 输出A之前检查是否含NaN这是系统未收敛的信号 if any(isnan(A(:))) error(A矩阵存在NaN检查潮流解或初始状态); end end function dx f_sys(x, sys) % 解析状态向量 ng 4; nx length(x); states_per nx/ng; dx zeros(size(x)); for k 1:ng delta x((k-1)*states_per 1); omega x((k-1)*states_per 2); Eq x((k-1)*states_per 3); Ed x((k-1)*states_per 4); Efd x((k-1)*states_per 5); % 由内部电势和网络方程求解端口电气量 [Pe, Qe, Id, Iq, Vt] solve_network_eq(delta, Eq, Ed, k, x, sys); % 转子运动方程 dx((k-1)*states_per 1) sys.omega0 * omega; dx((k-1)*states_per 2) (1/(2*sys.genData.H(k))) * ... (sys.genData.Tm(k) - Pe - sys.genData.D(k)*omega); % 励磁绕组动态 dx((k-1)*states_per 3) (1/sys.genData.Td0p(k)) * ... (Efd - Eq - (sys.genData.xd(k) - sys.genData.xdp(k))*Id); % 励磁系统一阶惯性模型: 输出电压跟随参考值 dx((k-1)*states_per 5) (1/sys.exc.Tr(k)) * ... (sys.exc.Ka(k)*(sys.exc.Vref(k) - Vt) - Efd); end end这段代码里 f_sys 的输入是全部发电机的状态向量输出是每个状态的导数。solve_network_eq 函数在每次调用时要解一个简化的网络方程把内部电势和功角映射到端口有功、无功和端电压。在 20 状态的小系统上这样嵌套求解的速度完全够用整段构建 A 矩阵跑几十次数值雅可比耗时不超过几秒。参数说明eps_diff 取 1e-6 是数值差分经验值。取得太大比如 1e-3差分截断误差会把低频模式的阻尼比算错取得太小比如 1e-10浮点舍入误差会主导结果。如果发现特征值虚部全是零或出现大量伪模式第一件事就检查这个差分步长。状态变量顺序必须和后面参与因子分析保持完全一致最好在建状态矩阵时同时生成一个状态名称列表后面做模态识别时对照使用。3.3 特征值分析与模式筛选把振荡模式从特征值海里捞出来A 矩阵建好后直接调用 MATLAB 自带的 eig 就能得到全部特征值。四机五阶模型加励磁系统总共 20 个特征值其中大部分是励磁系统和网络动态的快模研究低频振荡只需要盯着机电模式。识别方法是看频率范围电力系统低频振荡一般落在 0.1~2.0 Hz 之间对应的虚部在 0.6~12.6 rad/s。另外振荡模式必然是成对的共轭复数实特征值对应非振荡衰减模式。% mode_identify.m —— 特征值计算、排序、筛选振荡模式 [V, D] eig(A); lambda diag(D); % 计算每个模式的频率(Hz)和阻尼比 freq abs(imag(lambda)) / (2*pi); zeta -real(lambda) ./ abs(lambda); % 锁定机电振荡模式: 频率在0.1~2Hz范围的共轭对 candidate_idx find(freq 0.1 freq 2.0 abs(imag(lambda)) 1e-6); % 按阻尼比从低到高排序先看最危险的模式 [zeta_sorted, sort_idx] sort(zeta(candidate_idx)); freq_sorted freq(candidate_idx(sort_idx)); % 打印模式识别结果 for i 1:length(freq_sorted) fprintf(模式%d: 频率%.3f Hz, 阻尼比%.4f\n, ... i, freq_sorted(i), zeta_sorted(i)); end运行这段脚本典型输出会看到大约 0.54 Hz、阻尼比在 0.05 左右的区间振荡模式以及 1.1 Hz 附近的本地振荡模式。阻尼比为正说明小扰动稳定但如果阻尼比低于 0.03工程上就认为这个模式处于危险区必须加 PSS 或调整运行方式。如果算出来阻尼比是负的说明该模式增幅振荡现场叫“负阻尼”这就是要治理的目标。我在实际调试时一般会额外打印特征值的实部绝对值辅助判断模式的可辨识度。比 0.001 还小的实部意味着该模式几乎落在虚轴上数值误差可能让它时正时负这种模式计算结果的工程意义要谨慎对待不能直接下结论说不稳定。出现这种情况优先怀疑不是算法问题而是平衡点没算准。4. 模式识别与参与因子从“有振荡”到“谁在振荡”4.1 参与因子计算的正确姿势特征值告诉你系统有哪些振荡模式但没告诉你每台发电机在这些模式里参与了多少。参与因子就是从 A 矩阵的左、右特征向量推导出的量化指标它表示某个状态变量在某个模式里的相对参与程度。参与因子大的状态变量就是该模式的主导变量顺着它就能定位到具体机组和物理环节。参与因子的计算公式是左特征向量元素与右特征向量对应元素的乘积再归一化。这里有一个常见的隐蔽坑MATLAB 的 eig 函数返回的 V 是右特征向量矩阵左特征向量矩阵不是 inv(V)而应该是 inv(V).。写错这个转置参与因子的结果会完全对不上。% participation_factor.m —— 参与因子计算与可视化 % 输入: A矩阵, 已排序特征值lambda_sorted, 对应右特征向量V_sorted % 输出: 每个状态在指定模式中的参与度 function PF participation_factor(V_sorted, lambda_sorted, mode_id) % 左特征向量矩阵注意是inv(V).不是inv(V) W inv(V_sorted).; % 对该模式取左、右特征向量 right_vec V_sorted(:, mode_id); left_vec W(:, mode_id); % 参与因子 右特征向量元素 .* 左特征向量元素 PF right_vec .* left_vec; % 取模并归一化方便对比不同状态的参与程度 PF abs(PF); PF PF / sum(PF); end参数说明mode_id 是特征值排序后当前要分析的模式序号对应上一节筛选出的振荡模式。PF 向量长度等于状态数量在 Kundur 系统上就是 20 个元素。归一化后的 PF 可以跨模式比较。要注意的是参与因子对状态变量的缩放很敏感如果某个状态在状态方程里以角度为单位、另一个以标幺值 pu 为单位数值量级差异会影响 PF因此最好在构建 A 矩阵之前统一各状态变量到工程惯用单位。实际项目里我会把 20 个状态的 PF 画成柱状图横轴标注状态名纵轴是参与度。如果区间振荡模式对应的 PF 在四台发电机的 delta 状态上都明显偏大说明所有机组都在参与这是典型的区域间振荡特征如果 PF 集中在一两台发电机的 omega 和 Efd 状态上则说明是本地模式的励磁回路参与占主导。4.2 振荡模式与机组机电回路关联读懂参与因子的物理含义参与因子算出来只是第一步关键是读懂它。Kundur 系统因为结构对称分析结果非常干净区间振荡模式的参与因子在区域 1 的两台发电机和区域 2 的两台发电机之间呈“对顶”分布——区域 1 的机组相位一致区域 2 的机组相位一致两个区域之间相位相反。这种分布从参与因子上就能看出来同一个状态变量在两组机组上的 PF 值都大但因为特征向量对应元素的符号相反功角变化方向相反。“谁在振荡”的物理解读要落到机电回路上。我一般先看 delta 状态和 omega 状态的参与度这两者代表转子运动的参与强度再看 Efd 状态的参与度它代表励磁系统对模式的影响。Efd 参与度高意味着励磁放大倍数对该模式阻尼有显著影响这是后面设计 PSS 时的直接依据。还有一种常见误读是把参与因子当成灵敏度。参与因子大不代表调整这个状态对应的物理参数就一定能有效改变特征值这两个指标相关性存在但不是一回事。真正决定参数整定方向的是特征值灵敏度也就是特征值对某参数的偏导数这个指标要通过参数扰动或解析求导得到判断哪些参数是阻尼控制的有效杠杆。% sensitivity_scan.m —— 特征值对励磁增益Ka的灵敏度粗略扫描 % 通过小范围参数扰动观察目标模式阻尼比变化趋势 % 这里用控制变量法只改一台机的Ka扫典型值范围 Ka_range 20:5:200; Zeta_target zeros(length(Ka_range), 1); lambda_store zeros(length(Ka_range), 1); for i 1:length(Ka_range) sys_copy sys; sys_copy.exc.Ka(1) Ka_range(i); % 只动G1的励磁增益 A_new build_A_matrix(x0, sys_copy); [~, D_new] eig(A_new); lambda_new diag(D_new); % 取最接近目标模式频率的那对共轭特征值 omega_target 2*pi*0.54; [~, idx_match] min(abs(abs(imag(lambda_new)) - omega_target)); lambda_sel lambda_new(idx_match); lambda_store(i) lambda_sel; Zeta_target(i) -real(lambda_sel)/abs(lambda_sel); end plot(Ka_range, Zeta_target); xlabel(励磁增益 Ka (pu/pu)); ylabel(区间模式阻尼比);这段脚本能直观看到励磁增益升高会压缩阻尼比这是设计 PSS 前最重要的认识高放大倍数励磁在提升暂态稳定性的同时会牺牲小信号阻尼。阜实到 Kundur 系统上把某台机的 Ka 从 50 提到 150通常能看到区间模式阻尼比从 0.05 掉到 0.02 以下。这个敏感性结论也解释了为什么单纯调励磁参数很难同时满足暂态稳定和小信号稳定两方面的需求必须靠 PSS 提供附加阻尼。5. 避坑指南小信号分析落地中最容易翻车的五个细节5.1 现象特征值虚部全为零振荡模式一个都不剩第一回跑通代码时我盯着终端愣住了20 个特征值全部是实数一个共轭对都没有低频振荡分析无从谈起。原以为是系统的平衡点本来就是单调收敛的后来发现是转子运动方程里的 ω 状态被当成了常数。原因状态变量表里没有把 ω 纳入微分方程或者 f_sys 里转子方程用了稳态近似导致虚拟的“角速度”状态实际是静态变量整个系统退化成纯网络动态振荡自由度消失。解决检查状态清单里每台机是否包含 delta 和 omega 两个状态并在 f_sys 里确认 dx(ω) 的表达式是 (Tm − Pe − D·ω)/(2H)而不是直接把 ω 设成 0。把 ω 状态拿掉是二阶模型简化时常见的误操作在 Kundur 系统这种四机模型上一定不能省。5.2 现象重负荷工况下特征值判系统不稳定但非线性仿真却表现为稳定潮流解的平衡点本身是多解问题。重负荷下存在低压解和高压解两套平衡点潮流计算可能收敛到低压解。低压解对应的系统运行点功角远特征值计算出来常常是负阻尼。但实际非线性仿真是从一个更合理的运行点出发的所以两者结论冲突。原因潮流初值给得太差牛顿法收敛到了错误的平衡点分支。解决改变初值策略用平启动电压 1.0 pu、相角 0 度再做一次潮流或者用连续潮流沿着负荷增长路径追踪平衡点。工程上最简单的验证方式是把潮流解里的最低电压挑出来如果低于 0.85 pu就高度怀疑落入低压解分支重新设置初值。5.3 现象参与因子指向的“问题机组”与机电模态图对不上算出来的 PF 明明显示 G1 的主导度最高可画功率振荡曲线却发现振荡能量集中在 G3 和 G4 之间跟参与因子的结论完全相反。原因多数情况是左特征向量取错了。参与因子里要用的左特征向量是右特征向量矩阵逆的转置很多人只写了 inv(V)没有做转置。用错之后 PF 的数值就失去了物理意义。解决参与因子计算统一用 W inv(V).不要记“左特征向量等于逆矩阵”这种简化说法。同时用模态图交叉验证画出每个状态在该模式下的相对幅值和相位。正确的参与因子和模态图指向的机组必然一致不一致时优先检查特征向量矩阵的构造。5.4 现象调整发电机模型阶数后结果发生断崖式变化从三阶模型切换到五阶模型关键模式从阻尼比 0.08 变成了 −0.01振荡频率也漂了小半赫兹。一开始以为是模型阶数差异导致的正常现象后来发现是数据文件里五阶参数缺失对应的时间常数被 MATLAB 自动填成 0系统动态彻底变了。原因不同模型阶数需要不同组的参数xq、Tq0 这些参数在低阶模型里用不到高阶模型一旦引用到缺省值为 0 会让某些动态环节变成纯积分器A 矩阵的模态结构自然全面改变。解决切换模型阶数时把参数表也一起切换。用函数接口传参而不是散落在工作区的变量这样每次换模型必然换数据。我习惯在 sys 结构体里放一个 model_order 字段build_A_matrix 内部根据它的值决定调用哪套参数从机制上杜绝参数错配。5.5 现象自建代码计算结果和商用仿真软件相差超过 10%同一套 Kundur 系统自建脚本算出的区间振荡频率 0.54 Hz商用软件给出 0.59 Hz阻尼比差异也超出预期。排查了很久才发现问题不在算法在基准值。原因发电机参数如果以每台机额定容量为基准给出的而整个网络潮流方程以系统基准容量计算二者混用不做换算电磁功率和励磁环节的标幺值就全部失真。解决所有发电机阻抗、时间常数、增益在进入状态矩阵前统一用 S_gen / S_base 换算到系统基准容量下。换算规则是电抗的标幺值正比于容量时间的标幺值不变增益系数按容量比缩放。把这条换算写进 load_kundur.m不要靠手算。6. 从分析到治理用特征值结果反推 PSS 参数并验证6.1 谐振与相位补偿PSS 参数整定的底层逻辑特征值分析最终要为控制服务。Kundur 系统最经典的治理手段是电力系统稳定器 PSS它在励磁电压参考点叠加一个附加信号为低频振荡模式提供正阻尼。PSS 的核心设计任务是两件事相位补偿和增益整定。先看相位补偿的原理。励磁系统从电压偏差到电磁转矩之间存在相位滞后在振荡频率处可能接近 90 度甚至更多。PSS 的结构通常是一个隔直环节加两阶超前-滞后补偿隔直环节让稳态下 PSS 无输出超前-滞后环节在目标振荡频率处提供与励磁滞后相反的相位。目标频率在 Kundur 系统上就是区间振荡模式的大约 0.54 Hz。6.2 用特征值灵敏度扫描确定 PSS 增益相位环节设计好后剩下就是增益整定。增益太小阻尼提升不明显增益太大可能激发其它模式振铃甚至导致励磁限幅。我实践中用特征值的增益扫描代替现场试凑效率高很多。% pss_tuning.m —— 遍历PSS增益观察目标模式阻尼变化 % PSS传函: Kpss * sTw/(1sTw) * (1sT1)/(1sT2) * (1sT3)/(1sT4) Tw 10; % 隔直时间常数 T1 0.05; T2 0.02; % 超前-滞后环节 T3 0.05; T4 0.02; Kpss_range 0:2:50; zeta_best 0; K_best 0; for Kpss Kpss_range % 把PSS动态并入状态矩阵新增一个PSS状态变量 sys_pss sys; sys_pss.pss.K Kpss; sys_pss.pss.Tw Tw; sys_pss.pss.T1 T1; sys_pss.pss.T2 T2; sys_pss.pss.T3 T3; sys_pss.pss.T4 T4; A_pss build_A_matrix_with_pss(x0, sys_pss); lambda_pss eig(A_pss); % 追踪0.54Hz附近模式的阻尼比 freq_pss abs(imag(lambda_pss))/(2*pi); target_idx find(freq_pss 0.4 freq_pss 0.7); zeta_target max(-real(lambda_pss(target_idx))./abs(lambda_pss(target_idx))); if zeta_target zeta_best zeta_best zeta_target; K_best Kpss; end end fprintf(最佳PSS增益: Kpss%.1f, 阻尼比提升至%.4f\n, K_best, zeta_best);参数说明T1 到 T4 的取值来自相位补偿计算这里 0.05 和 0.02 的组合在 0.5 Hz 附近提供大约 25 度的超前角具体值要结合励磁模型相位滞后曲线微调。Tw 取 10 秒是为了让 PSS 只响应 0.1 Hz 以上的振荡不干扰稳态调压。Kpss 的搜索范围根据系统容量先给一个粗扫找到阻尼比峰值后再在峰值附近细扫避免粗扫间隔错过最佳点。6.3 时域验证用故障仿真检验阻尼效果特征值分析最终要回到时域验证。我惯用的验证手段是在区域间联络线中部设置一个三相短路故障100 毫秒后切除然后观察区域间传输功率的振荡曲线。没有 PSS 时功率曲线以大约 0.54 Hz 缓慢衰减甚至等幅振荡投入 PSS 后振荡应该明显衰减三五秒内趋于平稳。验证脚本的核心就是调用了前面所有函数先潮流求平衡点再分别构建无 PSS 和有 PSS 的状态矩阵对两种场景做时域数值积分。时域积分可以直接用 MATLAB 的 ode45但注意这是刚性系统我一般改用 ode23t 或直接按步长离散求解避免刚度问题把振荡曲线算得面目全非。这一整套流程走完你不仅拿到了 Kundur 系统的特征值和振荡模式还完成了一次完整的 PSS 设计与验证闭环。这里我想说句实在话小信号分析的本质是建立“参数变化→特征值移动→阻尼变化”的因果链很多人卡在中间某个环节就放弃了但真正做通了这一套再回去看自己手头项目的低频振荡问题视角会完全不一样。这套系统我前前后后跑了不下几十遍每遍都能找到新的参数细节值得琢磨。希望这点文字能帮你在同样的路上少绕几个弯直接站到能改代码、能调参数、能解释现象的位置上。本文还有配套的精品资源点击获取