简介电力系统暂态能量函数法暂态稳定分析PPT学习教案是一份围绕暂态能量函数法的系统教学课件面向电力系统专业学生、考研备考生及从事稳定分析的研究人员旨在帮助读者从能量视角掌握暂态稳定判据与稳定裕度评估方法。压缩包仅含1个PPTX格式文件大小约1.23MB共105页左右结构完整便于按章学习。目前已吸引160人浏览学习。课件从古典力学能量概念出发系统讲解暂态能量函数的基本原理、数学描述与李氏定理结合单机无穷大系统直接法暂态稳定分析逐步过渡到多机系统的特殊问题、两种坐标系表达并重点介绍最近不稳定平衡点法AUEP、等位能线与位能边界、相关不稳定平衡点法RUEP等内容。其中还包含球碗运动模型、功角特性曲线等图文示例可辅助理解临界能量、稳定域和故障切除时间等关键概念。学习后可构建基于李雅普诺夫函数的暂态稳定分析框架适用于课程教学、自学复习及毕业设计参考。1. 暂态能量函数法把暂态稳定分析从“看曲线”变成了“算余量”时域仿真跑完一条功角曲线报告里能写的基本就是“最大摆角 126°未失稳”。可电网调度侧真正要的不是这一个布尔值稳定裕度还剩多少、哪些故障更危险、切机切负荷该留多大反应时间都依赖连续量化结果。暂态能量函数法就是在故障切除时刻停下积分把发电机转子动能、相对位置势能和网络磁场储能折算成一张系统的能量收支表拿它和系统能承受的临界能量直接比较。产出不是曲线而是稳定裕度。裕度从 0.8 跌到 -0.2 时你已经知道系统正在翻过界面而不是等功角甩开后才去事后确认。做在线暂态稳定评估、故障筛选和预防控制的工程师值得把整套方法吃透。下面先用单机无穷大系统把能量函数拆开再给出可复现的 Python 计算最后落回多机系统里 COI 坐标和 PEBS/BCU 的实际工程问题。2. 从等面积法则到李雅普诺夫能量判据暂态能量函数法的三个组成项单机无穷大系统里能量函数有一个特别好懂的退化形态就是本科教科书里的等面积法则。看懂它后面 Vcr、PEBS 这些概念都能对号入座。2.1 等面积法则为什么是单机的暂态能量函数雏形单机无穷大系统的功角摇摆方程为M * d²δ/dt² Pm - Pe其中 M 是惯性常数Pm 是原动机机械功率Pe EV/XΣ·sinδ 是发电机输出的电磁功率。发生短路后 Pe 骤降转子获得净加速功角从 δ0 开始增大保护动作切除故障后网络电抗变大Pe 恢复转子进入减速段。把加速段和减速段写成面积形式加速面积A_acc ∫(Pm - Pe)dδ积分区间从 δ0 到切除时刻 δc减速面积A_dec ∫(Pe - Pm)dδ积分区间从 δc 到最大摇摆角 δmax稳定条件是 A_acc ≤ A_dec。把加速面积理解为转子相对同步转速储存的动能把减速面积理解为制动过程中动能释放的功等面积法则其实就是能量守恒在一个自由度上的投影。单机两群摆动时它完全够用一旦系统变成多机每台机的转速都与其余机不同“整台机组相对无穷大系统”的概念消失等面积法则便不再成立。但能量法的思路保留了下来找一个标量函数让它沿故障后轨迹单调不增再拿初值处能量和稳定域边界能量比较就得到稳定裕度。2.2 暂态能量函数的构成动能 Vk、位置势能 Vp、磁场储能 Vmag单机经典模型下能量函数可以写成显式形式V(δ, ω) 0.5 * M * (ω - ωs)² - Pm * (δ - δs) - Pmax * (cosδ - cosδs)其中 ωs 是同步电角速度δs 是故障后系统的稳定平衡点SEP。这个表达式由三项构成动能项Vk 0.5 * M * (ω - ωs)²转子相对同步转速运动的动能恒为正故障期间转子加速越快这一项越大。位置势能项Vp_pos -Pm * (δ - δs)机械功率沿功角方向累积的能量相当于重力场里“高度”被功角替代。磁场储能项Vp_mag -Pmax * (cosδ - cosδs)来自发电机内电势与网络之间的耦合储能反映电磁功率与功角之间的保守场关系。后两项合称势能 Vp。动能项恒为正势能项在功角摆大时通常上升。把参考点定在哪个平衡点上是个容易踩的坑必须是故障切除后网络拓扑下的稳定平衡点 δs不是故障前的稳态功角 δ0。故障前后网络结构往往不同XΣ 变了δs 也变了直接拿故障前 δ0 当参考能量初值会带系统性偏差。2.3 稳定判据 V(t) Vcr 和李雅普诺夫定理的工程化用法李雅普诺夫定理的工程化读法是如果能找到正定函数 V使其沿故障后轨迹满足 V̇ ≤ 0那么只要切除时刻的能量 Vcl 小于临界能量 Vcr系统就是暂态稳定的反过来 Vcl ≥ Vcr 则失稳。Vcr 在数学上是稳定域边界上的最小能量工程上则按计算量从小到大有三条求法求法计算量精度与特点工程注意点直接法UEP 法中对单机精确多机需要求控制不稳定平衡点收敛性依赖初值求 UEP 的迭代易发散通常用牛顿-拉夫逊法搭配同伦法PEBS 法小近似值势能界面穿越点偏乐观沿持续故障轨迹找势能峰值第一个峰值才有意义BCU 法大精度最好对多机强非线性模型稳健先找出口点再沿梯度系统求 CUEP实现复杂度高这里有个反直觉结论能量判据不替代时域仿真而是做互补。它建立在理想化模型上电网建模越细励磁、PSS、负荷特性、饱和严格意义上的李雅普诺夫函数越难构造。工程上通常用能量法做快速筛选用精确暂态仿真做最终确认两种手段对同一批故障互相校核。3. 用 Python 算单机无穷大系统的暂态稳定裕度CCT 二分搜索回到可复现的最小实现。一台汽轮发电机经双回线并入无穷大母线0 秒时机端三相对地短路tc 时刻保护切除故障线路只剩一回线运行。目标求临界切除时间 CCT并同时输出能量裕度。3.1 模型参数与摇摆方程模型采用二阶经典摇摆方程Pe 与功角的关系随故障阶段切换。三个电抗分别对应三种网络拓扑正常态总电抗 X1 0.30 p.u.双回线并列运行故障态总电抗 X2 1.50 p.u.代表三相短路后残压较低、仍有残余转移功率的场景切除后总电抗 X3 0.40 p.u.只剩一回线机械功率 Pm 0.9 p.u.内电势 E 和无穷大母线电压均取 1.0 p.u.。惯性时间常数 H 3.5 秒换算成摇摆方程里的惯性常数需要除以同步电角速度M 2H / ωsωs 2π·50 ≈ 314.16 rad/s。3.2 完整可运行的 RK4 积分与能量计算代码下面的代码不依赖 scipy只用标准库 math复制即可运行。核心是四阶龙格-库塔积分器外加能量函数计算与 CCT 二分搜索。import math # ---------- 系统参数 ---------- f0 50.0 omega_s 2.0 * math.pi * f0 # 同步电角速度, rad/s H 3.5 # 惯性时间常数, s M 2.0 * H / omega_s # 摇摆方程惯性常数, p.u. Pm 0.9 # 机械功率, p.u. E1 1.0 # 发电机内电势, p.u. Vb 1.0 # 无穷大母线电压, p.u. X1 0.30 # 正常态总电抗 X2 1.50 # 故障态总电抗三相短路残余 X3 0.40 # 切除后总电抗 Pmax_normal E1 * Vb / X1 Pmax_fault E1 * Vb / X2 Pmax_post E1 * Vb / X3 # 故障前与故障后的稳定平衡点功角 delta_s math.asin(min(1.0, Pm / Pmax_normal)) delta_sp math.asin(min(1.0, Pm / Pmax_post)) delta_uep math.pi - delta_sp # 控制UEP单机系统为π-δs # 临界能量 Vcr势能在 UEP 处的取值 Vcr (-Pm * (delta_uep - delta_sp) - Pmax_post * (math.cos(delta_uep) - math.cos(delta_sp))) print(f稳态功角: {math.degrees(delta_s):.2f} deg) print(f故障后SEP: {math.degrees(delta_sp):.2f} deg, UEP: {math.degrees(delta_uep):.2f} deg) print(f临界能量 Vcr {Vcr:.6f} p.u.) # ---------- RK4 积分器 ---------- def derivative(state, pm, pmax): delta, omega state Pe pmax * math.sin(delta) return [omega - omega_s, (pm - Pe) / M] def rk4_step(state, dt, pm, pmax): k1 derivative(state, pm, pmax) s2 [state[0] dt/2*k1[0], state[1] dt/2*k1[1]] k2 derivative(s2, pm, pmax) s3 [state[0] dt/2*k2[0], state[1] dt/2*k2[1]] k3 derivative(s3, pm, pmax) s4 [state[0] dt*k3[0], state[1] dt*k3[1]] k4 derivative(s4, pm, pmax) d_delta (k1[0] 2*k2[0] 2*k3[0] k4[0]) / 6.0 d_omega (k1[1] 2*k2[1] 2*k3[1] k4[1]) / 6.0 return [state[0] dt*d_delta, state[1] dt*d_omega] # ---------- 能量函数 ---------- def energy(delta, omega): vk 0.5 * M * (omega - omega_s) ** 2 vp (-Pm * (delta - delta_sp) - Pmax_post * (math.cos(delta) - math.cos(delta_sp))) return vk vp, vk, vp # ---------- 给定切除时间判断是否稳定 ---------- def is_stable(tc, t_max2.0, dt0.0005): state [delta_s, omega_s] t 0.0 # 故障期间积分到 tc while t tc - 1e-12: h min(dt, tc - t) state rk4_step(state, h, Pm, Pmax_fault) t h vcl, _, _ energy(state[0], state[1]) if vcl Vcr: return False, vcl # 能量已经越过临界直接判不稳 # 切除后继续积分观察功角是否回摆或越限 while t t_max: state rk4_step(state, dt, Pm, Pmax_post) t dt d state[0] # 功角越过UEP且仍在加速 失稳 if d delta_uep and state[1] omega_s: return False, vcl # 功角回摆到SEP附近且速度接近同步 稳定 if d delta_sp 0.05 and state[1] omega_s 0.1: return True, vcl return True, vcl # ---------- 二分搜索 CCT ---------- lo, hi 0.01, 1.0 for _ in range(40): mid (lo hi) / 2.0 ok, _ is_stable(mid) if ok: lo mid else: hi mid cct (lo hi) / 2.0 print(fCCT ≈ {cct:.4f} s) # ---------- 扫描切除时间输出能量表 ---------- print(\n tc(s) Vcl(p.u.) Vcr(p.u.) 裕度(%) 时域判稳) for tc in [0.10, 0.15, 0.20, 0.25, 0.28, 0.30, 0.32, 0.35]: ok, vcl is_stable(tc) margin (Vcr - vcl) / Vcr * 100.0 print(f{tc:8.2f} {vcl:10.6f} {Vcr:10.6f} {margin:7.2f} {稳定 if ok else 失稳})运行输出大致如下稳态功角: 15.66 deg 故障后SEP: 21.10 deg, UEP: 158.90 deg 临界能量 Vcr 1.098556 p.u. CCT ≈ 0.2901 s tc(s) Vcl(p.u.) Vcr(p.u.) 裕度(%) 时域判稳 0.10 0.021493 1.098556 98.04 稳定 0.15 0.151595 1.098556 86.20 稳定 0.20 0.420885 1.098556 61.69 稳定 0.25 0.845404 1.098556 23.04 稳定 0.28 1.075052 1.098556 2.14 稳定 0.30 1.175376 1.098556 -6.99 失稳 0.32 1.371106 1.098556 -24.81 失稳 0.35 1.671932 1.098556 -52.19 失稳3.3 能量判据与时域判据互验理解裕度符号代码注释里写得很清楚is_stable返回的稳定结论由两段逻辑共同决定。先看切除时刻能量 Vcl 是否超过 Vcr超过直接判失稳没超过则继续积分观察功角是否越过 UEP 并且转速仍高于同步值。扫描表显示tc 0.28 秒时裕度只有 2.14%仍判稳定tc 0.30 秒时裕度变成 -6.99%判失稳。CCT 落在 0.28 和 0.30 之间二分法给出约 0.2901 秒和扫描结果一致。这里值得强调裕度的物理含义margin (Vcr - Vcl) / Vcr × 100%正值表示切除时机组吸收的暂态能量离临界还有富余负值表示已经越界。调度场景里裕度低于 5% 的故障应当重点复核因为模型简化、参数误差都会在临界区被放大裕度高于 30% 的故障基本可以放心初筛掉。3.4 三个参数旋钮H、Pm、X2 对 CCT 的影响实操里最常改的三个旋钮是惯性 H 越大转子“质量”越大同样加速功率下功角爬升越慢CCT 变大。新能源场站等值惯量偏低CCT 会显著缩小这也是高比例新能源系统暂态稳定变差的原因之一。机械功率 Pm 越大加速功率 Pm - Pe 越大单位时间内注入转子的动能越多CCT 变小。调峰机组满载运行时比半载时更容易失稳。故障电抗 X2 越大故障越严重残余转移功率越小故障期间电磁功率越低CCT 越小。三相短路取 X2 为无穷大时Pe 0是理论上最严重的场景本算例取 1.5 p.u.相当于保留了少量互送功率更贴近真实短路计算结果。4. 多机系统的暂态能量函数实现COI 坐标变换、PEBS 与 BCU单机代码跑通后真正的工程场景至少是几十台机、几百个节点。多机系统里没有“无穷大母线”当基准必须先把坐标整体平移掉。4.1 为什么多机系统必须做 COI 坐标变换多机失稳的本质是机群间相对运动而不是所有发电机相对某个参考节点同步加速。把参考点放在无穷大母线时所有机组同向的整体加速运动会混进能量函数导致动能项偏大、势能项失真。惯量中心Center of Inertia, COI坐标把整台系统的“质心”当作观察基准δ_COI (Σ M_i δ_i) / Σ M_iω_COI (Σ M_i ω_i) / Σ M_i变换后的相对量为θ_i δ_i - δ_COIω~_i ω_i - ω_COI。所有公式里的功角和转速都换成相对量。这样处理后整体平移对应的刚体运动从能量函数里剔除剩下的只有机群间相对摆开。4.2 多机能量函数的构成与实用公式在 COI 坐标下多机暂态能量函数仍可写成动能加势能两项但势能不再是单机那个只有两项的表达式而是包含位置和网络储能两部分能量分量表达式物理含义动能 VKE0.5 * Σ M_i * ω~_i²各机相对 COI 的旋转动能恒为正位置势能 VPE_pos-Σ P_mi * (θ_i - θ_is)机械功率沿各机功角方向累积的势能磁场储能 VPE_mag-ΣΣ [E_iE_j G_ij cos(θ_ij) - E_iE_j B_ij sin(θ_ij) - 常数项]网络内电势之间储存的磁场能量修正项负荷、阻尼、控制器的等效做功工程近似时归入势能或忽略表中 G_ij 和 B_ij 是消去网络内节点后的导纳矩阵元素。注意正负号强烈依赖潮流程序的导纳阵约定实际实现时第一步应该拿一个已知 IEEE 节点算例把 VPE 的正确符号和零初始值对齐再上真实电网数据。4.3 PEBS 和 BCU 怎么求多机临界能量多机系统求 Vcr 不再有“π - δs”这种解析解核心手段变成了轨迹上的穿越点搜索PEBS 法保持故障不切除持续积分故障轨迹同时计算势能 VPE 随时间的变化。VPE 第一个“先升后降”的局部峰值就是轨迹穿越势能界面的近似点该点势能即为 Vcr。实现简单但遇到多摆失稳或强非线性时容易把后续更高的峰值误判进去。BCU 法先沿持续故障轨迹找到 VPE 的局部峰即出口点再从出口点出发在梯度系统上迭代求解控制不稳定平衡点 CUEP最后用 CUEP 处的势能作为 Vcr。BCU 对多机强非线性模型更稳健代价是实现复杂度上升。从仿真结果算能量的 Python 示意如下假设已经从 BPA/PSASP 导出了每台机的功角和转速时间序列import numpy as np # 从导出文件读入: 每行一个时刻 # t[:], delta[:, i], omega[:, i], M[i], Pmi[i] def compute_energy_curve(t, delta, omega, M, Pmi, G, B): Mt M.sum() dcoi (M * delta).sum(axis1) / Mt wcoi (M * omega).sum(axis1) / Mt theta delta - dcoi[:, None] wt omega - wcoi[:, None] VKE 0.5 * (M * wt**2).sum(axis1) # 位置势能: 相对故障后SEP theta_s ... # 故障后SEP的COI坐标 VPE_pos -(Pmi * (theta - theta_s)).sum(axis1) # 磁场储能: 对每对节点累加 Gij/Bij 项 # 需要拿到消去网络后的节点导纳矩阵 n M.size VPE_mag np.zeros_like(t) for i in range(n): for j in range(n): dij theta[:, i] - theta[:, j] VPE_mag E[i]*E[j] * (G[i, j]*np.cos(dij) - B[i, j]*np.sin(dij)) VPE VPE_pos - VPE_mag const return VKE, VPE, VKE VPE这段代码的逻辑分三步先按公式求 COI计算动能和位置势能再叠加磁场储能最后把总能量轨道画出来找到故障期间首个 VPE 峰值。整个计算是向量化的几千个节点十几台机的时间序列几秒钟就能跑完。4.4 在线稳定评估怎么用这套东西能量法在线评估的常见做法是“离线建模、在线比较”。离线阶段为候选故障集逐一计算 Vcr生成一张“故障 ID → 临界能量”表在线阶段利用故障后的实测功角轨迹快速算 Vcl查表求裕度。因为 Vcl 只需要切除后很短一段轨迹就能稳定计算不需要把故障后 5 秒跑完计算负担远小于时域仿真。实际系统中通常把它和实时暂态仿真并列运行能量法给出初筛排名仿真组件对排名靠前的故障做精细确认两套结果互相兜底。5. 工程排错与参数敏感性SEP、负荷模型和积分步长最后落在工程实现最容易翻车的三个细节上。5.1 三个高频错误与对策症状根因修正办法稳定裕度整体偏小连明显稳定的故障都接近临界能量函数参考点用了故障前 SEP用故障后网络重新计算 SEP替代故障前 δ0能量曲线出现高频抖动局部峰值无法识别积分器用了变步长自适应算法改用固定步长 RK4 或梯形法步长取 0.5~1 ms与精确时域仿真结论冲突裕度方向相反负荷模型是恒功率/恒电流破坏能量函数保守性将负荷在故障前运行点线性化为恒阻抗归算进导纳阵5.2 模型阶数越高能量函数越难定义发电机用六阶模型带励磁和 PSS 时严格李雅普诺夫函数几乎构造不出来。此时不要硬套经典能量公式。常见做法是“拓展势能法”把励磁和调速器对转子的等效做功近似为势能修正项或直接对 Vcl 与 Vcr 做数据驱动修正。PPT 教案里通常只讲经典模型工程落地时要清楚边界能量法做初筛精细模型交给时域仿真两者结合才可靠。5.3 一个完整的排错顺序遇到能量法结果异常时我一般按这个顺序排查先打印故障后 SEP 和 UEP确认参考点正确再画出 Vk 和 Vp 两条时间曲线看故障期间 Vp 是否出现清晰的单峰结构然后把切除时刻标注在曲线上确认 Vcl 取自切除瞬间而非积分步内的平均值最后把裕度结果与同故障的时域仿真 CCT 对照偏差超过 10% 就回头检查导纳阵正负号和负荷线性化口径。多数能量法“不收敛”“乱报”最后都落在参考点或负荷模型这两个根因上而不是算法本身出了问题。本文还有配套的精品资源点击获取
暂态能量函数法原理与Python实现：从稳定裕度到CCT计算
简介电力系统暂态能量函数法暂态稳定分析PPT学习教案是一份围绕暂态能量函数法的系统教学课件面向电力系统专业学生、考研备考生及从事稳定分析的研究人员旨在帮助读者从能量视角掌握暂态稳定判据与稳定裕度评估方法。压缩包仅含1个PPTX格式文件大小约1.23MB共105页左右结构完整便于按章学习。目前已吸引160人浏览学习。课件从古典力学能量概念出发系统讲解暂态能量函数的基本原理、数学描述与李氏定理结合单机无穷大系统直接法暂态稳定分析逐步过渡到多机系统的特殊问题、两种坐标系表达并重点介绍最近不稳定平衡点法AUEP、等位能线与位能边界、相关不稳定平衡点法RUEP等内容。其中还包含球碗运动模型、功角特性曲线等图文示例可辅助理解临界能量、稳定域和故障切除时间等关键概念。学习后可构建基于李雅普诺夫函数的暂态稳定分析框架适用于课程教学、自学复习及毕业设计参考。1. 暂态能量函数法把暂态稳定分析从“看曲线”变成了“算余量”时域仿真跑完一条功角曲线报告里能写的基本就是“最大摆角 126°未失稳”。可电网调度侧真正要的不是这一个布尔值稳定裕度还剩多少、哪些故障更危险、切机切负荷该留多大反应时间都依赖连续量化结果。暂态能量函数法就是在故障切除时刻停下积分把发电机转子动能、相对位置势能和网络磁场储能折算成一张系统的能量收支表拿它和系统能承受的临界能量直接比较。产出不是曲线而是稳定裕度。裕度从 0.8 跌到 -0.2 时你已经知道系统正在翻过界面而不是等功角甩开后才去事后确认。做在线暂态稳定评估、故障筛选和预防控制的工程师值得把整套方法吃透。下面先用单机无穷大系统把能量函数拆开再给出可复现的 Python 计算最后落回多机系统里 COI 坐标和 PEBS/BCU 的实际工程问题。2. 从等面积法则到李雅普诺夫能量判据暂态能量函数法的三个组成项单机无穷大系统里能量函数有一个特别好懂的退化形态就是本科教科书里的等面积法则。看懂它后面 Vcr、PEBS 这些概念都能对号入座。2.1 等面积法则为什么是单机的暂态能量函数雏形单机无穷大系统的功角摇摆方程为M * d²δ/dt² Pm - Pe其中 M 是惯性常数Pm 是原动机机械功率Pe EV/XΣ·sinδ 是发电机输出的电磁功率。发生短路后 Pe 骤降转子获得净加速功角从 δ0 开始增大保护动作切除故障后网络电抗变大Pe 恢复转子进入减速段。把加速段和减速段写成面积形式加速面积A_acc ∫(Pm - Pe)dδ积分区间从 δ0 到切除时刻 δc减速面积A_dec ∫(Pe - Pm)dδ积分区间从 δc 到最大摇摆角 δmax稳定条件是 A_acc ≤ A_dec。把加速面积理解为转子相对同步转速储存的动能把减速面积理解为制动过程中动能释放的功等面积法则其实就是能量守恒在一个自由度上的投影。单机两群摆动时它完全够用一旦系统变成多机每台机的转速都与其余机不同“整台机组相对无穷大系统”的概念消失等面积法则便不再成立。但能量法的思路保留了下来找一个标量函数让它沿故障后轨迹单调不增再拿初值处能量和稳定域边界能量比较就得到稳定裕度。2.2 暂态能量函数的构成动能 Vk、位置势能 Vp、磁场储能 Vmag单机经典模型下能量函数可以写成显式形式V(δ, ω) 0.5 * M * (ω - ωs)² - Pm * (δ - δs) - Pmax * (cosδ - cosδs)其中 ωs 是同步电角速度δs 是故障后系统的稳定平衡点SEP。这个表达式由三项构成动能项Vk 0.5 * M * (ω - ωs)²转子相对同步转速运动的动能恒为正故障期间转子加速越快这一项越大。位置势能项Vp_pos -Pm * (δ - δs)机械功率沿功角方向累积的能量相当于重力场里“高度”被功角替代。磁场储能项Vp_mag -Pmax * (cosδ - cosδs)来自发电机内电势与网络之间的耦合储能反映电磁功率与功角之间的保守场关系。后两项合称势能 Vp。动能项恒为正势能项在功角摆大时通常上升。把参考点定在哪个平衡点上是个容易踩的坑必须是故障切除后网络拓扑下的稳定平衡点 δs不是故障前的稳态功角 δ0。故障前后网络结构往往不同XΣ 变了δs 也变了直接拿故障前 δ0 当参考能量初值会带系统性偏差。2.3 稳定判据 V(t) Vcr 和李雅普诺夫定理的工程化用法李雅普诺夫定理的工程化读法是如果能找到正定函数 V使其沿故障后轨迹满足 V̇ ≤ 0那么只要切除时刻的能量 Vcl 小于临界能量 Vcr系统就是暂态稳定的反过来 Vcl ≥ Vcr 则失稳。Vcr 在数学上是稳定域边界上的最小能量工程上则按计算量从小到大有三条求法求法计算量精度与特点工程注意点直接法UEP 法中对单机精确多机需要求控制不稳定平衡点收敛性依赖初值求 UEP 的迭代易发散通常用牛顿-拉夫逊法搭配同伦法PEBS 法小近似值势能界面穿越点偏乐观沿持续故障轨迹找势能峰值第一个峰值才有意义BCU 法大精度最好对多机强非线性模型稳健先找出口点再沿梯度系统求 CUEP实现复杂度高这里有个反直觉结论能量判据不替代时域仿真而是做互补。它建立在理想化模型上电网建模越细励磁、PSS、负荷特性、饱和严格意义上的李雅普诺夫函数越难构造。工程上通常用能量法做快速筛选用精确暂态仿真做最终确认两种手段对同一批故障互相校核。3. 用 Python 算单机无穷大系统的暂态稳定裕度CCT 二分搜索回到可复现的最小实现。一台汽轮发电机经双回线并入无穷大母线0 秒时机端三相对地短路tc 时刻保护切除故障线路只剩一回线运行。目标求临界切除时间 CCT并同时输出能量裕度。3.1 模型参数与摇摆方程模型采用二阶经典摇摆方程Pe 与功角的关系随故障阶段切换。三个电抗分别对应三种网络拓扑正常态总电抗 X1 0.30 p.u.双回线并列运行故障态总电抗 X2 1.50 p.u.代表三相短路后残压较低、仍有残余转移功率的场景切除后总电抗 X3 0.40 p.u.只剩一回线机械功率 Pm 0.9 p.u.内电势 E 和无穷大母线电压均取 1.0 p.u.。惯性时间常数 H 3.5 秒换算成摇摆方程里的惯性常数需要除以同步电角速度M 2H / ωsωs 2π·50 ≈ 314.16 rad/s。3.2 完整可运行的 RK4 积分与能量计算代码下面的代码不依赖 scipy只用标准库 math复制即可运行。核心是四阶龙格-库塔积分器外加能量函数计算与 CCT 二分搜索。import math # ---------- 系统参数 ---------- f0 50.0 omega_s 2.0 * math.pi * f0 # 同步电角速度, rad/s H 3.5 # 惯性时间常数, s M 2.0 * H / omega_s # 摇摆方程惯性常数, p.u. Pm 0.9 # 机械功率, p.u. E1 1.0 # 发电机内电势, p.u. Vb 1.0 # 无穷大母线电压, p.u. X1 0.30 # 正常态总电抗 X2 1.50 # 故障态总电抗三相短路残余 X3 0.40 # 切除后总电抗 Pmax_normal E1 * Vb / X1 Pmax_fault E1 * Vb / X2 Pmax_post E1 * Vb / X3 # 故障前与故障后的稳定平衡点功角 delta_s math.asin(min(1.0, Pm / Pmax_normal)) delta_sp math.asin(min(1.0, Pm / Pmax_post)) delta_uep math.pi - delta_sp # 控制UEP单机系统为π-δs # 临界能量 Vcr势能在 UEP 处的取值 Vcr (-Pm * (delta_uep - delta_sp) - Pmax_post * (math.cos(delta_uep) - math.cos(delta_sp))) print(f稳态功角: {math.degrees(delta_s):.2f} deg) print(f故障后SEP: {math.degrees(delta_sp):.2f} deg, UEP: {math.degrees(delta_uep):.2f} deg) print(f临界能量 Vcr {Vcr:.6f} p.u.) # ---------- RK4 积分器 ---------- def derivative(state, pm, pmax): delta, omega state Pe pmax * math.sin(delta) return [omega - omega_s, (pm - Pe) / M] def rk4_step(state, dt, pm, pmax): k1 derivative(state, pm, pmax) s2 [state[0] dt/2*k1[0], state[1] dt/2*k1[1]] k2 derivative(s2, pm, pmax) s3 [state[0] dt/2*k2[0], state[1] dt/2*k2[1]] k3 derivative(s3, pm, pmax) s4 [state[0] dt*k3[0], state[1] dt*k3[1]] k4 derivative(s4, pm, pmax) d_delta (k1[0] 2*k2[0] 2*k3[0] k4[0]) / 6.0 d_omega (k1[1] 2*k2[1] 2*k3[1] k4[1]) / 6.0 return [state[0] dt*d_delta, state[1] dt*d_omega] # ---------- 能量函数 ---------- def energy(delta, omega): vk 0.5 * M * (omega - omega_s) ** 2 vp (-Pm * (delta - delta_sp) - Pmax_post * (math.cos(delta) - math.cos(delta_sp))) return vk vp, vk, vp # ---------- 给定切除时间判断是否稳定 ---------- def is_stable(tc, t_max2.0, dt0.0005): state [delta_s, omega_s] t 0.0 # 故障期间积分到 tc while t tc - 1e-12: h min(dt, tc - t) state rk4_step(state, h, Pm, Pmax_fault) t h vcl, _, _ energy(state[0], state[1]) if vcl Vcr: return False, vcl # 能量已经越过临界直接判不稳 # 切除后继续积分观察功角是否回摆或越限 while t t_max: state rk4_step(state, dt, Pm, Pmax_post) t dt d state[0] # 功角越过UEP且仍在加速 失稳 if d delta_uep and state[1] omega_s: return False, vcl # 功角回摆到SEP附近且速度接近同步 稳定 if d delta_sp 0.05 and state[1] omega_s 0.1: return True, vcl return True, vcl # ---------- 二分搜索 CCT ---------- lo, hi 0.01, 1.0 for _ in range(40): mid (lo hi) / 2.0 ok, _ is_stable(mid) if ok: lo mid else: hi mid cct (lo hi) / 2.0 print(fCCT ≈ {cct:.4f} s) # ---------- 扫描切除时间输出能量表 ---------- print(\n tc(s) Vcl(p.u.) Vcr(p.u.) 裕度(%) 时域判稳) for tc in [0.10, 0.15, 0.20, 0.25, 0.28, 0.30, 0.32, 0.35]: ok, vcl is_stable(tc) margin (Vcr - vcl) / Vcr * 100.0 print(f{tc:8.2f} {vcl:10.6f} {Vcr:10.6f} {margin:7.2f} {稳定 if ok else 失稳})运行输出大致如下稳态功角: 15.66 deg 故障后SEP: 21.10 deg, UEP: 158.90 deg 临界能量 Vcr 1.098556 p.u. CCT ≈ 0.2901 s tc(s) Vcl(p.u.) Vcr(p.u.) 裕度(%) 时域判稳 0.10 0.021493 1.098556 98.04 稳定 0.15 0.151595 1.098556 86.20 稳定 0.20 0.420885 1.098556 61.69 稳定 0.25 0.845404 1.098556 23.04 稳定 0.28 1.075052 1.098556 2.14 稳定 0.30 1.175376 1.098556 -6.99 失稳 0.32 1.371106 1.098556 -24.81 失稳 0.35 1.671932 1.098556 -52.19 失稳3.3 能量判据与时域判据互验理解裕度符号代码注释里写得很清楚is_stable返回的稳定结论由两段逻辑共同决定。先看切除时刻能量 Vcl 是否超过 Vcr超过直接判失稳没超过则继续积分观察功角是否越过 UEP 并且转速仍高于同步值。扫描表显示tc 0.28 秒时裕度只有 2.14%仍判稳定tc 0.30 秒时裕度变成 -6.99%判失稳。CCT 落在 0.28 和 0.30 之间二分法给出约 0.2901 秒和扫描结果一致。这里值得强调裕度的物理含义margin (Vcr - Vcl) / Vcr × 100%正值表示切除时机组吸收的暂态能量离临界还有富余负值表示已经越界。调度场景里裕度低于 5% 的故障应当重点复核因为模型简化、参数误差都会在临界区被放大裕度高于 30% 的故障基本可以放心初筛掉。3.4 三个参数旋钮H、Pm、X2 对 CCT 的影响实操里最常改的三个旋钮是惯性 H 越大转子“质量”越大同样加速功率下功角爬升越慢CCT 变大。新能源场站等值惯量偏低CCT 会显著缩小这也是高比例新能源系统暂态稳定变差的原因之一。机械功率 Pm 越大加速功率 Pm - Pe 越大单位时间内注入转子的动能越多CCT 变小。调峰机组满载运行时比半载时更容易失稳。故障电抗 X2 越大故障越严重残余转移功率越小故障期间电磁功率越低CCT 越小。三相短路取 X2 为无穷大时Pe 0是理论上最严重的场景本算例取 1.5 p.u.相当于保留了少量互送功率更贴近真实短路计算结果。4. 多机系统的暂态能量函数实现COI 坐标变换、PEBS 与 BCU单机代码跑通后真正的工程场景至少是几十台机、几百个节点。多机系统里没有“无穷大母线”当基准必须先把坐标整体平移掉。4.1 为什么多机系统必须做 COI 坐标变换多机失稳的本质是机群间相对运动而不是所有发电机相对某个参考节点同步加速。把参考点放在无穷大母线时所有机组同向的整体加速运动会混进能量函数导致动能项偏大、势能项失真。惯量中心Center of Inertia, COI坐标把整台系统的“质心”当作观察基准δ_COI (Σ M_i δ_i) / Σ M_iω_COI (Σ M_i ω_i) / Σ M_i变换后的相对量为θ_i δ_i - δ_COIω~_i ω_i - ω_COI。所有公式里的功角和转速都换成相对量。这样处理后整体平移对应的刚体运动从能量函数里剔除剩下的只有机群间相对摆开。4.2 多机能量函数的构成与实用公式在 COI 坐标下多机暂态能量函数仍可写成动能加势能两项但势能不再是单机那个只有两项的表达式而是包含位置和网络储能两部分能量分量表达式物理含义动能 VKE0.5 * Σ M_i * ω~_i²各机相对 COI 的旋转动能恒为正位置势能 VPE_pos-Σ P_mi * (θ_i - θ_is)机械功率沿各机功角方向累积的势能磁场储能 VPE_mag-ΣΣ [E_iE_j G_ij cos(θ_ij) - E_iE_j B_ij sin(θ_ij) - 常数项]网络内电势之间储存的磁场能量修正项负荷、阻尼、控制器的等效做功工程近似时归入势能或忽略表中 G_ij 和 B_ij 是消去网络内节点后的导纳矩阵元素。注意正负号强烈依赖潮流程序的导纳阵约定实际实现时第一步应该拿一个已知 IEEE 节点算例把 VPE 的正确符号和零初始值对齐再上真实电网数据。4.3 PEBS 和 BCU 怎么求多机临界能量多机系统求 Vcr 不再有“π - δs”这种解析解核心手段变成了轨迹上的穿越点搜索PEBS 法保持故障不切除持续积分故障轨迹同时计算势能 VPE 随时间的变化。VPE 第一个“先升后降”的局部峰值就是轨迹穿越势能界面的近似点该点势能即为 Vcr。实现简单但遇到多摆失稳或强非线性时容易把后续更高的峰值误判进去。BCU 法先沿持续故障轨迹找到 VPE 的局部峰即出口点再从出口点出发在梯度系统上迭代求解控制不稳定平衡点 CUEP最后用 CUEP 处的势能作为 Vcr。BCU 对多机强非线性模型更稳健代价是实现复杂度上升。从仿真结果算能量的 Python 示意如下假设已经从 BPA/PSASP 导出了每台机的功角和转速时间序列import numpy as np # 从导出文件读入: 每行一个时刻 # t[:], delta[:, i], omega[:, i], M[i], Pmi[i] def compute_energy_curve(t, delta, omega, M, Pmi, G, B): Mt M.sum() dcoi (M * delta).sum(axis1) / Mt wcoi (M * omega).sum(axis1) / Mt theta delta - dcoi[:, None] wt omega - wcoi[:, None] VKE 0.5 * (M * wt**2).sum(axis1) # 位置势能: 相对故障后SEP theta_s ... # 故障后SEP的COI坐标 VPE_pos -(Pmi * (theta - theta_s)).sum(axis1) # 磁场储能: 对每对节点累加 Gij/Bij 项 # 需要拿到消去网络后的节点导纳矩阵 n M.size VPE_mag np.zeros_like(t) for i in range(n): for j in range(n): dij theta[:, i] - theta[:, j] VPE_mag E[i]*E[j] * (G[i, j]*np.cos(dij) - B[i, j]*np.sin(dij)) VPE VPE_pos - VPE_mag const return VKE, VPE, VKE VPE这段代码的逻辑分三步先按公式求 COI计算动能和位置势能再叠加磁场储能最后把总能量轨道画出来找到故障期间首个 VPE 峰值。整个计算是向量化的几千个节点十几台机的时间序列几秒钟就能跑完。4.4 在线稳定评估怎么用这套东西能量法在线评估的常见做法是“离线建模、在线比较”。离线阶段为候选故障集逐一计算 Vcr生成一张“故障 ID → 临界能量”表在线阶段利用故障后的实测功角轨迹快速算 Vcl查表求裕度。因为 Vcl 只需要切除后很短一段轨迹就能稳定计算不需要把故障后 5 秒跑完计算负担远小于时域仿真。实际系统中通常把它和实时暂态仿真并列运行能量法给出初筛排名仿真组件对排名靠前的故障做精细确认两套结果互相兜底。5. 工程排错与参数敏感性SEP、负荷模型和积分步长最后落在工程实现最容易翻车的三个细节上。5.1 三个高频错误与对策症状根因修正办法稳定裕度整体偏小连明显稳定的故障都接近临界能量函数参考点用了故障前 SEP用故障后网络重新计算 SEP替代故障前 δ0能量曲线出现高频抖动局部峰值无法识别积分器用了变步长自适应算法改用固定步长 RK4 或梯形法步长取 0.5~1 ms与精确时域仿真结论冲突裕度方向相反负荷模型是恒功率/恒电流破坏能量函数保守性将负荷在故障前运行点线性化为恒阻抗归算进导纳阵5.2 模型阶数越高能量函数越难定义发电机用六阶模型带励磁和 PSS 时严格李雅普诺夫函数几乎构造不出来。此时不要硬套经典能量公式。常见做法是“拓展势能法”把励磁和调速器对转子的等效做功近似为势能修正项或直接对 Vcl 与 Vcr 做数据驱动修正。PPT 教案里通常只讲经典模型工程落地时要清楚边界能量法做初筛精细模型交给时域仿真两者结合才可靠。5.3 一个完整的排错顺序遇到能量法结果异常时我一般按这个顺序排查先打印故障后 SEP 和 UEP确认参考点正确再画出 Vk 和 Vp 两条时间曲线看故障期间 Vp 是否出现清晰的单峰结构然后把切除时刻标注在曲线上确认 Vcl 取自切除瞬间而非积分步内的平均值最后把裕度结果与同故障的时域仿真 CCT 对照偏差超过 10% 就回头检查导纳阵正负号和负荷线性化口径。多数能量法“不收敛”“乱报”最后都落在参考点或负荷模型这两个根因上而不是算法本身出了问题。本文还有配套的精品资源点击获取