单机无穷大系统暂态稳定仿真：SMIG_2.m脚本全解析

简介：单机无穷大系统是电力系统暂态稳定分析中的经典简化模型，这份资源面向电力系统专业学生、科研人员及MATLAB仿真学习者，提供该模型的脚本化仿真实现。压缩包内仅含1个m文件，大小约2KB，对应完整MATLAB源码，不依赖Simulink，可直接运行并深入阅读算法细节。代码涵盖发电机动态建模、励磁控制器逻辑、系统运动方程、初始条件设置、数值积分步进及结果绘图等多个模块，适合用于研究负荷突变、短路扰动等场景下的响应特性。资源已有263人学习，小巧精悍，特别适合希望从底层理解单机无穷大系统仿真流程，并在此基础上扩展自定义控制策略的读者。通过研读和运行这份代码，可直观掌握电力系统动态仿真的核心步骤，为后续开展更复杂的多机系统分析打下基础。

1. 单机无穷大系统仿真：从暂态稳定问题到SMIG_2.m脚本落地

无限大母线（infinite bus）这个词在电力系统教材里出现频率极高，但真正把它做成可复现代码的公开资源并不多。单机无穷大模型把电网等值成电压幅值稳定、频率恒定、容量无穷大的理想节点，研究焦点因此被压缩到发电机本体和它到母线之间的等值电抗上。无论是分析短路故障后的功角摇摆，还是校验励磁控制器（AVR）、电力系统稳定器（PSS）的参数，这个模型都比多机系统更干净、更易于定位问题，是暂态稳定分析的入门标配。

SMIG_2.rar 中的 SMIG_2.m 用纯 MATLAB 脚本（而不是 Simulink）实现了单机无穷大系统的完整仿真流程。脚本方式最大的好处在于每一行代码对应一个数学方程，从同步电机电压方程到励磁控制、再到数值积分全程透明，改参数、扫工况、做批量校核都很方便。适合电力系统分析与控制方向的学生做课程设计复现，也适合工程师在保护整定和稳定性预研时做快速验证。

2. 同步电机与无穷大母线建模：SMIG_2.m的动态方程内核

2.1 无穷大母线假设的边界在哪里

单机无穷大系统仿真的第一个关键动作不是写代码，而是确认模型假设能否成立。无穷大母线要求在研究的全部时间内满足三个条件：母线电压幅值恒定、频率恒定（相位以额定同步转速旋转）、等效短路容量远大于所接发电机的容量。第三点意味着当这台发电机投切或发生扰动时，不会对母线电压和频率造成可观测的影响。现实电网中，如果研究对象的容量占系统总容量比例低于3%到5%，并且选在强联系站点（短路比SCR大于等于3），这个假设通常可以被接受。

在某些场景下，这个假设会被滥用。比如分析次同步振荡（SSR）时，发电机与串补线路之间的机电扭振相互作用已经不能用等值电抗描述，无穷大母线模型会漏掉关键的谐振频段。再比如研究区域间低频振荡时，关心的正是多台发电机通过有限容量电网相互牵制的行为，此时把系统等值成无穷大母线等于直接删掉了要研究的问题。理解这个边界，才能正确解释SMIG_2.m仿真结果的有效范围。

2.2 派克变换后的电压与磁链方程

同步电机的三相定子方程在静止abc坐标下是时变系数的微分方程组，直接数值求解效率很低。派克变换把定子量投影到随转子旋转的dq0坐标之后，电感矩阵变成恒定值，这是所有同步电机数字仿真得以进行的前提。

在dq0坐标下，忽略定子磁链微分项（即不计定子电磁暂态，只保留基波分量），电压方程和磁链方程可以写成如下形式：

% 电枢电压方程（标幺值，忽略定子暂态） % vd = -rs*id + xq*iq % vq = -rs*iq - xd*id + eqp % 其中 eqp 为 q 轴暂态电动势

这里vd、vq是机端电压的dq轴分量，id、iq是定子电流的dq轴分量，rs是电枢电阻，xd是d轴电抗（在三阶模型中对应暂态电抗），xq是q轴同步电抗，eqp来自励磁绕组动态。省略定子磁链微分项，相当于假设定子磁链能瞬时跟随转子运动变化，这一近似把仿真步长从微秒级放宽到毫秒级，同时依然能正确反映机电暂态的主要特征，这也是绝大多数暂态稳定程序的做法。

2.3 转子运动方程的标幺化与离散化

发电机的转子运动方程，通常称为摇摆方程（swing equation），是暂态稳定分析的核心。标幺化之后的形式是：

% d(delta)/dt = omega_b * (omega - 1) % d(omega)/dt = (1/(2*H)) * (Pm - Pe - D*(omega - 1))

delta是功角（转子q轴与无穷大母线电压相量之间的夹角），omega是转子角速度的标幺值，omega_b是额定电角速度（50Hz系统对应100*pi rad/s约等于314.159 rad/s），H是机组惯性时间常数（单位秒），Pm和Pe分别是机械功率和电磁功率的标幺值，D是阻尼系数。第一个方程把角度变化与速度偏差联系起来，第二个方程描述转子动能的变化率等于加速功率。

一个常见的数值陷阱是阻尼系数D的量纲。在标幺值模型中，D通常定义为额定转矩基准下单位速度偏差对应的阻尼转矩，取值一般在1到5之间。但部分教材把D放在摇摆方程的另一个位置，写作d(omega)/dt = (Pm - Pe - D*omega)/(2H)，两者含义完全不同，代码直接照搬时很容易把阻尼作用放大或缩小一个数量级。看到仿真曲线上功角衰减过快或等幅振荡不止时，第一反应应该是核对D的写法而不是去调积分步长。

2.4 SMIG_2.m动态方程函数的代码骨架

把上述方程组合起来，SMIG_2.m中典型的动态模型函数可以写成：

function dx = smib_dynamics(t, x, u, p) % 单机无穷大系统三阶动态模型 % 状态向量 x = [delta; omega; Eqp] delta = x(1); omega = x(2); Eqp = x(3); % 电气回路等值电抗（发电机暂态电抗 + 外接电抗） XdS = p.Xdp + p.Xe; XqS = p.Xq + p.Xe; % 由网络代数方程解 dq 轴电流 Id = (Eqp - p.Vb * cos(delta)) / XdS; Iq = p.Vb * sin(delta) / XqS; % 电磁功率（含凸极效应项） Pe = Eqp * Iq + (XqS - XdS) * Id * Iq; % 转子运动方程 d_delta = p.omega_b * (omega - 1); d_omega = (p.Pm - Pe - p.D * (omega - 1)) / (2 * p.H); % 励磁绕组暂态方程，u.Efd 为励磁电动势（由控制器决定） d_Eqp = (u.Efd - Eqp - (p.Xd - p.Xdp) * Id) / p.Td0p; dx = [d_delta; d_omega; d_Eqp]; end

代码里的参数结构体p通常包含：无穷大母线电压幅值Vb、外接电抗（线路与变压器等值）Xe、发电机d轴同步电抗Xd、d轴暂态电抗Xdp、q轴同步电抗Xq、惯性时间常数H、阻尼系数D、d轴开路暂态时间常数Td0p、额定角速度omega_b。这个三阶模型保留了励磁绕组动态，因此可以自然接入励磁控制器和故障扰动分析，这是SMIG_2.m能支持闭环控制仿真的基础。

构造这段代码有两个关键点。第一，Id、Iq的求解除派克变换外还引入了外接电抗Xe，实际计算时需要注意基准容量必须与发电机额定容量一致。第二，Pe表达式中的凸极项(XqS - XdS)IdIq在XdS与XqS相等时自动消失，退化为Pe = EqpVbsin(delta)/XdS的经典模型；如果不需要考虑凸极效应，直接令Xq等于Xdp，既能少一个参数，也能减少一个潜在的错误源。

3. 励磁控制器与故障扰动注入：SMIG_2.m中的闭环仿真实现

3.1 励磁系统模型选择：从恒定励磁到AVR比例控制

三阶模型里的u.Efd是励磁电动势输入，最简单的做法是让Efd恒定，即发电机在故障期间保持恒定励磁，此时电磁暂态模型退化为二阶经典模型，可以用来研究失步的基本形态。但要贴近实际机组行为，就必须加入励磁调节器（AVR）模型。

常见的可控硅励磁方式可以简化为一阶惯性环节加限幅：

% 励磁系统一阶惯性模型 % d(Efd)/dt = (Ka * (Vref - Vt + Vpss) - Efd) / Ta

Ka是励磁增益，一般在100到400之间；Ta是励磁机时间常数，典型值为0.02到0.1秒；Vt是机端电压幅值，Vref是电压参考值，Vpss是PSS附加信号。限幅环节把Efd约束在Efd_min和Efd_max之间，这两个限幅值直接决定了故障后的强励能力，是电压恢复过程的关键参数。工程上发电机强励顶值为2倍额定励磁电压左右，如果仿真中电压迟迟不恢复，先检查限幅上限是否设置过低。

实际代码中，励磁系统多以独立函数或控制器结构体的形式实现，而不是混在主方程函数里。这样做的用意是，后续要对比恒定励磁、AVR、AVR加PSS三种策略时，只需要替换控制器回调函数，不改动机电动态核心代码。

3.2 故障场景注入：短路与负荷突变的时间事件处理

单机无穷大系统仿真最典型的扰动是机端附近的三相短路故障，故障在t_sw时刻发生，在t_cl时刻被切除。MATLAB脚本中，故障导致网络参数变化是用时间事件切换实现的，不需要把方程组重写一遍：

for step = 1:Nt t_now = (step - 1) * h; % 判断当前仿真时刻处于哪个阶段 if t_now < t_sw Xe_on = Xe_normal; elseif t_now < t_cl Xe_on = Xe_fault; % 故障期间母线电压被短路拉低 else Xe_on = Xe_normal; % 故障切除后系统拓扑恢复 end % 更新网络参数并做一步积分 p.Xe = Xe_on; x_next = rk4_step(@smib_dynamics, t_now, x_now, h, u_ctrl, p); x_now = x_next; end

故障期间外接电抗Xe_fault的取值是最容易出错的地方。如果是发电机机端三相短路，等效地认为母线电压被拉低到接近零，Xe_fault取一个很小的值（如0.001标幺值）就能模拟电气距离被短路旁路的效果；如果是输电线路某点故障，则要重新计算故障点到发电机之间的等值电抗。注意不能把Xe_fault直接设为0，否则电磁功率表达式除零，仿真直接发散。

3.3 闭环控制器与主仿真循环的整合

把励磁控制器和故障逻辑整合进主循环，就得到了SMIG_2.m的整体结构。实际实现中，每一步先根据当前机端电压Vt计算AVR输出的Efd，再做一次RK4步进：

% 主仿真循环：求解机端电压 -> AVR -> RK4 步进 for step = 1:Nt t_now = (step - 1) * h; % 网络参数切换（故障/正常） p.Xe = get_network_impedance(t_now, t_sw, t_cl, Xe_normal, Xe_fault); % 由当前状态求解 dq 轴电流和机端电压 XdS = p.Xdp + p.Xe; XqS = p.Xq + p.Xe; Id = (x_now(3) - p.Vb * cos(x_now(1))) / XdS; Iq = p.Vb * sin(x_now(1)) / XqS; % 机端电压幅值，q 轴为实轴、d 轴为虚轴的相量表示 Vt = sqrt((x_now(3) - p.Xdp * Id)^2 + (p.Xq * Iq)^2); % AVR 比例控制加限幅 Efd = max(Efd_min, min(Efd_max, Ka * (Vref - Vt))); u_ctrl.Efd = Efd; x_next = rk4_step(@smib_dynamics, t_now, x_now, h, u_ctrl, p); x_now = x_next; record(t_now, x_now, Vt, Efd); end

值得注意，饱和函数必须在RK4子步处理之外执行，否则控制器输出在每个子步内被重复限幅，尤其步长较大时会产生比实际励磁系统更多的非线性畸变。另一种做法是把限幅逻辑放进smib_dynamics内部，但那样相当于在每个子步反复计算限幅，语义不一致。建议在循环层调用saturate函数并保存未限幅的原始AVR输出，方便后处理时区分是电压调节作用还是励磁限幅起了作用。

3.4 扰动场景参数设置参考

负荷突增场景需要留意，在单机无穷大模型里负荷突变等价于发电机输出电功率变化，因为无穷大母线本身不会发生功率不平衡；典型实现会直接修改Pm或者修改外接电抗来改变发电机的输出功率，效果是等价的。

4. 数值积分与仿真参数整定：从欧拉法到四阶龙格-库塔迭代实现

4.1 稳态初始条件求解：先找潮流解，再启动动态仿真

一个常见错误是直接把功角初值设为0或随意设定转速初值，然后开始积分。单机无穷大系统的动态仿真要求从稳态运行点启动，否则一开始就会产生一段人为的过渡过程，污染故障后的响应波形。

标准做法是先指定发电机出口的有功功率Pg0和无功功率Qg0（或功率因数），反解出初始功角delta0、暂态电动势Eqp0和励磁电动势Efd0。对三阶模型，典型实现如下：

% 由 Pg0、Qg0 求初始功角和电动势 % 假设机端电压 Vt0 给定，电流相量 I0 = conj(S0 / Vt0) S0 = Pg0 + 1j * Qg0; I0 = conj(S0 / Vt0); Ep = Vt0 + 1j * p.Xdp * I0; % 暂态电动势相量 delta0 = angle(Ep); Eqp0 = abs(Ep); % 稳态励磁电动势：令 d_Eqp/dt = 0 反推 Id0 = real(I0 * exp(1j * (pi/2 - delta0))); % d 轴电流，相位参考与派克变换一致 Efd0 = Eqp0 + (p.Xd - p.Xdp) * Id0;

这里的Id0应通过派克反变换从相量电流中取得，不同教材的dq变换矩阵相差90度相位，导致Id、Iq的符号和大小不一致，这是跨教材复现代码时最容易出问题的地方。解决办法以代码里采用的派克变换矩阵为准，把稳态公式重新推导一遍，不要直接照抄其他文献里的Id0公式。初始转角错了，后续功角曲线的绝对值会整体偏移，但相对变化趋势可能看起来正常，这类错误很隐蔽。

4.2 四阶龙格-库塔法的MATLAB实现与步长选择

SMIG_2.m的时间推进部分，常见选择是定步长四阶龙格-库塔法（RK4）：在步长1到10毫秒范围内精度足够，实现简单、易调试。单步函数如下：

function x_next = rk4_step(fun, t, x, h, u, p) % 单步四阶龙格-库塔积分 k1 = fun(t, x, u, p); k2 = fun(t + h/2, x + h/2*k1, u, p); k3 = fun(t + h/2, x + h/2*k2, u, p); k4 = fun(t + h, x + h*k3, u, p); x_next = x + h * (k1 + 2*k2 + 2*k3 + k4) / 6; end

u是这一步起点处的控制量（如Efd）。RK4的局部截断误差为O(h^5)、总体误差O(h^4)，对机电暂态的0.1到10Hz频段，5毫秒步长可以把数值阻尼控制在很小范围。如果用隐式欧拉或改进欧拉法，同样步长下数值阻尼会明显偏大，功角曲线看起来“更平稳”，但这是一种假象，不要据此得出系统阻尼良好的结论。

4.3 积分方法与步长对比

若改用MATLAB内置ode45，虽然能自适应步长，但必须在t_sw和t_cl两个时间事件上用'Events'选项精确中断，否则积分器会跨过故障切点，结果依赖求解器容差设置。这个问题在手写RK4循环中不存在，因为时间事件显式放在循环层判断。

4.4 标幺值基准统一与参数换算

SMIG_2.m中还有一个容易被忽视的细节：标幺值基准的选取。发电机的容量和电压额定值来自铭牌，但线路、变压器电抗往往来自不同基准下的计算书。所有电抗数值在进入仿真前必须归算到同一个基准容量：

% 线路电抗从 100MVA 基准换算到发电机 250MVA 基准 X_line_100 = 0.12; % 以100MVA为基准的标幺值 S_base_gen = 250; % 发电机基准容量 MVA S_base_line = 100; % 原线路基准容量 MVA X_line_gen = X_line_100 * (S_base_gen / S_base_line); % 结果为 0.12 * 250 / 100 = 0.3 pu

换算关系是标幺值电抗与基准容量成正比。如果发电机容量600MW而线路电抗以100MVA基准给出，直接相加会导致等值电抗偏小、电磁功率偏大、故障时的加速面积被低估，最终CCT偏乐观。最稳妥的做法是在参数结构体p中统一记录基准容量S_base，每个电抗、时间常数都显式注明所属基准。

5. 从波形到决策：极限切除时间计算与稳定裕度验证

5.1 用等面积法则快速判读第一摆稳定性

拿到SMIG_2.m输出的功角曲线后，第一步不是看整段曲线是否收敛，而是看第一摆的走向。故障期间发电机加速、功角上升，切除后电磁功率跃升，如果切除时刻对应的角度小于临界切除角，转子有足够的减速面积吸收加速能量，第一摆稳定，之后功角在阻尼作用下衰减到新的平衡点。

数值仿真中，可以直接看功角峰值是否超过约120度，同时观察转速曲线是否越过同步速之后回落。功角和转速两者同时满足条件，基本可以判定第一摆稳定。但也不要忽视电压曲线：如果故障切除后机端电压长期低于0.75pu，说明励磁系统强励能力不足或外电抗偏大，即使功角稳住了，电压稳定性和恢复时间也值得怀疑。

5.2 极限切除时间CCT的二分搜索实现

故障切除时间t_cl越大，加速时间越长，临界失稳的切除时间就是极限切除时间（CCT）。这是衡量保护方案和运行方式的重要指标。二分法实现如下：

% 二分法搜索极限切除时间 CCT t_low = 0.05; % 下限：肯定稳定 t_high = 0.30; % 上限：肯定失稳 max_iter = 20; tol = 0.001; % 精度 0.001 秒 for k = 1:max_iter t_mid = (t_low + t_high) / 2; stable = run_smib_simulation(t_sw, t_mid); % 返回当前切除时间下是否稳定 if stable t_low = t_mid; else t_high = t_mid; end if (t_high - t_low) < tol break; end end CCT = (t_low + t_high) / 2; fprintf('CCT = %.4f s\n', CCT);

run_smib_simulation里判断“稳定”的标准一般是双判据：10秒后功角仍有界（如未超180度），并且转速偏差包络衰减到某一阈值（如小于0.001 pu）。注意单纯凭功角不超过180度来判定失稳有争议：实际工程中会同时校验转速和电磁功率趋势，两个判据同时满足才认为稳定。

5.3 批量工况扫描：把脚本封装成可复用函数

从单次仿真走向批量校核时，需要把SMIG_2.m的主流程封装成函数，输入t_sw、t_cl、参数结构体p，输出稳定状态、最大功角、最低电压等核心指标。随后可以对惯性常数H、阻尼D、励磁增益Ka做二维扫描，绘制等值线图，直接观察稳定域边界随参数的变化规律。

一个更实用的做法是把CCT二分搜索嵌套进参数扫描，输出“CCT随外电抗Xe变化”或“CCT随惯性常数H变化”的曲线。这类结果直接对应运行方式部门的决策需求：某台机组检修退出后系统等值电抗增大多少，CCT还能留多少余量，现有保护能否在极限时间内切除故障。整个扫描流程在MATLAB中运行时间通常在秒级到分钟级，比Simulink反复起停模型高效得多，这正是SMIG_2.m这类脚本资源在实际工作中的价值所在。

本文还有配套的精品资源，点击获取