简介面向自动控制、电子信息工程等专业学生及工程技术人员这份资源以Matlab/Simulink为平台实现了双罐系统的广义预测控制GPC仿真可用于课程设计、期末大作业及毕业设计中的控制算法验证与学习。压缩包共8个文件包含GPC系数计算、控制器及主程序等4个m脚本一个Simulink模型mdl文件以及说明文档txt和效果图jpg整体仅53KB轻量易用。代码采用参数化编程注释详细便于修改参数并观察系统响应新手也能通过readme快速上手并替换数据运行。已有58人学习下载适合希望理解GPC预测控制原理并快速搭建仿真实验的读者。1. 双罐系统用上 GPC从液位耦合到能“看几步”的控制器双罐系统是过程控制里最常见的二阶耦合对象上罐进水量同时影响上、下两罐液位下罐动态又反作用于流量差两个液位相互牵扯单纯 PID 很难在不牺牲响应速度的前提下同时压住耦合。GPC广义预测控制这一类基于模型预测的控制算法恰好用“未来一段时间的输出预测 滚动优化”来应对这种耦合和滞后这几年在 Simulink 仿真、课程设计、汽车热管理液位控制里出现频率都很高。把这套控制放在 Simulink 里实现核心不是拖一个 MPC 模块而是要把被控对象线性化、预测矩阵和在线求解三步理顺。下面按一条可复现的路径讲从双罐模型建立到 GPC 算法最小实现再到仿真回路和调参检查。2. 双罐对象的机理建模与可预测的线性化模型2.1 双罐系统的状态方程与参数表常见双罐结构是两个圆罐串联入水流量qin进 1 号罐1 号罐底部阀门连接 2 号罐2 号罐底部阀门出水。流量采用常用的平方根模型dh1/dt (qin - c1*sqrt(h1) - c12*sqrt(max(h1-h2,0))) / A1 dh2/dt (c12*sqrt(max(h1-h2,0)) - c2*sqrt(h2)) / A2sqrt(max(...))是为了防止仿真中出现负开方导致仿真中断这个写法在双罐 Simulink 模型中比单纯sqrt(h1-h2)更稳妥。参数含义本例取值A11 号罐截面积2.0 m²A22 号罐截面积1.5 m²c11 号罐出水阀系数1.2c22 号罐出水阀系数1.0c12连接管阀系数0.9h1_01 号罐操作点液位2.0 mh2_02 号罐操作点液位1.0 m在 Simulink 里搭双罐子系统时我一般不用 MATLAB Function 直接写微分方程而是用两个 Integrator 块构成状态h1、h2作为积分输出各自积分输入端用 Fcn 块按上面两式计算导数。这样后续做线性化、看状态变量、叠加扰动都更直观。给 h1、h2 的 Integrator 初始值分别填2.0和1.0输入口是qin输出口是两个液位。2.2 在工作点线性化得到 GPC 需要的离散状态空间GPC 的预测模型一般是线性模型直接在非线性模型上做预测不是不行但代码和计算量会明显上升。工程上更常见的做法是在工作点这里取 h12.0、h21.0附近线性化得到一个两输入两输出的线性状态空间模型再用 GPC 控制液位偏差。由于 GPC 只预测“相对操作点的增量”后续在回路里需要把绝对液位换算成偏差量。下面这段脚本可以直接放到模型目录下运行它给出了从机理参数到离散状态空间的完整转换。% two_tank_lin.m % 双罐系统在工作点附近的连续线性化 A12.0; A21.5; c11.2; c21.0; c120.9; h1_02.0; h2_01.0; dh_0h1_0-h2_0; % 平方根流量在工作点的导数 q1_der c1 / (2*sqrt(h1_0)); q12_der c12 / (2*sqrt(dh_0)); q2_der c2 / (2*sqrt(h2_0)); % 连续状态矩阵 dx Ac*x Bc*u % 两个状态为 dH1, dH2输入为 qin 的增量 Ac [-(q1_derq12_der)/A1, q12_der/A1; q12_der/A2, -(q2_derq12_der)/A2]; Bc [1/A1; 0]; Cc eye(2); Dc zeros(2,1); % 采样周期GPC 控制周期不是仿真步长 Ts 2.0; sys_d c2d(ss(Ac,Bc,Cc,Dc), Ts, zoh); Ad sys_d.A; Bd sys_d.B; Cd sys_d.C;这段脚本的关键是矩阵里的q12_der连接管流量对h1-h2求导后会同时出现在两个状态方程的耦合项里。很多初学者的错误是只在h1方程里加耦合项忽略h2方程里的同源项结果阶跃响应方向上就会出现偏差。2.3 用阶跃响应验证离散模型线性化结果合不合理不要等 GPC 跑完再看。最直接的做法是对连续非线性模型和离散线性模型施加同一个qin阶跃比较下罐液位h2的响应曲线。% 对离散线性模型做阶跃响应 sim_t 500; [y_lin, t_lin] step(sys_d, sim_t); % 对 Simulink 非线性模型做同样阶跃得到 y_nl % 两条曲线在操作点附近应基本重合如果两条曲线在最初几个采样周期内不重合优先检查c2d的采样时间Ts是否选得太大。双罐系统时间常数通常在几十秒到上百秒Ts2s是一个合理起步值若Ts超过系统最小时间常数的 1/5离散模型会丢失动态GPC 预测就会出现明显偏差。3. GPC 控制器核心算法从预测矩阵到滚动优化3.1 GPC 到底在“预测”什么GPC 与常规状态反馈控制的本质区别在于控制器不只依赖当前状态而是基于模型把未来 Np 步的输出预测出来然后找一串 Nc 步的控制增量使预测输出尽量靠近参考轨迹同时限制控制增量幅值。双罐系统的难点在于一个控制量qin要同时影响两个输出如果没有预测当前控制只能靠的“反馈误差”去校正有预测后控制器能提前知道这一步qin对h1和h2未来几步的影响从而在耦合通道之间做权衡。通常把预测时域 Np 取为 1025控制时域 Nc 取为 15。Np 太小控制器看不到闭环动态尾巴等效于一个近近视控制无法体现预测优势Np 太大计算量上升而且模型误差会被远期预测放大。3.2 用增量型状态空间表达 GPC 目标函数为了处理控制增量我把双罐线性模型扩写成带有上一拍输入的增广状态z(k) [dH1(k); dH2(k); du(k-1)]增广后的离散状态方程为z(k1) Aa*z(k) Ba*dU(k) y(k) Ca*z(k)其中Aa [Ad, Bd; zeros(1,2), eye(1)]; Ba [Bd; eye(1)]; Ca [Cd, zeros(2,1)];这样做的意义是让优化变量直接变成dU(k)而不再需要单独记录上一拍控制量。目标函数写成J sum_{j1}^{Np} || y(kj|k) - r(kj) ||^2_Q sum_{j0}^{Nc-1} || dU(kj) ||^2_R这个形式和无约束情况下 Clarke 的 CARIMA-GPC 等价但在 Simulink 里用状态空间矩阵实现更直接也不容易出现丢番图方程递推的编程错误。3.3 Gamma 预测矩阵和最小二乘解核心计算是构造预测矩阵Gamma。Gamma(i,j)表示第j个控制增量对未来第i步输出的影响。利用增广矩阵可以统一写成% build_gpc_matrices.m % 在模型初始化脚本中运行生成 P 结构体 Np 15; Nc 3; ny 2; % 两个液位输出 nu 1; % qin 一个输入 Aa [Ad, Bd; zeros(nu, size(Ad,2)), eye(nu)]; Ba [Bd; eye(nu)]; Ca [Cd, zeros(ny, nu)]; Gamma zeros(Np*ny, Nc*nu); for i 1:Np for j 1:Nc if i j Gamma((i-1)*ny1:i*ny, (j-1)*nu1:j*nu) Ca * (Aa^(i-j)) * Ba; end end end % 输出加权1 号罐液位权重低2 号罐液位权重高 Q_tmp diag([0.2, 1.0]); Qbar kron(eye(Np), Q_tmp); Rbar 0.5 * eye(Nc*nu); % 操作点 h0 [2.0; 1.0]; u0 c1*sqrt(h0(1)) c12*sqrt(h0(1)-h0(2));kron的作用是把单步输出权重展开到整个预测时域这样写比一层一层循环赋值更清楚。u0是操作点流量由两个罐在稳态时的流量平衡计算得到绝对流量必须从这个点往上加。每个控制周期在线只需要求解一个无约束线性最小二乘问题% 当前状态偏差 xk [h1_in; h2_in] - P.h0; zk [xk; u_km1 - P.u0]; % 自由响应假设未来控制增量全为 0 Y0 zeros(P.Np * P.ny, 1); Ak eye(size(P.Aa)); for k 1:P.Np Ak Ak * P.Aa; Y0((k-1)*P.ny1:k*P.ny) P.Ca * Ak * zk; end % 参考轨迹按两个输出分别给定 Rk zeros(P.Np * P.ny, 1); Rk(1:2:end) ref1 - P.h0(1); Rk(2:2:end) ref2 - P.h0(2); % 无约束解 H P.Gamma * P.Qbar * P.Gamma P.Rbar; dU H \ (P.Gamma * P.Qbar * (Rk - Y0)); u_next u_km1 dU(1);代码里Ak累乘得到Aa^k比每次都调用Aa^extra更省计算也更好调试。dU(1)只取第一个控制增量这一步就是滚动优化下一拍重新计算永远不会把整条控制序列直接发给被控对象。3.4 控制器参数表参数含义起步值取值范围Ts控制周期2 s系统时间常数 1/10 上下Np预测时域151025Nc控制时域315Q(1,1)h1 液位权重0.201Q(2,2)h2 液位权重1.0110R控制增量权重0.50.15R 越小控制增量越激进执行器阀门动作越大R 过大时系统会变得“不敢动”双罐液位跟踪速度明显下降。4. 在 Simulink 中集成双罐 GPC回路、函数块和采样设置4.1 为什么不用现成 MPC ToolboxSimulink 自带的 Model Predictive Control Toolbox 能直接识别状态空间对象但它在 GPC 细节上不一定完全匹配课程设计或论文里要求的 CARIMA 模型同时授权成本也高。实际做双罐 GPC 仿真时更多人会选择自己写 MATLAB Function 块因为双罐模型只有两个状态、一个输入实时计算量很小用最小二乘解完全够用。下面沿着这个路线搭回路。4.2 在 MATLAB Function 块中写在线 GPC 函数在 Simulink 库中放一个 MATLAB Function 块双击后把函数名改成gpc_step代码直接引用上一节的计算逻辑。注意函数入参里要包含结构体P这个P是在模型初始化脚本里生成的参数对象。function u_next gpc_step(h1, h2, u_km1, ref1, ref2, P) %#codegen % 双罐系统 GPC 控制器基于增量型状态空间模型 xk [h1; h2] - P.h0; zk [xk; u_km1 - P.u0]; Y0 zeros(P.Np * P.ny, 1); Ak eye(size(P.Aa)); for k 1:P.Np Ak Ak * P.Aa; Y0((k-1)*P.ny1:k*P.ny) P.Ca * Ak * zk; end Rk zeros(P.Np * P.ny, 1); Rk(1:2:end) ref1 - P.h0(1); Rk(2:2:end) ref2 - P.h0(2); H P.Gamma * P.Qbar * P.Gamma P.Rbar; dU H \ (P.Gamma * P.Qbar * (Rk - Y0)); u_next u_km1 dU(1); end函数内部没有persistent变量所有历史状态都通过外部输入u_km1带回。这样最便于调试也方便后续把u_km1换成阀门执行器反馈值。P.ny必须与实际输出个数一致这里两个液位都参与加权所以ny2。4.3 外围回路单位延迟、参考值和操作点叠加搭建控制器外围接线时按下面顺序操作从双罐非线性模型取出h1、h2两个输出分别接到gpc_step函数块的h1、h2输入。在控制器输出端接一个 Unit Delay 块输出作为u_km1反馈回控制器输入。Unit Delay 初始值填u0。将常数值u0与gpc_step的输出相加再加单位延迟前的信号不需要直接把qin u0 u_next送到双罐子系统。参考值ref2用 Step 或 Constant 生成ref1可以设成与h0(1)相同表示 1 号罐液位只需维持在当前操作点不要主动跟踪。需要特别注意的操作点是gpc_step输出的u_next是绝对流量还是控制增量我这里返回的是绝对流量因为u_km1是绝对流量。如果你把u_km1改成上一拍控制增量那么所有偏差换算都要重构最容易犯错的是u0被加了两遍。4.4 必须对齐的采样时间与求解器设置GPC 是离散控制器它的采样时间Ts只在控制器和 Unit Delay 上体现。Simulink 中被控对象是连续积分求解器应选择固定步长例如ode4步长取 0.1 s。不要用变步长求解器搭配离散控制器那会让 GPC 的控制间隔不是等间隔在线求解的预测模型时间基准就不成立。在 MATLAB Function 块的 Block Parameters 中Sample time 填2Unit Delay 的 Sample time 也填2。其余连续积分模块继承固定步长。检查方法是在仿真诊断窗口里看是否有黄色的“采样时间继承冲突”提示如果有优先检查 Unit Delay 是否误用了 Inherited。5. 参数调整与扰动抑制实验把双罐 GPC 调到可用状态5.1 先调 Np再调 Nc最后调 R很多人在 Simulink 里一上来就同时改四个参数结果曲线发散或抖动时根本不知道是哪一项造成的。常见做法是先固定Nc2、R0.5只改Np。对双罐系统你会发现 Np 太小比如 5时 GPC 和带前馈的 PID 差不多下罐超调明显Np 增大到 15 后控制器能看到完整的二阶响应超调显著下降。继续增大到 25改善很小但矩阵计算量线性上升。5.2 三组参数对比实验在同一个仿真模型里把ref2在 50 秒时从 1.0 m 阶跃到 1.2 m观察 h2 的响应和 qin 的变化得到下表。参数组合h2 超调调节时间控制量波动Np10, Nc2, R0.3约 12%约 160 s明显Np15, Nc3, R0.5约 4%约 210 s中等Np20, Nc5, R2.0几乎为 0约 320 s很小这三组数据不是理论公式而是仿真曲线上的典型读数。它们说明一个规律双罐 GPC 的“敏感旋钮”首先是 Np它决定控制器是否把耦合动态看全其次是 R它决定控制动作是否温和Nc 的作用更多是把控制轨迹的自由度打开Nc 从 2 变到 5 时响应形状会变化但不会出现 PID 那种临界振荡效应。5.3 液位耦合下的设定值跟踪与扰动实验为了验证耦合抑制能力建议在仿真 300 秒时给qin叠加入口扰动比如模拟进水管压力下降将 qin 直接减去 0.1 m³/s 持续 60 秒。观察两个液位的变化曲线。如果 GPC 调得合适h2 的最大偏差应小于设定值阶跃偏差的 1/3h1 因为加权较低会出现更大偏差但不会发散。实验中偶尔会遇到 h2 响应出现“先反方向走再回头”的现象。这不是控制算法错误而是耦合通道的动态下罐先被上罐液位抬高随后 GPC 降低进水量下罐又回落。若反方向幅度太大就把 Q 里 h1 的权重调大比如从 0.2 改到 0.5让控制器少用“抬高上罐”的方式去临时维持下罐流量。5.4 从仿真曲线反推参数原则参数调整的依据始终是曲线形状而不是误差积分指标。双罐系统里我一般看三个特征第一个是 h2 阶跃后是否有“回头”第二个是 qin 是否频繁触发饱和或速率限制第三个是稳态时 dU 是否在零附近抖动。如果 qin 持续高频抖动先增加 R而不是减小 Np。如果 qin 长时间顶在边界上说明 Nc 太小控制轨迹自由度不足。6. 把 GPC 仿真结果落到实时系统的检查方法Simulink 里跑通 GPC 只是第一步真正要把它用于外部硬件或自动化代码生成还需要做两件事第一是校验预测矩阵与离线仿真的一致性第二是确认采样时间和初始化逻辑在生成代码后仍然成立。校验预测矩阵用回归测试最省事。单独写一个脚本把dU设成一组固定随机值用gpc_step里相同逻辑算预测输出再用双罐离散模型实际迭代一遍两者误差小于1e-8才算通过。% verify_prediction.m dU_test [0.1; -0.02; 0.05]; u_km1 u0; zk [0.1; -0.05; 0]; % 初始偏差 Y_pred Y0_free Gamma * dU_test; x zk(1:2); u_cur u_km1; for i 1:Np if i Nc u_cur u_cur dU_test(i); end x Ad * x Bd * (u_cur - u0); Y_sim((i-1)*ny1:i*ny) Cd * x; end if norm(Y_pred - Y_sim, Inf) 1e-8 disp(预测矩阵校验通过); end这里最容易犯的错误是把u_cur - u0写成u_cur。因为Ad、Bd是从偏差量推导出来的输入必须是相对操作点的流量增量否则稳态会直接偏离。校验收敛后再看实时实现。如果目标平台支持 Simulink Coder常见做法是先把 GPC 函数块调整为“外部模式”在硬件上跑一圈验证通讯和采样节拍。重点检查三处MATLAB Function 块是否被识别为离散任务Unit Delay 初始值是否通过u0正确初始化H矩阵是否在初始化时预先求逆而不是每个周期重新分解。对于双罐系统这种小规模矩阵在线求逆没有问题但为了形成稳定的执行时间我通常把H的逆或者 Cholesky 分解结果缓存在函数块内部每个周期只做矩阵乘法。另一个实时化技巧是在生成代码前关闭 Simulink 的无限大矩阵检查并给gpc_step标注%#codegen确保所有变量尺寸固定。如果后续要封装成 FMU 或 DLL 供其他仿真工具调用同样可以沿用这套函数接口只把 Simulink 外围的 Unit Delay 和参数结构体一起封装进 S-Function 或导出的 C 接口里。本文还有配套的精品资源点击获取