尧图网站开发 
前阵子帮一个电网规划项目算静态稳定极限被同事问了一句“能不能一小时给我出一张功角曲线图”。说实话很多人一听电力系统静态稳定性仿真脑子里全是大型电网模型、时域故障扫描、一堆PSS参数越想越复杂。但我真正把这个问题做顺手之后发现Matlab加Simulink这一套组合完全可以用一条清晰的主线串起来先把物理过程翻译成数学方程再用脚本把特征值解出来最后用Simulink时域仿真去验证和展示。这篇文章我把这条路径完整走一遍从建模思路、Matlab脚本到Simulink模型搭建、PV曲线求取再到参数扫描和踩坑经验适合正在做课程设计、毕业设计或者刚接触电力系统稳定分析的工程师参考。1. 静态稳定性分析到底算的是什么账1.1 功角稳定一个“两个转子”的故事电力系统静态稳定说白了就是看系统在某个运行点上受到一个小扰动之后能不能自己回到原来的稳态。这里的关键词是“小扰动”比如负荷波动零点几个百分点、线路阻抗的缓慢变化而不是雷击短路那种大扰动。大扰动对应的是暂态稳定小扰动对应的是静态稳定这是两套完全不同的分析思路很多人一开始就把这两个概念混在一起导致仿真设计跑偏。功角稳定的物理图像其实特别简单。一台同步发电机经线路接到无穷大母线可以把发电机等效成一个旋转的电动势源无穷大母线看成电压幅值和频率都恒定的源。两个源之间通过电抗相连它们之间有一个相位差δ也就是功角。发电机原动机输入的机械功率和送出去的电磁功率必须平衡而电磁功率随功角变化的规律就是功角特性P_e (E * V / X_Σ) * sin(δ)E是发电机暂态电抗后的电势V是无穷大母线电压X_Σ是包括发电机暂态电抗、变压器电抗、线路电抗的总电抗。这条正弦曲线就是静态稳定分析的“作战地图”。稳定运行的直观条件是功角特性曲线关键点处的斜率必须为正也就是dP_e/dδ大于0。这就像推一个放在斜坡上的球位移和回复力的方向相反扰动之后球会自动滚回原来的位置如果斜率变成负的球就会继续滚远系统就失去了同步能力。静态稳定极限就是sin(δ)1也就是δ90度对应的最大传输功率P_maxEV/X_Σ。1.2 电压稳定鼻子尖上的临界点除了功角稳定静态稳定性还要看节点电压能不能维持在一个可运行的水平。研究输电系统时通常让负荷功率逐步增长跟踪某个节点的电压幅值变化画出一条鼻形曲线。曲线的上半支是稳定运行区负荷增加节点电压下降但系统能保持平衡曲线的下半支对应的平衡点本质上不稳定。到了曲线的“鼻子尖”就是电压崩溃的临界点这个点对应的负荷功率就是静态电压稳定极限。电压静态稳定和功角静态稳定不是一回事但在实际系统中经常耦合。功角稳定极限看的是转子动能能不能耗散掉电压稳定极限看的是无功功率能不能支撑住节点电压。做仿真时如果把这两条主线分清楚后续选择模型和工具就会自然很多。1.3 把问题翻译成状态方程现代稳定性分析的标准方法是小信号线性化。把电力系统的微分方程在运行点附近做泰勒展开只保留一阶项就得到一组线性化状态方程dΔx/dt A * Δx状态矩阵A的特征值直接决定了系统在小扰动下的动态行为。特征值实部为负扰动会衰减系统静态稳定实部为正扰动会发散系统不稳定实部为零则处于临界状态。对于一台发电机加无穷大母线的简单系统状态变量一般取功角偏差Δδ和转速偏差Δω只有二阶手算都能解但一旦加入励磁系统、调速器、多台机组状态矩阵的阶数立刻上升到十几阶甚至几十阶这时候必须依靠Matlab脚本批量组装。这里我想强调一个很实用的习惯不要一上来就打开Simulink拖模块而是先把这个线性化状态矩阵写出来、把特征值算出来。Simulink是用来做时域验证和非线性观察的Matlab脚本才是做机理判断的主力。两种工具配合才能形成一个完整的“理论—仿真—回归验证”闭环。2. Matlab脚本先行单机无穷大系统的数学化2.1 系统参数和标幺值基准做仿真第一步是确定参数。以最经典的单机无穷大系统SMIB为例我采用一组很常见的标幺值参数参数数值说明E1.05 pu发电机暂态电势V∞1.00 pu无穷大母线电压x_d0.30 pu发电机暂态电抗x_T0.15 pu变压器电抗x_L0.35 pu等效线路电抗X_Σ0.80 pu总电抗H4.00 s发电机惯量时间常数f_s50 Hz系统额定频率标幺值是一个很容易翻车的点。发电机、变压器、线路可能有名值差别巨大但换算到同一个功率基准之后并联串联关系就变得极其简单整个系统就剩下加法。这里所有电抗统一折算到同一个基准容量下直接相加得到X_Σ0.8 pu这个做法在实际工程中也很常用。转子运动方程写成M * d²δ/dt² P_m - P_e - D_d * dδ/dt其中M 2H/ω_sω_s 2πf_s 314.159 rad/s。代入H4.0M ≈ 0.0255。这个数值很小意味着同样的功率不平衡会产生很大的角加速度这也是电力系统机电暂态过程有时候看起来很“猛”的原因。2.2 同步功率系数、临界功角和静态稳定极限以一机出线带负荷P₀0.8 pu为例。稳态运行时机械功率等于电磁功率所以初始功角由下式决定sin(δ₀) P₀ * X_Σ / (E * V∞)代入数值sin(δ₀) 0.8 * 0.8 / (1.05 * 1.0) 0.6095δ₀ ≈ 37.6度。同步功率系数K_s dP_e/dδ (EV∞/X_Σ) * cos(δ₀) ≈ 1.04 pu/rad。从静态稳定极限看P_max EV∞/X_Σ 1.3125 pu当前运行点0.8 pu的稳定裕度大约39%。这个裕度数字直接回答了“系统还能带多少负荷”的问题工程上非常直观。可以再把几个典型运行点的参数列出来对比运行功率 P₀ (pu)初始功角 (deg)K_s (pu/rad)稳定裕度0.5022.41.21461.9%0.8037.61.04039.0%1.2066.10.5328.6%从表中能明显看出运行点越接近静态稳定极限同步功率系数越小系统抵抗扰动的“弹性”越差。2.3 状态矩阵A与特征值直接用eig看稳定性用Matlab脚本把上述计算串起来代码非常短%% SMIB系统静态稳定特征值计算 clear; clc; % 系统参数标幺值 E 1.05; % 暂态电势 V 1.00; % 无穷大母线电压 X 0.80; % 总电抗 H 4.00; % 惯量时间常数 s fs 50; omega_s 2*pi*fs; M 2*H/omega_s; % 转子动量系数 P0 0.80; % 初始输出功率 delta0 asin(P0*X/(E*V)); % 初始功角 rad Ks (E*V/X) * cos(delta0); % 同步功率系数 Pmax E*V/X; fprintf(delta0 %.2f deg\n, rad2deg(delta0)); fprintf(Ks %.3f pu/rad\n, Ks); fprintf(Pmax %.3f pu\n, Pmax); % 线性化状态方程 ddx A * dx Dd 0.02; % 阻尼转矩系数注意这里针对的是转速偏差量 A [0 1; -Ks/M -Dd/M]; lam eig(A); fprintf(特征值 %.4f - j%.4f\n, real(lam(1)), abs(imag(lam(1))));这段代码跑出来的特征值大概是-0.39 ± j6.38 rad/s。实部为负说明当前运行点静态稳定虚部对应振荡角频率约6.38 rad/s折算频率约1.02 Hz这正是典型的电力系统低频振荡频率范围。阻尼比约6%属于弱阻尼但尚能收敛的状态和工程上许多弱互联系统的实际表现很像。如果改变运行点让P₀逐渐升高K_s逐渐变小特征值的实部会从负值向零移动越过零以后系统失去稳定。这就是静态稳定极限的数学含义。2.4 再加入励磁和调速状态方程怎么长大经典二阶模型把发电机电动势E当成常数这只用于最基础的机理分析。实际系统里有励磁调节器、PSS、调速器每个环节都带来额外的状态变量。比如加一个一阶励磁惯性环节状态变量从2个变成3个加上自动电压调节器AVR的PID环节再增加1到2个调速器是机械环节又增加1个。多机系统每台机组的自由度还要翻倍。用Matlab脚本组装状态矩阵时我建议按“模块”来组装每个环节写成一个子函数输入是参数和运行点输出是它对应的状态矩阵分块然后按发电机、励磁、调速的顺序拼装而不是把整个矩阵一次性手写。这种模块化组装方式好处是系统规模扩大时不容易出错也方便单独投退某个环节做对比。Simulink模型本质上也是在做同样的事只是把矩阵拼装过程用图形化接口包装了。3. Simulink模型搭建先用模块把原理“摆”出来3.1 简化二阶摇摆模型留足观察振荡的空间特征值计算只是纸上谈兵最终还是要用Simulink把时域响应跑出来给人看。我推荐先搭一个纯数学的简化二阶模型占满足够空间让振荡曲线看得很清楚。这个模型的原理就是转子运动方程。输入是机械功率P_m和电磁功率P_e电磁功率通过功角正弦函数计算得到P_e EV/X * sin(δ)。在Simulink里用Integrator模块对M*dω/dt积分得到转速偏差再积分一次得到功角。功角反馈给sin函数形成闭环。这个模型用到的模块只有Gain、Sum、Integrator、Trigonometric Function、Scope和Constant半小时内就能搭完。搭建的时候注意把Dd设计成一个工作区变量方便后续用脚本批量修改。还要注意功角的单位是弧度Scope里显示的是时间曲线但纵坐标读出的是功角数值不要和度搞混。如果想直接看度就在功角输出后面加一个Gain模块系数取180/pi。纯数学模型的优势是求解快、逻辑透明任何仿真结果都能立刻对照解析公式。这也是我做参数扫描的主要工具。3.2 从原理模型升级到Simscape Electrical元件级模型简化模型验证了机理之后如果项目要求更贴近实际就需要用Simscape Electrical的元件库搭一个“电气级”模型。核心思路是用Synchronous Machine pu Standard模块代表发电机用Three-Phase Source模块代表无穷大母线中间用变压器和传输线模块连接再用Three-Phase Fault模块做故障注入。具体操作时我一般这样组织模型电源侧无穷大母线用Three-Phase Source线电压220kV频率50Hz短路容量设为大值使其近似为恒压源。发电机选Synchronous Machine pu Standard填入第2章的标幺值参数。输电线路用Three-Phase PI Section Line设置正序、零序阻抗。测量在关注位置加Three-Phase V-I Measurement把电压电流采样总线引到Scope或To Workspace。电气求解模型里必须有一个Powergui模块这是Simscape Electrical模型正常初始化和求解的前提。元件级模型的搭建确实比纯数学模型繁琐但它能自然处理暂态和稳态的统一比如短路故障时电流电压的畸变、励磁绕组的电磁暂态都是简化模型无法表现的。这里我要特别提醒两个模型之间的结果对比必须保证初始运行点一致。很多人把简化模型调到P₀0.8 pu又在元件模型里随手给机械功率0.8 pu结果初始功角差了一大截仿真曲线完全对不上最后怀疑人生。正确做法是让Powergui先执行一次“Machines and Load Flow Initialization”自动把发电机初始角度算出来再去和特征值计算时的δ₀做对比两者吻合再往下继续。3.3 励磁系统、调速器的接入位置元件级模型里发电机模块的电气端口和机械端口都要接外部环节。机械功率输入口接原动机模型或者先用一个常数代表稳态机械功率励磁电压输入口接励磁系统和AVR可以用Simscape Electrical自带的Exciter库也可以手动搭一个简单的一阶惯性环节。我的建议是第一版模型先不接AVR直接把励磁电压设为常数这样模型最简便于和二阶模型的特征值结果对照。等确认基本动态行为无误后再接入励磁系统观察加入AVR后特征值是往左半平面移动还是往右半平面移动。实际中AVR通常会增强同步转矩但同时可能削弱阻尼转矩是引发低频振荡的重要因素这个现象值得专门做一组对比仿真来感受。3.4 求解器与步长刚性系统别硬扛纯数学二阶模型用变步长ode45就能跑因为微分方程比较温和。但一旦加入励磁绕组电磁暂态、电力电子器件系统就变成刚性系统时间常数跨了好几个数量级再强行用ode45仿真步长会被最小时间常数卡死一个20秒的仿真能跑出天荒地老的感觉。我的经验是Simscape Electrical模型优先选ode23t或者ode15s相对容差设为1e-3到1e-4这样精度和速度都能兼顾。另外Max Step Size建议设一个上限比如0.01秒防止大跨度跳过振荡峰值。纯数学模型可以大胆用ode45但在对比曲线时两种模型的采样频率最好保持一致否则后续数据处理时会出现一堆麻烦。4. 小扰动时域仿真与特征值对账4.1 扰动注入阶跃、短路、负荷突变三种手段静态稳定研究的是小扰动所以仿真扰动也要体现“小”字。三种常见手段机械功率阶跃在原动机输入上加一个阶跃幅值取0.01到0.05 pu这是最干净的小扰动。负荷突变在负荷节点加一个短时间的负荷脉冲或者用Step模块直接改变负荷功率。三相短路故障持续0.1秒后切除虽然是暂态故障但可用于观察故障切除后系统能否恢复到原运行点。我做小信号验证时最喜欢用机械功率阶跃因为它的物理含义和线性化方程里的ΔP_m完全对应。幅度取0.02 pu既不会激发明显的非线性又能从功角曲线上数出清晰的振荡周期和衰减。4.2 从功角曲线里挖出振荡频率和阻尼比时域仿真结束后最直接的输出是功角随时间的变化曲线。把数据导出到工作区后用findpeaks函数提取前几个峰值的幅值和时间可以算出振荡频率和阻尼比。在Matlab中一次典型的信号处理流程是这样的% 仿真数据已经存为 time, delta % 取稳态之后的振荡段 idx time 5; sig delta(idx) - mean(delta(idx)); % 用FFT看主振荡频率 Fs 1 / (time(2) - time(1)); nfft 2^nextpow2(length(sig)); freq Fs/2 * linspace(0, 1, nfft/21); Y fft(sig, nfft); plot(freq, 2*abs(Y(1:nfft/21))); xlim([0 5]);从频谱里能明显看到大约1 Hz处有一个峰这就是机电振荡模式。再用findpeaks找功角峰值序列通过相邻峰值的对数衰减率估算阻尼比将这个阻尼比和第2章特征值计算得到的理论阻尼比做对比两者误差在5%以内就是合格的。4.3 仿真结果和特征根对不上的常见原因如果时域仿真的振荡频率和特征值计算结果差很多不用急着怀疑工具先按顺序排查初始功角是否一致。这是最高频的问题纯数学模型手动设定δ₀元件模型靠Powergui初始化两边不统一动态响应自然对不上。系统电抗是否一致。元件模型里变压器和线路有各自的阻抗如果某个参数没填对总电抗就和脚本里的X₀不一致。阻尼模型差异。特征值计算里D_d是对转速偏差量的阻尼系数时域仿真里可能需要换算尤其是使用元件模型时阻尼可能由绕组电阻和阻尼绕组自然提供数值上和手动设置的D_d不一定能直接对应。扰动幅度太大。小扰动线性化成立的前提是偏差足够小如果扰动幅度超过0.1 pu非线性效应就会让振荡频率和阻尼比显著偏移这时候不该拿线性特征值去苛求非线性时域仿真。这些细节看起来琐碎但项目经验恰恰体现在这些“对不上”时的排查效率上。5. 静态电压稳定性PV曲线与极限负荷5.1 连续潮流一步步逼出鼻形点功角稳定看的是转子电压稳定看的是节点电压。工程上最直观的做法是画PV曲线让系统负荷按某个功率因数逐步增长然后追踪某个关键节点的电压幅值变化直到潮流方程无解。这个“无解”的时刻就是电压稳定极限。暴力做法是把负荷从0开始每增加0.01 pu就用fsolve解一次潮流初值总是取上一级负荷的解。这种做法其实就是连续潮流的最简形式区别在于专业连续潮流软件会加入步长控制和参数化策略防止在临界点附近发散。Matlab里用fsolve逐步逼近代码骨架如下clear; clc; V_s 1.00; % 源电压 pu X 0.50; % 等效电抗 pu电压稳定分析时单独取线路部分 pf 0.95; % 负荷功率因数 phi acos(pf); P_list 0:0.01:1.2; V_l NaN(size(P_list)); theta NaN(size(P_list)); for k 1:length(P_list) P P_list(k); Q P * tan(phi); F (x) [ (V_s*x(1)/X)*sin(x(2)) - P; (V_s*x(1)/X)*cos(x(2)) - x(1)^2/X - Q ]; if k 1 x0 [1; 0.02]; else x0 [V_l(k-1); theta(k-1)]; % 用上一步解当初值 end [x, fval] fsolve(F, x0, optimoptions(fsolve, Display, off)); if norm(fval) 1e-6 x(1) 0 V_l(k) x(1); theta(k) x(2); else V_l(k) NaN; theta(k) NaN; end end plot(P_list, V_l, LineWidth, 1.5); grid on; xlabel(负荷有功 P (pu)); ylabel(节点电压 V_L (pu));5.2 用fsolve解两节点潮流并画PV曲线上面的代码里x(1)是负荷节点电压幅值x(2)是负荷节点相对于源的相位差。两个潮流方程分别是节点有功和无功平衡。接近临界点时雅可比矩阵变得奇异fsolve会开始发散这正是物理上电压失稳在数学上的表现。运行结果画出来后会看到典型的“鼻形曲线”先是随负荷增加电压平缓下降接近极限时电压下降速度明显加快曲线上翘的痕迹开始显现最后曲线在某个点中断或者跳变这个断点就是数值临界点对应的P值就是静态电压稳定极限。5.3 电压稳定裕度和崩溃点判据电压稳定裕度定义为裕度 (P_cr - P_0) / P_cr * 100%假设当前运行点P₀0.6临界P_cr约1.05裕度大约43%。这个数值可以作为规划设计时衡量“电网还有多少富余输送能力”的量化指标。除了看曲线形状工程上常用灵敏度判据“dQ/dV”作为辅助指标但这个指标需要潮流计算输出雅可比矩阵对初学者稍微复杂。先用PV曲线的直观方法把问题看明白再进一步学习灵敏度指标路径会更顺利。6. 参数扫描阻尼和电抗对稳定极限的影响6.1 阻尼系数D_d的影响弱阻尼是怎么毁掉一个运行点的只做一组仿真说明不了问题参数扫描才是仿真工具的价值所在。把阻尼系数D_d从0逐步增加到0.05每个值都算一遍特征值再跑一遍时域仿真就能清楚看到阻尼对振荡衰减的决定性作用。在Matlab脚本里简单的循环就能实现Dlist 0:0.005:0.05; lambda_real zeros(size(Dlist)); lambda_imag zeros(size(Dlist)); for i 1:length(Dlist) Dd Dlist(i); A [0 1; -Ks/M -Dd/M]; lam eig(A); lambda_real(i) real(lam(1)); lambda_imag(i) abs(imag(lam(1))); end plot(Dlist, lambda_real, o-); grid on; xlabel(阻尼系数 D_d); ylabel(特征值实部);会看到随着D_d增大特征值实部从0向负方向线性移动而虚部基本不变。这说明阻尼系数主要影响振荡的衰减速率对振荡频率影响很小。时域曲线上的表现就是D_d0时功角等幅振荡永不收敛D_d0.02时衰减得很慢振荡十几个周期才平息D_d0.05时振荡几个周期就消失了。把这个对比图放进项目报告里比单纯堆公式有说服力得多。6.2 总电抗X的影响机组离电源越远越脆弱另一个值得扫描的参数是总电抗X_Σ。保持其他条件不变将X_Σ从0.4增大到1.0观察静态稳定极限P_max的变化。P_max EV / X_Σ这是一个反比例关系。X_Σ翻倍理论极限功率直接减半。这个结论在工程上的意义非常直接发电机经长线路送电、或者经过多级变压器升压降压电抗累积效应会明显削减退电能力。有些地方电网“窝电”表面上看是通道容量不够深层次原因就是联络阻抗过大导致静态稳定极限过低。6.3 用parsim批量跑仿真并对比曲线如果需要把不同参数下的时域曲线画在同一张图上Simulink的批量仿真功能parsim可以帮忙。思路是用Simulink.SimulationInput对象携带不同的变量赋值然后并行仿真。mdl smib_math; % 换成你的模型名 load_system(mdl); Dlist [0 0.01 0.02 0.05]; for i 1:length(Dlist) simIn(i) Simulink.SimulationInput(mdl); simIn(i) simIn(i).setVariable(Dd, Dlist(i)); end simOut parsim(simIn, ShowProgress, on);随后用simOut(1).delta等字段把数据提出来统一画图比较。parsim的好处是多个参数可以同时跑时间省一大截尤其在仿真时长30秒、有几十组参数时体验明显。不过要注意模型里的Dd变量必须事先定义为工作区变量否则setVariable不起作用。7. 仿真过程中踩过的几个坑7.1 Simulink报错powergui初始化失败先查初值元件级模型最常见的报错就是Powergui初始化失败。大多数情况下问题出在发电机初始功角或励磁电压偏离真实平衡点太远。处理办法是先用简化计算得到合理的初始功角和励磁电压再填到发电机模块的参数里或者选择让Powergui自动初始化并手动检查结果。不要试图用“硬凑”的方式消除初始化报错因为电力系统仿真对初值极其敏感一开始就错可能导致整个动态过程全部跑偏。我的习惯是把简化模型算出来的稳态值作为初值写入然后在模型调试时先看稳态量再叠加扰动分两步走。7.2 标幺值基准不统一的连环炸标幺值这东西小系统里你感觉不到它的好一旦出错就是连环炸。发电机电抗和变压器电抗的基准容量不一样如果直接相加X_Σ就错得离谱后面所有特征值、功角、稳定极限全是错的。正确做法是先把所有元件电抗折算到统一基准容量通常取发电机额定容量作为系统基准。Matlab里这一步可以写成函数集中处理避免在Simulink模块里反复手算。7.3 输出数据的“套娃”问题别再拿游标数格子了刚开始做仿真时我习惯直接在Scope图上用游标读数但后来发现项目里动辄几十组曲线游标数量根本不够用而且容易看错。更可靠的做法是给关键信号单独加To Workspace模块把数据存成结构体或数组用统一的后处理脚本完成峰值提取、频率分析和曲线导出。每次仿真完成后直接运行一个post_process脚本把所有关键指标一次算出来既省时间又不会数错格子。7.4 模型导入导出和C代码生成的注意事项最后提一点和工程落地相关的内容。Simulink模型经过代码生成之后可以部署到RCP快速原型控制器或者硬件在环设备上但电力系统模型通常不是实时仿真友好的。模型里如果有纯Simscape电气元件代码生成之前往往要做降阶处理。另外在版本管理时模型文件是二进制格式没法像文本代码那样轻易diff我的建议是定期给模型截图并建立一个简单的参数变更记录文档。否则改了个参数一个月后根本想不起来哪一版是哪一版。在我实际做这套分析的时候最有价值的习惯还是“先算一笔再跑一圈”。先在Matlab脚本里把特征值、稳态初值、稳定极限这些解析量全部算出来再带着期望值去Simulink里做验证而不是反过来让仿真替你找答案。这套工作方式换到其他电力系统分析问题里同样适用。 
前阵子帮一个电网规划项目算静态稳定极限被同事问了一句“能不能一小时给我出一张功角曲线图”。说实话很多人一听电力系统静态稳定性仿真脑子里全是大型电网模型、时域故障扫描、一堆PSS参数越想越复杂。但我真正把这个问题做顺手之后发现Matlab加Simulink这一套组合完全可以用一条清晰的主线串起来先把物理过程翻译成数学方程再用脚本把特征值解出来最后用Simulink时域仿真去验证和展示。这篇文章我把这条路径完整走一遍从建模思路、Matlab脚本到Simulink模型搭建、PV曲线求取再到参数扫描和踩坑经验适合正在做课程设计、毕业设计或者刚接触电力系统稳定分析的工程师参考。1. 静态稳定性分析到底算的是什么账1.1 功角稳定一个“两个转子”的故事电力系统静态稳定说白了就是看系统在某个运行点上受到一个小扰动之后能不能自己回到原来的稳态。这里的关键词是“小扰动”比如负荷波动零点几个百分点、线路阻抗的缓慢变化而不是雷击短路那种大扰动。大扰动对应的是暂态稳定小扰动对应的是静态稳定这是两套完全不同的分析思路很多人一开始就把这两个概念混在一起导致仿真设计跑偏。功角稳定的物理图像其实特别简单。一台同步发电机经线路接到无穷大母线可以把发电机等效成一个旋转的电动势源无穷大母线看成电压幅值和频率都恒定的源。两个源之间通过电抗相连它们之间有一个相位差δ也就是功角。发电机原动机输入的机械功率和送出去的电磁功率必须平衡而电磁功率随功角变化的规律就是功角特性P_e (E * V / X_Σ) * sin(δ)E是发电机暂态电抗后的电势V是无穷大母线电压X_Σ是包括发电机暂态电抗、变压器电抗、线路电抗的总电抗。这条正弦曲线就是静态稳定分析的“作战地图”。稳定运行的直观条件是功角特性曲线关键点处的斜率必须为正也就是dP_e/dδ大于0。这就像推一个放在斜坡上的球位移和回复力的方向相反扰动之后球会自动滚回原来的位置如果斜率变成负的球就会继续滚远系统就失去了同步能力。静态稳定极限就是sin(δ)1也就是δ90度对应的最大传输功率P_maxEV/X_Σ。1.2 电压稳定鼻子尖上的临界点除了功角稳定静态稳定性还要看节点电压能不能维持在一个可运行的水平。研究输电系统时通常让负荷功率逐步增长跟踪某个节点的电压幅值变化画出一条鼻形曲线。曲线的上半支是稳定运行区负荷增加节点电压下降但系统能保持平衡曲线的下半支对应的平衡点本质上不稳定。到了曲线的“鼻子尖”就是电压崩溃的临界点这个点对应的负荷功率就是静态电压稳定极限。电压静态稳定和功角静态稳定不是一回事但在实际系统中经常耦合。功角稳定极限看的是转子动能能不能耗散掉电压稳定极限看的是无功功率能不能支撑住节点电压。做仿真时如果把这两条主线分清楚后续选择模型和工具就会自然很多。1.3 把问题翻译成状态方程现代稳定性分析的标准方法是小信号线性化。把电力系统的微分方程在运行点附近做泰勒展开只保留一阶项就得到一组线性化状态方程dΔx/dt A * Δx状态矩阵A的特征值直接决定了系统在小扰动下的动态行为。特征值实部为负扰动会衰减系统静态稳定实部为正扰动会发散系统不稳定实部为零则处于临界状态。对于一台发电机加无穷大母线的简单系统状态变量一般取功角偏差Δδ和转速偏差Δω只有二阶手算都能解但一旦加入励磁系统、调速器、多台机组状态矩阵的阶数立刻上升到十几阶甚至几十阶这时候必须依靠Matlab脚本批量组装。这里我想强调一个很实用的习惯不要一上来就打开Simulink拖模块而是先把这个线性化状态矩阵写出来、把特征值算出来。Simulink是用来做时域验证和非线性观察的Matlab脚本才是做机理判断的主力。两种工具配合才能形成一个完整的“理论—仿真—回归验证”闭环。2. Matlab脚本先行单机无穷大系统的数学化2.1 系统参数和标幺值基准做仿真第一步是确定参数。以最经典的单机无穷大系统SMIB为例我采用一组很常见的标幺值参数参数数值说明E1.05 pu发电机暂态电势V∞1.00 pu无穷大母线电压x_d0.30 pu发电机暂态电抗x_T0.15 pu变压器电抗x_L0.35 pu等效线路电抗X_Σ0.80 pu总电抗H4.00 s发电机惯量时间常数f_s50 Hz系统额定频率标幺值是一个很容易翻车的点。发电机、变压器、线路可能有名值差别巨大但换算到同一个功率基准之后并联串联关系就变得极其简单整个系统就剩下加法。这里所有电抗统一折算到同一个基准容量下直接相加得到X_Σ0.8 pu这个做法在实际工程中也很常用。转子运动方程写成M * d²δ/dt² P_m - P_e - D_d * dδ/dt其中M 2H/ω_sω_s 2πf_s 314.159 rad/s。代入H4.0M ≈ 0.0255。这个数值很小意味着同样的功率不平衡会产生很大的角加速度这也是电力系统机电暂态过程有时候看起来很“猛”的原因。2.2 同步功率系数、临界功角和静态稳定极限以一机出线带负荷P₀0.8 pu为例。稳态运行时机械功率等于电磁功率所以初始功角由下式决定sin(δ₀) P₀ * X_Σ / (E * V∞)代入数值sin(δ₀) 0.8 * 0.8 / (1.05 * 1.0) 0.6095δ₀ ≈ 37.6度。同步功率系数K_s dP_e/dδ (EV∞/X_Σ) * cos(δ₀) ≈ 1.04 pu/rad。从静态稳定极限看P_max EV∞/X_Σ 1.3125 pu当前运行点0.8 pu的稳定裕度大约39%。这个裕度数字直接回答了“系统还能带多少负荷”的问题工程上非常直观。可以再把几个典型运行点的参数列出来对比运行功率 P₀ (pu)初始功角 (deg)K_s (pu/rad)稳定裕度0.5022.41.21461.9%0.8037.61.04039.0%1.2066.10.5328.6%从表中能明显看出运行点越接近静态稳定极限同步功率系数越小系统抵抗扰动的“弹性”越差。2.3 状态矩阵A与特征值直接用eig看稳定性用Matlab脚本把上述计算串起来代码非常短%% SMIB系统静态稳定特征值计算 clear; clc; % 系统参数标幺值 E 1.05; % 暂态电势 V 1.00; % 无穷大母线电压 X 0.80; % 总电抗 H 4.00; % 惯量时间常数 s fs 50; omega_s 2*pi*fs; M 2*H/omega_s; % 转子动量系数 P0 0.80; % 初始输出功率 delta0 asin(P0*X/(E*V)); % 初始功角 rad Ks (E*V/X) * cos(delta0); % 同步功率系数 Pmax E*V/X; fprintf(delta0 %.2f deg\n, rad2deg(delta0)); fprintf(Ks %.3f pu/rad\n, Ks); fprintf(Pmax %.3f pu\n, Pmax); % 线性化状态方程 ddx A * dx Dd 0.02; % 阻尼转矩系数注意这里针对的是转速偏差量 A [0 1; -Ks/M -Dd/M]; lam eig(A); fprintf(特征值 %.4f - j%.4f\n, real(lam(1)), abs(imag(lam(1))));这段代码跑出来的特征值大概是-0.39 ± j6.38 rad/s。实部为负说明当前运行点静态稳定虚部对应振荡角频率约6.38 rad/s折算频率约1.02 Hz这正是典型的电力系统低频振荡频率范围。阻尼比约6%属于弱阻尼但尚能收敛的状态和工程上许多弱互联系统的实际表现很像。如果改变运行点让P₀逐渐升高K_s逐渐变小特征值的实部会从负值向零移动越过零以后系统失去稳定。这就是静态稳定极限的数学含义。2.4 再加入励磁和调速状态方程怎么长大经典二阶模型把发电机电动势E当成常数这只用于最基础的机理分析。实际系统里有励磁调节器、PSS、调速器每个环节都带来额外的状态变量。比如加一个一阶励磁惯性环节状态变量从2个变成3个加上自动电压调节器AVR的PID环节再增加1到2个调速器是机械环节又增加1个。多机系统每台机组的自由度还要翻倍。用Matlab脚本组装状态矩阵时我建议按“模块”来组装每个环节写成一个子函数输入是参数和运行点输出是它对应的状态矩阵分块然后按发电机、励磁、调速的顺序拼装而不是把整个矩阵一次性手写。这种模块化组装方式好处是系统规模扩大时不容易出错也方便单独投退某个环节做对比。Simulink模型本质上也是在做同样的事只是把矩阵拼装过程用图形化接口包装了。3. Simulink模型搭建先用模块把原理“摆”出来3.1 简化二阶摇摆模型留足观察振荡的空间特征值计算只是纸上谈兵最终还是要用Simulink把时域响应跑出来给人看。我推荐先搭一个纯数学的简化二阶模型占满足够空间让振荡曲线看得很清楚。这个模型的原理就是转子运动方程。输入是机械功率P_m和电磁功率P_e电磁功率通过功角正弦函数计算得到P_e EV/X * sin(δ)。在Simulink里用Integrator模块对M*dω/dt积分得到转速偏差再积分一次得到功角。功角反馈给sin函数形成闭环。这个模型用到的模块只有Gain、Sum、Integrator、Trigonometric Function、Scope和Constant半小时内就能搭完。搭建的时候注意把Dd设计成一个工作区变量方便后续用脚本批量修改。还要注意功角的单位是弧度Scope里显示的是时间曲线但纵坐标读出的是功角数值不要和度搞混。如果想直接看度就在功角输出后面加一个Gain模块系数取180/pi。纯数学模型的优势是求解快、逻辑透明任何仿真结果都能立刻对照解析公式。这也是我做参数扫描的主要工具。3.2 从原理模型升级到Simscape Electrical元件级模型简化模型验证了机理之后如果项目要求更贴近实际就需要用Simscape Electrical的元件库搭一个“电气级”模型。核心思路是用Synchronous Machine pu Standard模块代表发电机用Three-Phase Source模块代表无穷大母线中间用变压器和传输线模块连接再用Three-Phase Fault模块做故障注入。具体操作时我一般这样组织模型电源侧无穷大母线用Three-Phase Source线电压220kV频率50Hz短路容量设为大值使其近似为恒压源。发电机选Synchronous Machine pu Standard填入第2章的标幺值参数。输电线路用Three-Phase PI Section Line设置正序、零序阻抗。测量在关注位置加Three-Phase V-I Measurement把电压电流采样总线引到Scope或To Workspace。电气求解模型里必须有一个Powergui模块这是Simscape Electrical模型正常初始化和求解的前提。元件级模型的搭建确实比纯数学模型繁琐但它能自然处理暂态和稳态的统一比如短路故障时电流电压的畸变、励磁绕组的电磁暂态都是简化模型无法表现的。这里我要特别提醒两个模型之间的结果对比必须保证初始运行点一致。很多人把简化模型调到P₀0.8 pu又在元件模型里随手给机械功率0.8 pu结果初始功角差了一大截仿真曲线完全对不上最后怀疑人生。正确做法是让Powergui先执行一次“Machines and Load Flow Initialization”自动把发电机初始角度算出来再去和特征值计算时的δ₀做对比两者吻合再往下继续。3.3 励磁系统、调速器的接入位置元件级模型里发电机模块的电气端口和机械端口都要接外部环节。机械功率输入口接原动机模型或者先用一个常数代表稳态机械功率励磁电压输入口接励磁系统和AVR可以用Simscape Electrical自带的Exciter库也可以手动搭一个简单的一阶惯性环节。我的建议是第一版模型先不接AVR直接把励磁电压设为常数这样模型最简便于和二阶模型的特征值结果对照。等确认基本动态行为无误后再接入励磁系统观察加入AVR后特征值是往左半平面移动还是往右半平面移动。实际中AVR通常会增强同步转矩但同时可能削弱阻尼转矩是引发低频振荡的重要因素这个现象值得专门做一组对比仿真来感受。3.4 求解器与步长刚性系统别硬扛纯数学二阶模型用变步长ode45就能跑因为微分方程比较温和。但一旦加入励磁绕组电磁暂态、电力电子器件系统就变成刚性系统时间常数跨了好几个数量级再强行用ode45仿真步长会被最小时间常数卡死一个20秒的仿真能跑出天荒地老的感觉。我的经验是Simscape Electrical模型优先选ode23t或者ode15s相对容差设为1e-3到1e-4这样精度和速度都能兼顾。另外Max Step Size建议设一个上限比如0.01秒防止大跨度跳过振荡峰值。纯数学模型可以大胆用ode45但在对比曲线时两种模型的采样频率最好保持一致否则后续数据处理时会出现一堆麻烦。4. 小扰动时域仿真与特征值对账4.1 扰动注入阶跃、短路、负荷突变三种手段静态稳定研究的是小扰动所以仿真扰动也要体现“小”字。三种常见手段机械功率阶跃在原动机输入上加一个阶跃幅值取0.01到0.05 pu这是最干净的小扰动。负荷突变在负荷节点加一个短时间的负荷脉冲或者用Step模块直接改变负荷功率。三相短路故障持续0.1秒后切除虽然是暂态故障但可用于观察故障切除后系统能否恢复到原运行点。我做小信号验证时最喜欢用机械功率阶跃因为它的物理含义和线性化方程里的ΔP_m完全对应。幅度取0.02 pu既不会激发明显的非线性又能从功角曲线上数出清晰的振荡周期和衰减。4.2 从功角曲线里挖出振荡频率和阻尼比时域仿真结束后最直接的输出是功角随时间的变化曲线。把数据导出到工作区后用findpeaks函数提取前几个峰值的幅值和时间可以算出振荡频率和阻尼比。在Matlab中一次典型的信号处理流程是这样的% 仿真数据已经存为 time, delta % 取稳态之后的振荡段 idx time 5; sig delta(idx) - mean(delta(idx)); % 用FFT看主振荡频率 Fs 1 / (time(2) - time(1)); nfft 2^nextpow2(length(sig)); freq Fs/2 * linspace(0, 1, nfft/21); Y fft(sig, nfft); plot(freq, 2*abs(Y(1:nfft/21))); xlim([0 5]);从频谱里能明显看到大约1 Hz处有一个峰这就是机电振荡模式。再用findpeaks找功角峰值序列通过相邻峰值的对数衰减率估算阻尼比将这个阻尼比和第2章特征值计算得到的理论阻尼比做对比两者误差在5%以内就是合格的。4.3 仿真结果和特征根对不上的常见原因如果时域仿真的振荡频率和特征值计算结果差很多不用急着怀疑工具先按顺序排查初始功角是否一致。这是最高频的问题纯数学模型手动设定δ₀元件模型靠Powergui初始化两边不统一动态响应自然对不上。系统电抗是否一致。元件模型里变压器和线路有各自的阻抗如果某个参数没填对总电抗就和脚本里的X₀不一致。阻尼模型差异。特征值计算里D_d是对转速偏差量的阻尼系数时域仿真里可能需要换算尤其是使用元件模型时阻尼可能由绕组电阻和阻尼绕组自然提供数值上和手动设置的D_d不一定能直接对应。扰动幅度太大。小扰动线性化成立的前提是偏差足够小如果扰动幅度超过0.1 pu非线性效应就会让振荡频率和阻尼比显著偏移这时候不该拿线性特征值去苛求非线性时域仿真。这些细节看起来琐碎但项目经验恰恰体现在这些“对不上”时的排查效率上。5. 静态电压稳定性PV曲线与极限负荷5.1 连续潮流一步步逼出鼻形点功角稳定看的是转子电压稳定看的是节点电压。工程上最直观的做法是画PV曲线让系统负荷按某个功率因数逐步增长然后追踪某个关键节点的电压幅值变化直到潮流方程无解。这个“无解”的时刻就是电压稳定极限。暴力做法是把负荷从0开始每增加0.01 pu就用fsolve解一次潮流初值总是取上一级负荷的解。这种做法其实就是连续潮流的最简形式区别在于专业连续潮流软件会加入步长控制和参数化策略防止在临界点附近发散。Matlab里用fsolve逐步逼近代码骨架如下clear; clc; V_s 1.00; % 源电压 pu X 0.50; % 等效电抗 pu电压稳定分析时单独取线路部分 pf 0.95; % 负荷功率因数 phi acos(pf); P_list 0:0.01:1.2; V_l NaN(size(P_list)); theta NaN(size(P_list)); for k 1:length(P_list) P P_list(k); Q P * tan(phi); F (x) [ (V_s*x(1)/X)*sin(x(2)) - P; (V_s*x(1)/X)*cos(x(2)) - x(1)^2/X - Q ]; if k 1 x0 [1; 0.02]; else x0 [V_l(k-1); theta(k-1)]; % 用上一步解当初值 end [x, fval] fsolve(F, x0, optimoptions(fsolve, Display, off)); if norm(fval) 1e-6 x(1) 0 V_l(k) x(1); theta(k) x(2); else V_l(k) NaN; theta(k) NaN; end end plot(P_list, V_l, LineWidth, 1.5); grid on; xlabel(负荷有功 P (pu)); ylabel(节点电压 V_L (pu));5.2 用fsolve解两节点潮流并画PV曲线上面的代码里x(1)是负荷节点电压幅值x(2)是负荷节点相对于源的相位差。两个潮流方程分别是节点有功和无功平衡。接近临界点时雅可比矩阵变得奇异fsolve会开始发散这正是物理上电压失稳在数学上的表现。运行结果画出来后会看到典型的“鼻形曲线”先是随负荷增加电压平缓下降接近极限时电压下降速度明显加快曲线上翘的痕迹开始显现最后曲线在某个点中断或者跳变这个断点就是数值临界点对应的P值就是静态电压稳定极限。5.3 电压稳定裕度和崩溃点判据电压稳定裕度定义为裕度 (P_cr - P_0) / P_cr * 100%假设当前运行点P₀0.6临界P_cr约1.05裕度大约43%。这个数值可以作为规划设计时衡量“电网还有多少富余输送能力”的量化指标。除了看曲线形状工程上常用灵敏度判据“dQ/dV”作为辅助指标但这个指标需要潮流计算输出雅可比矩阵对初学者稍微复杂。先用PV曲线的直观方法把问题看明白再进一步学习灵敏度指标路径会更顺利。6. 参数扫描阻尼和电抗对稳定极限的影响6.1 阻尼系数D_d的影响弱阻尼是怎么毁掉一个运行点的只做一组仿真说明不了问题参数扫描才是仿真工具的价值所在。把阻尼系数D_d从0逐步增加到0.05每个值都算一遍特征值再跑一遍时域仿真就能清楚看到阻尼对振荡衰减的决定性作用。在Matlab脚本里简单的循环就能实现Dlist 0:0.005:0.05; lambda_real zeros(size(Dlist)); lambda_imag zeros(size(Dlist)); for i 1:length(Dlist) Dd Dlist(i); A [0 1; -Ks/M -Dd/M]; lam eig(A); lambda_real(i) real(lam(1)); lambda_imag(i) abs(imag(lam(1))); end plot(Dlist, lambda_real, o-); grid on; xlabel(阻尼系数 D_d); ylabel(特征值实部);会看到随着D_d增大特征值实部从0向负方向线性移动而虚部基本不变。这说明阻尼系数主要影响振荡的衰减速率对振荡频率影响很小。时域曲线上的表现就是D_d0时功角等幅振荡永不收敛D_d0.02时衰减得很慢振荡十几个周期才平息D_d0.05时振荡几个周期就消失了。把这个对比图放进项目报告里比单纯堆公式有说服力得多。6.2 总电抗X的影响机组离电源越远越脆弱另一个值得扫描的参数是总电抗X_Σ。保持其他条件不变将X_Σ从0.4增大到1.0观察静态稳定极限P_max的变化。P_max EV / X_Σ这是一个反比例关系。X_Σ翻倍理论极限功率直接减半。这个结论在工程上的意义非常直接发电机经长线路送电、或者经过多级变压器升压降压电抗累积效应会明显削减退电能力。有些地方电网“窝电”表面上看是通道容量不够深层次原因就是联络阻抗过大导致静态稳定极限过低。6.3 用parsim批量跑仿真并对比曲线如果需要把不同参数下的时域曲线画在同一张图上Simulink的批量仿真功能parsim可以帮忙。思路是用Simulink.SimulationInput对象携带不同的变量赋值然后并行仿真。mdl smib_math; % 换成你的模型名 load_system(mdl); Dlist [0 0.01 0.02 0.05]; for i 1:length(Dlist) simIn(i) Simulink.SimulationInput(mdl); simIn(i) simIn(i).setVariable(Dd, Dlist(i)); end simOut parsim(simIn, ShowProgress, on);随后用simOut(1).delta等字段把数据提出来统一画图比较。parsim的好处是多个参数可以同时跑时间省一大截尤其在仿真时长30秒、有几十组参数时体验明显。不过要注意模型里的Dd变量必须事先定义为工作区变量否则setVariable不起作用。7. 仿真过程中踩过的几个坑7.1 Simulink报错powergui初始化失败先查初值元件级模型最常见的报错就是Powergui初始化失败。大多数情况下问题出在发电机初始功角或励磁电压偏离真实平衡点太远。处理办法是先用简化计算得到合理的初始功角和励磁电压再填到发电机模块的参数里或者选择让Powergui自动初始化并手动检查结果。不要试图用“硬凑”的方式消除初始化报错因为电力系统仿真对初值极其敏感一开始就错可能导致整个动态过程全部跑偏。我的习惯是把简化模型算出来的稳态值作为初值写入然后在模型调试时先看稳态量再叠加扰动分两步走。7.2 标幺值基准不统一的连环炸标幺值这东西小系统里你感觉不到它的好一旦出错就是连环炸。发电机电抗和变压器电抗的基准容量不一样如果直接相加X_Σ就错得离谱后面所有特征值、功角、稳定极限全是错的。正确做法是先把所有元件电抗折算到统一基准容量通常取发电机额定容量作为系统基准。Matlab里这一步可以写成函数集中处理避免在Simulink模块里反复手算。7.3 输出数据的“套娃”问题别再拿游标数格子了刚开始做仿真时我习惯直接在Scope图上用游标读数但后来发现项目里动辄几十组曲线游标数量根本不够用而且容易看错。更可靠的做法是给关键信号单独加To Workspace模块把数据存成结构体或数组用统一的后处理脚本完成峰值提取、频率分析和曲线导出。每次仿真完成后直接运行一个post_process脚本把所有关键指标一次算出来既省时间又不会数错格子。7.4 模型导入导出和C代码生成的注意事项最后提一点和工程落地相关的内容。Simulink模型经过代码生成之后可以部署到RCP快速原型控制器或者硬件在环设备上但电力系统模型通常不是实时仿真友好的。模型里如果有纯Simscape电气元件代码生成之前往往要做降阶处理。另外在版本管理时模型文件是二进制格式没法像文本代码那样轻易diff我的建议是定期给模型截图并建立一个简单的参数变更记录文档。否则改了个参数一个月后根本想不起来哪一版是哪一版。在我实际做这套分析的时候最有价值的习惯还是“先算一笔再跑一圈”。先在Matlab脚本里把特征值、稳态初值、稳定极限这些解析量全部算出来再带着期望值去Simulink里做验证而不是反过来让仿真替你找答案。这套工作方式换到其他电力系统分析问题里同样适用。