§2  CARIMA 与丢番图递推

CARIMA 模型 · 丢番图恒等式 · F/G 递推公式 · na=1 具体形式 · 自由响应
严格对照教材《第五章 单变量广义预测控制》§5.2~§5.3 · 代码 gpc_na1.py
1. CARIMA 模型（教材式5.1~5.2）
★ 核心

GPC 基于 CARIMA（Controlled AutoRegressive Integrated Moving Average）模型（教材式5.1）：

$$ A(q^{-1})\, y(t) = B(q^{-1})\, u(t-k) + \frac{C(q^{-1})}{\Delta}\, \xi(t) \tag{5.1}$$

其中 $\Delta = 1 - q^{-1}$ 是差分算子，$C(q^{-1})=1$（简化）。两边乘 $\Delta$，得等价形式（教材式5.2）：

$$ A(q^{-1})\,\Delta y(t) = B(q^{-1})\,\Delta u(t-k) + \xi(t) \tag{5.2}$$
  • $A(q^{-1}) = 1 + a_1 q^{-1} + \cdots + a_{n_a} q^{-n_a}$，阶次 $n_a$
  • $B(q^{-1}) = b_0 + b_1 q^{-1} + \cdots + b_{n_b} q^{-n_b}$，阶次 $n_b$
  • $k$：纯时延步数
为什么引入 CARIMA
分母含 $1/\Delta$，即 $z=1$ 处有极点（积分器）。这使得 GPC 对阶跃扰动有零稳态误差，不需要额外加积分项。

水箱参数：$n_a=1,\ n_b=0,\ k=24$，$A=[1,\ a_1]$，$B=[b_0]$。

2. 丢番图恒等式（教材式5.3~5.4）
★ 核心

对每个预测步 $j = 1, 2, \ldots, p$，求多项式对 $(E_j,\, F_j)$ 满足（教材式5.3）：

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

各多项式的阶次（教材式5.4）：

多项式阶次项数
$E_j(q^{-1})$$j-1$$j$ 项
$F_j(q^{-1})$$n_a$$n_a+1$ 项

由此定义（教材式5.5）：

$$ G_j(q^{-1}) = E_j(q^{-1})\, B(q^{-1}) \tag{5.5}$$

$G_j$ 的阶次 $= (j-1) + n_b = n_b + j - 1$，长度 $= n_b + j$（$n_b=0$ 时长度 $= j$）。

关键性质（教材式5.16）
$G_j(q^{-1}) = g_0 + g_1 q^{-1} + \cdots + g_{j-1} q^{-(j-1)} + \cdots$，其中 $g_0, g_1, \ldots, g_{j-1}$ 是系统(5.1)的阶跃响应系数，与 DMC 的脉冲响应系数一致。
3. 递推公式（教材§5.3，式5.9~5.19）
★ 核心

起点 $j=1$（教材式5.9~5.11）：

$$ E_1 = 1, \quad G_1 = B(q^{-1}) \tag{5.11}$$ $$ F_1(q^{-1}) = (1-a_1) + (a_1-a_2)q^{-1} + \cdots + (a_{n_a-1}-a_{n_a})q^{-(n_a-1)} + a_{n_a} q^{-n_a} \tag{5.10}$$

即 $f_1^i = a_i - a_{i+1}$（$i=0,\ldots,n_a-1$，$a_0=1$），$f_1^{n_a} = a_{n_a}$。

递推 $j-1 \to j$（令 $e_{j-1} = f_{j-1}^0$，教材式5.17~5.19）：

$$ e_{j-1} = f_{j-1}^0 \tag{5.17}$$ $$ f_j^i = f_{j-1}^{i+1} - e_{j-1}(a_{i+1} - a_i), \quad i = 0, 1, \ldots, n_a-1 \quad (a_0 = 1) \tag{5.18a}$$ $$ f_j^{n_a} = e_{j-1} \cdot a_{n_a} \tag{5.18b}$$ $$ g_j^i = g_{j-1}^i, \quad 0 \leq i \leq j-2 \quad \text{（低段不变）} \tag{5.19a}$$ $$ g_j^{n_b+j-1} = e_{j-1} \cdot b_{n_b} \quad \text{（最高项，}n_b=0\text{ 时即 }e_{j-1} \cdot b_0\text{）} \tag{5.19c}$$
代码对应（gpc_na1.py）
def diophantine(A, b0, p):
    a = A   # a[0]=1, a[1]=a_1
    # j=1 初值（式5.10）
    F[i] = a[i] - a[i+1]   # i=0..na-1
    F[na] = a[na]
    # 递推（式5.18a, 5.18b）
    F_new[i] = F_prev[i+1] - e*(a[i+1]-a[i])
    F_new[na] = e * a[na]
    # G_j（式5.19，nb=0）
    G_new[:j-1] = G_prev   # 低段继承
    G_new[j-1]  = e * b0   # 最高项
4. na=1 时的具体数值
○ 数值示例

$A = [1,\ a_1] = [1,\ -0.9846]$，$b_0 = 0.024361$，$n_a=1,\ n_b=0$：

$j=1$ 初值（教材式5.10，$n_a=1$）：

$$ F_1 = [(1-a_1),\ a_1] = [1.9846,\ -0.9846] \tag{2.1}$$ $$ G_1 = [b_0] = [0.024361] \tag{2.2}$$

$j=2$（$e_1 = f_1^0 = 1.9846$）：

$$ f_2^0 = f_1^1 - e_1(a_1 - a_0) = -0.9846 - 1.9846 \times (-0.9846 - 1) = 2.9540 \tag{2.3}$$ $$ f_2^1 = e_1 \cdot a_1 = 1.9846 \times (-0.9846) = -1.9540 \tag{2.4}$$ $$ G_2 = [0.024361,\ e_1 \cdot b_0] = [0.024361,\ 0.048347] \tag{2.5}$$
代码验证
Python 中 F_list[0] = [1.9846, -0.9846]，G_list[1] = [0.024361, 0.048347]，与手算一致。
$G_j$ 的前 $j$ 项 = 系统(5.1)的阶跃响应系数（教材式5.16 的关键性质）。
5. 自由响应 Y1（教材式5.21~5.22）
★ 核心

$j$ 步最优预测（教材式5.8）：

$$ \hat y_{t+j|t} = F_j(q^{-1})\, y(t) + G_j(q^{-1})\, \Delta u(t+j-k) \tag{5.8}$$

将 $G_j \Delta u$ 分为已知（历史）和未知（未来）两部分，得预测分解（教材式5.21）：

$$ \hat{\mathbf{Y}} = \mathbf{Y}^1 + \mathbf{G}\,\Delta\mathbf{U} \tag{5.21}$$

自由响应 $\mathbf{Y}^1$（已知部分）：

$$ Y^1_j = F_j(q^{-1})\, y(t) + \sum_{m: j-k-m<0} G_j[m]\, \Delta u(t+j-k-m) \tag{2.6}$$

即 $Y^1_j = F_j \cdot y_t$（$F_j$ 项）加上 $G_j$ 中对应历史 $\Delta u$ 的贡献。

$n_a=1$ 时 $F_j$ 有 2 项：$Y^1_j \supseteq f_j^0 y(t) + f_j^1 y(t-1)$。

$\mathbf{G}$ 矩阵（教材式5.22）：$p \times p$ 下三角，$G[i,j] = g_{i-j}$（$i \geq j$，否则0），$g_i = G_{i+1}[i]$。

代码对应（gpc_na1.py）
def free_response(F_list, G_list, y_buf, du_buf, p, k):
    for j in range(p):
        # F_j * y_t（式5.18，F_j有2项）
        fy = f_j^0*y(t) + f_j^1*y(t-1)
        # G_j 中历史 Delta*u 的贡献
        # lag = k + m - j，lag>0 时是历史值
        gu = sum(G_j[m]*du_buf[-lag] for m if lag>0)
        Y1[j] = fy + gu
6. 踩坑：F_j 维度错误（已修复）
⚠ 实际 bug

错误现象：仿真中 $y$ 完全不跟踪 SP，控制量几乎不动。

根本原因：第一版代码把 $\tilde A = A\Delta$ 当成新的 $A$ 来递推，$F_j$ 有3项（$n_a+2$），而教材式(5.4b) 明确 $F_j$ 阶次 $= n_a$，有 $n_a+1$ 项。

错误版本正确版本（教材）
递推对象$\tilde A = A\Delta$（2阶）原始 $A$（1阶）
$F_j$ 项数$n_a + 2 = 3$$n_a + 1 = 2$（式5.4b）
$F_1$ 初值$[\tilde A_0-\tilde A_1,\ \tilde A_1-\tilde A_2,\ \tilde A_2]$$[(1-a_1),\ a_1]$（式5.10）
修复方法（gpc_na1.py）
# 正确：用原始 A 的系数递推，F_j 有 na+1=2 项
F[i] = a[i] - a[i+1]   # i=0..na-1，a_0=1
F[na] = a[na]           # 式5.10
7. 踩坑：G_mat 索引偏移（C++ 版，已修复）
⚠ 实际 bug

错误现象：C++ 版 G2 从第 3 项开始与 Python 不一致，denom 差约 4 倍。

根本原因：Python 中 $G_j$ 长度为 $j$，最高项在索引 $j-1$；C++ 递推时多分配了一位，最高项写到了索引 $j$ 而不是 $j-1$，导致 $G\_\text{mat}[j][i]$ 的索引偏移了 1。

$j=2$ 时Python $G_\text{list}[1]$C++ 错误版
最高项位置索引 1（$G[1]=0.04835$）索引 2（错位）
修复方法（gpc_na1.cpp）
// 正确：G_new 长度 j+1，最高项在索引 j（对应 Python G_j[j-1]）
double G_new[P] = {0};
for (int kk = 0; kk < jj; kk++) G_new[kk] = G_cur[kk];
G_new[jj] = e * b0;   // 最高项在 jj

// G_mat[j][i] = G_new[j-i]
G_mat[jj][ii] = (idx >= 0) ? G_new[idx] : 0.0;