简介基于特征值分析法的电力系统稳定性研究是一份面向电力系统研究人员、电气工程师及高校相关专业学生的技术文献聚焦次同步谐振问题的分析判断与仿真验证。资源系统阐述了特征值分析法的理论优势、模型建立过程以及次同步谐振产生机理并对发电机电磁回路、汽轮机和SVC等环节的线性化模型进行了推导展示结合具体算例给出SSO特性分析与对比结论。整包共1个PDF文件压缩包大小约244KB内容精炼属于参考文献与专业指导类资料。已有206人学习这份资源适合需要借助严格线性化方法判断串补输电系统稳定性、开展振荡事故机理分析的读者。阅读后可系统掌握特征值分析法在电力系统振荡研究中的完整应用链路为稳定性建模与仿真工作提供理论支撑。1. 特征值分析法为什么是电力系统稳定性研究的主力工具电力系统最令人头疼的运行状态不是瞬时短路而是“看着还在运行、其实正在失步”一台发电机组的功角在 0.8 Hz 附近小幅往复摆动幅值一次比一次大十几秒后保护动作跳闸运行人员事后复盘时往往只能看到功率波动和电压抖动这两个表象。而小信号稳定分析在故障发生之前就已经给出了答案——系统状态矩阵的某个特征值实部从负变正。特征值分析法就是把“看不见的负阻尼”变成一个可计算、可排序、可优化的数学量这正是它在电力系统稳定性研究里沿用数十年、至今仍是工业级校核标准的原因。这篇内容面向做并网仿真、机组控制参数设计与稳定性校核的工程师和研究生覆盖从单机无穷大系统到多机系统的状态矩阵搭建、参与因子定位振荡源、PSS 参数整定到最终校核的完整管线。需要先说明边界这里讨论的是小信号稳定性即系统在某个运行点附近受小扰动后的动态行为大扰动暂态稳定不在本篇范围内。2. 从 DAE 到状态矩阵特征值分析法的建模与计算管线2.1 电力系统模型为什么写成微分代数方程电力系统的完整动态描述是一组微分代数方程Differential Algebraic EquationsDAE结构为dx/dt f(x, y, u) 0 g(x, y)其中 x 是动态状态变量包括发电机功角 δ、转速 ω、暂态电动势 Eq、励磁电压 Efd、调速器阀门开度等y 是代数变量主要是节点电压幅值与相角由潮流方程约束。g(x, y) 0 代表网络方程它没有导数项因为电磁暂态比机电暂态快得多工程上将其视为瞬时平衡。小信号分析的核心步骤是在运行点 (x₀, y₀) 处做一阶泰勒展开把 f 和 g 线性化然后消去代数变量 Δy最终得到dΔx/dt A · Δx这个 A 就叫状态矩阵维度等于动态状态变量的个数。特征值分析法的所有结论都来自对 A 的特征分解。需要强调的是线性化之后的结论只在运行点附近成立运行点变了A 就变特征值也要重新算。常见的做法是对一批典型运行方式分别计算而不是只算一个工况。2.2 单机无穷大系统的 Phillips-Heffron 模型要理解 A 矩阵的物理来源最经典的载体是单机无穷大SMIB系统的 Phillips-Heffron 模型。它把一台同步发电机、励磁系统和无穷大母线等效为一个四阶线性系统状态变量取 Δδ、Δω、ΔEq、ΔEfd。模型中的六个常数 K1~K6 是从潮流解和电机参数中推导出来的每个常数都有明确的物理含义。常数物理含义典型数值示例工况K1同步转矩系数功角变化引起的电磁转矩变化0.5 ~ 1.5K2暂态电动势变化引起的电磁转矩变化0.8 ~ 1.8K3与 d 轴暂态电抗相关的阻抗系数0.2 ~ 0.5K4去磁效应系数功角变化引起的励磁绕组去磁1.0 ~ 2.5K5电压变化对功角的敏感度-0.2 ~ 0.2K6电压变化对暂态电动势的敏感度0.3 ~ 0.7注意 K5 的符号在不同教材里可能有差异因为它代表负荷和网络结构对电压支撑特性的综合影响系统重负荷时 K5 为负的情况在工程中很常见而这正是负阻尼振荡的根源之一。建模时务必以自己采用的电机模型和基准方向为准不要跨文献混用符号约定。2.3 组装状态矩阵并计算特征值的 MATLAB 实现拿到 K1~K6 之后状态矩阵 A 可以直接按模型结构组装。下面是一段可运行的最小 MATLAB 脚本展示从参数到特征值的完整计算%% SMIB 四阶 Phillips-Heffron 模型的特征值计算 % 参数来自某典型重负荷运行点单位均为标幺值或秒 K1 1.0755; K2 1.2578; K3 0.3072; K4 1.7124; K5 0.0405; K6 0.4971; M 10.0; % 惯性时间常数按 2H 计算单位秒 D 0.02; % 机械阻尼系数 Tdo 6.0; % d 轴开路暂态时间常数单位秒 KA 200; % 励磁系统增益 TA 0.05; % 励磁系统时间常数单位秒 ws 314.159; % 同步角速度50 Hz 系统 % 状态变量顺序: [Δδ, Δω, ΔEq, ΔEfd] A [0 ws 0 0; -K1/M -D/M -K2/M 0; -K4/Tdo 0 -1/(K3*Tdo) 1/Tdo; -KA*K5/TA 0 -KA*K6/TA -1/TA]; lambda eig(A); % 求全部特征值 fprintf(特征值实部 虚部 频率(Hz) 阻尼比\n); for i 1:length(lambda) sigma real(lambda(i)); omega imag(lambda(i)); f omega / (2*pi); zeta -sigma / sqrt(sigma^2 omega^2); fprintf(%12.4f %12.4f %10.4f %10.4f\n, sigma, omega, f, zeta); end代码的逻辑并不复杂第一行到第六行是模型参数K1~K6 取一组能跑通教学示例的数值状态矩阵 A 的每一行对应一个状态变量的导数方程比如第一行表示 Δδ 的导数等于 ws 乘以 Δω这是功角与转速的固有运动关系第四行把励磁电压的动态和 Δδ、ΔEq 的反馈耦合起来。eig(A)一次返回全部特征值然后用实部、虚部换算频率和阻尼比。参数说明M 取 10 秒大约是 200 MW 级汽轮发电机的典型值TA 取 0.05 秒代表现代快速励磁系统KA 取 200 说明励磁增益较高这种参数组合下系统通常会出现弱阻尼的机电振荡模式。如果换成一个励磁增益只有 50 的旧式励磁系统特征值实部会明显更负这就是后面 PSS 整定章节的伏笔。2.4 特征值怎么读实部、虚部、阻尼比与频率特征值是一个复数 λ σ jω工程上不直接看实部和虚部而是换成阻尼比和振荡频率。阻尼比的定义是 ζ -σ / √(σ² ω²)。σ 为正时特征值在复平面右半平面系统不稳定任何小扰动都会发散σ 为负但绝对值很小时系统虽然稳定动态过程却拖得很长运行人员看到的就是“功率来回摆、半天停不下来”。模式类型频率范围主要参与对象工程关注度区间振荡模式0.1 ~ 0.3 Hz区域间联络线上的功率摆动极高本地机电模式0.4 ~ 0.8 Hz单台机组或一个厂站的功角摆动高控制相关模式1 ~ 2 Hz励磁、调速器与机组相互耦合中次同步模式10 ~ 50 Hz轴系扭振与串联补偿线路视场景而定判断标准方面国内工程校核一般要求机电模式阻尼比不低于 0.03主流设计规范更倾向 0.05 以上弱阻尼阻尼比小于 0.03的区域间模式需要重点关注。另有一个高频陷阱特征值实部为负但接近零时线性化模型说“稳定”实际系统中因为噪声和参数漂移可能已经临界振荡所以留足裕度是工程常规做法。3. 参与因子与模态识别在几十台机组里定位振荡源3.1 特征值告诉你有振荡但没告诉振荡在哪里多机系统的状态矩阵维度动辄几十甚至上百特征值求出来之后你只会面对一串复数某个 0.5 Hz 的弱阻尼模式确实存在但它到底是哪台机组和哪台机组在对着摆主导状态变量是功角还是励磁电压这些问题靠特征值本身回答不了需要引入参与因子Participation Factor。参与因子的数学定义是 P_ki v_ki · w_ik其中 v_ki 是第 k 个特征值对应右特征向量的第 i 个分量w_ik 是左特征向量对应分量。左右特征向量共同起作用的意义在于右特征向量衡量该模态在这个状态上“可以被看到”的程度左特征向量衡量外部激励对这个模态的“激发能力”两者相乘才代表该状态与模态之间的总体关联强度。如果只取其中一个会漏掉不可观或不可控的情况。3.2 用左右特征向量计算参与因子的 MATLAB 实现MATLAB 的 eig 函数支持同时返回左右特征向量代码路径很直接%% 计算参与因子矩阵 [V, D, W] eig(A); n size(A, 1); P zeros(n, n); for k 1:n for i 1:n % W 的列是左特征向量V 的列是右特征向量 P(i, k) abs(V(i, k) * W(i, k)); end % 每列归一化使参与因子总和为 1 P(:, k) P(:, k) / sum(P(:, k)); end lambda diag(D); % 假设我们关心频率在 0.5 Hz 附近的那个机电模式 target_f 0.5; [~, idx] min(abs(imag(lambda) / (2*pi) - target_f)); % 对该模态的参与因子排序显示前 5 个主导状态 [sortP, sortIdx] sort(P(:, idx), descend); fprintf(模态频率 %.2f Hz 的主导参与因子:\n, imag(lambda(idx)) / (2*pi)); for i 1:min(5, n) fprintf( 状态 %d: %.4f\n, sortIdx(i), sortP(i)); end逻辑说明[V, D, W] eig(A)返回的 W 列是左特征向量满足 W·A D·W与右特征向量之间具备双正交性因此参与因子的绝对值可以直接用对应分量相乘得到。归一化不是必须的但把所有值缩放到总和为 1 之后更容易横向比较不同模态之间的主导程度。参数说明target_f要根据实际系统选择单机系统看 0.4~0.8 Hz 的本地模式区域联网系统还要额外检查 0.1~0.3 Hz。排序后如果某个模态的参与因子集中在少数几个状态上说明模态局部化程度高治理手段也相对明确。3.3 多机系统里怎么用参与因子选机组参与因子的工程价值体现在多机系统里。以三机九节点系统为例通常的做法是先把每台发电机的 Δδ 和 Δω 都作为独立状态引入计算完成后把每个模态按机组聚合参与因子。下表是一个典型结果的示意数值为演示数据排序状态变量模态 A0.42 Hz参与因子1G2 转速 Δω0.4522G2 功角 Δδ0.2213G1 转速 Δω0.1834G1 功角 Δδ0.094解读规则如果某个模态的最大参与因子落在 G2 的转速和功角上且 G1 的参与因子次之这个模态就是 G1 与 G2 之间的相对摆动模式G2 是主导机组。PSS 装在这台机组上效果最直接装在其他机组上可能绕一大圈还压不住。如果参与因子分散在五六台机组上且数值相近说明这是一个区域级模式单靠一台机的 PSS 很难解决需要考虑区域间的协调控制策略。3.4 留数从参与因子到控制点选择参与因子描述的是“状态与模态的关系”而控制设计还需要知道“哪个输入输出通道对这个模态影响最大”这要用到留数Residue。对线性系统从输入 u_j 到输出 y_i 的传递函数可以按模态展开G(s) Σ R_k / (s - λ_k)留数 R_k C_i · V_k · W_k · B_j其中 C、B 是输出、输入矩阵。工程上最常用的做法是对候选机组逐一计算目标模态的留数留数幅值越大说明该通道对模态的能控能观性越强留数相角则决定 PSS 需要提供的相位补偿量。这个指标与参与因子配合使用一个用来选状态一个用来选通道是特征值分析在控制设计里的标准组合拳。4. 特征值指导 PSS 参数整定根轨迹与阻尼比校核4.1 励磁系统为什么可能带来负阻尼快速励磁系统是电压控制的好帮手但它有一个副作用高增益励磁会削弱机电模式的阻尼严重时直接把特征值推到右半平面。物理上的解释是励磁电压调节引入的附加转矩中有一部分与转速变化方向相反构成负阻尼分量。特征值分析能把这个过程量化成一条一目了然的轨迹——励磁增益 KA 从 50 加到 400机电模式的特征值实部如何从 -0.4 逐渐靠近 0 再变正这在复平面上是一段随时可能失控的位移。电力系统稳定器PSS的使命就是提供一个附加励磁信号抵消励磁系统带来的负阻尼。PSS 输出的信号叠加到励磁电压参考点上相位经过校正使附加转矩与转速变化同相。特征值分析法在这个过程中承担的角色是确定目标模式、计算留数相角、扫描增益找到最优补偿、校核整定结果。4.2 PSS 的标准结构与参数窗口工程中绝大多数 PSS 是超前滞后结构典型传递函数为G_pss(s) Kpss · [s·Tw / (1 s·Tw)] · [(1 s·T1)/(1 s·T2)] · [(1 s·T3)/(1 s·T4)]各环节作用和典型参数范围如下环节典型参数作用与整定要点增益 Kpss2 ~ 30太小补偿不足太大会激发高频模态隔直环节 Tw5 ~ 20 s滤除稳态直流分量只对动态分量起作用超前滞后 T1/T20.02 ~ 0.2 s补偿励磁系统与发电机转子回路的相位滞后输出限幅±0.05 ~ 0.1 p.u.防止 PSS 在扰动时过度调节整定的目标非常明确使目标机电模式的阻尼比大于等于 0.05同时不恶化其他模态。实际操作时 T1~T4 由留数相位决定Kpss 由根轨迹扫描决定Tw 一般取 10 秒附近是一个对结果不太敏感的常规值。4.3 用根轨迹扫描增益的最小代码下面这段 MATLAB 脚本演示了最核心的整定动作固定相位补偿参数扫描 Kpss观察目标模式阻尼比的变化。为演示逻辑这里把 PSS 简化为纯增益反馈实际工程中必须带超前滞后环节。%% PSS 增益扫描跟踪机电模式阻尼比 Kpss_list 0:1:20; zeta_list zeros(size(Kpss_list)); for idx 1:length(Kpss_list) Kpss Kpss_list(idx); % 基础 A 矩阵同样来自 SMIB 四阶模型 A [... % 与第 2 节相同的矩阵赋值 0 ws 0 0; -K1/M -D/M -K2/M 0; -K4/Tdo 0 -1/(K3*Tdo) 1/Tdo; -KA*K5/TA -KA*Kpss/TA -KA*K6/TA -1/TA]; lambda eig(A); % 筛选机电模式虚部对应频率 0.1~2 Hz 的共轭对 f imag(lambda) / (2*pi); mech_idx find(abs(f) 0.1 abs(f) 2.0); [~, j] max(abs(imag(lambda(mech_idx)))); lam_mech lambda(mech_idx(j)); sigma real(lam_mech); omega imag(lam_mech); zeta_list(idx) -sigma / sqrt(sigma^2 omega^2); end % 找到阻尼比首次超过 0.05 的最小增益 k_need Kpss_list(find(zeta_list 0.05, 1)); fprintf(阻尼比 0.05 所需的最小 Kpss: %d\n, k_need); % 绘制根轨迹观察增益增大时特征值方向 figure; for idx 1:length(Kpss_list) Kpss Kpss_list(idx); A(4,2) -KA*Kpss/TA; % 只更新反馈项 lambda eig(A); plot(real(lambda), imag(lambda), x); hold on; end grid on; xlabel(实部); ylabel(虚部);代码逻辑的关键在A(4,2) -KA*Kpss/TA这一行第 4 行第 2 列对应励磁电压方程里来自 Δω 的反馈项PSS 的增益通过这个位置进入系统。每次扫描只需要更新这一个元素不需要重新组装整个矩阵。筛选机电模式时限定虚部对应的频率范围确保跟踪的是目标模式而不是励磁模式。参数说明Kpss 从 0 扫到 20步长 1已经覆盖多数同步机组的实用范围。扫描结果通常呈现这样的趋势Kpss 从 0 到某一点阻尼比先升后降超过某个临界值后原本的高频励磁模式开始向虚轴移动甚至出现新的不稳定。这就是“PSS 增益不是越大越好”的数学解释。4.4 整定后的四步校验增益定下来之后不要急着交付按下面的顺序复核一遍计算完整状态矩阵的全部特征值确认 0.1~2 Hz 范围内所有模式的阻尼比都大于 0.03而不是只盯着目标模式。在时域仿真中对参考电压施加 1% 阶跃扰动观察功率和转速的响应曲线确认振荡在几次摆动内明显衰减。检查 PSS 输出是否频繁触及限幅值。如果限幅动作过于频繁说明 Kpss 偏大或者输入信号中包含过多噪声需要回退增益。换一个运行点重新计算特征值。重负荷方式下阻尼通常变差如果新运行点下阻尼比接近 0.02说明整定余量不足。这四步里最容易被忽略的是第 4 步。实际电网考核的是“最严重方式下仍然稳定”而不是“典型方式下稳定”。用特征值分析做 PSS 整定一定要把整个运行方式集合作为参数空间的一部分。5. 工具链与最小复现流程从潮流结果到特征值5.1 常用工具与选型特征值分析在不同场景下有不同的工具选择。做研究原型验证和教学演示时我一般用 MATLAB 或 Python 自己组装矩阵灵活且可控做工程交付时商用软件是主流它们的模态分析模块把这些流程固化成了标准操作。工具类型优势边界MATLAB PST学术免费代码开放便于教学和二次开发需要自己搭数据文件PSAT学术免费潮流与小信号分析一体文档和社区规模有限Python scipy开源适合全流程自动化与批处理多机系统建模工作量较大PSS/E、DIgSILENT商业模态分析成熟工程认可度高授权成本高模型可解释性弱这里不展开具体商业软件的操作因为它与版本强相关。下面用一个 Python 最小示例说明从模型到特征值的完整路线这段代码在普通笔记本上可以直接运行。5.2 一个可以直接跑的 Python 双机摆方程示例双机系统是理解多机振荡的最小载体两台发电机通过等效电抗相连每台机用经典二阶模型描述状态变量取 Δδ₁、Δω₁、Δδ₂、Δω₂。其线性化状态矩阵可以直接写出来import numpy as np ws 2 * np.pi * 50 # 50 Hz 系统的同步角速度 # 两台机组的惯性时间常数和阻尼系数 M np.array([6.0, 8.0]) # 2H单位秒 D np.array([0.02, 0.03]) # 阻尼系数 Ks 0.8 # 等效同步转矩系数标幺值 # 状态顺序: [Δδ1, Δω1, Δδ2, Δω2] A np.array([ [0, ws, 0, 0], [-Ks/M[0], -D[0]/M[0], Ks/M[0], 0], [0, 0, 0, ws], [Ks/M[1], 0, -Ks/M[1], -D[1]/M[1]] ]) lam np.linalg.eigvals(A) print(特征值实部 虚部 频率(Hz) 阻尼比) for l in lam: sigma, omega l.real, l.imag f omega / (2 * np.pi) zeta -sigma / np.sqrt(sigma**2 omega**2) print(f{sigma:10.4f} {omega:10.4f} {f:10.4f} {zeta:10.4f})逻辑说明第一行到第三行定义基础参数惯性时间常数 M 的单位是秒阻尼系数是无量纲的标幺值Ks 代表两台机之间的同步功率系数。状态矩阵的结构与单机系统一脉相承每个功角方程都乘以 ws 连接到转速每个转速方程中自己机组的功角项是负反馈另一台机组的功角项是正耦合这正是相对振荡的数学来源。运行后会看到一对共轭特征值其虚部对应双机系统的固有振荡频率实部对应阻尼水平。参数说明Ks 取 0.8 代表中等强度的机间电气耦合实际系统中联络线越弱 Ks 越小振荡频率也随之降低。把这段代码里的 M、D、Ks 换成具体工程数据就是最简化的区域间振荡模型。5.3 从潮流到特征值的完整复现步骤上面的双机例子跳过了潮流计算因为参数是直接给定的。真实工程流程中特征值分析的上游永远是潮流结果。完整链路如下用潮流计算工具求取系统运行点得到各节点电压幅值、相角和发电机注入功率。从潮流解中提取换流变压器变比、发电机端电压等参数代入同步电机模型计算 Phillips-Heffron 常数或等效的状态矩阵元素。把所有动态元件的线性化方程按状态变量顺序拼装成全局 A 矩阵节点网络方程以代数约束形式消去。调用特征值求解器得到全部特征值按频率段归类为区间模式、本地模式、控制模式。对每个关注模态计算参与因子与留数输出机组排序和控制点推荐。这个流程里最常出问题的环节是第 2 步和第 3 步之间的单位制转换。发电机参数一般以自身容量为基准而网络方程以系统基准容量为基准两者不经换算直接拼装会导致特征值数量级完全失真。不同文献对 M 的定义也有差异有的用 2H有的直接写 H换算因子漏乘一个 2 就会让振荡频率偏差 40%这类错误在特征值分析里极难靠肉眼发现。6. 用时域仿真校准特征值Prony 分析与小技巧特征值分析算出来的阻尼比本质上是线性化模型的推论模型参数不准或者线性化点偏离实际结果就会有偏差。因此工程上有一个交叉验证的常规动作在完整时域仿真模型里施加一个小扰动记录功率或转速响应然后从这段响应曲线里反向提取振荡频率和阻尼比与特征值计算结果对照。误差在 5% 以内说明小信号模型是可信的对不上就说明建模时丢掉了某个动态环节。提取振荡参数最常用的方法是 Prony 分析它能把时域信号分解为若干指数衰减正弦分量的叠加。下面给出用 Python 实现的简化思路import numpy as np def prony_analysis(signal, dt, n_terms4): 从时域信号提取主导振荡频率和阻尼比的简化实现 signal: 一维数组扰动后的功率或转速响应 dt: 采样间隔单位秒 N len(signal) # 构造线性预测方程阶数为 2 * n_terms L 2 * n_terms if N 2 * L: raise ValueError(信号长度不足以支持 Prony 分析) # 最小二乘求解自回归系数 A np.array([signal[i:iL] for i in range(N - L)]) b -signal[L:] coeff np.linalg.lstsq(A, b, rcondNone)[0] # 特征多项式求根 roots np.roots(np.concatenate([[1], coeff[::-1]])) dt_vals np.log(np.abs(roots)) / dt f_vals np.arctan2(roots.imag, roots.real) / (2 * np.pi * dt) % 50 return dt_vals, f_vals逻辑说明Prony 分析的核心思路是先拟合信号的自回归模型再从特征多项式根的位置换算成频率和衰减系数。上述实现是教学级简化版工程信号中的噪声和高阶模态会让结果发散建议用 MATLAB 的 prony 函数或专门的系统辨识工具箱。使用时的检验标准是提取出的主导模式数量要大于等于特征值分析中关注的目标模态数这样才不会漏振荡。还有一个容易被忽略的排错技巧在特征值分析和小信号建模阶段把时域仿真里施加的扰动幅度设置在额定值的 1% 以内可以保证响应基本落在线性化假设的范围内如果扰动加得太大仿真结果里混入了非线性畸变成分Prony 分析提取出的阻尼比会偏离特征值计算结果这种偏差是方法本身带来的不是模型错误。反过来如果扰动很小但两边结果仍然差异明显就得回头检查动态参数的基准值换算。最后一个推荐做法把特征值分析结果和 Prony 辨识结果放到同一张波特图或复平面上对比既能看到频率是否吻合也能直观看到阻尼比的差距。这套交叉验证流程一旦跑通后续每次更换运行点或调整 PSS 参数都可以快速回归一遍比单纯依赖特征值计算要可靠得多。本文还有配套的精品资源点击获取