一文理清最小二乘法估计

大家好，欢迎来到IT知识分享网。

1 最小二乘法估计(LS)

1.1 原理与推导

最小二乘法最早是高斯在预估星体轨道时提出来的，后来成为了估计理论的奠基石。考虑如下CAR模型：

其中：

 参数估计的任务就是根据输入和输出，估计出a1,a2,—-,ana,b1,b2,…,bnb这na+nb+1个参数。

将1-1式改成差分方程形式：

 对于L组输入{y(k),u(k),k=1,2,…,L}，系统参数的最小二乘估计为：

其中：

上式推导过程为：

对于第k次测量，其残差为：

构造如下性能指标J:

 当J取极小值时，将有：

此时有：

 1.2 例子

考虑如下系统：

 Matlab程序为：

%最小二乘参数估计（LS） clear all; a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数 na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次 L=400; %数据长度 uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值 x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值 xi=sqrt(1)*randn(L,1); %白噪声序列 theta=[a(2:na+1);b]; %对象参数真值 for k=1:L phi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据 M=xor(x3,x4); %产生M序列 IM=xor(M,S); %产生逆M序列 if IM==0 u(k)=-1; else u(k)=1; end S=not(S); %产生方波 %更新数据 x4=x3; x3=x2; x2=x1; x1=M; for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); end thetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae（结果见MATLAB命令窗口）

 运行结果：

thetae =

 2 递推最小二乘法(RLS)

 递推最小二乘法RLS是最小二乘法LS的改进版，它可以根据采样的实时数据来不断修正估计参数，从而做到参数的实时估计。其基本思想是：

RLS的实施步骤为：

对于上节中的例子，Matlab程序：

%递推最小二乘参数估计（RLS） clear all; close all; a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数 na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次 L=400; %仿真长度 uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %白噪声序列 theta=[a(2:na+1);b]; %对象参数真值 thetae_1=zeros(na+nb+1,1); %thetae初值 P=10^6*eye(na+nb+1); for k=1:L phi=[-yk;uk(d:d+nb)]; %此处phi为列向量 y(k)=phi'*theta+xi(k); %采集输出数据 %递推最小二乘法 K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(na+nb+1)-K*phi')*P; %更新数据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); end plot([1:L],thetae); %line([1,L],[theta,theta]); xlabel('k'); ylabel('参数估计a、b'); legend('a_1','a_2','b_0','b_1'); axis([0 L -2 1.5]);

 运行结果：

3 遗忘因子最小二乘法(FFRLS)

3.1 FFRLS原理

对于参数时变系统，RLS易出现数据饱和现象，就是随着时间的变化，矫正矩阵P(k)和K(k)越来越小，从而其矫正作用越来越弱，这将导致较大的估计误差。遗忘因子最小二乘法可以较好的处理这个问题。

FFRLS的实施步骤为：

3.2 例子

 考虑系统

Matlab程序为：

%遗忘因子递推最小二乘参数估计（FFRLS） clear all; close all; a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数 na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次 L=2000; %仿真长度 uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %白噪声序列 thetae_1=zeros(na+nb+1,1); %thetae初值 P=10^6*eye(na+nb+1); lambda=1; %遗忘因子范围[0.9 1] for k=1:L if k==501 a=[1 -1 0.4]';b=[1.5 0.2]'; %对象参数突变 end theta(:,k)=[a(2:na+1);b]; %对象参数真值 phi=[-yk;uk(d:d+nb)]; y(k)=phi'*theta(:,k)+xi(k); %采集输出数据 %遗忘因子递推最小二乘法 K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(na+nb+1)-K*phi')*P/lambda; %更新数据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); end subplot(1,2,1) plot([1:L],thetae(1:na,:)); hold on; plot([1:L],theta(1:na,:),'k:'); xlabel('k'); ylabel('参数估计a'); legend('a_1','a_2'); axis([0 L -2 2]); subplot(1,2,2) plot([1:L],thetae(na+1:na+nb+1,:)); hold on; plot([1:L],theta(na+1:na+nb+1,:),'k:'); xlabel('k'); ylabel('参数估计b'); legend('b_0','b_1'); axis([0 L -0.5 2]);

 当设置遗忘因子为1时，参数估计如图所示有比较大的误差：

当设置遗忘因子为0.99时，参数估计如图所示：

4 递推增广最小二乘法(RELS)

4.1 RELS原理

当系统模型的干扰为有色噪声时，系统方程表示为

 写成最小二乘形式为：

增广最小二乘法 实施步骤：

 4.2 例子

.考虑系统

 Matlab程序：

%递推增广最小二乘参数估计（RELS） clear all; close all; a=[1 -1.5 0.7]'; b=[1 0.5]'; c=[1 -1 0.2]'; d=3; %对象参数 na=length(a)-1; nb=length(b)-1; nc=length(c)-1; %na、nb、nc为A、B、C阶次 L=1000; %仿真长度 uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值 xik=zeros(nc,1); %噪声初值 xiek=zeros(nc,1); %噪声估计初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %白噪声序列 theta=[a(2:na+1);b;c(2:nc+1)]; %对象参数 thetae_1=zeros(na+nb+1+nc,1); %na+nb+1+nc为辨识参数个数 P=10^6*eye(na+nb+1+nc); for k=1:L phi=[-yk;uk(d:d+nb);xik]; y(k)=phi'*theta+xi(k); %采集输出数据 phie=[-yk;uk(d:d+nb);xiek]; %组建phie %递推增广最小二乘法 K=P*phie/(1+phie'*P*phie); thetae(:,k)=thetae_1+K*(y(k)-phie'*thetae_1); P=(eye(na+nb+1+nc)-K*phie')*P; xie=y(k)-phie'*thetae(:,k); %白噪声的估计值 %更新数据 thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end yk(1)=y(k); for i=nc:-1:2 xik(i)=xik(i-1); xiek(i)=xiek(i-1); end xik(1)=xi(k); xiek(1)=xie; end figure(1) plot([1:L],thetae(1:na,:)); xlabel('k'); ylabel('参数估计a'); legend('a_1','a_2'); axis([0 L -2 2]); figure(2) plot([1:L],thetae(na+1:na+nb+1,:)); xlabel('k'); ylabel('参数估计b'); legend('b_0','b_1'); axis([0 L 0 1.5]); figure(3) plot([1:L],thetae(na+nb+2:na+nb+nc+1,:)); xlabel('k'); ylabel('参数估计c'); legend('c_1','c_2'); axis([0 L -2 2]);

 运行结果

参考：系统辨识与自适应控制Matlab仿真（第3版），庞中华，崔红编著