简介：一份以Python代码为驱动的电力系统暂态稳定性分析资料，围绕三机九节点系统完整覆盖“建模—潮流—仿真—评估”主线：发电机经典二阶模型、负荷恒阻抗建模、导纳矩阵构建、牛顿-拉夫逊法潮流计算、改进欧拉法故障过程模拟，以及临界切除时间计算。面向具备电力系统基础、有一定编程经验的研究人员、工程师和高校学生，可用于实际工程校验与教学实践。压缩包共1个docx文档，大小约54KB，内含可运行代码、逐段中文解释、结果可视化图表和Simulink验证说明，便于对照调试。目前已有96人学习该资源。读者能借此掌握暂态稳定判据的具体实现，理解阻尼、故障切除时间对功角曲线的影响，并参考文中的智能算法、实时稳定判别等优化方向进行扩展。
1. 暂态稳定性仿真不是“跑通代码”：为什么是三机九节点
很多初学者拿到“电力稳定性分析”的题目，第一反应是找个现成仿真工具，把模型一拖、按钮一按，看功角曲线不发散就算完事。但真正到要写报告、要解释临界切除时间的时候，才发现自己手里是个黑匣子。三机九节点系统的价值就在于：它刚好复杂到能承载暂态稳定性仿真的完整链路——系统建模、网络化简、潮流初值、时域积分、故障设置、临界切除时间搜索，又简单到每一行代码都能和教材公式对上。这套体系吃透了，换39节点、IEEE 118节点只是改数据和矩阵规模的事。适合电力系统专业学生、刚入门的继保或稳控工程师，以及想把“暂态稳定”从概念变成可复现结果的研究者。
2. 系统建模：从三机九节点参数到可计算的导纳矩阵
2.1 参数别凭感觉填：发电机、变压器、线路数据的标幺化
暂态稳定性仿真里，单位错误是最隐蔽的翻车点。所有模型参数统一用标幺值，基准功率一般是100 MVA，电压基准取各电压等级的额定值。发电机侧需要两个核心参数：暂态电抗xd'和惯性时间常数H。H的单位是秒，物理含义是转子在额定转矩下从静止加速到额定转速所需的时间。H越大，转子越“沉”，功角变化越慢，系统越不容易失稳。
这个经典测试系统的数据在工程里已经用得比较固定，我一般直接写成一个 Python 字典：
import numpy as np # ---- 三机九节点系统参数（标幺值，基准功率 100 MVA） ---- gen_data = { 1: {"bus": 1, "xd": 0.0608, "H": 23.64}, 2: {"bus": 2, "xd": 0.1198, "H": 6.4}, 3: {"bus": 3, "xd": 0.1813, "H": 3.01}, } # 支路：母线i, 母线j, 串联电阻, 电抗, 对地电纳 branches = [ (1, 4, 0.0000, 0.0576, 0.000), (2, 7, 0.0000, 0.0625, 0.000), (3, 9, 0.0000, 0.0586, 0.000), (4, 5, 0.0100, 0.0850, 0.088), (4, 6, 0.0170, 0.0920, 0.079), (5, 7, 0.0320, 0.1610, 0.153), (6, 9, 0.0390, 0.1700, 0.179), (7, 8, 0.0085, 0.0720, 0.0745), (8, 9, 0.0119, 0.1008, 0.1045), ] # 负荷：母线 -> (有功P, 无功Q)，单位 pu load_data = {5: (1.25, 0.50), 6: (0.90, 0.30), 8: (1.00, 0.35)}第一条到第三条支路是升压变压器，电阻为 0，只有漏抗。第四条以后才是输电线路，用 π 形等值电路，对地电纳b/2加在线路两端。注意线路 5-7 和 6-9 的电抗比较大，这两条是主要功率传输通道，后面做故障仿真时最容易让系统失稳。
2.2 构建9母线导纳矩阵：负荷用恒阻抗还是恒功率
导纳矩阵是暂态仿真的地基。构建时有个选择：负荷到底按恒功率还是恒阻抗处理。在经典暂态稳定仿真里，工程上最常见的做法是把负荷近似成恒定阻抗，因为它可以直接并入导纳矩阵，避免在每个时步反复迭代潮流。代价是负荷吸收功率会随电压平方变化，和实际恒功率负荷有偏差。本文统一用恒阻抗模型，代码里把这个选择固化下来。
def build_Y9(fault=None, trip=None): """构建 9 母线导纳矩阵。 fault=(i,j): 线路 i-j 中点发生三相短路。 trip=(i,j) : 故障后跳开线路 i-j。 """ Y = np.zeros((9, 9), dtype=complex) def is_target(a, b, target): if target is None: return False return (a == target[0] and b == target[1]) or (a == target[1] and b == target[0]) for k, m, r, x, b in branches: if is_target(k, m, trip): # 跳闸线路直接不加入网络 continue y = 1.0 / (r + 1j * x) Y[k - 1, k - 1] += y + 1j * b / 2 Y[m - 1, m - 1] += y + 1j * b / 2 Y[k - 1, m - 1] -= y Y[m - 1, k - 1] -= y if is_target(k, m, fault): # 三相短路：把线路中点看成接地母线 # 等效做法是两侧母线各并入一条半段线路的对地导纳 y_ground = 2.0 / (r + 1j * x) Y[k - 1, k - 1] += y_ground Y[m - 1, m - 1] += y_ground # 负荷按恒定阻抗并入导纳阵，取额定电压 1.0 pu 估算 for bus, (P, Q) in load_data.items(): Y[bus - 1, bus - 1] += (P - 1j * Q) return Y这里的故障等效要重点说：线路中点三相短路，故障点电压近似为零，原来的一条线路被拆成两段，每一段末端都接地。所以程序先把故障线路从导纳矩阵里拿掉，再在两端母线上各自并联一个2/(r+jx)的对地导纳。这是把短路故障“嵌入”网络矩阵的常见技巧，很多教材里不会写这么细，但自己实现时最容易在这里错。
2.3 发电机内节点与Kron约化：为什么最后只需要3×3矩阵
有了9母线导纳矩阵还不行，因为暂态仿真里发电机不是一个简单功率注入源，而是一个“内电势 + 暂态电抗”的戴维南等值电路。每台发电机要增加一个内节点，通过jxd'连接到机端母线。这样网络从9节点变成12节点，然后再把非发电机节点全部消掉，得到3×3的约化导纳矩阵。Kron约化就是分块矩阵消去法。
def make_Yred(fault=None, trip=None): """构造只含发电机内节点的 3x3 约化导纳矩阵。""" Y9 = build_Y9(fault, trip) Y = np.zeros((12, 12), dtype=complex) Y[:9, :9] = Y9 # 发电机内节点编号：9, 10, 11（对应发电机 1, 2, 3） for g in range(1, 4): bus_idx = gen_data[g]["bus"] - 1 inner_idx = 8 + g yg = 1.0 / (1j * gen_data[g]["xd"]) Y[bus_idx, bus_idx] += yg Y[inner_idx, inner_idx] += yg Y[bus_idx, inner_idx] -= yg Y[inner_idx, bus_idx] -= yg keep = [9, 10, 11] remove = [i for i in range(12) if i not in keep] Ykk = Y[np.ix_(keep, keep)] Ykr = Y[np.ix_(keep, remove)] Yrk = Y[np.ix_(remove, keep)] Yrr = Y[np.ix_(remove, remove)] Yred = Ykk - Ykr @ np.linalg.inv(Yrr) @ Yrk return Yred约化之后，每台发电机的电功率可以直接用Pe_i = Re(E_i * conj(I_i))计算，其中I = Yred @ E。这里E_i是复数内电势，幅值由初始潮流决定，相角就是功角delta。这一步做完，网络方程就退化成3个节点，后面时域积分每步只算3×3的复数乘法，速度快得多。
3. 时域仿真实现：让摇摆方程真正跑起来
3.1 先求稳态运行点：数值潮流反推暂态初值
暂态仿真的初值不是随便给的。先要做一次潮流计算，得到各母线电压和发电机注入功率，再反推每台发电机的内电势幅值和初始功角。这个系统只有14个未知量，我用scipy.optimize.root直接解潮流方程，省去手写雅可比矩阵的繁琐。
from scipy.optimize import root def solve_initial_state(): """潮流求解：返回 Vph（9母线复数电压）和 Pm（发电机机械功率）。""" Y9 = build_Y9() def pf_residual(x): theta = np.zeros(9) theta[1:] = x[0:8] V = np.ones(9) V[0] = 1.04 # 平衡机母线电压 V[1] = 1.025 # 发电机2母线电压 V[2] = 1.025 # 发电机3母线电压 V[3:] = x[8:14] # PQ母线电压幅值 E = V * np.exp(1j * theta) S = E * np.conj(Y9 @ E) P = S.real Q = S.imag P_target = np.zeros(9) P_target[1] = 1.63 # 发电机2有功注入 P_target[2] = 0.85 # 发电机3有功注入 res = [] for i in range(1, 9): res.append(P[i] - P_target[i]) for i in range(3, 9): res.append(Q[i]) # PQ母线无功平衡 return np.array(res) x0 = np.zeros(14) x0[8:] = 1.0 sol = root(pf_residual, x0, method="hybr") if not sol.success: raise RuntimeError(f"潮流不收敛: {sol.message}") theta = np.zeros(9) theta[1:] = sol.x[0:8] V = np.ones(9) V[0] = 1.04 V[1] = 1.025 V[2] = 1.025 V[3:] = sol.x[8:14] return V * np.exp(1j * theta)潮流解出来之后，暂态初值就定了：Pm取三台发电机注入有功，内电势按E_g = V + jxd' * I_g反推，这里I_g = conj(S_g / V)。注意Pm不能用节点注入功率代替，必须取发电机母线上的复功率，因为变压器支路还有无功损耗。
3.2 摇摆方程怎么离散：RHS、积分器与步长参数
三机九节点的经典暂态模型用二阶摇摆方程描述，每台发电机两个状态量：功角delta和标幺转速omega。时间常数上，H是秒，omega_s = 2*pi*60 ≈ 377。微分方程为：
d(delta)/dt = omega_s * (omega - 1)d(omega)/dt = (Pm - Pe - D*(omega-1)) / (2*H) * omega_s
Pe是电磁功率，和导纳矩阵、内电势相角直接相关。积分器我一般选RK45，但必须限制最大步长，否则在故障切除瞬间容易跳过暂态峰值。
from scipy.integrate import solve_ivp omega_s = 2 * np.pi * 60 def make_rhs(Yred, Pm, E0, H, D=0.0): """构造摇摆方程右端函数。""" def rhs(t, z): delta = z[0:3] omega = z[3:6] E = E0 * np.exp(1j * delta) # 内电势相量 I = Yred @ E # 网络方程 Pe = np.real(E * np.conj(I)) # 电磁功率 ddelta = omega_s * (omega - 1.0) domega = (Pm - Pe - D * (omega - 1.0)) / (2.0 * H) * omega_s return np.concatenate([ddelta, domega]) return rhs def simulate(tcl, T=2.0, D=0.0): """tcl：故障切除时间；T：仿真时长。""" Vph = solve_initial_state() Y9 = build_Y9() S_net = Vph * np.conj(Y9 @ Vph) Pm = np.zeros(3) E0 = np.zeros(3, dtype=complex) H = np.zeros(3) for g in range(1, 4): idx = gen_data[g]["bus"] - 1 Sg = S_net[idx] Ig = np.conj(Sg / Vph[idx]) E0[g - 1] = Vph[idx] + 1j * gen_data[g]["xd"] * Ig Pm[g - 1] = Sg.real H[g - 1] = gen_data[g]["H"] delta0 = np.angle(E0) omega0 = np.ones(3) z0 = np.concatenate([delta0, omega0]) Yf = make_Yred(fault=(5, 7), trip=None) # 故障中网络 Yp = make_Yred(fault=None, trip=(5, 7)) # 故障后跳开5-7 def rhs(t, z): Yred = Yf if t < tcl else Yp return make_rhs(Yred, Pm, E0, H, D)(t, z) sol = solve_ivp(rhs, [0, T], z0, method="RK45", max_step=0.005, rtol=1e-6, atol=1e-8) return sol, E0, H故障设置在这里是线路 5-7 中点三相短路，故障切除的同时跳开这条线路。tcl之前用Yf，之后用Yp。solve_ivp内部虽然自适应步长，但max_step=0.005s必须给，否则切除瞬间的数值跳变可能被大步长跨过去。
3.3 结果判断：相对功角和失稳判据
仿真结束不能只看曲线“大概没飞”，要有一个可量化的失稳判据。工程里最常用的是相对功角判据：以惯量加权的惯性中心角为基准，看任意一台发电机相对中心角是否超过180度。超过180度意味着两台发电机之间已经失去同步，系统基本判失稳。
def is_stable(sol, H): """基于惯性中心相对功角判据判断稳定性。""" z = sol.y delta = z[0:3] H_sum = np.sum(H) delta_coi = (H @ delta) / H_sum rel = delta - delta_coi return np.max(np.abs(rel)) < np.pi这个判据比“最大功角差小于某值”更稳。因为三机系统里功角整体偏移是允许的，只有相对运动表征同步能力。如果你想看得更细致，可以打印最后一秒的相对功角变化率，如果持续朝一个方向增大，即使还没到180度，潜在失稳也已经发生。
4. 临界切除时间确定：二分搜索与评估口径
4.1 为什么要二分而不是“跑一堆 tcl”
临界切除时间（CCT）是暂态稳定性仿真里最有工程价值的输出：故障必须在这么短时间内被切除，否则系统失稳。直接枚举所有可能切除时间当然可行，但每跑一次仿真都要做潮流和数值积分，效率太低。CCT对切除时间的变化是单调的：切得越晚越容易失稳，所以可以用二分搜索，把区间稳定地缩小到几十毫秒以内。
4.2 CCT搜索代码与收敛边界
def find_cct(tlo=0.01, thi=0.5, tol=0.001): """二分搜索临界切除时间，精度 1ms。""" while thi - tlo > tol: tmid = 0.5 * (tlo + thi) sol, _, H = simulate(tmid) if is_stable(sol, H): tlo = tmid else: thi = tmid return 0.5 * (tlo + thi)注意我把下界设成 0.01s 而不是 0。因为tcl=0在物理上代表故障一开始就被切除，可以直接用故障后网络跑，不需要构造故障矩阵；让二分从很小的值开始更稳妥。搜索精度tol=0.001s对工程报告够用，如果要做更精细的稳控策略整定，可以提到 0.0005s，代价是搜索次数增加约 10 次仿真。
这个场景下线路 5-7 三相短路的临界切除时间通常在 0.2s 附近。你跑出来的具体值会因为负荷导纳的估计方式、阻尼系数和积分容差有少量浮动。二分法给出的是一个置信区间：tlo是能稳定恢复的最晚切除时间，thi是系统开始失稳的最早切除时间，最后返回二者中点。
4.3 CCT结果与影响参数
| 参数变化 | 对CCT的影响 | 原因 |
|---|---|---|
| 阻尼系数 D 从 0 加到 3 | CCT 可能增加 10~30ms | 阻尼帮助衰减功角振荡，延缓失稳 |
| 负荷从恒阻抗改为恒功率 | CCT 可能减小 30~50ms | 恒功率负荷在低压时吸收功率不减，加重功率缺额 |
| 故障点从 5-7 移到 8-9 | CCT 明显变化 | 不同通道传输功率不同，故障严重度不同 |
| 发电机 H 减半 | CCT 大幅下降 | 转子惯量小，功角加速度快 |
做过暂态稳定整定的工程师都会有这个感觉：CCT 不是一个精确到小数点后三位的物理常数，而是和负荷模型、阻尼假设强相关的评估口径。写报告时最好把工况条件写清楚，否则单给一个 CCT 数字没有复现意义。
5. 避坑与常见问题：五条拿命换来的暂态仿真经验
5.1 母线编号从1起、numpy从0起，错位半天查不出
现象：潮流不收敛，或者收敛结果里发电机母线电压奇低。 原因：数据文件里母线编号从 1 开始，build_Y9里虽然做了k-1转换，但如果在solve_initial_state里直接拿母线编号当索引，就会错位。 解决：所有取网络元素的地方统一走bus_idx = bus - 1。写代码时不要省这一步，宁可多写一个变量也不要直接Vph[gen_data[g]["bus"]]。
5.2 发电机内电抗没加进导纳阵，暂态功角直接“飞”
现象：故障前潮流完全正常，但一进入暂态仿真，功角曲线立刻以接近线性的速度飞散。 原因：make_Yred里只取了9母线网络，发电机被当成理想电压源接到母线，没有经过jxd'。这样网络戴维南等值阻抗偏小，电磁功率传输能力虚高，故障期间减速面积被低估。 解决：检查约化矩阵是不是3×3复数矩阵，且对角线实部不为正得离谱。一个简单自检：把故障前Yred代入稳态Pe，如果和潮流算出的Pm差超过 5%，说明网络等值有问题。
5.3 故障与跳闸用同一套导纳阵，得到假CCT
现象：CCT 比文献值大一倍，怎么调都不对。 原因：故障期间用故障后导纳阵跑，相当于故障在线路中点但又被瞬间隔离，物理上少了故障持续期间的能量积累。 解决：simulate里必须同时准备Yf和Yp两套约化矩阵，以t < tcl为界切换。我一般还会把故障支路打印出来，确认fault=(5,7)和trip=(5,7)在build_Y9里确实走了不同分支。
5.4 积分步长与失稳判据配合不好，出现假稳定
现象：某次把max_step从 0.005 改成 0.05，CCT 增大明显。 原因：RK45 在故障切除瞬间有快速暂态，步长太粗会漏掉功角峰值，让失稳系统看起来“稳定”了。 解决：max_step不要大于 0.01s，最好保持 0.005s。同时把rtol保持在 1e-6 以下。如果做批量 CCT 搜索，固定max_step后还要固定求解器，不同求解器之间的结果对比要特别小心。
5.5 初始功率不对齐，仿真曲线起点就在漂
现象：故障前功角曲线不是平的，而是从t=0就开始缓慢变化。 原因：Pm用的是给定功率注入，而Pe是由潮流网络算出来的，两者在初始点不相等，系统一开始就有加速功率。 解决：Pm必须取潮流解算出的发电机注入有功，而不是数据表里的目标值。尤其平衡机的Pm一定是潮流算出来的，因为平衡机有功是自由的，数据表里根本没有预设值。这个细节很多参考代码都含糊带过，但它是暂态仿真初值自洽的底线。
6. 进阶与验证：把CCT从单个数字变成可信区间
仿真做完了，如果只交一个 CCT 数字，答辩或审查时很容易被动。我习惯再做两个交叉验证：第一，换一条故障线路重跑整套流程，对比不同通道的临界切除时间变化趋势，确认系统薄弱环节和理论分析一致；第二，检查D=0和D=3两组阻尼下的 CCT 差异，如果差异超过 30ms 就顺便写进报告，说明结果对阻尼模型的敏感性，有理有据。
一个更轻量的自检方法是换积分器：把method="RK45"改成method="LSODA"或"BDF"，重新跑一次二分搜索。如果两个积分器给出的 CCT 差超过 10ms，说明你的max_step或容差设置不合理，而不是积分器本身的问题。这比对着曲线“看感觉”靠谱得多。
# 验证脚本：打印稳定/失稳边界附近的相对功角最大值 for tcl in [0.20, 0.21, 0.22, 0.23, 0.24]: sol, _, H = simulate(tcl) z = sol.y delta = z[0:3] rel = delta - (H @ delta) / np.sum(H) max_rel = np.max(np.abs(rel), axis=0) stable = "稳定" if np.max(max_rel) < np.pi else "失稳" print(f"tcl={tcl:.2f}s 最大相对功角={max_rel[-1]*180/np.pi:.1f}deg {stable}")这段脚本会把二分搜索的边界摊开给你看，比黑盒的“稳定/失稳”结果更有说服力。我自己每次跑完 CCT，都会把这个边界打印出来留档，因为后面调负荷模型或改造网络时，回来比较“哪些工况进入了失稳区”比只盯着一个临界值有用得多。
暂态稳定性仿真的坑，大多数不在算法理论，而在网络矩阵的构造和初值自洽。把这两点守住，三机九节点这套流程就能稳稳迁移到更大系统上。希望帮到你。
本文还有配套的精品资源，点击获取
