MATLAB三机系统暂态稳定时域仿真：从摇摆方程到故障切除分析

简介这是一份面向电力系统暂态稳定分析学习者的多机系统短路故障时域仿真MATLAB代码包适用于理解三机系统两相短路接地故障后的动态行为与稳定性判断。压缩包仅4KB包含5个m文件覆盖主程序、改进欧拉法数值积分、故障前后模型计算、案例数据读取等核心环节结构清晰便于直接运行和二次修改。目前已吸引395人学习下载。通过该代码包读者可以掌握时域仿真在暂态稳定分析中的实现流程获得发电机转速、电压、电流等关键量的时间序列结果进而分析故障切除时间对系统恢复能力的影响为后续开展更复杂的电力系统保护策略研究提供可复用的仿真基础。1. 多机系统暂态稳定时域仿真这个MATLAB包里到底有什么做电力系统故障分析的工程师对“暂态稳定、多机系统、电力系统故障、时域仿真”这四个词组合在一起应该不陌生一台发电机失稳往往不是它自己的事而是整个网络重新分配功率的问题。上个月我在整理一批旧项目的仿真资料时翻到一个三机系统短路故障时域仿真包里面只有五个 MATLAB 脚本——main.m、Improved_euler.m、calculate.m、readcase.m、data1.m但跑通它能完整模拟一次两相短路接地故障从发生到切除再到判稳的全程。这个包适合两类人一类是刚接触暂态稳定计算、想拿现成代码理解摇摆方程怎么离散求解的学生另一类是手里有实际系统参数、想快速验证故障切除时间对系统功角曲线影响的工程师。它就是一套能直接复现、参数能改、曲线能看的时域仿真框架。2. 从摇摆方程到时域离散暂态稳定仿真前要搞清楚的模型基础2.1 三机系统中经典的发电机二阶模型在多机系统暂态分析里最常用的是经典二阶模型也就是把每台发电机表示成暂态电抗后的恒定电势源加一个惯性转子。这个模型舍弃了励磁调节、调速器和凸极效应换来的是计算量小、物理图像清晰尤其适合短路故障后第一个摆动周期的稳定性分析。它的核心是摇摆方程M_i * (dω_i/dt) P_mi - P_ei dδ_i/dt ω_i - ω_s第一个式子里ω_i是第 i 台发电机的转子角速度ω_s是同步转速P_mi是机械功率P_ei是电磁功率M_i是惯性时间常数。这个方程描述的是发电机转子上不平衡转矩与功角变化的关系故障期间电磁功率骤降而机械功率来不及变转子加速切除故障后如果制动能量能大于加速能量系统就能稳定否则功角持续拉大直到失步。这个模型里最关键的是P_ei的计算它不是孤立的而是所有发电机共同作用的结果。对于 n 机系统第 i 台机的电磁功率为P_ei sum_{j1..n} E_i E_j (G_ij cos(δ_i - δ_j) B_ij sin(δ_i - δ_j))其中E_i是第 i 台机暂态电抗后的电势幅值G_ij和B_ij分别是化简到发电机内节点的节点导纳矩阵第 i 行第 j 列元素的实部和虚部。这就是为什么时域仿真绕不开网络方程的求解——功角一变导纳矩阵对应关系一变注入功率就要重新算。2.2 网络方程与节点导纳矩阵计算故障前后的系统状态时域仿真的每个积分步里都要根据当前网络拓扑重新生成节点导纳矩阵然后消去非发电机节点得到只保留发电机内节点的化简导纳矩阵Y_reduced。calculate.m大概率就在干这件事构建原始导纳矩阵做高斯消去或 Kron 化简最后输出G和B矩阵供摇摆方程使用。故障期间和切除故障后导纳矩阵不一样原因在于网络结构变了。两相短路接地发生在 AB 线路首端时故障点额外引入一条对地支路0.1 秒后保护动作切除故障线路相当于从网络中拿掉一条边。如果不改导纳矩阵算出来的电磁功率和故障实际值完全对不上功角曲线直接失真。我自己在做这类仿真时一般不会在每步都重新化简整个矩阵因为在固定拓扑下导纳矩阵是常数矩阵只在故障发生和故障切除这两个时刻重新生成一次就好。效率高也符合物理过程。2.3 改进欧拉法为什么是这里的主力求解器摇摆方程是一组常微分方程初值问题理论上有四阶龙格库塔等精度更高的方法但这个包用了改进欧拉法也就是二阶龙格-库塔法以Improved_euler.m单独成文件。原因是暂态稳定仿真通常步长取 0.01 秒二阶方法在这个步长下已经有不错的精度而且每步只需两次右端函数求值程序好调试迭代中出现发散也更容易定位是数值问题还是物理失稳。改进欧拉法的标准做法是先用欧拉公式计算预报值再用梯形公式计算校正值。写成 MATLAB 就是function [t, x] improved_euler(f, x0, tspan, h, netSwitch) % f: 状态导数函数 handle, 形参 (t, x, topology_tag) % x0: 初始状态向量 [delta_1..delta_n, omega_1..omega_n] % tspan: [t_start t_end] % h: 积分步长 % netSwitch: 用于在故障前/故障中/故障后切换拓扑的函数 handle t tspan(1):h:tspan(2); x zeros(length(t), length(x0)); x(1,:) x0(:).; for k 1:length(t)-1 % 根据当前时刻选择拓扑状态 tag netSwitch(t(k)); % 预报步骤 k1 f(t(k), x(k,:), tag); % 校正步骤 tag2 netSwitch(t(k)h); k2 f(t(k)h, x(k,:) h*k1, tag2); x(k1,:) x(k,:) (h/2)*(k1 k2); end end这里f是摇摆方程的右端函数返回[d_delta, d_omega]netSwitch根据时间返回拓扑编号比如 0 表示故障前1 表示故障中2 表示故障切除后。这样改进欧拉法的积分过程就和事件切换解耦了主程序里只需要维护一个时间到拓扑的映射。参数h一般取 0.005 到 0.02 秒之间初学建议先用 0.01。注意x矩阵的行是时间列是状态量。排列顺序务必和摇摆方程定义的维度一致否则算出来全是错位数据。3. 把五个脚本串起来从data1.m到main.m的仿真流程3.1 data1.m和readcase.m系统参数到底怎么存data1.m不是函数而是一个被读取的脚本存放整个三机系统的原始参数。常见的内容包括发电机总数、每台机的暂态电抗Xd_prime、惯性常数、机械功率、内电势幅值、节点导纳矩阵原始数据、线路拓扑以及初始功角。以三机系统为例数据结构大致是这样% data1.m 三机系统参数定义 n_engine 3; % 发电机参数 [暂态电抗, 惯性时间常数, 机械功率, 电势幅值] gen_param [ 0.2 30.0 0.9 1.05; 0.25 25.0 0.8 1.02; 0.15 35.0 1.0 1.08 ]; % 母线导纳原始数据格式: [首端节点, 末端节点, 电阻, 电抗, 对地电纳] line [ 1 2 0.02 0.15 0.01; 2 3 0.025 0.18 0.012; 1 3 0.018 0.12 0.008 ]; % 故障线路编号: 第1条(1-2)首端发生两相短路接地 fault_line_idx 1; fault_pos 0.0; % 0线路首端, 1末端 % 初始功角和角速度偏差 delta_init [0.28; 0.42; 0.65]; % rad omega_init [0; 0; 0]; % 偏差值readcase.m的工作就是把这份脚本里的变量刷到仿真程序的空间里。更规范的做法是把它封装成一个返回结构体的函数function sys readcase(fname) % 执行数据脚本并把需要的变量打包到结构体 sys 中 run(fname); sys.n n_engine; sys.gen gen_param; sys.line line; sys.fault_line fault_line_idx; sys.delta_init delta_init; sys.omega_init omega_init; % 这里顺手完成一些派生计算导纳矩阵、化简等 sys.Y build_y(sys); end之所以用脚本存数据而不是用.mat、Excel 或纯文件是因为脚本能直接写表达式比如把标幺值统一成一个数组改动时一眼能看到原始计算。缺点是数据文件不能脱离 MATLAB 环境不过对这个规模的工程完全够用。3.2 calculate.m一次要解哪些方程calculate.m是整个仿真的右端函数输入是当前时刻、当前状态和拓扑标记输出是状态量的导数。它的职责有两个一是根据拓扑标记生成该时刻的化简导纳矩阵二是按摇摆方程计算每台发电机的电磁功率和角速度变化率。逻辑上很像这样function dx calculate(t, x, tag, sys) % tag: 0故障前稳态, 1故障中, 2切除后 n sys.n; delta x(1:n); omega x(n1:end); % 按拓扑选择导纳矩阵, 这里只做示意 switch tag case 0 Yred sys.Y_normal; case 1 Yred sys.Y_fault; case 2 Yred sys.Y_post; end G real(Yred); B imag(Yred); E sys.gen(:, 4); Pm sys.gen(:, 3); M sys.gen(:, 2); Pe zeros(n, 1); for i 1:n for j 1:n Pe(i) Pe(i) E(i)*E(j) * (G(i,j)*cos(delta(i)-delta(j)) ... B(i,j)*sin(delta(i)-delta(j))); end end d_delta omega - sys.omega_s; % omega_s 为同步角速度 d_omega (Pm - Pe) ./ M; dx [d_delta; d_omega]; end这段代码里的d_omega是实现了两段机械功率Pm在仿真过程中保持常数这是经典模型的假设电磁功率Pe从化简导纳矩阵和当前功角算出。故障期间由于导纳矩阵中增加了故障支路Pe大幅下降d_omega变大转子加速。代码逻辑本身不复杂但要注意矩阵维度G和B是n x nE是n维列向量cos和sin里的角度差如果用弧度制初始功角也要是弧度。用角度制也行但所有公式内部都得统一否则算出来的功率会周期性跳变这是新手最容易翻车的地方。3.3 故障模拟与0.1s切除逻辑时域仿真中的事件切换故障事件在仿真里的处理方式不是“到点改一个开关量”而是重新生成网络结构。两相短路接地发生在 AB 线路首端意味着故障点在节点 A 附近额外引入一条接地支路这条支路的阻抗按两相短路接地的复合序网计算正序网络中等效为短路阻抗并联接到大地。对应到程序里就是变更原始导纳矩阵的故障相关节点对角元然后再次化简得到Y_fault。切除逻辑同样是一个时刻触发的事件。0.1 秒后保护动作切除故障线路 AB相当于把这条线路的串联阻抗从导纳矩阵中剔除。注意这里“0.1秒后”指的是故障开始后的延时。如果仿真在 0 秒启动、0.1 秒发生故障那么切除时刻一般是t_fault t_clear即 0.2 秒。很多初学朋友会在main.m里把切除时间直接写成 0.1那实际成了故障刚发生就切除跟题目设定不一致结果曲线自然不同。事件切换的实现方式在main.m里常见是写一个匿名时间映射函数t_fault 0.1; % 故障发生时刻 t_clear 0.2; % 故障切除时刻 topologyTag (t) (t t_fault)*0 ... (t t_fault t t_clear)*1 ... (t t_clear)*2;这个映射式的写法很紧凑数值上用逻辑表达式返回 0、1、2 三种状态积分器每走一步都会调用它一次来决定用哪套导纳矩阵。比堆if-else清晰也方便改成多级故障序列比如重合闸、切机等等。4. 跑起来与参数调整短路故障仿真的可复现实验4.1 第一次运行main.m的完整步骤拿到这个包第一件事是确认 MATLAB 当前路径在解压目录下然后把五个脚本放进同一个文件夹。打开main.m依次检查三组变量仿真总时长t_end、步长h、故障与切除时间。第一次跑建议保持原始参数不动直接执行% main.m 中的调用过程示意 clear; clc; % 读取系统数据 sys readcase(data1.m); % 仿真控制参数 t_end 5.0; % 总仿真时长, 单位秒 h 0.01; % 积分步长, 单位秒 t_fault 0.1; % 故障发生时刻 t_clear 0.2; % 保护动作切除时刻 % 初始状态 x0 [sys.delta_init; sys.omega_init]; % 定义拓扑切换逻辑 netSwitch (t) (t t_fault)*0 ... (t t_fault t t_clear)*1 ... (t t_clear)*2; % 调用改进欧拉法 [t, x] improved_euler((t, x, tag) calculate(t, x, tag, sys), ... x0, [0 t_end], h, netSwitch); % 输出结果 delta x(:, 1:sys.n); figure; plot(t, delta*180/pi); xlabel(时间/s); ylabel(功角/deg); legend(G1,G2,G3); grid on;这段代码里有几个点值得说明。readcase(data1.m)返回结构体sys里面已经包含化简导纳矩阵和初始功角所以main.m本身不用关心导纳矩阵怎么形成。improved_euler接收的右端函数用了匿名函数把sys捕获进去省得全局变量满天飞。输出里的x是(时间点数) x (2*n)矩阵前半列是功角后半列是角速度偏差画图时只取前n列即可。4.2 时间步长、故障时间和切除时间这三个参数怎么定时间步长h决定仿真的离散分辨力。三机系统稳定后转子相对功角摆动频率大概是 0.5 到 2 赫兹一个周期对应 0.5 到 2 秒步长取 0.01 秒可以保证每周期有 50 到 200 个采样点画出曲线的毛刺很小。如果想把计算速度提上来可以放宽到 0.02 秒如果功角曲线首摆细节存在明显尖锐转折就缩到 0.005 秒。故障时刻t_fault的取值影响的是故障前稳态阶段的长度。你希望系统先在原平衡点稳定运行一小段时间再跳变通常取 0.05 到 0.3 秒。取值太小会让初始冲击和故障冲击叠加曲线起始段不好看取值太大会白白增加计算量并且如果初始功角不是精确平衡点积累的偏移会影响故障后结果。切除时间t_clear是关键中的关键。题目里明确是 0.1 秒后切除故障线路那么折算到仿真时间轴上是t_fault 0.1。如果要探索系统稳定边界就把t_clear从t_fault开始逐步后移每步增加 0.02 秒观察功角曲线什么时候从“摆动后衰减”变成“单调增大”。这就是求临界切除时间的笨办法但也最可靠。4.3 判稳不只看功角需要输出的关键曲线很多初学朋友画完功角曲线看到三条线在抖动就一脸懵。判稳的正确做法不是看绝对功角而是看相对功角差。三机系统一般选 1 号机做参考机计算delta_2 - delta_1和delta_3 - delta_1观察这两个差值是否在某一时刻开始单调增大并超过 180 度。% 相对功角差 delta_rel2 (delta(:,2) - delta(:,1)) * 180/pi; delta_rel3 (delta(:,3) - delta(:,1)) * 180/pi; figure; plot(t, delta_rel2, t, delta_rel3); legend(delta2-delta1, delta3-delta1); xlabel(时间/s); ylabel(相对功角/deg); grid on;如果相对功角差曲线在振荡几次后趋于恒定或衰减到同步值系统暂态稳定如果没过多久就直奔 180 度甚至超过 360 度说明发电机已经失步需要更早切除故障或调整线路参数。角速度偏差曲线也可以看失稳时对应的那台发电机角速度会持续偏离同步转速。对于两相短路接地这种非对称故障仿真的物理过程比三相短路温和切除时间窗口相对宽裕但分析手段完全一样。5. 避坑指南时域仿真里最常见的问题与排查思路5.1 功角曲线发散不是数值发散而是模型失稳现象仿真跑到最后某台发电机的功角相对参考机单调增大一路冲到上千度曲线看起来像是螺旋上升。原因这大概率不是数值方法发散而是切除故障时间过晚系统真的失稳。解决先缩短t_clear再看曲线如果步长减半之后曲线几乎不变基本判断为物理失稳要缩短故障切除时间或修改线路参数。我判断这个有个习惯把t_clear设成等于t_fault即故障刚发生就切除如果依然发散才考虑数值问题。5.2 改进欧拉法步长太大导致振铃现象功角曲线在稳态阶段出现高频等幅振荡看起来像叠加了一个几十赫兹的正弦波而不是平滑的摆动。原因步长太大二阶 Runge-Kutta 的频率响应在中高频段产生了数值振荡三机系统模型的自然振荡频率在暂态过程中瞬时可能很高。解决把h从 0.01 降到 0.005 再跑一遍如果振荡幅度随步长减半明显减小就是步长问题。另外检查一下d_delta是否用了omega与参考机角速度的差值如果直接用绝对角速度也会导致高频分量被放大。5.3 readcase读取失败路径与数据格式问题现象执行readcase(data1.m)报错Undefined function or variable或Invalid file name。原因最常见是 MATLAB 当前路径不在脚本所在目录或者data1.m里定义的变量名和readcase.m中run之后引用的变量名拼写不一致。解决先cd到资源目录或者在main.m第一行相对路径设置好然后打开data1.m和readcase.m逐个对变量名重点是n_engine和n、gen_param和gen这类的缩写映射。我的经验是把数据文件名和内部变量名固定下来以后换系统只需要改数据文件readcase一行不动。5.4 节点编号和参考机约定容易搞混现象仿真结果里 1 号机和 2 号机的功角差初始值不为常数而是在故障前阶段就在漂移。原因初始功角量测的参考节点和你编程时取的参考机不一致。比如数据里delta_init可能是相对于无穷大母线的绝对角度但程序里相对角计算用的是 1 号机作为参考导致初始时刻相对功角差就不是平衡值。解决在readcase阶段统一约定——所有初始功角都转换为相对 1 号机的差值并将 1 号机的初始角设为 0。三机系统里参考机的选择不影响稳定性结论但影响曲线形态务必让数据文件和计算逻辑使用同一个参考机。5.5 故障前初值算错潮流初值对暂态结果的影响现象故障发生前功角曲线在 0 到 0.1 秒内就有明显波动甚至向下漂移。原因初始功角不是真正稳态平衡点仿真程序把初始状态当成平衡点但右端函数在那个点导数为非零于是故障还没来系统自己先动起来了。解决在启动暂态仿真前先做一次潮流或稳态求解把delta_init收敛到满足P_mi P_ei的点。如果包里没有潮流函数可以自己写几行牛顿迭代或者直接把数据里的初始功角代入calculate.m看看d_omega是否接近 0不为 0 就要修正。注意这种“自然漂移”不是仿真步长能解决的它反映的是初值不匹配属于建模问题。6. 把仿真结果用起来从功角曲线到临界切除时间定位前几章把仿真跑通后最直接的应用是找到临界切除时间 CCT——这是暂态稳定分析里最常用的一个量化指标。做法很简单固定t_fault从很小的t_clear开始逐步增大切除时间每跑一次记录系统是否稳定稳定与失稳的分界点就是 CCT。手动人肉试太慢可以用二分法自动找% 二分法搜索临界切除时间 t_low 0.1; % 初值取故障时刻, 必然稳定 t_high 0.5; % 一个足够大且肯定不稳定的时间 for iter 1:15 t_test (t_low t_high) / 2; sys.clearingTime t_test; [t, x] run_simulation(sys); % 重新运行完整仿真 delta_rel max(abs(x(:,[2 3]) - x(:,1)), [], 2); if max(delta_rel) pi % 设定不超过180度作为稳定判据 t_low t_test; else t_high t_test; end end cct (t_low t_high) / 2;这段二分法的稳定判据用的是“相对功角差最大值是否超过 180 度”这个阈值对应静态稳定极限是经典模型下最常用的判据。15 次迭代可以把区间缩窄到 0.00003 秒量级足够工程使用。如果你在算一个保护定值拿到 CCT 之后再乘一个 0.8 的安全系数就是建议的最大切除时间。我还喜欢在这个基础上再画一张“切除时间-最大相对功角”的曲线。横轴是t_clear纵轴是摆开过程中的相对功角峰值可以看到曲线在 CCT 附近有一个接近垂直的上升段这个陡峭拐点比单纯一个数值更能说明系统稳定裕度。从那次翻出这个包之后我每次做故障仿真都强制走一遍流程先用原始参数确认代码能跑再把t_clear缩短一半确认曲线形态正确然后才放心去改系统参数。这套习惯帮我避开了好几回初值错误和数据格式错误。希望这个资源也能帮你在暂态稳定分析上少走些弯路。本文还有配套的精品资源点击获取