§7 代码对照

gpc_na1.py 每个函数 ↔ 教材公式 ↔ 笔记章节

代码文件：mpc_split/gpc_na1.py · 教材：第五章 §5.1~§5.5

0. 文件结构总览

代码函数/类	教材对应	笔记章节	作用
`norm_y / norm_u`	—	§1-1	物理量 ↔ 归一化 [0,1]
`diophantine(A, b0, p)`	式5.10~5.19	§2-3	递推求 $F_j, G_j$
`free_response(...)`	式5.21	§2-5	计算自由响应 $Y^1$
`soften_setpoint(...)`	式5.20	§3-2	柔化设定值轨迹 $W$
`stepped_gpc(...)`	式5.33~5.36	§3-4	阶梯式 GPC 控制律
`class RLS`	§5.1	§4-2	递推最小二乘在线辨识
主循环 `for t in range(N_sim)`	§5.1 图1	§4-4	自校正 GPC 每步执行

1. 归一化 ↔ §1 系统辨识

所有计算在归一化偏差空间进行，物理量仅在输入/输出端转换。

教材公式

$y_\text{norm} = (y_V - 1.0) / 3.634$
$u_\text{norm} = (u_\text{mA} - 4.0) / 16.0$
偏差量：$y_\text{dev} = y_\text{norm} - y_0$

代码（gpc_na1.py）

Y_MIN, Y_RANGE = 1.0, 3.634
U_MIN, U_RANGE = 4.0, 16.0

def norm_y(y):
    return (y - Y_MIN) / Y_RANGE

# 偏差量（主循环中）
y_dev = y_norm - y0n
        

为什么必须在归一化空间辨识

原始空间 $b_0^\text{raw}=0.0055$，稳态增益 $K_{ss}=0.36$，到达 3V 需要 28.6 mA（超限）。归一化后 $b_0=0.02436$，$K_{ss}=1.58$，9.6 mA 即可到达 3V。

2. 丢番图递推 ↔ §2 CARIMA 与丢番图

教材公式（§5.3）

起点 $j=1$（式5.10）
$f_1^i = a_i - a_{i+1}$，$i=0..n_a-1$，$a_0=1$
$f_1^{n_a} = a_{n_a}$
$G_1 = [b_0]$（式5.11）

递推 $j-1\to j$（式5.17~5.19）
$e_{j-1} = f_{j-1}^0$
$f_j^i = f_{j-1}^{i+1} - e_{j-1}(a_{i+1}-a_i)$
$f_j^{n_a} = e_{j-1} \cdot a_{n_a}$
$G_j[0..j-2] = G_{j-1}[0..j-2]$
$G_j[j-1] = e_{j-1} \cdot b_0$

代码（diophantine 函数）

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]
    G = [b0]  # 式5.11

    for j in range(2, p+1):
        e = F_prev[0]  # 式5.17
        # 式5.18a
        F_new[i] = F_prev[i+1] - e*(a[i+1]-a[i])
        F_new[na] = e * a[na]  # 式5.18b
        # 式5.19
        G_new[:j-1] = G_prev
        G_new[j-1]  = e * b0
        

na=1 时的具体值

$F_1 = [1.9846,\ -0.9846]$（2项），$G_1 = [0.02436]$，$G_2 = [0.02436,\ 0.04835]$
$G_j$ 前 $j$ 项 = 系统阶跃响应系数（教材式5.16 关键性质）

3. 自由响应 ↔ §2-5 自由响应 Y1

教材公式（式5.8, 5.21）

$\hat y_{t+j|t} = F_j(q^{-1})y_t + G_j(q^{-1})\Delta u_{t+j-k}$

自由响应（已知部分）：
$Y^1_j = F_j \cdot y_t + \sum_{m:\,j-k-m<0} G_j[m]\,\Delta u_{t+j-k-m}$

$n_a=1$：$F_j$ 有2项
$F_j \cdot y_t = f_j^0 y(t) + f_j^1 y(t-1)$

历史 $\Delta u$ 贡献：
$\text{lag} = k + m - j$，$\text{lag}>0$ 时是历史值

代码（free_response 函数）

def free_response(F_list, G_list,
                    y_buf, du_buf, p, k):
    for j in range(p):
        Fj = F_list[j]  # 长度2
        Gj = G_list[j]  # 长度j

        # F_j * y_t（式5.8）
        fy = Fj[0]*y(t) + Fj[1]*y(t-1)

        # 历史 Δu 贡献
        for m in range(len(Gj)):
            lag = k + m - j
            if lag > 0:
                gu += Gj[m] * du_buf[-lag]

        Y1[j] = fy + gu
        

4. 柔化设定值 ↔ §3-2 柔化设定值轨迹

教材公式（式5.20）

$w_{t+k-1} \approx y(t)$（用当前值近似）
$w_{t+k+i} = \alpha\, w_{t+k+i-1} + (1-\alpha)\, SP$
$i = 0, 1, \ldots, p-1$

$\alpha \in [0,1)$：柔化因子
$\alpha=0$：$W = SP$（无柔化）
$\alpha=0.9$：约30步逼近 SP

代码（soften_setpoint 函数）

def soften_setpoint(y_now, SP, alpha, p):
    W = zeros(p)
    # i=0（式5.20，用 y(t) 近似）
    W[0] = alpha*y_now + (1-alpha)*SP
    # i=1..p-1
    for i in range(1, p):
        W[i] = alpha*W[i-1] + (1-alpha)*SP
    return W
        

5. 阶梯式 GPC ↔ §3-4 阶梯式 GPC

教材公式（式5.28, 5.33~5.36）

$G_1 \in \mathbb{R}^{p \times p_u}$：$G_1[j,i] = g_{j-i}$（式5.28）

阶梯约束（式5.33）：$\Delta u_{t+i} = \beta^i \delta$
$\mathbf{G}_2 = G_1 \cdot [1,\beta,\ldots,\beta^{p_u-1}]^T$（式5.34）

控制律（式5.36）：
$\delta = \dfrac{\mathbf{G}_2^T(\mathbf{W}-\mathbf{Y}^1)}{\mathbf{G}_2^T\mathbf{G}_2 + \lambda\sum_{i=0}^{p_u-1}\beta^{2i}}$

代码（stepped_gpc 函数）

def stepped_gpc(G_list, Y1, W, p, pu, lam, beta):
    # G_1 矩阵（式5.28）
    for j in range(p):
        G_mat[j,i] = G_list[j][j-i]

    # 阶梯向量（式5.33）
    beta_vec = [beta**i for i in range(pu)]

    # G_2（式5.34）
    G2 = G_mat @ beta_vec

    # 分母（式5.35）
    denom = G2@G2 + lam*sum(beta_vec**2)

    # 控制律（式5.36）
    return G2 @ (W - Y1) / denom
        

6. RLS 更新 ↔ §4-2 RLS 算法

教材公式（§5.1）

$\varphi(t) = [y(t-1),\ u(t-k)]^T$
$\mathbf{K} = \mathbf{P}\varphi / (\lambda + \varphi^T\mathbf{P}\varphi)$
$\hat\theta \leftarrow \hat\theta + \mathbf{K}(y - \varphi^T\hat\theta)$
$\mathbf{P} \leftarrow (\mathbf{P} - \mathbf{K}\varphi^T\mathbf{P})/\lambda$

符号约定：
$\hat\theta[0] = -a_1 > 0$（正值）
$A = [1, -\hat\theta[0]]$（$a_1 < 0$）

代码（class RLS + 主循环）

# 主循环中构造 A
a1_rls = rls.theta[0]  # = -a_1 > 0
A_cur = [1.0, -a1_rls]  # a_1 < 0

# RLS 更新（控制律之后）
phi = [y_buf[-2], u_buf[-(k+1)]]
rls.update(phi, y_t)

# class RLS.update
Pp = P @ phi
K  = Pp / (lam + phi @ Pp)
theta += K * (y - phi @ theta)
P = (P - outer(K, phi@P)) / lam
        

7. 主循环逐行注解

○ 每步执行顺序（对应 §4-4 自校正 GPC 流程）

for t in range(N_sim):
    SP = SP_seq[t]

    # ① RLS 参数 → A → 丢番图（§2-3）
    a1_rls = rls.theta[0]          # θ_0 = -a_1 > 0
    A_cur  = [1.0, -a1_rls]        # A = [1, a_1]，a_1 < 0
    F_list, G_list = diophantine(A_cur, b0_rls, p)

    # ② 自由响应（§2-5）
    Y1 = free_response(F_list, G_list, y_buf, du_buf, p, k)

    # ③ 柔化设定值（§3-2）
    W  = soften_setpoint(y_buf[-1], SP, alpha, p)

    # ④ 阶梯式 GPC 控制律（§3-4）
    du_t = stepped_gpc(G_list, Y1, W, p, pu, lam, beta)
    du_t = clip(du_t, -DU_MAX, DU_MAX)   # 增量限幅
    u_t  = clip(u_prev + du_t, 0, 1)     # 绝对值限幅
    du_t = u_t - u_prev                  # 实际执行增量

    # ⑤ 更新历史缓冲
    du_buf.append(du_t)   # Δu 历史，供自由响应使用
    u_buf.append(u_t)
    y_t = a1_rls*y_buf[-1] + b0_rls*u_buf[-(k+1)]  # ARX 仿真
    y_buf.append(y_t)

    # ⑥ RLS 更新（控制律之后，§4-2）
    phi = [y_buf[-2], u_buf[-(k+1)]]
    rls.update(phi, y_t)
    

顺序的关键

RLS 更新必须在控制律计算之后：先用旧参数算控制量并施加，再用新观测 $y(t)$ 更新参数，供下一步使用。如果先更新参数再算控制律，会用到还未发生的观测。