简介面向自动化、控制工程方向的学习者与工程调试人员这份增量式PID控制算法MATLAB仿真文档围绕PID调节原理、增量式算法实现与仿真程序编写展开。文档以传递函数G(s)5/(s^22s10)为被控对象覆盖离散化、Z传递函数推导、差分方程建立以及单位阶跃、正弦信号输入下的PID程序实现采样时间设为1ms控制器输出限幅[-5,5]并给出系统输出与误差曲线绘制方法。压缩包内仅含1个docx文件大小约162KB以文字、公式与程序清单为主便于直接查阅与复现。作者还记录了kp、ki、kd由小到大逐项整定、观察临界振荡与微调增益的过程并附带限幅与选择条件对波形影响的讨论。该资源已有2088人学习下载适合作为课程实验、毕业设计或控制算法入门阶段的参考范例帮助读者理解增量式PID公式与MATLAB仿真代码之间的对应关系并掌握参数整定思路。1. 增量式 PID 的 MATLAB 仿真同一组 kp、ki、kd 为什么跑出两条曲线同一组 kp6、ki45、kd5脚本画出来的波形和 Simulink 画出来的波形对不上误差曲线一直压不下去限幅削顶之后更是各走各的。很多人第一次做增量式 PID 仿真就卡在这里然后开始怀疑算法、怀疑参数、怀疑自己的数学。实际上卡点通常不在算法本身而在三个极容易被跳过的细节被控对象的分母系数被 MATLAB 解析成了另一个阶数控制器参数在两套定义之间没有按采样周期折算差分方程只取了前两拍历史值导致模型被截断。这篇围绕 G(s)5/(s²2s10) 这个二阶弱阻尼对象把 c2d 离散化、增量式差分方程推导、主循环更新顺序、输出限幅与抗积分饱和、参数整定与发散排查一路走完。脚本可以直接复制运行参数怎么改、改完看哪条曲线、曲线不对时先查什么都落在具体代码上。适合正在做课程设计、准备把 PID 从 MATLAB 搬到嵌入式代码里的工程师。2. 增量式 PID 的差分方程推导与 Kp、Ki、Kd 折算2.1 位置式与增量式的区别到底在哪一行代码位置式 PID 把控制量写成误差的绝对表达u(k) Kp·e(k) Ki·T·Σe(j) (Kd/T)·[e(k)−e(k−1)]。求和项要一直保存量程会随时长增长输出一旦被限幅削平积分仍在内部继续累积退饱和时就会有一段明显过冲这就是经典的积分饱和windup。手自动切换时位置式还需要把输出值强行搬到积分累加器里否则切换瞬间会跳变。增量式不再算绝对量而是算这一拍该加多少。把 u(k) 和 u(k−1) 两个式子相减求和符号正好抵消Δ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)三项分别对应比例、积分、微分的增量贡献这也正是原程序里 x(1)、x(2)、x(3) 三行的来历。要注意增量式并不是不累积累积被搬到了输出侧u(k)u(k−1)Δu(k)。限幅只改 u不改内部的历史状态所以增量式天然比位置式更容易做抗饱和。工程上还有一种做法是只在未饱和时更新 u(k−1)效果更干净5.2 节会给出代码。2.2 把 G(s)5/(s²2s10) 离散成可迭代的差分方程2.2.1 c2d 与 tfdata 的正确用法MATLAB 提供 c2d 直接把连续传递函数转成 Z 传递函数默认方法是零阶保持。老代码里常见 z 这个写法新版本建议显式写 zoh避免不同版本下的默认行为差异ts 0.001; % 采样周期 1 ms sys tf(5,[1,2,10]); % 被控对象 G(s)5/(s^22s10) dsys c2d(sys,ts,zoh); % 零阶保持离散化 [num,den] tfdata(dsys,v); % v 返回行向量而不是 cell fprintf(num %s\n, mat2str(num,6)); fprintf(den %s\n, mat2str(den,6));tfdata 第二个参数给 v返回的是数值行向量写循环时比 cell 取值方便得多。den 的首元素恒为 1说明返回的是 z 的降幂系数num 会补齐到与 den 同长度所以严格真传递函数的 num(1) 往往是 0。这一步打印出来的 den 有几个元素后面的差分方程就必须取几拍历史这是后面排查的核心依据。2.2.2 从 Z 传递函数到差分方程设 den [1, a1, a2]、num [0, b1, b2]Z 传递函数为 (b1·z b2)/(z² a1·z a2)。由 z 的位移定理 Z[e(t−kT)] z^(−k)·E(z)交叉相乘后做逆变换得到可直接迭代的差分方程% y(k) -a1*y(k-1) - a2*y(k-2) b1*u(k-1) b2*u(k-2) yout(k) -den(2)*y_1 - den(3)*y_2 num(2)*u_1 num(3)*u_2;系数下标与历史值拍数是严格绑定的den(2) 配 y(k−1)den(3) 配 y(k−2)依此类推。原程序里 y_1、y_2 只取两拍对二阶对象是够的对被控对象是三阶的情况还缺 y(k−3) 和 u(k−3)这一项缺失会让脚本里的被控对象和 Simulink 里的完全不是同一个系统。3.1 节会把这个坑挖开讲。2.3 三套 PID 参数写法的换算表调参翻车最常见的原因不是参数不对而是三套定义混着用。连续域的并联式、串级式Ti/Td 形式、以及离散增量式里 Δu 前面的系数量纲和数量级都不一样。下表按采样周期 T 给出折算关系T1 ms 时差异是 1000 倍量级混用必然发散。写法连续域表达式增量式 Δu 的三个系数并联式C Kp Ki/s Kd·skpKpkiKi·TkdKd/T串级式C Kp(1 1/(Ti·s) Td·s)kpKpkiKp·T/TikdKp·Td/TSimulink PID 模块C Kp Ki/s Kd·s同并联式需乘/除 T 后再填脚本这张表解释了 5.3 节里那个疑问Simulink 里填 Kp6、Ki45、Kd5脚本里也写 kp6、ki45、kd5看起来一模一样实际上脚本把 Ki 放大了 1/T1000 倍、把 Kd 缩小了 1000 倍。脚本里 ki45 对应连续域积分增益 45000任何执行器都跟不上而 kd5 对应连续域微分增益 0.005阻尼几乎没加上。两条曲线当然对不上。3. MATLAB 主循环实现离散模型迭代、增量计算与输出限幅3.1 第一个坑tf(5,[1,2,1 0]) 被解析成三阶系统原程序第一行写的是 systf(5,[1,2,1 0])。MATLAB 里空格就是列分隔符[1,2,1 0] 等价于 [1 2 1 0]是四个系数对应分母 s³2s²s也就是 G(s)5/(s(s1)²)一个三阶系统跟题目给的 s²2s10 不是一回事。花十秒验证一下sysW tf(5,[1,2,1 0]); [~,denW] tfdata(sysW,v); fprintf(denW %s系数个数 %d\n, mat2str(denW), numel(denW)); % 输出denW [1 2 1 0]系数个数 4 - 三阶三阶对象离散化后 den 同样是四个元素差分方程需要 y(k−3)。脚本只写了 y_1、y_2 两拍等于把模型硬截断跑出来的曲线自然和 Simulink 对不上。可以再用 Routh 判据核对一遍三阶对象加比例增益后特征方程为 s³2s²s5Kp0临界比例增益只有 0.4而脚本里 Kp6 却能画出一条看起来还行的曲线——这种与理论明显矛盾的现象本身就是模型与差分方程不匹配的信号。所以第一件事是把分母改回 [1,2,10]并且用 numel(den) 确认阶数。3.2 主循环的更新顺序循环里最容易写乱的是顺序。合理的顺序是先算对象输出、再算误差、再算增量、再累积输出、最后滚动历史值这样物理上对应采样—计算—输出for k 1:N t(k) k*ts; rin(k) (mode1) * 1 (mode2) * 0.5*sin(2*pi*1*t(k)); % 阶跃或 1 Hz 正弦 y -den(2)*y1 - den(3)*y2 num(2)*u1 num(3)*u2; % 1) 对象输出 y(k) e rin(k) - y; % 2) 本拍误差 du kp*(e - e1) ki*e kd*(e - 2*e1 e2); % 3) PID 增量 u u1 du; % 4) 累积成绝对输出 u max(min(u, u_max), u_min); % 5) 输出限幅 yout(k)y; err(k)e; uo(k)u; u2u1; u1u; y2y1; y1y; e2e1; e1e; % 6) 滚动历史 end顺序反了会出现两类问题。先算 u 再算 y等于对象用了本拍的控制量构成代数环在离散脚本里表现为凭空快了一拍先滚动 e1、e2 再算 du增量式就退化成了位置式的近似误差曲线会多出一个稳态偏差。原程序把 x(1)、x(2)、x(3) 的更新放在循环末尾du 用的是上一拍的误差增量整体等效滞后一拍在 1 ms 采样下影响很小但概念上不如上面的写法干净。3.3 输出限幅与抗积分饱和限幅本身只是一行 max/min真正要处理的是限幅之后的积分状态。原程序的做法是先算 u再钳到 ±5然后把钳位后的值赋给 u_1这叫钳位法能挡住大部分饱和但积分项仍在内部累积长时间贴着限幅运行后误差会反向滞后。更彻底的做法是条件积分只有未饱和时才更新 u(k−1) 和积分历史。u_raw u1 du; % 未限幅的期望输出 u_sat max(min(u_raw, u_max), u_min); % 限幅后的实际输出 if u_raw u_sat u1 u_sat; % 未饱和正常累积 e2 e1; e1 e; % 未饱和才滚动误差历史 end % 饱和时冻结 u1 与误差历史误差继续反向时能立刻退出饱和对 G(s)5/(s²2s10) 这种直流增益只有 0.5 的对象单位阶跃下稳态输出是 0.5加入积分后才会爬到 1所以限幅区间 [−5,5] 相对宽松一般只在启动瞬间碰到边界。但如果把限幅收紧到 [−1,1]或者把正弦幅值调大饱和时间就会拉长这时候条件积分的差别非常明显。3.4 完整可运行脚本%% 增量式 PID 仿真G(s)5/(s^22s10)阶跃与正弦双工况 clear; clc; close all; ts 0.001; Tend 10; N round(Tend/ts); % 采样 1 ms仿真 10 s sys tf(5,[1,2,10]); % 注意分母是 [1 2 10] dsys c2d(sys,ts,zoh); [num,den] tfdata(dsys,v); fprintf(对象阶数校验den 系数个数 %d\n, numel(den)); Kp 2; Ki 10; Kd 1; % 连续域并联式增益 kp Kp; ki Ki*ts; kd Kd/ts; % 折算为差分系数关键一步 u_max 5; u_min -5; for mode 1:2 % 1阶跃 21 Hz 正弦 u10; u20; y10; y20; e10; e20; tzeros(1,N); rint; youtt; errt; uot; for k 1:N t(k) k*ts; if mode1, rin(k) 1; else, rin(k) 0.5*sin(2*pi*1*t(k)); end y -den(2)*y1 - den(3)*y2 num(2)*u1 num(3)*u2; e rin(k) - y; du kp*(e-e1) ki*e kd*(e-2*e1e2); u max(min(u1du, u_max), u_min); yout(k)y; err(k)e; uo(k)u; u2u1; u1u; y2y1; y1y; e2e1; e1e; end figure; plot(t,rin,b,t,yout,r,LineWidth,1.2); grid on; xlabel(time(s)); ylabel(rin, yout); title(sprintf(mode%d Kp%g Ki%g Kd%g, mode, Kp, Ki, Kd)); legend(rin,yout); figure; plot(t,err,r); grid on; xlabel(time(s)); ylabel(error); end脚本里 numel(den) 那句校验建议长期保留它是判断对象是不是你以为的那个对象的唯一可靠手段。kp、ki、kd 三个系数是从连续域增益折算来的改 Kp、Ki、Kd 就够不要直接改 kp、ki、kd否则下次换采样周期时全部要重算。4. 参数整定二阶弱阻尼对象上 Kp、Ki、Kd 的取法4.1 为什么从小到大调 Kp 找临界振荡在这个对象上不好用原程序里写的整定步骤是先让 kikd0把 kp 从小到大加到临界稳定再加积分、再加微分。这套流程来自经典临界比例度法前提是对象阶数较高、纯比例存在临界增益。G(s)5/(s²2s10) 是二阶对象加比例后闭环特征方程是 s²2s(105Kp)0根的实部恒为 −1与 Kp 无关也就是说纯比例在这个对象上永远临界不了只会越来越振荡。用阶跃响应指标来定 Kp 更实际for Kp [0.5 1 2 3 5] Gcl feedback(Kp*5, [1 2 10]); % 纯比例闭环 S stepinfo(Gcl); fprintf(Kp%4.1f 超调%5.1f%% 上升时间%.3fs 稳态值%.3f\n, ... Kp, S.Overshoot, S.RiseTime, dcgain(Gcl)); end跑出来会看到超调始终在 40% 以上因为对象阻尼比只有 0.316纯比例改不了阻尼。这组数据正好说明了为什么必须上微分项D 项提供的是相位超前作用是补阻尼不是锦上添花。原程序里 kp6、kd0 时曲线振荡不止本质原因就在这里。4.2 用零极点对消直接给出一组基准参数更省事的办法是让控制器零点去吃掉对象的极点。并联式 PID 写成 C(s) (Kd·s² Kp·s Ki)/s让分子等于 Kd·(s²2s10)即取 Kp2Kd、Ki10Kd。取 Kd1得到 Kp2、Ki10、Kd1闭环开环传递函数化简为 5/s闭环等效为一阶系统 5/(s5)理论无超调时间常数 0.2 s。对应的差分系数是 kp2、ki0.01、kd1000。注意 kd1000 这个数字。它之所以大是因为 Kd/T1/0.001纯粹是采样周期折算出来的不是参数调猛了。这也提示一个现实约束微分项对测量噪声的放大倍数正比于 1/T1 ms 采样下这个倍数很高仿真里没有噪声所以看不出问题一旦接实物就要给微分加一阶低通滤波。4.3 不完全微分把 kd 拉回可用范围工程上很少用理想微分常见做法是把 D 项换成带滤波的形式用一个时间常数 Tf 把高频段压下去。在增量式里最省事的实现是把二阶差分替换成一阶差分加滤波Tf 5*ts; % 滤波时间常数取 3~10 倍采样周期 alpha Tf/(Tf ts); % 一阶低通系数 ed_f alpha*ed_f (1-alpha)*(e - e1); % 滤波后的误差变化率 du kp*(e-e1) ki*e kd*ed_f; % 用滤波值替代 kd*(e-2e1e2)Tf 越大滤波越强噪声压得越狠但微分作用也越滞后通常会取 3~10 倍采样周期。改完之后 kd 可以适当放大而不抖代价是超调抑制能力下降一点需要重新微调 Kp。原程序试出来的那组 kp150、ki0.132、kd2400换算成连续域大约是 Kp150、Ki132、Kd2.4方向是对的高比例增益压静差、微分补阻尼只是因为限幅和阶数问题没能收敛。4.4 整定记录表调参一定要留记录否则改到第五轮就忘了哪组最好。建议按下面的字段记录每次只动一个参数轮次KpKiKd超调调节时间稳态误差备注1200大长0.5纯比例稳态是 0.522100减小中0积分消静差振荡仍在32101≈00.7 s0零极点对消基准43121.5小0.5 s0加快响应注意限幅超调、调节时间可以用 stepinfo 直接取也可以从 err 曲线上读最大值。判据我一般看三条阶跃下超调小于 10%、调节时间满足指标、正弦跟踪时误差幅值衰减在可接受范围。第三条对被控对象带宽有硬约束闭环等效 5/(s5) 的带宽约 5 rad/s跟踪 1 Hz6.28 rad/s正弦时幅值衰减到 62% 左右误差必然很大这是对象和参数共同决定的不是 bug。5. 仿真发散与脚本、Simulink 结果不一致的排查路径5.1 三个自动校验点遇到曲线不对先跑这三段校验比盯着参数看快得多。第一段确认对象阶数第二段确认折算关系第三段确认差分方程拍数与 den 长度匹配% 校验一对象阶数 sys tf(5,[1,2,10]); [~,den] tfdata(c2d(sys,0.001,zoh),v); assert(numel(den)3, 对象不是二阶检查分母系数); % 校验二参数折算从连续域反推回来看是否等值 Kp2; Ki10; Kd1; ts0.001; fprintf(脚本系数 kp%.4g ki%.4g kd%.4g\n, Kp, Ki*ts, Kd/ts); % 校验三差分方程拍数 fprintf(den 长度%d需要的历史拍数%d\n, numel(den), numel(den)-1);assert 那句很值钱它把分母写错这类低级问题挡在跑图之前。校验三打印出来的历史拍数如果大于脚本里保存的 y 历史个数说明模型被截断必须补 y_3、u_3或者回头把对象阶数改对。5.2 限幅引起的低频振荡现象是输出贴着 ±5 走一段然后突然反向误差长期同号大值。这不是参数问题是积分饱和。判别方法很简单把限幅区间临时放开到 [−100,100] 再跑一次如果曲线立刻变好就确认是饱和问题。处理办法在 3.3 节已经给出就是条件积分。另一个更粗暴但有效的做法是积分分离误差绝对值大于阈值时直接令 ki0小误差时才启用积分if abs(e) 0.2 du kp*(e-e1) kd*(e-2*e1e2); % 大误差段暂时关掉积分 else du kp*(e-e1) ki*e kd*(e-2*e1e2); end阈值取阶跃幅值的 10%~30% 之间比较合适太小起不到作用太大就变成纯 PD 了。5.3 采样周期与数值精度采样周期 1 ms 相对这个对象是绰绰有余的对象自然频率约 3.16 rad/s一个周期内采了 2000 个点离散化误差可以忽略。但采样周期直接决定了 kdKd/T 的大小T 从 1 ms 改成 0.1 mskd 就放大 10 倍定点实现时很容易溢出。反过来说把 T 从 1 ms 放宽到 10 mskd 缩小 10 倍更好实现代价是相位裕度下降需要重新校核。5.4 脚本与 Simulink 对齐的操作清单要让两边跑出同一条线不只是把数字抄过去。Simulink 侧要确认求解器类型固定步长、步长等于 ts、PID 模块是否设了采样时间、被控对象用的是连续传递函数还是离散模块。如果对象用的是连续 tf 模块而求解器又是变步长那和脚本里 ZOH 离散的对象天生不同误差曲线必然有差异。下表的排查顺序基本能覆盖九成情况。现象最可能原因验证动作输出始终为 0 或极小分母写成 [1,2,1 0] 且 num 前几项为 0numel(den)、mat2str(den)一开始就发散ki 或 kd 没做 T 折算打印 kp、ki、kd 三个系数波形形状对但整体偏移增益量纲不一致用 dcgain 核对稳态值输出贴限幅不退出积分饱和临时放开限幅复跑曲线有小幅毛刺理想微分放大数值误差换成不完全微分6. 把脚本封装成函数并用频域指标校核闭环脚本跑通之后下一步是做成能复用的函数不然每换一组参数就要复制一百行代码。函数签名建议把采样周期和工况一起传进去输出保持固定顺序方便批量跑参数扫描function [t, rin, yout, err, uo] incpid_sim(Kp, Ki, Kd, ts, Tend, mode) % 增量式 PID 仿真通用函数 % mode1单位阶跃 21 Hz 正弦 Kp/Ki/Kd 为连续域并联式增益 N round(Tend/ts); sys tf(5,[1,2,10]); [num,den] tfdata(c2d(sys,ts,zoh),v); kp Kp; ki Ki*ts; kd Kd/ts; % 统一在这里做折算 u10; u20; y10; y20; e10; e20; tzeros(1,N); rint; youtt; errt; uot; for k 1:N t(k) k*ts; if mode1, rin(k)1; else, rin(k)0.5*sin(2*pi*1*t(k)); end y -den(2)*y1 - den(3)*y2 num(2)*u1 num(3)*u2; e rin(k) - y; du kp*(e-e1) ki*e kd*(e-2*e1e2); u max(min(u1du,5),-5); yout(k)y; err(k)e; uo(k)u; u2u1; u1u; y2y1; y1y; e2e1; e1e; end end折算只在函数里做一次调用方永远填连续域增益这样换采样周期时上层代码不用动。批量扫描参数时可以直接套两层循环用 max(abs(err)) 作为代价函数挑最优组。真正省时间的是先把参数在频域上圈定范围再去跑时域仿真。用 margin 一次性给出幅值裕度、相位裕度和穿越频率能提前发现时域看着还行、其实裕度只有 3 dB的假稳定s tf(s); G 5/(s^2 2*s 10); C 2 10/s 1*s; % 零极点对消那组参数 [Gm, Pm, Wcg, Wcp] margin(C*G); fprintf(幅值裕度 %.1f dB相位裕度 %.1f 度穿越频率 %.2f rad/s\n, ... 20*log10(Gm), Pm, Wcp);这组参数下 LC·G 化简为 5/s输出是幅值裕度无穷大、相位裕度 90 度、穿越频率 5 rad/s和 4.2 节的理论推导吻合。如果换成 Kp6、Ki45、Kd5 直接填进脚本不做 T 折算margin 会给出相位裕度为负这就是它在时域上必然发散的原因。把频域校核放进调参流程能省掉大量改一个参数跑一次图的来回。本文还有配套的精品资源点击获取