最小二乘系统辨识：从ARX模型到参数估计的工程实践

1. 从“黑箱”到“白箱”系统辨识的工程价值在工业控制、信号处理乃至经济建模的日常工作中我们常常面对一个“黑箱”你给它一个输入信号它吐出一个输出信号但箱子内部的结构、参数你一无所知。比如你想预测一个化学反应器的温度变化或者想设计一个无人机的飞控算法你首先得知道这个系统“听不听话”——给它一个推力它会以多快的速度响应这个响应过程是平滑的还是震荡的系统辨识就是把这“黑箱”变成“白箱”的科学与艺术。它不要求你拆开物理设备而是通过输入输出数据用数学方法“猜”出系统内部的动态规律也就是数学模型。而“最小二乘”无疑是打开这扇门最经典、最实用的一把钥匙。它不是什么高深莫测的理论其核心思想朴素得惊人找一条线或一个曲面让所有数据点到这条线的“距离”平方和最小。在系统辨识里这条“线”就是我们的候选模型那些“距离”就是模型预测输出与实际观测输出之间的误差。最小二乘法告诉我们哪个模型参数能让预测和现实最贴合。这门“最小二乘系统辨识课”的上篇我们就聚焦于最基础、也最核心的部分如何为你的系统选择一个合适的“辨识模型”以及最小二乘法是如何在这个模型框架下大显身手的。无论你是自动化专业的学生还是从事算法开发的工程师理解这套基础范式都是后续处理更复杂非线性、时变系统问题的基石。2. 模型选择给系统“画像”的第一步在动手算参数之前选对模型结构至关重要。这就像画画你得先决定是用素描还是油画是写实还是抽象。选错了后续再精妙的算法也是徒劳。在经典的系统辨识中主要有两大类模型结构方程误差模型和输出误差模型。理解它们的区别是避免后续踩坑的关键。2.1 方程误差模型最直接的线性回归视角方程误差模型也叫ARX模型是入门系统辨识最先接触的结构。它的形式非常直观A(z)y(k) B(z)u(k) e(k)这里y(k)是系统输出u(k)是系统输入e(k)是一个白噪声序列。A(z)和B(z)是后移算子z^{-1}的多项式。展开写就是一个差分方程y(k) a1*y(k-1) ... ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k)这个方程告诉我们当前的输出y(k)是由过去若干时刻的输出y(k-1)...、过去若干时刻的输入u(k-1)...以及一个当前的噪声e(k)共同决定的。为什么它如此受欢迎核心原因在于这个模型关于待估参数a1, a2, ..., b1, b2, ...是线性的。我们把方程稍微变形一下y(k) -a1*y(k-1) - ... - ana*y(k-na) b1*u(k-1) ... bnb*u(k-nb) e(k)现在我们把等号右边除了e(k)之外的所有项看作是一组已知的“特征”-y(k-1),-y(k-2), ...,u(k-1),u(k-2)...。而a1, a2, ..., b1, b2, ...就是这些特征对应的权重系数。看这完美地契合了多元线性回归y θ1*x1 θ2*x2 ... θn*xn e的形式这意味着我们可以直接套用成熟、高效、解析解明确的最小二乘法来估计参数计算简单全局最优。它的“阿喀琉斯之踵”噪声假设方程误差模型有一个很强的假设噪声e(k)是直接加在方程等式两端的。在实际的物理系统中噪声往往不是这样作用的。更常见的场景是噪声影响了系统内部状态然后经过系统自身的动力学特性后才体现在输出上。这就引出了方程误差模型的一个主要缺点当真实系统的噪声通道与模型假设不符时即使数据量无限其参数估计值也会存在偏差无法收敛到真实值。这种现象在辨识理论中称为“有偏估计”。因此ARX模型更适合对模型精度要求不是极端苛刻或者噪声影响较小的初步建模场景。2.2 输出误差模型更贴近物理现实的描述为了克服方程误差模型的偏差问题输出误差模型OE模型被提出。它的结构更贴近我们对许多物理系统的认知y(k) [B(z)/F(z)] * u(k) e(k)这里[B(z)/F(z)]代表一个传递函数描述了输入u到系统“无噪声输出”x(k)的动态过程。而最终的观测输出y(k)是这个无噪声输出x(k)加上一个独立的观测噪声e(k)。用差分方程写出来是x(k) f1*x(k-1) ... fnf*x(k-nf) b1*u(k-1) ... bnb*u(k-nb)y(k) x(k) e(k)模型优势与代价OE模型清晰地分离了系统的确定性动态B/F和随机性干扰e(k)。只要e(k)是白噪声理论上通过一些迭代算法如预报误差法可以得到无偏的、一致的参数估计。这听起来很完美对吧但代价是模型关于参数b1, b2, ..., f1, f2, ...是非线性的。因为无噪声输出x(k)本身依赖于参数f这使得模型无法写成关于所有参数的线性回归形式。我们不能再用简单的最小二乘直接求解而必须借助非线性优化算法如高斯-牛顿法、梯度下降法进行迭代搜索计算复杂且可能陷入局部最优。注意在实际工程中模型结构的选择ARX还是OE以及阶次na, nb, nf如何确定本身就是一个重要课题。通常建议从简单的ARX模型开始利用其计算快的优势进行模型阶次的试探和初步分析如果残差序列预测误差表现出明显的相关性说明ARX的噪声假设不成立再考虑使用OE等更复杂的模型。不要一开始就追求复杂模型简单有效永远是第一原则。3. 最小二乘原理误差平方和最小化的几何与统计意义现在我们以最经典的ARX模型为例深入最小二乘法的内核。假设我们通过实验采集到了一组长度为N的输入输出数据{u(1), y(1)}, {u(2), y(2)}, ..., {u(N), y(N)}。我们选定模型阶次na和nb想要估计参数向量θ [a1, a2, ..., ana, b1, b2, ..., bnb]^T。根据ARX模型对于第k个采样时刻k max(na, nb)我们可以写出y(k) φ(k)^T θ e(k)其中φ(k)称为回归向量或数据向量φ(k) [-y(k-1), -y(k-2), ..., -y(k-na), u(k-1), u(k-2), ..., u(k-nb)]^T你看这个式子把非线性关于y和u的系统动力学巧妙地转化成了关于参数θ的线性形式。对于所有k L1 到 NL max(na, nb)我们可以构建一个庞大的线性方程组y(L1) φ(L1)^T θ e(L1) y(L2) φ(L2)^T θ e(L2) ... y(N) φ(N)^T θ e(N)写成矩阵形式简洁有力Y Φ θ E其中Y [y(L1), y(L2), ..., y(N)]^T是输出向量。Φ [φ(L1), φ(L2), ..., φ(N)]^T是数据矩阵也叫信息矩阵或回归矩阵。E [e(L1), e(L2), ..., e(N)]^T是噪声向量。我们的目标是找到一组参数θ使得模型预测值Φθ尽可能接近真实观测值Y。最小二乘准则规定这个“接近”的程度用所有误差的平方和来衡量即代价函数J(θ) E^T E (Y - Φθ)^T (Y - Φθ)。我们要找到使J(θ)最小的那个θ。从几何角度理解向量Y存在于一个N维空间中。矩阵Φ的列向量张成了一个子空间列空间。模型预测值Φθ是这个子空间中的一个点。最小二乘解的意义在于在子空间中寻找一个点Φθ使得它到真实点Y的欧几里得距离最短。根据几何知识这个最短距离是通过Y向子空间做正交投影得到的。因此最优的Φθ是Y在Φ列空间上的投影而误差向量E Y - Φθ垂直于该列空间。从统计角度理解如果我们假设噪声e(k)是零均值、同方差且互不相关的白噪声那么最小二乘估计量具有一系列优良性质它是无偏的估计值的期望等于真值并且是所有无偏估计中方差最小的有效估计。这使得最小二乘在统计意义上也是最优的。求解这个最优化问题可以通过对代价函数J(θ)求关于θ的梯度并令其为零∇J(θ) -2Φ^T (Y - Φθ) 0由此得到著名的正规方程(Φ^T Φ) θ Φ^T Y只要矩阵Φ^T Φ是可逆的这要求数据足够丰富且输入信号具有一定的激励性即持续激励条件我们就可以得到最小二乘参数估计的解析解θ_LS (Φ^T Φ)^{-1} Φ^T Y这个公式干净利落是所有系统辨识、机器学习线性回归问题的基石。4. 实操核心数据矩阵构建与持续激励条件理论很优美但落地到代码和实验有两个细节决定成败如何正确构建数据矩阵Φ以及如何保证Φ^T Φ可逆。很多初学者在这里栽跟头。4.1 数据矩阵构建的“边界”问题构建Φ矩阵时第一个实际问题是数据索引的起始点。注意我们的回归向量φ(k)包含了y(k-1), ..., y(k-na)和u(k-1), ..., u(k-nb)。这意味着要构造φ(L1)我们需要y(L)和u(L)的数据而L max(na, nb)。因此有效的数据段是从k L1开始到k N结束。你用于拟合的数据长度实际上是N - L而不是N。在编程时一个常见的错误是直接从k1开始循环构造φ(k)导致数组越界。正确的做法是import numpy as np # 假设 y_data, u_data 是长度为 N 的数组 na, nb 2, 2 L max(na, nb) N len(y_data) Y y_data[L:] # 从索引L到末尾 Phi [] for k in range(L, N): phi_k np.concatenate([-y_data[k-1:k-na-1:-1], u_data[k-1:k-nb-1:-1]]) Phi.append(phi_k) Phi np.array(Phi) # 然后求解 theta np.linalg.inv(Phi.T Phi) Phi.T Y这里y_data[k-1:k-na-1:-1]利用了Python切片来获取[y(k-1), y(k-2), ..., y(k-na)]注意顺序。4.2 持续激励让数据“会说话”第二个也是更本质的问题是什么样的输入数据u(k)才能让我们唯一地、可靠地辨识出所有参数答案就是输入信号必须满足“持续激励”条件。直观理解如果你用一个恒定值比如u(k)1去激励系统你只能得到系统在某个静态工作点附近的信息无法分辨出系统动态a1, a2,...中不同模式的影响。这就像你想了解一个弹簧的质量和阻尼系数只把它压住不动是测不出来的必须用不同频率的力去推拉它。从数学上看Φ^T Φ的可逆性要求数据矩阵Φ是列满秩的。对于ARX模型Φ的列由过去的输出和输入数据组成。而过去的输出y(k-i)本身又依赖于过去的输入和参数。因此归根结底要求输入信号u(k)能充分激发系统所有模态。一个经典且实用的选择是伪随机二进制序列。它看起来像随机的0/1跳变但具有周期性和良好的自相关特性能在一个宽频带内提供近似白噪声的激励是系统辨识实验中常用的输入信号。实操心得在仿真中你可以用PRBS或白噪声作为输入。但在实际物理系统测试中输入信号必须考虑系统的安全限幅和执行器的物理限制。一个折中的好办法是使用幅值受限的、不同频率的正弦扫频信号或者幅值随机变化的阶跃信号序列。关键是要让输入有足够丰富的变化覆盖你关心的频率范围。记录数据时务必确保输入输出数据是同步采集的并注意剔除明显的野值。5. 评估与诊断你的模型“合格”了吗参数θ算出来了模型就有了。但模型质量如何不能只靠感觉需要有量化的评估和诊断工具。这里介绍三个最实用的方法。5.1 拟合优度量化匹配程度最直接的指标是拟合优度通常用归一化的均方误差来表示例如FIT (1 - norm(Y - Y_pred) / norm(Y - mean(Y))) * 100%其中Y_pred Φ θ_LS是模型的预测输出。FIT越接近100%说明模型对这段训练数据的解释能力越强。但要注意高拟合优度可能意味着过拟合——模型不仅拟合了系统动态还拟合了数据中特定的噪声。因此它通常用于同一组数据上不同模型结构的横向比较而不是绝对质量的评判。5.2 残差分析检验模型假设的“试金石”残差就是预测误差序列ε(k) y(k) - y_pred(k)。如果我们的模型包括其噪声假设是完美的那么残差序列应该是一个白噪声序列——均值为零序列自身不同时刻的值互不相关。如何检验绘制残差序列图肉眼观察是否围绕零均值线随机波动有无明显的趋势或周期性。计算残差的自相关函数对于一个白噪声其自相关函数在时滞τ≠0时应接近于零。我们可以计算残差的自相关函数并观察其是否落在95%的置信区间内。如果很多点落在区间外特别是前几个时滞的相关性显著不为零则说明残差中存在未建模的动态信息模型结构如阶次可能选择不当。残差与输入信号的互相关函数一个理想的模型其残差应与过去的输入信号不相关。如果互相关函数显著不为零说明模型未能完全捕获输入到输出的动态关系或者存在非线性未建模。残差分析是系统辨识中极其重要的一步它比单纯的拟合优度更能揭示模型的本质问题。我个人的经验是一个FIT只有85%但残差接近白噪声的模型通常比一个FIT高达95%但残差自相关严重的模型更可靠、更具泛化能力。5.3 交叉验证防范过拟合的黄金准则最终极的检验是将模型用在它“没见过”的数据上。这就是交叉验证。具体做法将你的数据集分为两部分一部分用于参数估计训练集另一部分用于模型验证测试集。只用训练集的数据来构建Φ_train和Y_train并计算参数θ_LS。锁定这个θ_LS将其应用到测试集上。用测试集的输入数据和模型公式注意此时计算预测输出y_pred_test时使用的过去输出值应是测试集真实的y值而不是模型自己递归预测的值这称为“仿真验证”模式得到测试集上的预测输出。计算模型在测试集上的拟合优度或误差指标。如果模型在训练集上表现很好但在测试集上表现大幅下降这就是典型的过拟合。交叉验证能最真实地反映模型对新数据的预测能力是评估模型泛化性能的黄金标准。在实际项目中我强烈建议至少保留30%的数据作为测试集。6. 一个完整的仿真案例二阶离散系统辨识让我们用一个具体的MATLAB/Python仿真例子把上述所有步骤串起来。假设真实系统是一个二阶离散系统y(k) - 1.5y(k-1) 0.7y(k-2) u(k-1) 0.5u(k-2) e(k)其中e(k)是方差为0.1的高斯白噪声。步骤1数据生成% MATLAB N 1000; u randn(N, 1); % 使用白噪声作为输入满足持续激励 e 0.1 * randn(N, 1); y zeros(N,1); y(1:2) [0; 0]; % 初始条件 for k 3:N y(k) 1.5*y(k-1) - 0.7*y(k-2) u(k-1) 0.5*u(k-2) e(k); end % 划分训练集和测试集 train_ratio 0.7; N_train floor(N * train_ratio); u_train u(1:N_train); y_train y(1:N_train); u_test u(N_train1:end); y_test y(N_train1:end);# Python import numpy as np N 1000 np.random.seed(42) u np.random.randn(N) # 白噪声输入 e 0.1 * np.random.randn(N) y np.zeros(N) y[:2] [0, 0] for k in range(2, N): y[k] 1.5*y[k-1] - 0.7*y[k-2] u[k-1] 0.5*u[k-2] e[k] # 划分数据集 train_ratio 0.7 N_train int(N * train_ratio) u_train, y_train u[:N_train], y[:N_train] u_test, y_test u[N_train:], y[N_train:]步骤2模型阶次选择与数据矩阵构建我们猜测系统可能是二阶的na2, nb2。用训练集数据构建ARX模型的数据矩阵。na, nb 2, 2 L max(na, nb) # 构建训练集数据矩阵 Y_train y_train[L:] Phi_train [] for k in range(L, len(y_train)): phi_k np.concatenate([-y_train[k-1:k-na-1:-1], u_train[k-1:k-nb-1:-1]]) Phi_train.append(phi_k) Phi_train np.array(Phi_train)步骤3最小二乘参数估计# 求解正规方程 theta_LS np.linalg.inv(Phi_train.T Phi_train) Phi_train.T Y_train print(f估计参数: {theta_LS}) print(f真实参数: a1-1.5, a20.7, b11.0, b20.5) # 注意我们回归方程写的是 y(k) -a1*y(k-1) -a2*y(k-2) b1*u(k-1) b2*u(k-2) e(k) # 所以 theta_LS 的前两个元素对应 -a1, -a2后两个对应 b1, b2 a1_est, a2_est, b1_est, b2_est -theta_LS[0], -theta_LS[1], theta_LS[2], theta_LS[3] print(f估计的差分方程: y(k) {a1_est:.4f}*y(k-1) {a2_est:.4f}*y(k-2) {b1_est:.4f}*u(k-1) {b2_est:.4f}*u(k-2))运行后估计参数应该非常接近真实值[-1.5, 0.7, 1.0, 0.5]但由于噪声存在会有微小偏差。步骤4模型评估与诊断# 1. 训练集拟合 y_pred_train Phi_train theta_LS fit_train 100 * (1 - np.linalg.norm(y_train[L:] - y_pred_train) / np.linalg.norm(y_train[L:] - np.mean(y_train[L:]))) print(f训练集拟合优度: {fit_train:.2f}%) # 2. 测试集验证 (仿真模式) # 注意测试集仿真需要递归计算因为每一步的预测输出会作为下一步的输入 y_pred_test np.zeros_like(y_test) y_pred_test[:L] y_test[:L] # 用真实值初始化 for k in range(L, len(y_test)): # 使用模型和过去的“预测输出”及“真实输入” phi_k_test np.concatenate([-y_pred_test[k-1:k-na-1:-1], u_test[k-1:k-nb-1:-1]]) y_pred_test[k] phi_k_test theta_LS fit_test 100 * (1 - np.linalg.norm(y_test[L:] - y_pred_test[L:]) / np.linalg.norm(y_test[L:] - np.mean(y_test[L:]))) print(f测试集拟合优度: {fit_test:.2f}%) # 3. 残差分析 (以训练集残差为例) residual y_train[L:] - y_pred_train # 计算残差自相关函数 (这里简化为计算前20个时滞) max_lag 20 acf np.correlate(residual, residual, modefull) acf acf[len(acf)//2 : len(acf)//2 max_lag 1] acf acf / acf[0] # 归一化 # 绘制自相关图 (略)观察是否在置信带内。通过这个完整流程你不仅得到了模型参数还通过测试集验证和残差分析对模型质量有了定量和定性的评估。如果测试集拟合度尚可且残差接近白噪声那么这个模型就可以用于后续的控制器设计或系统分析了。7. 局限与进阶思考最小二乘辨识的边界通过上篇的讨论我们掌握了基于最小二乘和ARX模型进行系统辨识的完整流程。然而在实际工程中你会很快遇到它的边界。认识到这些局限正是迈向更高级辨识方法如递推最小二乘、广义最小二乘、预报误差法等的起点。首先是噪声模型的局限。ARX模型假设噪声直接加在方程输出端这在实际中往往不成立。当噪声通过系统动力学环节即有色噪声时标准最小二乘估计是有偏的。这时就需要引入更复杂的模型如ARMAX方程误差模型但噪声为移动平均过程或BJBox-Jenkins模型并使用广义最小二乘或极大似然等方法来获得无偏估计。其次是时变系统的挑战。我们假设系统参数是定常的。但对于缓慢变化或突变的系统如化学反应器催化剂活性衰减、飞机在不同空速下的气动参数需要能够在线跟踪参数变化的算法这就是递推最小二乘RLS及其变种如带遗忘因子的RLS的用武之地。再者是关于线性与静态的假设。最小二乘辨识本质上是线性回归它只能辨识线性系统的参数。对于非线性系统你需要选择非线性模型结构如Hammerstein模型、Wiener模型、神经网络等并使用非线性优化方法进行参数估计计算复杂度和初始值敏感性会急剧增加。最后我想分享一个最容易被忽视的心得系统辨识不仅仅是数学和算法更是对物理对象的理解与实验设计的艺术。再精巧的算法如果输入数据质量差激励不足、测量噪声大、采样不同步也得不到好模型。在启动辨识算法之前请务必花时间思考我的输入信号能激发系统的所有重要模式吗我的采样频率足够高吗满足香农定理我的传感器数据可靠吗很多时候花在改进实验设计和数据预处理上的时间远比调试算法参数带来的回报大得多。最小二乘给了我们一个强大的工具但用好这个工具的前提是深刻理解你的系统和你的数据。