简介面向电力系统暂态稳定分析与多机系统仿真方向的初学者与研究人员这份资源以三机系统线路AB段首端发生两相短路接地故障为场景演示了从故障发生到0.1秒后切除故障线路的完整时域仿真流程。压缩包内共5个m文件大小仅4KB包含主程序main.m、改进欧拉法数值积分Improved_euler.m、核心计算逻辑calculate.m、系统参数数据data1.m以及案例读取函数readcase.m便于快速理解各模块分工。目前已有395人学习下载。通过运行这些脚本可获得发电机转速、电压、电流等关键量随时间变化的曲线直观观察暂态过程中转子角度摆动与系统恢复能力进而判断多机系统能否重新回到稳定运行状态。对于电力系统继电保护策略评估、稳定性机理研究及MATLAB仿真实践均具有参考价值尤其适合作为课程设计或科研入门的起步模板。1. 短路故障后的暂态稳定为什么只能靠时域仿真做电力系统分析的工程师十有八九都遇过这样的场景送出线路发生三相短路保护在 120 ms 内切除故障运行人员想知道系统会不会失稳。发电机一堆线路几十条负荷还在波动任何解析公式都算不出答案。多机系统短路故障后的暂态稳定分析本质上是在回答一个非常具体的问题故障切除后每台发电机的转子能不能重新回到同步转速附近彼此功角差能不能收敛。时域仿真就是把这场“机械惯性 电磁功率”的博弈一步步算给你看。它从潮流断面出发逐时刻求解发电机的转子运动方程和电网的代数方程直到我们能够判断系统稳住还是失稳。单机无穷大系统的等面积法则在多机系统里基本失效因为每台机的电磁功率会受到其他机组功角摆开的非线性影响这种影响没有闭合表达式只能靠数值积分推进。适合它的用户很明确——调度机构做安全稳定校核的工程师、设计院做送出工程论证的规划人员、研究暂态稳定控制策略的学生。门槛不算低但上手路径非常固定选模型、搭网络、设故障、跑仿真、判结果每一步都有章可循。2. 多机系统时域仿真的建模从摇摆方程到增广导纳矩阵时域仿真的第一步不是写代码而是把“多机系统”翻译成一套可以计算的数学结构。这里面有两个核心问题发电机怎么描述网络和故障怎么进入方程。建模深度决定了仿真结果的可用范围也决定了代码的复杂程度。2.1 发电机模型选到什么程度先古典模型再逐步加细节暂态稳定最常用的出发点是古典二阶模型。每台发电机用一个恒定内电势和暂态电抗表示转子方程只有两个状态变量功角δ和转速偏差。方程写出来就是dδ/dt ω0(ω_pu - 1)dω_pu/dt (Pm - Pe - D(ω_pu - 1)) / (2H)其中H是发电机惯性时间常数单位秒Pm是机械功率标幺值Pe是电磁功率标幺值D是阻尼系数。古典模型的核心假设是励磁系统能够维持内电势幅值恒定暂态过程中原动机功率基本不变。对大多数故障后第一摆的暂态稳定问题这个假设够用而且状态变量少仿真速度快便于大量扫描故障场景。什么时候不能用古典模型当故障点离机端很近机端电压跌到接近零励磁系统的强励动作会显著改变内电势这时候古典模型会偏乐观或偏悲观取决于系统条件。工程中常见的做法是先用古典模型做初步筛选把失稳风险最高的断面挑出来再换上更详细的模型重新验证。详细模型会在发电机转子上增加励磁绕组、阻尼绕组等状态同时把励磁系统、调速器、PSS 都写进微分方程。多机系统时域仿真中这些模型并不是越多越好而是要和你的数据条件匹配——没有调速器参数硬加调速器只会让结果更虚。选择模型时可以参考这张表模型层次状态变量适用场景数据要求古典二阶模型δ、ω第一摆暂态稳定、CCT 扫描、故障筛选H、Xd、E详细机电暂态模型δ、ω、Eq、Ed、励磁/调速/PSS动态行为验证、阻尼分析、控制措施校验励磁参数、调速器参数、PSS 参数电磁暂态模型三相瞬时值、含网络暂态换流器暂态、过电压、次同步振荡线路分布参数、控制框图仿真步长微秒级多机系统短路故障后时域仿真通常落在第二层区间。第一层用于快速判断第三层用于特殊问题。你如果刚接手一个算例数据里只有惯性和电抗那就先按古典模型跑通流程等拿到详细参数后再升级不要一上来就追求状态变量的齐全。2.2 网络方程与增广导纳矩阵把故障“装进”矩阵里发电机模型定了接下来要回答的问题是电磁功率 Pe 怎么算答案是通过网络方程。时域仿真中电网本身被认为是瞬时平衡的不需要微分方程只需用代数方程描述。多机系统的网络节点电压和注入电流关系用节点导纳矩阵表示为Y * V I对于经典模型发电机内电势节点通过暂态电抗 Xd 连接到机端母线。因此常规做法是把发电机内电势节点也放进导纳矩阵形成增广节点导纳矩阵维度是“发电机数 母线数”。构造规则是发电机内节点 k 与对应机端母线 j 之间加入一条支路导纳为 y 1/(jXd_k)内节点和机端母线的自导纳分别累加 y互导纳置为 -y母线部分的原始导纳来自线路、变压器和负荷等效阻抗。负荷在暂态稳定仿真中常按恒定阻抗处理即在负荷母线处并联一个导纳由潮流断面中的电压和有功无功反推。这是第一版仿真最容易忽略的点负荷模型选得不对后面看什么都别扭。故障怎么进矩阵三相短路相当于在故障母线处接入一条接地支路导纳为一个很大的数线路跳闸相当于把该线路的支路导纳从 Ybus 中删除。用两个不同矩阵的比较来表示situation 1 用短路支路处理给故障母线增加接地导纳 Yf 1/ZfZf 取一个很小的过渡电阻值比如 0.0001 pu 到 0.001 pu。situation 2 用切除处理将故障线路两端母线间的导纳置零。实际仿真程序中会用三段网络矩阵故障前矩阵、故障中矩阵、故障后矩阵。在时间轴上切换这三个矩阵就搭建了完整的仿真场景。故障前运行点由潮流计算结果给出再通过潮流计算反推每台发电机的内电势幅值和初始功角保证仿真初始时刻系统处于平衡状态。这一步如果直接拿任意初值开始跑前面几几十毫秒会出现严重的功率震荡那不是物理现象是初始条件没对齐。3. 故障与切除工况怎么设短路类型、故障位置与 CCT 扫描设计模型搭好之后仿真场景的设计往往决定了结论的工程价值。很多新手拿到系统就随机挑一条线路短路跑一条曲线出来看着稳定就说“没问题”。这对于多机系统暂态稳定分析来说基本等于白做。故障位置、故障类型、切除时间三个要素必须成套定义。3.1 故障位置与类型机端短路、近端短路和远端短路的差别短路故障的严重程度跟离发电机的电气距离紧密相关。机端三相短路时该发电机端电压接近于零电磁功率几乎掉到零转子在机械功率作用下猛烈加速远端线路短路时系统中各台发电机的电磁功率跌落幅度小功角摆动也小。做稳定校核时最忌讳的是“平均用力”把几百条线路逐一跑一遍。节省时间的做法是先在拓扑图上找出送端电源集中、负荷中心薄弱的断面重点关注这些断面上的机端母线和送出线路首端。故障类型方面工程校核通常分两层第一层用三相短路评估最严重工况第二层用单相接地故障评估概率最高的事故场景。三相短路属于对称故障可以直接用正序增广导纳矩阵处理单相接地属于不对称故障严格来说要用序网法求故障支路导纳再折算到正序网络。很多商用稳定程序内部已经做了这一层转换界面里只要求你选“三相短路”“单相接地”“两相短路接地”即可。但自编代码时要注意不要把一个不对称故障简单套进正序 Ybus 的对称短路处理里那样会低估或高估故障严重程度。故障位置还有一点值得强调线路中间短路比母线短路更考验系统吗不一定。线路中间短路时故障分量对两侧系统的影响相对均衡母线短路则把全部短路电流集中在一个节点上对相邻机组冲击更集中。实际操作中两者都要扫重点看机端母线、变压器高压侧母线和长线路首端。3.2 故障时序怎么定故障时刻、切除时刻和仿真总时长一条完整的时域仿真时序包括三个时刻故障发生时刻 t0、故障切除时刻 tc、仿真终止时刻 T_end。故障发生前系统处于稳态潮流断面就是初始状态t0 到 tc 是故障持续期tc 之后是故障后网络系统开始进入转子摇摆过程。如果考虑重合闸还要把重合时刻加入变成一个更长的时序链。切除时刻怎么取稳定校核中不取保护实际动作时间而是扫“极限切除时间”也就是最大能保持暂态稳定的故障持续时长。打个比方某线路首端三相短路的 CCT 是 0.28 s而该线路主保护动作时间是 0.1 s那么裕度足够系统风险低如果 CCT 只有 0.11 s保护动作稍慢就失稳这时候就要考虑切机、快关汽门等控制措施。仿真总时长也不是随便定的。古典模型下观察首摆稳定仿真 5 s 基本够带励磁和调速器、存在弱阻尼模式的系统区域间振荡周期可能到 3~8 s最好仿真到 10~20 s否则你会看到系统在一两秒内“回稳”却在五六秒后慢慢摆开。一个常见的工程序列是0 s 开始1 s 加故障用 1 s 作为故障前稳定运行段避开初始时刻的数值过渡故障持续 tc-1 stc 时刻切除仿真到 15 s。代码里把 t0 放在 1 s而不是 0 s目的是让数值积分器先跑一段“纯稳态”预热这在自编代码时能少很多麻烦。3.3 用二分法扫出极限切除时间 CCTCCT 不是数值积分直接给出的它是一组时序仿真结果的边界。同一故障位置切除时间短则系统稳定切除时间拉长到某个临界值之后系统失稳这个临界值就是 CCT。多机系统没有解析表达式只能用数值搜索。最稳定的搜索方法是二分法设定下限 t_low 0.05 s已知大概率稳定设定上限 t_high 0.6 s已知大概率失稳在仿真实测中调整取 tc (t_low t_high) / 2跑一次仿真如果稳定把 t_low 更新为 tc如果失稳把 t_high 更新为 tc重复 12~15 次把边界收敛到你想要的精度。这个搜索每轮都要跑一次完整时域仿真所以计算量不小。对于多机大系统可以先在古典模型下用大步长快速扫锁定关键故障再对重要故障用详细模型做精扫。二分法对“稳定/失稳判据”的依赖很强判据不清晰时二分结果会来回跳这一块留到第五章详细说。4. 跑通一个三机多母线算例核心代码与关键参数理论讲清楚之后最直接的上手方式是搭一个三机系统跑通流程。这里我以一台“三机、三条母线的自定义测试系统”为例重点展示核心解算逻辑如果你手头有 IEEE 3 机 9 节点系统或 Kundur 四机两区系统的数据只需要替换 Ybus 和发电机参数数组代码无需改动。4.1 算例的数据组织母线导纳矩阵和发电机参数仿真程序里数据通常按三类组织母线数据、支路数据、发电机数据。母线数据给出编号和负荷支路数据给出线路/变压器的阻抗和导纳发电机数据给出机端母线、惯性常数、暂态电抗和机械功率。组织好之后第一步是用潮流计算结果构建故障前的母线导纳矩阵 Ybus_pre。这一步可以直接从潮流程序导出也可以自编潮流计算得到。用 Python 做这个事最清晰的结构是把 Ybus 和发电机参数封装成两个对象import numpy as np from scipy.integrate import solve_ivp # 发电机参数: H(s), Xd(pu), Pm(pu), 机端母线编号(0-based) # 这里用一组示意值, 实际应替换为潮流/铭牌数据 gen_data [ {bus: 0, H: 5.0, Xd: 0.25, Pm: 0.6}, {bus: 1, H: 4.0, Xd: 0.30, Pm: 0.4}, {bus: 2, H: 4.5, Xd: 0.28, Pm: 0.5}, ] # Ybus_pre: 故障前母线导纳矩阵, 来自潮流数据 # Ybus_fault: 故障中矩阵(在故障母线增加接地导纳后) # Ybus_post: 故障后矩阵(切除线路后) Ybus_pre np.array([[5-20j, -28j, -14j], [-28j, 4-15j, -39j], [-14j, -39j, 4-14j]]) Ybus_fault Ybus_pre.copy() Ybus_fault[0, 0] 1 / (1e-4 1j*1e-4) # 母线0三相短路接地 Ybus_post Ybus_pre.copy() Ybus_post[0, 1] 0 Ybus_post[1, 0] 0 # 切除母线0-1之间的线路这段代码中gen_data 每台机组都带上了自己所属的机端母线编号。注意 H、Xd 和 Pm 都是用系统基准容量归一的标幺值或秒值做多机系统时如果各发电机额定容量不同必须先归一到统一基准否则后面每台机的加速度会失真。Ybus_fault 里给故障母线加接地导纳这是三相短路近似处理Zf 取 1e-4 pu 级别的过渡电阻是为了避免导纳矩阵奇异也是商用程序里常见做法。4.2 核心代码增广导纳矩阵与转子运动方程有了网络矩阵和发电机参数下一步是把发电机内节点扩展进去得到增广导纳矩阵。然后在每一个积分步中求解网络方程计算出电磁功率 Pe再喂给转子运动方程n_gen len(gen_data) n_bus Ybus_pre.shape[0] N n_gen n_bus def build_extended(Ybus): 把发电机内节点扩展到母线导纳矩阵中 Y np.zeros((N, N), dtypecomplex) Y[n_gen:, n_gen:] Ybus.copy() for k, g in enumerate(gen_data): y 1 / (1j * g[Xd]) # 内节点到机端母线的导纳 b n_gen g[bus] # 机端母线在增广矩阵中的位置 Y[k, k] y Y[b, b] y Y[k, b] - y Y[b, k] - y return Y Y_aug_pre build_extended(Ybus_pre) Y_aug_fault build_extended(Ybus_fault) Y_aug_post build_extended(Ybus_post) omega_s 2 * np.pi * 50.0 # 工频角速度, rad/s H np.array([g[H] for g in gen_data]) Pm np.array([g[Pm] for g in gen_data]) D np.zeros(n_gen) # 先不加阻尼, 便于观察振荡 # 初始功角由潮流结果提供, 这里以0.2/0.35/0.5 rad示意 delta0 np.array([0.2, 0.35, 0.50]) omega0 np.ones(n_gen) # 标幺转速初值1 def compute_pe(delta, Y_aug): 由内电势求电磁功率 E np.exp(1j * delta) # 假设内电势幅值标幺值1.0 Iinj E / (1j * np.array([g[Xd] for g in gen_data])) V np.linalg.solve(Y_aug, np.concatenate([Iinj, np.zeros(n_bus)])) Pe np.real(V[:n_gen] * np.conj(Iinj)) return Pe def rhs(t, x, Y_aug): delta, omega x[:n_gen], x[n_gen:] Pe compute_pe(delta, Y_aug) d_delta omega_s * (omega - 1.0) d_omega (Pm - Pe - D * (omega - 1.0)) / (2 * H) return np.concatenate([d_delta, d_omega])这里最核心的部分是 compute_pe。它把每台发电机内节点看作电流源求解包含发电机内节点的增广网络方程得到内节点电压再计算出每台机的电磁功率。这样做的好处是故障前、故障中、故障后只需要切换 Y_aug不用改变微分方程结构。机械功率 Pm 是常数这是古典模型的设定如果你之后加入调速器Pm 就要变成状态变量此处留好了扩展位。4.3 按时间段积分故障、切除两次切换不能少仿真过程需要分三段推进稳态段、故障段、故障后段。我用 solve_ivp 分别积分并把上一段的终点状态作为下一段初值def simulate(t_fault, t_clear, t_end, x0): 经典模型三段式时域仿真 t_seg1 (0, t_fault) t_seg2 (t_fault, t_clear) t_seg3 (t_clear, t_end) sol1 solve_ivp(rhs, t_seg1, x0, args(Y_aug_pre,), methodRK45, rtol1e-6, atol1e-8, max_step0.02) sol2 solve_ivp(rhs, t_seg2, sol1.y[:, -1], args(Y_aug_fault,), methodRK45, rtol1e-6, atol1e-8, max_step0.005) sol3 solve_ivp(rhs, t_seg3, sol2.y[:, -1], args(Y_aug_post,), methodRK45, rtol1e-6, atol1e-8, max_step0.02) ts np.concatenate([sol1.t, sol2.t, sol3.t]) xs np.concatenate([sol1.y, sol2.y, sol3.y], axis1) return ts, xs ts, xs simulate(t_fault1.0, t_clear1.2, t_end10.0, x0np.concatenate([delta0, omega0]))注意这里 max_step 的设置故障段步长收紧到 0.005 s稳态和故障后段放宽到 0.02 s。原因是短路瞬间电磁功率变化剧烈显式 RK45 如果步长太大会在故障切除时刻附近产生明显误差。t_fault 放在 1.0 s刚建好 1 s 稳态输出曲线里也容易区分时间段。切出故障后的结果取出相对功角曲线就能判断稳定与否相对功角摆开一定角度后回落并持续衰减说明系统稳定如果某个相对功角单调增大且不再回头说明失稳。5. 时域仿真常见的 5 个坑现象、原因与处理自编时域仿真程序的人几乎都要经历一样的过程模型看着没错参数没什么异常结果却完全不像物理。下面这五条是从实际操作中整理出来的高频问题每一条都按“现象 → 原因 → 解决”写直接对照排查。仿真一开始就出 NaN或者前 0.1 s 功角剧烈震荡。现象是电压或功角瞬间出现非数值或者功角在 0.05 s 内来回甩了几十度。原因大概率是初始条件没对齐发电机内电势幅值、初始功角与潮流结果不一致或者增广导纳矩阵在构造时有符号错误。解决方法是先把故障前状态跑一遍纯稳态仿真功角、电压应该是一条水平线如果不是就先别做故障场景回头查 Ybus 和 E 的计算。故障切除时刻附近结果对步长极其敏感换个 max_step 结论就反转。这是显式积分器处理刚性系统的经典翻车。短路的过渡电阻很小增广导纳矩阵条件数很大方程偏刚性如果在切除时刻附近步长太粗电磁功率跌落和恢复的细节被抹掉CCT 可能差出几十毫秒。解决方法是故障段用 0.005 s 之外还要保证切换时刻本身落在积分采样点上solve_ivp 分段调用能够保证这一点。自写 RK4 循环时要专门做“跨切换点步长截断”不要让积分器一步跨过切除时刻。功角曲线第一摆看着稳了仿真拉到 10 s 却失稳。这种问题在弱阻尼多机系统里很典型。前一两秒相对功角可能只是小幅摆动但区域间振荡模式周期较长负阻尼或零阻尼时振荡会缓慢增长直到失去同步。这不是仿真 bug是仿真时长不够。古典模型至少仿真到 5~8 s加详细励磁和调速器后建议 10~20 s。如果只是为了找 CCT失稳通常发生在首摆可以缩短时长但做控制措施评估必须拉长观察窗口。二分法扫 CCT 不收敛相同切除时间有时判稳定、有时判失稳。很多时候不是数值随机性而是稳定判据写得模糊。比如只看“最大相对功角是否超过 180°”而 180° 并不是严格判据系统可能在 170° 以后回摆也可能在 200° 时触发保护切机而回到稳定。解决方法是定义一个统一判据仿真结束时任意两台发电机的相对功角差小于设定阈值且角速度差趋近于零才判稳相对功角差单调增加且超过工程阈值判失稳。再用角速度符号变化辅助判断“首摆是否回头”二分结果就会稳定很多。负荷模型用恒定功率结果出现莫名其妙的“低压失稳”。恒定功率负荷在母线电压跌落时仍保持功率不变相当于等效导纳增大故障期间加重了发电机负担会让结果偏保守。实际配电网中负荷含大量恒阻抗分量电压低时负荷功率会自然下降。解决方法是先用恒定阻抗负荷模型跑基准工况观察摇摆形态是否正常再换成 ZIP 模型做敏感性分析。多机系统稳定结论如果对负荷模型特别敏感说明系统稳定裕度本身不高这一点恰恰是分析中要重点写出来的结论。6. 把仿真再往前推一步自动判稳与二分搜索CCT 扫描是最值得自动化的一步。人工一个个改切除时间跑仿真不仅慢而且容易误操作。我一般把仿真函数封装成一个返回布尔值的接口参数是故障位置和切除时间内部完成三段积分和稳定判据判断然后直接套二分循环。def is_steady(tc, t_fault1.0, t_end10.0, max_rel_angle2.5): ts, xs simulate(t_faultt_fault, t_cleartc, t_endt_end, x0np.concatenate([delta0, omega0])) n n_gen delta xs[:n, :] # 各发电机功角逐时刻曲线 rel np.max(delta, axis0) - np.min(delta, axis0) return np.all(rel max_rel_angle) # 全程最大相对功角差阈值判断 lo, hi 0.05, 0.50 for _ in range(15): tc 0.5 * (lo hi) if is_steady(tc): lo tc else: hi tc print(ftc{tc:.4f}s - {稳定 if is_steady(tc) else 失稳})阈值 2.5 rad约 143°只是示意值工程上需要根据系统同步稳定标准调整。更严谨的做法是同时检查最后一小段时间内相对功角差的斜率如果斜率持续为正说明还在拉开不能判稳。二分循环跑通之后下一步就是并行扫描。多机系统的故障场景非常多每个场景一次时域仿真耗时从几秒到几分钟不等按故障位置和故障类型做两层循环然后用多进程并行能省下大量等待时间。我自己的习惯是先做一波粗略扫描用 2 s 仿真时长、粗阈值筛掉明显稳定的故障对落在边界附近的故障再用 10 s 时长、严格阈值精扫。两轮筛选下来计算量通常能少一半以上。验证仿真结果的方法也很重要。如果有现场故障录波或 PMU 数据把录波的功角差、电压曲线和仿真曲线叠在一起看趋势这比任何“理论正确”都有说服力。模型参数不准确的时候曲线首摆可能吻合后续摆动却逐渐漂移这说明阻尼参数或励磁模型还不够准需要继续校模型。这一步做得好时域仿真就不只是“算一个数”而是能真正支撑定值整定和稳定控制决策。我这些年做多机系统时域仿真吃过最大的亏就是忽略稳定判据的一致性同一套数据在两个人手里判出相反结论。后来把所有仿真统一封装成“输入故障、输出是否稳定”的标准函数才把重复劳动变成批量扫描。希望这些踩坑经验能帮你少跑几次数小时的无效仿真。本文还有配套的精品资源点击获取