系统辨识

系统辨识（system identification）根据系统输入与输出的记录，建立或估计描述动态关系的模型。它包括选择模型结构、设计输入、估计参数和验证预测；拟合一组数值只是其中一步。不同结构可能在有限数据上给出相近拟合，却有不同的外推行为。

一阶动态模型中的两个参数

考虑采样模型 $$y_{k+1}=a{y}_k+b{u}_k+{\epsilon }_{k+1}.$$ $u_k$ 为第 $k$ 个时段的已知输入，$y_k$ 为输出，$a$ 描述无输入时的延续，$b$ 描述输入作用，$\epsilon$ 为未建模部分。采样周期固定；若物理单位不同，$b$ 的单位应为输出单位除以输入单位。

收集 $N$ 个转移，令设计矩阵第 $k$ 行为 $(y_k,u_k)$，目标向量为 $Y=(y_1,\dots ,y_N{)}^{𝖳}$，参数 $\theta =(a,b{)}^{𝖳}$。最小二乘估计最小化 $\Vert Y-X\theta {\Vert }^2$。展开平方并求导得 $$X^{𝖳}X\hat{\theta }=X^{𝖳}Y.$$ 若 $X$ 列满秩，唯一解为 $\hat{\theta }=(X^{𝖳}X{)}^{-1}{X}^{𝖳}Y$。令残差 $r=Y-X\hat{\theta }$，则 $X^{𝖳}r=0$；任意参数偏移 $h$ 使平方误差变为 $\Vert r{\Vert }^2+\Vert Xh{\Vert }^2$，证明该解确实最小，而不仅是驻点。

可复算的输入输出记录

构造无误差模型 $a=1/2,b=2,y_0=0$，输入序列为 $(1,0,1,0)$，逐步输出为 $(0,2,1,5/2,5/4)$。因此 $$X=\left(\begin{array}{cc}0& 1\\ 2& 0\\ 1& 1\\ 5/2& 0\end{array}\right),X^{𝖳}X=\left(\begin{array}{cc}45/4& 1\\ 1& 2\end{array}\right),X^{𝖳}Y=\left(\begin{array}{c}61/8\\ 9/2\end{array}\right).$$ 行列式为 $43/2$，解出 $\hat{a}=1/2,\hat{b}=2$，残差全部为零。改变输入的时段让数据包含不同方向的信息；若从零初态始终施加零输入，输出恒为零，任意 $a,b$ 都符合，完全无法辨识。

折线仅连接离散采样时刻。输入与下一时刻输出之间有一个动态转移，不能把同一横坐标的输入与输出当成静态比例关系。

参数计算与验证范围

下列代码直接累计两列正规方程，适合复算这个小例子；大规模或病态数据通常使用QR/SVD，避免正规方程放大条件数。输入必须有限，输出长度比输入多一，不改写数据；时间 $O(N)$、额外空间 $O(N)$ 用于输入副本，累计求解本身只需常数空间。

Python 3

import math

def fit_arx1(outputs, inputs):
    y, u = tuple(outputs), tuple(inputs)
    if len(u) < 2 or len(y) != len(u) + 1:
        raise ValueError("need at least two transitions")
    if not all(math.isfinite(x) for x in y + u):
        raise ValueError("finite observations required")
    yy = yu = uu = y_next = u_next = 0.0
    for k, uk in enumerate(u):
        yy += y[k] ** 2
        yu += y[k] * uk
        uu += uk ** 2
        y_next += y[k] * y[k + 1]
        u_next += uk * y[k + 1]
    det = yy * uu - yu ** 2
    if det <= 1e-12 * max(yy * uu, 1.0):
        raise ValueError("rank deficient or poorly scaled design")
    a = (uu * y_next - yu * u_next) / det
    b = (yy * u_next - yu * y_next) / det
    return a, b

print(fit_arx1([0, 2, 1, 2.5, 1.25], [1, 0, 1, 0]))

它返回 $(0.5,2.0)$。变式。 改成两次记录 $y=(0,3,3/2)$、$u=(1,0)$，可独立求出 $a=1/2,b=3$。

有噪声时，最小平方残差不自动保证参数代表真实机制。回归量可能与误差相关，闭环控制也可能使输入和扰动相关；一致性需要说明创新噪声、信息时序和激励等条件。验证应使用未参与估计的时段，分别检查一步预测、自由运行和残差结构。模型阶数、时间延迟和采样周期选错，不能仅靠增加拟合精度补救。持续激励进一步说明输入中需要哪些独立信息。