加权递推最小二乘法仿真

6
2
4
1.5
Y
2
Y
0 5 X 10
1
0
0.5
-2
Hale Waihona Puke Baidu
0
0
5 X
10
图 2 加权递推最小二乘算法
表 1 两种算法的误差对比 X 1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10
e1（一般） 0.886 3 e2（加权 递推） 2.541 6
0.857 5 0.030 4
四、仿真程序
clc;clear; x=[1.01,2.03,3.02,4.01,5,6.02,7.03,8.04,9.03,10]; y=[9.6,4.1,1.3,0.4,0.05,0.1,0.7,1.8,3.8,9.0]; l=ones(1,10); h=[x.^2;x;l]; z=y; %一般最小二乘法 theta1=polyfit(x,y,2) q1=polyval(theta1,x); e1=abs(q1-y) %拟合值与实际值之差的绝对值 figure(1); subplot(1,2,1);
三、加权递推最小二乘算法仿真结果和分析
这里使用 MATLAB 软件进行仿真。首先采用了一般最小二乘算法进行仿真， 仿真结果如图 1 所示。其次采用了加权递推最小二乘算法进行仿真，权阵值由
1/ y Y ，初始值的确定采用了上述的第一种方法，仿真结果如图 2 所示。表
1 是两种算法的误差值的比较，表 2 是两种算法的辨识结果。
所以有P 1 k = i hT i h i k hT k h k
i 1
k 1
P
1
k-1 k hT k h k
k 1 ˆ k 1 ( H T H ) 1 H T Z =P k 1 i h i z i 由 k 1 k 1 k 1 k 1 k 1 k 1 i 1
9.0
首先写出系统的最小二乘表达式，分析 Y 与 X 的关系式 Y aX 2 bX c ， 与式 z (k ) hT (k ) e(k ) 对比可得，取 hT (k ) ( x2
x 1) ， (a b c) ，采用
加权递推最小二乘辨识对参数 a、b 与 c 进行估计。给出加权递推最小二乘算法 和仿真程序及仿真结果。
1 T 可得P k = P k-1 k h k h k T 1 1
1 = P k 1 P k 1 h k h k P k 1 hT k P k 1 h k k P k 1 h k hT k = I T P k 1 1 h k P k 1 k 1 所以K k P k 1 h k hT k P k 1 h k k 综上所述，加权递推最小二乘算法（RWLS）的递推过程如下： K k P k 1 h k hT k P k 1 h k 1 k 1 T P k = I K k h k P k 1 ˆ k ˆ k 1 K k z k hT k ˆ k 1 ˆ k ˆ k 1 递推过程结束条件为 max , 是适当小的数。 ˆ k 1 初始值的确定方法分为两种：1.利用一次完成算法预先求出
ˆ k 1 = i h i z i 可得P 1 k 1
i 1 k ˆ k =( H T H ) 1 H T Z =P k i h i z i 进一步可得 k k k k k k i 1 1 ˆ P k P k 1 k 1 k h k z k 1 T ˆ P k P k k h k h k k 1 k h k z k ˆ k 1 P k h k k z k hT k ˆ k 1
ˆ) (Z H ˆT ˆ 加权最小二乘算法。一般最小二乘算法可由 J ( L L ) ( Z L H L ) min ˆ) (Z H ˆT ˆ 表示，加权最小二乘算法可由 J ( L L ) L ( Z L H L ) min 表示，其中
L diag[(1), (2), , (L)] 为权阵，则在加权最小二乘算法下的 的估计值为
ˆ ( H T H )1 H T Z 。 WLS L L L L L L
2.加权递推最小二乘法介绍 一般最小二乘或加权最小二乘为一次完成算法或批处理算法。 但是计算量大、 存储大、不适合在线辨识。因此采用一种参数递推估计——递推最小二乘算法。 在递推的过程中加入权阵，就成为加权递推最小二乘法。 其基本思想是当前估计
二、加权递推最小二乘算法介绍
1.加权最小二乘算法简要介绍 一般最小二乘估计精度不高的原因之一是对测量数据同等对待， 但是实际情 况是各次测量数据很难在相同的条件下获得。有的测量值置信度高，有的测量值 置信度低。因此提出对不同置信度的测量值采用加权的办法分别对待。置信度高 的，权重取得大些；置信度低的，权重取的小些。这种改进的最小二乘算法成为
2
一般最小二乘拟合曲线 10 实际曲线 拟合曲线 8 1.4
一般最小二乘拟合误差曲线
1.2
1 6 0.8 4
Y
Y
0.6 0.4 0.2 0 5 X 10 0 0
2
0
-2
5 X
10
图 1 一般最小二乘算法
加权递推最小二乘拟合曲线 10 实际曲线 拟合曲线 8 2.5 加权递推最小二乘拟合误差曲线 3
结果分析：与一般最小二乘算法相比，加权递推最小二乘算法的辨识后 的误差值在 x=1.01 和 x=10 两个点明显变大， x=4.01 这个点的误差值略微变大， 其他点均有不同程度的减少。 出现这种结果的原因可能就是在这两个点的 Y 的观 测值误差较大，置信度较低，采用一般最小二乘算法辨识的误差不明显，而采用 加权递推算法后，加入了权阵，减少了这两个置信度低的点的权值，增加了其他 置信度高的点的权值。
plot(x,y,'r',x,q1);%拟合曲线与实际曲线对比 xlabel('X'); ylabel('Y'); title('一般最小二乘拟合曲线'); legend('实际曲线','拟合曲线',1); subplot(1,2,2); plot(x,e1);%误差曲线 xlabel('X'); ylabel('Y'); title('一般最小二乘拟合误差曲线'); %加权递推最小二乘法 a=1./((y-sum(y)/10).^2); w=diag(a);%定义加权矩阵 p1=inv(h*w*h');%定义 P 矩阵的初始值 theta2=inv(h*w*h')*h*w*z';%定义估计值 theta 的初始值 i=1; m=1; epsilon=0.002;%定义足够小的值 while((m>epsilon)&&(i<11)) %加权递推过程 k=p1*h(:,i)*inv(h(:,i)'*p1*h(:,i)+1/w(i,i)); p2=(eye(3,3)-k*h(:,i)')*p1; theta3=theta2+k*(z(:,i)-h(:,i)'*theta2); m=max(abs((theta3-theta2)./theta2)); i=i+1; p1=p2; theta2=theta3; end theta2 q2=polyval(theta2,x); e2=abs(q2-y) %拟合值与实际值之差的绝对值 figure(2); subplot(1,2,1); plot(x,y,'r',x,q2);%拟合曲线与原曲线对比 xlabel('X'); ylabel('Y'); title('加权递推最小二乘拟合曲线'); legend('实际曲线','拟合曲线',1); subplot(1,2,2); plot(x,e2);%误差曲线 xlabel('X'); ylabel('Y'); title('加权递推最小二乘拟合误差曲线');
一、问题描述和作业要求
设某物理量 Y 与 X 满足关系式 Y aX 2 bX c ,实验获得一批数据如下表， 试辨识模型 a，b 和 c。 X 1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10
Y
9.6
4.1
1.3
0.4
0.05
0.1
0.7
1.8
3.8
0.922 1 0.590 2
0.016 5 0.021 3
0.608 2 0.386 9
0.690 5 0.474 8
0.384 5 0.369 0
0.354 9 0.023 6
1.063 8 0.111 0
0.612 3 2.308 2
表 2 两种算法的辨识结果 a theta1(一般) theta2（加权递推） 0.4575 0.3625 b -5.0734 -4.0323 c 13.3711 10.7612
1
k 1 P k i hT i h i H kT k H k i 1 综合上述式子可得 k 1 P 1 k-1 i hT i h i H T H k 1 k 1 k 1 i 1
ˆ(k 1) ＋修正项。每获得一次新的观测数据就修正一次参 ˆ(k ) ＝上次估计值 值
数估计值，随着时间的推移，可得满意的辨识。它是由一次完成算法推导得出。 具体推导过程如下：
最小二乘算法表达式为Z L H L VL hT (1) z (1) T z (2) h (2) 其中Z L , HL T z ( L) h ( L) ˆ ( H T H ) 1 H T Z =P L H T Z 加权最小二乘估计值 WLS L L L L L L L L L L L = i hT i h i i h T i z i i 1 i 1 T h (1) z (1) 1 T z (2) L h (2) T ， 其中P L = i h i h i ，H k ，Z k i 1 T z (k) h (k) k diag 1 2 1 2 k k
T 1 P L0 =( H L0 L0 H L0 ) ，L0是数据长度 ˆ T L P L H Z 0 0 L0 L0 L0 ˆ 0 ˆL 令 P 0 P L , 0 0 1
2. 直接令
2 P 0 = I ， 是充分大的实数， 是充分小的实向量 ˆ 0 =
k 1


定义增益矩阵K k P k h k k ˆk ˆ k 1 K k z k hT k ˆ k 1 则 由P 1 k =P 1 k-1 k hT k h k