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

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

1、2008级硕士研究生系统建模理论试卷已知一个三阶线性离散系统的输入、输出数据，共有40个采样值，试分别 用：最小二乘法(LS)、递推最小二乘法(RLS)、广义最小二乘法(GLS)进 行参数估计，给出相应的数学模型，并阐述相应的辨识原理。5.9004 3.98702.2486 0.9525 0.5325 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.066

3、60.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(z-1)u(k) + n(k)(1.1.1)其中，u(k)为过程的输入量，z(k)为过程的输出量，n(k)是噪声，多项式A(z-1)和B(z-1)为：f A( z-1) = 1 + a z-1 + a z -2 HF a z -

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

5、 a , b , b , b t对于k = 1,2,L ,123123方程式(1.1.2)构成一个线性方程组，可以把它写成-hT(1)-h =hT11 IhT (i)利用数据序列-z (0)-1)使 J (0) = min 的 0来计算0的方法称值记，称为参数的最z = z (1), z ,z (l )tln = n (1), n (2), n(l )ti00-z (0)0u(0) u(13)0u(0)T (k)02(1.1.4)zh（ k） 极小化准则e未知模型参数最到的，这种模型输出能最好地接近实际过程的输2、辨识原理-2)u(l - 3)(1.1.3)计值。通过极小化（1.1.4）式）都

6、是可观测的数据次型的形式Q作来衡生实际观测值与计算值之累次误差的平方和考虑模型（1.1.2）式的辨识问题，其中准则函数取（1.1.4）根据（1.1.3）的定义，准则函数J（0）可写J (0) = (z - H 9 口 (z - H 0)i ii i(1.2.1)(1.2.2)0 = （HtH ）-1 Ht zi ii i将使模型的输出最好的预报过程的输出。3、辨识结果0 = aa a b b b = -1.0000 0.0000 0.0000 0.0000 0.0000 0.0000123123二、递推最小二乘法（RLS）1、数学模型在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但

7、是由于具体使用时占用 内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。 递推算法的基本思想可以概括成：新的估计值0 (k)=老的估计值。(k -1)+修正项(2.1.1)在此算法中，时不变SISO动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格 式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的辨识原理中描述。2、辨识原理首先将1.2.2式一次完成算法写成0 二(Ht H )-1Ht z 2 P(l)Ht z 二工h(i)d (i)卜 1工h(i)z(i)1 11 11 1(2.2.1)i=1i=1定义P -1(k)=工 h(i)hT (i) = Ht

8、Hk k5、(2.2.2)1 h(i)hT (i) = Ht HIi=1k-1k-1hT (1)hT (1)其中：Hk =hT (2)，H k-1=hT (2)hT (k)hT (k -1)由(2.2.2)式可得:P-1 (k) = t1 h(i)hT (i) + h(k)d (k) = P-1(k -1) + h(k(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 - 1)0(k -1)=艺h(i)z(i)(2.2

9、.4)i=1令 z = z (1),z (2),z (k )tk利用(2.2.3)和(2.2.4)式，可得0(k) = (HtH )-1 Htz = P(k)f 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)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) = 6 (k -1) + K (k) z (k) - h

10、i (k )0 (k -1)(2.2.7)进一步把(2.2.3)式写成P(k) = P-i (k -1) + h(k)hi(k)-i(2.2.8)为了避免矩阵求逆运算，利用矩阵反演公式可将(2.2.8)式演变成P(k) = P-i (k -1) + h(k )hi (k )-i = I - K (k )h 丫 (k )P(k -1)(2.2.9)将(2.2.9)式代入(2.2.6)式，整理后有K (k) = P(k - 1)h(k )hi (k )P(k - 1)h(k) + 1-i(2.2.10)综合(2.2.7)、(2.2.9)、(2.2.10)式便得到最小二乘参数估计递推算法。6 (k)

11、 = 6 (k -1) + K (k) z (k) - hi (k )6 (k -1)P(k)二P-1 (k -1) + h(k)hi (k)-1 二I - K(k)h工(k)P(k -1)K(k)二 P(k - 1)h(k)hi (k)P(k - 1)h(k) +1-1利用上述公式即可求得参数的估计值63、辨识结果6 = aa a b b b = -1.0005 -0.0007 - 0.0002 0.0000 - 0.0005 -0.0001123123三、广义最小二乘法(GLS)1、数学模型设时不变SISO动态过程的数学模型为1A( z -1) z (k) = B (z -1)u (k)

12、+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 Hb a z-役12n B(z-1) = b z-1 + b z-2 Hb b z -nb12nbC (z -1) = 1 + c z -1 + c z -2 Hb c z -nc12nc在本题中，在本题中，n = n = n =3.即a b cA( z t) = 1 + a z-i + a z -2 + a z -3123B (z -i) = b z -i +

13、b z -2 + b z -3123C (z -1)二 1 + c z -1 + c z -2 + c z -3 123(3.1.2)广义最小二乘发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。2、辨识原理由(3.1.1)式可得 A(z-1)C(z-1)z(k)二 B(z-1)C(z-1)u(k) + v(k)(3.2.1)(k)二 C(z -1)z (k)(k) = C (z -1)u (k)(3.2.2)(k) = -z (k -1),-z (k - 2),-z (k - 3), u (k -1), u (k - 2), u (k - 3幵fff

14、fff6 = a , a , a , b , b , b b123123(3.2.3)可将模型(3.1.1)式化成最小二乘格式zf(k)= h:(k)6 + v(k)(3.2.4)由于上式v(k)是白噪声，所以利用最小二乘法即可获得参数的无偏估计。但是，数据向量h (k)中的变量均需要按照(3.2.2)式计算，然而噪声模型C(z-1)并不知道。为此需要 f用迭代的方法来估计C(z-1)。令 e(k)=1C (z -1)v(k)(3.2.5)(3.2.6)m h (k) = -e(k - 1),-e(k - 2),-e(k - 3汗 置S ec6 (k) = c , c , c 卞Je123就把

15、噪声模型(3.2.5)也化成最小二乘格式e(k) =(k)6 + 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)t(3.2.7)e其中 e(k) = 0,k 0 时，按照式子e(k) = z(k) 一 h(k)9 (k)计算，式子中 h(k) = z (k 1),-z (k 2),-z (k 3), u (k 1), u (k 2), u (k 3”

16、综上分析，广义最小二乘递推算法可归纳成：八6(k)二 e(k -1) + K (k)z (k) -hT (k)6(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)6 (k -1)eeeeeP (k) = I K (k)hT (k)P (k 1)eeeeK (k) = P (k 1)h (k)hT (k)P (k 1)h (k) + 1-ieeeeee利用上述公式即可求得参数

17、 的估计值6并能求出噪声模型的估计。3、估计结果e = aa abb 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.06801.0666-0.5284 -0.58353.1471 -3.71856.2149 -6.30267.2705. .-9.05528.1735 -5.90043.9870 -2.2486

18、0.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.7599 -0.4739 -0.1784 -

19、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,:) = -z(i) -z(i-1) -z(i-2) u(i)

20、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.05528.1735 -5.90043.9870 -2.24860.9525 -0.5325 -1.52270.42001.0786.-1.5579 0.6640 -1.4222 2.

21、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.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591.-1.0576 -1.0071 1.1342 -0.0740

22、 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,:) = -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,

23、:)*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);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,:);e

24、nd%辨识结果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 -5.90043.9870 -2.24860.9525 -0.5325 -1.52270.42001.0786.-1.5579 0.6640 -1.4222 2.6444 -2.9

25、572 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.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591.-1.0576 -1.0071 1.1342 -0.0740 0.6759 0.

26、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;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

27、);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,40);pe(:,:,1) = eye(3)-ke(1,:)*he(1,:)*pe0;oe=zeros(3,40);oe(1,:)=oe0+ke(1,:)*e(1)-he(1)*oe0

28、;%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(1) 0 0;ke(2, :) = pe(:,:,2-1)*he(2,:)*inv(he(2,:)*pe(:,:,2-1)*he(2,:)+1);pe(:,:,2) = eye(3

29、)-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,:) = of(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;

30、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)*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,

31、:) = -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);