【附代码获取链接】递推最小二乘参数辨识：原理、实现与代码详解-编程阁

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

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

在控制系统、信号处理和系统建模领域，我们经常面临这样的问题：给定一个"黑箱"系统，我们能够测量其输入和输出信号，但不知道系统内部的具体数学模型。参数辨识技术就是解决这类问题的关键工具。递推最小二乘法（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 代码文件组织

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

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

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

%% Example Recursive Least Squares %%clc close all clear all%% 系统定义Gs=tf([0.5],[111]);% 实际连续系统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]% 处理延迟ifdelay>0Bz=[zeros(1,delay),cell2mat(num)]./Az(1);elseBz=cell2mat(num)./Az(1);end

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

2.4 延迟检测算法

% 检测延迟d（Bz中前导零的数量）ind=find(Bz==0);test=isempty(ind);iftest==1d=0;elseifind(1)==1d=1;elseifind(1)~=1d=0;end% 处理多个连续零的情况iflength(ind)>1fori=1:length(ind)-1ifind(i+1)-ind(i)==1d=i+1;elsebreakendendend

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

三、RLS核心算法实现

3.1 数据向量构建（Regressor Vector）

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

数据向量φ(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的核心递推公式：

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

四、递推流程分析

4.1 时间递推结构

% 在主程序中的递推循环Theta=zeros(nu,1);% 初始参数P=10^6*eye(nu,nu);% 初始协方差矩阵Phi=zeros(1,nu);% 初始数据向量fori=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 onform=1:Sa(1)plot(t,a(m,:),'color',colors(m),'LineWidth',1.5)plot(t,A(m)*ones(length(t),1),'--','color','k')% 真实值参考线grid ontitle('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 onplot(t,Y(1:end),'--','LineWidth',1.5)% 实际系统输出plot(t,y(1:end),'o','LineWidth',1.5,'color','r')% 估计系统输出