企业官网建设流程全解析

简介：面向自动控制、电子信息工程等专业学生及工程技术人员，这份资源以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)) / A2

sqrt(max(...))是为了防止仿真中出现负开方导致仿真中断，这个写法在双罐 Simulink 模型中比单纯sqrt(h1-h2)更稳妥。

在 Simulink 里搭双罐子系统时，我一般不用 MATLAB Function 直接写微分方程，而是用两个 Integrator 块构成状态：h1、h2作为积分输出，各自积分输入端用 Fcn 块按上面两式计算导数。这样后续做线性化、看状态变量、叠加扰动都更直观。给 h1、h2 的 Integrator 初始值分别填2.0和1.0，输入口是qin，输出口是两个液位。

2.2 在工作点线性化，得到 GPC 需要的离散状态空间

GPC 的预测模型一般是线性模型，直接在非线性模型上做预测不是不行，但代码和计算量会明显上升。工程上更常见的做法是在工作点（这里取 h1=2.0、h2=1.0）附近线性化，得到一个两输入两输出的线性状态空间模型，再用 GPC 控制液位偏差。由于 GPC 只预测“相对操作点的增量”，后续在回路里需要把绝对液位换算成偏差量。

下面这段脚本可以直接放到模型目录下运行，它给出了从机理参数到离散状态空间的完整转换。

% two_tank_lin.m % 双罐系统在工作点附近的连续线性化 A1=2.0; A2=1.5; c1=1.2; c2=1.0; c12=0.9; h1_0=2.0; h2_0=1.0; dh_0=h1_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_der+q12_der)/A1, q12_der/A1; q12_der/A2, -(q2_der+q12_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是否选得太大。双罐系统时间常数通常在几十秒到上百秒，Ts=2s是一个合理起步值；若Ts超过系统最小时间常数的 1/5，离散模型会丢失动态，GPC 预测就会出现明显偏差。

3. GPC 控制器核心算法：从预测矩阵到滚动优化

3.1 GPC 到底在“预测”什么

GPC 与常规状态反馈控制的本质区别在于：控制器不只依赖当前状态，而是基于模型把未来 Np 步的输出预测出来，然后找一串 Nc 步的控制增量，使预测输出尽量靠近参考轨迹，同时限制控制增量幅值。双罐系统的难点在于一个控制量qin要同时影响两个输出，如果没有预测，当前控制只能靠的“反馈误差”去校正；有预测后，控制器能提前知道这一步qin对h1和h2未来几步的影响，从而在耦合通道之间做权衡。

通常把预测时域 Np 取为 10～25，控制时域 Nc 取为 1～5。Np 太小，控制器看不到闭环动态尾巴，等效于一个近近视控制，无法体现预测优势；Np 太大，计算量上升，而且模型误差会被远期预测放大。

3.2 用增量型状态空间表达 GPC 目标函数

为了处理控制增量，我把双罐线性模型扩写成带有上一拍输入的增广状态：

z(k) = [dH1(k); dH2(k); du(k-1)]

增广后的离散状态方程为：

z(k+1) = 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_{j=1}^{Np} || y(k+j|k) - r(k+j) ||^2_Q + sum_{j=0}^{Nc-1} || dU(k+j) ||^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)*ny+1:i*ny, (j-1)*nu+1: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.ny+1: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 控制器参数表

R 越小，控制增量越激进，执行器阀门动作越大；R 过大时系统会变得“不敢动”，双罐液位跟踪速度明显下降。

4. 在 Simulink 中集成双罐 GPC：回路、函数块和采样设置

4.1 为什么不用现成 MPC Toolbox

Simulink 自带的 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.ny+1: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必须与实际输出个数一致，这里两个液位都参与加权，所以ny=2。

4.3 外围回路：单位延迟、参考值和操作点叠加

搭建控制器外围接线时，按下面顺序操作：

1. 从双罐非线性模型取出h1、h2两个输出，分别接到gpc_step函数块的h1、h2输入。
2. 在控制器输出端接一个 Unit Delay 块，输出作为u_km1反馈回控制器输入。Unit Delay 初始值填u0。
3. 将常数值u0与gpc_step的输出相加，再加单位延迟前的信号？不需要，直接把qin = u0 + u_next送到双罐子系统。
4. 参考值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 填2，Unit Delay 的 Sample time 也填2。其余连续积分模块继承固定步长。检查方法是在仿真诊断窗口里看是否有黄色的“采样时间继承冲突”提示，如果有，优先检查 Unit Delay 是否误用了 Inherited。

5. 参数调整与扰动抑制实验：把双罐 GPC 调到可用状态

5.1 先调 Np，再调 Nc，最后调 R

很多人在 Simulink 里一上来就同时改四个参数，结果曲线发散或抖动时根本不知道是哪一项造成的。常见做法是先固定Nc=2、R=0.5，只改Np。对双罐系统，你会发现 Np 太小（比如 5）时 GPC 和带前馈的 PID 差不多，下罐超调明显；Np 增大到 15 后，控制器能看到完整的二阶响应，超调显著下降。继续增大到 25，改善很小，但矩阵计算量线性上升。

5.2 三组参数对比实验

在同一个仿真模型里，把ref2在 50 秒时从 1.0 m 阶跃到 1.2 m，观察 h2 的响应和 qin 的变化，得到下表。

这三组数据不是理论公式，而是仿真曲线上的典型读数。它们说明一个规律：双罐 GPC 的“敏感旋钮”首先是 Np，它决定控制器是否把耦合动态看全；其次是 R，它决定控制动作是否温和；Nc 的作用更多是把控制轨迹的自由度打开，Nc 从 2 变到 5 时响应形状会变化，但不会出现 PID 那种临界振荡效应。

5.3 液位耦合下的设定值跟踪与扰动实验

为了验证耦合抑制能力，建议在仿真 300 秒时给qin叠加入口扰动，比如模拟进水管压力下降，将 qin 直接减去 0.1 m³/s 持续 60 秒。观察两个液位的变化曲线。如果 GPC 调得合适，h2 的最大偏差应小于设定值阶跃偏差的 1/3，h1 因为加权较低会出现更大偏差，但不会发散。

实验中偶尔会遇到 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)*ny+1: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 接口里。

本文还有配套的精品资源，点击获取