递推最小二乘法算法

题目：
（递推最小二乘法）
考虑如下系统：
)()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-=-+--
式中，)(k ξ为方差为0.1的白噪声。

取初值I P 610)0(=、00=∧
）（θ。

选择方差为1的白噪声作为输入信号)(k u ，采用PLS 法进行参数估计。

Matlab 代码如下：
clear all
close all
L=400; %仿真长度
uk=zeros(4,1); %输入初值：uk(i)表示u(k -i)
yk=zeros(2,1); %输出初值
u=randn(L,1); %输入采用白噪声序列
xi=sqrt(0.1)*randn(L,1); %方差为0.1的白噪声序列
theta=[-1.5;0.7;1.0;0.5]; %对象参数真值
thetae_1=zeros(4,1); %（）θ初值
P=10^6*eye(4); %题目要求的初值
for k=1:L
phi=[-yk;uk(3:4)]; %400×4矩阵phi 第k 行对应的y(k-1),y(k-2),u(k-3), u(k-4) 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(4)-K*phi')*P;
%更新数据
thetae_1=thetae(:,k);
for i=4:-1:2
uk(i)=uk(i -1);
end
uk(1)=u(k);
for i=2:-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 2]);
结果
分析如下
系统方程为
)()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-=-+--
由CAR 模型：
)()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-+---=
可以得到：
2117.05.11)(----+=z z z A
)5.01()(131---+=z z z B
我们可以知道纯迟延为3T ，na=2，nb=1，d=3，
phi(k,:)=[-yk;uk(3:4)]'; 400×4矩阵phi第k行对应的y(k-1),y(k-2),u(k-3), u(k-4)
y(k)=phi*theta+xi(k);输出等于矩阵phi与真值矩阵相乘加上白噪声。

在这里phi是列向量不需要用phi(k,:)形式
这是循环的主体，接下来
一开始我们赋给初值给输入以及输出
初始输入u(0)=u(-1)=u(-2)=u(-3)=0；初始输出y(0)=y(-1)=0；
所以
=u
-
+
u
+
y
y
y
+
)2-(
5.0
)3-(
)1-(
7.0
)1(
5.1
)0(
)1(ξ
经过白噪声的更新代码后，以及输入输出值得更新后我们有，也就是把计算出来的y(1)给到下一组y(2)的计算中区去：
-
+
+
=u
y
+
y
y
u
)2(
)1-(
5.0
)2-(
)2(ξ
)0(
5.1
)1(
7.0
于是我让程序中显示phi以便只观的表现递推最小二乘法和批处理最小二乘法的不同：
>> phi
phi =
-1.9333
-3.3952
-0.4448
-0.4710
递推最小二乘法phi截图如下：
相比于第一次作业的批处理最小二乘法的phi：
可以发现递推最小二乘法和批处理最小二乘法的运算结果在L数据长度是一样的情况下是相同的，两者的实质原理相同，仅仅是数学原理上的差别。

但是优点在于后一种办法，取一组数据用一组，避免了矩阵求逆的过程，减轻了计算机软件的负担，因为在matlab中求逆是应该避免的。