简介广义预测控制算法程序压缩包面向自动化、控制工程领域的研究者与工程师聚焦多输入多输出及单输入单输出系统的预测控制实现。资源包含两个MATLAB脚本文件分别对应多变量与单变量系统的标准预测控制程序整体仅4KB代码精炼便于阅读和二次开发。广义预测控制采用基于模型的多步预测与滚动优化策略能够处理非线性、时变等复杂被控对象程序特意针对无扰动情况做了简化便于初学者抓住算法主线。已有337人学习下载是入门广义预测控制并对照算法流程的实用参考。通过运行这些程序读者可以直观理解预测模型构建、控制律计算、滚动优化等核心环节也可基于现有框架修改模型参数或扩展到带扰动、非线性场景为进一步研究和工程应用提供起点。1. 当 consSISO.m 和 consMIMO.m 同时出现时GPC 的约束边界就不只是“预测误差”解开 GPC.zip里面是两个 MATLAB 入口文件consSISO.m 和 consMIMO.m。这不是随手整理的一堆脚本而是广义预测控制Generalized Predictive Control, GPC可复现的标准骨架cons 前缀明确定义了控制增量和输出幅值的约束。它适合三类人做运动控制或过程控制、想把 PID 参数换成 MPC 框架的工程师需要可运行 baseline 做对比的研究者以及准备把老版本 Simulink 代码改造成独立控制器的维护者。所谓“无扰动情况”并不是为了省事而是为了让预测方程的因果边界先暴露出来——扰动补偿该放在哪个环节后续才有对比基准。2. 从 CARIMA 模型开始GPC 预测方程与 Diophantine 恒等式的数值落地2.1 为什么 GPC 必须用 CARIMA 而不是裸传递函数普通传递函数模型 A(z^{-1}) y(k) B(z^{-1}) u(k-1) 不携带积分因子。模型失配时预测残差会直接沉淀成稳态误差用 GPC 处理这类问题会在模型分母上强制乘上 Δ1-z^{-1}形成 CARIMA 模型A(z^{-1}) y(k) B(z^{-1}) u(k-1) C(z^{-1}) ξ(k)/ΔΔ 放进分母相当于在控制器内部嵌入一个积分器。即使像 GPC.zip 描述那样取 C1 且 ξ 是零均值白噪声控制律依然具备消除稳态误差的能力。换句话说“无扰动”不是把噪声项删掉而是让 C 和 Δ 的配合退化成最简单的积分预测结构方便单独验证预测矩阵和约束书写是否正确。从代码实现看CARIMA 把扰动模型直接放进被控对象模型里后面构造自由响应时不需要额外状态估计器。2.2 预测方程 YGΔUF 的展开标准 GPC 预测通过 Diophantine 恒等式 1E_j(z^{-1})A(z^{-1})Δz^{-j}F_j(z^{-1}) 把未来输出拆成两块一块由当前时刻以前已知的 y 和 Δu 决定叫自由响应 F一块由未来控制增量 ΔU 决定叫受迫响应 GΔU。最终预测方程是Y GΔU F其中 Y 是 Np×1 的预测输出向量ΔU 是 Nu×1 的控制增量向量G 是 Np×Nu 下三角 Toeplitz 矩阵。矩阵每一列就是模型对一次单位脉冲的响应系数。很多人纠结为什么叫 Diophantine 恒等式其实它只是在求一个因果分解把当前已知信息和未来未知输入分离。监视数值稳定性有个简单办法E_j 的长度随预测步 j 线性增长F_j 的长度保持在 max(na)-1如果递推过程出现跳变先检查是否把 C 多项式当作 1 时丢掉了扰动整形。2.3 用 MATLAB 递推方式求 G 矩阵最好不要把 tf2ss 的结果拿来构造预测矩阵。常见做法是直接递推 CARIMA 差分方程的脉冲响应系数再通过 tril 和 repmat 扩展成 Toeplitz 矩阵。我把这段函数称为“G 矩阵构造器”它也是 consSISO.m 里最值得单独抽出来测试的部分function [G, g] carima_step_response(A, B, Np) % A, B: CARIMA 模型多项式系数按 z^{-1} 升幂排列A(1)1 % 求单位脉冲响应前 Np 个系数并组装成预测矩阵 G Ad conv(A, [1, -1]); % A(z^{-1}) * (1-z^{-1}) na length(Ad) - 1; nb length(B) - 1; g zeros(1, Np); for i 1:Np past 0; for m 1:min(i-1, na) past past Ad(m1) * g(i-m); % 过去响应的加权和 end bi 0; if i-1 nb bi B(i); % 当前时刻激励的分子系数 end g(i) bi - past; end G tril(repmat(g, Np, 1)); % 下三角 Toeplitz G G(:, 1:Np); end代码里的 conv(A,[1,-1]) 完成 Δ 卷积得到的 Ad 是 CARIMA 分母。循环里的 past 用已经算出的 g 反推当前项这是差分方程的标准长除法。B(i) 只在 i-1nb 时才有效这是为了防止数组越界如果对象带纯延迟B 前几项为零G 前几行也会是零这属于正常现象。提示如果对象纯延迟大于一个采样周期预测矩阵会出现一整段零行QP 的 Hessian 会接近半正定。工程上先估计最小纯滞后把预测起点后移到首个非零响应处再进 QP比调整权重更干净。2.4 自由响应 F 的差分方程推进自由响应的含义是“假设未来 Δu 全为 0模型自己会怎么走”。它完全由当前和过去的输出、过去控制增量决定。我习惯用一个按时间倒序存历史值的函数做推进逐拍外推 Np 步function F free_response(A, B, y_past, du_past, Np) % y_past: 顺序是 y(k), y(k-1), ...; du_past: 顺序是 du(k-1), du(k-2), ... % 返回未来 Np 步自由响应列向量 Ad conv(A, [1, -1]); na length(Ad) - 1; nb length(B) - 1; F zeros(Np, 1); y_hist zeros(na, 1); y_hist(1:min(na, length(y_past))) y_past(1:min(na, length(y_past))); for i 1:Np y_next -Ad(2:end) * y_hist; % 历史输出对当前步的贡献 for r 1:nb1 idx i - r; % 这个系数对应增量 du(kidx) if idx 0 past_idx -idx; if past_idx length(du_past) y_next y_next B(r) * du_past(past_idx); end end end F(i) y_next; y_hist [y_next; y_hist(1:na-1)]; % 滚动更新输出历史 end endy_next 先用分母多项式对历史输出加权再对过去控制增量逐项累加。idxi-r 判断的是“第 r 个分子系数对应哪一拍的增量”idx0 表示这一拍还在过去idx0 表示当前拍idx0 表示未来后两者在自由响应里都取零所以直接跳过。这样得到的 F 和 G 放在一起才能构成完整的预测方程。表 1 列出需要核对的矩阵维度调试时按这个表逐步加断点。矩阵/向量维度来源GNp×Nucarima_step_responseFNp×1free_responseΔUNu×1QP 输出YNp×1GΔUFWNp×1设定点轨迹3. 约束 SISO GPC 的 QP 构造从解析解到 consSISO.m 的核心循环3.1 无约束解析解与约束问题的分界线当没有输入输出约束时最小化 J(Y-W)^T(Y-W)λΔU^TΔU 有解析解ΔU (G^T G λI)^{-1} G^T (W-F)consSISO.m 和大部分教学代码一样会先写这个解再用 quadprog 替代。边界在哪里解析解只有在所有约束都处于非激活状态时成立一旦执行器饱和或输出越限解析解就不再是最优解。因此看 consSISO.m第一件事不是读滚动优化循环而是看它对 quadprog 输入的 Hessian 和线性项怎么写这决定了约束是否真正生效。3.2 把执行器饱和和输出约束写成线性不等式SISO 系统的约束一共三类控制增量约束、控制量约束、输出约束。控制增量直接就是 ΔU 本身控制量约束用下三角矩阵 L 把 ΔU 累加成真实控制量输出约束则必须经过预测矩阵 G 转换因为输出是 ΔU 的线性函数。把它们堆成一个不等式组L ΔU ≤ U_max - u(k-1)-L ΔU ≤ u(k-1) - U_min G ΔU ≤ Y_max - F-G ΔU ≤ F - Y_min。增量约束不需要 L。写代码时最容易错的是符号方向quadprog 约束一律是 Aineq·x ≤ bineq所以“≥”形式的约束要乘负号再放进去。下面这段代码可以直接对齐 consSISO.m 的 QP 构造部分。3.3 一个可直接对准 consSISO.m 的 QP 求解骨架function du solve_gpc_qp(G, F, W, lambda, u_prev, u_min, u_max, ... du_min, du_max, y_min, y_max) Np size(G, 1); Nu size(G, 2); H 2 * (G * G lambda * eye(Nu)); f 2 * G * (F - W); % quadprog 目标: 0.5*x*H*x f*x L tril(ones(Nu)); % 累积矩阵 Aineq [ eye(Nu); -eye(Nu); ... L; -L; ... G; -G ]; bineq [ du_max*ones(Nu,1); -du_min*ones(Nu,1); ... (u_max-u_prev)*ones(Nu,1); (u_prev-u_min)*ones(Nu,1); ... y_max-F; F-y_min ]; options optimoptions(quadprog, Display, off, ... Algorithm, interior-point-convex); du quadprog(H, f, Aineq, bineq, [], [], [], [], [], options); if isempty(du) error(GPC QP infeasible: 检查 y_min/y_max 是否过紧); end endH 前乘 2 是为了配合 quadprog 的二次项定义。f 写成 2G(F-W) 而不是 (F-W)*G目的是避免产生行向量导致维度报错。bineq 里的控制量约束用 ones(Nu,1) 扩展成向量因为每一个未来控制步都必须落在同一个上下限之间。输出约束直接用 y_max-F 和 F-y_min这是从 YGΔUF 推出来的。只要看到 quadprog 报 dimension mismatch优先检查 G 的列数确认它设置成了 Nu 而不是 Np。3.4 约束饱和判断与滚动优化入口consSISO.m 的主循环按“预测 → 优化 → 取第一步 → 更新状态”的顺序推进。每次只取 du(1)因为下一个采样周期要用新的量测重新预测这就是滚动优化的核心。约束是否真的参与了解可以用一个非常简单的触发器判断如果 du 落在 du_max 或 du_min 边界上或与边界之差小于 1e-6说明 QP 的约束集里有成员被激活。表 2 给出各变量的维度关系便于和工作区逐项比对。QP 输入说明常见错误HNu×Nu对称正定忘记乘 2fNu×1线性项转置导致维度错误Aineq(2Nu2Nu2Np)×Nu约束顺序与 bineq 不一致bineq(2Nu2Nu2Np)×1标量直接放进来未扩展整个 SISO 部分算清楚以后MIMO 只是把这里的 G、Q、R 变成分块矩阵QP 求解器的调用方式可以保持不变。4. MIMO 广义预测控制consMIMO.m 的多变量解耦与分块预测4.1 输入输出配对与分块预测矩阵MIMO 广义预测控制的第一步不是解耦而是把多输入多输出关系组织成“通道对通道”的模型网格。假设被控对象有 nu 个输入、ny 个输出每个输出 y_i 都受所有输入 u_j 影响。对每一对 (i,j) 构建一个 CARIMA 模型预测方程为Y_i Σ_{j1}^{nu} G_{ij} ΔU_j F_i把所有输出通道竖向堆叠得到整体预测方程。这样 QP 的决策变量仍然是所有输入增量只是 G 从 Np×Nu 变成 (ny·Np)×(nu·Nu)。这里最容易被忽略的是G 的列顺序必须和 ΔU 的排列顺序一致否则约束张量写对了也会解出错误通道的控制量。4.2 consMIMO.m 的分块矩阵构建用一个二维结构体存通道模型再循环调用单变量函数就可以把 SISO 代码直接升级为 MIMO 预测矩阵。这段逻辑在代码里通常体现为“块 Toeplitz 拼接”function Gt mimo_prediction_matrix(mods, ny, nu, Np, Nu) % mods{i,j}.A 和 mods{i,j}.B: 输出 i 对输入 j 的 CARIMA 模型 % Gt 的尺寸: (ny*Np) x (nu*Nu) Gt zeros(ny*Np, nu*Nu); for i 1:ny for j 1:nu [Gij, ~] carima_step_response(mods{i,j}.A, mods{i,j}.B, Np); row (i-1)*Np 1 : i*Np; col (j-1)*Nu 1 : j*Nu; Gt(row, col) Gij(:, 1:Nu); end end end外层循环遍历输出 i内层遍历输入 j。row 和 col 的作用是把单通道的 Np×Nu 矩阵放进整个分块矩阵的对应位置。Gij(:,1:Nu) 只取前 Nu 列因为预测时域 Np 一般大于控制时域 Nu后段 Δu 已经固定为 0不应该再出现在决策变量里。这个细节如果写错MIMO 代码会在第一个采样周期就把控制器输出跳到上限因为 QP 会错误地认为未来控制增量还能继续变化。4.3 权重矩阵与约束张量的分块堆叠MIMO 的目标函数需要把每个输出通道的权重独立设置。常见做法是用 Kronecker 积构造块对角权重qy [1; 2; 0.8]; % 每个输出通道的权重 Qu kron(diag(qy), eye(Np)); % 块对角: 每个通道一个 Np×Np 块 R lambda * kron(eye(nu), eye(Nu)); H 2 * (Gt * Qu * Gt R); f 2 * Gt * Qu * (Ft - Wt);Qu 的维度是 (ny·Np)×(ny·Np)它把不同输出通道按重要程度放大误差R 的维度是 (nu·Nu)×(nu·Nu)控制每个输入通道的增量变化速度。Ft 和 Wt 分别是把每个通道自由响应和设定点轨迹按输出顺序堆叠出的列向量。表 3 总结了块矩阵的尺寸来源按这张表能快速发现 MIMO 代码里最常见的问题某一行忘了乘 ny导致二维索引错位。矩阵维度来源Gt(ny·Np)×(nu·Nu)每个子块循环调用Ft(ny·Np)×1各输出通道自由响应堆叠Wt(ny·Np)×1各输出通道参考轨迹Qu(ny·Np)×(ny·Np)kron(diag(qy), eye(Np))R(nu·Nu)×(nu·Nu)kron(eye(nu), lambda·I)4.4 无扰动假定在多变量场景下的工程代价GPC.zip 中的两个文件都标记了“无扰动情况”。在 SISO 里这意味着 C1自由响应不需要白噪声估计在 MIMO 里代价会被交互耦合放大。任何一个通道模型失配预测误差都会通过 G_{ij} 的影响传递到所有输出通道如果自由响应不做修正闭环会呈现出“看似稳定但总是差一点”的残差。工程上最常见的补救是在每个采样周期计算输出预测残差 d_i y_i(k) - ŷ_i(k)然后直接把 d_i 加到该通道未来 Np 步的自由响应上。这一步不属于基础 GPC却是从 consMIMO.m 走向真实对象时最值得补的功能。5. 参数整定与约束激活验证跑闭环前先把这两组测试做完5.1 四个参数的建议起点预测时域 Np、控制时域 Nu、控制加权 lambda、约束边界这四者的影响并不独立。表 4 是我对这类程序一直采用的起点和调法。参数推荐起点系统表现调整方向Np4~7 倍上升时间响应过早/超调增大 NpNumax(5, Np/3)计算量过大减小 Nulambda0.1控制量抖动增大 lambday_min/y_max物理执行边界求解失败先放松再逐步回紧Np 太小会让预测看不到系统主惯性约束解会过于激进Np 太大虽然稳定但计算量随 (ny·Np)×(nu·Nu) 快速膨胀。lambda 不是越大越好因为它不对输出误差起作用只能压低增量变化会在约束边界附近使控制器变得“迟钝”。5.2 确定性验证脚本把对象设成一阶惯性模型Np30、Nu10、lambda0.1跑一个 200 步的阶跃响应观察 du 是否在 du_max 处出现平台。出现平台说明约束被激活这是正常的不是数值问题。A [1, -0.85]; B [0, 0.15]; % 对象真实模型: y(k)0.85y(k-1)0.15u(k-1) Np 30; Nu 10; lambda 0.1; du_max 0.2; y_max 1.2; y_min -1.2; u_prev 0; y_past 0; du_past 0; for k 1:200 [G, ~] carima_step_response(A, B, Np); F free_response(A, B, y_past, du_past, Np); du solve_gpc_qp(G, F, 1.0*ones(Np,1), lambda, u_prev, ... -0.5, 0.5, -du_max, du_max, y_min, y_max); u_curr u_prev du(1); y_curr -A(2) * y_past(1) B(2) * u_prev; % 真实对象差分方程 u_past [u_curr; u_prev]; y_past [y_curr; y_past]; du_past [du(1); du_past]; u_prev u_curr; end这段脚本里的 free_response 接收的 du_past 顺序是 du(k-1), du(k-2)……与实现保持一致。solve_gpc_qp 的返回值在约束激活时会正好落在 du_max 或 y_max 边界因此可以把 du(1) 直接打印出来肉眼观察是否存在越过边界的瞬间。5.3 约束激活的三个检查点增量子落在边界上说明 QP 的增量约束生效如果所有步都落在边界说明 lambda 太小或 Np 太短。输出轨迹贴着 y_max 走说明输出约束生效若轨迹在边界附近来回振荡需要给输出约束增加软化带而不是继续放大 QP 里的松弛变量。quadprog 返回空解绝大多数不是算法问题而是 y_min/y_max 设置得比执行机构物理能力更紧。先把边界放宽到执行器极限验证控制器内核再逐步逼近真实约束。最后一个检查点是最容易被忽视的“逻辑约束”和“物理约束”之间的冲突。调节器内部用数学边界做 QP被控对象用物理边界做硬件保护两者相差过大会直接导致 infeasible。把 y_min/y_max 先放回执行机构的物理极限确认闭环能跑通再把约束一步步收紧这个过程比任何调参公式都可靠。本文还有配套的精品资源点击获取 GPC广义预测控制核心解析：CARIMA模型、预测矩阵与约束QP实现
简介广义预测控制算法程序压缩包面向自动化、控制工程领域的研究者与工程师聚焦多输入多输出及单输入单输出系统的预测控制实现。资源包含两个MATLAB脚本文件分别对应多变量与单变量系统的标准预测控制程序整体仅4KB代码精炼便于阅读和二次开发。广义预测控制采用基于模型的多步预测与滚动优化策略能够处理非线性、时变等复杂被控对象程序特意针对无扰动情况做了简化便于初学者抓住算法主线。已有337人学习下载是入门广义预测控制并对照算法流程的实用参考。通过运行这些程序读者可以直观理解预测模型构建、控制律计算、滚动优化等核心环节也可基于现有框架修改模型参数或扩展到带扰动、非线性场景为进一步研究和工程应用提供起点。1. 当 consSISO.m 和 consMIMO.m 同时出现时GPC 的约束边界就不只是“预测误差”解开 GPC.zip里面是两个 MATLAB 入口文件consSISO.m 和 consMIMO.m。这不是随手整理的一堆脚本而是广义预测控制Generalized Predictive Control, GPC可复现的标准骨架cons 前缀明确定义了控制增量和输出幅值的约束。它适合三类人做运动控制或过程控制、想把 PID 参数换成 MPC 框架的工程师需要可运行 baseline 做对比的研究者以及准备把老版本 Simulink 代码改造成独立控制器的维护者。所谓“无扰动情况”并不是为了省事而是为了让预测方程的因果边界先暴露出来——扰动补偿该放在哪个环节后续才有对比基准。2. 从 CARIMA 模型开始GPC 预测方程与 Diophantine 恒等式的数值落地2.1 为什么 GPC 必须用 CARIMA 而不是裸传递函数普通传递函数模型 A(z^{-1}) y(k) B(z^{-1}) u(k-1) 不携带积分因子。模型失配时预测残差会直接沉淀成稳态误差用 GPC 处理这类问题会在模型分母上强制乘上 Δ1-z^{-1}形成 CARIMA 模型A(z^{-1}) y(k) B(z^{-1}) u(k-1) C(z^{-1}) ξ(k)/ΔΔ 放进分母相当于在控制器内部嵌入一个积分器。即使像 GPC.zip 描述那样取 C1 且 ξ 是零均值白噪声控制律依然具备消除稳态误差的能力。换句话说“无扰动”不是把噪声项删掉而是让 C 和 Δ 的配合退化成最简单的积分预测结构方便单独验证预测矩阵和约束书写是否正确。从代码实现看CARIMA 把扰动模型直接放进被控对象模型里后面构造自由响应时不需要额外状态估计器。2.2 预测方程 YGΔUF 的展开标准 GPC 预测通过 Diophantine 恒等式 1E_j(z^{-1})A(z^{-1})Δz^{-j}F_j(z^{-1}) 把未来输出拆成两块一块由当前时刻以前已知的 y 和 Δu 决定叫自由响应 F一块由未来控制增量 ΔU 决定叫受迫响应 GΔU。最终预测方程是Y GΔU F其中 Y 是 Np×1 的预测输出向量ΔU 是 Nu×1 的控制增量向量G 是 Np×Nu 下三角 Toeplitz 矩阵。矩阵每一列就是模型对一次单位脉冲的响应系数。很多人纠结为什么叫 Diophantine 恒等式其实它只是在求一个因果分解把当前已知信息和未来未知输入分离。监视数值稳定性有个简单办法E_j 的长度随预测步 j 线性增长F_j 的长度保持在 max(na)-1如果递推过程出现跳变先检查是否把 C 多项式当作 1 时丢掉了扰动整形。2.3 用 MATLAB 递推方式求 G 矩阵最好不要把 tf2ss 的结果拿来构造预测矩阵。常见做法是直接递推 CARIMA 差分方程的脉冲响应系数再通过 tril 和 repmat 扩展成 Toeplitz 矩阵。我把这段函数称为“G 矩阵构造器”它也是 consSISO.m 里最值得单独抽出来测试的部分function [G, g] carima_step_response(A, B, Np) % A, B: CARIMA 模型多项式系数按 z^{-1} 升幂排列A(1)1 % 求单位脉冲响应前 Np 个系数并组装成预测矩阵 G Ad conv(A, [1, -1]); % A(z^{-1}) * (1-z^{-1}) na length(Ad) - 1; nb length(B) - 1; g zeros(1, Np); for i 1:Np past 0; for m 1:min(i-1, na) past past Ad(m1) * g(i-m); % 过去响应的加权和 end bi 0; if i-1 nb bi B(i); % 当前时刻激励的分子系数 end g(i) bi - past; end G tril(repmat(g, Np, 1)); % 下三角 Toeplitz G G(:, 1:Np); end代码里的 conv(A,[1,-1]) 完成 Δ 卷积得到的 Ad 是 CARIMA 分母。循环里的 past 用已经算出的 g 反推当前项这是差分方程的标准长除法。B(i) 只在 i-1nb 时才有效这是为了防止数组越界如果对象带纯延迟B 前几项为零G 前几行也会是零这属于正常现象。提示如果对象纯延迟大于一个采样周期预测矩阵会出现一整段零行QP 的 Hessian 会接近半正定。工程上先估计最小纯滞后把预测起点后移到首个非零响应处再进 QP比调整权重更干净。2.4 自由响应 F 的差分方程推进自由响应的含义是“假设未来 Δu 全为 0模型自己会怎么走”。它完全由当前和过去的输出、过去控制增量决定。我习惯用一个按时间倒序存历史值的函数做推进逐拍外推 Np 步function F free_response(A, B, y_past, du_past, Np) % y_past: 顺序是 y(k), y(k-1), ...; du_past: 顺序是 du(k-1), du(k-2), ... % 返回未来 Np 步自由响应列向量 Ad conv(A, [1, -1]); na length(Ad) - 1; nb length(B) - 1; F zeros(Np, 1); y_hist zeros(na, 1); y_hist(1:min(na, length(y_past))) y_past(1:min(na, length(y_past))); for i 1:Np y_next -Ad(2:end) * y_hist; % 历史输出对当前步的贡献 for r 1:nb1 idx i - r; % 这个系数对应增量 du(kidx) if idx 0 past_idx -idx; if past_idx length(du_past) y_next y_next B(r) * du_past(past_idx); end end end F(i) y_next; y_hist [y_next; y_hist(1:na-1)]; % 滚动更新输出历史 end endy_next 先用分母多项式对历史输出加权再对过去控制增量逐项累加。idxi-r 判断的是“第 r 个分子系数对应哪一拍的增量”idx0 表示这一拍还在过去idx0 表示当前拍idx0 表示未来后两者在自由响应里都取零所以直接跳过。这样得到的 F 和 G 放在一起才能构成完整的预测方程。表 1 列出需要核对的矩阵维度调试时按这个表逐步加断点。矩阵/向量维度来源GNp×Nucarima_step_responseFNp×1free_responseΔUNu×1QP 输出YNp×1GΔUFWNp×1设定点轨迹3. 约束 SISO GPC 的 QP 构造从解析解到 consSISO.m 的核心循环3.1 无约束解析解与约束问题的分界线当没有输入输出约束时最小化 J(Y-W)^T(Y-W)λΔU^TΔU 有解析解ΔU (G^T G λI)^{-1} G^T (W-F)consSISO.m 和大部分教学代码一样会先写这个解再用 quadprog 替代。边界在哪里解析解只有在所有约束都处于非激活状态时成立一旦执行器饱和或输出越限解析解就不再是最优解。因此看 consSISO.m第一件事不是读滚动优化循环而是看它对 quadprog 输入的 Hessian 和线性项怎么写这决定了约束是否真正生效。3.2 把执行器饱和和输出约束写成线性不等式SISO 系统的约束一共三类控制增量约束、控制量约束、输出约束。控制增量直接就是 ΔU 本身控制量约束用下三角矩阵 L 把 ΔU 累加成真实控制量输出约束则必须经过预测矩阵 G 转换因为输出是 ΔU 的线性函数。把它们堆成一个不等式组L ΔU ≤ U_max - u(k-1)-L ΔU ≤ u(k-1) - U_min G ΔU ≤ Y_max - F-G ΔU ≤ F - Y_min。增量约束不需要 L。写代码时最容易错的是符号方向quadprog 约束一律是 Aineq·x ≤ bineq所以“≥”形式的约束要乘负号再放进去。下面这段代码可以直接对齐 consSISO.m 的 QP 构造部分。3.3 一个可直接对准 consSISO.m 的 QP 求解骨架function du solve_gpc_qp(G, F, W, lambda, u_prev, u_min, u_max, ... du_min, du_max, y_min, y_max) Np size(G, 1); Nu size(G, 2); H 2 * (G * G lambda * eye(Nu)); f 2 * G * (F - W); % quadprog 目标: 0.5*x*H*x f*x L tril(ones(Nu)); % 累积矩阵 Aineq [ eye(Nu); -eye(Nu); ... L; -L; ... G; -G ]; bineq [ du_max*ones(Nu,1); -du_min*ones(Nu,1); ... (u_max-u_prev)*ones(Nu,1); (u_prev-u_min)*ones(Nu,1); ... y_max-F; F-y_min ]; options optimoptions(quadprog, Display, off, ... Algorithm, interior-point-convex); du quadprog(H, f, Aineq, bineq, [], [], [], [], [], options); if isempty(du) error(GPC QP infeasible: 检查 y_min/y_max 是否过紧); end endH 前乘 2 是为了配合 quadprog 的二次项定义。f 写成 2G(F-W) 而不是 (F-W)*G目的是避免产生行向量导致维度报错。bineq 里的控制量约束用 ones(Nu,1) 扩展成向量因为每一个未来控制步都必须落在同一个上下限之间。输出约束直接用 y_max-F 和 F-y_min这是从 YGΔUF 推出来的。只要看到 quadprog 报 dimension mismatch优先检查 G 的列数确认它设置成了 Nu 而不是 Np。3.4 约束饱和判断与滚动优化入口consSISO.m 的主循环按“预测 → 优化 → 取第一步 → 更新状态”的顺序推进。每次只取 du(1)因为下一个采样周期要用新的量测重新预测这就是滚动优化的核心。约束是否真的参与了解可以用一个非常简单的触发器判断如果 du 落在 du_max 或 du_min 边界上或与边界之差小于 1e-6说明 QP 的约束集里有成员被激活。表 2 给出各变量的维度关系便于和工作区逐项比对。QP 输入说明常见错误HNu×Nu对称正定忘记乘 2fNu×1线性项转置导致维度错误Aineq(2Nu2Nu2Np)×Nu约束顺序与 bineq 不一致bineq(2Nu2Nu2Np)×1标量直接放进来未扩展整个 SISO 部分算清楚以后MIMO 只是把这里的 G、Q、R 变成分块矩阵QP 求解器的调用方式可以保持不变。4. MIMO 广义预测控制consMIMO.m 的多变量解耦与分块预测4.1 输入输出配对与分块预测矩阵MIMO 广义预测控制的第一步不是解耦而是把多输入多输出关系组织成“通道对通道”的模型网格。假设被控对象有 nu 个输入、ny 个输出每个输出 y_i 都受所有输入 u_j 影响。对每一对 (i,j) 构建一个 CARIMA 模型预测方程为Y_i Σ_{j1}^{nu} G_{ij} ΔU_j F_i把所有输出通道竖向堆叠得到整体预测方程。这样 QP 的决策变量仍然是所有输入增量只是 G 从 Np×Nu 变成 (ny·Np)×(nu·Nu)。这里最容易被忽略的是G 的列顺序必须和 ΔU 的排列顺序一致否则约束张量写对了也会解出错误通道的控制量。4.2 consMIMO.m 的分块矩阵构建用一个二维结构体存通道模型再循环调用单变量函数就可以把 SISO 代码直接升级为 MIMO 预测矩阵。这段逻辑在代码里通常体现为“块 Toeplitz 拼接”function Gt mimo_prediction_matrix(mods, ny, nu, Np, Nu) % mods{i,j}.A 和 mods{i,j}.B: 输出 i 对输入 j 的 CARIMA 模型 % Gt 的尺寸: (ny*Np) x (nu*Nu) Gt zeros(ny*Np, nu*Nu); for i 1:ny for j 1:nu [Gij, ~] carima_step_response(mods{i,j}.A, mods{i,j}.B, Np); row (i-1)*Np 1 : i*Np; col (j-1)*Nu 1 : j*Nu; Gt(row, col) Gij(:, 1:Nu); end end end外层循环遍历输出 i内层遍历输入 j。row 和 col 的作用是把单通道的 Np×Nu 矩阵放进整个分块矩阵的对应位置。Gij(:,1:Nu) 只取前 Nu 列因为预测时域 Np 一般大于控制时域 Nu后段 Δu 已经固定为 0不应该再出现在决策变量里。这个细节如果写错MIMO 代码会在第一个采样周期就把控制器输出跳到上限因为 QP 会错误地认为未来控制增量还能继续变化。4.3 权重矩阵与约束张量的分块堆叠MIMO 的目标函数需要把每个输出通道的权重独立设置。常见做法是用 Kronecker 积构造块对角权重qy [1; 2; 0.8]; % 每个输出通道的权重 Qu kron(diag(qy), eye(Np)); % 块对角: 每个通道一个 Np×Np 块 R lambda * kron(eye(nu), eye(Nu)); H 2 * (Gt * Qu * Gt R); f 2 * Gt * Qu * (Ft - Wt);Qu 的维度是 (ny·Np)×(ny·Np)它把不同输出通道按重要程度放大误差R 的维度是 (nu·Nu)×(nu·Nu)控制每个输入通道的增量变化速度。Ft 和 Wt 分别是把每个通道自由响应和设定点轨迹按输出顺序堆叠出的列向量。表 3 总结了块矩阵的尺寸来源按这张表能快速发现 MIMO 代码里最常见的问题某一行忘了乘 ny导致二维索引错位。矩阵维度来源Gt(ny·Np)×(nu·Nu)每个子块循环调用Ft(ny·Np)×1各输出通道自由响应堆叠Wt(ny·Np)×1各输出通道参考轨迹Qu(ny·Np)×(ny·Np)kron(diag(qy), eye(Np))R(nu·Nu)×(nu·Nu)kron(eye(nu), lambda·I)4.4 无扰动假定在多变量场景下的工程代价GPC.zip 中的两个文件都标记了“无扰动情况”。在 SISO 里这意味着 C1自由响应不需要白噪声估计在 MIMO 里代价会被交互耦合放大。任何一个通道模型失配预测误差都会通过 G_{ij} 的影响传递到所有输出通道如果自由响应不做修正闭环会呈现出“看似稳定但总是差一点”的残差。工程上最常见的补救是在每个采样周期计算输出预测残差 d_i y_i(k) - ŷ_i(k)然后直接把 d_i 加到该通道未来 Np 步的自由响应上。这一步不属于基础 GPC却是从 consMIMO.m 走向真实对象时最值得补的功能。5. 参数整定与约束激活验证跑闭环前先把这两组测试做完5.1 四个参数的建议起点预测时域 Np、控制时域 Nu、控制加权 lambda、约束边界这四者的影响并不独立。表 4 是我对这类程序一直采用的起点和调法。参数推荐起点系统表现调整方向Np4~7 倍上升时间响应过早/超调增大 NpNumax(5, Np/3)计算量过大减小 Nulambda0.1控制量抖动增大 lambday_min/y_max物理执行边界求解失败先放松再逐步回紧Np 太小会让预测看不到系统主惯性约束解会过于激进Np 太大虽然稳定但计算量随 (ny·Np)×(nu·Nu) 快速膨胀。lambda 不是越大越好因为它不对输出误差起作用只能压低增量变化会在约束边界附近使控制器变得“迟钝”。5.2 确定性验证脚本把对象设成一阶惯性模型Np30、Nu10、lambda0.1跑一个 200 步的阶跃响应观察 du 是否在 du_max 处出现平台。出现平台说明约束被激活这是正常的不是数值问题。A [1, -0.85]; B [0, 0.15]; % 对象真实模型: y(k)0.85y(k-1)0.15u(k-1) Np 30; Nu 10; lambda 0.1; du_max 0.2; y_max 1.2; y_min -1.2; u_prev 0; y_past 0; du_past 0; for k 1:200 [G, ~] carima_step_response(A, B, Np); F free_response(A, B, y_past, du_past, Np); du solve_gpc_qp(G, F, 1.0*ones(Np,1), lambda, u_prev, ... -0.5, 0.5, -du_max, du_max, y_min, y_max); u_curr u_prev du(1); y_curr -A(2) * y_past(1) B(2) * u_prev; % 真实对象差分方程 u_past [u_curr; u_prev]; y_past [y_curr; y_past]; du_past [du(1); du_past]; u_prev u_curr; end这段脚本里的 free_response 接收的 du_past 顺序是 du(k-1), du(k-2)……与实现保持一致。solve_gpc_qp 的返回值在约束激活时会正好落在 du_max 或 y_max 边界因此可以把 du(1) 直接打印出来肉眼观察是否存在越过边界的瞬间。5.3 约束激活的三个检查点增量子落在边界上说明 QP 的增量约束生效如果所有步都落在边界说明 lambda 太小或 Np 太短。输出轨迹贴着 y_max 走说明输出约束生效若轨迹在边界附近来回振荡需要给输出约束增加软化带而不是继续放大 QP 里的松弛变量。quadprog 返回空解绝大多数不是算法问题而是 y_min/y_max 设置得比执行机构物理能力更紧。先把边界放宽到执行器极限验证控制器内核再逐步逼近真实约束。最后一个检查点是最容易被忽视的“逻辑约束”和“物理约束”之间的冲突。调节器内部用数学边界做 QP被控对象用物理边界做硬件保护两者相差过大会直接导致 infeasible。把 y_min/y_max 先放回执行机构的物理极限确认闭环能跑通再把约束一步步收紧这个过程比任何调参公式都可靠。本文还有配套的精品资源点击获取