系统辨识各类最小二乘法汇总.doc

最小二乘方法的总结与比较
Y30150638 控制二班 王思明（作业二）
最小二乘法大约是1795年高斯在他那著名的星体运动轨迹研究工作中提出来的。后来，最小二乘法就成了理论的奠基石。由于最小二乘原理简单，编程程序也不困难，所以它颇受人们重视，应用相当广泛。最小二乘的基本结果有两种形式，一种是经典的一次完成算法，另一种是现代的递推算法。本次作业只讨论递推算法和其他扩展的最小二乘方法比较。

一、 最小二乘法
最小二乘法具有两个缺陷：
1、当噪声模型是有色噪声时，最小二乘参数估计不是你无偏、一致估计。
2、随着数据的增长，最小二乘法出现数据饱和现象，以致算法慢慢失去修正能力。
下面是其分别在白噪声和有色噪声时的参数估计曲线。最小二乘的程序见附录。

图1-1 RLS白噪声时参数估计变化曲线

图1-2 RLS有色噪声时参数估计变化曲线

从图中可以很明显的看出，当噪声变成有色噪声时，RLS算法是有偏的，始终不能得出参数的真值，而且可以看出随着数据的增多，算法对参数的修正能力明显下降，到最后参数值几乎不变，这说明RLS算法进入了数据饱和状态。

二、 遗忘因子法和限定记忆法
这两种算法（RFF,RFM）是基于RLS的数据饱和现象而提出来的一种辨识算法。
RFF很简单，在老数据加上遗忘因子以降低老数据提供的信息量，增加新数据的信息量。仅仅在RLS算法内部加入遗忘因子系数，这里因为篇幅有限，不给出其程序以及最后的参数变化曲线。
RFM就比较特殊，其限定有效数据的长度L，即始终只有长度L的数据在影响参数估计，可以从根本上解决数据饱和现象。下面是其实验的参数估计变化曲线。

图2-1 RFM白噪声时的参数估计曲线
从图中可以看出，相比于RLS的白噪声时的参数估计，变化曲线围绕真值震荡变化，这说明没有数据饱和现象，辨识结果良好。

三、 偏差补偿法
当噪声是有色噪声时，我们提出偏差补偿最小二乘法，因为从图1-2我们可以看出，最后的参数估计是有偏的，所以我们在算法中每次都减去这个偏差量，使得最后的结果能够逼近真值。我们称这种算法为RCLS。
但是，既然我们加入了偏差量，如果我们的噪声是白噪声时，那这个偏差量就变成了误差量，使得参数估计轨迹偏真值。
下图是其参数估计曲线，程序见附录。

图3-1 RCLS有色噪声时的参数变化曲线

图3-2 RCLS白噪声时的参数估计曲线
从图中可以很清楚的看出，RCLS可以对有色噪声情况下，对参数进行无偏一致估计，但是在白噪声的情况下，却不能得到无偏估计，其原因上面也提到了，这里就不在赘述了。
四、 增广最小二乘和广义最小二乘
针对不同的有色噪声模型，我们用不同的辨识方法，并能顺带求出噪声的模型模型，通过扩充参数向量和数据向量的维数，把噪声模型的辨识同时考虑进去。
当噪声模型： ek=DZ-1*vk ( v(k) 为白噪声 )时，我们采用增广最小二乘方法。其实验结果如下，其程序见附录。

图4-1 RELS参数估计曲线
当噪声模型： ek=1C（z-1）*vk ( v(k) 为白噪声 )时，我们采用广义最小二乘方法。其实验结果如下，其程序见附录。

图4-2 RGLS参数估计曲线
从上图可以看出，针对不同的参数模型，RELS和RGLS分别能发挥不同的作用，并且给出参数的无偏一致估计，并顺带求出了噪声参数，其辨识结果良好，符合要求.
五、 辅助变量法
在有些问题中，我们并不能事先知道参数的模型，而且我们也不需要求出噪声参数，这时普通的RLS算法已经不能满足我们的要求，我们提出了辅助变量法，RIV。下图给出其实验结果，具体程序见附录。

图5-1 RIV有色噪声时参数估计变化曲线
从图中可以看出，RIV算法对有色噪声的情况，也能给出参数的无偏估计，其结果与RLS（图1-2）比较可以看出区别。
六、 总结
RLS：
当噪声模型是有色噪声时，最小二乘参数估计不是你无偏、一致估计。而且随着数据的增长，最小二乘法出现数据饱和现象，以致算法慢慢失去修正能力。
RFF.RFM：
RFF很简单，在老数据加上遗忘因子以降低老数据提供的信息量，增加新数据的信息量，RFM就比较特殊，其限定有效数据的长度L，即始终只有长度L的数据在影响参数估计，可以从根本上解决数据饱和现象
RCLS：
当噪声是有色噪声时，我们提出偏差补偿最小二乘法，在算法中每次都减去这个偏差量，使得最后的结果能够逼近真值。但却不适合噪声为白噪声的时候。
RELS：
当噪声模型： ek=DZ-1*vk ( v(k) 为白噪声 )时，我们采用增广最小二乘方法。能辨识出参数（包括噪声参数）的无偏估计。
RGLS：
当噪声模型： ek=1C（z-1）*vk ( v(k) 为白噪声 )时，我们采用广义最小二乘方法。能辨识出参数（包括噪声参数）的无偏估计。
RIV：
当不确定参数模型时，我们可以直接采用RIV模型得到参数的无偏一致估计。

七、 附录
1、 RLS实验程序：
clc
clear

% y1=randn(1,2500); %产生高斯白噪声
% y1=y1/std(y1);
% y1=y1-mean(y1);
% a=0;
% b=sqrt(0.1);
% y1=a+b*y1;
% figure(1);
% plot(y1);
% title('高斯白噪声');

for i=1:2500 %产生有色噪声
y1(i)=1;
end

N=400;
uk=rand(1,N);
for i=1:N %产生随机输入
uk(i)=uk(i)*(-1)^(i-1);
end
figure(2);
plot(uk);
title('随机输入曲线');

yk=zeros(1,N);
for k=3:N %定义输出
yk(k)=1.5*yk(k-1)-0.7*yk(k-2)+uk(k-1)+0.5*uk(k-2)+y1(k);
end
figure(3);
plot(yk);
title('对应输出曲线');

theta=[0;0;0;0];
p=10^6*eye(4);

for t=3:N %迭代计算 （ls算法）
h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2)]);
x=1+h'*p*h;
p=(p-p*h*1/x*h'*p);
theta=theta+p*h*(yk(t)-h'*theta);

a1t(t)=theta(1);
a2t(t)=theta(2);
b1t(t)=theta(3);
b2t(t)=theta(4);

end

theta
t=1:N; %绘制参数变化曲线
figure(4);
plot(t,a1t(t),t,-1.5,t,a2t(t),t,0.7,t,b1t(t),t,1,t,b2t(t),t,0.5);
title('参数变化曲线');
text(300,-1.5,'a1');
text(300,0.7,'a2');
text(300,1,'b1');
text(300,0.5,'b2');
2、 RFM实验程序（部分）
只给出算法更改部分：
for s=3:100;
h=([yk(s-1);yk(s-2);uk(s-1);uk(s-2)]);
p=p-p*h*(inv(1+h'*p*h))*h'*p;
theta=theta+p*h*(yk(s)-h'*theta);
end
s
for t=3:N-100
m=s+t;
h=([yk(m-1);yk(m-2);uk(m-1);uk(m-2)]);
p=p-p*h*(inv(1+h'*p*h))*h'p;
theta=theta+p*h*(yk(m)-h'*theta);
h=([yk(t-1);yk(t-2);uk(t-1);uk(t-2)]);
p=p+p*h*(inv(1-h'*p*h))*h'*p;
theta=theta+p*h*(yk(t)-h'*theta);

a1t(t)=theta(1);
a2t(t)=theta(2);
b1t(t)=theta(3);
b2t(t)=theta(4);

end
3、 RCLS实验程序（部分）
theta=[0;0;0;0];
p=10^6*eye(4);
j=0;
thetc=[0;0;0;0];

for t=3:N
h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2)]);
x=1+h'*p*h;
j=j+((yk(t)-h'*theta)^2)/(1+h'*p*h);
p=(p-p*h*1/x*h'*p);
theta=theta+p*h*(yk(t)-h'*theta)
d=[1 0 0 0;0 1 0 0;0 0 0 0; 0 0 0 0];
w=j/(t*(1+thetc'*d*theta));
thetc=theta-t*w*p*d*thetc;

a1t(t)=thetc(1);
a2t(t)=thetc(2);
b1t(t)=thetc(3);
b2t(t)=thetc(4);

end
4、 RELS实验程序（部分）
for t=3:N
h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2);y1(t-1);y1(t-2)]);
x=1+h'*p*h;
p=(p-p*h*1/x*h'*p);
theta=theta+p*h*(yk(t)-h'*theta);

a1t(t)=theta(1);
a2t(t)=theta(2);
b1t(t)=theta(3);
b2t(t)=theta(4);
d1t(t)=theta(5);
d2t(t)=theta(6);

end
5、 RGLS试验程序（部分）
for t=3:N

he=([-e(t-1);-e(t-2)]);
xe=1+he'*pe*he;
pe=(pe-pe*he*1/xe*he'*pe);
thete=thete+pe*he*(e(t)-he'*thete);

c1t(t)=thete(1);
c2t(t)=thete(2);

end

for t=3:N
ykf(t)=yk(t)+c1t(t)*yk(t-1)+c2t(t)*yk(t-2);
ukf(t)=uk(t)+c1t(t)*uk(t-1)+c2t(t)*uk(t-2);

h=([-ykf(t-1);-ykf(t-2);ukf(t-1);ukf(t-2)]);
x=1+h'*p*h;
p=(p-p*h*1/x*h'*p);
theta=theta+p*h*(ykf(t)-h'*theta);

a1t(t)=theta(1);
a2t(t)=theta(2);
b1t(t)=theta(3);
b2t(t)=theta(4);

end
6、 RIV试验程序（部分）
for t=5:N
hx=([-yk(t-3);-yk(t-4);uk(t-1);uk(t-2)]);
h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2)]);
x=1+h'*p*hx;
p=(p-p*hx*1/x*h'*p);
theta=theta+p*hx*(yk(t)-h'*theta);

a1t(t)=theta(1);
a2t(t)=theta(2);
b1t(t)=theta(3);
b2t(t)=theta(4);

end
15