用MATLAB搞懂计算机控制：离散化、数字PID与状态空间实战

《计算机控制理论与应用》这门课，大概是自动化及相关专业里"公式占比最高、直观感最低"的一门课。很多同学学到z变换和差分方程时，还是习惯用连续系统的直觉去推离散系统，结果上课点头、做题摇头、MATLAB仿真直接懵。我见过不少做课程设计的朋友，手算传递函数头头是道，一到MATLAB里画阶跃响应就发现曲线方向不对、稳定性判断对不上。其实问题往往不在MATLAB，而在"连续"和"离散"的思维切换没有完成。这篇文章想做的，就是把这门课里最容易被忽略、又最值得用MATLAB验证的几块内容展开讲：离散化怎么用代码验证、数字PID调参的仿真思路、状态空间设计从Bode图思维切换到状态反馈思维，以及仿真跑通之后还要补哪些坑。适合正在学这门课、准备课程设计，或者复习时想自己写一遍代码的人。

1. 学计算机控制理论为什么总卡在"看得懂但算不对"

1.1 连续域那套直觉，在离散域会反过来坑你

连续系统的分析方法，很多人脑海中已经形成了一套比较固定的路径：拉普拉斯变换、传递函数、Bode图、稳定判据，靠这些工具可以处理绝大多数线性时不变系统。离散系统看上去公式长得差不多：z变换、脉冲传递函数、单位圆判稳，也是一条平行路径。但两者一旦混在一起，问题就出现了。最常见的第一类错误，是把连续传递函数里的s直接换成z，比如典型被控对象G(s) = 1 / (s(s+2))，有同学会写成G(z) = 1 / (z(z+2))，然后用zplane判断稳定性，这一步的物理意义就完全错了。连续域的s和离散域的z之间不是替换关系，而是经过"采样 + 保持 + 离散化"的映射关系，只有通过G(z) = Z[G_h(s)G(s)]这个环节，才能得到真正的脉冲传递函数。

第二类常见错误，是忽略终值定理和初值定理的使用条件。离散系统的终值定理要求(1 - z^(-1))F(z)在单位圆上以及单位圆外没有极点，而很多同学连续系统用习惯了，直接套用，算出来的稳态误差经常和仿真对不上。这些错误在MATLAB里暴露得非常快：只要把手算结果写成代码，再和c2d的答案对比，立刻就知道哪一步想错了。我带课程设计的时候，看到最多的问题就是"手算结果和仿真曲线稳态值差几个数量级"，最后排查下来基本都是z变换表查错、或者保持器环节漏加。这类问题如果只靠手算，可能反复检查两小时都发现不了，但用MATLAB对照，定位就是几秒钟的事。

1.2 验证手算、搭实验平台、快速探索参数：MATLAB的三个用法

MATLAB在这门课里的用途，我大致归成三类。

第一类是验证手算。教科书上的课后习题大量涉及z变换、差分方程、脉冲传递函数化简，手算完不确定对错，用符号计算和c2d核对一下，效率极高。这里的建议是先手算、再验证，不要反过来。否则很容易形成"会用代码但不懂原理"的假熟练，考试时还是要吃亏。

第二类是搭实验平台。simulink可以用离散模块搭出完整的采样控制系统，也可以直接在M脚本里用for循环按拍迭代。M脚本的好处是每一步都透明，能清清楚楚看到u(k)、y(k)、e(k)的变化，对理解"拍"这个概念特别有帮助。很多同学对"系统是一拍一拍运行的"没有体感，直到自己在循环里写出y(k+1)的递推式，才真正明白上一拍的控制量如何影响下一拍的输出。

第三类是快速探索参数对性能的影响。把采样周期、PID参数、极点位置写成变量，跑一个参数扫描，画一组曲线对比，就能直观感受到课本里"采样周期越小，数字控制越接近连续控制"这句话到底意味着什么。这三个用途也对应三个层次：会验证、会搭平台、会探索。把这三件事做完，这门课的很多结论就不再是背下来的，而是自己"看见"的。

2. 离散化这一段：用MATLAB把"近似"这件事讲透

2.1 z变换算不准的常见原因：用c2d一行代码对答案

先说离散化。书上的公式一般长这样：G(z) = Z[G_h(s)G(s)]，其中G_h(s) = (1 - e^(-sT)) / s是零阶保持器的传递函数。手动做这个z变换并不轻松，要用到z变换表和部分分式展开，一个符号不小心就前功尽弃。MATLAB里只需要一行：

s = tf('s'); G = 1/(s*(s+2)); % 连续被控对象 Ts = 0.1; % 采样周期 Gd = c2d(G, Ts, 'zoh'); % 零阶保持器离散化

得到的结果是一个z传递函数，分子分母的系数可以直接用[numd, dend] = tfdata(Gd, 'v');取出来，再和手算的结果逐项比较。如果手算对了，这段代码几秒钟就能给出确认；如果手算错了，也能通过对比快速定位是分子错还是分母错、是变换表用错还是保持器环节漏掉。

这里特别解释一下c2d的第三个参数。'zoh'、'foh'、'tustin'分别对应零阶保持器、一阶保持器、双线性变换。三者的差异，说白了是"用什么形状去逼近连续信号"：零阶保持器认为一个采样周期内信号保持不变，像一个台阶；一阶保持器用斜线连接相邻两点；双线性变换是在频率域做一种代数近似。实际课程设计里选'zoh'就好，因为和DAC输出行为最接近，也最符合课本定义。它们的特点可以整理成一张表：

2.2 采样周期是离散系统最大的设计变量

离散系统和连续系统最大的区别，就是凭空多出一个采样周期T。T选多大，直接影响系统的稳定性、动态性能和可控性。MATLAB非常适合做这件事：把同一个连续被控对象分别用T=0.01、T=0.1、T=1离散化，画三条阶跃响应曲线，对比一目了然。

s = tf('s'); G = 1/(s*(s+2)); for Ts = [0.01 0.1 1] Gd = c2d(G, Ts, 'zoh'); step(Gd); hold on; end legend({'T=0.01', 'T=0.1', 'T=1'});

实际跑下来你会发现，T=0.01时曲线跟连续系统几乎重合；T=0.1时曲线开始出现可察觉的延迟感；T=1时可能已经变了味道，稳定裕度明显下降。这说明采样周期的选取有实际约束，不是越小越好，太小计算量增大、数值敏感性增强；也不是越大越省，性能会劣化。工程上常用的经验是取系统上升时间的1/10到1/20作为初始值，再根据仿真结果微调。这个"经验值"在MATLAB里可以被量化验证，比死记数字有用得多。以G(s)=1/(s(s+2))为例，上升时间大约1秒左右，那么T=0.1就落在合理区间，T=0.01偏保守，T=1就明显偏大了。

2.3 零阶保持器不是玄学，是一个sinc

很多同学不理解为什么离散化要在系统前面加一个保持器，觉得那是书本上的形式主义。其实零阶保持器的物理含义很朴素：计算机算出一个控制量u(k)后，DAC要把数字量变成模拟量送给执行机构，在下一个采样时刻到来之前，电平保持不变。这个"阶梯波形"就是零阶保持器的输出。它的传递函数H_0(s) = (1 - e^(-sT)) / s，从频域看，幅频特性是|sinc(ωT/2)|的形状，本质上是一个带低通性质的重构滤波器。

在MATLAB里可以直接画离散化后系统的Bode图，直观看到高频段增益衰减和相位延迟：

Ts = 0.1; Gd = c2d(1/(s*(s+2)), Ts, 'zoh'); bode(Gd); grid on;

理解了这一点，再回头看课本推导就不会觉得是硬背公式了。离散化就是从"脉冲序列"到"阶梯信号"再到"被控对象"的完整链路，每一个环节都有对应的传递函数。如果把保持器环节漏掉，你就等于默认DAC输出是一串理想脉冲，这在实际系统里是不存在的。

3. 数字PID：MATLAB里调参，等于把试错成本压到最低

3.1 位置式PID和增量式PID：代码差别只有五行，行为差很多

PID是这门课里最接近工程的内容，也是课程设计里最常被要求实现的控制器。MATLAB仿真的价值在于，可以把试错成本从实验台上的反复烧程序，压到"改一个变量重新跑一下"。

位置式PID的差分方程是u(k) = Kp·e(k) + Ki·T·Σe(j) + Kd·(e(k) - e(k-1)) / T，增量式PID输出的是控制量增量Δu(k) = u(k) - u(k-1)，系数算出来的是差值。两者在无扰动、无饱和的理想情况下，性能几乎一样，但面对执行器饱和时行为差异很大。增量式的输出天然不累积历史误差，饱和时自动停止增长，所以工程上更常用；位置式则一定需要额外的抗积分饱和逻辑。

下面这段代码是我常用的位置式PID仿真模板，为了把循环写得短一点，被控对象选了一阶惯性环节1/(s+1)的零阶保持离散模型，二阶对象的逻辑完全一样，只是递推式多一两个历史项：

% 位置式PID仿真，一阶对象G(s)=1/(s+1)，零阶保持离散化 Ts = 0.01; t = 0:Ts:2; r = ones(size(t)); % 单位阶跃给定 a = exp(-Ts); b = 1 - exp(-Ts); % 离散模型 y(k+1)=a*y(k)+b*u(k) Kp = 1.5; Ki = 3.0; Kd = 0.05; u = zeros(size(t)); y = zeros(size(t)); e_int = 0; e_prev = 0; for k = 1:length(t)-1 e = r(k) - y(k); e_int = e_int + e*Ts; u(k+1) = Kp*e + Ki*e_int + Kd*(e - e_prev)/Ts; e_prev = e; y(k+1) = a*y(k) + b*u(k+1); end

增量式的核心改动只有几行，把u(k+1)换成Δu，再累加：

du = Kp*(e - e_prev) + Ki*e*Ts + Kd*(e - 2*e_prev + e_prev2)/Ts; u(k+1) = u(k) + du;

这里e_prev2是再前一次的误差，需要多存一个变量。

3.2 积分饱和与微分噪声：仿真里一眼就能看见的隐患

很多教材会在末尾提一句"PID参数过大容易饱和"，但只有仿真才能让你真正"看见"积分饱和有多严重。把前面的代码加上限幅u = max(min(u, 1), -1);，然后刻意把Ki调大，你会发现阶跃响应出现一个很大的超调，而且回稳很慢。原因很直观：误差持续为正时，积分项不停累积，就算实际输出已经顶到执行器上限，积分项还在暗地里增加。等误差反向时，积分项需要很长时间才能"消化"掉积累的值，超调就被严重放大。

微分项的问题也类似。实际系统的测量信号总是带噪声的，微分算子对噪声极其敏感，一个幅值只有0.01的随机噪声，差分后可能被放大成幅值几倍的波动。在仿真里加一行：

y_meas = y + 0.01*randn(size(y)); e = r - y_meas; % 用带噪声的测量值计算误差

再看控制量曲线，就知道为什么工业上要给微分项加一阶低通滤波器，或者干脆只用PI。仿真能帮你提前建立这种"敬畏感"，在实验课上往往是要踩了坑才能体会到的事。

3.3 一个"先比例后积分再微分"的调参顺序

如果只是想让仿真尽快跑通，调参顺序可以固定成一套路子。先把Ki和Kd设成0，只留Kp，从小到大试Kp，直到阶跃响应出现轻微振荡；然后加入Ki，用小一点的Ki慢慢消除稳态误差，注意观察超调量；最后加入Kd，用很小的Kd抑制振荡，超调减小后再微调Kp。这个过程在M脚本里就是一个for循环套一个参数数组，每次画一条曲线，比在simulink里手动拖滑块更可控。

也可以用MATLAB自带的pidtune做自动化设计，拿到初值后再手工微调。但经验是先相信仿真曲线，再去理解参数为什么这样取，而不是背"越大越灵敏"之类的口诀。同一套被控对象，换个采样周期，PID参数很可能就要全部重调，因为离散化后的模型变了，控制器的相位裕度、增益裕度都跟着变。这种情况在课堂上容易被忽略，但做一次仿真对比就有体感了。

4. 状态空间设计：从Bode图思维切换到状态反馈思维

4.1 极点配置不是"放几个极点"，而是重写特征方程

到了状态空间这部分，很多人的思维还停留在传递函数和Bode图上，看到矩阵就头疼。但实际上状态空间设计的逻辑比频域法更直白：想让系统快一点，就把闭环极点放到离原点更远的地方；想让系统稳一点，就把极点放到实轴左侧更远的位置。极点配置要做的，就是找一个状态反馈增益矩阵K，使得闭环系统矩阵A - BK的特征值等于你指定的极点。

MATLAB里这个操作被压缩成了一行：

A = [0 1; -2 -3]; B = [0; 1]; p = [-5 -6]; % 期望闭环极点 K = place(A, B, p); eig(A - B*K) % 验证闭环特征值

但先别急着背命令，要想清楚两个前提。第一个前提是系统必须能控，也就是可控性矩阵[B AB A^2B ...]必须满秩，否则有些极点根本"够不着"，place会直接报错。这个报错不是代码问题，而是系统的结构性限制，值得停下来想一想。第二个前提是极点的位置和控制器增益密切相关，极点放得越远，K的数值通常越大，意味着控制量越大、对噪声越敏感。这意味着"快"是有代价的。仿真中建议把期望极点依次取为[-2, -3]、[-5, -6]、[-10, -12]，画出控制量曲线，直观感受"快"和"用力猛"的关系。

4.2 状态观测器：另一个极点配置问题

真正常见的情况是，状态变量不能全部直接测量。这时候要设计状态观测器，用输入u和输出y重构出状态的估计值。观测器的核心是让估计误差的动态矩阵A - LC具有快速衰减的特征值，本质上就是第二个极点配置问题。

C = [1 0]; p_obs = [-15 -16]; % 观测器极点通常比控制器极点快3到5倍 L = place(A', C', p_obs)';

设计完K和L后，完整的闭环仿真才是重点。可以让真实系统从x0 = [1; 0]出发，观测器从x_hat0 = [0; 0]出发，状态反馈用的是估计状态而不是真实状态。把真实状态和估计状态画在同一个图上，你会看到估计值在几个采样周期内追上真实值，之后系统才真正体现"极点配置后的动态"。仿真闭环的核心片段大致长这样：

sys = ss(A, B, C, 0); sysd = c2d(sys, Ts, 'zoh'); [Ad, Bd, Cd, Dd] = ssdata(sysd); x = [1; 0]; xh = [0; 0]; for k = 1:length(t)-1 u_k = -K*xh; % 用估计状态做反馈 x = Ad*x + Bd*u_k; % 真实系统更新 y = Cd*x; xh = Ad*xh + Bd*u_k + L*(y - Cd*xh); % 观测器更新 end

这个过程非常直观地展示了分离原理：状态反馈和观测器可以分开设计，组合在一起后，整个闭环系统的极点近似等于控制器极点和观测器极点的并集。在仿真里把观测器初值故意设偏，再看估计误差收敛，比背十遍原理都管用。

4.3 LQR与极点配置：一个管位置，一个管代价

极点配置的问题在于，它只告诉你极点放哪里，不告诉你这个"快"要花多少控制能量。LQR则把性能指标写成J = ∫(x'Qx + u'Ru)dt，通过调整Q和R的比值，在"状态跟踪得准"和"控制量不浪费"之间找平衡。MATLAB里的lqr命令：

Q = [10 0; 0 1]; % 加大Q对角元，让对应状态更快收敛 R = 1; K_lqr = lqr(A, B, Q, R);

一个比较实用的体会是，极点配置适合教学和概念理解，它直白、可控、容易验证；LQR适合工程初设，因为权重矩阵反映了你的真实诉求——到底更在乎快速还是更在乎节能。两者不是对立的，LQR得到的闭环极点也可以反过来近似用place实现，但空间里多了一层"代价"的视角，这是状态空间方法相比频域法的核心优势。

5. 仿真通过不等于实验通过：那些MATLAB没告诉你的坑

5.1 反因果错误：仿真里最隐蔽、实验中最致命的坑

这是这门课最值得单独讲的一个坑。很多同学在仿真里写控制律u(k) = K·e(k)，用的e(k)是当前拍的误差，在MATLAB里这只是数组索引的问题，跑起来完全正常。但实际系统里，y(k)这一拍的数据往往是这一拍结束才测量并传到控制器的，控制器算完再输出，执行器再动作，整个链路已经引入了至少一拍延迟。如果仿真用的是"零拍延迟"模型，现场程序几乎必然会振荡。

正确的做法是在仿真里主动引入一拍延迟：让控制律使用y(k-1)或者u(k-1)对应的误差，把系统行为改成"测量-计算-执行"的真实时序。具体实现就是在迭代循环里把状态更新和控制量计算错开一拍。很多教材的simulink示例里，模型本身就有一个Unit Delay模块，很多同学不知道为什么要加，其实就是因为真实系统做不到"这一拍测量、同一拍输出控制量"。你在MATLAB里随手把延迟去掉，绝对是最容易犯、也最不容易被发现的错误。

5.2 连续+离散混合仿真的数值积分设置

当被控对象是连续的、控制器是离散的，simulink默认的变步长求解器有时会带来奇怪的结果——连续部分的积分步长和离散采样周期之间没有整数倍数关系，导致离散环节的触发时刻和积分步长对不上，曲线抖动甚至出现假的不稳定。建议把求解器设为固定步长，步长取采样周期的1/10到1/100，再运行一遍对比。

这个经验在M脚本里不容易遇到，但在simulink里很常见。很多人以为是控制器参数不对，调了半天Kp、Ki，最后发现只是数值积分设置的问题。检查顺序一般是：先固定步长，步长设为采样周期的整数分之一，再回看仿真曲线是否恢复平滑。如果连续部分本身有快动态（比如高频振荡模态），还需要把固定步长取得足够小，否则会引入数值阻尼或虚假振荡。

5.3 饱和、死区、纯延时：把硬件的脾气加进模型里

课本上的仿真几乎都是线性理想模型，但实验台上的执行器有行程限制，传感器有死区，通信和计算有纯延时。把这三样加进仿真，你会立刻看到性能下降甚至发散。饱和限幅的代码很简单：

u = max(min(u, u_max), u_min);

但把它放进仿真后，PID参数可能又要重新调一遍。死区可以建模为：

if abs(u) < delta u = 0; end

纯延时可以用m函数的延迟缓存，或者simulink的Transport Delay模块。我的建议是课程设计的仿真分两步走：先跑线性理想模型验证控制算法，再把饱和、死区和一拍延迟加进去，检查控制器在这种"更真实"的系统里还能不能稳住。这一步做完，才算真正把理论变成工程，也才能在答辩时讲清楚"为什么仿真结果和理论计算有差别"。

我自己这些年用MATLAB研究这门课的体会是，仿真的价值不只是帮你把作业做完，而是逼着你把每一个公式都变成可以交互的东西。如果你现在正在学这门课，建议挑一个具体的例子，从头到尾走一遍完整的链路——连续模型、离散化、数字PID或状态反馈、加饱和、加噪声、加一拍延迟，每一步都停下来和书本上的公式对照。这个过程跑完，你会发现自己对"计算机控制"的理解完全不同了。最后再分享一个很实用的小技巧：所有MATLAB脚本，一开始就把采样周期T写成一个变量，所有参数都用变量表达。后续扫描参数时只需要改一行，整个链路都会跟着重算，能省下大量来回改代码的时间。