简介计算机控制仿真选择计算.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/s1~10 s0.1~1 sT 过大导致稳态值偏离电机速度环10~100 rad/s1~10 ms0.2~1 ms保持器等效延时吃相位裕度电流环、开关电源1k~10k rad/s10~50 µs1~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, n400): 以 T 为采样周期用零阶保持器离散化后递推单位阶跃响应 Ad, Bd, Cd, Dd, _ signal.cont2discrete( (plant.A, plant.B, plant.C, plant.D), T, methodzoh) 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(fT{T:5} ISE{ise:.5f} 末值{y[-1]:.5f})代码逻辑是标准的离散状态递推x[k1] Ad·x[k] Bd·u[k]y[k] Cd·x[k] Dd·u[k]。cont2discrete的methodzoh保证离散模型的阶跃响应在采样点上与连续系统完全一致这正是零阶保持器的物理含义。参数上步数 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 保持器的真实回路输出恒滞后一拍高频段相位不准双线性Tustins (2/T)·(z-1)/(z1)左半平面映射到单位圆内无混叠需要保留频域形状的滤波器、控制器频率轴畸变需按 ω (2/T)·tan(ωT/2) 预畸变零极点匹配z_i exp(s_i·T)零极点位置与增益低阶、零极点明确的被控对象阶数高时要补零点否则高频增益不对零阶保持器法可以手算一阶惯性环节 G(s) b/(sa) 的离散结果是 G(z) b(1-e^{-aT}) / [a(z-e^{-aT})]。代入 a 2、b 1、T 0.05 se^{-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, methodzoh) 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(fh{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 调用次数全局误差阶线性问题发散边界典型用途前向欧拉1O(h)|hλ| 2快速原型、h 极小改进欧拉2O(h²)|hλ| 2教学、中等精度经典 RK44O(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 仿真发散了怎么查五个高频原因现象常见原因排查动作几步内数值冲到 1e30h 超出显式方法的收敛边界步长减半重跑若消失即为步长问题控制量抖振、幅度与步长相关保持器延时叠加微分项放大噪声减小 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。本文还有配套的精品资源点击获取