定步长积分方法对电磁暂态实时仿真的优势

RLC侠收录于

2026-06-17 2026-07-08 约 7700 字预计阅读 16 分钟 - 次阅读

前面几篇文章一直在讲欧拉、梯形法、theta 积分和 Dommel 形式。读到这里，很多人自然会问一个问题：

既然数值分析里有那么多更高级的积分器，比如 ode4、ode45、ode15s，为什么电磁暂态实时仿真软件还经常围着固定步长梯形法、后向欧拉、theta 积分打转？

高阶积分器为什么诱人

先说更高级的算法为什么好。

以经典四阶 Runge-Kutta 为例，工程软件里常把固定步长 RK4 叫作 ode4。对常微分方程：

\frac{\mathrm{d}x}{\mathrm{d}t}=f(t,x)

RK4 在一个步长内会计算四个斜率：

k_1=f(t_k,x_k)

k_2=f\left(t_k+\frac{T}{2},x_k+\frac{T}{2}k_1\right)

k_3=f\left(t_k+\frac{T}{2},x_k+\frac{T}{2}k_2\right)

k_4=f(t_k+T,x_k+Tk_3)

然后组合：

x_{k+1}=x_k+\frac{T}{6}(k_1+2k_2+2k_3+k_4)

它比简单欧拉法精细得多，因为它不是只看一步开始或结束的斜率，而是在一步内部取了多个采样点。

再进一步，像 MATLAB 里的 ode45 这类变步长算法，还会估计局部误差，然后自动调节步长：

误差小，就放大步长，少算几步；
误差大，就缩小步长，保证精度；
遇到平滑段，效率很高；
遇到剧烈变化，也能自动加密计算。

下面给一个很小的 Python 例子。它不模拟某一个具体设备，只是构造了一个带开关边沿的二阶系统。这个系统在 $t=30\,ms$ 附近受到一个快速阶跃激励，类似一个开关动作把系统从平滑段推入暂态段。脚本用 Dormand-Prince 5(4) 公式做了一个简化版 RK45，也就是 ode45 常用的那类自适应思想：同时算一个五阶解和四阶解，用二者差值估计局部误差，再决定下一步走多大。

网页正文把这个例子的完整 Python 脚本放在下面。它包括教学用二阶系统、自适应 RK45、固定步长欧拉、固定步长梯形预测校正、固定步长 RK4，以及两张图的生成代码：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
import numpy as np
import matplotlib.pyplot as plt


def switched_rlc(t, x):
    """A small second-order system with a switching event.

    It is not meant to model a specific device. The point is to show how an
    adaptive RK45-style solver shrinks its step size around a fast switching
    edge.
    """
    y, dy = x
    t_sw = 0.03
    source = 0.5 * (1.0 + np.tanh((t - t_sw) / 1.2e-4))
    omega = 2.0 * np.pi * 420.0
    zeta = 0.08
    ddy = omega * omega * (source - y) - 2.0 * zeta * omega * dy
    return np.array([dy, ddy])


def rk45_step(fun, t, x, h):
    # Dormand-Prince 5(4) coefficients, the same family behind ode45.
    k1 = fun(t, x)
    k2 = fun(t + h * 1 / 5, x + h * (1 / 5) * k1)
    k3 = fun(t + h * 3 / 10, x + h * (3 / 40 * k1 + 9 / 40 * k2))
    k4 = fun(t + h * 4 / 5, x + h * (44 / 45 * k1 - 56 / 15 * k2 + 32 / 9 * k3))
    k5 = fun(
        t + h * 8 / 9,
        x
        + h
        * (
            19372 / 6561 * k1
            - 25360 / 2187 * k2
            + 64448 / 6561 * k3
            - 212 / 729 * k4
        ),
    )
    k6 = fun(
        t + h,
        x
        + h
        * (
            9017 / 3168 * k1
            - 355 / 33 * k2
            + 46732 / 5247 * k3
            + 49 / 176 * k4
            - 5103 / 18656 * k5
        ),
    )

    x5 = x + h * (
        35 / 384 * k1
        + 500 / 1113 * k3
        + 125 / 192 * k4
        - 2187 / 6784 * k5
        + 11 / 84 * k6
    )

    k7 = fun(t + h, x5)
    x4 = x + h * (
        5179 / 57600 * k1
        + 7571 / 16695 * k3
        + 393 / 640 * k4
        - 92097 / 339200 * k5
        + 187 / 2100 * k6
        + 1 / 40 * k7
    )
    err = np.linalg.norm(x5 - x4, ord=np.inf)
    return x5, err


def solve_adaptive(fun, t0, t1, x0, h0=4e-4, rtol=2e-4, atol=1e-6):
    t = t0
    x = np.array(x0, dtype=float)
    h = h0
    ts = [t]
    xs = [x.copy()]
    hs = []

    while t < t1:
        h = min(h, t1 - t)
        x_new, err = rk45_step(fun, t, x, h)
        scale = atol + rtol * max(np.linalg.norm(x, ord=np.inf), np.linalg.norm(x_new, ord=np.inf))
        ratio = err / scale

        if ratio <= 1.0:
            t += h
            x = x_new
            ts.append(t)
            xs.append(x.copy())
            hs.append(h)

        if ratio == 0.0:
            factor = 2.0
        else:
            factor = 0.9 * ratio ** (-0.2)
        factor = min(3.0, max(0.15, factor))
        h = min(1.5e-3, max(2e-7, h * factor))

    return np.array(ts), np.array(xs), np.array(hs)


def solve_fixed_euler(fun, t0, t1, x0, h):
    ts = np.arange(t0, t1 + 0.5 * h, h)
    xs = np.zeros((len(ts), len(x0)))
    xs[0] = np.array(x0, dtype=float)
    for k in range(len(ts) - 1):
        xs[k + 1] = xs[k] + h * fun(ts[k], xs[k])
    return ts, xs


def solve_fixed_trapezoidal_predictor_corrector(fun, t0, t1, x0, h):
    """Explicit predictor + trapezoidal corrector.

    This keeps the example lightweight. It is not a full implicit Newton
    trapezoidal solve, but it is enough to contrast fixed-step behavior with
    adaptive RK45 for this introductory figure.
    """
    ts = np.arange(t0, t1 + 0.5 * h, h)
    xs = np.zeros((len(ts), len(x0)))
    xs[0] = np.array(x0, dtype=float)
    for k in range(len(ts) - 1):
        f0 = fun(ts[k], xs[k])
        pred = xs[k] + h * f0
        f1 = fun(ts[k + 1], pred)
        xs[k + 1] = xs[k] + 0.5 * h * (f0 + f1)
    return ts, xs


def solve_fixed_rk4(fun, t0, t1, x0, h):
    ts = np.arange(t0, t1 + 0.5 * h, h)
    xs = np.zeros((len(ts), len(x0)))
    xs[0] = np.array(x0, dtype=float)
    for k in range(len(ts) - 1):
        t = ts[k]
        x = xs[k]
        k1 = fun(t, x)
        k2 = fun(t + 0.5 * h, x + 0.5 * h * k1)
        k3 = fun(t + 0.5 * h, x + 0.5 * h * k2)
        k4 = fun(t + h, x + h * k3)
        xs[k + 1] = x + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    return ts, xs


ts, xs, hs = solve_adaptive(switched_rlc, 0.0, 0.06, [0.0, 0.0])
fixed_dt = 50e-6
te, xe = solve_fixed_euler(switched_rlc, 0.0, 0.06, [0.0, 0.0], fixed_dt)
tt, xt = solve_fixed_trapezoidal_predictor_corrector(switched_rlc, 0.0, 0.06, [0.0, 0.0], fixed_dt)
tr, xr = solve_fixed_rk4(switched_rlc, 0.0, 0.06, [0.0, 0.0], fixed_dt)

fig, (ax1, ax2) = plt.subplots(
    2,
    1,
    figsize=(8.8, 5.2),
    sharex=True,
    gridspec_kw={"height_ratios": [2.1, 1.2]},
)

ax1.plot(ts * 1000.0, xs[:, 0], color="#2563eb", linewidth=1.9, label="adaptive RK45")
ax1.plot(tr * 1000.0, xr[:, 0], color="#7c3aed", linewidth=1.35, linestyle="-.", label="fixed-step RK4")
ax1.plot(tt * 1000.0, xt[:, 0], color="#16a34a", linewidth=1.25, linestyle="--", label="fixed-step trapezoidal")
ax1.plot(te * 1000.0, xe[:, 0], color="#f97316", linewidth=1.05, linestyle=":", label="fixed-step Euler")
ax1.axvline(30.0, color="#ef4444", linestyle="--", linewidth=1.2)
ax1.set_ylabel("state y")
ax1.set_title("Switching response: adaptive RK45 vs fixed-step methods")
ax1.grid(True, alpha=0.28)
ax1.legend(loc="upper right", frameon=False)
ax1.annotate(
    "switching event",
    xy=(30.0, 0.05),
    xytext=(47.0, 0.48),
    arrowprops={"arrowstyle": "->", "color": "#ef4444"},
    color="#991b1b",
)

step_t = ts[:-1] * 1000.0
ax2.step(step_t, hs * 1e6, where="post", color="#16a34a", linewidth=1.6)
ax2.axvline(30.0, color="#ef4444", linestyle="--", linewidth=1.2)
ax2.set_xlabel("time / ms")
ax2.set_ylabel("accepted step / us")
ax2.grid(True, alpha=0.28)

fig.tight_layout()
fig.savefig("ode45-adaptive-step-switching.svg", format="svg")


ref_ts, ref_xs, _ = solve_adaptive(
    switched_rlc,
    0.0,
    0.06,
    [0.0, 0.0],
    h0=1e-4,
    rtol=1e-8,
    atol=1e-10,
)
ref_y = np.interp(te, ref_ts, ref_xs[:, 0])
err_euler = np.abs(xe[:, 0] - ref_y)
err_trap = np.abs(xt[:, 0] - ref_y)
err_rk4 = np.abs(xr[:, 0] - ref_y)

fig2, bx = plt.subplots(figsize=(8.8, 3.7))

eps = 1e-12
bx.semilogy(te * 1000.0, err_euler + eps, color="#f97316", linewidth=1.15, linestyle=":", label="Euler error")
bx.semilogy(tt * 1000.0, err_trap + eps, color="#16a34a", linewidth=1.25, linestyle="--", label="trapezoidal PC error")
bx.semilogy(tr * 1000.0, err_rk4 + eps, color="#7c3aed", linewidth=1.3, linestyle="-.", label="RK4 error")
bx.axvline(30.0, color="#ef4444", linestyle="--", linewidth=1.1)
bx.annotate(
    "switching event",
    xy=(30.0, 1e-5),
    xytext=(39.0, 2e-3),
    arrowprops={"arrowstyle": "->", "color": "#ef4444"},
    color="#991b1b",
)
bx.set_xlabel("time / ms")
bx.set_ylabel("absolute error")
bx.set_title("Fixed-step method error at the same dt")
bx.grid(True, which="both", alpha=0.28)
bx.legend(loc="upper right", frameon=False, fontsize=8.8)

fig2.tight_layout()
fig2.savefig("fixed-step-accuracy-comparison.svg", format="svg")

自适应 RK45 在开关事件附近自动缩小步长 — 自适应 RK45 在开关边沿附近自动缩小步长

上图里，上半部分给出了四条波形：自适应 RK45、固定步长 RK4、固定步长梯形预测校正和固定步长欧拉法。下半部分给出了 RK45 实际接受的步长。可以看到，平滑段里 RK45 的步长可以变大；到了开关边沿和振荡较明显的区域，步长会自动缩小。这正是变步长算法在离线仿真里诱人的地方：它把计算资源集中到“更难算”的时间段。至于程序具体怎样根据误差估计放大或缩小步长，后面会单独用一篇文章拆开讲；这里先抓住现象和工程含义。

这里还有一个容易误会的点：图中的振荡不是 RK45 算法凭空造出来的数值振荡，而是这个二阶系统在开关激励后的欠阻尼响应。换句话说，物理模型本身就会振。RK45 做的是在振荡更剧烈的地方自动缩小步长，让这段响应被更细地解析出来；它并不会把系统本来存在的物理振荡消掉。

这和前面讨论过的梯形积分数值振荡不是一回事。数值振荡通常是离散算法、开关突变和历史项不协调共同引出的寄生模态；而这里的振荡来自方程本身。判断二者的一个简单办法是：如果换成足够小步长的不同算法，振荡仍然存在，而且频率和衰减趋势一致，那它更可能是模型响应；如果振荡形态强烈依赖积分器和步长，甚至表现为一正一负交替的离散毛刺，那就要警惕数值振荡。

图里固定步长欧拉法的偏差很大，也不是偶然。显式欧拉只用步长起点的斜率外推下一拍：

x_{k+1}=x_k+T f(t_k,x_k)

对于开关后的欠阻尼振荡，状态在一个步长内的斜率变化很快。欧拉法仍然拿“起点斜率”代表整个步长的平均变化，于是每一步都会积累相位误差和幅值误差。对振荡系统来说，这种误差还可能表现得像给系统额外注入或抽走能量，所以波形会明显偏离参考解。RK4 在同一个步长内取了多个斜率，等于更好地估计了这一步的平均变化，因此误差通常小很多。

为了单独看数值精度，再把同一个系统用固定步长欧拉、固定步长梯形预测校正、固定步长 RK4 各算一遍。这里统一取 $\Delta t=50\,\mu s$ ，再用一个更严格容差的自适应 RK45 结果作为参考解。这样比较并不是说参考解就等于解析解，而是为了让读者看到同一条件下不同积分公式的误差量级差异。

固定步长欧拉、梯形预测校正和 RK4 的数值精度对比 — 同一步长下，RK4 的数值误差通常明显小于欧拉法

从误差图可以看出，高阶方法确实有它的价值。RK4 在同样步长下通常比欧拉法精确得多，也常常比简单的梯形预测校正更贴近参考解。这也是高阶积分器在离线仿真里非常受欢迎的原因。

但实时仿真关心的不只是“这一步算得多准”。RK4 每一步需要多次评估 $f(t,x)$ ；如果 $f(t,x)$ 背后是一整套电力网络方程、控制器和开关模型，多评估几次就可能意味着更多网络求解、更多状态更新和更重的实时调度压力。也就是说，高阶方法在精度上有优势，但它并不自动满足实时仿真的固定计算节拍。

如果只是做离线仿真，这些能力非常诱人。你关心的是最后波形准不准、总计算时间能不能接受。程序中间某一步多算一点、少算一点，通常不是大问题。

实时仿真的目标不一样

实时仿真的约束更硬。

假设实时仿真步长为 $T=50\,\mu s$ 。这意味着软件或硬件每隔 $50\,\mu s$ 就必须交出一次结果。如果某一步用了 $70\,\mu s$ 才算完，即使这一小步的数学误差很低，实时系统也已经失败了。

所以实时仿真首先追求：

固定周期；
固定计算负载；
固定输入输出延迟；
最坏情况也能在一个步长内算完；
适合并行化、流水线和硬件映射。

高阶和变步长算法的问题，就出现在这里。

ode4/RK4 每一步要多次调用 $f(t,x)$ 。如果 $f$ 背后是一个大规模电力网络、控制器和电力电子开关模型，多算几次斜率就意味着多次网络求解或多次状态评估。

变步长算法更麻烦。它每一步到底走多大，取决于误差估计；如果误差不满足要求，还可能拒绝当前步，缩小步长重算。这在离线仿真里很合理，但在实时仿真里很危险：真实时间不会等你重算。

固定步长 Dommel 的隐藏优势

对一个线性电感，梯形法给出的等效导纳为：

G_L=\frac{T}{2L}

对一个线性电容，梯形法给出的等效导纳为：

G_C=\frac{2C}{T}

如果步长 $T$ 固定，元件参数 $L$ 、 $C$ 固定，网络拓扑也不变，那么这些等效导纳就是常数。抽象地说，这会让节点方程里的左端矩阵尽量稳定：

Gv_k=b_k

其中 $G$ 是节点导纳矩阵， $v_k$ 是当前步节点电压未知量， $b_k$ 是由电流源、历史源、控制注入等组成的右端项。

只说这句话还是有点空，下面给一个具体数值例子。

假设某个局部网络经过节点消元、等值支路和控制接口装配以后，保留下来 3 个节点。静态电阻、线路等值和已经消元的内部节点，会给出一个稠密的静态导纳矩阵：

G_{\mathrm{static}}= \begin{bmatrix} 0.42 & -0.18 & -0.07\\ -0.18 & 0.36 & -0.11\\ -0.07 & -0.11 & 0.31 \end{bmatrix}

单位可以理解为西门子。这个矩阵之所以是稠密的，不是因为物理上每两个节点都一定有一根线直连，而是因为消元、等效和耦合支路会把原来局部的连接“折叠”到保留下来的节点之间。

再加上几个动态元件：

节点 1 对地有电感 $L_1=5\,\mathrm{mH}$ ；
节点 2 对地有电容 $C_2=20\,\mu\mathrm{F}$ ；
节点 2 和节点 3 之间有电感 $L_{23}=2\,\mathrm{mH}$ ；
节点 3 对地有电容 $C_3=5\,\mu\mathrm{F}$ 。

若固定步长取：

T=50\,\mu s

那么：

G_{L1}=\frac{T}{2L_1}=0.005

G_{C2}=\frac{2C_2}{T}=0.8

G_{L23}=\frac{T}{2L_{23}}=0.0125

G_{C3}=\frac{2C_3}{T}=0.2

装配以后，当前步左端矩阵就是：

G_{50\mu s}= \begin{bmatrix} 0.425 & -0.180 & -0.070\\ -0.180 & 1.1725 & -0.1225\\ -0.070 & -0.1225 & 0.5225 \end{bmatrix}

这里每一个数都可以追踪。例如节点 2 的对角项从 $0.36$ 变成 $1.1725$ ，是因为加上了对地电容 $0.8$ 和节点 2-3 电感支路贡献 $0.0125$ ；节点 2 和节点 3 之间的互导纳从 $-0.11$ 变成 $-0.1225$ ，是因为支路 $L_{23}$ 又贡献了 $-0.0125$ 。

固定步长的好处在这里就非常具体了。只要拓扑和参数不变， $G_{50\mu s}$ 每一拍都一样。程序可以在实时循环开始前做一次：

确定矩阵非零结构；
做节点重排序，减少分解填充；
做 LU 分解，得到 $G=LU$ ；
在每个实时步里只更新 $b_k$ ；
然后解两次三角方程 $Ly=b_k$ 、 $Uv_k=y$ 。

也就是说，实时循环里不是每一步真的去算 $G^{-1}$ ，也不是每一步重新从头分解矩阵，而是尽量把昂贵、不可预测的部分放到循环外面。

如果改成变步长，情况就变了。设下一步误差估计要求把步长缩小到：

T_k=25\,\mu s

那么电感等效导纳会减半，电容等效导纳会加倍：

G_{L1}=0.0025,\qquad G_{C2}=1.6,\qquad G_{L23}=0.00625,\qquad G_{C3}=0.4

此时左端矩阵变成：

G_{25\mu s}= \begin{bmatrix} 0.4225 & -0.180 & -0.070\\ -0.180 & 1.96625 & -0.11625\\ -0.070 & -0.11625 & 0.71625 \end{bmatrix}

如果某一步又把步长放大到：

T_k=100\,\mu s

则：

G_{100\mu s}= \begin{bmatrix} 0.430 & -0.180 & -0.070\\ -0.180 & 0.785 & -0.135\\ -0.070 & -0.135 & 0.435 \end{bmatrix}

注意，这三个矩阵的非零位置看起来没有变，但数值已经变了。对于这个 3 阶小矩阵，重新分解当然不算什么；但在实时 EMT 里，矩阵可能是几百阶、几千阶，甚至更多。稠密矩阵做一次 LU 分解的计算量大约是 $O(n^3)$ ，而分解以后每步做前代、回代大约是 $O(n^2)$ 。如果 $n=1000$ ，一次稠密 LU 分解大约是 $\frac{2}{3}n^3\approx 6.7\times 10^8$ 量级的浮点操作，而两次三角回代大约是 $2n^2\approx 2\times 10^6$ 量级。两者不是一个档位。

真实电力网络常常是稀疏矩阵，不能直接拿稠密矩阵的 $O(n^3)$ 当精确耗时，但结论仍然成立：能复用分解，就尽量复用分解；只更新右端项并回代，通常比每步重建、重排、重分解可控得多。

这就是固定步长 Dommel 方法在实时仿真里特别吃香的原因之一。它不是单纯在“积分公式”层面好用，而是在“计算结构”层面也很规整。

固定步长更容易对齐实时硬件时序

矩阵复用只是固定步长适合实时仿真的一部分原因。另一个很实际的原因是：实时仿真不是单纯在电脑里算完一条曲线，它还要和外部硬件按固定节拍交换数据。

以硬件在环 HIL 或实时 EMT 平台为例，一个仿真步通常不是“求解器自己想算多久就算多久”，而是被一个严格的时钟框住。假设步长仍然是：

T=50\,\mu s

那么每一拍大致都要完成类似这样的流程：

采样外部模拟量、数字量或通信量；
锁存上一拍的输入，避免本拍计算过程中输入乱跳；
根据当前网络状态和控制状态装配右端项 $b_k$ ；
回代求出节点电压、电流和内部状态；
更新开关逻辑、保护逻辑和控制器；
把输出写到 D/A、PWM、光纤接口、I/O 板卡或通信缓存；
等待下一个同步时钟边沿。

如果固定步长是 $50\,\mu s$ ，平台就可以把这段时间切成相对稳定的时序窗口。例如前 $5\,\mu s$ 采样和锁存输入，中间 $35\,\mu s$ 做网络求解和控制更新，最后 $10\,\mu s$ 写输出、留保护裕量和处理同步。真实平台的分配当然会更复杂，但核心思想就是：每一拍的输入、计算、输出都能围绕同一个固定节拍排队。

这对 FPGA 尤其重要。FPGA 擅长做固定流水线：第一级读数据，第二级做乘加，第三级做局部求和，后面再做回代、限幅、接口输出。只要每一拍的计算路径稳定，流水线就可以提前安排好，延迟也可以提前数清楚。比如某个求解流水线固定需要 300 个 FPGA 时钟周期，I/O 同步需要 80 个周期，控制计算需要 120 个周期，只要总周期数小于一个仿真步允许的时钟周期数，就可以稳定运行。

变步长会让这个问题麻烦很多。它不只是让数值积分的 $T_k$ 改变，还会让实时调度变得尴尬：

这一拍如果误差太大，需要缩小步长重算，但外部硬件的真实时间不会倒回去；
下一拍如果允许放大步长，输出采样点又不再和固定 I/O 节拍自然对齐；
FPGA 流水线希望每拍走同样的路径，但变步长可能改变矩阵数值、迭代次数或接受/拒绝逻辑；
HIL 接口、保护装置和控制器通常期待固定采样周期，频繁变化的仿真步长会让接口协议和延迟补偿更难处理。

所以实时仿真里的“固定步长”不只是数值算法偏好，也是系统工程约束。它让求解器、I/O、通信、控制器和 FPGA/多核 CPU 的调度都围绕同一个节拍工作。这样做未必在每一个局部时刻都是数学上最省步数的，但它更容易满足实时系统最看重的要求：固定延迟、确定时序、最坏情况可预估。

小结

这一篇我们没有深入推导某一种高阶算法，而是讨论了算法选择背后的工程原因：

ode4/RK4、ode45 等方法精度高、适应性强，适合很多离线仿真场景。
实时仿真要求固定周期、确定延迟和最坏情况可控。
固定步长 Dommel 方法让 $G_L=T/(2L)$ 、 $G_C=2C/T$ 这类等效导纳保持常数。
当拓扑和参数不变时，节点导纳矩阵 $G$ 可以预处理、预分解和复用。
变步长会让等效导纳和矩阵结构更难固定，从而削弱实时调度和加速能力。
梯形法和 theta 积分之所以常用，不是因为它们最高级，而是因为它们在精度、稳定性和实时可实现性之间取得了很好的平衡。

以后再看实时仿真软件的积分方法，就不要只问“它为什么不用更高阶算法”。更准确的问题是：这个算法能不能在每一个固定节拍里，把网络、控制和接口全部准时算完。

定步长积分方法对电磁暂态实时仿真的优势

高阶积分器为什么诱人

实时仿真的目标不一样

固定步长 Dommel 的隐藏优势

固定步长更容易对齐实时硬件时序

小结

相关内容

目录