简介这份PDF是面向自动化、控制工程等专业学生的计算机控制仿真复习资料侧重选择题与计算题训练可用于课程期末、考研专业课或仿真实验前的查漏补缺。资源包仅含1个PDF文件约981KB内容以题库、概念辨析和计算题型为主便于打印或电子设备随时翻阅。已有43人学习下载。资料覆盖系统定义与分类、模型校核与验证、连续与离散系统仿真、保持器、数值积分及截断误差、龙格-库塔法、离散相似法、双线性变换、采样控制系统、MATLAB中c2d函数等核心考点并延伸至参数优化、目标函数与单纯形寻优等内容。读者可通过题目练习梳理仿真建模流程强化公式推导、稳定性判断和误差分析能力适合作为知识点自测与考前速查的配套材料。1. 控制系统仿真到底在算什么从一次数值发散说起很多人第一次用二阶龙格-库塔法解 $y -y$看到步长取 2.5 时结果直接炸掉才意识到仿真不是把方程丢给计算机这么简单。计算机控制仿真这道门槛本质上是三件事的叠加把物理系统抽象成数学模型一次模型化把数学模型转成计算机能跑的差分递推式二次模型化再在稳定性、精度、速度之间做妥协。这份题库覆盖的正是这三个环节的考点与算例——系统边界怎么划、保持器怎么选、龙格-库塔法步长如何限制、双线性变换为什么绝对稳定、根匹配法在纯延迟环节上怎么处理。它适合两类人一类是正在准备控制系统仿真、计算机仿真课程考试的学生需要把填空和计算题背后的推导链条补齐另一类是已经上手 MATLAB 做数字控制器实现的工程师想回头确认 c2d 的 tustin 和 matched 两种 method 到底差在哪、为什么离散模型会引入附加零点。下面不复述题库而是把每类题背后的机制拆开配上能直接跑的代码和参数边界。2. 系统建模与二次模型化边界、模型分类与状态空间落地2.1 从系统边界到三类模型描述的映射关系定义一个系统时最先要确定的是边界边界以内是系统本体边界以外对系统的作用记为输入系统对边界以外的作用记为输出。这个定义看起来是术语题但它直接决定了后面建模时状态变量的个数输入输出一旦确定微分方程的阶次和传递函数的分子分母阶次就固定了。系统的三大要素是实体、属性和活动描述系统时还会用到事件这一术语。按属性分系统分工程系统和非工程系统按信号特征分分连续系统、离散系统、采样数据系统和离散-连续系统。采样控制系统就是典型的混合系统它既有连续信号又有离散信号按采样周期 T 重复工作这也是后面离散相似法必须满足采样定理的原因。模型层面有两级划分。按表达形式分物理模型和数学模型按数学表达形式又把数学模型分成静态模型和动态模型静态模型的表达形式是代数方程和逻辑关系表达式动态模型是微分方程和差分方程。再按变量间函数关系分线性模型和非线性模型。这里有个容易混淆的点——校核和验证不是一回事校核是检验数字仿真模型和数学模型是否一致验证是检验数字仿真模型和实际系统是否一致。前者查的是转换有没有错后者查的是模型本身对不对。2.2 一次模型化与二次模型化的衔接把实际系统抽象为数学模型叫一次模型化把数学模型转化为可在计算机上运行的仿真模型叫二次模型化。这两步之间的桥对连续系统来说通常是离散化。以状态方程 $\dot{x} Ax Bu$、$y Cx$ 为例二次模型化的目标是拿到 $x(k1) \Phi x(k) \Gamma u(k)$ 形式。常见做法是先求状态转移矩阵 $\Phi e^{AT}$再用 $\Gamma \int_0^T e^{A\tau} B , d\tau$ 求输入矩阵。MATLAB 里可以直接用expm或c2d完成这一步。% 连续系统状态空间模型 - 离散化仿真模型 A [0 1; 0 -2]; % 状态矩阵特征值 0 和 -2 B [0; 1]; % 输入矩阵 C [1 0]; % 输出矩阵 D 0; sys_c ss(A, B, C, D); T 0.1; % 仿真步长/虚拟采样周期 sys_d_zoh c2d(sys_c, T, zoh); % 零阶保持器离散化 sys_d_tust c2d(sys_c, T, tustin);% 双线性变换离散化 disp(sys_d_zoh.A); % 即 Phi disp(sys_d_zoh.B); % 即 Gammac2d的第一个参数是连续系统对象第二个是采样周期第三个 method 决定离散化方式。zoh对应加虚拟采样开关和零阶保持器是最贴近实际采样控制系统物理过程的做法tustin走双线性替换matched走零极点根匹配。这里要提醒一点连续系统仿真中计算速度和计算精度是一对固有矛盾。步长 h 越小截断误差越小但每步计算量固定总耗时线性上升步长放大能提速代价是稳定性边界被逼近。2.3 保持器与采样定理的约束保持器是把离散时间信号恢复成连续信号的装置。零阶保持器能较好地再现阶跃信号一阶保持器能较好地再现斜坡信号这是选型时最直接的判据。离散相似法在采样周期上必须满足采样定理即采样频率不低于信号最高频率分量的两倍。工程上一般取 5 到 10 倍留裕量因为保持器本身会引入相位滞后零阶保持器在频率 $\omega$ 处引入约 $\omega T/2$ 的相位滞后这个滞后会吃掉闭环系统的相位裕度裕量不够就直接振荡。仿真步长选完之后一定要回头检查相位裕度损失别只看时域波形看起来还行。3. 数值积分算法龙格-库塔法步长边界与稳定性判定3.1 单步法、多步法与显隐式的分类依据数值积分算法有三组常见分类维度。第一组按计算稳定性能否对步长 h 提出限制分条件稳定算法和绝对稳定算法第二组按本次计算用到前一次还是更前面的多次结果分单步法和多步法——常见的 RK 法是单步法Adams 法是多步法第三组按本次计算用到的数据是否全部已知分显式算法和隐式算法。三阶隐式 Adams 算法的截断误差是 $O(h^4)$。龙格-库塔法的基本思想是用几个点上函数值的线性组合来避免计算函数的高阶导数同时把精度做上去。二阶 RK 法的局部截断误差是 $O(h^3)$四阶 RK 法是 $O(h^5)$。3.2 二阶 RK 法解 y-y 的步长上限推导以 $y -y$即 $\tau 1$为例二阶 RK 迭代式为$$y_{k1} \left(1 - h \frac{h^2}{2}\right) y_k$$令放大因子 $R(h) 1 - h h^2/2$数值稳定的条件是 $|R(h)| 1$。解这个不等式得到 $0 h 2$。也就是步长必须小于两倍时间常数否则离散递推会发散。import numpy as np def rk2_solve(h, t_end10.0, y01.0): 二阶RK法解 y-y返回时间序列和数值解 t, y 0.0, y0 ts, ys [t], [y] while t t_end - 1e-12: k1 -y k2 -(y h * k1) # 二级RK的中间斜率 y y h * (k1 k2) / 2.0 # 加权平均更新 t h ts.append(t); ys.append(y) return np.array(ts), np.array(ys) for h in (0.5, 1.5, 1.9, 2.1): _, ys rk2_solve(h) print(fh{h:4} 末值{ys[-1]:12.4e} {发散 if abs(ys[-1])1e3 else 收敛})这段代码的核心是k2 -(y h*k1)它表示先用 k1 向前试走一步再取斜率。放大因子随 h 增长先减后增h 接近 2 时 $R(h)$ 趋于 1数值解出现持续振荡而不衰减h 超过 2 后 $|R(h)|1$直接发散。跑一遍上面的输出就能看到 1.9 收敛但振荡、2.1 发散的边界现象。3.3 四阶 RK 法与欧拉法的精度差对照经典四阶 RK 公式是$$y_{n1} y_n \frac{h}{6}(k_1 2k_2 2k_3 k_4)$$四个斜率分别取自区间起点、中点、中点、终点。以 $y t y$、$y(0)1$、$h0.1$ 为例精确解是 $y(t) 2e^t - t - 1$。def euler_and_rk4(f, h, y0): 同一初值下欧拉法与RK4法的单步结果对比 k1 f(0.0, y0) y_euler y0 h * k1 k2 f(0.5*h, y0 0.5*h*k1) k3 f(0.5*h, y0 0.5*h*k2) k4 f(h, y0 h*k3) y_rk4 y0 h*(k1 2*k2 2*k3 k4)/6 return y_euler, y_rk4 f lambda t, y: t y ye, yr euler_and_rk4(f, 0.1, 1.0) print(f欧拉法 y(0.1){ye:.4f}) # 1.1000 print(fRK4法 y(0.1){yr:.4f}) # 1.1103 print(f精确解 y(0.1){2*np.exp(0.1)-0.1-1:.4f})欧拉法只保留一阶项误差是 $O(h^2)$ 量级RK4 用四次函数求值换到 $O(h^5)$ 局部截断误差同样的 h 下精度高一个数量级以上。两者差异的根源就在泰勒展开保留的项数——欧拉法扔掉的高阶项RK4 用中点斜率把它们近似补回来了。另一处常考误差是 $y (y - k)/T$ 形式的系统。欧拉法递推为 $y_{m1} (1 - h/T)y_m (h/T)k$RK2 递推为带二次项的版本。系统特征值为 $-1/T$当 $h 2T$ 时数值解发散——判断依据不是看方程本身而是看 $|h\lambda|$ 是否超出该算法的绝对稳定域。3.4 算法稳定性对照表算法类型稳定条件说明欧拉法单步、显式$h 2/\lambda二阶 RK单步、显式$0 h 2/\lambda四阶 RK单步、显式约 $h 2.78/\lambda三阶隐式 Adams多步、隐式更大稳定域截断误差 $O(h^4)$双线性替换法离散化无条件稳定恒映射左半 s 平面到单位圆内根匹配法离散化无条件稳定极零点逐一映射这张表里最值得记住的结论是常见数值积分法几乎都是条件稳定的而双线性替换法和根匹配法是无条件稳定的。这也是快速数字仿真算法增广矩阵法、时域矩阵法、替换法、根匹配法被独立拎出来讲的原因——它们要求每步计算量小且具有良好的计算稳定性二者缺一不可。4. 离散化实战c2d 三种 method、双线性变换与根匹配法4.1 双线性变换把左半 s 平面映射到哪双线性变换的定义式是$$s \frac{2}{T}\cdot\frac{z-1}{z1}$$反过来 $z \dfrac{1 sT/2}{1 - sT/2}$。令 $s \sigma j\omega$ 代入并取模可得$$|z|^2 \frac{(2/T \sigma)^2 \omega^2}{(2/T - \sigma)^2 \omega^2}$$只要 $\sigma 0$即 s 在左半平面分子恒小于分母$|z| 1$即 z 落在单位圆内。结论双线性变换把左半 s 平面映射到 z 平面的单位圆内。这个映射的代价是频率畸变warpings 域虚轴上 $\omega$ 趋于无穷时z 域相位只到 $\pi$低频段近似线性、高频段被压缩。所以双线性变换适合低通和带通设计做宽带系统时要先做预畸变修正。4.2 c2d 的 method 参数与传递函数离散化c2d(sys, T, method)里 method 的取值直接决定离散化路径。当 method 为tustin时用的是双线性变换法当 method 为matched时用的是零极点根匹配法。% 双线性变换离散化 G(s)2/(s^22s1) G tf(2, [1 2 1]); T 0.1; Gz_tustin c2d(G, T, tustin); Gz_zoh c2d(G, T, zoh); Gz_matched c2d(G, T, matched); [nt, dt] tfdata(Gz_tustin, v); % 提取分子分母系数 fprintf(差分方程: y(k) -%f*y(k-1) - %f*y(k-2) %f*u(k) ...\n, ... dt(2), dt(3), nt(1));提取出 G(z) 的分子分母系数之后差分方程就是直接的系数搬运分子系数对应 $u(k), u(k-1), \dots$分母系数对应 $y(k-1), y(k-2), \dots$注意归一化到 $y(k)$ 系数为 1。例如数字校正环节 $D(z) \dfrac{z^2 - 0.04}{z^2 - 0.3z}$ 在采样周期 T0.02 秒下的仿真模型是 $y_n 0.3y_{n-1} u_n - 0.04u_{n-2}$这类题本质上就是对比系数。4.3 根匹配法与纯延迟环节的处理根匹配法的依据映射是 $z e^{Ts}$。操作步骤把 G(s) 的极点和零点逐个按 $z_i e^{s_i T}$ 映射到 z 平面增益按终值定理或稳态增益相等来配。如果 G(s) 的分母阶次 n 高于分子阶次 m那么在 G(z) 的分子上还需要配上 $n - m$ 个附加零点。这一步最容易被漏掉——G(s) 是严格真分式但 G(z) 分子分母阶次必须相同缺的零点要补上位置一般取在 z 平面原点附近或按特定规则放置。纯延迟环节 $G(s) e^{-0.45s}$ 在 T0.2 时延迟时间除以采样周期得 2.25。整数部分 C₀2小数部分 C₁0.25。整数部分直接用 z 的负幂次移位实现小数部分用线性插值在相邻两个采样点之间分配权重最终得到形如 $y_k 0.75u_{k-2} 0.25u_{k-3}$ 的仿真模型。# 纯延迟环节的线性插值仿真模型 def pure_delay(u_hist, delay_samples): u_hist: 输入历史序列最近在前; delay_samples: 延迟采样数(可为小数) c0 int(np.floor(delay_samples)) # 整数部分 c1 delay_samples - c0 # 小数部分 return (1 - c1) * u_hist[c0] c1 * u_hist[c0 1] buf [1.0, 2.0, 3.0, 4.0, 5.0] # u_k, u_{k-1}, u_{k-2}, ... print(pure_delay(buf, 2.25)) # 0.75*u_{k-2} 0.25*u_{k-3} 2.75delay_samples拆成整数与小数两部分是关键整数部分决定取哪两个历史值小数部分决定两者的插值权重。这个方法在通信仿真和过程控制里都很常用误差集中在延迟的分数部分。4.4 离散化结果的三个校验点拿到 G(z) 之后别急着写控制器先过三道校验极点在不在单位圆内稳定性、稳态增益和连续模型是否一致终值、分子分母阶次关系是否符合根匹配的补零规则。双线性变换会保持阶次关系所以 G(s) 分子一阶分母二阶时G(z) 分子分母都是二阶稳态增益方面G(s) 稳态增益为 0 则 G(z) 也为 0。根匹配法则要保持极点数量和位置对应附加零点单独确认。5. 参数优化与仿真步长自动控制目标函数与误差估计5.1 IAE、ISE、ITAE 这几类误差积分的选型差别控制系统参数优化设计中目标函数一般分两类加权性能指标型目标函数和误差积分型目标函数。后者常用的有误差绝对值积分IAE、误差平方积分ISE、时间乘误差绝对值积分ITAE、时间乘误差平方积分ITSE还有时间平方乘误差绝对值和时间平方乘误差平方两种扩展形式。选型上的经验ISE 对大误差惩罚重收敛快但超调可能偏大IAE 对大小误差权重均衡ITAE 因为乘了时间因子会压制长时间拖尾的残差阶跃响应更干净是工业整定里用得最多的一种。% 以二阶系统为例对比不同误差积分指标下的最优阻尼比 T_end 10; t 0:0.01:T_end; w0 1; for zeta [0.4 0.707 1.0] G tf(w0^2, [1 2*zeta*w0 w0^2]); [y, ~] step(G, t); e 1 - y; IAE trapz(t, abs(e)); ISE trapz(t, e.^2); ITAE trapz(t, t .* abs(e)); fprintf(zeta%.3f IAE%.4f ISE%.4f ITAE%.4f\n, zeta, IAE, ISE, ITAE); end代码用trapz做梯形数值积分e 1 - y是单位阶跃下的误差序列。改zeta就能看到三种指标给出的排序不一定一致——这正是目标函数选哪个必须结合具体性能要求来定的原因没有普适最优。参数优化问题也叫静态优化问题寻优途径分间接寻优法和直接寻优法两类。单纯形法属于直接寻优二维情况下正规单纯形是正三角形三维是正三棱锥它不需要求梯度靠反射、扩张、收缩几步迭代逼近最优点。5.2 目标函数梯度与单纯形法落地以 $Q(\alpha) \alpha_1^2 \frac{3}{2}\alpha_2^2 - \frac{1}{2}\alpha_1\alpha_2$ 为例梯度为$$\nabla Q \begin{bmatrix} 2\alpha_1 - \frac{1}{2}\alpha_2 \ 3\alpha_2 - \frac{1}{2}\alpha_1 \end{bmatrix}$$在 $\alpha_0 (2, 4)^T$ 处代入得梯度方向 $[1, 11]^T$ 方向的分量负梯度方向 $[-1, -11]^T$ 即为下降方向。from scipy.optimize import minimize def Q(a): a1, a2 a return a1**2 1.5*a2**2 - 0.5*a1*a2 res minimize(Q, x0[2.0, 4.0], methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-8}) print(res.x, res.fun) # 最优解接近 [0, 0]Q0methodNelder-Mead就是单纯形法的工程实现无需梯度信息适合目标函数不可导或求导代价高的场景。xatol/fatol控制单纯形的收敛判据设得过松会提前停在山脊上。5.3 步长自动控制的误差估计前提控制系统仿真过程中实现步长自动控制的前提是先做误差估计。常见做法是用两个不同阶次的算法同时算一步比较结果差异作为局部误差估计再按当前误差与容差的比值缩放步长。可变步长求解器如 MATLAB 的 ode45、ode15s内部就是这个逻辑。手动实现一个简单的步长自适应外层循环思路是每次用 h 和 h/2 各积分一步比较局部误差误差偏大就减半、偏小就放大。def adaptive_rk4(f, t0, y0, t_end, tol1e-5): 基于步长加倍法的自适应RK4误差超限则回退减半 t, y, h t0, y0, 0.1 while t t_end: y1 rk4_step(f, t, y, h) y2a rk4_step(f, t, y, h/2) y2 rk4_step(f, t h/2, y2a, h/2) err abs(y2 - y1) / 15.0 # RK4 误差估计因子 if err tol: t, y t h, y2 h min(h * 1.5, 0.5) # 误差达标则放大步长 else: h h / 2.0 # 误差超限则回退重算 return t, yerr |y2 - y1| / 15里的 15 来自 RK4 与两级半步法之间的误差阶数差是经典的经验因子。步长放大用 1.5 而不是 2是为了避免连续放大后立刻又超限造成抖动。生产环境里更稳妥的做法是直接用scipy.integrate.solve_ivp的rtol/atol参数让求解器自己管步长。参数作用典型取值rtol相对误差容差1e-3 ~ 1e-6atol绝对误差容差1e-6 ~ 1e-9max_step步长上限信号最小时间常数/10first_step初始步长由系统特征值估算采样控制系统的数字仿真一般方法有两种差分方程递推求解法和双重循环方法。前者适合控制器局部模型直接按差分方程逐拍递推后者按外层的采样周期循环和内层的仿真步长循环嵌套适合采样周期远大于仿真步长的场景。6. 仿真发散排查从特征值到步长的定位链路仿真跑出 NaN 或者数值幅度指数上升第一件事是算系统特征值而不是急着调小步长。定位链路的顺序是求连续模型特征值 $\lambda \sigma j\omega$用 $|h\lambda|$ 对照所选算法的绝对稳定域判断是算法本身发散还是模型本身就不稳定。几个常见坑的对照现象可能原因验证方法数值解指数增大h 超出条件稳定算法边界减半 h 观察是否收敛持续等幅振荡h 接近稳定边界放大因子接近 1打印 $R(h)$ 或 $双线性变换后频率响应偏移高频段频率畸变对比离散前后 bode 相位根匹配后直流增益不对缺少附加零点或增益未配用终值定理核对 G(1)Simulink 模型跑不动代数环或零交叉过多加单位延迟打断代数环Simulink 是 MATLAB 下的数字仿真工具文件类型是 .mdl用鼠标画出系统框图的方式建模主要用于动态系统建模。离散化模型要真正验证最直接的办法是把连续模型和离散模型放在同一个输入下逐拍对比。% 连续模型与离散模型在同一阶跃输入下的逐拍对比 sys_c tf(2, [1 2 1]); sys_d c2d(sys_c, 0.05, zoh); t 0:0.05:2; [y_c, ~] step(sys_c, t); [y_d, ~] step(sys_d, t); err max(abs(y_c - y_d)); fprintf(最大偏差 %.4f\n, err); % 步长越小偏差越小err是离散模型对连续模型的复现误差它随采样周期减小而单调下降。如果某次缩小步长后 err 反而变大说明离散化方法选错了——比如对含高频振荡的系统用了双线性变换却没做预畸变。这套先特征值、再步长、最后方法的排查顺序比无脑调参省时间。至于采样周期取值工程上我一般按闭环带宽的 10 到 20 倍来定再回头用上面这段误差对比确认精度够用。本文还有配套的精品资源点击获取 控制系统仿真：数值积分、离散化与步长稳定性解析
简介这份PDF是面向自动化、控制工程等专业学生的计算机控制仿真复习资料侧重选择题与计算题训练可用于课程期末、考研专业课或仿真实验前的查漏补缺。资源包仅含1个PDF文件约981KB内容以题库、概念辨析和计算题型为主便于打印或电子设备随时翻阅。已有43人学习下载。资料覆盖系统定义与分类、模型校核与验证、连续与离散系统仿真、保持器、数值积分及截断误差、龙格-库塔法、离散相似法、双线性变换、采样控制系统、MATLAB中c2d函数等核心考点并延伸至参数优化、目标函数与单纯形寻优等内容。读者可通过题目练习梳理仿真建模流程强化公式推导、稳定性判断和误差分析能力适合作为知识点自测与考前速查的配套材料。1. 控制系统仿真到底在算什么从一次数值发散说起很多人第一次用二阶龙格-库塔法解 $y -y$看到步长取 2.5 时结果直接炸掉才意识到仿真不是把方程丢给计算机这么简单。计算机控制仿真这道门槛本质上是三件事的叠加把物理系统抽象成数学模型一次模型化把数学模型转成计算机能跑的差分递推式二次模型化再在稳定性、精度、速度之间做妥协。这份题库覆盖的正是这三个环节的考点与算例——系统边界怎么划、保持器怎么选、龙格-库塔法步长如何限制、双线性变换为什么绝对稳定、根匹配法在纯延迟环节上怎么处理。它适合两类人一类是正在准备控制系统仿真、计算机仿真课程考试的学生需要把填空和计算题背后的推导链条补齐另一类是已经上手 MATLAB 做数字控制器实现的工程师想回头确认 c2d 的 tustin 和 matched 两种 method 到底差在哪、为什么离散模型会引入附加零点。下面不复述题库而是把每类题背后的机制拆开配上能直接跑的代码和参数边界。2. 系统建模与二次模型化边界、模型分类与状态空间落地2.1 从系统边界到三类模型描述的映射关系定义一个系统时最先要确定的是边界边界以内是系统本体边界以外对系统的作用记为输入系统对边界以外的作用记为输出。这个定义看起来是术语题但它直接决定了后面建模时状态变量的个数输入输出一旦确定微分方程的阶次和传递函数的分子分母阶次就固定了。系统的三大要素是实体、属性和活动描述系统时还会用到事件这一术语。按属性分系统分工程系统和非工程系统按信号特征分分连续系统、离散系统、采样数据系统和离散-连续系统。采样控制系统就是典型的混合系统它既有连续信号又有离散信号按采样周期 T 重复工作这也是后面离散相似法必须满足采样定理的原因。模型层面有两级划分。按表达形式分物理模型和数学模型按数学表达形式又把数学模型分成静态模型和动态模型静态模型的表达形式是代数方程和逻辑关系表达式动态模型是微分方程和差分方程。再按变量间函数关系分线性模型和非线性模型。这里有个容易混淆的点——校核和验证不是一回事校核是检验数字仿真模型和数学模型是否一致验证是检验数字仿真模型和实际系统是否一致。前者查的是转换有没有错后者查的是模型本身对不对。2.2 一次模型化与二次模型化的衔接把实际系统抽象为数学模型叫一次模型化把数学模型转化为可在计算机上运行的仿真模型叫二次模型化。这两步之间的桥对连续系统来说通常是离散化。以状态方程 $\dot{x} Ax Bu$、$y Cx$ 为例二次模型化的目标是拿到 $x(k1) \Phi x(k) \Gamma u(k)$ 形式。常见做法是先求状态转移矩阵 $\Phi e^{AT}$再用 $\Gamma \int_0^T e^{A\tau} B , d\tau$ 求输入矩阵。MATLAB 里可以直接用expm或c2d完成这一步。% 连续系统状态空间模型 - 离散化仿真模型 A [0 1; 0 -2]; % 状态矩阵特征值 0 和 -2 B [0; 1]; % 输入矩阵 C [1 0]; % 输出矩阵 D 0; sys_c ss(A, B, C, D); T 0.1; % 仿真步长/虚拟采样周期 sys_d_zoh c2d(sys_c, T, zoh); % 零阶保持器离散化 sys_d_tust c2d(sys_c, T, tustin);% 双线性变换离散化 disp(sys_d_zoh.A); % 即 Phi disp(sys_d_zoh.B); % 即 Gammac2d的第一个参数是连续系统对象第二个是采样周期第三个 method 决定离散化方式。zoh对应加虚拟采样开关和零阶保持器是最贴近实际采样控制系统物理过程的做法tustin走双线性替换matched走零极点根匹配。这里要提醒一点连续系统仿真中计算速度和计算精度是一对固有矛盾。步长 h 越小截断误差越小但每步计算量固定总耗时线性上升步长放大能提速代价是稳定性边界被逼近。2.3 保持器与采样定理的约束保持器是把离散时间信号恢复成连续信号的装置。零阶保持器能较好地再现阶跃信号一阶保持器能较好地再现斜坡信号这是选型时最直接的判据。离散相似法在采样周期上必须满足采样定理即采样频率不低于信号最高频率分量的两倍。工程上一般取 5 到 10 倍留裕量因为保持器本身会引入相位滞后零阶保持器在频率 $\omega$ 处引入约 $\omega T/2$ 的相位滞后这个滞后会吃掉闭环系统的相位裕度裕量不够就直接振荡。仿真步长选完之后一定要回头检查相位裕度损失别只看时域波形看起来还行。3. 数值积分算法龙格-库塔法步长边界与稳定性判定3.1 单步法、多步法与显隐式的分类依据数值积分算法有三组常见分类维度。第一组按计算稳定性能否对步长 h 提出限制分条件稳定算法和绝对稳定算法第二组按本次计算用到前一次还是更前面的多次结果分单步法和多步法——常见的 RK 法是单步法Adams 法是多步法第三组按本次计算用到的数据是否全部已知分显式算法和隐式算法。三阶隐式 Adams 算法的截断误差是 $O(h^4)$。龙格-库塔法的基本思想是用几个点上函数值的线性组合来避免计算函数的高阶导数同时把精度做上去。二阶 RK 法的局部截断误差是 $O(h^3)$四阶 RK 法是 $O(h^5)$。3.2 二阶 RK 法解 y-y 的步长上限推导以 $y -y$即 $\tau 1$为例二阶 RK 迭代式为$$y_{k1} \left(1 - h \frac{h^2}{2}\right) y_k$$令放大因子 $R(h) 1 - h h^2/2$数值稳定的条件是 $|R(h)| 1$。解这个不等式得到 $0 h 2$。也就是步长必须小于两倍时间常数否则离散递推会发散。import numpy as np def rk2_solve(h, t_end10.0, y01.0): 二阶RK法解 y-y返回时间序列和数值解 t, y 0.0, y0 ts, ys [t], [y] while t t_end - 1e-12: k1 -y k2 -(y h * k1) # 二级RK的中间斜率 y y h * (k1 k2) / 2.0 # 加权平均更新 t h ts.append(t); ys.append(y) return np.array(ts), np.array(ys) for h in (0.5, 1.5, 1.9, 2.1): _, ys rk2_solve(h) print(fh{h:4} 末值{ys[-1]:12.4e} {发散 if abs(ys[-1])1e3 else 收敛})这段代码的核心是k2 -(y h*k1)它表示先用 k1 向前试走一步再取斜率。放大因子随 h 增长先减后增h 接近 2 时 $R(h)$ 趋于 1数值解出现持续振荡而不衰减h 超过 2 后 $|R(h)|1$直接发散。跑一遍上面的输出就能看到 1.9 收敛但振荡、2.1 发散的边界现象。3.3 四阶 RK 法与欧拉法的精度差对照经典四阶 RK 公式是$$y_{n1} y_n \frac{h}{6}(k_1 2k_2 2k_3 k_4)$$四个斜率分别取自区间起点、中点、中点、终点。以 $y t y$、$y(0)1$、$h0.1$ 为例精确解是 $y(t) 2e^t - t - 1$。def euler_and_rk4(f, h, y0): 同一初值下欧拉法与RK4法的单步结果对比 k1 f(0.0, y0) y_euler y0 h * k1 k2 f(0.5*h, y0 0.5*h*k1) k3 f(0.5*h, y0 0.5*h*k2) k4 f(h, y0 h*k3) y_rk4 y0 h*(k1 2*k2 2*k3 k4)/6 return y_euler, y_rk4 f lambda t, y: t y ye, yr euler_and_rk4(f, 0.1, 1.0) print(f欧拉法 y(0.1){ye:.4f}) # 1.1000 print(fRK4法 y(0.1){yr:.4f}) # 1.1103 print(f精确解 y(0.1){2*np.exp(0.1)-0.1-1:.4f})欧拉法只保留一阶项误差是 $O(h^2)$ 量级RK4 用四次函数求值换到 $O(h^5)$ 局部截断误差同样的 h 下精度高一个数量级以上。两者差异的根源就在泰勒展开保留的项数——欧拉法扔掉的高阶项RK4 用中点斜率把它们近似补回来了。另一处常考误差是 $y (y - k)/T$ 形式的系统。欧拉法递推为 $y_{m1} (1 - h/T)y_m (h/T)k$RK2 递推为带二次项的版本。系统特征值为 $-1/T$当 $h 2T$ 时数值解发散——判断依据不是看方程本身而是看 $|h\lambda|$ 是否超出该算法的绝对稳定域。3.4 算法稳定性对照表算法类型稳定条件说明欧拉法单步、显式$h 2/\lambda二阶 RK单步、显式$0 h 2/\lambda四阶 RK单步、显式约 $h 2.78/\lambda三阶隐式 Adams多步、隐式更大稳定域截断误差 $O(h^4)$双线性替换法离散化无条件稳定恒映射左半 s 平面到单位圆内根匹配法离散化无条件稳定极零点逐一映射这张表里最值得记住的结论是常见数值积分法几乎都是条件稳定的而双线性替换法和根匹配法是无条件稳定的。这也是快速数字仿真算法增广矩阵法、时域矩阵法、替换法、根匹配法被独立拎出来讲的原因——它们要求每步计算量小且具有良好的计算稳定性二者缺一不可。4. 离散化实战c2d 三种 method、双线性变换与根匹配法4.1 双线性变换把左半 s 平面映射到哪双线性变换的定义式是$$s \frac{2}{T}\cdot\frac{z-1}{z1}$$反过来 $z \dfrac{1 sT/2}{1 - sT/2}$。令 $s \sigma j\omega$ 代入并取模可得$$|z|^2 \frac{(2/T \sigma)^2 \omega^2}{(2/T - \sigma)^2 \omega^2}$$只要 $\sigma 0$即 s 在左半平面分子恒小于分母$|z| 1$即 z 落在单位圆内。结论双线性变换把左半 s 平面映射到 z 平面的单位圆内。这个映射的代价是频率畸变warpings 域虚轴上 $\omega$ 趋于无穷时z 域相位只到 $\pi$低频段近似线性、高频段被压缩。所以双线性变换适合低通和带通设计做宽带系统时要先做预畸变修正。4.2 c2d 的 method 参数与传递函数离散化c2d(sys, T, method)里 method 的取值直接决定离散化路径。当 method 为tustin时用的是双线性变换法当 method 为matched时用的是零极点根匹配法。% 双线性变换离散化 G(s)2/(s^22s1) G tf(2, [1 2 1]); T 0.1; Gz_tustin c2d(G, T, tustin); Gz_zoh c2d(G, T, zoh); Gz_matched c2d(G, T, matched); [nt, dt] tfdata(Gz_tustin, v); % 提取分子分母系数 fprintf(差分方程: y(k) -%f*y(k-1) - %f*y(k-2) %f*u(k) ...\n, ... dt(2), dt(3), nt(1));提取出 G(z) 的分子分母系数之后差分方程就是直接的系数搬运分子系数对应 $u(k), u(k-1), \dots$分母系数对应 $y(k-1), y(k-2), \dots$注意归一化到 $y(k)$ 系数为 1。例如数字校正环节 $D(z) \dfrac{z^2 - 0.04}{z^2 - 0.3z}$ 在采样周期 T0.02 秒下的仿真模型是 $y_n 0.3y_{n-1} u_n - 0.04u_{n-2}$这类题本质上就是对比系数。4.3 根匹配法与纯延迟环节的处理根匹配法的依据映射是 $z e^{Ts}$。操作步骤把 G(s) 的极点和零点逐个按 $z_i e^{s_i T}$ 映射到 z 平面增益按终值定理或稳态增益相等来配。如果 G(s) 的分母阶次 n 高于分子阶次 m那么在 G(z) 的分子上还需要配上 $n - m$ 个附加零点。这一步最容易被漏掉——G(s) 是严格真分式但 G(z) 分子分母阶次必须相同缺的零点要补上位置一般取在 z 平面原点附近或按特定规则放置。纯延迟环节 $G(s) e^{-0.45s}$ 在 T0.2 时延迟时间除以采样周期得 2.25。整数部分 C₀2小数部分 C₁0.25。整数部分直接用 z 的负幂次移位实现小数部分用线性插值在相邻两个采样点之间分配权重最终得到形如 $y_k 0.75u_{k-2} 0.25u_{k-3}$ 的仿真模型。# 纯延迟环节的线性插值仿真模型 def pure_delay(u_hist, delay_samples): u_hist: 输入历史序列最近在前; delay_samples: 延迟采样数(可为小数) c0 int(np.floor(delay_samples)) # 整数部分 c1 delay_samples - c0 # 小数部分 return (1 - c1) * u_hist[c0] c1 * u_hist[c0 1] buf [1.0, 2.0, 3.0, 4.0, 5.0] # u_k, u_{k-1}, u_{k-2}, ... print(pure_delay(buf, 2.25)) # 0.75*u_{k-2} 0.25*u_{k-3} 2.75delay_samples拆成整数与小数两部分是关键整数部分决定取哪两个历史值小数部分决定两者的插值权重。这个方法在通信仿真和过程控制里都很常用误差集中在延迟的分数部分。4.4 离散化结果的三个校验点拿到 G(z) 之后别急着写控制器先过三道校验极点在不在单位圆内稳定性、稳态增益和连续模型是否一致终值、分子分母阶次关系是否符合根匹配的补零规则。双线性变换会保持阶次关系所以 G(s) 分子一阶分母二阶时G(z) 分子分母都是二阶稳态增益方面G(s) 稳态增益为 0 则 G(z) 也为 0。根匹配法则要保持极点数量和位置对应附加零点单独确认。5. 参数优化与仿真步长自动控制目标函数与误差估计5.1 IAE、ISE、ITAE 这几类误差积分的选型差别控制系统参数优化设计中目标函数一般分两类加权性能指标型目标函数和误差积分型目标函数。后者常用的有误差绝对值积分IAE、误差平方积分ISE、时间乘误差绝对值积分ITAE、时间乘误差平方积分ITSE还有时间平方乘误差绝对值和时间平方乘误差平方两种扩展形式。选型上的经验ISE 对大误差惩罚重收敛快但超调可能偏大IAE 对大小误差权重均衡ITAE 因为乘了时间因子会压制长时间拖尾的残差阶跃响应更干净是工业整定里用得最多的一种。% 以二阶系统为例对比不同误差积分指标下的最优阻尼比 T_end 10; t 0:0.01:T_end; w0 1; for zeta [0.4 0.707 1.0] G tf(w0^2, [1 2*zeta*w0 w0^2]); [y, ~] step(G, t); e 1 - y; IAE trapz(t, abs(e)); ISE trapz(t, e.^2); ITAE trapz(t, t .* abs(e)); fprintf(zeta%.3f IAE%.4f ISE%.4f ITAE%.4f\n, zeta, IAE, ISE, ITAE); end代码用trapz做梯形数值积分e 1 - y是单位阶跃下的误差序列。改zeta就能看到三种指标给出的排序不一定一致——这正是目标函数选哪个必须结合具体性能要求来定的原因没有普适最优。参数优化问题也叫静态优化问题寻优途径分间接寻优法和直接寻优法两类。单纯形法属于直接寻优二维情况下正规单纯形是正三角形三维是正三棱锥它不需要求梯度靠反射、扩张、收缩几步迭代逼近最优点。5.2 目标函数梯度与单纯形法落地以 $Q(\alpha) \alpha_1^2 \frac{3}{2}\alpha_2^2 - \frac{1}{2}\alpha_1\alpha_2$ 为例梯度为$$\nabla Q \begin{bmatrix} 2\alpha_1 - \frac{1}{2}\alpha_2 \ 3\alpha_2 - \frac{1}{2}\alpha_1 \end{bmatrix}$$在 $\alpha_0 (2, 4)^T$ 处代入得梯度方向 $[1, 11]^T$ 方向的分量负梯度方向 $[-1, -11]^T$ 即为下降方向。from scipy.optimize import minimize def Q(a): a1, a2 a return a1**2 1.5*a2**2 - 0.5*a1*a2 res minimize(Q, x0[2.0, 4.0], methodNelder-Mead, options{xatol: 1e-6, fatol: 1e-8}) print(res.x, res.fun) # 最优解接近 [0, 0]Q0methodNelder-Mead就是单纯形法的工程实现无需梯度信息适合目标函数不可导或求导代价高的场景。xatol/fatol控制单纯形的收敛判据设得过松会提前停在山脊上。5.3 步长自动控制的误差估计前提控制系统仿真过程中实现步长自动控制的前提是先做误差估计。常见做法是用两个不同阶次的算法同时算一步比较结果差异作为局部误差估计再按当前误差与容差的比值缩放步长。可变步长求解器如 MATLAB 的 ode45、ode15s内部就是这个逻辑。手动实现一个简单的步长自适应外层循环思路是每次用 h 和 h/2 各积分一步比较局部误差误差偏大就减半、偏小就放大。def adaptive_rk4(f, t0, y0, t_end, tol1e-5): 基于步长加倍法的自适应RK4误差超限则回退减半 t, y, h t0, y0, 0.1 while t t_end: y1 rk4_step(f, t, y, h) y2a rk4_step(f, t, y, h/2) y2 rk4_step(f, t h/2, y2a, h/2) err abs(y2 - y1) / 15.0 # RK4 误差估计因子 if err tol: t, y t h, y2 h min(h * 1.5, 0.5) # 误差达标则放大步长 else: h h / 2.0 # 误差超限则回退重算 return t, yerr |y2 - y1| / 15里的 15 来自 RK4 与两级半步法之间的误差阶数差是经典的经验因子。步长放大用 1.5 而不是 2是为了避免连续放大后立刻又超限造成抖动。生产环境里更稳妥的做法是直接用scipy.integrate.solve_ivp的rtol/atol参数让求解器自己管步长。参数作用典型取值rtol相对误差容差1e-3 ~ 1e-6atol绝对误差容差1e-6 ~ 1e-9max_step步长上限信号最小时间常数/10first_step初始步长由系统特征值估算采样控制系统的数字仿真一般方法有两种差分方程递推求解法和双重循环方法。前者适合控制器局部模型直接按差分方程逐拍递推后者按外层的采样周期循环和内层的仿真步长循环嵌套适合采样周期远大于仿真步长的场景。6. 仿真发散排查从特征值到步长的定位链路仿真跑出 NaN 或者数值幅度指数上升第一件事是算系统特征值而不是急着调小步长。定位链路的顺序是求连续模型特征值 $\lambda \sigma j\omega$用 $|h\lambda|$ 对照所选算法的绝对稳定域判断是算法本身发散还是模型本身就不稳定。几个常见坑的对照现象可能原因验证方法数值解指数增大h 超出条件稳定算法边界减半 h 观察是否收敛持续等幅振荡h 接近稳定边界放大因子接近 1打印 $R(h)$ 或 $双线性变换后频率响应偏移高频段频率畸变对比离散前后 bode 相位根匹配后直流增益不对缺少附加零点或增益未配用终值定理核对 G(1)Simulink 模型跑不动代数环或零交叉过多加单位延迟打断代数环Simulink 是 MATLAB 下的数字仿真工具文件类型是 .mdl用鼠标画出系统框图的方式建模主要用于动态系统建模。离散化模型要真正验证最直接的办法是把连续模型和离散模型放在同一个输入下逐拍对比。% 连续模型与离散模型在同一阶跃输入下的逐拍对比 sys_c tf(2, [1 2 1]); sys_d c2d(sys_c, 0.05, zoh); t 0:0.05:2; [y_c, ~] step(sys_c, t); [y_d, ~] step(sys_d, t); err max(abs(y_c - y_d)); fprintf(最大偏差 %.4f\n, err); % 步长越小偏差越小err是离散模型对连续模型的复现误差它随采样周期减小而单调下降。如果某次缩小步长后 err 反而变大说明离散化方法选错了——比如对含高频振荡的系统用了双线性变换却没做预畸变。这套先特征值、再步长、最后方法的排查顺序比无脑调参省时间。至于采样周期取值工程上我一般按闭环带宽的 10 到 20 倍来定再回头用上面这段误差对比确认精度够用。本文还有配套的精品资源点击获取