简介面向电力系统研究人员与调度自动化工程师的暂态稳定分析资源包聚焦李雅普诺夫直接法在多机系统在线稳定评估中的应用。内容在传统能量函数基础上引入转移电导与阻尼修正降低保守性并提升计算效率以辽南系统算例验证改进效果同时给出完整Python实现涵盖转子运动方程建模、能量函数计算、不稳定平衡点搜索、故障场景仿真与稳定性判据。包体为1份docx文档压缩包大小约50KB文档内含理论推导、公式解析、可运行代码及逐步解释便于按图索骥复现实验。现有103人学习适合具备电力系统基础和编程能力的读者用于算法验证、在线评估工具开发或调度自动化实时性改进。通过该资源可掌握从机理建模到代码落地的完整链路并理解参数调优与计算效率优化思路。1. 为什么在线暂态稳定分析要绕开数值积分电力系统暂态稳定分析的经典做法是逐步积分法即对描述发电机转子运动的微分方程组做时域数值积分故障前、故障中、故障后三段各自建模步长取得越密结果越准。这个方法在离线分析中没有任何问题但放到在线动态安全评估场景里就会遇到一个尴尬的矛盾系统规模越大需要积分的微分方程数量越多计算时间随节点数非线性增长而调度自动化给稳定评估留出的时间窗口往往只有秒级。李雅普诺夫直接法的思路完全不同它不跟踪故障后的完整运动轨迹而是构造一个能量型函数比较故障切除时刻的系统能量与系统能容纳的临界能量大小关系直接决定系统是否稳定。这个判断过程只需要一次函数值计算计算量几乎不随故障持续时间的长度变化。下面把这套方法的数学构造、代码实现和故障场景仿真拆开讲清楚使用的代码直接复现了论文中考虑转移电导的能量型 Lyapunov 函数并包含了不稳定平衡点的搜索逻辑。2. 能量型 Lyapunov 函数的构造从转子运动方程到可计算代码2.1 多机系统经典模型与转子运动方程的离散化多机系统的暂态稳定分析通常采用经典二阶模型每台发电机用一个转子运动方程描述。对第 i 台发电机转子运动方程写作M_i * d(omega_i)/dt P_mi - P_ei - D_i * omega_i d(delta_i)/dt omega_i其中 M_i 是惯性时间常数D_i 是阻尼系数P_mi 是机械功率P_ei 是电磁功率delta_i 是转子功角omega_i 是相对同步转速的角速度偏差。电磁功率 P_ei 不是独立的它由全部发电机的功角差和网络导纳矩阵共同决定。这里的关键是把网络方程压缩到发电机内节点上得到一个 n 阶的导纳矩阵 Y然后通过下式计算每台发电机的电磁功率P_ei sum_j V_i * V_j * (G_ij * cos(delta_i - delta_j) B_ij * sin(delta_i - delta_j))G 是电导矩阵实部B 是电纳矩阵虚部。经典模型中电压幅值 V 视为恒定所以 P_ei 只依赖于功角差。下面这段代码实现了这个计算过程def swing_equation(self, x, t): delta x[:self.n] omega x[self.n:] Pe np.zeros(self.n) for i in range(self.n): for j in range(self.n): Pe[i] self.V[i] * self.V[j] * ( self.G[i,j] * np.cos(delta[i] - delta[j]) self.B[i,j] * np.sin(delta[i] - delta[j]) ) ddelta omega domega (self.P_m - Pe - self.D * omega) / self.M return np.concatenate((ddelta, domega))这个方法是转子运动方程的右侧函数直接传给 SciPy 的odeint做数值积分。注意它的输出格式是状态导数的拼接向量前 n 个元素是功角导数后 n 个元素是角速度导数这个顺序要求调用方严格对应。P_m 与 Pe 的差值代表加速功率阻尼项 D*omega 的作用是让系统在稳态扰动后逐渐收敛如果阻尼设置为 0系统会保持等幅振荡不衰减这和实际电网中励磁系统及原动机提供的阻尼效果不符合所以 D 不能省略。2.2 Lyapunov 函数的三段式表达式及其物理含义能量型 Lyapunov 函数在直接法中扮演的角色是系统的能量计。它的构造逻辑是把系统总能量拆成动能和势能两部分动能项对应转子相对运动的动能势能项对应功角偏离稳定平衡点所储存的位能。论文中给出的函数形式如下def lyapunov_function(self, delta, omega): V_kin 0.5 * np.sum(self.M * omega**2) V_pot 0.0 for i in range(self.n): for j in range(i1, self.n): V_pot self.V[i] * self.V[j] * ( -self.G[i,j] * (np.cos(delta[i] - delta[j]) - np.cos(self.delta0[i] - self.delta0[j])) self.B[i,j] * (np.sin(delta[i] - delta[j]) - np.sin(self.delta0[i] - self.delta0[j])) ) V_trans 0.0 for i in range(self.n): V_trans self.P_m[i] * (delta[i] - self.delta0[i]) return V_kin V_pot V_trans动能项的计算是逐台发电机累加的M 是惯性常数omega 是角速度偏差这一项在故障切除瞬间最大随着系统振荡逐渐转化为势能。势能项的累加只遍历 ij 的上三角部分这是因为互导纳项在交换下标后对称重复累加会引入 2 倍系数错误。注意势能项里同时出现了 G 和 BB 项来自正常电磁功率路径G 项来自转移电导的贡献。经典李雅普诺夫函数在处理转移电导时做了忽略处理这里保留 G 项是为了更精确地反映能量在电阻性支路上的耗散路径。V_trans 是机械功率的修正项它保证函数在稳定平衡点处取值为 0否则能量函数会有一个常数偏移导致临界能量比较失真。2.3 故障切除时刻能量与临界能量的比较逻辑直接法的稳定性判据可以写成一个非常直接的不等式如果故障切除时刻的系统能量 V_fault 小于临界能量 V_critical系统是暂态稳定的反之则不稳定。这个判据不需要继续积分故障后的轨迹这正是直接法计算速度快的根源。实际应用中的关键问题是 V_critical 如何确定它在理论上等于系统稳定边界上某个不稳定平衡点对应的 Lyapunov 函数值。找到一个合理的不稳定平衡点是直接法应用中最容易出错的环节下一章专门处理这个问题。3. 不稳定平衡点搜索与临界能量直接法判断稳定性的关键3.1 不稳定平衡点的理论数量与工程近似理论上n 机系统有 2^(n-1)-1 个不稳定平衡点每个点对应一种失步模式。以三机系统为例理论上有 3 个不稳定平衡点分别对应 1 号机失步、2 号机失步和 1、2 号机同时相对 3 号机失步这三种模式。随着发电机数量增加不稳定平衡点数量指数增长逐一搜索在在线计算中根本不现实。工程上的常见做法是只搜索主导失步模式对应的不稳定平衡点即某一台发电机相对其余机组失步的情况忽略多机同时失步的组合模式。这种近似会引入一定的保守性但计算量大幅下降论文中辽南系统的计算也验证了这种做法在工程上是可接受的。def find_unstable_equilibrium(self, methodapproximate): if method approximate: n self.n ref_gen n - 1 unstable_points [] for i in range(n - 1): delta_ue np.zeros(n) delta_ue[i] np.pi unstable_points.append(delta_ue) V_values [] for delta in unstable_points: V self.lyapunov_function(delta, np.zeros(self.n)) V_values.append(V) min_idx np.argmin(V_values) return unstable_points[min_idx], V_values[min_idx]这段近似搜索的逻辑是逐个假设第 i 台发电机相对参考机失步把它的功角置为 π其余发电机的功角保持初始值然后计算对应的 Lyapunov 函数值。取所有失步模式中能量最小的那个作为临界能量的候选这是偏保守的处理方式。功角置 π 的物理含义是失步发电机与参考机完全反相此时系统处于最接近稳定边界的状态。如果换一个更大的角度如 2π会得到完全不同的能量值那个值对应的是多转一圈后的不稳定平衡点不是首摆稳定边界不能用于暂态稳定判断。3.2 牛顿法精化从近似解到精确不稳定平衡点近似方法得到的不稳定平衡点精度有限因为真实的不稳定平衡点需要满足所有发电机的电磁功率与机械功率相等即 RHS 全部为 0。功角直接置 π 的做法通常不能满足这个条件。更严谨的做法是以近似解为初值用牛顿法在功率平衡方程上迭代求解。下面是用 SciPy 的newton函数做精化的实现def find_unstable_equilibrium(self, methodnewton): from scipy.optimize import newton def equilibrium_eq(delta): Pe np.zeros(self.n) for i in range(self.n): for j in range(self.n): Pe[i] self.V[i] * self.V[j] * ( self.G[i,j] * np.cos(delta[i] - delta[j]) self.B[i,j] * np.sin(delta[i] - delta[j]) ) return self.P_m - Pe delta_approx, _ self.find_unstable_equilibrium(methodapproximate) delta_ue newton(equilibrium_eq, delta_approx) V_ue self.lyapunov_function(delta_ue, np.zeros(self.n)) return delta_ue, V_ueequilibrium_eq返回的是机械功率与电磁功率的差值向量牛顿法迭代的目标是让这个差值逼近 0。这里有一个参数选择的细节newton默认使用一阶差分估计雅可比矩阵对 n 机系统需要 n1 次函数求值精度比近似方法高一个量级。但牛顿法对初值敏感如果初始猜测距离真实解太远可能收敛到平凡的稳定平衡点而非不稳定平衡点。验证的方法是检查解对应的势能是不是局部的鞍点值如果不稳定平衡点搜索正确对应的 V_ue 应该明显大于稳定平衡点处的 0同时小于相邻其他失步模式的能量值。如果牛顿法收敛到 0 附近的解说明初值有问题需要换一个失步模式作为初始猜测。3.3 临界能量计算的稳定性讨论临界能量 V_critical 的取值直接决定了稳定性判据的保守程度。V_critical 取的是多个失步模式中能量的最小值这意味着系统能量还没达到某些模式的临界值时就已经被判为不稳定结果偏保守。反过来如果只取最大值的失步模式判据会偏乐观可能给出错误的稳定结论。工程上折中的做法是取最小值作为 V_critical宁保守勿冒进因为暂态稳定判断失误的代价远高于保守带来的误切。另外需要注意阻尼系数对临界能量的影响阻尼项的物理作用是耗散能量理论上会让系统的真实稳定域比无阻尼时略大但直接法构造的 Lyapunov 函数在数学上要求沿系统轨迹单调不增有阻尼项时这一条件仍然成立V_critical 的数值变化不大。4. 故障场景仿真从故障切除时刻到稳定性判定4.1 故障仿真的流程设计故障场景仿真遵循一个固定的三段式流程故障前系统运行在稳定平衡点故障发生时网络导纳矩阵改变电磁功率传输路径被削弱故障切除后系统恢复到原始网络拓扑。代码实现的简化处理方式是不区分故障中和故障后的网络参数而是在故障时间段内用修改后的导纳矩阵积分切除时刻之后切回原矩阵。更精细的做法是给出三个不同的导纳矩阵 Y_pre、Y_fault、Y_post在每个时间区间用对应的参数积分。下面给出包含故障仿真的稳定性检查方法def stability_check(self, fault_duration, t_end10, dt0.01): delta_ue, V_critical self.find_unstable_equilibrium() t, sol, V self.simulate(t_end, dt) idx np.argmin(np.abs(t - fault_duration)) V_fault V[idx] stable V_fault V_critical plt.figure(figsize(10, 6)) plt.plot(t, V, labelLyapunov Function) plt.axhline(V_critical, colorr, linestyle--, labelCritical Value) plt.axvline(fault_duration, colorg, linestyle:, labelFault Cleared) plt.xlabel(Time (s)) plt.ylabel(Lyapunov Function Value) plt.title(fSystem is {stable if stable else unstable} for fault duration {fault_duration}s) plt.legend() plt.grid(True) plt.show() return stable故障切除时刻的能量 V_fault 取自 Lyapunov 函数随时间变化的曲线。idx是时间轴上最接近故障持续时间的下标dt 取 0.01 秒时时间量化误差不超过 5 毫秒对稳定性判断的影响很小。但如果故障持续时间很短比如 0.05 秒量化误差占比就到了 10%这时需要把 dt 缩小到 0.001 量级重新仿真。图上那条绿色虚线标出的位置很关键故障切除时刻的系统能量如果低于红色虚线代表的临界值判为稳定否则判为不稳定。4.2 稳定性判据的物理解读与陷阱比较 V_fault 与 V_critical 时有一个隐含假设故障切除后系统能量沿 Lyapunov 函数单调不增。这个假设在理想的能量型函数下成立但实际构造的函数由于转移电导项的存在并不严格满足单调递减条件。也就是说V 在故障切除后的某个时刻可能短暂升高然后在阻尼作用下再回落。直接比较切除时刻的能量会忽略这种先升后降的过渡行为造成偏保守的判断。处理方式有两种一种是把比较点向后推移到 V 达到峰值的时刻另一种是引入能量裕度系数即要求 V_fault 小于 V_critical 乘以一个小于 1 的系数比如 0.9才判为稳定。后者在工程中更常用计算量没有增加但稳定性结论更可靠。4.3 仿真结果的可视化解读仿真的输出是三个子图功角轨迹、角速度轨迹、Lyapunov 函数轨迹。功角曲线能直观看出各发电机是否同步摇摆如果两条功角曲线持续分开且没有收敛趋势对应系统失去同步。角速度轨迹的衰减速度反映阻尼的效果阻尼系数越大角速度振荡衰减越明显。Lyapunov 函数曲线在故障期间迅速上升故障切除后如果系统稳定V 应该缓慢下降或保持平稳如果 V 持续攀升说明系统正在失去稳定。判读 Lyapunov 函数曲线时要注意它在振荡过程中可能出现的上下起伏不应因为单次下降就认定稳定要看整体的趋势是否收敛到一个小于 V_critical 的水平。5. 转移电导的三种改进方法减少保守性的工程取舍5.1 转移电导对能量函数构造的影响转移电导就是导纳矩阵中的电导实部 G_ij它对应输电线路和变压器的电阻损耗。经典构造 Lyapunov 函数时会忽略 G_ij理由是能量函数在存在转移电导时难以满足严格的单调性条件。但实际电网中转移电导不小忽略它会让临界能量计算结果与真实稳定域偏差明显。论文中给出了三种考虑转移电导的改进方法第一种方法直接把电纳 B_ij 替换为导纳幅值 sqrt(G_ij^2 B_ij^2)第二种方法将势能项的系数调整为 G_ij 与导纳幅值的比值第三种方法在势能项中引入与功角差线性相关的修正项。下面这段代码实现了三种改进的势能项计算def improved_lyapunov_function(self, delta, omega, method1): V_kin 0.5 * np.sum(self.M * omega**2) if method 1: V_pot 0.0 for i in range(self.n): for j in range(i1, self.n): Y_mag np.sqrt(self.G[i,j]**2 self.B[i,j]**2) V_pot self.V[i] * self.V[j] * ( -self.G[i,j] * (np.cos(delta[i] - delta[j]) - np.cos(self.delta0[i] - self.delta0[j])) Y_mag * (np.sin(delta[i] - delta[j]) - np.sin(self.delta0[i] - self.delta0[j])) ) elif method 2: V_pot 0.0 for i in range(self.n): for j in range(i1, self.n): Y_mag np.sqrt(self.G[i,j]**2 self.B[i,j]**2) V_pot self.V[i] * self.V[j] * ( -self.G[i,j]/Y_mag * (np.cos(delta[i] - delta[j]) - np.cos(self.delta0[i] - self.delta0[j])) (np.sin(delta[i] - delta[j]) - np.sin(self.delta0[i] - self.delta0[j])) ) elif method 3: V_pot 0.0 for i in range(self.n): for j in range(i1, self.n): V_pot self.V[i] * self.V[j] * ( -self.G[i,j] * (np.cos(delta[i] - delta[j]) - np.cos(self.delta0[i] - self.delta0[j]) self.B[i,j] * (np.sin(delta[i] - delta[j]) - np.sin(self.delta0[i] - self.delta0[j])) self.G[i,j] * (delta[i] - delta[j] - (self.delta0[i] - self.delta0[j])) * np.sin(self.delta0[i] - self.delta0[j]) )) V_trans 0.0 for i in range(self.n): V_trans self.P_m[i] * (delta[i] - self.delta0[i]) return V_kin V_pot V_trans三种方法的核心差异在于如何处理 G_ij 项与前一项的耦合关系。方法 1 用导纳幅值替换电纳相当于把转移电导的损耗折算到无功功率路径上实现最直接但物理意义的对应关系相对模糊。方法 2 的构造方式是让势能项的幅值系数归一化G_ij 与 Y_mag 的比值在 G_ij 越大时越接近 1势能项受转移电导的调制更平滑。方法 3 最精细在线性化项里保留了 delta_i - delta_j 与初始功角差的差这个差在小扰动下趋于 0避免了大角度偏差时修正项过度放大势能值。三种方法对临界能量的影响程度可以用下面的对比关系来说明改进方法势能修正逻辑临界能量变化趋势适用场景方法 1B_ij 替换为导纳幅值临界能量提高保守性下降线路电阻占比高的网络方法 2系数归一化处理临界能量变化平缓多种运行方式切换频繁方法 3线性化修正项临界能量最接近真实稳定域高精度离线校核表格中的临界能量变化趋势是从数值实验和理论推导两方面得出的。方法 1 的修正幅度较大在重负荷场景下可能把临界能量抬高 5% 到 10%方法 2 的修正相对温和适合作为默认选项方法 3 的精度最高但计算量最大在线应用时需要评估时间预算。实际使用中我一般先用方法 1 快速扫描失步风险对风险较高的场景再用方法 3 精确计算。5.2 改进方法在故障场景中的验证验证三种改进方法的实际效果需要在同一故障场景下分别计算临界能量 V_critical 和故障切除时刻能量 V_fault然后比较稳定性判断结果。如果某种方法判稳定而参考时域仿真结果判不稳定说明该方法偏乐观临界能量取值过高反之判不稳定而时域仿真判稳定说明该方法偏保守。在示例三机系统中方法 1 计算出的临界能量最高方法 3 的最低方法 2 居中。这个排序在物理上可以这样理解方法 1 把电纳替换为更大的导纳幅值等效地增大了发电机间的电气距离临界能量自然偏高方法 3 的线性化修正抵消了一部分 G_ij 带来的势能增长临界能量更接近真实值。三者的比较还受初始功角差的影响初始功角差越大方法间的差异越明显。6. 从仿真到在线应用参数调优与四个落地细节直接法在线应用的瓶颈不在算法本身而在参数维护和计算流程的工程化。以下四个细节是实际落地时最容易踩坑的地方。第一个细节是阻尼系数 D 的整定。论文中的阻尼系数设定为固定值但实际系统中阻尼与运行方式强相关励磁系统的 PSS 投入情况会影响等效阻尼。常见做法是先按离线方式刷一遍小扰动数据用 Prony 分析提取振荡衰减时间常数反推等效阻尼系数。在线运行时如果阻尼估计偏大Lyapunov 函数中的动能项衰减速度被高估稳定性结论可能偏乐观所以阻尼系数建议取下限值。第二个细节是导纳矩阵的实时更新。电网拓扑变化频繁每次开断线路都要重新计算导纳矩阵否则 Lyapunov 函数计算使用的 G 和 B 与实际网络不匹配。工程实现时建议保留一份节点支路模型拓扑变化时只重新计算受影响区域的导纳子矩阵不必全矩阵重建。三机示例中的矩阵维度低刷新时间不是问题但实际大网中导纳矩阵刷新耗时占了在线计算的相当比例值得优化。第三个细节是故障切除时刻能量的采样方式。仿真步长 dt 的取值直接影响 V_fault 的精度但盲目缩小步长会增加计算量。折中方案分两段处理故障期间用小步长精细积分故障切除后放大步长。对于切除时刻的 V 值用步长 0.005 秒能比 0.01 秒减少约一倍的量化误差时间开销增加也可控。第四个细节是临界能量的匹配更新。不稳定平衡点取决于网络拓扑和运行方式不随故障位置变化所以可以离线预先计算多个典型运行方式下的 V_critical在线运行时按当前运行方式查表取值。运行方式剧烈变化后需要重新计算时用上一章的牛顿法以历史解为初值做少量迭代即可收敛计算量远小于从零搜索。临界能量的值会受故障切除角度分布影响换算出标幺值时要同步检查基准功率是否一致。本文还有配套的精品资源点击获取