增广递推最小二乘参数估计

%增广递推最小二乘参数估计(ERLS)

考虑如下系统：

y(k)1.5y(k1)0.7y(k2)u(k3)0.5u(k4)(k)(k1)0.2(k2)

式中(k)为方差为0.1的白噪声。

选择方差为1的白噪声作为输入信号u(k).

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; %计算阶次

L=1000; %数据长度

uk=zeros(d+nb,1);yk=zeros(na,1); %输入输出初值

xik=zeros(nc,1); %噪声初值

xiek=zeros(nc,1); %噪声估计初值

u=randn(L,1); %输入采用方差为1的白噪声序列

xi=sqrt(0.1)*randn(L,1); % 方差为0.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);

lambda=0.98; %遗忘因子范围[0.9 1]

for k=1:L

phi=[-yk;uk(d:d+nb);xik]; %测量向量

y(k)=phi'*theta+xi(k); %采样输出数据

phie=[-yk;uk(d:d+nb);xiek]; %估计的测量向量

%增广递推最小二乘公式

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]);