1、2008級碩士研究生系統(tǒng)建模理論試卷已知一個(gè)三階線性離散系統(tǒng)的輸入、輸出數(shù)據(jù)，共有40個(gè)采樣值，試分別用：最小二乘法（LS）、遞推最小二乘法（RLS）、廣義最小二乘法（GLS）進(jìn)行參數(shù)估計(jì)，給出相應(yīng)的數(shù)學(xué)模型，并闡述相應(yīng)的辨識原理。k 1 2 3 4 5 6 7 8 9 10u(k) 0.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035y(k) 1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705k 11 12

2、 13 14 15 16 17 18 19 20u(k) -0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0.0591y(k) -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786k 21 22 23 24 25 26 27 28 29 30u(k) -1.0576 -1.0071 1.1342 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936y(k) -1

3、.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8334 -3.2542 1.1568k 31 32 33 34 35 36 37 38 39 40u(k) 1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.1157y(k) 0.0615 0.9120 -0.0692 -3.2731 3.7486 -4.3194 4.7230 -5.2781 5.1507 -2.7235一、最小二乘法（LS）1、數(shù)學(xué)模型設(shè)時(shí)不變SISO動態(tài)過程的數(shù)學(xué)模型為 （）其

4、中，為過程的輸入量，為過程的輸出量，是噪聲，多項(xiàng)式和為：在本題中，=3.即將此模型寫成最小二乘格式 （）其中，是過程的輸出量；是可觀測的數(shù)據(jù)向量；是均值為零的隨機(jī)噪聲。式中對于，方程式（）構(gòu)成一個(gè)線性方程組，可以把它寫成（）利用數(shù)據(jù)序列和，極小化準(zhǔn)則函數(shù) （）使的估計(jì)值記作，稱為參數(shù)的最小二乘估計(jì)值。通過極小化（）式來計(jì)算的方法稱作最小二乘法，未知模型參數(shù)最可能的值是在實(shí)際觀測值與計(jì)算值之累次誤差的平方和達(dá)到最小處所得到的，這種模型輸出能最好地接近實(shí)際過程的輸出。2、辨識原理考慮模型（）式的辨識問題，其中和都是可觀測的數(shù)據(jù)，是待估計(jì)參數(shù)，準(zhǔn)則函數(shù)取（）根據(jù)（）的定義，準(zhǔn)則函數(shù)可寫成二次型的形式

5、 （）顯然上式中的代表模型的輸出，或者說是過程的輸出預(yù)報(bào)值。因此可以看作來衡量模型輸出與實(shí)際過程輸出的接近情況。極小化，求得參數(shù)的估計(jì)值 （）將使模型的輸出最好的預(yù)報(bào)過程的輸出。3、 辨識結(jié)果二、遞推最小二乘法（RLS）1、數(shù)學(xué)模型在第一部分中建立了最小二乘法，并用一次完成算法進(jìn)行了計(jì)算。但是由于具體使用時(shí)占用內(nèi)存量大，而且不能用于在線辨識。所以引入了最小二乘參數(shù)估計(jì)的遞推算法。遞推算法的基本思想可以概括成：新的估計(jì)值=老的估計(jì)值+修正項(xiàng) （）在此算法中，時(shí)不變SISO動態(tài)過程的數(shù)學(xué)模型仍與最小二乘法的一樣。模型的最小二乘格式也相同，只是計(jì)算方法不同，具體計(jì)算方法，在遞推最小二乘法的辨識原理中

6、描述。2、辨識原理首先將式一次完成算法寫成 （）定義 （）其中：，由（）式可得：（）令：則：于是有 （）令利用（）和（）式，可得 （）引進(jìn)增益矩陣，定義為 （）則（）式寫成 （）進(jìn)一步把（）式寫成 （）為了避免矩陣求逆運(yùn)算，利用矩陣反演公式可將（）式演變成 （）將（）式代入（）式，整理后有 （）綜合（）、（）、（）式便得到最小二乘參數(shù)估計(jì)遞推算法。利用上述公式即可求得參數(shù)的估計(jì)值3、辨識結(jié)果三、廣義最小二乘法（GLS）1、數(shù)學(xué)模型設(shè)時(shí)不變SISO動態(tài)過程的數(shù)學(xué)模型為 （）其中，和 分別表示過程的輸入和輸出，是均值為零的不相關(guān)隨機(jī)噪聲，多項(xiàng)式、和為：在本題中，在本題中，=3.即 （）廣義最小二乘

7、發(fā)的基本思想是基于對數(shù)據(jù)先進(jìn)行一次濾波預(yù)處理，然后利用普通最小二乘法對濾波后的數(shù)據(jù)進(jìn)行辨識。2、辨識原理由（）式可得 （）令， （）及（）可將模型（）式化成最小二乘格式 （）由于上式是白噪聲，所以利用最小二乘法即可獲得參數(shù)的無偏估計(jì)。但是，數(shù)據(jù)向量中的變量均需要按照（）式計(jì)算，然而噪聲模型并不知道。為此需要用迭代的方法來估計(jì)。令 （）置 （）就把噪聲模型（）也化成最小二乘格式由于上式的噪聲已是白噪聲，所以再次利用最小二乘法可獲得噪聲模型參數(shù)的無偏估計(jì)。但是，數(shù)據(jù)向量依然包含著不可測得噪聲量，它可用相應(yīng)的估計(jì)值代替，置 （）其中，當(dāng)時(shí)，按照式子計(jì)算，式子中綜上分析，廣義最小二乘遞推算法可歸納成：

8、利用上述公式即可求得參數(shù)的估計(jì)值并能求出噪聲模型的估計(jì)。3、估計(jì)結(jié)果噪聲模型的估計(jì)值為1.0e-003 * -0.0000 0.4115 0.0001 即附錄：MATLAB程序1、一次完成的最小二乘法%一次完成的最小二乘法clears=(1:40);z=1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705. -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786.-1.5579 0.6640 -1.4222

9、 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.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.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.13

10、42 -0.0740 0.6759 0.5221 0.9954 0.5271 -1.7656 0.4936.1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.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);endo=inv(h'*h)*h'*z'%誤差分析J

11、=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.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705. -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0786.-1.5579 0.6640 -1.4222 2.6444 -2.9572 3.6340 -3.1281 3.8

12、334 -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.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.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.

13、5271 -1.7656 0.4936.1.4810 0.9591 -3.1293 -0.3604 -0.4251 0.4185 -0.6728 -0.0027 2.1145 1.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,:)*p0*h(1,:)'+1)

14、9;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:40 k(i,:)=p(:,:,i-1)*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)

15、-h(i,:)*o(i-1,:)''end%辨識結(jié)果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、 廣義最小二乘法%廣義預(yù)測的最小二乘法clears=(1:40);z=1.5333 -1.0680 1.0666 -0.5284 -0.5835 3.1471 -3.7185 6.2149 -6.3026 7.2705. -9.0552 8.1735 -5.9004 3.9870 -2.2486 0.9525 -0.5325 -1.5227 0.4200 1.0

16、786.-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.8251 0.0988 0.4628 -0.9168 2.2325 0.0777 2.3654 0.3476 1.1473 -1.9035.-0.9229 1.6400 -0.8410 0.7599 -0.4739 -0.1784 -1.7760 -1.6722 1.2959 -0

17、.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.4251 0.4185 -0.6728 -0.0027 2.1145 1.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時(shí)kf(1,:)=pf0*hf(1,:)

18、'*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,40);pe(:,

19、:,1)=eye(3)-ke(1,:)'*he(1,:)*pe0;oe=zeros(3,40)'oe(1,:)=oe0+ke(1,:)'*e(1)-he(1)*oe0;%k=2時(shí)kf(2,:)=pf(:,:,2-1)*hf(2,:)'*inv(hf(2,:)*pf(:,:,2-1)*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)

20、=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)-ke(2,:)'*he(2,:)*pe(:,:,2-1);oe(2,:)=oe(2-1,:)'+ke(2,:)'*e(2)-he(2,:)*oe(2-1,:)''%k=3時(shí)kf(3,:)=pf(:,:,3-1)*hf(3,:)'*inv(hf(3,:)*pf(:,:,3-1)

21、*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;ke(3,:)=pe(:,:,3-1)*he(3,:)'*inv(he(3,:)*pe(:,:,3-1)*he(3,:)'+1)'pe(:,:,3)=eye(3)-ke(3,:)'*he(3,:)*pe(:,:,3-1);oe(3,:)=oe(3-1,:)'+ke(3,:)'*e(3)-he(3,:)*oe(3-1,:)