递推最小二乘参数辨识：原理、实现与代码详解

递推最小二乘参数辨识：原理、实现与代码详解

引言：为什么需要参数辨识？

在控制系统、信号处理和系统建模领域，我们经常面临这样的问题：给定一个"黑箱"系统，我们能够测量其输入和输出信号，但不知道系统内部的具体数学模型。参数辨识技术就是解决这类问题的关键工具。递推最小二乘法（Recursive Least Squares, RLS）作为最经典的在线辨识算法之一，以其计算效率高、适用于实时系统等优点，在工程实践中得到广泛应用。

本文将深入解析递推最小二乘法的原理，并通过三个完整的MATLAB代码示例，详细展示其实现过程。

一、递推最小二乘法基本原理

1.1 系统模型表示

考虑离散时间系统，通常可用以下差分方程描述：

y(k) + a₁y(k-1) + ... + a_na y(k-na) = 
b₀u(k-d) + b₁u(k-d-1) + ... + b_nb u(k-d-nb) + e(k)

其中：

• y(k) 是系统在时刻k的输出
• u(k) 是系统在时刻k的输入
• d 是系统延迟
• e(k) 是测量噪声
• a₁,..., a_na 和 b₀,..., b_nb 是待辨识的系统参数

写成向量形式：

y(k) = φ(k)ᵀθ + e(k)

其中：

φ(k) = [-y(k-1), ..., -y(k-na), u(k-d), ..., u(k-d-nb)]ᵀ
θ = [a₁, ..., a_na, b₀, ..., b_nb]ᵀ

1.2 最小二乘法的递推形式

传统的最小二乘法需要处理所有历史数据，计算量大且不适合在线应用。递推最小二乘法通过迭代更新，每次只处理新数据，大大提高了计算效率。

递推公式如下：

K(k) = P(k-1)φ(k)[1 + φ(k)ᵀP(k-1)φ(k)]⁻¹
θ̂(k) = θ̂(k-1) + K(k)[y(k) - φ(k)ᵀθ̂(k-1)]
P(k) = P(k-1) - P(k-1)φ(k)[1 + φ(k)ᵀP(k-1)φ(k)]⁻¹φ(k)ᵀP(k-1)

其中：

• K(k) 是增益矩阵
• P(k) 是协方差矩阵
• θ̂(k) 是k时刻的参数估计值

二、代码结构解析

2.1 代码文件组织

提供的代码包含三个文件：

1. example.m - 示例程序，展示RLS在仿真系统中的应用
2. main.m - 主程序，适用于实际数据处理
3. RecursiveLeastSquares.m - RLS核心算法函数

2.2 系统建模与初始化（example.m分析）

%% Example Recursive Least Squares %%
clc
close all
clear all

%% 系统定义
Gs = tf([0.5],[1 1 1]);                     % 实际连续系统
Ts = 0.2;                                   % 采样时间
Gz = c2d(Gs,Ts);                           % 离散化

这部分代码创建了一个二阶连续系统，并将其离散化。离散化是数字控制系统分析和设计的关键步骤。

2.3 阶次确定与参数提取

%% 阶次计算
num = get(Gz,'num');
den = get(Gz,'den');
delay = get(Gz,'iodelay');

% 提取多项式系数
Az = cell2mat(den);                  
Az = Az./Az(1);                             % 归一化：[1 a1 a2 ... ana]

% 处理延迟
if delay >0
    Bz = [zeros(1,delay),cell2mat(num)]./Az(1);
else
    Bz = cell2mat(num)./Az(1);
end

这段代码展示了如何从传递函数对象中提取分子、分母多项式和延迟信息，为参数辨识做准备。

2.4 延迟检测算法

% 检测延迟d（Bz中前导零的数量）
ind = find(Bz == 0);
test = isempty(ind);

if test == 1
    d = 0;
elseif ind(1) == 1
    d = 1;
elseif ind(1)~= 1 
    d = 0;
end

% 处理多个连续零的情况
if length(ind)>1
    for i=1:length(ind)-1
        if ind(i+1)-ind(i) == 1
           d = i+1;
        else
            break
        end
    end
end

这是一个实用的延迟检测算法，通过查找B(z)多项式中的前导零来确定系统延迟。

三、RLS核心算法实现

3.1 数据向量构建（Regressor Vector）

function [a,b,P,Theta,phi,K] = RecursiveLeastSquares(U,Y,d,nb,na,P,Theta,phi,n)
% 构建数据向量φ
for j = 1:nu
    if j <= na % y的项
        if (n-j)<=0
            phi(n,j) = 0;           % 初始时刻，历史数据不足时补零
        else
            phi(n,j) = -Y(n-j);     % 输出历史值
        end
    else       % u的项
        if (n-d-(j-(na+1)))<=0
            phi(n,j) = 0;           % 考虑延迟，历史数据不足时补零
        else
            phi(n,j) = U(n-d-(j-(na+1)));  % 输入历史值
        end
    end
end

数据向量φ(k)的构建是RLS算法的关键步骤之一。它包含了系统的历史输入输出信息，按照模型的阶次和延迟排列。

3.2 RLS递推更新核心

% RLS递推更新
K = P*phi(n,:)'*inv(1+phi(n,:)*P*phi(n,:)');      % 计算增益矩阵
Theta = Theta + K*(Y(n)-phi(n,:)*Theta);          % 参数更新
P = P - P*phi(n,:)'*inv(1+phi(n,:)*P*phi(n,:)')*phi(n,:)*P;  % 协方差矩阵更新

% 提取参数
a = Theta(1:na);      % 分母多项式系数
b = Theta(na+1:end);  % 分子多项式系数

这三行代码实现了RLS的核心递推公式：

1. 增益计算：确定新数据对参数更新的"权重"
2. 参数更新：基于预测误差修正参数估计
3. 协方差更新：反映参数估计的不确定性

四、递推流程分析

4.1 时间递推结构

% 在主程序中的递推循环
Theta = zeros(nu,1);                      % 初始参数
P = 10^6 * eye(nu,nu);                    % 初始协方差矩阵
Phi = zeros(1,nu);                        % 初始数据向量

for i = 1 : length(U)  % 时间递推循环
    [a(:,i),b(:,i),P,Theta,Phi,K(:,i)] = ...
        RecursiveLeastSquares(U(1:i),Y(1:i),d,nb,na,P,Theta,Phi,i);
end

这个for循环体现了RLS的时间递推特性：

• 每个时间步处理一组新的输入输出数据
• 参数估计不断更新，无需保存所有历史数据
• 适合在线、实时应用

4.2 协方差矩阵的初始化

P = 10^6 * eye(nu,nu);  % 初始协方差矩阵

这里将协方差矩阵初始化为一个很大的单位矩阵，反映了对初始参数估计的高度不确定性。随着数据积累，P矩阵会逐渐减小，表明参数估计越来越确定。

五、参数辨识结果与验证

5.1 参数收敛性分析

% 参数收敛图
figure
hold on
for m = 1:Sa(1)
    plot(t,a(m,:),'color',colors(m),'LineWidth',1.5)
    plot(t,A(m)*ones(length(t),1),'--','color','k')  % 真实值参考线
    grid on
    title('Parameters a_i Convergence')
    xlabel('time (sec.)')
    ylabel('Parameter Value')
end

代码绘制了参数a_i随时间的变化曲线，并与真实值（虚线）对比，直观展示RLS算法的收敛性能。

5.2 模型验证

% 使用辨识得到的参数构建估计系统
y = lsim(tf(b(:,end)',[1 a(:,end)'],Ts,'iodelay',d,'variable','z^-1'),U,t);

% 对比实际输出与估计输出
figure
grid on 
hold on
plot(t,Y(1:end),'--','LineWidth',1.5)      % 实际系统输出
plot(t,y(1:end),'o','LineWidth',1.5,'color','r')  % 估计系统输出

通过比较实际系统输出和基于辨识参数的系统输出，可以验证辨识结果的质量。

六、实际应用注意事项

6.1 数据预处理

在实际应用（main.m）中，需要注意：

1. 数据质量：确保输入输出数据无异常值
2. 持续激励：输入信号应包含足够频率成分以激励所有模态
3. 噪声处理：考虑测量噪声的影响，可能需要改进算法

6.2 算法改进方向

基本RLS算法可能存在以下问题：

1. 数据饱和：老数据权重与新数据相同，可能不适应时变系统
2. 数值稳定性：协方差矩阵可能失去正定性

改进方法包括：

• 带遗忘因子的RLS：λRLS，给新数据更高权重
• 平方根RLS：提高数值稳定性
• 正则化：防止过拟合

七、RLS在控制系统中的应用价值

7.1 自适应控制

RLS常用于自适应控制系统，如模型参考自适应控制（MRAC）和自校正调节器。通过在线辨识系统参数，控制器可以实时调整以应对系统变化。

7.2 故障检测与诊断

通过监控辨识参数的变化，可以检测系统故障。例如，机械系统的刚度或阻尼参数异常变化可能预示故障。

7.3 系统健康监测

在关键设备（如航空发动机、电力变压器）中，RLS可用于监测系统参数漂移，实现预测性维护。

结论

递推最小二乘法是一种强大而灵活的系统辨识工具。通过本文的代码解析，我们可以看到：

1. 原理与实现对应：RLS的数学公式直接对应代码实现
2. 递推特性明确：通过时间循环实现参数在线更新
3. 实用性高：代码结构清晰，易于修改适应不同应用

提供的三个MATLAB文件构成了一个完整的RLS辨识框架：

• example.m 提供了仿真示例和验证方法
• main.m 展示了实际数据处理流程
• RecursiveLeastSquares.m 是核心算法，可独立使用

理解这些代码不仅有助于掌握RLS算法，也为进一步学习更先进的自适应滤波和系统辨识技术奠定了基础。在实际工程应用中，可以根据具体需求调整算法参数和改进实现方式，以获得更好的辨识性能。

---

版权声明：本文代码示例仅供学习使用，实际应用时请根据具体需求进行修改和优化。递推最小二乘法是系统辨识领域的经典方法，掌握其原理和实现对控制工程师至关重要。