ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

大林控制算法Matlab仿真:从原理推导到控制器实现

大林控制算法Matlab仿真:从原理推导到控制器实现 简介面向本硕博及教研人员基于MATLAB的Dalin大林控制算法仿真资源配套操作视频解决Dalin算法从理论到代码实现的编程学习问题。资源共2个文件包含1个M脚本与1个AVI操作录像整个压缩包仅159KB轻量易用。目前已有2519人浏览学习。M脚本为主程序Runme_.m可在MATLAB2021a及以上版本中直接运行通过仿真展示大林控制算法的设计流程与效果操作录像则逐步演示运行环境配置、当前目录设置及执行步骤帮助学习者快速上手。除代码和视频外资源还提示了常见运行注意事项如需在工程路径下运行、避免直接运行子函数等便于排错。适合控制类课程设计、算法仿真入门及毕业设计参考。1. 大林控制算法的仿真为什么绕不开 Matlab大林控制算法Dahlin Control Algorithm是处理纯滞后过程的经典方法在温度、流量、浓度等工业对象里滞后往往导致 PID 参数不敢调大而大林算法通过直接指定闭环脉冲传递函数来补偿滞后。真正要验证这个算法是否适合你的对象光靠手算传递函数是不够的必须在采样周期、量化效应和模型误差都被模拟的状态下观察输出曲线。Matlab 同时具备符号运算、离散化工具和 Simulink 环境是复现大林算法最顺手的地方。这篇文章按我实际做仿真的路径来写先讲清楚设计公式里每一项的来历再给出一套可以直接改参数运行的 Matlab 脚本最后聊振铃抑制和模型失配。配套的代码操作视频演示了从打开脚本到调出曲线的完整过程文字部分则补上视频里不太方便展开的参数解释。2. 大林控制算法的原理与设计参数2.1 从纯滞后对象到闭环期望传递函数大林算法解决问题的核心思路不是去反复调整 PID 参数而是把闭环系统的期望传递函数直接指定成一阶惯性加纯滞后。被控对象通常是带有纯滞后的一阶或二阶惯性环节比如温度对象常常写成[ G_0(s) \frac{K e^{-\tau s}}{T_1 s 1} ]其中 (K) 是增益(T_1) 是时间常数(\tau) 是纯滞后。数字控制器的设计目标是让整个闭环系统的传递函数 (W(z)) 近似为一个时间常数为 (T_0) 的一阶惯性加同样的滞后 (\tau)。选择一阶惯性而不是更高阶是因为一阶模型在离散化后形式简单而且对被控对象的参数变化不那么敏感。这里要注意“期望闭环”并不等同于“控制器”我们需要根据对象的离散模型反推出控制器的脉冲传递函数 (D(z))。这一步需要用到零阶保持器ZOH的离散化结果因为实际的计算机控制系统输出是阶梯波不是连续信号。所以第一步是把连续对象加上零阶保持器之后进行 z 变换得到 (G(z))然后再从期望闭环传递函数 (W(z)) 和 (G(z)) 的关系中解出 (D(z))。2.2 数字控制器 D(z) 的推导步骤最基本的闭环结构是前向通道 (D(z)) 与 (G(z)) 串联反馈为单位反馈于是闭环脉冲传递函数为[ \Phi(z) \frac{D(z)G(z)}{1 D(z)G(z)} ]假设我们要求的期望闭环传递函数 (W(z))并令 (\Phi(z)W(z))就可以解出控制器[ D(z) \frac{W(z)}{G(z) \left[ 1 - W(z) \right]} ]这个公式看起来简单但每一项都需要精确离散化。假设对象纯滞后 (\tau) 是采样周期 (T) 的整数倍记 (\tau L \cdot T)(L) 为整数且 (T_0) 为期望闭环时间常数则连续期望闭环 (W_0(s) e^{-s\tau}/(T_0 s 1)) 加上零阶保持器离散化后为[ W(z) \frac{1 - a}{1 - a z^{-1}} , z^{-(L1)}, \quad a e^{-T/T_0} ]而被控对象 (G_0(s) K e^{-s\tau}/(T_1 s 1)) 的零阶保持离散化结果为[ G(z) \frac{K(1 - b)}{1 - b z^{-1}} , z^{-(L1)}, \quad b e^{-T/T_1} ]注意这里 (z^{-(L1)}) 中的 (L1)其中 (L) 来自纯滞后额外的 (1) 来自零阶保持器本身引入的一个采样周期延迟。把这两个式子代入 (D(z)) 公式经过代数化简后得到控制器的解析表达式。这个推导过程手工做很容易在符号化简时出错我习惯先在 Matlab 里用符号工具箱验证一遍。下面是我常用的校验脚本% 符号推导大林控制器 D(z) syms z a b K real syms L integer % 期望闭环离散传递函数 W (1 - a) / (1 - a * z^(-1)) * z^(-(L1)); % 被控对象离散模型一阶惯性加纯滞后 G K * (1 - b) / (1 - b * z^(-1)) * z^(-(L1)); % D(z) W / (G * (1 - W)) D W / (G * (1 - W)); D simplify(D); % 打印化简后的结果 pretty(D)代码里用syms声明了符号变量simplify会把分子分母公因子约掉。运行后你会看到 D(z) 的分子和分母都是关于 (z^{-1}) 的多项式并没有出现 (z^{-(L1)}) 项因为 (W/G) 中的滞后因子互相抵消了。也就是说大林控制器的阶次只和对象阶次相关和滞后步数无关这是它区别于史密斯预估器的一个重要特点。2.3 需要手工指定的三个关键参数使用大林算法时有三个参数直接决定控制效果。下表总结了它们的含义、推荐范围和调整倾向。参数含义推荐范围调整说明采样周期 (T)计算机控制器的采样间隔满足香农定理且建议 (T \tau / 4)(T) 过大会导致离散模型失真过小会放大控制器输出抖动期望闭环时间常数 (T_0)闭环系统希望达到的惯性时间一般为对象时间常数 (T_1) 的 0.30.5 倍(T_0) 越小响应越快但对模型误差和噪声越敏感滞后步数 (L)纯滞后除以采样周期后取整必须满足 (L \ge 1)尽量使 (\tau/T) 接近整数若 (\tau/T) 非整数离散模型会产生额外误差需要近似处理这里特别强调采样周期 (T) 的选择。大林算法的基础是滞后时间恰好等于采样周期的整数倍当 (\tau / T) 不是整数时离散模型会产生额外的高阶项直接套公式会导致控制器复杂化。遇到这种情况我一般先尝试把采样周期微调一下让 (L \text{round}(\tau / T))然后在仿真中检查闭环是否稳定。如果不行就需要在期望闭环中添加一个额外的零点来补偿小数滞后但这已经属于进阶话题。提示三个参数不是独立调整的。(T) 变小会让 (L) 变大而 (T_0) 与 (T) 的关系又影响 (a) 的取值所以改参数后一定要重新计算 (a)、(b) 和控制器系数不要只在仿真里改某个数字。3. 用 Matlab 实现大林控制算法的代码骨架3.1 建立被控对象并离散化先用传递函数对象建立被控对象然后用c2d函数加零阶保持器离散化。下面这段代码是完整脚本的第一部分% 被控对象参数 K 2.0; % 增益 T1 5.0; % 对象时间常数 tau 1.0; % 纯滞后时间 Ts 0.5; % 采样周期 % 连续对象模型G(s) K*exp(-tau*s)/(T1*s1) s tf(s); G0 K * exp(-tau*s) / (T1*s 1); % 零阶保持器离散化 Gz c2d(G0, Ts, zoh); disp(Gz)代码说明c2d是 Control System Toolbox 中的离散化函数第三个参数zoh表示使用零阶保持器这是数字控制系统最常用的离散化方式因为它模拟了 DAC 和保持电路的实际行为。运行后会在命令行打印 Gz 的表达式例如0.09516 z^-6 / (z - 0.9048)其中z^-6对应滞后步数 (L2) 以及零阶保持器带来的额外一拍。注意这里exp(-tau*s)在连续模型里可以被c2d识别为输入延迟为了更严谨也可以写成set(G0,InputDelay,tau)效果相同。如果不想让c2d的处理变成黑箱也可以自己写出离散模型。假设滞后步数 (L \tau / T_s)一阶惯性离散化后是[ G(z) \frac{K(1 - b)}{1 - b z^{-1}} , z^{-(L1)} ]其中 (b e^{-T_s/T_1})。你可以用这个公式和c2d的输出进行比较确认两者一致后再继续向下设计控制器这样后面每一步模型都是可追踪的。3.2 将符号结果转换为数值传递函数上一章的符号推导给出了 D(z) 的通用形式但在仿真中我们需要一个具体的数值离散传递函数。下面这段代码直接按公式计算并利用 Control System Toolbox 的tf对象搭建控制器。% 设定参数 T0 2.0; % 期望闭环时间常数取 0.4*T1 a exp(-Ts / T0); b exp(-Ts / T1); L round(tau / Ts); % 创建离散时间变量 z z tf(z, Ts); % 对象 G(z) Gz K * (1 - b) / (1 - b * z^(-1)) * z^(-(L1)); % 期望闭环 W(z) Wz (1 - a) / (1 - a * z^(-1)) * z^(-(L1)); % 控制器 D(z) Dz (Wz / Gz) / (1 - Wz); Dz minreal(Dz); % 查看控制器表达式 disp(Dz)minreal用来消除分子分母中的公因子。由于第 2 章符号推导提到 W/G 会消去滞后项但这里用tf对象计算时数值舍入可能导致零极点没有严格对消minreal可以强制约分。如果不调用后续仿真可能出现状态空间不满秩的问题。3.3 闭环仿真的脚本式实现得到控制器 Dz 后可以用feedback直接构造闭环并绘制阶跃响应。完整代码如下% 前向通路Dz * Gz Go Dz * Gz; % 单位负反馈闭环 closed_loop feedback(Go, 1); % 仿真时间序列 t 0:Ts:20; % 阶跃响应 [y, t_out] step(closed_loop, t); % 绘图 figure; plot(t_out, y, LineWidth, 1.5); grid on; xlabel(时间 (s)); ylabel(输出); title(大林控制算法闭环阶跃响应);代码里的feedback(Go, 1)是对前向通路做单位负反馈等价于 (G_o/(1G_o))。step在使用离散传递函数时自动按采样时刻输出时间向量t的步长要和Ts保持一致。运行这段代码得到一条平滑的、近似一阶惯性加滞后的曲线超调量很小响应时间接近 (T_0 \tau)。如果想在 Simulink 中复现同样的闭环可以在模型中放置两个Discrete Transfer Fcn模块分别填入 Dz 和 Gz 的分子分母系数。系数可以从工作区变量提取Dz.Numerator{1}、Dz.Denominator{1}。然后使用Sum和Gain搭出负反馈结构。这个搭建过程在配套的代码操作视频里有完整演示视频中还演示了如何用simplot将 Simulink 的仿真结果与脚本代码的阶跃响应画在同一张图上。4. 大林控制算法仿真结果分析与边界4.1 振铃现象及其消除方法大林算法在实际应用中最大的问题是振铃。振铃是指控制器的输出在采样时刻之间快速正负波动而系统的实际输出可能看起来正常或者只是轻微振荡。振铃通常由控制器 D(z) 中位于 z 平面负实轴附近的极点引起。当期望闭环时间常数 (T_0) 取取得太小或者采样周期 (T_s) 相对较大时这些极点可能接近 -1导致控制器输出剧烈摆动。判断是否存在振铃可以在 Matlab 里绘制控制器的单位脉冲响应或者直接计算参考输入到控制器的传递函数再观察阶跃输入下的控制器输出。下面的代码提取并绘制控制信号% 参考输入到控制器输出的传递函数 U/R R_to_U Dz / (1 Dz * Gz); R_to_U minreal(R_to_U); % 阶跃输入下的控制器输出 [u, t_u] step(R_to_U, t); figure; plot(t_u, u, LineWidth, 1.5); title(控制器输出观察振铃);如果控制器输出波形在相邻采样点之间以接近采样频率的频率上下跳动说明存在振铃。消除振铃的标准方法是修正控制器把控制器中引起振铃的因子替换为它的稳态值。具体来说假设 D(z) 中可以分解出一个接近 -1 的极点 (z_p)则将因子 ((z - z_p)) 中的 (z) 替换为 1得到 ((1 - z_p))。这个小操作会让控制器在低频段的增益基本不变但会去除高频振荡分量。典型一阶对象修正后的控制器形式为[ D(z) D(z) \cdot \frac{(1 - z_p)}{z - z_p} ]上式在理论上抵消了振铃极点但工程上更稳妥的做法是让 (D(z)) 与原始 (D(z)) 在 (z1) 处的增益保持一致也就是对剩余部分乘以一个比例系数修正。我在实际仿真中常把振铃因子消除后的控制器再代入闭环对比修正前后的控制器输出确认振铃幅度降低了 80% 以上。4.2 参数失配与仿真发散排查仿真发散通常不是算法本身的问题而是模型离散化出错或参数不合理。常见的原因有三个采样周期过大导致 (b) 或 (a) 过小离散模型失真滞后步数 (L) 计算错误比如 (\tau / T_s) 不是整数时直接取整导致 (G(z)) 的滞后步数与 (W(z)) 不一致控制器 D(z) 中的零极点未约分或者反馈路径符号错误。遇到发散时我的排查顺序是先打印 Gz 和 Dz 的零点极点看看是否都在单位圆内然后在 Simulink 中逐步检查每个模块的采样时间是否一致最后缩小仿真步长和采样周期。下面这个代码可以快速检查闭环极点% 检查闭环极点 pzmap(closed_loop); % 计算最大极点模值 max_abs_pole max(abs(pole(closed_loop))); disp([最大极点模值 , num2str(max_abs_pole)]);如果最大极点模值略大于 1需要减小 (T_0) 或增大 (T_s)然后重新设计控制器。另外提醒一下round(tau/Ts)的取整方式会引入模型误差如果取整后误差超过 5%建议改用更精确的离散化方法或者将期望闭环的滞后步数设置成与对象完全相同的值。4.3 和 PID 控制的对比在相同对象上做一组对比仿真能让你快速理解大林算法的适用场景。我通常使用pidtune自动整定一个 PID 控制器然后与修正后的大林算法放在同一张图里比较。以下是对比代码的骨架% 设计用于对比的 PID Gc_pid pidtune(G0, PID); loop_pid feedback(Gc_pid * G0, 1); % 大林闭环使用修正后的 Dz Dz_corrected ...; % 将修正后的控制器赋值给这个变量 loop_dalin feedback(Dz_corrected * Gz, 1); % 统一时间轴并绘图 t_cont 0:0.01:20; step(loop_pid, t_cont); hold on; step(loop_dalin, t_cont); legend(PID, 大林);对比结果通常显示大林算法对滞后系统的响应更平滑、超调更小但对模型增益和时间常数的变化更敏感。PID 调得好的话可以接近大林的效果但需要更多适配。因此大林算法适合对象参数比较稳定、滞后明确的场合而 PID 适合参数摄动大的场合。这个结论在我的仿真中多次重复出现。5. 将大林控制算法封装成可复用函数并验证5.1 封装设计函数把大林控制器设计封装成函数可以避免每次仿真都复制粘贴一段脚本。下面这个函数接收对象参数和期望参数返回控制器、对象离散模型和闭环模型。function [Dz, Gz, closed_loop] design_dahlin(K, T1, tau, Ts, T0) % DESIGN_DAHLIN 设计大林控制器 % 输入 % K 对象增益 % T1 对象时间常数 % tau 纯滞后时间 % Ts 采样周期 % T0 期望闭环时间常数 % 输出 % Dz 大林控制器离散传递函数 % Gz 被控对象零阶保持离散模型 % closed_loop 闭环传递函数 a exp(-Ts / T0); b exp(-Ts / T1); L round(tau / Ts); z tf(z, Ts); Gz K * (1 - b) / (1 - b * z^(-1)) * z^(-(L1)); Wz (1 - a) / (1 - a * z^(-1)) * z^(-(L1)); Dz minreal((Wz / Gz) / (1 - Wz)); closed_loop feedback(Dz * Gz, 1); end使用时只需调用[D, G, cl] design_dahlin(2, 5, 1, 0.5, 2); step(cl);这个函数把第 2 章的公式全部翻译成了离散传递函数注意z^(-(L1))中的 L1 是因为零阶保持器带来的一个采样周期延迟和公式推导保持一致。5.2 验证零极点位置封装后一定要做一次验证。我习惯用一组随机参数跑批量仿真检查闭环极点是否都在单位圆内以及阶跃响应的稳态值是否为 1。下面的验证代码基于assert一旦发现不稳定会立刻报错param_list [ 2 5 1 0.5 2 1 3 0.8 0.2 1 0.5 10 2 1 4]; for i 1:size(param_list, 1) [Dz, Gz, cl] design_dahlin(... param_list(i,1), param_list(i,2), param_list(i,3), ... param_list(i,4), param_list(i,5)); assert(max(abs(pole(cl))) 1, 闭环不稳定参数组 %d, i); end这段代码中pole(cl)返回闭环传递函数的所有极点只要有一个极点模值大于等于 1assert就会中止脚本。批量验证的意义在于大林算法的参数空间里存在稳定域边界单纯一组参数看不出来多组才能暴露问题。5.3 配合操作视频的复现要点配套的代码操作视频基于上面这份函数封装演示了三个操作第一如何在 Matlab 命令行中调用design_dahlin函数并查看输出变量第二如何把 Gz 和 Dz 导出到 Simulink 模型进行协同仿真第三如何利用legend或compare观察不同 T0 下的响应差异。我这里要给出的复现建议是视频中使用的 Matlab 版本建议在 R2019a 以上因为老版本的tf对离散时间变量z的处理方式稍有不同。如果你用的是 R2018a 或更老版本需要先用c2d显式转换。关于代码操作视频本身我的观点是它更适合用来确认“应该看到什么结果”而不是替代自己敲代码。因为大林算法仿真最常踩的坑不是公式推导而是离散化方式不一致。视频里展示了从 G0 到 Gz 的转换过程如果你想复现最好手打一遍代码而不是直接粘贴这样能体会到InputDelay和exp(-tau*s)在c2d中的差异。对于想要快速上手的人我这里提供一个最小复现清单Matlab 环境变量和当前路径要保持干净避免变量名冲突先运行clear all; close all;清空工作区用format long检查控制器系数的小数位避免中间舍入误差如果发现仿真发散优先检查L round(tau/Ts)这一行是否与对象的InputDelay属性一致。以上细节在文字版里都已展开视频里只是快速操作两者搭配使用效果最好。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进