最小二乘法在系统辨识中的实战应用与优化策略-优快云博客

1. 从“猜”到“算”：最小二乘法如何成为系统辨识的基石

我刚入行做控制系统的时候，最头疼的就是给一个黑盒子似的设备建数学模型。比如，给你一个温控箱，你知道输入电流能改变温度，但具体电流变化1安培，温度会怎么变、多久能变到位？这中间的数学关系就是“系统模型”。最开始，老师傅们靠经验“猜”参数，调试全靠手感，一个项目耗上几个月是常事。直到我真正用上了最小二乘法，才明白什么叫“数据驱动”的建模——让数据自己说话，把“猜”变成“算”。

最小二乘法的核心思想，其实特别生活化。想象一下，你在纸上随手画了一堆散点，想找一条最合适的直线穿过它们，用来预测新的点。你怎么判断这条线“最合适”？最小二乘法说：找那条让所有散点到这条直线垂直距离的平方和最小的线。为什么是平方和？因为距离有正有负，直接加可能会抵消，平方一下就能保证都是正数，而且数学上求导找最小值特别方便。这个“距离平方和最小”的原则，就是“二乘”的由来。

在系统辨识里，这个“散点”就是你的实验数据：输入信号和对应的输出信号。那个“黑盒子”系统内部的结构（比如，是惯性环节还是振荡环节？时间常数多大？）我们用一个带未知参数的方程来假设，比如一个差分方程。最小二乘法的任务，就是根据输入输出的数据，反推出那些未知参数最可能的值，使得模型计算出来的输出和实际测到的输出，两者的误差平方和最小。

我最早用它是在一个伺服电机的参数辨识项目上。电机模型里有几个关键参数，像转动惯量、阻尼系数，手册给的值和实际装配后总有偏差。我们给电机输入一系列已知的电流信号，同时用编码器记录它的位置变化，采集了几千组数据。然后，按照最小二乘的公式一套，用MATLAB一算，几分钟就把参数估计出来了。把参数代入模型，仿真出来的位置曲线和实际记录的重合度高达95%以上，那一刻的成就感，比盲调一周都强。这方法之所以成为基石，就是因为它原理直观、数学严谨、实现简单，为后续所有优化算法提供了一个可靠的起点。

2. 实战第一步：手把手完成你的第一次最小二乘辨识

光说原理太抽象，咱们直接看一个能跑通的例子。假设我们有一个简单的直流电机速度控制系统，我们认为它的一阶离散模型可以写成这样： y(k) = a * y(k-1) + b * u(k-1) 其中，y(k)是当前时刻的速度，u(k)是当前时刻的电压输入，a和b就是我们待辨识的两个参数。

我们的任务就是通过实验数据，把a和b找出来。怎么做呢？

第一步，设计实验并采集数据。 这是最关键的一步，数据质量直接决定辨识结果。你不能只输入恒定的电压，那样激不出系统的全部动态。我常用的方法是输入一个伪随机二进制序列（PRBS），它看起来像随机的0/1跳变，但具有近似白噪声的频谱，能很好地激励系统在各种频率下的响应。我们用这个信号作为电压u，同时以固定采样周期（比如10ms）记录电机的实际速度y。假设我们采集了N=1000组数据。

第二步，构造数据矩阵。 根据我们的模型，可以把每一组数据的关系写成： y(k) = [y(k-1), u(k-1)] * [a; b] 看，这就变成了我们熟悉的线性方程形式：z = H * θ。其中：

z 是我们的观测输出向量，z = [y(2), y(3), ..., y(1000)]^T
H 是数据矩阵，每一行是[y(k-1), u(k-1)]，所以H是一个999行2列的矩阵。
θ 就是待求参数向量 [a; b]

第三步，套用公式进行计算。 最小二乘的标准解是：θ_hat = inv(H^T * H) * H^T * z 这里inv()表示求逆，^T表示转置。在MATLAB或Python里实现起来非常简单。

import numpy as np
import matplotlib.pyplot as plt

# 假设我们已经有了数据数组 u 和 y
# u = [...], y = [...]， 长度均为1000

# 构造数据矩阵H和观测向量z
N = len(y)
H = np.column_stack((y[:-1], u[:-1]))  # 第k行是 [y(k-1), u(k-1)]
z = y[1:].reshape(-1, 1)  # 从y(2)到y(N)

# 最小二乘估计
theta_hat = np.linalg.inv(H.T @ H) @ H.T @ z
a_est, b_est = theta_hat.flatten()

print(f"估计的参数 a = {a_est:.4f}, b = {b_est:.4f}")

# 验证：用估计的参数和输入u，模拟计算输出y_sim
y_sim = np.zeros(N)
y_sim[0] = y[0]  # 初始值
for k in range(1, N):
    y_sim[k] = a_est * y_sim[k-1] + b_est * u[k-1]

# 绘制对比图
plt.figure(figsize=(10, 6))
plt.plot(y, 'b-', label='实际测量速度', linewidth=2)
plt.plot(y_sim, 'r--', label='模型仿真速度', linewidth=1.5)
plt.xlabel('采样点')
plt.ylabel('速度')
plt.legend()
plt.title('最小二乘法辨识结果验证')
plt.grid(True)
plt.show()

跑完这段代码，你就能看到模型输出和实际数据的拟合情况。如果曲线贴合得好，说明辨识成功；如果差别大，可能模型结构没选对（比如该用二阶模型），或者数据噪声太大。这就是一次完整的离线批处理最小二乘辨识。它把问题转化成了一个线性代数求解，不需要迭代，一次计算就得到全局最优解，对于新手来说非常友好，能快速建立信心。

3. 当数据“厚此薄彼”：加权最小二乘法的精妙之处

用上面基本方法做完第一个项目后，我很快遇到了新问题。在一次无人机飞控系统的传感器参数辨识中，我发现不同时间段的传感器数据质量差异很大。比如，在电机刚启动的瞬间，振动大，电流传感器读数噪声也大；而在平稳飞行阶段，数据就干净得多。如果对所有数据“一视同仁”，那些噪声大的“坏数据”就会把我们的参数估计“带偏”。

这就引出了加权最小二乘法的核心思想：给不同的测量数据赋予不同的“信任度”或权重。你更信任的数据，就在计算中给它更高的话语权；你不信任的数据，就降低它的影响力。这就像是一个专家评审团，资深专家的票权重要高于新人的票。

数学上，这体现在目标函数里。普通最小二乘是最小化误差平方和 J = Σ e(k)^2。加权最小二乘则把它变成了 J = Σ w(k) * e(k)^2。这里的 w(k) 就是第k个数据点的权重。如果某个数据点我们非常确信，就设一个很大的 w(k)（比如10）；如果某个数据点我们知道它噪声很大，就设一个很小的 w(k)（比如0.1）。

那么，权重矩阵 W 怎么选呢？最常用也最合理的选择是取误差协方差矩阵的逆。假设我们知道不同时刻测量噪声的方差，记作 σ^2(k)。如果这些噪声是相互独立的，那么权重矩阵 W 就是一个对角阵，对角线上的元素就是 w(k) = 1 / σ^2(k)。这意味着，噪声方差越大的数据点，其权重越小。这个特例被称为马尔可夫估计，理论证明它能在所有线性无偏估计中，给出方差最小的估计结果，也就是最优的。

在实际工程中，我们往往不知道精确的 σ^2(k)。我的经验是，可以通过以下方法动态确定权重：

残差分析法：先用普通最小二乘做一次初步辨识，计算每个数据点的残差（模型预测值与实际值之差）。残差大的点，通常意味着该点拟合差，可能受噪声或异常值影响大，在下一次迭代中就给一个较小的权重。
分段加权：根据先验知识对数据分段。比如，我知道系统启动阶段数据不可靠，就给前100个数据点统一设一个很小的权重（如0.01）；在稳态阶段，权重设为1。
自适应加权：在递推算法中（下一节会讲），可以根据新息（预测误差）的大小实时调整权重因子。

使用加权最小二乘后，上面那个无人机传感器辨识的结果得到了显著改善。模型在平稳飞行阶段的预测精度提高了，因为那些阶段的“好数据”得到了应有的重视。加权，本质上是一种利用先验知识对数据进行预处理的手段，它让算法变得更智能，也更贴近工程实际——毕竟，现实世界的数据从来都不是同等质量的。

4. 让模型“活”起来：递推最小二乘与在线实时辨识

批处理最小二乘法虽然强大，但它有个硬伤：必须等所有数据都采集完了才能计算。这在很多场景下是不可接受的。比如，我想让一个机械臂在运行过程中，实时地更新自身的动力学模型参数，以适应抓取不同负载的变化。这时候，就需要递推最小二乘法。

递推算法的精髓是“利用旧结果，结合新数据，快速得到新结果”。它不需要保存和处理全部历史数据，只需要记住上一时刻的参数估计值和一个特定的矩阵（通常记为 P 矩阵，代表估计的协方差或信息矩阵），当新的输入输出数据对到来时，通过一套固定的更新公式，就能立即算出新的参数估计。

公式看起来有点复杂，但理解其思想后很简单： 新的估计 = 旧的估计 + 增益矩阵 * (新观测值 - 基于旧估计的预测值) 后面括号里的 (新观测值 - 基于旧估计的预测值) 被称为“新息”，它代表了新数据带来的信息量。增益矩阵 K 则决定了我们有多信任这个新信息。P 矩阵在这个过程中不断更新，它反映了我们当前参数估计的不确定度。随着数据越来越多，P 矩阵通常会越来越小，意味着我们的估计越来越确定。

# 递推最小二乘法 (RLS) 的Python示例
def rls_identification(u, y, n_params, lambda_=1.0):
    """
    u: 输入信号序列
    y: 输出信号序列
    n_params: 待辨识参数个数 (例如，对于一阶系统为2)
    lambda_: 遗忘因子，通常取0.95~1.0
    """
    N = len(y)
    theta = np.zeros((n_params, 1))  # 参数初始估计，可以设为小随机数或零
    P = 1e6 * np.eye(n_params)  # 初始协方差矩阵，取一个很大的数，表示初始不确定性很大
    theta_history = np.zeros((N, n_params))  # 记录参数变化历史

    for k in range(n_params, N):  # 需要足够的数据构造初始回归向量
        # 构造当前时刻的回归向量 h(k)
        # 对于模型 y(k) = a1*y(k-1) + b1*u(k-1)， h(k) = [y(k-1), u(k-1)]
        h = np.array([[y[k-1], u[k-1]]]).T  # 列向量

        # 计算先验误差（新息）
        y_pred = h.T @ theta
        error = y[k] - y_pred

        # 计算增益矩阵 K
        K = (P @ h) / (lambda_ + h.T @ P @ h)

        # 更新参数估计
        theta = theta + K * error

        # 更新协方差矩阵 P
        P = (1/lambda_) * (P - K @ h.T @ P)

        theta_history[k, :] = theta.flatten()

    return theta, theta_history

这段代码实现了一个带遗忘因子的递推最小二乘。lambda_（λ）就是遗忘因子，它是应对“数据饱和”和“时变系统”的关键。当 λ=1 时，算法记得所有历史数据，适用于时不变系统。但当系统参数随时间缓慢变化时（比如电池内阻随老化增加），旧数据的信息会过时，反而会阻碍算法跟踪新参数。这时，我们设置 λ < 1（比如0.99），让算法逐渐“遗忘”旧数据，更关注新数据，从而实现对时变参数的跟踪。

我在一个自适应温度控制项目中就用了这个。加热炉的散热特性会随着环境温度和炉内物料变化。我们采用RLS算法，在线实时更新模型参数，控制器再根据新模型调整功率输出。整个过程完全自动，无需人工干预，系统始终能保持最佳控制性能。从“离线批处理”到“在线递推”，最小二乘法从一种建模工具，进化成了让系统具备“自学习”能力的核心引擎。

5. 直面噪声：增广最小二乘法与更精确的建模

前几种方法都有一个重要的前提假设：系统中的噪声是白噪声，并且只加在输出端。这个假设在很多时候是为了数学推导方便。但现实中，噪声可能存在于各个环节，并且常常是有色的（即前后时刻的噪声相关）。忽略噪声特性或者假设不当，会导致参数估计产生偏差，即使数据量无限增大，估计值也不会收敛到真实值。

为了同时估计系统参数和噪声模型的参数，增广最小二乘法 被提了出来。它的思路非常巧妙：把噪声模型的参数也当作待辨识的系统参数，一起放进参数向量里。具体来说，如果我们认为噪声是一个自回归模型，比如 v(k) = c1*v(k-1) + ... + e(k)（其中e(k)是白噪声），那么原来的模型方程就会变成一个更高阶的方程，其中既包含原来的系统输入输出，也包含不可测的噪声项 v(k)。

直接处理这个方程很麻烦，因为 v(k) 未知。ELS的诀窍在于，它用残差（实际输出与模型预测输出之差）来近似代替噪声。在递推过程中，我们把历史的残差也当作一种“虚拟的输入”，和真实的输入输出数据一起，构造出一个增广的回归向量 h_aug(k)。

# 增广最小二乘法 (ELS) 简要思路示意
# 假设系统模型为： y(k) = a*y(k-1) + b*u(k-1) + v(k)
# 噪声模型为： v(k) = c*v(k-1) + e(k) （一阶噪声）
# 则增广参数向量为： theta_aug = [a, b, c]^T
# 增广回归向量为： h_aug(k) = [y(k-1), u(k-1), epsilon(k-1)]^T
# 其中 epsilon(k-1) 是 k-1 时刻的残差（噪声估计）

# 在递推公式中，我们使用 h_aug(k) 代替原来的 h(k)
# 初始残差可以设为0，随着递推进行，用 (y(k) - h_aug(k).T @ theta) 来更新当前残差，并用于构造下一时刻的 h_aug

这样做的好处是，我们用一个统一的算法，同时辨识出了系统的动态特性和噪声的统计特性。得到的模型不仅能更准确地预测系统输出，还能对不确定性（噪声）进行量化，这对于后续设计鲁棒控制器至关重要。

当然，ELS也有它的代价。首先，待估参数变多了，需要更长的数据才能获得好的估计。其次，因为用残差近似噪声，在算法初期或噪声较大时，估计可能会有波动。我在实际应用中发现，将ELS与递推算法结合，并设置一个谨慎的遗忘因子，效果最好。它特别适用于对模型精度要求高、且噪声特性明显的场合，比如高精度伺服系统、振动信号分析等。

6. 避坑指南：工程实战中的优化策略与经验谈

理论很完美，但一脚踩进工程现实，坑可不少。这里分享几个我用最小二乘法做系统辨识时，总结出的关键优化策略和避坑经验。

策略一：激励信号的设计是成败的一半。 你不能用幅度太小或频率成分太单一的信号去激励系统。这就好比你想了解一个人的全部性格，却只跟他聊天气，得到的信息肯定片面。我强烈推荐使用伪随机二进制序列（PRBS）或扫频正弦信号作为输入u(k)。PRBS能覆盖很宽的频带，能量分布均匀，能同时激发出系统各种模态的响应。在设计时，要确保信号的幅值能覆盖系统预期的正常工作范围，但不要大到引起非线性饱和。信号的周期要足够长，以保证能观察到系统的完整暂态过程。

策略二：采样周期的选择需要权衡。 采样太快（周期太短），相邻数据间相关性太强，信息冗余，还可能引入高频测量噪声；采样太慢（周期太长），会丢失系统快速动态的信息，导致混叠。一个实用的经验法则是：采样频率至少应为系统预估最高工作频率的10倍以上。如果不确定，可以先快速采样，然后对数据进行低通滤波和下采样处理。

策略三：模型结构的先验知识至关重要。 最小二乘法只能帮你估计给定模型结构下的参数。如果模型结构选错了，比如一个二阶系统你用一阶模型去拟合，那参数估计得再准也没用。在动手辨识前，一定要通过机理分析、阶跃响应观察、频域分析等方法，对系统的阶次、有无纯延时等有一个大致的判断。可以尝试不同的模型阶次，用损失函数值（残差平方和）和AIC/BIC信息准则来辅助判断哪个模型更合适。

策略四：处理异常值与数据预处理。 传感器偶尔的跳变、通信瞬间中断，都会产生异常值。这些“野点”对基于平方误差的最小二乘法影响巨大。在构造数据矩阵H和向量z之前，一定要先做数据清洗。简单的阈值法、滑动平均滤波，或者更复杂的统计离群点检测方法，都能有效提升数据质量。别忘了，加权最小二乘本身也是一种通过权重来抑制异常值影响的高级预处理手段。

策略五：验证、验证、再验证。 辨识出的模型，绝不能只在一个数据集上看起来好看。一定要用另一组独立的、未参与辨识的数据进行验证。计算模型在这组验证数据上的预测误差。如果训练误差小，验证误差大，说明模型可能“过拟合”了，即过分迎合了训练数据中的噪声。这时需要考虑简化模型结构，或增加正则化项。

最后，工具的选择上，MATLAB的System Identification Toolbox功能非常强大，适合快速原型验证。而在嵌入式或需要在线运行的场合，用C/C++或Python手动实现递推算法（RLS或ELS）是必须的。记住，系统辨识是一个“实验-建模-验证-再实验”的迭代过程，没有一蹴而就的完美模型。每一次调试，你对系统和最小二乘法的理解都会更深一层。