3. 传递函数模型阶次选择与参数估计

问题设定与数据生成

本节使用同一组离散输入输出数据，对传递函数模型 G(z) = B(z⁻¹)/A(z⁻¹) 分别拟合 1～4 阶，比较各阶模型的拟合优度、AIC 值与残差白噪声检验结果，最终选出合适阶次并给出参数估计。为避免与一阶、二阶案例重复，本节重点不在单点参数反算，而在多阶候选模型的横向比较与阶次判决。

构造一个已知真值的二阶系统作为教学示例：采样周期 Ts = 0.1 s，连续域 G(s) = 2/(s² + 1.2s + 1)，离散化后叠加标准差 0.02 的测量噪声，共 400 个采样点。真值仅用于最后对照，辨识过程只使用输入输出数据。

import numpy as np

np.random.seed(7)
Ts = 0.1
N = 400
# 真值二阶离散系统（零阶保持离散化）
from scipy.signal import cont2discrete
num = [2.0]
den = [1.0, 1.2, 1.0]
sys_d = cont2discrete((num, den), Ts, method='zoh')
b_true = sys_d[0][0][1:]   # 去掉首项 0
a_true = sys_d[1][0][1:]   # 去掉首项 1

u = np.random.randn(N)
y = np.zeros(N)
for k in range(2, N):
    y[k] = (-a_true[0]*y[k-1] - a_true[1]*y[k-2]
            + b_true[0]*u[k-1] + b_true[1]*u[k-2])
y_meas = y + 0.02*np.random.randn(N)
print("真值 a:", a_true, " b:", b_true)

多阶模型拟合与 AIC 计算

对每个候选阶次 n = 1,2,3,4，构造回归矩阵 Φ，用最小二乘求解参数向量 θ = [a₁…aₙ, b₁…bₙ]ᵀ，再计算残差平方和 V 与 AIC：

AIC = N·ln(V/N) + 2·p，其中 p = 2n 为待估参数个数。AIC 越小越好，惩罚项 2p 抑制过拟合。

def fit_arx(y, u, n):
    N = len(y)
    Phi = np.zeros((N-n, 2*n))
    Y = y[n:]
    for i in range(n):
        Phi[:, i]     = -y[n-1-i:N-1-i]   # a 项
        Phi[:, n+i]   =  u[n-1-i:N-1-i]   # b 项
    theta, *_ = np.linalg.lstsq(Phi, Y, rcond=None)
    resid = Y - Phi @ theta
    V = resid @ resid
    p = 2*n
    aic = N*np.log(V/N) + 2*p
    return theta, V, aic, resid

results = {}
for n in [1, 2, 3, 4]:
    theta, V, aic, resid = fit_arx(y_meas, u, n)
    results[n] = dict(theta=theta, V=V, aic=aic, resid=resid)
    print(f"n={n}  V={V:.4f}  AIC={aic:.2f}")

典型输出（示例数据下的预期量级）：n=1 时 AIC 明显偏高；n=2 时 AIC 最低；n=3、4 时 V 略降但惩罚项增大，AIC 回升。这正是“拟合优度 vs 复杂度”的权衡点。

残差白噪声检验

AIC 只反映拟合与复杂度的折中，还需确认残差是否近似白噪声。对残差做自相关，检查是否落在 ±1.96/√N 置信带内；若某阶残差仍有显著自相关，说明模型未充分提取动态信息。

def residual_acf(resid, max_lag=20):
    r = resid - resid.mean()
    acf = np.array([np.sum(r[:len(r)-k]*r[k:])/np.sum(r*r)
                    for k in range(max_lag+1)])
    bound = 1.96/np.sqrt(len(r))
    return acf, bound

for n in [1, 2, 3, 4]:
    acf, bound = residual_acf(results[n]['resid'])
    exceed = np.sum(np.abs(acf[1:]) > bound)
    print(f"n={n}  超界滞后数={exceed}  置信带=±{bound:.3f}")

判读逻辑：n=1 残差自相关通常显著超界（模型欠拟合）；n=2 残差基本落入置信带；n≥3 虽也落入带内，但 AIC 已回升，且参数方差增大。综合 AIC 最小 + 残差白噪声通过，选 n=2。

最终模型与验证

选定 n=2 后，用全部数据重新估计参数，并做一步预测交叉验证：将后 100 个点作为验证集，用前 300 点估计的模型预测验证集输出，计算拟合度 fit = 100·(1 − ‖y−ŷ‖/‖y−ȳ‖)。若拟合度高于 85% 且残差 ACF 无超界，即可认为模型可用。

# 用前 300 点估计，后 100 点验证
theta2, V2, aic2, resid2 = fit_arx(y_meas[:300], u[:300], 2)
a1, a2, b1, b2 = theta2

y_pred = np.zeros(100)
for k in range(2, 100):
    idx = 300 + k
    y_pred[k] = (-a1*y_meas[idx-1] - a2*y_meas[idx-2]
                 + b1*u[idx-1] + b2*u[idx-2])

y_val = y_meas[300:]
fit = 100*(1 - np.linalg.norm(y_val[2:]-y_pred[2:])
             / np.linalg.norm(y_val[2:]-y_val[2:].mean()))
print(f"验证集拟合度 = {fit:.1f}%")
print(f"最终参数 a1={a1:.4f}, a2={a2:.4f}, b1={b1:.4f}, b2={b2:.4f}")