遗忘因子递推最小二乘参数估计

%遗忘因子递推最小二乘参数估计(FFRLS)
考虑如下系统：
1201()(1)(2)(3)(4)()y k a y k a y k b u k b u k k ξ+-+-=-+-+
式中ξ(k )为均值为0、方差为0.1的白噪声，
对象时变参数T 1201()[,,,]k a a b b θ=为：
[ 1.5,0.7,1,0.5],500()[1,0.4,1.5,0.2],500
T T k k k θ⎧-≤⎪=⎨->⎪⎩ 取遗忘因子λ=0.98，
采用方差为1的白噪声序列作为输入信号u (k ):
clear all; close all;
a=[1 -1.5 0.7]';b=[1 0.5]';d=3; %对象参数
na=length(a)-1;nb=length(b)-1; %计算阶次
L=1000;%数据长度
uk=zeros(d+nb,1);yk=zeros(na,1); %输入输出初值 u=randn(L,1); %输入采用方差为1的白噪声序列
xi=sqrt(0.1)*randn(L,1); % 方差为0.1的白噪声干扰序列
%theta=[a(2:na+1);b]; %对象参数真值
thetae_1=zeros(na+nb+1,1); %参数初值
P=10^6*eye(na+nb+1);
lambda=0.98; %遗忘因子范围[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(2,1,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(2,1,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]);
thetae(:,400)=
-1.5334
0.7266
0.9841
0.4743 thetae(:,900)=
-0.9731
0.3607
1.5207
0.2374。