rls.m

利用RLS算法实现参数的自适应估计 估计y(n) = 1.5*y(n-1)-0.7*y(n-2)+0.3*u(n-2)中的参数

%参数估计习题y(n) = 1.5*y(n-1)-0.7*y(n-2)+0.3*u(n-2)
%author:杨德伟
%Email：[e-mail]
%date:2008.3.24
%version2:由刘畅增加了简化算法

clc;close all;clear;
N_FIR = 500;%有限脉冲响应的长度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%画出标准的模型输出曲线
u(1:2) = 0;
u(3:N_FIR) = 1.0;%u=1
%tn = [0,0,0.2:0.2:0.2*(N_FIR-2)];%u = 1.0+sin(0.2*t)
%u = sin(tn)+1;
%u=sin(tn);
plot(u,'b');hold on;
irpt = normrnd(0,0.5,1,N_FIR);%干扰信号e(t)=N(0,1)
y(1:2) = 0;
ExOrd = 1;%选择练习的题号
switch ExOrd
    case 1
        for k = 3:N_FIR
            y(k) = 1.5*y(k-1)-0.7*y(k-2)+0.3*u(k-2);
        end
    case 2
        for k = 3:N_FIR
            y(k) = 1.5*y(k-1)-0.7*y(k-2)+0.3*u(k-2)+irpt(k-2);
        end
end
plot(y,'g');hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算回归向量
N_Count = 500;
phy_S = zeros(3,N_Count);%phy_S的各列元素为小φ'（t）;即phy_S为大Φ’（t）
for k = 3:N_Count; %用于参数估计的数据点
    phy_S(:,k) = [-y(k-1),-y(k-2),u(k-2)]';%phy_S(:,k-1)就不能收敛到指定的值。why?
end
%phy = zeros(3,N_Count-2);
phy = phy_S(:,3:N_Count);%取有用的N-2列元素进行运算。维数3*N_Count-2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1:离线方式估计参数
Y = y(:,3:N_Count)';%维数(N_Count-2)*1
theta_Offline = inv(phy*phy')*phy*Y;
%theta_Offline = inv(phy_S(:,10:N_Count-10)*(phy_S(:,10:N_Count-10))')*phy_S(:,10:N_Count-10)*Y(10:N_Count-10,:);%选取大Φ（t）中的部分行进行离线运算
for k = 3:N_FIR   %绘制估计参数的系统响应曲线
    y1(k) =-theta_Offline(1,1)*y(k-1)-theta_Offline(2,1)*y(k-2)+theta_Offline(3,1)*u(k-2);
end
%plot(y1,'c');hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2：利用迭代算法计算theta_Cap
P = zeros(3,3,N_Count);%初始化P（t）矩阵3*3维
theta_Cap = zeros(3,N_Count);%初始化θcap（t）矩阵3*1维
%theta_Cap(:,2) = [5,7,3]';%设计θcap（t）的初值对其收敛结果无影响。
K = zeros(3,N_Count);%初始化K（t）矩阵3*1维
P(:,:,2) = 1000*eye(3); %选定一初值P，为正定的。

for k = 3:N_Count; 
    K(:,k) = P(:,:,k-1)*phy_S(:,k)/(1+(phy_S(:,k))'*P(:,:,k-1)*phy_S(:,k));
    P(:,:,k) =  (eye(3)-K(:,k)*(phy_S(:,k))')*P(:,:,k-1);
    theta_Cap(:,k) = theta_Cap(:,k-1)+K(:,k)*(y(k)-(phy_S(:,k))'*theta_Cap(:,k-1));
    y2(k) = -theta_Cap(1,k-1)*y(k-1)-theta_Cap(2,k-1)*y(k-2)+theta_Cap(3,k-2);
end
plot(y2,'r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3:投影算法
theta_Cap1 = zeros(3,N_Count);%初始化θcap（t）1矩阵3*1维
K = zeros(3,N_Count);%初始化K（t）矩阵3*1维
a=0.1;b=0.4;%a的取值影响收敛速度，在无干扰情况下，a=0时收敛较快，且较精确。a取0.2时收敛过程就比较震荡了，慢，但最终可以收敛。a取0.5收敛性就已经不好了。同时，b的取值也影响收敛速度
for k=3:N_Count
    K(:,k) = b*phy_S(:,k)/(a+(phy_S(:,k))'*phy_S(:,k));
    theta_Cap1(:,k) = theta_Cap1(:,k-1)+K(:,k)*(y(k)-(phy_S(:,k))'*theta_Cap1(:,k-1));
    y3(k) = -theta_Cap1(1,k-1)*y(k-1)-theta_Cap1(2,k-1)*y(k-2)+theta_Cap1(3,k-2);
end
plot(y3,'y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%4:最均方差算法（LMS，Least Mean Square算法）
theta_Cap2 = zeros(3,N_Count);%初始化θcap（t）1矩阵3*1维
gama_LMS=0.15;
for k=3:N_Count
    theta_Cap2(:,k) = theta_Cap2(:,k-1)+gama_LMS*phy_S(:,k)*(y(k)-(phy_S(:,k))'*theta_Cap2(:,k-1));
    y4(k) = -theta_Cap2(1,k-1)*y(k-1)-theta_Cap2(2,k-1)*y(k-2)+theta_Cap2(3,k-2);
end
%plot(y4,'m');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘制出参数的估计实时更新过程；
theta_1 = -1.5*ones(N_Count,1);
theta_2 = 0.7*ones(N_Count,1);
theta_3 = 0.3*ones(N_Count,1);
figure(2);hold on;
subplot(3,1,1),plot(theta_1,'--');hold on;plot(theta_Cap(1,:),'r');
title('参数的实时估计曲线—迭代最小二乘估计(RLS)');
axis([0,N_Count,-2,1]);xlabel('\theta(1)的估计值');
subplot(3,1,2),plot(theta_2,'--');hold on;plot(theta_Cap(2,:),'r');
xlabel('\theta(2)的估计值');
subplot(3,1,3),plot(theta_3,'--');hold on;plot(theta_Cap(3,:),'r');%axis([0,N_Count,-1,1]);
xlabel('\theta(3)的估计值');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘制出参数的估计实时更新过程；
theta_1 = -1.5*ones(N_Count,1);
theta_2 = 0.7*ones(N_Count,1);
theta_3 = 0.3*ones(N_Count,1);
figure(3);hold on;
subplot(3,1,1),plot(theta_1,'--');hold on;plot(theta_Cap1(1,:),'r');
title('参数的实时估计曲线-Kaczmarz’s 投影算法');
axis([0,N_Count,-2,1]);xlabel('\theta(1)的估计值');
subplot(3,1,2),plot(theta_2,'--');hold on;plot(theta_Cap1(2,:),'r');
xlabel('\theta(2)的估计值');
subplot(3,1,3),plot(theta_3,'--');hold on;plot(theta_Cap1(3,:),'r');%axis([0,N_Count,-1,1]);
xlabel('\theta(3)的估计值');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘制出参数的估计实时更新过程；
theta_1 = -1.5*ones(N_Count,1);
theta_2 = 0.7*ones(N_Count,1);
theta_3 = 0.3*ones(N_Count,1);
figure(4);hold on;
subplot(3,1,1),plot(theta_1,'--');hold on;plot(theta_Cap2(1,:),'r');
title('参数的实时估计曲线-最小均方差算法');
axis([0,N_Count,-2,1]);xlabel('\theta(1)的估计值');
subplot(3,1,2),plot(theta_2,'--');hold on;plot(theta_Cap2(2,:),'r');
xlabel('\theta(2)的估计值');
subplot(3,1,3),plot(theta_3,'--');hold on;plot(theta_Cap2(3,:),'r');%axis([0,N_Count,-1,1]);
xlabel('\theta(3)的估计值');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%