2. 电力系统动态仿真数学基础

电力系统动态仿真的本质，是把发电机转子运动、励磁调节、原动机调速、负荷动态等连续过程，与网络节点功率平衡这一瞬时约束联立求解。前者是微分方程，后者是代数方程，二者耦合形成微分代数方程（DAE）。理解这一数学结构，是正确选择MATLAB求解器、判断仿真结果可信度的前提。

一、从物理对象到微分代数方程

以单机无穷大系统为例，设发电机采用二阶经典模型，状态量为功角 δ 与转速 ω，惯性时间常数 M、阻尼系数 D、机械功率 P_m、电磁功率 P_e，则转子运动方程为：

dδ/dt = ω - ω_s
M·dω/dt = P_m - P_e - D·(ω - ω_s)

其中电磁功率 P_e 由网络潮流决定，而网络节点功率平衡方程不含时间导数，属于代数约束：

P_i = Σ U_i·U_j·(G_ij·cosθ_ij + B_ij·sinθ_ij)
Q_i = Σ U_i·U_j·(G_ij·sinθ_ij - B_ij·cosθ_ij)

把状态变量记为 x（δ、ω 及励磁、调速等内部变量），代数变量记为 y（节点电压幅值与相角），DAE 的标准形式为：

x' = f(x, y, t)
0  = g(x, y, t)

这里 f 描述动态元件的演化，g 描述网络的瞬时平衡。需要特别注意的是：y 不是由微分方程积分得到的，而是在每个时刻由 g=0 隐式求解出来的。因此仿真推进的每一步，都要在给定 x 的条件下解一次代数方程组，这就是电力系统仿真比一般ODE问题更复杂的原因。

二、状态空间模型与DAE的相互转化

当网络规模较小、且所有动态元件都能显式写出导数时，可把代数约束消去，得到纯状态空间形式 x' = Ax + Bu。线性化后，A 矩阵的特征值直接决定系统的小干扰稳定性：实部为负表示该模式衰减，实部为正表示失稳，虚部对应振荡频率。这是特征值分析的基础。

但实际电网中，负荷的恒功率特性、ZIP负荷、FACTS装置、以及大量节点电压约束，使代数方程难以解析消去。此时保留 DAE 形式更自然，也更接近工程实际。DAE 的另一个重要概念是指数（index）：若对代数约束求一次时间导数即可解出 y'，则称为 index-1 系统；电力系统机电暂态模型通常属于 index-1，但电磁暂态模型可能达到更高指数，求解难度显著上升。

三、刚性：电力系统仿真的核心难点

电力系统动态过程的时间尺度跨度极大：发电机转子运动的时间常数约为秒级，励磁绕组为0.1～1 s，而网络电磁暂态、线路分布参数过程可达微秒级。当这些快慢过程同时出现在一个模型中，系统就表现出刚性（stiffness）。

刚性的直接后果是：显式方法（如 ode45、RK4）为保持数值稳定，步长必须小于最快时间常数的量级，即使我们只关心秒级的功角摇摆，也不得不以微秒级步长推进，计算量无法接受。刚性比可用最大与最小时间常数之比衡量，电力系统中该比值常达 10³～10⁶。

四、数值解法与MATLAB求解器选择

MATLAB 的 ode 系列求解器按适用场景可分为三类：

对于 DAE 形式，MATLAB 提供两种处理方式。第一种是质量矩阵形式：把方程写成 M(t,x)·x' = f(t,x)，其中代数方程对应的行质量矩阵为奇异（全零），调用 ode15s 时通过 odeset('Mass', M) 传入。第二种是半显式形式：直接使用 ode15s 配合 'Mass' 或 ode15i 求解完全隐式 F(t,x,x')=0。

% 质量矩阵形式的DAE示例（示意）
M = [1 0 0; 0 1 0; 0 0 0];   % 第三行为代数约束
opts = odeset('Mass', M, 'RelTol', 1e-6, 'AbsTol', 1e-8);
[t, x] = ode15s(@dae_fun, [0 10], x0, opts);

function dx = dae_fun(t, x)
    delta = x(1); omega = x(2); U = x(3);
    Pe = U * sin(delta) / 0.5;        % 简化网络关系
    dx = zeros(3,1);
    dx(1) = omega - 1;
    dx(2) = (0.8 - Pe - 0.1*(omega-1)) / 6;
    dx(3) = 1.0 - U - 0.5*cos(delta); % 代数约束 g=0
end

选择求解器时，建议按以下顺序判断：先估计系统刚性比，若大于 10³ 优先 ode15s；若模型含开关、限幅等不连续环节，可考虑 ode23tb 或 ode23s；若为纯非刚性小系统，ode45 更高效。容差设置上，机电暂态通常取 RelTol=1e-4～1e-6、AbsTol=1e-6～1e-8，过松会掩盖振荡衰减，过紧则显著增加计算量。

五、结果判读与常见问题

仿真结束后，应重点观察三类信息：一是功角曲线是否单调或呈衰减振荡，若持续等幅或发散，说明阻尼不足或失稳；二是求解器返回的步长与事件信息，若步长长期处于下限，提示刚性过强或容差过严；三是代数残差 g(x,y) 是否始终接近零，若残差增大，说明代数方程未被正确满足，结果不可信。

常见问题包括：把刚性系统交给 ode45 导致计算停滞；质量矩阵维度或奇异行设置错误，使代数约束被当作微分方程处理；以及忽略初值一致性，导致仿真起始段出现非物理尖峰。排查时，可先用固定步长小算例验证方程本身，再切换到变步长刚性求解器对比。

综上，电力系统动态仿真的数学基础可概括为：用 DAE 描述“动态元件 + 网络约束”的耦合系统，用刚性理论解释为什么必须采用隐式方法，用 ode15s 等求解器完成数值推进，并通过功角曲线、步长行为和代数残差三方面判断结果是否可信。掌握这条主线，后续建模与仿真调试才有稳固的落脚点。