📚 课程目录第3章最小二乘法基础

3. 最小二乘法基础

各位同学，今天我们来聊聊系统辨识里最基础、也最实用的一个工具——最小二乘法。说实话，我入行做控制那会儿，第一个接触的辨识方法就是它。为什么？因为它简单、直观，而且效果往往还不错。你想想看，一堆带噪声的数据摆在你面前，怎么找到那条最靠谱的线？最小二乘法就是干这个的。

3.1 最小二乘原理：直觉与数学

最小二乘法的核心思想，说白了就是四个字：误差平方和最小。

我习惯这样理解：假设我们有一个系统，输入是 \( u \)，输出是 \( y \)。我们想用一个模型 \( \hat{y} = \theta^T \phi \) 来近似它。这里的 \( \theta \) 是参数，\( \phi \) 是输入或过去输出的组合。那么，模型预测值和真实值之间肯定有误差：

\[ e(k) = y(k) - \hat{y}(k) \]

最小二乘法要做的，就是找到一组参数 \( \theta \)，让所有时间点上的误差平方和最小：

\[ J(\theta) = \sum_{k=1}^{N} e(k)^2 = \sum_{k=1}^{N} (y(k) - \theta^T \phi(k))^2 \]

嗯，这里要注意：为什么是平方？而不是值？我在项目中遇到过这个问题。值函数在零点不可导，求极值很麻烦。平方函数光滑可导，求导后直接得到解析解。这就是数学上的便利性。

  核心要点： 最小二乘法不要求数据分布满足任何特定条件，它只是一个优化准则。但如果你假设噪声是高斯白噪声，那么最小二乘估计就是最优的（最大似然估计）。

3.2 批量最小二乘（BLS）公式推导

好，现在我们来做公式推导。这部分我建议你亲手推一遍，别只看。

假设我们有 N 组观测数据。把所有的输出写成向量形式：

\[ Y = \begin{bmatrix} y(1) \\ y(2) \\ \vdots \\ y(N) \end{bmatrix} \]

把所有回归向量写成矩阵形式（这叫信息矩阵或设计矩阵）：

\[ \Phi = \begin{bmatrix} \phi^T(1) \\ \phi^T(2) \\ \vdots \\ \phi^T(N) \end{bmatrix} \]

那么误差向量就是：

\[ E = Y - \Phi \theta \]

代价函数写成矩阵形式：

\[ J(\theta) = E^T E = (Y - \Phi \theta)^T (Y - \Phi \theta) \]

展开它：

\[ J(\theta) = Y^T Y - 2 \theta^T \Phi^T Y + \theta^T \Phi^T \Phi \theta \]

对 \( \theta \) 求导，令导数为零：

\[ \frac{\partial J}{\partial \theta} = -2 \Phi^T Y + 2 \Phi^T \Phi \theta = 0 \]

得到正规方程：

\[ \Phi^T \Phi \theta = \Phi^T Y \]

如果 \( \Phi^T \Phi \) 可逆，那么：

\[ \boxed{\hat{\theta}_{BLS} = (\Phi^T \Phi)^{-1} \Phi^T Y} \]

这就是批量最小二乘的估计公式。

我的经验： 在实际工程中，\( \Phi^T \Phi \) 可能接近奇异（比如数据量太少或输入信号不够丰富）。我曾经遇到过这种情况，算出来的参数跳来跳去，完全没法用。后来加了激励信号（比如PRBS），问题就解决了。所以，输入信号一定要充分激励系统。

3.3 在ARX模型中的应用

ARX模型是啥？AutoRegressive with eXogenous input，自回归外生输入模型。它的形式是：

\[ y(k) + a_1 y(k-1) + \dots + a_{n_a} y(k-n_a) = b_1 u(k-1) + \dots + b_{n_b} u(k-n_b) + e(k) \]

写成预测形式：

\[ \hat{y}(k) = -a_1 y(k-1) - \dots - a_{n_a} y(k-n_a) + b_1 u(k-1) + \dots + b_{n_b} u(k-n_b) \]

你看，这正好可以写成 \( \hat{y}(k) = \theta^T \phi(k) \) 的形式。其中：

\[ \theta = [a_1, a_2, \dots, a_{n_a}, b_1, b_2, \dots, b_{n_b}]^T \] \[ \phi(k) = [-y(k-1), -y(k-2), \dots, -y(k-n_a), u(k-1), u(k-2), \dots, u(k-n_b)]^T \]

然后直接用BLS公式就能算出参数。就这么简单。

注意事项：

ARX模型假设噪声是白噪声。如果噪声是有色的，最小二乘估计会有偏。这时候要用辅助变量法或广义最小二乘法。
模型阶次 \( n_a, n_b \) 的选择很关键。阶次太低，模型欠拟合；阶次太高，过拟合。我一般用AIC或BIC准则来选阶。
数据预处理很重要：去均值、去趋势、滤波。别小看这一步，有时候数据里有个直流偏置，辨识出来的参数就全偏了。

3.4 知识体系图

下面我用一张SVG图来总结本章的核心逻辑，帮你理清思路：

3.5 代码示例：用Python实现BLS辨识ARX模型

光说不练假把式。下面我给一个简单的Python实现，帮你把公式落地：

import numpy as np
import matplotlib.pyplot as plt

# 生成仿真数据：一个二阶ARX系统
# y(k) = 1.5*y(k-1) - 0.7*y(k-2) + 0.5*u(k-1) + 0.3*u(k-2) + e(k)
true_theta = np.array([1.5, -0.7, 0.5, 0.3])
na, nb = 2, 2
N = 500

# 输入：PRBS信号
u = np.random.choice([-1, 1], size=N)
# 噪声
e = 0.1 * np.random.randn(N)
# 输出
y = np.zeros(N)
for k in range(2, N):
    y[k] = 1.5*y[k-1] - 0.7*y[k-2] + 0.5*u[k-1] + 0.3*u[k-2] + e[k]

# 构造回归矩阵 Phi 和输出向量 Y
# 注意：从 max(na, nb) 开始取数据
start = max(na, nb)
Y = y[start:].reshape(-1, 1)

Phi = []
for k in range(start, N):
    row = [y[k-1], y[k-2], u[k-1], u[k-2]]
    Phi.append(row)
Phi = np.array(Phi)

# 批量最小二乘
theta_hat = np.linalg.inv(Phi.T @ Phi) @ Phi.T @ Y
theta_hat = theta_hat.flatten()

print("真实参数:", true_theta)
print("估计参数:", np.round(theta_hat, 4))

# 验证：用估计参数做一步预测
y_pred = Phi @ theta_hat
mse = np.mean((Y.flatten() - y_pred)**2)
print("一步预测MSE:", round(mse, 6))

我的建议： 跑完这段代码后，试试改变噪声强度 \( e \) 的方差，看看参数估计的精度变化。再试试减少数据量 N，看看什么时候 \( \Phi^T \Phi \) 开始接近奇异。这些实验比单纯看公式更能帮你理解最小二乘法的本质。

好了，这一章的内容就到这里。最小二乘法是系统辨识的基石，后面的递推最小二乘、带遗忘因子的最小二乘、辅助变量法等等，都是在这个基础上发展起来的。把这一章吃透，后面的路就好走了。