4. ARX模型结构确定与参数辨识

问题设定与数据准备

ARX（AutoRegressive with eXogenous input）模型把系统输出写成“历史输出 + 历史输入 + 白噪声”的线性组合：

A(q⁻¹)·y(k) = B(q⁻¹)·u(k−nk) + e(k)

其中 A(q⁻¹)=1+a₁q⁻¹+…+a_na·q⁻¹，B(q⁻¹)=b₁+b₂q⁻¹+…+b_nb·q⁻¹，nk 为纯延迟（以采样周期为单位）。本节任务：给定一段离散输入输出序列，先确定 na、nb、nk，再用最小二乘估计系数。

为便于复现，构造一个已知真值的教学示例系统（采样周期 0.1 s）：

y(k) = 1.5·y(k−1) − 0.7·y(k−2) + 0.8·u(k−3) + 0.4·u(k−4) + e(k)
真值：na = 2, nb = 2, nk = 3, a = [−1.5, 0.7], b = [0.8, 0.4]

输入 u(k) 取零均值、方差 1 的伪随机序列，e(k) 取方差 0.01 的高斯白噪声，共 400 个采样点。前 300 点用于辨识，后 100 点用于交叉验证。

候选结构搜索：损失函数与残差分析

对每一组候选 (na, nb, nk)，构造回归向量：

φ(k) = [−y(k−1), …, −y(k−na), u(k−nk), …, u(k−nk−nb+1)]ᵀ

最小二乘解为 θ̂ = (ΦᵀΦ)⁻¹ΦᵀY，损失函数取残差平方和 V = Σe²(k)。仅看 V 会随阶次单调下降，因此必须配合惩罚项与残差检验：

• 损失函数：V 随 na、nb 增大而下降，但 nk 取错时 V 会明显偏高，这是识别延迟最直接的信号。
• 惩罚准则：用 AIC = N·ln(V/N) + 2p（p = na + nb）抑制过拟合。
• 残差白噪声检验：计算残差自相关 R_ε(τ)，若在 τ ≥ 1 处显著超出 ±1.96/√N 置信带，说明结构欠拟合（常见于 na 偏小或 nk 偏小）。

关键实现：回归矩阵与最小二乘

下面给出可直接运行的 Python 实现，覆盖数据生成、候选结构遍历、AIC 计算、残差自相关检验与交叉验证。所有参数在本例中显式定义，不依赖前文代码。

import numpy as np

rng = np.random.default_rng(2024)
N = 400
u = rng.standard_normal(N)
e = 0.1 * rng.standard_normal(N)          # 方差 0.01

# 真值系统：na=2, nb=2, nk=3
y = np.zeros(N)
for k in range(N):
    y[k] = (1.5 * (y[k-1] if k-1 >= 0 else 0.0)
            - 0.7 * (y[k-2] if k-2 >= 0 else 0.0)
            + 0.8 * (u[k-3] if k-3 >= 0 else 0.0)
            + 0.4 * (u[k-4] if k-4 >= 0 else 0.0)
            + e[k])

Ntr = 300
def build_regressors(y, u, na, nb, nk, start, stop):
    rows, Y = [], []
    k0 = max(na, nk + nb - 1)
    for k in range(max(start, k0), stop):
        phi = [-y[k-i] for i in range(1, na+1)]
        phi += [u[k-nk-j] for j in range(nb)]
        rows.append(phi)
        Y.append(y[k])
    return np.array(rows), np.array(Y)

def fit_arx(na, nb, nk):
    Phi, Y = build_regressors(y, u, na, nb, nk, 0, Ntr)
    theta, *_ = np.linalg.lstsq(Phi, Y, rcond=None)
    resid = Y - Phi @ theta
    V = float(resid @ resid)
    p = na + nb
    aic = len(Y) * np.log(V / len(Y)) + 2 * p
    return theta, V, aic, resid

def residual_whiteness(resid, max_lag=10):
    r = resid - resid.mean()
    denom = r @ r
    acf = [float(r[:-t] @ r[t:] / denom) for t in range(1, max_lag+1)]
    bound = 1.96 / np.sqrt(len(r))
    return acf, bound

results = []
for na in range(1, 5):
    for nb in range(1, 4):
        for nk in range(0, 5):
            theta, V, aic, resid = fit_arx(na, nb, nk)
            acf, bound = residual_whiteness(resid)
            white = all(abs(rho) < bound for rho in acf)
            results.append((na, nb, nk, V, aic, white, theta))

# 先按残差白噪声过滤，再取 AIC 最小
candidates = [r for r in results if r[5]]
best = min(candidates, key=lambda r: r[4])
print("选定结构 na, nb, nk =", best[0], best[1], best[2])
print("AIC =", round(best[4], 2), " V =", round(best[3], 4))
print("theta =", np.round(best[6], 3))

运行后预期输出（教学示例，非实测）：

选定结构 na, nb, nk = 2 2 3
AIC = -1180.6  V = 2.87
theta = [-1.502  0.701  0.798  0.402]

系数与真值 [−1.5, 0.7, 0.8, 0.4] 接近，说明结构选择正确。若把 nk 误设为 2，V 会明显抬升（示例中约 3.9），且残差在 τ=1 处出现显著负相关，这是延迟错位的典型指纹。

结果判读与交叉验证

结构确定后，用留出的后 100 点做一步预测，计算拟合优度：

def fit_percent(y_true, y_pred):
    return 100 * (1 - np.linalg.norm(y_true - y_pred) / np.linalg.norm(y_true - y_true.mean()))

na, nb, nk = best[0], best[1], best[2]
theta = best[6]
Phi_va, Y_va = build_regressors(y, u, na, nb, nk, Ntr, N)
Y_hat = Phi_va @ theta
print("验证段 FIT =", round(fit_percent(Y_va, Y_hat), 2), "%")

判读要点：

• FIT ≥ 90%：结构合理，可用于后续控制器设计或仿真。
• FIT 偏低但训练段 V 很小：典型过拟合，应回退阶次或检查数据段是否含突变。
• 残差 ACF 在 τ=1 显著：优先怀疑 nk 偏小或 na 偏小；若在 τ 较大处仍有峰，考虑噪声为有色噪声，应转向 ARMAX/OE 结构。
• a₁ 接近 ±1：系统接近积分或临界稳定，需检查采样周期与数据长度是否足够。