简介面向电力系统暂态稳定研究的一份MATLAB仿真资源聚焦三机系统线路AB段首端两相短路接地故障后的时域动态过程。资源针对多机系统故障分析需求给出了从故障设定到0.1秒后切除故障线路的完整仿真流程适合电力系统方向学生、研究生及入门工程师用于暂态稳定机理学习与仿真验证。压缩包内共5个m文件整体仅4KBmain.m负责仿真控制与参数设置calculate.m完成故障前后系统模型与核心计算Improved_euler.m采用改进欧拉法迭代求解微分方程data1.m存储发电机、线路及故障参数readcase.m用于读取案例数据。脚本结构清晰可快速了解多机暂态仿真的程序组织方式。目前已有395人学习下载。通过运行与研读代码可掌握短路故障建模、数值积分步进、故障切除时间设定等关键环节直观观察发电机转速、电压电流等动态响应并据此分析系统能否恢复稳定为保护策略优化提供参考。1. 多机系统短路故障后时域仿真从一次 10ms 的切除时间偏差说起同样是母线三相短路故障切除时间差 10ms多机系统的功角曲线可能从收敛变成发散。多机系统短路故障后时域仿真就是把电网故障前、故障中、故障后三个阶段统一写成微分-代数方程组从 0 时刻开始一步步积分出每台发电机的功角、转速响应回答一个问题这种故障、这种切除方式、这个运行方式下系统到底稳不稳哪台机先失稳。做继电保护定值校核、调度运行方式分析、新能源并网稳定评估的人几乎每天都离不开这套流程。这篇文章按建模、故障注入、参数设定、避坑、指标输出一条线走完新手能跟着跑通熟手可以对照检查自己在判稳和故障阻抗处理上的习惯是否有问题。2. 暂态稳定为什么必须用时域仿真模型、判据与适用边界2.1 多机系统里“暂态稳定”到底看什么单机无穷大系统里功角只有一个稳定判据也直接功角越过不稳定平衡点就失步。多机系统完全不是这个逻辑。三台以上发电机彼此通过线路和负荷耦合短路故障时每台机的加速功率不一样故障切除后它们之间会出现相对摆动。失稳模式可能是单机相对其余机组甩开也可能是两群机组对甩极端情况下还有连锁失步。只看某一台机的绝对功角没有意义因为参考机自身也在摆。实际工程中判断多机系统暂态稳定看的是任意两台机之间的最大相对功角差或者各机组相对惯性中心COI的偏移量。用 COI 坐标是比较稳的做法。先算惯性中心的加权角δ_coi (Σ H_i * δ_i) / (Σ H_i)然后把每台机换成相对 COI 的角差 Δδ_i δ_i - δ_coi。失稳的核心指标是任意 Δδ_i 和其余机组的偏差是否持续扩大并越过 180°。需要注意的是多机系统首摆稳定不代表后面几摆稳定。有的工况第一摆最大相对角差 120°看着安全第二摆直接冲到 200° 失步这就是多机多摆失稳只截取前两秒曲线根本看不出来。2.2 直接用等面积法则行不行时域仿真不可替代的原因等面积法则在单机无穷大系统里非常漂亮加速面积小于减速面积就稳定。这套理论用到多机上工程上最常见的处理是“同调等值”——先把多台机聚合成两群再对两群系统用扩展等面积法则判断。PSASP 里的 EEAC、DEEAC 模块就是这个思路。问题是分群结果对故障点位置和切除时间极敏感同一个故障G1 和 G3 在第一摆同调第二摆可能就不同调了分错群以后等值参数全错结论自然翻车。时域仿真完全不依赖分群假设。它把每台发电机的摇摆方程和网络方程一起积分功角曲线是“积”出来的不是“套”出来的。哪怕故障场景再复杂——不对称短路、断路器单相拒动、重合闸时序、新能源变流器低压穿越只要数学模型和参数给得对仿真结果就是基准。调度部门做安稳校核、保护部门做断路器失灵校验最终认的都是时域仿真曲线而不是某个解析判据。所以时域仿真在暂态稳定分析里的地位就是“黑匣子但这个黑匣子最可信”。2.3 发电机模型选型经典二阶模型和详细模型的取舍多机时域仿真里发电机模型按详细程度分三档模型描述适用场景需要参数结论偏差经典二阶恒定内电势 E暂态电抗 xd 后电压恒定快速扫描、方案初筛、故障点扫描H、xd、E偏保守三相短路近区会夸大失稳计及励磁的详细模型E 随励磁调节变化带励磁机模型最终校核、投运前分析励磁系统参数、调压器增益更接近实际结论偏乐观计及励磁调速器PSS功频调节和阻尼控制全部参与实际电网校核、多摆稳定分析PSS 整定值、调速器参数最贴近现场参数错误也会带来新偏差我一般的习惯是研究阶段用经典模型扫故障集把临界切除时间摸出大致范围再对边界工况换详细模型复核一遍。经典模型因为忽略励磁强行增磁效应普遍低估系统稳定能力所以它判失稳的工况详细模型下大概率还是失稳它判稳定的工况才需要进一步确认。反过来直接把详细模型套到每个故障点上光收敛和调参数就够折腾一晚上。2.4 大扰动和小扰动仿真为什么线性化方法在这里失效有人会问静态稳定分析和特征值分析也能判断系统稳不稳何必积分几十秒。区分在扰动大小。特征值分析处理的是运行点附近的小扰动把方程在工作点线性化看主导特征值实部是否小于零。三相短路是拓扑级突变功角被推到远离平衡点的位置线路电抗、发电机内电势、负荷电压特性全部偏离线性化假设线性化方程在这个尺度下失真严重。时域仿真没有这个限制。故障就是修改导纳矩阵切除就是再改一次导纳矩阵方程结构和参数突变都能直接表达唯一要求是数值积分方法足够稳、步长足够小。这也是为什么时域仿真是整个电力系统稳定分析体系里的“裁判”能量函数法、特征值法、等面积法则的结果边界最终都要拿时域仿真来校验。3. 多机系统时域仿真落地三机系统建模与短路故障注入全过程3.1 从数据到可仿真系统的五个准备步骤拿到一个多机系统不管是有 3 台机还是 30 台机时域仿真建模流程都是固定的五步。第一步整理发电机参数至少要有惯性时间常数 H、暂态电抗 xd、次暂态电抗 xd经典模型只需要前两个。第二步整理网络参数线路电阻、电抗、对地充电电容变压器漏抗和变比。第三步确定负荷模型恒阻抗、恒电流还是恒功率对暂态稳定结果影响最大的就是负荷电压特性。第四步做潮流计算得到各母线电压幅值相角和各发电机出力这是时域仿真的初值初值不对后面全是白算。第五步把发电机动态方程和网络方程拼接形成微分-代数方程组准备积分。日常工程里初值这步最常见的做法是用 BPA 或 PSASP 先算一个潮流导出母线电压、发电机出力再填进时域仿真里。自己写 Python 做研究时也可以用 pypower 的 runpf 接口直接取初值。下面给出的示例是一个自洽的三机三机端系统代码里用初始功角对应的电磁功率作为机械功率初值保证仿真从 t0 开始就是平衡的不需要额外潮流软件介入。3.2 三机系统参数与导纳矩阵构造为了把核心的时域仿真流程讲清楚这里用一个三机系统示例。三台机经机端母线构成三角形环网线路参数和发电机参数如下表。单位全部为标幺值基准容量取 100MVA。对象参数数值线路 L12电抗j0.20线路 L23电抗j0.25线路 L31电抗j0.30G1H / xd / E6.0s / 0.20 / 1.05G2H / xd / E4.5s / 0.25 / 1.03G3H / xd / E5.0s / 0.22 / 1.04初始功角δ1 / δ2 / δ30° / 8° / 15°构建含发电机内节点的增广导纳矩阵前 3 个节点是机端母线后 3 个节点是发电机内电势节点每一台发电机通过暂态电抗 xd 接到对应机端母线。线路参数已知直接组母线导纳矩阵再用 xd 把内节点耦合进来。import numpy as np # 三机系统数据标幺值 X12, X23, X31 0.20, 0.25, 0.30 # 线路电抗 y12, y23, y31 1/(1j*X12), 1/(1j*X23), 1/(1j*X31) H np.array([6.0, 4.5, 5.0]) # 惯性时间常数秒 Xd np.array([0.20, 0.25, 0.22]) # 暂态电抗 E np.array([1.05, 1.03, 1.04]) # 内电势幅值 delta0 np.deg2rad([0.0, 8.0, 15.0]) # 初始功角弧度 def build_base_ybus(): 机端母线导纳矩阵三角形环网忽略电阻和对地电容 Y np.zeros((3, 3), dtypecomplex) for i, j, y_line in [(0, 1, y12), (1, 2, y23), (2, 0, y31)]: Y[i, i] y_line Y[j, j] y_line Y[i, j] - y_line Y[j, i] - y_line return Y def build_full_ybus(Ym): 把发电机内节点纳入形成6节点增广导纳矩阵 n len(Ym) Y np.zeros((2*n, 2*n), dtypecomplex) Y[:n, :n] Ym for g in range(n): y_xd 1 / (1j * Xd[g]) Y[g, g] y_xd # 机端节点贡献 Y[ng, ng] y_xd # 内节点自导纳 Y[g, ng] -y_xd Y[ng, g] -y_xd return Y这段代码里build_base_ybus 组装的是线路网络的母线导纳矩阵环网结构每一条线路在两个节点上各自产生自导纳和互导纳。build_full_ybus 再把发电机暂态电抗以串联支路形式接入发电机内节点作为独立的恒压源节点这个增广矩阵是后续网络求解的基础。注意 xd 支路只有内节点到机端节点一条路径内节点之间没有直接耦合它们的耦合完全通过机端母线网络实现。3.3 短路故障如何“注入”节点三相短路的导纳修改法时域仿真里短路故障不是“删掉节点”而是“修改导纳”。金属性三相短路时故障点对地阻抗趋近于零相当于在故障节点上并联一个极大导纳。数值上不能让导纳直接取无穷大常见做法是取一个很小的阻抗 Zf 1e-4 j1e-6然后向故障节点的自导纳上加上 1/Zf。这样故障节点的电压被压到接近零电磁功率降到接近零发电机进入加速状态。线路中点故障需要多一步处理。以线路 2-3 中点三相短路为例把线路拆成两段每段电抗取原线路的一半中间新增一个故障节点 4。正常运行时节点 4 是隐藏的故障时刻给节点 4 加接地导纳切除时刻不仅移除该接地导纳还要把整条线路 2-3 从网络中摘除对应 Y23 和 Y32 清零同时从 Y22 和 Y33 的自导纳中减去原线路导纳。这一步最容易出错只清互导纳、不清自导纳导纳矩阵不对称仿真结果会直接飞掉。def apply_fault(Ym, fault_node, Zf1e-4 1j*1e-6): 故障注入在指定机端母线并联小阻抗接地 Yf Ym.copy() Yf[fault_node, fault_node] 1 / Zf return Yf def remove_line(Ym, i, j, y_line): 切除线路互导纳清零自导纳相应修正 Yp Ym.copy() Yp[i, i] - y_line Yp[j, j] - y_line Yp[i, j] 0.0 Yp[j, i] 0.0 return Ypapply_fault 处理的是机端母线短路Zf 取极小值模拟金属性短路的极强接地效果。remove_line 对应保护动作切除故障线路注意自导纳修正和互导纳清零必须成对出现。如果仿真场景是母线故障且保护跳开母联那在 Ypost 里还要额外处理母线所带电源/负荷的转移这里不过度展开工程分析里由故障设置模块统一管理。3.4 电磁功率求解与摇摆方程积分主循环有了导纳矩阵每步积分先求机端母线电压再反推发电机电磁功率。内节点电压已知为 E∠δ机端母线电压通过解网络方程 Ymm·Vm -Ymg·Vg 得到然后算发电机注入网络的电流和电磁功率。def solve_network(Yfull, delta): 给定内电势相角求解机端电压和发电机电磁功率 n len(delta) Ymm Yfull[:n, :n] Ymg Yfull[:n, n:] Ygm Yfull[n:, :n] Ygg Yfull[n:, n:] Vg E * np.exp(1j * delta) # 内电势相量 Vm np.linalg.solve(Ymm, -Ymg Vg) # 机端母线电压 Ig Ygm Vm Ygg Vg # 流入网络的发电机电流 Pe np.real(Vg * np.conj(Ig)) # 电磁功率 return Vm, Pesolver 的核心是五步内电势相量由幅值和功角合成机端电压通过线性方程组求得发电机电流由增广导纳矩阵的耦合分块计算电磁功率取内电势电压乘电流共轭的实部。这个功率表达式已经包含了各机之间的相互影响不需要再单独列出功率角方程。摇摆方程的积分采用经典二阶模型每台机两个状态量功角和角速度。def build_rhs(Y0, Yf, Ypost, t_fault, t_clear, D0.0): 构造微分方程右侧函数按时间阶段切换导纳矩阵 n len(H) omega0 2 * np.pi * 50.0 def ybus_at(t): if t t_fault: return Y0 elif t t_clear: return Yf else: return Ypost def rhs(t, x): delta x[:n] omega x[n:] Y ybus_at(t) _, Pe solve_network(Y, delta) ddelta omega - omega0 domega (omega0 / (2 * H)) * (Pm - Pe - D * (omega - omega0)) return np.concatenate([ddelta, domega]) return rhs def rk4_step(rhs, t, x, dt): 经典四阶龙格库塔单步积分 k1 rhs(t, x) k2 rhs(t 0.5*dt, x 0.5*dt*k1) k3 rhs(t 0.5*dt, x 0.5*dt*k2) k4 rhs(t dt, x dt*k3) return x (dt/6.0) * (k1 2*k2 2*k3 k4)rhs 函数内部按 t 自动选择当前阶段的导纳矩阵故障发生在 t_fault切除发生在 t_clear三个阶段各自对应 Y0、Yf、Ypost。注意初始状态下把 Pm 初始化为故障前的电磁功率这是保证仿真开始不出现虚假功率冲击的关键。rk4_step 是四阶龙格库塔的标准写法每一步调用四次右侧函数步长建议取 0.001 到 0.005 秒。整个主循环模拟故障前的稳态段、故障段和故障后段完整记录每台机的功角轨迹最后按最大相对功角差判稳。Y0 build_full_ybus(build_base_ybus()) Yf build_full_ybus(apply_fault(build_base_ybus(), fault_node2)) Ypost Y0 # 母线故障清除线路保持完整 n len(H) Pm solve_network(Y0, delta0)[1] # 取初始电磁功率为机械功率 x0 np.concatenate([delta0, 2*np.pi*50.0*np.ones(n)]) t_fault, t_clear 0.10, 0.25 dt, t_end 0.002, 20.0 rhs build_rhs(Y0, Yf, Ypost, t_fault, t_clear, D0.5) t_list, angle_list [0.0], [delta0.copy()] x x0 t 0.0 while t t_end: x rk4_step(rhs, t, x, dt) t dt t_list.append(t) angle_list.append(x[:n].copy()) angle_array np.rad2deg(np.array(angle_list)) max_dev angle_array.max(axis1) - angle_array.min(axis1) stable max_dev.max() 180.0 print(失稳 if not stable else 稳定, 最大功角差, round(max_dev.max(), 1), 度)主循环从 0 到 20 秒步长 0.002 秒一共 1 万步三机系统跑起来非常快。故障时刻设在 0.1 秒0.25 秒切除模拟保护装置动作时间 100ms 加断路器开断时间 50ms 的典型组合。判稳逻辑取任意两台机最大功角差的最大值小于 180 度判稳超过则判失稳。3.5 换用标准节点系统从示例走向工程模型上面的三机系统用于理解原理够用真正工程分析一般不手写这个网络。日常工作是拿 BPA 或 PSASP 的卡片数据把网络和发电机参数直接读进来形成大导纳矩阵。标准算例里 IEEE 3 机 9 节点、IEEE 39 节点系统是校验新算法、新脚本最常见的两个网架。手写代码对接标准算例时关键在于把潮流初值吃准。BPA 潮流卡片里的节点电压、发电机出力直接映射成时域仿真初值不要自己算功角更不要凭空拍一组初值。pypower 用户可以用 loadcase 读 IEEE 9 节点数据然后 runpf 取潮流结果再把它转换成语义清晰的字典结构传给时域仿真函数网络求解和积分主循环完全复用上面这套框架。参数表一换代码主体一行不用动这就是分层建模的好处。4. 时域仿真关键参数怎么调切除时间、积分步长与故障阻抗的取值4.1 故障切除时间影响稳定结论最直接的参数短路故障时发电机输出电磁功率接近零机械功率不变转子持续加速。切除越晚加速时间越长功角冲得越高。同一个故障点0.2 秒切除时最大功角差 150 度0.25 秒切除就变成 210 度稳定结论直接翻转。切除时间的工程来源是保护动作时间加断路器固有开断时间一般 220kV 线路是 100ms 左右如果带断路器失灵保护二次跳闸要按 250ms 到 300ms 校核。做运行方式分析时通常不直接把切除时间写死而是扫出一个临界切除时间CCTCritical Clearing Time再对比保护实际动作时间。CCT 大于实际切除时间留有余量CCT 小于实际切除时间说明该故障下系统承受不住这个切除速度要调整保护定值或改变运行方式。扫 CCT 的标准做法是二分搜索下一章给完整脚本。4.2 故障阻抗与故障类型为什么不能直接填 0金属性短路的物理本质是故障点电压为零但数值仿真里给故障节点并联一个阻抗为零的支路导纳直接是无穷大线性方程组求解直接崩。实际金属性短路存在弧光电阻标幺值通常在 0.001 到 0.01 之间所以工程习惯取 Zf 1e-3 j1e-4 或更小既足够接近金属性短路又不至于让矩阵病态。故障类型处理方法备注三相短路节点并联小阻抗接地最严重工况稳定校核首选项两相相间短路复合序网或相域模型需负序网络单相接地复合序网中性点接地方式影响大发展性故障分时段切换故障类型需事件序列表三相短路对系统冲击最大也是暂态稳定标准校核的默认工况。两相和单相故障要构造负序、零序网络手写代码工作量陡增日常直接用商业软件的事件卡片配置原理上都是修改故障节点导纳只是对称分量法要把正负零序网络分开求解再合成。4.3 积分步长和仿真时长数值稳定性与多摆失稳的权衡显式 RK4 在刚性系统里对步长很挑剔。电力系统网络方程本身是代数方程发电机转子方程时间常数在秒级但励磁绕组和网络暂态会把系统变刚。只考虑机电暂态、用经典模型时0.005 秒步长通常没问题一旦引入励磁模型或双馈风机变流器模型建议直接降到 0.001 秒否则高频振荡一旦激发结果直接废掉。仿真时长要覆盖到多摆过程。首摆峰值一般出现在故障切除后 0.3 到 1.0 秒但第二摆、第三摆可能在 2 到 5 秒才达到峰值。阻尼较弱的系统每隔几秒摆一次摆幅缓慢增大20 秒内必然越过 180 度。我一般设 20 秒如果 20 秒内相对角差持续衰减且最大值低于 150 度判稳定的信心才足够。4.4 阻尼系数 D 的影响与整定经典模型里的阻尼系数 D 代表机械阻尼和电气阻尼的合计效果取值很玄学。汽轮机组一般取 1 到 2水轮机组取 0.5 到 1研究性算例经常直接取 0相当于最保守地忽略一切阻尼因素。D 取大了摇摆曲线衰减快容易掩盖系统真实的负阻尼问题取小了曲线长期不衰减又会把本来稳定的系统判成临界。我的习惯是先用 D0 扫 CCT得到保守下界再按实际机组阻尼经验值取 D1 复核几个关键故障两次结果差异大说明结论对阻尼敏感这种情况下要谨慎下结论。5. 多机时域仿真避坑指南五个让结果翻车的常见工程错误5.1 故障阻抗设成 0矩阵奇异与数值爆炸现象仿真跑到故障时刻瞬间电压、功率输出 NaN或者功角曲线出现不合理的跳变日志里能看到 LinAlgError 奇异矩阵报错。原因故障节点并联阻抗为 0对应导纳无穷大导纳矩阵在该节点出现溢出或严重病态。金属性短路在物理上确实是零阻抗但数值计算需要一个小阻抗来代表弧光电阻和数值容差。解决故障阻抗取 Zf 1e-3 j1e-4敏感问题可以再压到 1e-4 j1e-6。不要取 0也不要用 1e-8 这种极端小值否则矩阵条件数过差结果仍然不可信。做灵敏度分析时Zf 变化一个数量级对功角曲线的影响极小说明这个取值是合理的。5.2 初值不自洽仿真开始就出现虚假摇摆现象t0 到 t0.02 秒内功角就快速偏移然后在 0.1 秒故障前已经摆出十几度故障时刻的状态完全偏离工作点。原因机械功率 Pm 和初始电磁功率 Pe 不相等。时域仿真从静态工作点出发发电机转子初始加速度必须为零。很多手写脚本直接拍一个 Pm没考虑初始潮流对应的 Pe仿真开始时转子就带上了初速后面所有结论全部失真。解决初始化时先解一次网络方程把故障前状态的 Pe 算出来再令 Pm Pe。BPA 和 PSASP 里这个动作是自动完成的但在自写代码里是显式步骤。代码写法是Pm solve_network(Y0, delta0)[1]一行代码解决。5.3 切除线路时只改互导纳、忘改自导纳现象故障切除后功角曲线出现异常阶跃或者切除前后功率不平衡曲线明显不连续甚至出现功率凭空增加的情况。原因线路切除涉及两个修改线路两端节点间的互导纳清零同时两个节点的自导纳要扣除该线路贡献的那部分导纳。只清互导纳自导纳还保留着线路的并联支路网络拓扑自相矛盾。解决把线路切除封装成独立函数按“自导纳减、互导纳清”的顺序成对修改每次切线后检查导纳矩阵是否对称。代码里 remove_line 已经体现了这个顺序实际工程里切双回线中的一回、切变压器支路逻辑完全一致。5.4 把参考机当“绝对参考”忽略参考机自身摇摆现象同一套仿真数据用 G1 做参考显示相对功角差最大 160 度判稳换个参考机显示超过 200 度判失稳。两边都有道理但结论打架。原因多机系统没有绝对静止的参考机。任何一台机都在摆动以单机无穷大的习惯选参考机参考机自身的大幅摆动会把真实相对运动淹没。解决改用 COI 参考系计算功角偏差或者直接用任意两机功角差的最大值作为判稳指标。代码里max_dev angle_array.max(axis1) - angle_array.min(axis1)就是取每时刻最大和最小功角的差这个差值不依赖具体参考机结论客观稳定。5.5 仿真时长太短只看到首摆漏掉多摆失稳现象5 秒仿真显示最大功角差 140 度曲线没有越过 180 度结论写“稳定”。把仿真延长到 20 秒发现第 9 秒第二摆冲到 230 度实际是失稳的。原因部分系统故障切除后首摆看起来收敛但故障期间注入的不平衡能量没有完全耗散后续摆动逐渐放大。尤其是阻尼弱、多机振荡模式耦合的系统失稳往往发生在第二、第三摆。解决仿真时长至少取 15 到 20 秒判稳前同时看最大功角差随时间的变化趋势——是持续减小还是波动增大。末尾阶段的转速偏差也要看转速不归零、相对角差振荡不衰减都不能轻易判稳。这个规则对任何多机系统都适用也是最容易因“跑得快”而踩的坑。6. 从功角曲线到稳定裕度临界切除时间的二分搜索与仿真检验6.1 用二分法扫出临界切除时间故障切除时间是运行方式分析里最关心的稳定边界。二分搜索比线性扫描效率高得多先定一个稳定肯定成立的切除时间下界再定一个必失稳的上界反复折半直到区间宽度小于允许误差。def is_stable(t_clear_candidate): t_fault 0.10 Yf build_full_ybus(apply_fault(build_base_ybus(), fault_node2)) x0 np.concatenate([delta0, 2*np.pi*50.0*np.ones(n)]) rhs build_rhs(Y0, Yf, Y0, t_fault, t_clear_candidate, D0.5) x, t x0, 0.0 max_angle 0.0 while t 20.0: x rk4_step(rhs, t, x, 0.002) t 0.002 rel max(x[:n]) - min(x[:n]) max_angle max(max_angle, np.rad2deg(rel)) if max_angle 180.0: return False return True lo, hi 0.10, 0.60 # 下界稳定上界失稳 while hi - lo 0.005: mid (lo hi) / 2 if is_stable(mid): lo mid else: hi mid print(临界切除时间约, round((lo hi) / 2, 3), 秒)二分搜索里每次调用 is_stable 都要完整跑一遍 20 秒仿真三机系统耗时很小几十次循环完全可接受大型系统建议把步长放宽到 0.005 秒做初扫锁定区间后再精扫。出来 CCT 以后拿保护实际切除时间做对比比值大于 1.3 一般认为是可靠接近 1.0 就要警惕。6.2 仿真结果的三条自检路径第一步长收敛性。把步长从 0.002 换成 0.001 和 0.0005 各跑一遍三条功角曲线主峰差在 1 度以内说明积分可靠。第二能量校核。故障切除后的转子动能增量应近似等于加速功率对时间的积分偏差超过 5% 说明积分或功率计算有逻辑错误。第三商业软件对照。同一工况的曲线走向和 CCT 结果和 BPA、PSASP 或 PSS/E 对比功角峰值差在几度内属于正常差十几度要检查模型参数和故障事件时序。我的个人习惯是 CCT 二分搜索永远配一条详细模型复核曲线两边对上了才敢写进分析报告。这套流程踩过的坑不少但把自检步骤固定下来以后翻车的概率低了很多希望帮到你。本文还有配套的精品资源点击获取
简介面向电力系统暂态稳定研究的一份MATLAB仿真资源聚焦三机系统线路AB段首端两相短路接地故障后的时域动态过程。资源针对多机系统故障分析需求给出了从故障设定到0.1秒后切除故障线路的完整仿真流程适合电力系统方向学生、研究生及入门工程师用于暂态稳定机理学习与仿真验证。压缩包内共5个m文件整体仅4KBmain.m负责仿真控制与参数设置calculate.m完成故障前后系统模型与核心计算Improved_euler.m采用改进欧拉法迭代求解微分方程data1.m存储发电机、线路及故障参数readcase.m用于读取案例数据。脚本结构清晰可快速了解多机暂态仿真的程序组织方式。目前已有395人学习下载。通过运行与研读代码可掌握短路故障建模、数值积分步进、故障切除时间设定等关键环节直观观察发电机转速、电压电流等动态响应并据此分析系统能否恢复稳定为保护策略优化提供参考。1. 多机系统短路故障后时域仿真从一次 10ms 的切除时间偏差说起同样是母线三相短路故障切除时间差 10ms多机系统的功角曲线可能从收敛变成发散。多机系统短路故障后时域仿真就是把电网故障前、故障中、故障后三个阶段统一写成微分-代数方程组从 0 时刻开始一步步积分出每台发电机的功角、转速响应回答一个问题这种故障、这种切除方式、这个运行方式下系统到底稳不稳哪台机先失稳。做继电保护定值校核、调度运行方式分析、新能源并网稳定评估的人几乎每天都离不开这套流程。这篇文章按建模、故障注入、参数设定、避坑、指标输出一条线走完新手能跟着跑通熟手可以对照检查自己在判稳和故障阻抗处理上的习惯是否有问题。2. 暂态稳定为什么必须用时域仿真模型、判据与适用边界2.1 多机系统里“暂态稳定”到底看什么单机无穷大系统里功角只有一个稳定判据也直接功角越过不稳定平衡点就失步。多机系统完全不是这个逻辑。三台以上发电机彼此通过线路和负荷耦合短路故障时每台机的加速功率不一样故障切除后它们之间会出现相对摆动。失稳模式可能是单机相对其余机组甩开也可能是两群机组对甩极端情况下还有连锁失步。只看某一台机的绝对功角没有意义因为参考机自身也在摆。实际工程中判断多机系统暂态稳定看的是任意两台机之间的最大相对功角差或者各机组相对惯性中心COI的偏移量。用 COI 坐标是比较稳的做法。先算惯性中心的加权角δ_coi (Σ H_i * δ_i) / (Σ H_i)然后把每台机换成相对 COI 的角差 Δδ_i δ_i - δ_coi。失稳的核心指标是任意 Δδ_i 和其余机组的偏差是否持续扩大并越过 180°。需要注意的是多机系统首摆稳定不代表后面几摆稳定。有的工况第一摆最大相对角差 120°看着安全第二摆直接冲到 200° 失步这就是多机多摆失稳只截取前两秒曲线根本看不出来。2.2 直接用等面积法则行不行时域仿真不可替代的原因等面积法则在单机无穷大系统里非常漂亮加速面积小于减速面积就稳定。这套理论用到多机上工程上最常见的处理是“同调等值”——先把多台机聚合成两群再对两群系统用扩展等面积法则判断。PSASP 里的 EEAC、DEEAC 模块就是这个思路。问题是分群结果对故障点位置和切除时间极敏感同一个故障G1 和 G3 在第一摆同调第二摆可能就不同调了分错群以后等值参数全错结论自然翻车。时域仿真完全不依赖分群假设。它把每台发电机的摇摆方程和网络方程一起积分功角曲线是“积”出来的不是“套”出来的。哪怕故障场景再复杂——不对称短路、断路器单相拒动、重合闸时序、新能源变流器低压穿越只要数学模型和参数给得对仿真结果就是基准。调度部门做安稳校核、保护部门做断路器失灵校验最终认的都是时域仿真曲线而不是某个解析判据。所以时域仿真在暂态稳定分析里的地位就是“黑匣子但这个黑匣子最可信”。2.3 发电机模型选型经典二阶模型和详细模型的取舍多机时域仿真里发电机模型按详细程度分三档模型描述适用场景需要参数结论偏差经典二阶恒定内电势 E暂态电抗 xd 后电压恒定快速扫描、方案初筛、故障点扫描H、xd、E偏保守三相短路近区会夸大失稳计及励磁的详细模型E 随励磁调节变化带励磁机模型最终校核、投运前分析励磁系统参数、调压器增益更接近实际结论偏乐观计及励磁调速器PSS功频调节和阻尼控制全部参与实际电网校核、多摆稳定分析PSS 整定值、调速器参数最贴近现场参数错误也会带来新偏差我一般的习惯是研究阶段用经典模型扫故障集把临界切除时间摸出大致范围再对边界工况换详细模型复核一遍。经典模型因为忽略励磁强行增磁效应普遍低估系统稳定能力所以它判失稳的工况详细模型下大概率还是失稳它判稳定的工况才需要进一步确认。反过来直接把详细模型套到每个故障点上光收敛和调参数就够折腾一晚上。2.4 大扰动和小扰动仿真为什么线性化方法在这里失效有人会问静态稳定分析和特征值分析也能判断系统稳不稳何必积分几十秒。区分在扰动大小。特征值分析处理的是运行点附近的小扰动把方程在工作点线性化看主导特征值实部是否小于零。三相短路是拓扑级突变功角被推到远离平衡点的位置线路电抗、发电机内电势、负荷电压特性全部偏离线性化假设线性化方程在这个尺度下失真严重。时域仿真没有这个限制。故障就是修改导纳矩阵切除就是再改一次导纳矩阵方程结构和参数突变都能直接表达唯一要求是数值积分方法足够稳、步长足够小。这也是为什么时域仿真是整个电力系统稳定分析体系里的“裁判”能量函数法、特征值法、等面积法则的结果边界最终都要拿时域仿真来校验。3. 多机系统时域仿真落地三机系统建模与短路故障注入全过程3.1 从数据到可仿真系统的五个准备步骤拿到一个多机系统不管是有 3 台机还是 30 台机时域仿真建模流程都是固定的五步。第一步整理发电机参数至少要有惯性时间常数 H、暂态电抗 xd、次暂态电抗 xd经典模型只需要前两个。第二步整理网络参数线路电阻、电抗、对地充电电容变压器漏抗和变比。第三步确定负荷模型恒阻抗、恒电流还是恒功率对暂态稳定结果影响最大的就是负荷电压特性。第四步做潮流计算得到各母线电压幅值相角和各发电机出力这是时域仿真的初值初值不对后面全是白算。第五步把发电机动态方程和网络方程拼接形成微分-代数方程组准备积分。日常工程里初值这步最常见的做法是用 BPA 或 PSASP 先算一个潮流导出母线电压、发电机出力再填进时域仿真里。自己写 Python 做研究时也可以用 pypower 的 runpf 接口直接取初值。下面给出的示例是一个自洽的三机三机端系统代码里用初始功角对应的电磁功率作为机械功率初值保证仿真从 t0 开始就是平衡的不需要额外潮流软件介入。3.2 三机系统参数与导纳矩阵构造为了把核心的时域仿真流程讲清楚这里用一个三机系统示例。三台机经机端母线构成三角形环网线路参数和发电机参数如下表。单位全部为标幺值基准容量取 100MVA。对象参数数值线路 L12电抗j0.20线路 L23电抗j0.25线路 L31电抗j0.30G1H / xd / E6.0s / 0.20 / 1.05G2H / xd / E4.5s / 0.25 / 1.03G3H / xd / E5.0s / 0.22 / 1.04初始功角δ1 / δ2 / δ30° / 8° / 15°构建含发电机内节点的增广导纳矩阵前 3 个节点是机端母线后 3 个节点是发电机内电势节点每一台发电机通过暂态电抗 xd 接到对应机端母线。线路参数已知直接组母线导纳矩阵再用 xd 把内节点耦合进来。import numpy as np # 三机系统数据标幺值 X12, X23, X31 0.20, 0.25, 0.30 # 线路电抗 y12, y23, y31 1/(1j*X12), 1/(1j*X23), 1/(1j*X31) H np.array([6.0, 4.5, 5.0]) # 惯性时间常数秒 Xd np.array([0.20, 0.25, 0.22]) # 暂态电抗 E np.array([1.05, 1.03, 1.04]) # 内电势幅值 delta0 np.deg2rad([0.0, 8.0, 15.0]) # 初始功角弧度 def build_base_ybus(): 机端母线导纳矩阵三角形环网忽略电阻和对地电容 Y np.zeros((3, 3), dtypecomplex) for i, j, y_line in [(0, 1, y12), (1, 2, y23), (2, 0, y31)]: Y[i, i] y_line Y[j, j] y_line Y[i, j] - y_line Y[j, i] - y_line return Y def build_full_ybus(Ym): 把发电机内节点纳入形成6节点增广导纳矩阵 n len(Ym) Y np.zeros((2*n, 2*n), dtypecomplex) Y[:n, :n] Ym for g in range(n): y_xd 1 / (1j * Xd[g]) Y[g, g] y_xd # 机端节点贡献 Y[ng, ng] y_xd # 内节点自导纳 Y[g, ng] -y_xd Y[ng, g] -y_xd return Y这段代码里build_base_ybus 组装的是线路网络的母线导纳矩阵环网结构每一条线路在两个节点上各自产生自导纳和互导纳。build_full_ybus 再把发电机暂态电抗以串联支路形式接入发电机内节点作为独立的恒压源节点这个增广矩阵是后续网络求解的基础。注意 xd 支路只有内节点到机端节点一条路径内节点之间没有直接耦合它们的耦合完全通过机端母线网络实现。3.3 短路故障如何“注入”节点三相短路的导纳修改法时域仿真里短路故障不是“删掉节点”而是“修改导纳”。金属性三相短路时故障点对地阻抗趋近于零相当于在故障节点上并联一个极大导纳。数值上不能让导纳直接取无穷大常见做法是取一个很小的阻抗 Zf 1e-4 j1e-6然后向故障节点的自导纳上加上 1/Zf。这样故障节点的电压被压到接近零电磁功率降到接近零发电机进入加速状态。线路中点故障需要多一步处理。以线路 2-3 中点三相短路为例把线路拆成两段每段电抗取原线路的一半中间新增一个故障节点 4。正常运行时节点 4 是隐藏的故障时刻给节点 4 加接地导纳切除时刻不仅移除该接地导纳还要把整条线路 2-3 从网络中摘除对应 Y23 和 Y32 清零同时从 Y22 和 Y33 的自导纳中减去原线路导纳。这一步最容易出错只清互导纳、不清自导纳导纳矩阵不对称仿真结果会直接飞掉。def apply_fault(Ym, fault_node, Zf1e-4 1j*1e-6): 故障注入在指定机端母线并联小阻抗接地 Yf Ym.copy() Yf[fault_node, fault_node] 1 / Zf return Yf def remove_line(Ym, i, j, y_line): 切除线路互导纳清零自导纳相应修正 Yp Ym.copy() Yp[i, i] - y_line Yp[j, j] - y_line Yp[i, j] 0.0 Yp[j, i] 0.0 return Ypapply_fault 处理的是机端母线短路Zf 取极小值模拟金属性短路的极强接地效果。remove_line 对应保护动作切除故障线路注意自导纳修正和互导纳清零必须成对出现。如果仿真场景是母线故障且保护跳开母联那在 Ypost 里还要额外处理母线所带电源/负荷的转移这里不过度展开工程分析里由故障设置模块统一管理。3.4 电磁功率求解与摇摆方程积分主循环有了导纳矩阵每步积分先求机端母线电压再反推发电机电磁功率。内节点电压已知为 E∠δ机端母线电压通过解网络方程 Ymm·Vm -Ymg·Vg 得到然后算发电机注入网络的电流和电磁功率。def solve_network(Yfull, delta): 给定内电势相角求解机端电压和发电机电磁功率 n len(delta) Ymm Yfull[:n, :n] Ymg Yfull[:n, n:] Ygm Yfull[n:, :n] Ygg Yfull[n:, n:] Vg E * np.exp(1j * delta) # 内电势相量 Vm np.linalg.solve(Ymm, -Ymg Vg) # 机端母线电压 Ig Ygm Vm Ygg Vg # 流入网络的发电机电流 Pe np.real(Vg * np.conj(Ig)) # 电磁功率 return Vm, Pesolver 的核心是五步内电势相量由幅值和功角合成机端电压通过线性方程组求得发电机电流由增广导纳矩阵的耦合分块计算电磁功率取内电势电压乘电流共轭的实部。这个功率表达式已经包含了各机之间的相互影响不需要再单独列出功率角方程。摇摆方程的积分采用经典二阶模型每台机两个状态量功角和角速度。def build_rhs(Y0, Yf, Ypost, t_fault, t_clear, D0.0): 构造微分方程右侧函数按时间阶段切换导纳矩阵 n len(H) omega0 2 * np.pi * 50.0 def ybus_at(t): if t t_fault: return Y0 elif t t_clear: return Yf else: return Ypost def rhs(t, x): delta x[:n] omega x[n:] Y ybus_at(t) _, Pe solve_network(Y, delta) ddelta omega - omega0 domega (omega0 / (2 * H)) * (Pm - Pe - D * (omega - omega0)) return np.concatenate([ddelta, domega]) return rhs def rk4_step(rhs, t, x, dt): 经典四阶龙格库塔单步积分 k1 rhs(t, x) k2 rhs(t 0.5*dt, x 0.5*dt*k1) k3 rhs(t 0.5*dt, x 0.5*dt*k2) k4 rhs(t dt, x dt*k3) return x (dt/6.0) * (k1 2*k2 2*k3 k4)rhs 函数内部按 t 自动选择当前阶段的导纳矩阵故障发生在 t_fault切除发生在 t_clear三个阶段各自对应 Y0、Yf、Ypost。注意初始状态下把 Pm 初始化为故障前的电磁功率这是保证仿真开始不出现虚假功率冲击的关键。rk4_step 是四阶龙格库塔的标准写法每一步调用四次右侧函数步长建议取 0.001 到 0.005 秒。整个主循环模拟故障前的稳态段、故障段和故障后段完整记录每台机的功角轨迹最后按最大相对功角差判稳。Y0 build_full_ybus(build_base_ybus()) Yf build_full_ybus(apply_fault(build_base_ybus(), fault_node2)) Ypost Y0 # 母线故障清除线路保持完整 n len(H) Pm solve_network(Y0, delta0)[1] # 取初始电磁功率为机械功率 x0 np.concatenate([delta0, 2*np.pi*50.0*np.ones(n)]) t_fault, t_clear 0.10, 0.25 dt, t_end 0.002, 20.0 rhs build_rhs(Y0, Yf, Ypost, t_fault, t_clear, D0.5) t_list, angle_list [0.0], [delta0.copy()] x x0 t 0.0 while t t_end: x rk4_step(rhs, t, x, dt) t dt t_list.append(t) angle_list.append(x[:n].copy()) angle_array np.rad2deg(np.array(angle_list)) max_dev angle_array.max(axis1) - angle_array.min(axis1) stable max_dev.max() 180.0 print(失稳 if not stable else 稳定, 最大功角差, round(max_dev.max(), 1), 度)主循环从 0 到 20 秒步长 0.002 秒一共 1 万步三机系统跑起来非常快。故障时刻设在 0.1 秒0.25 秒切除模拟保护装置动作时间 100ms 加断路器开断时间 50ms 的典型组合。判稳逻辑取任意两台机最大功角差的最大值小于 180 度判稳超过则判失稳。3.5 换用标准节点系统从示例走向工程模型上面的三机系统用于理解原理够用真正工程分析一般不手写这个网络。日常工作是拿 BPA 或 PSASP 的卡片数据把网络和发电机参数直接读进来形成大导纳矩阵。标准算例里 IEEE 3 机 9 节点、IEEE 39 节点系统是校验新算法、新脚本最常见的两个网架。手写代码对接标准算例时关键在于把潮流初值吃准。BPA 潮流卡片里的节点电压、发电机出力直接映射成时域仿真初值不要自己算功角更不要凭空拍一组初值。pypower 用户可以用 loadcase 读 IEEE 9 节点数据然后 runpf 取潮流结果再把它转换成语义清晰的字典结构传给时域仿真函数网络求解和积分主循环完全复用上面这套框架。参数表一换代码主体一行不用动这就是分层建模的好处。4. 时域仿真关键参数怎么调切除时间、积分步长与故障阻抗的取值4.1 故障切除时间影响稳定结论最直接的参数短路故障时发电机输出电磁功率接近零机械功率不变转子持续加速。切除越晚加速时间越长功角冲得越高。同一个故障点0.2 秒切除时最大功角差 150 度0.25 秒切除就变成 210 度稳定结论直接翻转。切除时间的工程来源是保护动作时间加断路器固有开断时间一般 220kV 线路是 100ms 左右如果带断路器失灵保护二次跳闸要按 250ms 到 300ms 校核。做运行方式分析时通常不直接把切除时间写死而是扫出一个临界切除时间CCTCritical Clearing Time再对比保护实际动作时间。CCT 大于实际切除时间留有余量CCT 小于实际切除时间说明该故障下系统承受不住这个切除速度要调整保护定值或改变运行方式。扫 CCT 的标准做法是二分搜索下一章给完整脚本。4.2 故障阻抗与故障类型为什么不能直接填 0金属性短路的物理本质是故障点电压为零但数值仿真里给故障节点并联一个阻抗为零的支路导纳直接是无穷大线性方程组求解直接崩。实际金属性短路存在弧光电阻标幺值通常在 0.001 到 0.01 之间所以工程习惯取 Zf 1e-3 j1e-4 或更小既足够接近金属性短路又不至于让矩阵病态。故障类型处理方法备注三相短路节点并联小阻抗接地最严重工况稳定校核首选项两相相间短路复合序网或相域模型需负序网络单相接地复合序网中性点接地方式影响大发展性故障分时段切换故障类型需事件序列表三相短路对系统冲击最大也是暂态稳定标准校核的默认工况。两相和单相故障要构造负序、零序网络手写代码工作量陡增日常直接用商业软件的事件卡片配置原理上都是修改故障节点导纳只是对称分量法要把正负零序网络分开求解再合成。4.3 积分步长和仿真时长数值稳定性与多摆失稳的权衡显式 RK4 在刚性系统里对步长很挑剔。电力系统网络方程本身是代数方程发电机转子方程时间常数在秒级但励磁绕组和网络暂态会把系统变刚。只考虑机电暂态、用经典模型时0.005 秒步长通常没问题一旦引入励磁模型或双馈风机变流器模型建议直接降到 0.001 秒否则高频振荡一旦激发结果直接废掉。仿真时长要覆盖到多摆过程。首摆峰值一般出现在故障切除后 0.3 到 1.0 秒但第二摆、第三摆可能在 2 到 5 秒才达到峰值。阻尼较弱的系统每隔几秒摆一次摆幅缓慢增大20 秒内必然越过 180 度。我一般设 20 秒如果 20 秒内相对角差持续衰减且最大值低于 150 度判稳定的信心才足够。4.4 阻尼系数 D 的影响与整定经典模型里的阻尼系数 D 代表机械阻尼和电气阻尼的合计效果取值很玄学。汽轮机组一般取 1 到 2水轮机组取 0.5 到 1研究性算例经常直接取 0相当于最保守地忽略一切阻尼因素。D 取大了摇摆曲线衰减快容易掩盖系统真实的负阻尼问题取小了曲线长期不衰减又会把本来稳定的系统判成临界。我的习惯是先用 D0 扫 CCT得到保守下界再按实际机组阻尼经验值取 D1 复核几个关键故障两次结果差异大说明结论对阻尼敏感这种情况下要谨慎下结论。5. 多机时域仿真避坑指南五个让结果翻车的常见工程错误5.1 故障阻抗设成 0矩阵奇异与数值爆炸现象仿真跑到故障时刻瞬间电压、功率输出 NaN或者功角曲线出现不合理的跳变日志里能看到 LinAlgError 奇异矩阵报错。原因故障节点并联阻抗为 0对应导纳无穷大导纳矩阵在该节点出现溢出或严重病态。金属性短路在物理上确实是零阻抗但数值计算需要一个小阻抗来代表弧光电阻和数值容差。解决故障阻抗取 Zf 1e-3 j1e-4敏感问题可以再压到 1e-4 j1e-6。不要取 0也不要用 1e-8 这种极端小值否则矩阵条件数过差结果仍然不可信。做灵敏度分析时Zf 变化一个数量级对功角曲线的影响极小说明这个取值是合理的。5.2 初值不自洽仿真开始就出现虚假摇摆现象t0 到 t0.02 秒内功角就快速偏移然后在 0.1 秒故障前已经摆出十几度故障时刻的状态完全偏离工作点。原因机械功率 Pm 和初始电磁功率 Pe 不相等。时域仿真从静态工作点出发发电机转子初始加速度必须为零。很多手写脚本直接拍一个 Pm没考虑初始潮流对应的 Pe仿真开始时转子就带上了初速后面所有结论全部失真。解决初始化时先解一次网络方程把故障前状态的 Pe 算出来再令 Pm Pe。BPA 和 PSASP 里这个动作是自动完成的但在自写代码里是显式步骤。代码写法是Pm solve_network(Y0, delta0)[1]一行代码解决。5.3 切除线路时只改互导纳、忘改自导纳现象故障切除后功角曲线出现异常阶跃或者切除前后功率不平衡曲线明显不连续甚至出现功率凭空增加的情况。原因线路切除涉及两个修改线路两端节点间的互导纳清零同时两个节点的自导纳要扣除该线路贡献的那部分导纳。只清互导纳自导纳还保留着线路的并联支路网络拓扑自相矛盾。解决把线路切除封装成独立函数按“自导纳减、互导纳清”的顺序成对修改每次切线后检查导纳矩阵是否对称。代码里 remove_line 已经体现了这个顺序实际工程里切双回线中的一回、切变压器支路逻辑完全一致。5.4 把参考机当“绝对参考”忽略参考机自身摇摆现象同一套仿真数据用 G1 做参考显示相对功角差最大 160 度判稳换个参考机显示超过 200 度判失稳。两边都有道理但结论打架。原因多机系统没有绝对静止的参考机。任何一台机都在摆动以单机无穷大的习惯选参考机参考机自身的大幅摆动会把真实相对运动淹没。解决改用 COI 参考系计算功角偏差或者直接用任意两机功角差的最大值作为判稳指标。代码里max_dev angle_array.max(axis1) - angle_array.min(axis1)就是取每时刻最大和最小功角的差这个差值不依赖具体参考机结论客观稳定。5.5 仿真时长太短只看到首摆漏掉多摆失稳现象5 秒仿真显示最大功角差 140 度曲线没有越过 180 度结论写“稳定”。把仿真延长到 20 秒发现第 9 秒第二摆冲到 230 度实际是失稳的。原因部分系统故障切除后首摆看起来收敛但故障期间注入的不平衡能量没有完全耗散后续摆动逐渐放大。尤其是阻尼弱、多机振荡模式耦合的系统失稳往往发生在第二、第三摆。解决仿真时长至少取 15 到 20 秒判稳前同时看最大功角差随时间的变化趋势——是持续减小还是波动增大。末尾阶段的转速偏差也要看转速不归零、相对角差振荡不衰减都不能轻易判稳。这个规则对任何多机系统都适用也是最容易因“跑得快”而踩的坑。6. 从功角曲线到稳定裕度临界切除时间的二分搜索与仿真检验6.1 用二分法扫出临界切除时间故障切除时间是运行方式分析里最关心的稳定边界。二分搜索比线性扫描效率高得多先定一个稳定肯定成立的切除时间下界再定一个必失稳的上界反复折半直到区间宽度小于允许误差。def is_stable(t_clear_candidate): t_fault 0.10 Yf build_full_ybus(apply_fault(build_base_ybus(), fault_node2)) x0 np.concatenate([delta0, 2*np.pi*50.0*np.ones(n)]) rhs build_rhs(Y0, Yf, Y0, t_fault, t_clear_candidate, D0.5) x, t x0, 0.0 max_angle 0.0 while t 20.0: x rk4_step(rhs, t, x, 0.002) t 0.002 rel max(x[:n]) - min(x[:n]) max_angle max(max_angle, np.rad2deg(rel)) if max_angle 180.0: return False return True lo, hi 0.10, 0.60 # 下界稳定上界失稳 while hi - lo 0.005: mid (lo hi) / 2 if is_stable(mid): lo mid else: hi mid print(临界切除时间约, round((lo hi) / 2, 3), 秒)二分搜索里每次调用 is_stable 都要完整跑一遍 20 秒仿真三机系统耗时很小几十次循环完全可接受大型系统建议把步长放宽到 0.005 秒做初扫锁定区间后再精扫。出来 CCT 以后拿保护实际切除时间做对比比值大于 1.3 一般认为是可靠接近 1.0 就要警惕。6.2 仿真结果的三条自检路径第一步长收敛性。把步长从 0.002 换成 0.001 和 0.0005 各跑一遍三条功角曲线主峰差在 1 度以内说明积分可靠。第二能量校核。故障切除后的转子动能增量应近似等于加速功率对时间的积分偏差超过 5% 说明积分或功率计算有逻辑错误。第三商业软件对照。同一工况的曲线走向和 CCT 结果和 BPA、PSASP 或 PSS/E 对比功角峰值差在几度内属于正常差十几度要检查模型参数和故障事件时序。我的个人习惯是 CCT 二分搜索永远配一条详细模型复核曲线两边对上了才敢写进分析报告。这套流程踩过的坑不少但把自检步骤固定下来以后翻车的概率低了很多希望帮到你。本文还有配套的精品资源点击获取
