预测控制 · 经典论文精读

广义预测控制 Part I — 基础算法

Generalized Predictive Control — Part I. The Basic Algorithm

原著：D. W. Clarke, C. Mohtadi, P. S. Tuffs · Automatica, Vol. 23, No. 2, 1987

笔记日期：2026-06-06 关键词：GPC, CARIMA, Diophantine, 预测控制, 自适应控制, 滚动时域

§0摘要

现有自整定算法对死时间和模型阶次的先验选择缺乏鲁棒性。本文提出一种新方法——广义预测控制（GPC），仿真研究表明其优于广义最小方差（GMV）和极点配置等公认技术。

该滚动时域方法通过对未来控制动作的假设，在多步内预测被控对象输出。其核心假设是：在"控制域"之后所有控制增量置零——这对鲁棒性和简化计算均有益。通过选取特定的输出域与控制域参数，GPC 可退化为 GMV、EPSAC、Peterka 预测控制器（1984）和 Ydstie 扩展域设计（1984）等子集。

核心能力

GPC 可控制：非最小相位对象、开环不稳定对象、可变/未知死时间对象、模型过参数化对象，且消除稳态偏差（CARIMA 模型的积分特性）。

§1引言

一个通用自适应控制算法需能处理以下四类挑战：

非最小相位对象：连续时间传递函数在足够快采样速率下会产生单位圆外的离散零点。
开环不稳定或阻尼不足对象：如柔性航天器、机器人。
可变/未知死时间：最小方差自整定器对死时间假设高度敏感。
未知阶次对象：极点配置和 LQG 自整定器在模型阶次过估计时因极零点对消而性能恶化。

GPC 在单一算法中克服了这些问题，可在参数、死时间和模型阶次均发生变化时稳定控制被控对象，对同时具有非最小相位和开环不稳定特性的对象亦有效，且无需特殊预防措施处理过参数化估计。

工业过程扰动通常表现为随机时刻的随机阶跃（确定性情形）或布朗运动（随机情形）。GPC 通过 CARIMA 模型的积分特性自然地具有积分作用，无需人为添加，从而实现无偏差调节。

§2CARIMA 被控对象模型与输出预测

在某操作点附近，非线性被控对象通常可以局部线性化为：

$$A(q^{-1})y(t) = B(q^{-1})u(t-1) + x(t) \tag{1}$$

其中 $q^{-1}$ 为后移算子，$A$、$B$ 为多项式：

$$A(q^{-1}) = 1 + a_1 q^{-1} + \cdots + a_{n_a} q^{-n_a}, \quad B(q^{-1}) = b_0 + b_1 q^{-1} + \cdots + b_{n_b} q^{-n_b}$$

$A$ 为首一多项式；若 $b_0 = 0$，则对象含纯滞后（dead-time）。

工业中扰动常为非平稳型，适合建模为 $x(t) = C(q^{-1})\xi(t)/\Delta$，其中 $\Delta = 1-q^{-1}$。合并后得 CARIMA 模型：

$$A(q^{-1})y(t) = B(q^{-1})u(t-1) + C(q^{-1})\xi(t)/\Delta$$

为什么 CARIMA 而非 CARMA？

CARMA 适合平稳扰动；工业扰动多为非平稳型。CARIMA 中 $1/\Delta$ 引入积分，使控制器自然具备积分作用，消除稳态偏差，无需人为添加积分环节。

为简化推导，取 $C(q^{-1}) = 1$，得到本文基础模型：

$$A(q^{-1})y(t) = B(q^{-1})u(t-1) + \xi(t)/\Delta \tag{5}$$

等价写法：$\Delta A(q^{-1})y(t) = B(q^{-1})\Delta u(t-1) + \xi(t)$，控制增量 $\Delta u$ 自然出现。

推导 $j$ 步最优预测器

引入 Diophantine 恒等式：对每个预测步长 $j$，存在唯一多项式对 $(E_j, F_j)$ 满足：

$$1 = E_j(q^{-1})\Delta A(q^{-1}) + q^{-j}F_j(q^{-1}) \tag{6}$$

其中 $E_j$ 阶次为 $j-1$，$F_j$ 阶次为 $n_a$。

Diophantine 方程的作用

式(6)将"1"分解为两部分：$E_j \Delta A$（已知量相关）和 $q^{-j}F_j$（延迟项）。这一分解是提取最优预测器的关键，使未来噪声项与已知量完全分离。

将式(5)乘以 $E_j \Delta q^j$ 并代入式(6)，得 $j$ 步最优预测器：

$$\hat{y}(t+j \mid t) = G_j(q^{-1})\Delta u(t+j-1) + F_j(q^{-1})y(t) \tag{8}$$

其中 $G_j(q^{-1}) = E_j(q^{-1})B(q^{-1})$。由于 $E_j$ 阶次为 $j-1$，噪声项 $E_j\xi(t+j)$ 只含 $\xi(t+1),\ldots,\xi(t+j)$（全为未来量），期望为零。$G_j$ 的系数恰为对象阶跃响应的前 $j$ 项，与 $j$ 无关：$g_{ij} = g_i$（$i < j$）。

2.1 Diophantine 方程的递推计算

对每个 $j$ 分别数值求解 Diophantine 方程代价过高。更高效的方案是利用相邻步长间的递推关系：给定 $(E_j, F_j)$ 直接得到 $(E_{j+1}, F_{j+1})$。

设 $\tilde{A} = \Delta A$，两步 Diophantine 方程相减并分析多项式结构，得递推关系：

$$r_j = f_0 \tag{11a}$$

$$s_i = f_{i+1} - \tilde{a}_{i+1} r_j \quad (i = 0, 1, \ldots) \tag{11b}$$

$$E_{j+1}(q^{-1}) = E_j(q^{-1}) + q^{-j} r_j \tag{12}$$

$$G_{j+1}(q^{-1}) = B(q^{-1}) E_{j+1}(q^{-1}) \tag{13}$$

递推含义逐条解释

式(11a)：递推系数 $r_j$ = 当前 $F_j$ 的常数项 $f_0$（因 $\tilde{A}$ 首项为1）
式(11b)：$F_{j+1}$ 各系数由 $F_j$ 系数向右移位一步并减去修正量 $\tilde{a}_{i+1}r_j$ 得到
式(12)：$E_{j+1}$ = $E_j$ 末尾追加一项 $r_j q^{-j}$，长度每步增加1
式(13)：下一步预测响应多项式直接由 $B \cdot E_{j+1}$ 给出

整个递推仅需少量加减乘法，无需重解 Diophantine 方程，对自适应版本尤为重要。

初始化（$j=1$）：由 $1 = E_1\tilde{A} + q^{-1}F_1$ 且 $\tilde{A}$ 首项为1，得 $E_1 = 1$，$F_1(q^{-1}) = q(1 - \tilde{A}(q^{-1}))$。

§3预测控制律

设未来设定值序列 $\{w(t+j)\}$ 已知（通常取当前设定值 $w$ 的常数），滚动时域控制律每个采样步执行以下四步：

计算未来设定值序列 $w(t+j)$
用预测模型(8)生成预测输出 $\hat{y}(t+j|t)$，注意 $j>k$ 时依赖于待确定的未来控制信号
最小化未来误差和控制量的二次代价函数，得到未来控制序列
施加第一个控制量 $u(t)$，下一采样步重复上述过程（滚动时域）

代价函数为：

$$J(N_1, N_2, N_U, \lambda) = \sum_{j=N_1}^{N_2}\left[\hat{y}(t+j|t)-w(t+j)\right]^2 + \sum_{j=1}^{N_U}\lambda(j)\left[\Delta u(t+j-1)\right]^2 \tag{14}$$

其中 $N_1$ 为最小输出域，$N_2$ 为最大输出域，$\lambda(j)$ 为控制加权序列。第一项惩罚输出跟踪误差，第二项惩罚控制增量幅度，平衡跟踪精度与控制能量。

取 $N_1=1$，$N_2=N$，$\lambda(j)=\lambda$（常数），将 $N$ 步预测写成向量形式：

$$\hat{\mathbf{y}} = G\tilde{\mathbf{u}} + \mathbf{f} \tag{15}$$

其中 $\hat{\mathbf{y}} = [\hat{y}(t+1),\ldots,\hat{y}(t+N)]^T$，$\tilde{\mathbf{u}} = [\Delta u(t),\ldots,\Delta u(t+N-1)]^T$，$\mathbf{f}$ 为 $t$ 时刻已知量（历史控制和输出）决定的"自由响应"向量。

矩阵 $G$ 是下三角 Toeplitz 矩阵（阶跃响应系数构成）：

$$G = \begin{pmatrix} g_0 & 0 & \cdots & 0 \\ g_1 & g_0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ g_{N-1} & g_{N-2} & \cdots & g_0 \end{pmatrix}$$

下三角结构反映因果性：$t+j$ 时刻的输出只受 $\Delta u(t),\ldots,\Delta u(t+j-1)$ 影响。若死时间为 $k>1$，$G$ 前 $k-1$ 行全为零，但 GPC 在此情形下仍能给出稳定解（这是 GPC 对未知死时间鲁棒的关键）。

代价函数展开为：

$$J_1 = (G\tilde{\mathbf{u}} + \mathbf{f} - \mathbf{w})^T(G\tilde{\mathbf{u}} + \mathbf{f} - \mathbf{w}) + \lambda\,\tilde{\mathbf{u}}^T\tilde{\mathbf{u}} \tag{16}$$

对 $\tilde{\mathbf{u}}$ 求偏导令其为零（$\partial J_1/\partial\tilde{\mathbf{u}} = 0$）：

$$2G^T(G\tilde{\mathbf{u}} + \mathbf{f} - \mathbf{w}) + 2\lambda\tilde{\mathbf{u}} = 0$$

解得 GPC 最优控制增量向量：

$$\hat{\tilde{\mathbf{u}}} = (G^T G + \lambda I)^{-1} G^T (\mathbf{w} - \mathbf{f}) \tag{17}$$

式(17)关键性质

这是标准正则化最小二乘解；$(G^TG + \lambda I)$ 正定，解唯一存在
$\lambda > 0$ 时即使 $G^TG$ 奇异（如死时间估计错误）也有唯一解，增强鲁棒性
只取第一个元素 $\Delta u(t)$ 施加，即 $u(t) = u(t-1) + \mathbf{g}^T(\mathbf{w}-\mathbf{f})$

实际控制律（取 $\hat{\tilde{\mathbf{u}}}$ 第一行）：

$$u(t) = u(t-1) + \mathbf{g}^T(\mathbf{w} - \mathbf{f}) \tag{18}$$

其中 $\mathbf{g}^T$ 为 $(G^TG+\lambda I)^{-1}G^T$ 的第一行。控制律包含 $u(t-1)$，即对 $\Delta u$ 积分，自然具备积分作用。

可验证：由 Diophantine 方程(6)在 $q=1$ 处，$F_j(1)=1$（因 $\Delta A(1)=0$），故自由响应 $f(t+j)$ 的稳态值等于 $y(t)$，从而保证常值设定点下无稳态偏差。

3.1 控制域（Control Horizon）$N_U$

在 $N_U < N$ 之后假设控制增量为零：

$$\Delta u(t+j-1) = 0, \quad j > N_U \tag{19}$$

此假设等价于对后续控制变化施加无穷大权重。$\tilde{\mathbf{u}}$ 维度降为 $N_U$，矩阵 $G$ 变为 $N \times N_U$ 矩阵 $G_1$（右侧列截断），控制律更新为：

$$\hat{\tilde{\mathbf{u}}} = (G_1^T G_1 + \lambda I)^{-1} G_1^T (\mathbf{w} - \mathbf{f}) \tag{20}$$

求逆矩阵维度从 $N \times N$ 降为 $N_U \times N_U$。$N_U=1$ 时退化为标量运算，计算极为高效。

引入 $N_U$ 的三重作用

稳定非最小相位对象：$\lambda=0$ 时的最小方差律引入与非最小相位零点对应的控制信号增长模式；通过限制控制域可约束这些模式，即使 $\lambda=0$ 也能稳定闭环
降低计算量：矩阵求逆维度从 $N^2$ 降为 $N_U^2$，$N_U=1$ 时为标量
控制闭环极点：$N_U$ 至少等于不稳定或阻尼不足极点数目时可获得良好控制

3.2 输出域与控制域的选择

$N_1$（最小输出域）：若死时间 $k$ 精确已知，取 $N_1 = k$（$j < k$ 的预测不受 $u(t)$ 影响，计算多余）；死时间未知时取 $N_1=1$，并扩大 $B$ 的阶数覆盖可能的死时间范围。

$N_2$（最大输出域）：对含非最小相位性（初始响应为负）的对象，$N_2$ 需超过 $\deg B(q^{-1})$，使后续正向响应样本包含在代价中。实践中建议 $N_2$ 接近对象上升时间。

$N_U$（控制域）：最关键的设计参数：

对简单稳定对象（含死时间和非最小相位性），$N_U=1$ 通常给出可接受控制
对复杂系统，$N_U$ 至少等于不稳定或阻尼不足极点数目
增大 $N_U$ 使控制更激进；超过某阈值后继续增大几乎无改善
从零开始增大 $\lambda$ 可进一步阻尼控制动作

默认参数推荐

对工业过程控制：$N_1=1$，$N_2 \approx$ 对象上升时间，$N_U=1$，$\lambda$ 从0开始调整。大多数场合无需精细整定即可获得合理性能。

§4与其他方法的关系

GPC 融合五个核心思想：①CARIMA 模型；②超过死时间的有限域长程预测；③递推 Diophantine 方程；④控制增量加权；⑤控制域（域后增量置零）。以下各方法均为 GPC 的特例：

方法	GPC 退化条件	局限性
IDCOM (Richalet 1978)	无精确对应	仅限全零模型；不适用于开环不稳定/非最小相位对象
DMC (Cutler 1980)	阶跃响应模型 + 控制域	不适用开环不稳定；启发式处理偏差
GMV (Clarke 1975)	$N_1=N_2=k$，$N_U=1$	对变死时间敏感；只能稳定特定非最小相位对象
Ydstie (1984)	$N_1=N_2=d>k$，$N_U=1$，$\lambda=0$，CARMA	使用 CARMA 模型；阻尼不足极点时不稳定
EPSAC (De Keyser 1985)	$N_1=1$，$N_2=d$，$N_U=1$，$\lambda=0$，CARIMA	与默认参数 GPC 最接近
Peterka (1984)	哲学相近，但集中于 $N_2\to\infty$	无限域无法处理约束；矩阵分解推导较复杂

GPC 本质上是有限域方法，带受限控制域，因而可通过修改将约束（控制饱和、输出界限）纳入计算——这是无限域设计无法实现的推广。

§5仿真研究

在一个死时间、阶次和参数均发生大幅变化的被控对象上，对比固定 PID、自适应 GMV、自适应极点配置和自适应 GPC 四种方法（均使用标准递推最小二乘，遗忘因子0.9，无噪声）。

仿真被控对象按每 80 步切换一次（共400步）：

编号	采样区间	传递函数（连续时间）
1	1–79	$\dfrac{1}{1+10s+40s^2}$（二阶稳定）
2	80–159	$\dfrac{e^{-2s}}{1+10s+40s^2}$（含死时间）
3	160–239	$\dfrac{e^{-0.7s}}{1+10s}$（一阶+小死时间）
4	240–319	$\dfrac{1}{1+10s}$（一阶稳定）
5	320–400	$\dfrac{1}{10s(1+2.5s)}$（积分+一阶，开环不稳定）

结果汇总：

固定 PID：模型1和5尚可，其余因增益过大而性能不佳；模型3/4出现稳态偏差
自适应 GMV：模型2下弱阻尼振荡；整体优于固定PID
自适应极点配置：模型3/4（一阶）因 Diophantine 解奇异而完全失效；验证了对模型阶次变化敏感的问题
自适应 GPC（$N_1=1$，$N_2=10$，$N_U=1$，固定不变）：所有模型下均表现优异，每次换模型最多两步设定值跳变即可完全适应，无任何不稳定迹象

GPC 的关键优势

即使面对从稳定二阶到开环不稳定积分对象的剧烈切换，GPC 使用相同的默认参数均能稳定运行，而其他三种方法均在某些对象上失效。

§6结论

GPC 是一种新型鲁棒算法，适用于具有挑战性的自适应控制应用。推导和实现均简便；短控制域时可在微计算机上运行。

GPC 是 GMV 的进一步推广，可配备设计多项式和传递函数，继承 GMV 的理论框架。Part II 论文进一步探讨这些思想，给出 GPC 参数选择的稳定性理论依据，并展示更复杂控制任务上的仿真。

附录AGPC 最小化性质

考虑代价函数（含随机噪声）：

$$J_1 = E\left\{(\mathbf{y}-\mathbf{w})^T(\mathbf{y}-\mathbf{w}) + \lambda\tilde{\mathbf{u}}^T\tilde{\mathbf{u}}\right\}, \quad \mathbf{y} = G\tilde{\mathbf{u}} + \mathbf{f} + \mathbf{e}$$

其中 $\mathbf{e} = [E_1\xi(t+1), E_2\xi(t+2), \ldots, E_N\xi(t+N)]^T$。

A.1 开环反馈最优（OLFO）

假设未来控制序列与未来测量无关（开环施加），则 $E\{\tilde{\mathbf{u}}^T G^T \mathbf{e}\} = 0$，对 $\tilde{\mathbf{u}}$ 求导置零得：

$$\tilde{\mathbf{u}} = (G^TG + \lambda I)^{-1}G^T(\mathbf{w} - \mathbf{f})$$

这正是主文的控制律(17)。

A.2 闭环反馈最优（CLFO）等价性

对于 CARIMA（自回归型）扰动模型，可以证明：在因果性约束（$\Delta u(t+i)$ 只依赖 $t+i$ 时刻及之前的数据）下，CLFO 最优解与 OLFO 解等价。

证明思路：从最后一个控制量 $\Delta u(t+N-1)$ 反向递推，每步利用未来噪声的不可预见性（$E\{\Delta u(t+i)\xi(t+j)\}=0$，$j>i$），说明每步最优决策均等价于假设 $\xi(t+j)=0$（$j>i$）时的结果，最终第一步 $\Delta u(t)$ 的最优解与 OLFO 相同。

意义

GPC 的开环假设不是一种近似——对自回归型扰动模型，它恰好给出与闭环最优策略相同的第一步控制量。这为 GPC 的最优性提供了理论保证。

附录B一阶系统示例

考虑含非最小相位零点的一阶系统（含分数死时间）：

$$(1 + a_1 q^{-1})y(t) = (b_0 + b_1 q^{-1})u(t-1)$$

取数值 $a_1 = -0.9$，$b_0 = 1$，$b_1 = 2$，即：

$$(1 - 0.9q^{-1})y(t) = (1 + 2q^{-1})u(t-1)$$

注意 $B$ 的根为 $q^{-1}=-0.5$，即 $q=−2$（单位圆外零点），是非最小相位系统。

B.1 递推 Diophantine 计算示例

由 $\tilde{A} = \Delta A = (1-q^{-1})(1-0.9q^{-1}) = 1 - 1.9q^{-1} + 0.9q^{-2}$，递推结果：

$$E_1=1,\ f_0^{(1)}=1.9,\ f_1^{(1)}=-0.9$$

$$E_2=1+1.9q^{-1},\ f_0^{(2)}=2.71,\ f_1^{(2)}=-1.71$$

$$E_3=1+1.9q^{-1}+2.71q^{-2},\ f_0^{(3)}=3.439,\ f_1^{(3)}=-2.439$$

阶跃响应系数：$g_0=1,\ g_1=3.9,\ g_2=6.51$（与 $G_j = E_j B$ 一致）。

B.2 不同 $N_2$ 下的闭环极点（$N_U=1$，$\lambda=0$）

$N_2$	闭环极点	稳定性
1	$1-0.09q^{-1}$（对消非最小相位零点）	不稳定（零点对消律）
2（$> \deg B$）	$1-0.416q^{-1}$	稳定（极点在单位圆内）
3	$1-0.416q^{-1}$	稳定，极点趋向开环极点

关键结论

$N_2=1$（$\le \deg B$）时控制律对消非最小相位零点，导致闭环不稳定。$N_2 > \deg B$ 即可稳定。对一般一阶系统，可证明 $N_2 > \deg B(q^{-1}) + k - 1$ 即可稳定所有开环稳定对象（前提是直流增益符号正确）。

参考主要参考文献

Clarke, D.W., Mohtadi, C., & Tuffs, P.S. (1987). Generalized predictive control — Part I. The basic algorithm. Automatica, 23(2), 137–148.
Clarke, D.W. & Gawthrop, P.J. (1975). Self-tuning controller. Proc. IEE, 122, 929–934.
Cutler, C.R. & Ramaker, B.L. (1980). Dynamic Matrix Control. JACC, San Francisco.
Richalet, J. et al. (1978). Model predictive heuristic control. Automatica, 14, 413–428.
Peterka, V. (1984). Predictor-based self-tuning control. Automatica, 20, 39–50.
Ydstie, B.E. (1984). Extended horizon adaptive control. IFAC 9th World Congress.
De Keyser, R.M.C. & Van Cauwenberghe, A.R. (1985). Extended prediction self-adaptive control. IFAC Symp.