刚接触电力系统暂态稳定分析时我最大的感觉是课本上的公式看得懂但真正要算一个具体系统、验证一个故障场景时手算几乎无从下手。功角摇摆方程是非线性的故障瞬间电气参数跳变靠解析方法只能处理极简单的理想模型。后来我把重心转到 Matlab 编程加 Simulink 仿真才感觉真正打开了暂态稳定分析的大门。这套组合的好处很直接Matlab 负责把数学模型算清楚Simulink 负责把物理系统搭出来看动态过程两者互相验证既能深入理解理论又能贴近工程实际。这篇文章就围绕“用 Matlab Simulink 做电力系统暂态稳定性分析”展开从数学模型推导、编程计算、仿真建模到临界切除时间验证把完整流程走一遍。适合正在学习电力系统分析的在校学生、刚入行的电力工程师以及任何对机电暂态仿真感兴趣、想用工具验证理论的人。我尽量把每一步的原理和实操都讲透包括我踩过的坑和调试经验方便你直接照着做。1. 暂态稳定分析到底在分析什么1.1 从“发电机能不能稳住”说起电力系统暂态稳定通俗讲就是系统遭遇大扰动最常见的是短路故障、大容量负荷切除或发电机跳闸之后系统中各发电机转子能不能继续保持同步运行。如果能系统是暂态稳定的如果某些发电机相对其他发电机越摆越大最终失步系统就暂态失稳了。这种稳定问题的时间尺度通常在零点几秒到几秒之间正好是机电暂态过程的时间范围所以也常称为功角稳定或第一摆稳定。很多人会把暂态稳定和静态稳定搞混。静态稳定研究的是系统在某一运行点受到微小扰动后能否回到原来的平衡状态用的是线性化方法本质上是看雅可比矩阵的特征值暂态稳定研究的是大扰动下的非线性动态行为发电机功角可能摆出很大的角度线性化不再适用。另一个容易混淆的概念是动态稳定它关注的往往是加了各种控制器励磁调节器、PSS、调速器等之后系统在小扰动或中等扰动下的阻尼特性。三种稳定的时间尺度、分析方法和数学工具都不一样做暂态稳定时要用的是非线性微分方程数值积分。暂态稳定分析为什么重要因为电力系统的安全防线里继电保护动作切除故障通常在几十到几百毫秒内完成而发电机的功角摆开则需要一秒甚至更久。这个时间差里系统能不能扛得住第一摆、能不能在故障切除后逐渐衰减振荡稳定下来直接决定事故是止步于“跳闸”还是升级为“大面积停电”。所以工程上做安全稳定分析时暂态稳定校核是必做项目而且要遍历各种预想故障。1.2 为什么选单机无穷大系统作为切入点真实电力系统动辄几十台、上百台发电机分析起来极其复杂。学习阶段最好的切入点是单机无穷大系统Single Machine Infinite BusSMIB。它由一台发电机经过升压变压器、双回输电线路连接到无穷大母线构成无穷大母线代表一个容量无限大、电压和频率都恒定的理想系统。SMIB 系统虽然结构简单但包含了暂态稳定分析的全部核心要素发电机有惯性、有机械功率输入和电磁功率输出、输电网络有电抗、故障会改变网络拓扑。它足够简单可以用解析方法等面积法则求出理论解又足够通用多机系统的许多结论可以从 SMIB 系统推广而来。换句话说把 SMIB 系统彻底搞懂再去做多机系统就会轻松很多。用 Matlab 和 Simulink 来做这个课题正好发挥两者的特长。Simulink 里搭一个 SMIB 系统模型非常直观三相电源、变压器、输电线路、故障模块都是现成的图形化搭建能避免很多编程细节的干扰Matlab 脚本则用来做批量参数计算、数值积分、批量绘图和结果后处理。另外Matlab 的数值计算能力很强等面积法则中需要的非线性方程求解、临界切除时间计算几行代码就能搞定这是手算完全做不到的效率。2. Matlab 编程实现暂态稳定计算2.1 摇摆方程暂态稳定的数学内核SMIB 系统的暂态稳定分析建立在发电机的转子运动方程上工程上通常叫摇摆方程Swing Equation。写成标幺值形式2H/ωs · d²δ/dt² Pm - Pe其中 H 是发电机组的惯性时间常数秒ωs 是同步角速度δ 是发电机转子角功角Pm 是原动机输入的机械功率Pe 是发电机输出的电磁功率。这个方程的本质是牛顿第二定律在旋转机械上的应用转子加速功率等于机械功率和电磁功率之差差值大于零则功角增大小于零则功角减小。电磁功率 Pe 在理想情况下满足Pe (EV/X) · sinδ Pmax · sinδ其中 E 是发电机暂态电动势近似恒定V 是无穷大母线电压X 是发电机到无穷大母线之间的总电抗。注意故障发生前、故障期间、故障切除后网络拓扑不同对应的总电抗 X 不同所以 Pmax 有三个不同的值。这是暂态稳定分析最核心的要点之一系统在故障过程中经历了三个阶段每个阶段的功率特性曲线都不一样。一阶微分方程组形式为dδ/dt ω - ωs dω/dt (ωs/2H) · (Pm - Pe)这样就把二阶方程转化成了两个一阶方程方便用数值积分方法求解。2.2 用 ODE45 求解功角摇摆曲线求解摇摆方程最直接的方式是用 Matlab 内置的 ODE45 函数。ODE45 是基于龙格-库塔法的自适应步长求解器对这类非线性常微分方程有很好的适应性不需要手动选择步长精度也有保障。我先写一个简单的脚本定义状态方程函数然后调用 ODE45 进行积分。先定义系统参数。假设一个典型的 SMIB 模型发电机暂态电抗 Xd 0.3标幺值变压器电抗 Xt 0.1双回线每回电抗 Xl 0.4总电抗在故障前为 X1 Xd Xt Xl/2 0.3 0.1 0.2 0.6。系统额定运行时发电机向无穷大系统输送有功功率 Pm 0.8无穷大母线电压 V 1.0暂态电动势 E 1.2假设恒定。机械功率在暂态过程中保持恒定这是经典暂态稳定分析的常用假设。正常运行时最大电磁功率Pmax1 EV/X1 1.2×1.0/0.6 2.0初始功角δ0 arcsin(Pm/Pmax1) arcsin(0.8/2.0) arcsin(0.4) ≈ 23.58°换算成弧度约 0.4115 rad接下来定义三阶段的 Pmax。假设 t0.1s 发生三相短路故障故障点在输电线路首端。故障期间发电机到无穷大系统的电气联系几乎被切断电磁功率 Pe ≈ 0严格说是经故障点接地的少量功率工程上常近似为 0。t0.25s 时保护动作切除故障线路切除一回线此时系统变为单回线运行总电抗 X3 Xd Xt Xl 0.3 0.1 0.4 0.8所以Pmax3 EV/X3 1.2/0.8 1.5用 Matlab 求解时我采用分段积分的方式第一段从 0 到 0.1sPmax 2.0第二段从 0.1s 到 0.25sPmax 0第三段从 0.25s 到 3sPmax 1.5。每一段的起始状态是上一段的结束状态这样逐段积分就能得到完整的功角曲线。关键代码如下function dydt swing_equation(t, y, Pmax, Pm, H, ws) delta y(1); omega y(2); Pe Pmax * sin(delta); dydt [omega - ws; (ws / (2*H)) * (Pm - Pe)]; end主脚本中分段调用% 系统参数标幺值 H 5.0; % 惯性时间常数秒 ws 2*pi*50; % 同步角速度rad/s Pm 0.8; % 机械功率 Pmax1 2.0; % 故障前最大电磁功率 Pmax0 0.0; % 故障期间最大电磁功率 Pmax2 1.5; % 故障切除后最大电磁功率 % 初始功角 delta0 asin(Pm / Pmax1); y0 [delta0; ws]; % 分段积分 t1 [0 0.1]; [t1_out, y1_out] ode45((t,y) swing_equation(t, y, Pmax1, Pm, H, ws), t1, y0); ...积分完成后将功角 δ 从弧度转换为角度用 plot 函数绘制功角随时间变化的曲线。观察曲线就能判断系统是否稳定如果功角在达到最大值后开始回摆说明系统能保持稳定如果功角单调递增穿越 180°则系统失稳。2.3 计算过程的关键细节数值积分时有一个容易被忽略的细节故障瞬间和切除瞬间Pmax 发生跳变微分方程的右侧不连续。ODE45 是自适应步长算法在间断点附近可能步长取得不合适导致结果偏离真实解。解决方法是强制在事件发生时刻把积分分段也就是上面代码里把 [0, 0.1]、[0.1, 0.25]、[0.25, 3] 分开积分。这样每个区间的方程都是光滑的误差可控。另一个细节是初始功角的选取。很多人直接取 δ0 0 或者随便给个值这是错误的。初始功角必须满足故障前稳态运行条件即 Pm Pmax1 · sinδ0所以 δ0 arcsin(Pm/Pmax1)。如果初值不对仿真开始阶段就会产生虚假的振荡影响对稳定性的判断。还有一点功角 δ 的单位。状态方程里角速度用 rad/sPmax 用标幺值角度自然用弧度。绘制曲线时要注意单位换算否则图像上的角度值看着很奇怪。建议统一用弧度计算绘图时再乘以 180/pi。3. Simulink 搭建 SMIB 暂态稳定仿真模型3.1 为什么这一步要换用 SimulinkMatlab 编程解决的是“纯数学”的暂态过程求解但它有一个比较大的局限模型是抽象的。当我想观察故障瞬间三相电压电流波形或者想加入更复杂的保护逻辑、励磁系统模型时纯编程的建模成本会变得很高。Simulink 的优势在于可视化建模尤其是电力系统领域Simscape Electrical 组件库提供了完整的电力元件模块三相电源、变压器、输电线路、故障发生器等直接拖拽就能用。Simulink 模型还能直观地展示系统拓扑故障发生在哪里、切除哪条线路一目了然。更重要的是Simulink 里发电机的内部模型可以选得比较复杂考虑励磁、PSS、调速器这在纯编程里很难实现。因此我习惯的组合方式是Matlab 编程做理论验证和参数计算Simulink 搭更接近实际的系统模型做仿真验证。3.2 在 Simulink 中从零搭建 SMIB 模型打开 Matlab在命令行输入 simulink 回车新建一个空白模型。然后在库浏览器中找到 Simscape Electrical Specialized Power Systems 库这里就是电力系统仿真所需模块的大本营。先搭主干网络。从 Electrical Sources 库拖入一个 Three-Phase Source 作为无穷大母线这时把它的容量设得非常大比如短路容量 10000 MVA电压设为 220kV频率 50Hz就相当于一个理想的无穷大母线。从 Elements 库拖入 Three-Phase Transformer (Two Windings) 作为升压变压器把变比设置为发电机出口电压和输电电压的比值比如 13.8kV/220kV。输电线路用 Three-Phase Series RLC Branch 表示在参数里设置正序电阻和电抗值。发电机用 Machines 库里的 Synchronous Machine 模块选择“Fundamental”模型设置额定功率、额定电压、惯性常数、电抗等参数。为了模拟故障从 Elements 库拖入 Three-Phase Fault 模块放在输电线路首端。双击模块设置故障类型为 Three-Phase Fault三相短路故障时序用外部控制信号控制这样就可以在仿真中用阶跃信号精确控制故障发生和切除时刻。连接好电路之后还要在模型中加入 Powergui 模块拖入一个 Powergui 就行不用连到电路里。Powergui 是电力系统仿真的核心用来设置仿真类型连续/离散、求解器、潮流初值等。这里选 Continuous求解器选 ode23tb这个求解器对电力电子和开关突变场景更稳定。3.3 发电机功角的提取方法要观察暂态稳定分析中最关心的量——发电机功角可以有两种做法。一种是用同步电机模块自带的测量输出。Synchronous Machine 模块的测量端口输出一系列信号包括转子角、转速、电磁功率等。但要注意它输出的转子角是相对于转子初始位置的电角度不一定是和无穷大母线之间的功角差。严格的功角定义是发电机 q 轴电动势和无穷大母线电压之间的相位差需要经过换算。另一种更直观的方法是利用 Three-Phase V-I Measurement 模块分别测量发电机端电压和无穷大母线电压的相角然后用加减运算得到两者的相角差再减去变压器和线路的固有相移这个在仿真结果中自然包含就能得到功角曲线。实际操作中我通常用 PLL锁相环模块分别追踪两个电压的相位再求差值。这种方式更接近实际工程中通过 PMU同步相量测量装置测量功角的方式概念上很直观。Simulink 中测量到的功角信号需要转成角度制再送到 Scope 观察。仿真的时间步长设置要注意如果设得太大开关事件故障发生和切除时刻的响应会失真如果设得太小仿真速度会很慢。一般暂态稳定分析关注的是 0 到几秒的机电暂态最大步长设置在 0.0001s 到 0.001s 之间比较合适。3.4 仿真波形解读模型搭建完成后设置故障时序t0.1s 故障发生t0.25s 故障切除仿真时长 3s。运行仿真后观察功角曲线。正常情况下会看到功角从初始值 23.58° 开始故障发生后快速增大因为电磁功率突降到几乎为零机械功率大于电磁功率转子加速切除故障后电磁功率恢复但线路少了一回Pmax 下降功角继续增大一段时间后达到最大值然后在减速功率作用下回摆并围绕新的平衡点做衰减振荡最终稳定。这个波形特征和第二节里 Matlab 编程计算的结果应该一致。如果仿真参数设置得合理两条曲线几乎可以重合。这套“编程 仿真”互验的方法正是这个项目最有价值的地方Matlab 算出来的结果和 Simulink 仿出来的结果互相印证说明模型搭建正确、编程正确、理论理解正确。三者对上了你才可以放心去做更复杂的系统分析。4. 极限切除角与临界切除时间理论计算与仿真验证4.1 等面积法则的应用暂态稳定分析的一个典型工程问题是故障必须在多长时间内切除才能保证系统不失稳这个时间就是临界切除时间Critical Clearing TimeCCT对应的故障切除角度叫极限切除角。工程上继电保护的动作时间必须小于临界切除时间系统才是安全的。等面积法则Equal Area CriterionEAC是求解单机无穷大系统极限切除角的理论工具。它的物理基础是只要加速面积等于减速面积系统就能恢复到稳定运行状态如果可用减速面积小于加速面积系统将失稳。用公式表示极限切除角满足∫(Pm - Pe_fault)dδ ∫(Pe_after - Pm)dδ其中左边是故障期间的加速面积从初始功角 δ0 积分到极限切除角 δcr右边是故障切除后的减速面积从 δcr 积分到最大摇摆角 δmax。而 δmax 又取决于故障切除后功率曲线上能与 Pm 相交的最大角度也就是 π - arcsin(Pm/Pmax_after)。把电磁功率表达式代入可以写出具体积分表达式。以我用的参数为例Pm 0.8故障期间 Pmax0 0切除后 Pmax2 1.5初始功角 δ0 0.4115 rad。等面积方程变成一个包含三角函数的非线性方程手算比较麻烦但用 Matlab 的 fsolve 函数几行代码就能解出来。4.2 用 Matlab 计算临界切除时间计算临界切除时间还有一个更通用的方法二分搜索。思路是给定一个切除时间 tc仿真得到最终的功角曲线判断系统是否稳定功角是否越过失稳门槛。然后不断调整 tc用二分法逼近临界值。具体实现思路% 二分搜索临界切除时间 tc_low 0.1; % 确信稳定的切除时间 tc_high 0.5; % 确信失稳的切除时间 for i 1:20 tc_mid (tc_low tc_high) / 2; % 用 tc_mid 作为故障切除时间运行暂态仿真 stable simulate_stability(tc_mid); if stable tc_low tc_mid; else tc_high tc_mid; end end这里的 simulate_stability 函数可以直接调用第二节的 ODE45 分段积分代码把故障切除时间作为参数传入稳定判据是功角最大值是否超过某阈值通常取 120°180° 之间的值工程上常用 100° 作为摆开角的警示值。也可以用等面积法则先求出极限切除角再通过数值积分故障期间的摇摆方程得到临界切除时间。实际操作中两种方法可以对比验证。等面积法属于理论解计算速度快但只适用于单机系统二分搜索法本质上是数值实验思路可以推广到多机系统。两种方法结果互相吻合说明计算流程没问题。4.3 Simulink 中的对应验证在 Simulink 模型里把 Three-Phase Fault 模块的切除时刻设成不同值运行仿真并观察功角波形可以直观看到切除时间小于临界值时功角摆到某个最大值后回摆系统稳定切除时间大于临界值时功角单调增大系统失稳。这个“分子与分母”的边界状态就是临界切除时间。我试过用 Simulink 来标定 CCT分别设切除时间为 0.30s、0.32s、0.34s、0.36s 等运行仿真看功角曲线。比如在本例参数下Matlab 编程计算的临界切除时间约为 0.365s。Simulink 里设 0.36s 时功角回摆设 0.38s 时功角发散临界值就在这个区间内。两者误差在 0.02s 以内对于工程精度已经非常可观。这一步很有意义因为它让你把三种能力串起来了理论分析等面积法则给出可解析求解的极限切除角、编程计算ODE45 数值积分和二分搜索给出临界切除时间、仿真验证Simulink 可视化的功角摇摆过程。无论哪一种单独用都很难建立对暂态稳定的直观理解三者结合才能真正形成闭环。5. 常见问题与排查技巧实录5.1 仿真发散问题的处理做电力系统仿真仿真结果发散是新手最常遇到的问题。现象是波形直接飞掉数值变成 NaN 或者 InfScope 里的曲线一路冲上天。原因不外乎几种求解器不合适、步长过大、模型本身存在代数环问题。我建议优先检查求解器设置。对于包含开关事件故障投入/切除的电力系统模型ode23tb 或者 ode15s 这类刚性求解器通常比 ode45 更合适。然后把最大步长调小比如从默认值改成 0.0001s。如果原来是离散仿真试着改成连续仿真试试。还有一种情况是 Powergui 的仿真类型选错了建议选 Continuous。如果仿真发散发生在故障切除瞬间大概率是代数环或者数值振荡问题。可以在故障模块的控制信号路径中加入一个 Memory 模块或者 Unit Delay 模块来打断代数环或者在故障切除后用一个小的时间常数模拟断路器弧隙动态过程避免瞬间开路带来的数值冲击。5.2 初始工作点对仿真结果的影响这是一个比较容易踩坑的地方。Simulink 中同步电机模块默认的初始状态可能是零功率运行如果直接开始仿真从空载状态突然加上机械功率会产生非常大的暂态振荡这和我们要分析的暂态稳定问题完全不是一回事。解决办法是用 Powergui 的 Load Flow 工具做潮流初始化。在 Powergui 模块上右键选择 Powergui → Tools → Load Flow设置发电机出力、母线电压幅值等运行条件让仿真从正确的稳态工作点开始。这样仿真一开始各电气量就处于故障前的稳态值故障发生后的动态过程才是真正要分析的暂态过程。Matlab 编程计算的初值也要注意。前面说了 δ0 必须由 Pm Pmax1 · sinδ0 反解这是一个经常被忽略但极其关键的细节。如果 δ0 给错后面的加速面积、减速面积全都不对临界切除时间算出来自然也是错的。5.3 标幺制和单位制混乱电力系统分析领域习惯用标幺值per unitSimulink 的同步电机模块里面也能设置基值。最怕的就是混用发电机用标幺值线路用欧姆电压用 kV结果算出来电磁功率和机械功率对不上功角自然算不对。我习惯的做法是所有电气参数在 Matlab 脚本里先统一换算成标幺值Simulink 里能填标幺值的模块直接填标幺值比如同步电机的电抗、惯性常数需要填国际单位欧姆、亨利的模块比如 RLC 线路再用基值换算回去。基值换算的公式是Zbase Vbase²/SbaseLbase Zbase/ωsCbase 1/(ωs·Zbase)。把这些写成一个小工具函数每次建模时调用能避免大量低级错误。5.4 常用故障排查速查表现象可能原因解决办法功角曲线高频振荡且幅度很大初始工作点不对仿真从非稳态启动用 Powergui 潮流初始化故障切除瞬间波形跳变严重求解器步长过大或代数环问题减小步长、用 ode23tb、在控制回路加 Memory功角看起来很小且不摆动Pmax 计算错误总电抗设置不对检查各阶段电抗值确认线路模型是否接入正确仿真结果和 Matlab 编程结果偏差大Simulink 中考虑了励磁或阻尼绕组而编程模型没考虑对比时统一模型复杂度或编程中也加入阻尼项断路器切除后电流不为零三相 Fault 模块的切换逻辑有误检查 Fault 模块 Transition 时间和控制信号时序这些问题的排查思路有一条主线先确认稳态工作点正确再检查动态模型参数最后查看故障时序。按照这个顺序逐项排除大部分仿真问题都能解决。5.5 一些实操心得我做这个课题最大的体会是Matlab 编程和 Simulink 仿真不是二选一的关系而是互相成就。纯编程时数学模型很清晰但容易忽略物理意义纯仿真时模型很直观但容易变成“黑盒实验”不知道为什么这个参数要这样设。两者结合才会逼迫你同时理解数学和物理。另外建议仿真时长不要设太短。初做时我只仿真 1s结果只看到了首摆没看到功角衰减振荡的收敛过程。把仿真时长延长到 3~5s你才能真正看到阻尼效应也才能判断这个系统是单调稳定还是振荡稳定的。功角摇摆曲线不是一个简单的上升下降而是一个包含多个振荡周期的衰减过程看到完整的曲线你对暂态稳定动态过程的理解会上一个台阶。最后给一个小技巧无论 Matlab 还是 Simulink参数不要硬编码在模型里尽量在脚本中定义变量传入工作区Simulink 模型通过变量名引用。这样批量扫参数时比如扫不同切除时间、不同出力水平只需在脚本里改循环模型自动跟着变。这个习惯能为后续做多机系统、做参数灵敏度分析省下大量时间。