最小二乘法在系统辨识中的应用包含相关的三种算法

《最小二乘法在系统辨识中的应用包含相关的三种算法》由会员分享，可在线阅读，更多相关《最小二乘法在系统辨识中的应用包含相关的三种算法（11页珍藏版）》请在装配图网上搜索。

1、2008级硕士研究生系统建模理论试卷已知一个三阶线性离散系统的输入、输出数据，共有40个采样值，试分别 用：最小二乘法(LS)、递推最小二乘法(RLS)、广义最小二乘法(GLS)进 行参数估计，给出相应的数学模型，并阐述相应的辨识原理。5.9004 3.98702.2486 0.9525 05325 1.5227 0.42001.0786k21222324252627282930u(k) 1.0576 1.0071 1 13420.07400.67590.5221099540.5271 1.76560.4936y(k) T55790.6640 1.42222.64442.95723.6340-

2、3.128138334 3.25421.1568k31323334353637383940u(k) 1.48100.9591 3.12930.36040.42510.41850.6728 0.00272.11451.1157y(k) 0.06150.9120 0.06923.27313.74864.31944.72305.27815.1507 2.7235y(k) 9.0552 8.1735k12345678910u(k) 0.82510.09880.46280.91682.23250.07772.36540.34761.1473-1.9035y(k) 1.5333-1.06801.06660

3、.5284 -0.58353.1471 3.71856.21496.30267.2705k11121314151617181920u(k) 0.92291.6400 0.84100.7599 0.4739 0.1784 1.7760 1.67221.2959 0.0591一、最小二乘法(LS)1、数学模型设时不变SISO动态过程的数学模型为A(z-1)z(k)二 B(zt)u(k) + n(k)(1.1.1)其中，u(k)为过程的输入量，z(k)为过程的输出量,n(k)是噪声，多项式A(z-1)和B(z-1)为：A( z-1) = 1 + a z-1 + a z-2 + b a z -na丨

4、B (z-1) = b z-1 + b z-2 + + b z-nbl12nh在本题中，n = n =3.即a bf A( z-1)二 1 + a z-1 + a z -2 + a z -3123丨 B(z-1)二 bz-1 + b z-2 + b z-3J123将此模型写成最小二乘格式z(k)二he(k) 6+ n(k)(1.1.2) 其中，z(k)是过程的输出量；he(k)是可观测的数据向量；(k)是均值为零的随机噪声。h(k) = z (k 1),z (k 2),-z (k 3), u (k 1), u (k 2), u (k 3)卞0 二a , a , a , b , b , b 卡1

5、23123对于k = 1,2,L,方程式(1.1.2)构成一个线性方程组，可以把它写成z = z (1), z (2),z (l )t ihT (1)-z (0)l00u(0)00 _hT=zz (0)0u(1)u(0)0_hT(l)_ z (l 1)z (l 2)z (l 3)u(l 1)u(l 2)u(l 3)n = n(1), n (2), n(l )tVh =i(1.1.3)利用数据序列如和h(k)，极小化准则函数Jz(k) hT(k)02(1.1.4)k=1使J(0) = min的0估计值记作0，称为参数0的最小二乘估计值。通过极小化(1.1.4)式来计算0的方法称作最小二乘法，未知

6、模型参数0最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得 到的，这种模型输出能最好地接近实际过程的输出。2、辨识原理考虑模型(1.1.2)式的辨识问题，其中z(k)和h(k)都是可观测的数据，0是待估计参数，准则函数取(1.1.4)根据(1.1.3)的定义，准则函数J(0)可写成二次型的形式J =(z H 9 )T (z H 0)(1.2.1)IlII显然上式中的H/0代表模型的输出，或者说是过程的输出预报值。因此J(0)可以看作来衡量模型输出与实际过程输出的接近情况。极小化J(0)，求得参数0的估计值0 = (HtH )-1 Htz(1.2.2)l ll l将使模型的输出最

7、好的预报过程的输出。3、辨识结果0 = a a a b b b = -1.0000 0.0000 0.0000 0.0000 0.0000 0.0000123123二、递推最小二乘法(RLS)1、数学模型在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但是由于具体使用时占用内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。 递推算法的基本思想可以概括成：新的估计值9 (k)=老的估计值0 (k -1)+修正项(2.1.1)在此算法中，时不变SISO动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格 式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的

8、辨识原理中描述。2、辨识原理首先将1.2.2式一次完成算法写成(2.2.1)(2.2.2)0 二(HtH )-iHtz 2 P(/)Ht z 二工h(i)5 (i)-i工h(i)z(i)l ll ll li=1i=1定义P-1(k)=工 h(i)b (i) = Ht Hk kP-1 (k -1)=丈 h(i)hT (i) = Ht Hk-1k-1hT (1)- hT (2)，H k-1 =-hT (1)hT (2)hT (k)k (k -1)i=1其中：Hk由(2.2.2)式可得:P-1 (k) = t1 h(i)hT (i) + h(k)hT (k) = P-1 (k -1) + h(k血(

9、k) (2.2.3)i =1令: z = z(1), z(2),z(k - 1)tk-1则：0(k 1) = (Ht H)-1 Ht z = P(k 1)刃 h(i)z(i)k-1k -1k-1 k -1i=1于是有P-1 (k - l)0(k -1)=乞h(i)z(i)(2.2.4)i=1令 z = z (1), z (2),z (k )t k利用(2.2.3)和(2.2.4)式，可得0(k) = (HtH )-1 Htz = P(k)h h(i)z(i)k kk ki=1=P(k)P-1 (k - 1)0(k -1) + h(k)z(k)(2.2.5)=P(k )P-1 (k) - h(k

10、 )hT (k )0(k -1) + h(k) z (k)=0(k -1) + p(k )h(k) z (k) - hT (k )0(k -1)引进增益矩阵K (k)，定义为K(k)二 P(k)h(k)(2.2.6)则(2.2.5)式写成0(k) = 0(k -1) + K(k)z(k) -d (k)0(k -1)(2.2.7)进一步把(2.2.3)式写成P(k)二P-i(k -1) + h(k)hT(k)-i(2.2.8)为了避免矩阵求逆运算，利用矩阵反演公式可将(2.2.8)式演变成P(k)二P-1 (k -1) + h(k)k (k)-i 二I - K(k)d (k)P(k -1)(2.

11、2.9)将(2.2.9)式代入(2.2.6)式，整理后有K(k)二 P(k - 1)h(k)d (k)P(k - 1)h(k) +1-1(2.2.10)综合(2.2.7)、(2.2.9)、(2.2.10)式便得到最小二乘参数估计递推算法。0(k) = 0(k -1) + K(k)z(k) - d (k)0(k -1)P(k)二P-1 (k -1) + h(k)k (k)-1 二I - K(k)d (k)P(k -1)K(k)二 P(k - 1)h(k)d (k)P(k - 1)h(k) +1-1利用上述公式即可求得参数0的估计值03、辨识结果0 = a a a b b b = -1.0005

12、- 0.0007 - 0.0002 0.0000 - 0.0005 -0.0001123123三、广义最小二乘法(GLS)1、数学模型设时不变SISO动态过程的数学模型为1A(z -1)z(k)二 B(z -1)u(k) +_- v(k)C (z -1)(3.1.1)其中，u(k)和z(k)分别表示过程的输入和输出，v(k)是均值为零的不相关随机噪声，多项式 A(z-1)、B(z-1)和 C(z-1)为：A(z-1) = 1 + a z-1 + a z -2 + a z -na12n B(z-1) = b z-1 + b z-2 + b z-nb12nbC (z -1) = 1 + c z -

13、1 + c z -2 + c z -nc12nc在本题中，在本题中，n = n = n =3.即a b cA( z -i) = 1 + a z-i + a z -2 + a z -3123 B(z-i) = b z-i + b z-2 + b z-31123(3.1.2)C (z -1) = 1 + c z -1 + c z -2 + c z -3J123广义最小二乘发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法 对滤波后的数据进行辨识。2、辨识原理(k)二 C(z -1)z(k)(k)二 C(z-Qu(k)由(3.1.1)式可得 A(zt)C(z-1)z(k)二 B(z

14、t)C(z-1)u(k) + v(k)(3.2.1)(3.2.2)(3.2.3)(k) = z (k 1),-z (k 2),-z (k 3), u (k 1), u (k 2), u (k 3”ffffff0 二a , a , a , b , b , b 卡123123(3.2.4)可将模型(3.1.1)式化成最小二乘格式zf (k)= hf T (k)0 + v(k)由于上式v(k)是白噪声，所以利用最小二乘法即可获得参数0的无偏估计。但是，数据向量h (k)中的变量均需要按照(3.2.2)式计算，然而噪声模型C(z-1)并不知道。为此需要 f用迭代的方法来估计C(z-1)。令 e(k)=

15、1C (z -1)(3.2.5)(326)fh (k) = -e(k - 1),-e(k - 2),-e(k - 3片 置 e0 (k) = c , c , c tJe123就把噪声模型(3.2.5)也化成最小二乘格式e(k)=血(k)0 + v(k)ee由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数0的无偏估计。e但是，数据向量h (k)依然包含着不可测得噪声量e(k -1),e(k - 2),e(k - 3)，它可用相应e的估计值代替，置h (k)= -e(k - 1),-e(k - 2),-e(k -3”(3.2.7)e其中 e(k)= 0, k 0 时，按照式子 e

16、(k)= z (k)- h (k)o (k)计算，式子中 h(k) = z (k 1),z (k 2),z (k 3), u (k 1), u (k 2), u (k 3)F 综上分析，广义最小二乘递推算法可归纳成：八0(k)二 e(k -1) + K (k)z (k) - hT (k)0(k -1)f ffP (k) = I - K (k)hT (k)P (k -1)ff f fK (k) = P (k - 1)h (k)hT (k)P (k - 1)h (k) + 1-ifff f ff八e (k) = e (k -1) + K (k)e(k) -hT(k)0 (k -1)eeeeeP (

17、k) = I K (k )hT (k )P (k 1)eeeeK (k) = P (k 1)h (k)hT (k)P (k 1)h (k) +1-1eeeeee利用上述公式即可求得参数e的估计值0并能求出噪声模型的估计。3、估计结果0 = aa a b b b = -1.0005 -0.0008 0.0002 -0.0000 -0.0005 -0.0001123123噪声模型的估计值为 1.0e-003 * -0.00000.41150.0001 即 C(z-1)沁 1附录：MATLAB程序1、一次完成的最小二乘法%次完成的最小二乘法clears=(l:40);z=1.5333 -1.0680

18、1.0666-0.5284 -0.58353.1471 -3.71856.2149 -6.30267.2705. .-9.05528.1735 -5.90043.9870 -2.24860.9525 -0.5325 -1.52270.42001.0786.-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568.0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235; u=0.82510.09880.4628

19、-0.91682.23250.07772.3654 0.34761.1473-1.9035.-0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591.-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936.1.4810 0.9591 -3.1293-0.3604 -0.42510.4185 -0.6728 -0.00272.11451.1157;h=zeros(40,6);h(1,:) = -z(1) 0 0

20、 u(1) 0 0;h(2,:) = -z(2) -z(1) 0 u(2) u(1) 0;for i=3:1:40h(i,:) = -z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2);endo=inv(h*h)*h*z%误差分析J=0;for i=1:1:40e(i) = z(i)-h(i,:)*o;J=J+e(i)2;endJplot( s,e);2、递推的最小二乘法%递推的最小二乘法clears=(1:40);z=1.5333 -1.06801.0666-0.5284 -0.58353.1471 -3.71856.2149 -6.30267.2705. .-9

21、.05528.1735 -5.90043.9870 -2.24860.9525 -0.5325 -1.52270.42001.0786.-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568.0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235; u=0.82510.09880.4628 -0.91682.23250.07772.3654 0.34761.1473-19035-0.9229 1.6400 -

22、0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591.-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936.1.4810 0.9591 -3.1293-0.3604 -0.42510.4185 -0.6728 -0.00272.11451.1157;h=zeros(40,6);h(1,:) = -z(1) 0 0 u(1) 0 0;h(2,:) = -z(2) -z(1) 0 u(2) u(1) 0;for i=3:1:40h(i,:

23、) = -z(i) -z(i-1) -z(i-2) u(i) u(i-1) u(i-2);endp0=106;o0=0.001;k(1,:) = p0*h(1,:)*inv(h(1,:)*p0*h(1,:)+1); K=zeros(6,40);K(:,1)=k(1,:);p=zeros(6,6,40);p(:,:,1) = eye(6)-K(:,1)*h(1,:)*p0;o=zeros(6,40);o(1,:)=o0+K(:,1)*z(1)-h(1)*o0;for i=2:1:40k(i,:) = p(:,:,iT)*h(i,:)*inv(h(i,:)*p(:,:,i-1)*h(i,:)+1)

24、;p(:,:,i) = eye(6)-k(i,:)*h(i,:)*p(:,:,i-1);o(i,:) = o(i-1,:)+k(i,:)*z(i)-h(i,:)*o(i-1,:);end%辨识结果o(40,:)%误差分析J=0;for i=2:1:40e(i) = z(i)-h(i,:)*o(i-1,:);J=J+e(i)2;endplot( s,e);J3、广义最小二乘法%广义预测的最小二乘法clears=(1:40);z=1.5333 -1.06801.0666-0.5284 -0.58353.1471 -3.71856.2149 -6.30267.2705. .-9.05528.1735

25、 -5.90043.9870 -2.24860.9525 -0.5325 -1.52270.42001.0786.-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568.0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235;u=0.82510.09880.4628 -0.91682.23250.07772.3654 0.34761.1473-1.9035.-0.9229 1.6400 -0.8410 0.7

26、599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591.-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936.1.4810 0.9591 -3.1293-0.3604 -0.42510.4185 -0.6728 -0.00272.11451.1157;zf=z;uf=u;hf=zeros(40,6);hf(1,:) = -zf(1) 0 0 uf(1) 0 0;hf(2,:) = -zf(2) -zf(1) 0 uf(2) uf(1) 0;o0=0.001

27、;pf0=106;oe0=0;pe0=1;%k=1 时kf(1,:) = pf0*hf(1,:)*inv(hf(1,:)*pf0*hf(1,:)+1);pf=zeros(6,6,40);pf(:,:,1) = eye(6)-kf(1,:)*hf(1,:)*pf0;of=zeros(6,40);of(1,:)=o0+kf(1,:)*z(1)-hf(1)*o0;ke=zeros(40,3);e(1)=z(1)-hf(1,:)*of(1,:);he(1,:) = 0 0 0;ke(1,:) = pe0*he(1,:)*inv(he(1,:)*pe0*he(1,:)+1);pe=zeros(3,3,4

28、0);pe(:,:,1) = eye(3)-ke(1,:)*he(1,:)*pe0;oe=zeros(3,40);oe(1,:)=oe0+ke(1,:)*e(1)-he(1)*oe0;%k=2 时kf(2, :) = pf(:,:,2-1)*hf(2,:)*inv(hf(2,:)*pf(:,:,2T)*hf(2,:)+1);pf(:,:,2) = eye(6)-kf(2,:)*hf(2,:)*pf(:,:,2-1);of(2,:) = of(2-1,:)+kf(2,:)*z(2)-hf(2,:)*of(2-1,:); e(2)=z(2)-hf(2,:)*of(2,:);he(2,:) = -e

29、(1) 0 0;ke(2, :) = pe(:,:,2-1)*he(2,:)*inv(he(2,:)*pe(:,:,2-1)*he(2,:)+1);pe(:,:,2) = eye(3)-ke(2,:)*he(2,:)*pe(:,:,2-1);oe(2,:) = oe(2-1,:)+ke(2,:)*e(2)-he(2,:)*oe(2T,:);%k=3 时kf(3,:) = pf(:,:,3-1)*hf(3,:)*inv(hf(3,:)*pf(:,:,3-1)*hf(3,:)+1);pf(:,:,3) = eye(6)-kf(3,:)*hf(3,:)*pf(:,:,3-1);of(3,:) = o

30、f(3-1,:)+kf(3,:)*z(3)-hf(3,:)*of(3-1,:); e(3)=z(3)-hf(3,:)*of(3,:);he(3,:) = -e(2) -e(1) 0;ke(3,:) = pe(:,:,3T)*he(3,:)*inv(he(3,:)*pe(:,:,3T)*he(3,:)+1);pe(:,:,3) = eye(3)-ke(3,:)*he(3,:)*pe(:,:,3T); oe(3,:) = oe(3T,:)+ke(3,:)*e(3)-he(3,:)*oe(3T,:);%k4 时for i=4:1:40 zf(i)=z(i)+z(iT)*oe(iT,l)+z(i-2)

31、*oe(iT,2)+z(i-3)*oe(iT,3); uf(i)=u(i)+u(iT)*oe(iT,l)+u(i-2)*oe(iT,2)+u(i-3)*oe(iT,3); hf(i,:) = -zf(i) -zf(i-1) -zf(i-2) uf(i) uf(iT) uf(i-2); kf(i,:) = pf(:,:,il)*hf(i,:)*inv(hf(i,:)*pf(:,:,il)*hf(i,:)+l) pf(:,:,i) = eye(6)-kf(i,:)*hf(i,:)*pf(:,:,i-l);of(i,:) = of(iT,:)+kf(i,:)*z(i)hf(i,:)*of(iT,:); e(i)=z(i)hf(i,:)*of(i,:);he(i,:) = -e(i-l) -e(i-2) -e(i-3);ke(i,:)二pe(:,:,i_1)*he(i,:)，*inv(he(i,:)*pe(:,:,i_1)*he(i,:)，+)， pe(:,:,i) = eye _ke(i,:)*he(i,:)*pe(:,:,i_l);oe(i,:) = oe(i_1,:)+ke(i,:)*e(i)_he(i,:)*oe(i_1,:);end%辨识结果of(40,:)oe(40,:)%误差分析plot( s,e);