简介：计算机控制仿真选择+计算.pdf 是一份面向自动化、控制工程及相关专业学生与备考人员的仿真知识点梳理资料，围绕控制系统数字仿真的题库式选择与计算内容展开，适合课堂复习、期末突击与考研复试前查漏补缺。压缩包内共1个PDF文件，约981KB，体量轻便，便于手机或电脑随时查阅。内容以填空题形式串联系统定义、实体属性活动、连续与离散系统分类、物理模型与数学模型、线性与非线性模型、仿真模型校核与验证等基础概念，并延伸至零阶/一阶保持器、龙格-库塔法局部截断误差、条件稳定与绝对稳定、单步法与多步法、离散相似法与香农定理、增广矩阵法等快速数字仿真算法，以及MATLAB中c2d函数、双线性替换法、采样控制系统数字仿真模型和差分方程递推求解。已有43人学习，可帮助读者快速定位概念盲区，形成从建模、离散化到数值积分与稳定性分析的完整复习脉络。
1. 一份「计算机控制仿真选择+计算」的 PDF，真正要解决的是什么
很多人硬盘里都存着几份名为「计算机控制仿真选择+计算」的 PDF 或课件，公式看着眼熟，抄下来一跑却对不上：连续域里调得很听话的 PID，按采样周期一离散就开始抖；仿真步长从 1 ms 改成 10 ms，超调凭空多出十几个百分点；同一套参数在 Python 里跑得好好的，烧进单片机就变成了小幅自持振荡。
问题通常不在公式抄错了，而在标题里那两个动词。「选择」指的是选仿真域（连续近似、采样混合还是全离散）、选离散化变换（零阶保持器、双线性还是零极点匹配）、选积分方法（欧拉、RK4 还是隐式变步长）；「计算」指的是把 G(s) 换成能一行行递推的差分方程，把每个采样时刻的状态更新、控制量和限幅都算准，并且知道浮点舍入、步长和量化各自贡献了多少误差。这份内容写给已经会写控制器、但仿真曲线和上机曲线对不上的工程师，也写给需要从课件文档里把公式还原成代码的学生。
2. 计算机控制仿真里的「选择」：连续域、离散域与采样周期
真实现场里，被控对象是连续的，控制器是离散的，两者之间夹着采样器和零阶保持器。仿真时把这三段怎么摆，决定了你能看到什么、看不到什么。选错域，后面所有参数整定都是白费力气。
2.1 连续域仿真与离散域仿真的差别到底在哪
最常见的三条路线，各自的取舍差别很大。
第一条是全连续近似：把数字控制器用一个等效连续传递函数代替，直接用 ode45 之类的连续求解器跑。工具现成、上手快，代价是把采样器和保持器藏起来了，采样引入的等效延时完全暴露不出来。它只适合做早期控制律验证。
第二条是采样混合仿真：被控对象保留连续状态方程，用固定步长 h 积分，每到 T 时刻采一次样、算一次控制量，再用零阶保持器把控制量保持到下一个采样点。它最接近真实系统，也最容易暴露相位滞后和采样引起的振荡，是我在整定阶段默认选的一条路。
第三条是全离散仿真：把对象也用零阶保持器离散掉，整套系统退化成纯差分递推，没有积分、没有连续求解器。计算量最小，也正好是最终要写进 MCU 的形式。开发后期我会切到这条路，用同一组参数复核一遍。
2.2 采样周期 T 与仿真步长 h：两个量，别混着用
T 是物理意义上的控制周期，由传感器刷新率、执行器响应和 CPU 预算决定；h 是数值积分步长，只影响解算精度。两者混用是仿真翻车的头号原因：把 h 直接取成 T，等于用一个采样周期去做一次积分，精度和收敛性都没有余量。
经验取值是 h ≤ T/10，用欧拉法时收紧到 T/20，用 RK4 可以放宽到 T/5。对连续部分还有一条独立约束：h·ωmax ≤ 0.1~0.2，其中 ωmax 取闭环带宽的 5~10 倍，用来覆盖高次模态。
| 被控对象类型 | 典型闭环带宽 | 建议 T | 建议 h | 主要风险 |
|---|---|---|---|---|
| 温度、液位等慢过程 | 0.1~1 rad/s | 1~10 s | 0.1~1 s | T 过大导致稳态值偏离 |
| 电机速度环 | 10~100 rad/s | 1~10 ms | 0.2~1 ms | 保持器等效延时吃相位裕度 |
| 电流环、开关电源 | 1k~10k rad/s | 10~50 µs | 1~5 µs | 测量噪声被微分项放大 |
还有一条容易忽略的账：零阶保持器带来的等效延时约为 T/2，对应相位损失约 ωc·T/2。带宽 100 rad/s、T 取 5 ms 时相位损失 0.25 rad，也就是 14 度左右，这个数字必须从相位裕度里预先扣掉。
2.3 用一段代码看采样周期对同一对象的影响
下面这段代码用零阶保持器把同一个连续对象按不同 T 离散化，递推单位阶跃响应并算误差平方积分，观察 T 的影响。
import numpy as np from scipy import signal # 被控对象 G(s) = 1 / (0.5s + 1)，写成状态空间便于离散化 plant = signal.StateSpace([[-2.0]], [[1.0]], [[1.0]], [[0.0]]) def step_response(T, n=400): """以 T 为采样周期，用零阶保持器离散化后递推单位阶跃响应""" Ad, Bd, Cd, Dd, _ = signal.cont2discrete( (plant.A, plant.B, plant.C, plant.D), T, method='zoh') x = np.zeros(1) y = np.zeros(n) for k in range(n): y[k] = (Cd @ x + Dd * 1.0)[0] # 当前输出，输入恒为 1 x = Ad @ x + Bd * 1.0 # 状态推进一个采样周期 return y for T in (0.5, 0.2, 0.05, 0.01): y = step_response(T) e = 1.0 - y ise = float(np.sum(e ** 2) * T) # 手写矩形累加，避开版本差异 print(f"T={T:<5} ISE={ise:.5f} 末值={y[-1]:.5f}")代码逻辑是标准的离散状态递推：x[k+1] = Ad·x[k] + Bd·u[k]，y[k] = Cd·x[k] + Dd·u[k]。cont2discrete的method='zoh'保证离散模型的阶跃响应在采样点上与连续系统完全一致，这正是零阶保持器的物理含义。参数上，步数 n 要满足 n·T 远大于 5 倍时间常数（本例 τ = 0.5 s，所以 n·T 至少 2.5 s），否则记录不到稳态；ISE 用矩形累加近似 ∫e²dt，和梯形积分的差别在小 T 下可以忽略。
结果会显示 T = 0.5 s 时 ISE 明显偏大、末值也偏离 1，因为此时 T 和 τ 同量级，保持器延时已经主导了动态；T 降到 0.05 s 以下基本收敛。想快速判断自己手上的 T 是否合适，把 T 减半再跑一次，曲线肉眼可分就说明还没收敛。
3. 从 G(s) 到差分方程：计算机控制仿真的「计算」主线
选好了域和步长，剩下的全是计算。这一步最容易出错的不是数学，而是「哪一拍用哪个量」——零阶保持器天然带一拍延时，写错下标，仿真出来的相位裕度会比实际乐观一大截。
3.1 三种离散化变换怎么选：零阶保持器、双线性、零极点匹配
| 变换方法 | 变换关系 | 保住的特性 | 适用场景 | 注意点 |
|---|---|---|---|---|
| 零阶保持器 | G(z) = (1-z⁻¹)·Z[G(s)/s] | 采样点阶跃响应一致 | 有 DAC + 保持器的真实回路 | 输出恒滞后一拍，高频段相位不准 |
| 双线性（Tustin） | s = (2/T)·(z-1)/(z+1) | 左半平面映射到单位圆内，无混叠 | 需要保留频域形状的滤波器、控制器 | 频率轴畸变，需按 ω = (2/T)·tan(ωT/2) 预畸变 |
| 零极点匹配 | z_i = exp(s_i·T) | 零极点位置与增益 | 低阶、零极点明确的被控对象 | 阶数高时要补零点，否则高频增益不对 |
零阶保持器法可以手算，一阶惯性环节 G(s) = b/(s+a) 的离散结果是 G(z) = b(1-e^{-aT}) / [a(z-e^{-aT})]。代入 a = 2、b = 1、T = 0.05 s：e^{-0.1} = 0.9048，增益为 1×(1-0.9048)/2 = 0.0476，于是 G(z) = 0.0476/(z - 0.9048)，对应的差分方程是 y[k] = 0.9048·y[k-1] + 0.0476·u[k-1]。
注意：这里的输入是 u[k-1] 而不是 u[k]，这一个下标就是零阶保持器的物理延时。很多资料为了书写方便把它略掉，抄进代码就会凭空多出或少掉一拍延时。
3.2 PID 离散化的两种写法与积分饱和处理
位置式和增量式在数学上等价，在工程上差别很大。
位置式写作 u(k) = Kp·e(k) + Ki·T·Σe(j) + (Kd/T)·[e(k) - e(k-1)] + u₀，需要一个持续累加的积分器。它的好处是抗饱和好做，限幅时把超出部分回算进积分器即可；坏处是手动/自动切换要重置累加器，浮点累加误差也会随时间漂移。
增量式写作 Δu(k) = Kp·[e(k)-e(k-1)] + Ki·T·e(k) + (Kd/T)·[e(k)-2e(k-1)+e(k-2)]，u(k) = u(k-1) + Δu(k)。它只需要 u(k-1) 和两拍误差历史，天然支持无扰切换，适合定点运算和资源受限的 MCU；代价是限幅时积分仍在继续累积，需要额外做限幅回算。
另一个细节是微分先行：把微分项作用在测量值 y 上而不是误差 e 上，写成 -(Kd/T)·[y(k)-y(k-1)]，这样设定值阶跃时不会产生微分冲击。
3.3 一个可复现的增量式 PID 闭环仿真
import numpy as np from scipy import signal # 被控对象 G(s) = 1 / (s^2 + 1.2s + 1)，欠阻尼二阶 plant = signal.StateSpace([[0.0, 1.0], [-1.0, -1.2]], [[0.0], [1.0]], [[1.0, 0.0]], [[0.0]]) T = 0.02 # 控制周期 20 ms Ad, Bd, Cd, Dd, _ = signal.cont2discrete( (plant.A, plant.B, plant.C, plant.D), T, method='zoh') Kp, Ki, Kd = 1.2, 8.0, 0.05 # 增量式参数，Ki 项里再乘 T u_max, u_min = 5.0, -5.0 r = 1.0 x = np.zeros(2) u, e1, e2 = 0.0, 0.0, 0.0 ys, us = [], [] for k in range(1500): y = float((Cd @ x + Dd * u)[0]) # 先算输出 e = r - y # 再算误差 du = Kp * (e - e1) + Ki * T * e + (Kd / T) * (e - 2 * e1 + e2) u = float(np.clip(u + du, u_min, u_max)) # 位置限幅，最简抗饱和 ys.append(y); us.append(u) x = Ad @ x + Bd * u # 最后推进状态 e2, e1 = e1, e # 误差历史整体后移 ys = np.array(ys) print(f"超调={ys.max() - 1.0:.4f} 末值误差={1.0 - ys[-1]:.2e} 峰值控制量={max(us):.3f}")递推顺序必须严格是「算 y → 算 e → 算 u → 推进状态 → 移误差历史」。顺序写反会引入额外一拍延时，测出来的振荡频率和裕度都会偏离真实系统，这是纯数字仿真最常见的一类隐性错误。
参数方面，Ki 项已经显式乘了 T，所以 Ki 本身可以按连续域的量级给（这里 8.0），不需要再折算；Kd/T 在大 T 下会被放大，T = 0.02 s 时 Kd/T = 2.5，当 T 进一步缩小这个增益会继续变大，测量噪声会被指数级放大。常见做法是给微分项加一阶低通，形式为 D(z) = Kd·(1-z⁻¹) / [T·(1-αz⁻¹)]，α 取 0.8~0.95。np.clip只做了位置限幅，如果发现限幅期间超调仍然变大，说明积分还在累积，需要加上限幅回算：把被削掉的 du 从积分累加量里扣除。
4. 数值积分法的取舍：欧拉、RK4 与刚性问题
被控对象里只要出现多个时间尺度相差很大的模态，积分方法的选择就会直接决定仿真能不能跑完。这一节把三种常用方法的实现、误差阶和发散边界放在一起对比。
4.1 三种积分法的最小实现与误差量级
import numpy as np lam = 10.0 # dy/dt = -lam*y，解析解 y(t) = exp(-lam*t) f = lambda t, y: -lam * y exact = lambda t: np.exp(-lam * t) def euler(f, y0, t0, tf, n): h = (tf - t0) / n y, t = y0, t0 for _ in range(n): y = y + h * f(t, y); t += h return y def rk4(f, y0, t0, tf, n): h = (tf - t0) / n y, t = y0, t0 for _ in range(n): k1 = f(t, y) k2 = f(t + h / 2, y + h * k1 / 2) k3 = f(t + h / 2, y + h * k2 / 2) k4 = f(t + h, y + h * k3) y = y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6 t += h return y for n in (10, 50, 200, 1000): ye = euler(f, 1.0, 0.0, 1.0, n) yr = rk4(f, 1.0, 0.0, 1.0, n) print(f"h={1.0/n:<8.4f} 欧拉误差={abs(ye-exact(1)):.2e} RK4误差={abs(yr-exact(1)):.2e}")两个函数都是显式单步法，接口一致，唯一区别是每步调用 f 的次数：欧拉一次，RK4 四次。欧拉的局部截断误差是 O(h²)、全局 O(h)，RK4 分别是 O(h⁵) 和 O(h⁴)。打印出来的结果会显示，h 缩小到 1/1000 时欧拉误差才降到 1e-3 量级，而 RK4 在 h = 1/50 时就已经到了 1e-6。
| 方法 | 每步 f 调用次数 | 全局误差阶 | 线性问题发散边界 | 典型用途 |
|---|---|---|---|---|
| 前向欧拉 | 1 | O(h) | |hλ| < 2 | 快速原型、h 极小 |
| 改进欧拉 | 2 | O(h²) | |hλ| < 2 | 教学、中等精度 |
| 经典 RK4 | 4 | O(h⁴) | 约 |hλ| < 2.78 | 通用非刚性系统 |
| 后向欧拉 / 梯形 | 需迭代 | O(h) / O(h²) | 对线性刚性问题无步长上界 | 刚性系统 |
4.2 步长、计算复杂度和实时预算怎么权衡
RK4 每步贵 4 倍，但达到同样精度所需的步数往往少几十倍，总计算量反而更省。真正的约束来自实时仿真：一个控制周期 T 内必须完成 m 个积分步，于是有 h = T/m，并且要满足 n_f × m × cost(f) < 0.5 × CPU 预算，留一半余量给通信、保护和调度开销。
空间复杂度怎么计算也不难估：显式方法只需要状态向量 y ∈ Rⁿ 和几个中间向量 k，占用是 O(n)；隐式方法要组装并分解雅可比矩阵，占用 O(n²)。阶数上去以后，隐式方法往往先撞内存再撞算力，这就是嵌入式里优先用显式方法的原因之一。
4.3 仿真发散了怎么查：五个高频原因
| 现象 | 常见原因 | 排查动作 |
|---|---|---|
| 几步内数值冲到 1e30 | h 超出显式方法的收敛边界 | 步长减半重跑，若消失即为步长问题 |
| 控制量抖振、幅度与步长相关 | 保持器延时叠加微分项放大噪声 | 减小 T，给微分项加一阶低通 |
| 稳态附近小幅自持振荡 | 量化误差或限幅引起的极限环 | 放宽限幅值，检查整型溢出 |
| 慢变量发散、快变量正常 | 刚性系统，显式方法步长被迫极小 | 换隐式或变步长求解器 |
| 结果与解析解差固定比例 | 单位不统一（ms 与 s）或漏乘 T | 逐项核对量纲和增益 |
第一条最值得展开。对象特征值为 λ = -100 时，前向欧拉要求 h < 2/100 = 0.02 s。如果采样周期 T 恰好取 0.02 s，又把仿真步长也取成 0.02 s，就等于正好压在发散边界上：参数稍微一动就跑飞。这类问题在把 T 当步长用的代码里非常普遍。
5. 用高精度计算与交叉验证确认仿真结果可信
5.1 误差是从算法来的还是从浮点来的
单精度浮点只有 24 位有效位，约合 7 位十进制。位置式 PID 的积分累加器在 T = 1 ms、连续跑 1 小时的情况下要累加 3.6×10⁶ 次，低位舍入会被逐步放大，这就是单精度 MCU 上积分项缓慢漂移的来源。常见做法是把累加器提升到 float64 或 32 位定点，并周期性重整。
二进制计算里 0.1 这类十进制小数本身就没有精确表示，离散化系数 e^{-aT} 每算一次都带舍入。想判断误差归属，用 decimal 或 fractions 对关键系数做一次高精度复算，再和浮点结果比对，就能区分是算法引入的还是舍入引入的。
矩阵分析也要看一眼。Ad = exp(A·T) 在 scipy 里默认用 Padé 近似，当 A 的特征值实部很大而 T 又不小时，Ad 的条件数会迅速变差，用np.linalg.cond(Ad)看，超过 1e10 就要警惕，此时的状态递推对舍入非常敏感。
5.2 三条交叉验证路线
import numpy as np from scipy import signal from scipy.linalg import expm A = np.array([[-2.0]]); B = np.array([[1.0]]) T = 0.05 Ad_num, Bd_num, _, _, _ = signal.cont2discrete( (A, B, np.array([[1.0]]), np.array([[0.0]])), T) Ad_ref = expm(A * T) # 矩阵指数参考值 Bd_ref = np.linalg.solve(A, Ad_ref - np.eye(1)) @ B # ZOH 精确积分式 print("Ad 绝对误差:", np.abs(Ad_num - Ad_ref).max()) print("Bd 相对误差:", np.abs(Bd_num - Bd_ref).max() / abs(Bd_ref).max()) print("Ad 条件数:", np.linalg.cond(Ad_ref)) print("离散直流增益:", float(np.linalg.inv(np.eye(1) - Ad_ref) @ Bd_ref)) # 应为 1/2第一路线是解析对照：Ad 用矩阵指数的参考实现复算，Bd 用 A⁻¹(Ad - I)B 这个精确积分式复算，两者与库函数结果的差值就是库内部近似的量级。第二路线是不变量校验：离散系统的直流增益 inv(I - Ad)·Bd 必须等于连续系统的 G(0) = 1/a = 0.5，稳态增益对不上，说明离散化或量纲出了问题。第三路线是步长减半收敛验证：同一算例用 h 和 h/2 各跑一次，把两条曲线的最大差值打印出来，只有小于 1e-6 才能认为数值已经收敛。
步长减半还能白拿一阶精度：误差按 O(h^p) 衰减时，Richardson 外推 y ≈ (2^p·y_{h/2} - y_h)/(2^p - 1) 可以把结果再推高一阶。把 h/2 与 h 的差值直接打印并设成 1e-6 的门限，比对着曲线目测收敛可靠得多——曲线看着平滑的时候，误差可能还停在 1e-2。
本文还有配套的精品资源，点击获取
