2008級碩士研究生《系統建模理論》試卷已知一個三階線性離散系統的輸入、輸出數據，共有40個采樣值，試分別用：最小二乘法(LS)、遞推最小二乘法(RLS)、廣義最小二乘法(61，)進行參數估計，給出相應的數學模型，并闡述相應的辨識原理。u(k)-0.9229ju(k)-0.9229j(k)-9.05521.6400-0.84100.75998.1735-5.90043.9870-0.4739-0.1784-1.7760-1.6722-2.24860.9525-0.5325-1.52271.2959-0.05910.42001.0786k1234 56 78910u(k)0.82510.09880.4628-0.91682.23250.07772.36540.34761.1473-1.9035y(k)1.5333-1.06801.0666-0.5284-0.58353.1471-3.71856.2149-6.30267.2705k11 121314 15 16 17 18 19 20k2122 2324252627282930u(k)-1.0576-1.00711.1342-0.07400.67590.52210.99540.5271-1.76560.4936y(k)-1.55790.6640-1.42222.6444-2.95723.6340-3.12813.8334-3.25421.1568k3132 3334353637383940u(k)1.48100.9591-3.1293-0.3604-0.42510.4185-0.6728-0.00272.11451.1157y(k)0.06150.9120-0.0692-3.27313.7486-4.31944.7230-5.27815.1507-2.7235一、最小二乘法(LS)1、數學模型設時不變SISO動態過程的數學模型為A(z-1)z(k)=B(z-i)u(k)+n(k) (1.1.1)其中，u(k)為過程的輸入量，z(k)為過程的輸出量，n(k)是噪聲，多項式A(z-1)和B(z-1)為：JA(z-1)=1+az-1+az-2+—Faz-naib(5=u+—+『在本題中，n=〃廣3.即JA(z-1)=1+匕z-1+a2z-2+a3z-3iB(z-1)=bz-1+bz-2+bz-3v 1 2 3(1.1.2)將此模型寫成最小二乘格式z(k)=hT(k)e+n((1.1.2)其中，z(k)是過程的輸出量；hT(k)是可觀測的數據向量；n(k)是均值為零的隨機噪聲。

式中h(k)=[—z(k—1),_z(k—2),-z(k—3),u(k—1),u(k—2),u(k—3)]t式中0=[a,a,a,b,b,bR偵 12 3 12 3對于k=1,2,…L，方程式(1.1.2)構成一個線性方程組，可以把它寫成匕=[z匕=[z(1),z(2),...,z(l計n=[n(1),n⑵,,n(l)k10 0 u(0)—z(0) 0 u(1)…—z(1—2) —z(1—3)u(1—1)0 0u(0) 0… (1.1.3)u(1—2)u(1—3)利用數據序列(z(k)｝和｛h(k)｝，極小化準則函數J(0)=W[z(k)—hT(k)0]2 (1.1.4)k=1.一_— . - . ..A . . 一一. .. .. .使J(0)=min的0估計值記作0，稱為參數0的取小二乘估計值。通過極小化(1.1.4)式來計算0的方法稱作最小二乘法，未知模型參數0最可能的值是在實際觀測值與計算值之累次誤差的平方和達到最小處所得到的，這種模型輸出能最好地接近實際過程的輸出。2、 辨識原理考慮模型(1.1.2)式的辨識問題，其中z(k)和h(k)都是可觀測的數據，0是待估計參數，準則函數取(1.1.4)根據(1.1.3)的定義，準則函數J(0)可寫成二次型的形式J(0)=(z1-HQ)T(z1-H10) (1.2.1)顯然上式中的H0代表模型的輸出，或者說是過程的輸出預報值。因此J(0)可以看作來衡量模型輸出與實際過程輸出的接近情況。極小化J(0)，求得參數0的估計值TOC\o"1-5"\h\z△ 一0=(HtH)-1Htz (1.2.2)1 1 11將使模型的輸出最好的預報過程的輸出。3、 辨識結果人 ............... ..... .....0=[aaabbb]=[-1.00000.00000.00000.00000.00000.0000]12 3 12 3二、遞推取小二乘法(RLS)1、數學模型在第一部分中建立了最小二乘法，并用一次完成算法進行了計算。但是由于具體使用時占用

內存量大，而且不能用于在線辨識。所以引入了最小二乘參數估計的遞推算法。遞推算法的基本思想可以概括成：'一、一..新的估計值。(k)=老的估計值。(k-1)+修正項在此算法中，時不變SISO動態過程的數學模型仍與最小二乘法的一樣。模型的最小二乘格式也相同，只是計算方法不同，具體計算方法，在遞推最小二乘法的辨識原理中描述。2、辨識原理首先將1.2.2式一次完成算法寫成0 ———— —— — 0 ———— —— — ——e=(HTH)-1HTz二P(l)HTziiii ii=[￡h(i)ht(i)]-1[￡h(i)z(i)]i=1 i=1(2.2.1)定義(2.2.2)P-1(k)=￡h(i)hT(i)=HtHkk(2.2.2)P-1(k-1)=吏h(i)hT(i)=HtH-hT⑴--hT⑴一hT⑵，H=hT⑵:.k-1:|_hT(k)JLhT(k-1)Ji=1其中：Hk由(2.2.2)式可得:P-1(k)=歹1h(i)ht(i)+h(k)ht(k)=P-1(k-1)+h(k)ht(k)(2.2.3)i=1令：zk-1=[z⑴，z(2),???z(k-1)]T則：e(k-1)=(HTH )-1HTz=P(k-1)[蕓h(i)z(i)]k-1k-1 k-1k-1i=1于是有P-1(k-1)6(k-1)=￡h(i)z(i) (2.2.4)i=1令zk=[z⑴，z(2), z(k)]T利用(2.2.3)和(2.2.4)式，可得e(k)=(HTH)-1HTz=P(k)芭h(i)z(i)]kkkki=1_ _ _ _ 0_ - _ _=P(k)[P-1(k-1)6(k-1)+h(k)z(k)] (2.2.5)=P(k){[P-1(k)-h(k)hT(k)]6(k-1)+h(k)z(k)}=e(k-1)+p(k)h(k)[z(k)-hT(k)6(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)-h「(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-i(k-1)+h(k)ht(k)]-i=[I-K(k)ht(k)]P(k-1) (2.2.9)將(2.2.9)式代入(2.2.6)式，整理后有K(k)=P(k-1)h(k)[ht(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)-ht(k)0(k-1)]P(k)=[P-1(k-1)+h(k)ht(k)]-1=[I-K(k)ht(k)]P(k-1)K(k)=P(k-1)h(k)[hT(k)P(k-1)h(k)+1]-1 .入利用上述公式即可求得參數0的估計值03、辨識結果△................0=[a aabbb]=[-1.0005-0.0007-0.00020.0000-0.0005-0.0001]12 3 12 3三、廣義最小二乘法(GLS)1、數學模型設時不變SISO動態過程的數學模型為A(z-1)z(k)=B(z-1)u(k)+ 1 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+az-1+az-2+ +az-匕<B(z-1)=bz-1+bz-2+ +bz-金1 2 駐C(z-1)=1+C]z-1+C2z-2+ +cz-nc

在本題中，在本題中，n=n.=n=3.即A(z-i在本題中，在本題中，n=n.=n=3.即A(z-i)=1+az-i+az-2+az-31 2 3<B(z-1)=bz-1+bz-2+bz-3「1 2 3 (3.1.2)C(z-1)=1+cz-1+cz-2+cz-3I 1 2 3廣義最小二乘發的基本思想是基于對數據先進行一次濾波預處理，然后利用普通最小二乘法對濾波后的數據進行辨識。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)|z(k)=C(z-1)z(k)u(k)=C(z-1)u(k)，1f(3.2.2)h(k)=[—z(k—1),-z(k—2),-z(k—3),u(k—1),u(k—2),u(k—3)]tf f f f ff f (3.2.3)0=[a,a,a,b,b,b卜- 12 3 12 3可將模型(3.1.1)式化成最小二乘格式zf(kNhfT(k)0+v(k)(3.2.4)由于上式v(k)是白噪聲，所以利用最小二乘法即可獲得參數0的無偏估計。但是，數據向量hf(k)中的變量均需要按照(322)式計算，然而噪聲模型C(z-1)并不知道。為此需要用迭代的方法來估計C(z-1)。令e(k)=—-——v(k)

C(z-1)(3.2.5)(3.2.6)h(k)=[-e(k-1),-e(k-2),-e(k-3)卡(3.2.6)0(k)=[匕,c2,c3k就把噪聲模型(3.2.5)也化成最小二乘格式e(k)=he(k)0°+v(k)由于上式的噪聲已是白噪聲，所以再次利用最小二乘法可獲得噪聲模型參數七的無偏估計。但是，數據向量h^(k)依然包含著不可測得噪聲量e(k-1),e(k-2),e(k-3)，它可用相應(3.2.7)的估計值代替，置h^(k)=[-e(k-1),-e(k-2),-e(k-(3.2.7),._ . _ _ _ FC .. 人其中e(k)=0,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)]e綜上分析，廣義最小二乘遞推算法可歸納成：0(k)=e(k-1)+Kf(k)[%(k)-hf(k渺(k-1)]P(k)=[I-Kf(k)hf(k)]p(k-1)Kf(k)=P(k-1)hf(k)[hf(k)p(k-1)hf(k)+1]-iTOC\o"1-5"\h\z\o"CurrentDocument"△ 一一. _ .fe(k)=e(k-1)+k(k)[e(k)-hf(k)e(k-1)]e e e ee\o"CurrentDocument"P(k)=[I-K(k)hf(k)]P(k-1)e ee e\o"CurrentDocument"K(k)=P(k-1)h(k)[hf(k)P(k-1)h(k)+1]-1e e e ee e ..一A... . ._一一利用上述公式即可求得參數e的估計值e并能求出噪聲模型的估計。3、估計結果0=[a aabbb]=[-1.0005-0.00080.0002-0.0000-0.0005-0.0001]12 3 12 3一 一~噪聲模型的估計值為1.0e-003*[-0.0000 0.4115 0.0001]即C(z-1)就1附錄：MATLAB程序1、 一次完成的最小二乘法%%一次完成的最小二乘法clears=(1:40);z=[1.5333-1.06801.0666-0.5284-0.58353.1471-3.71856.2149-6.3026.-9.05528.1735-5.90043.9870-2.24860.9525-0.5325-1.52270.42001.0786...-1.55790.6640-1.42222.6444-2.95723.6340-3.12813.8334-3.25421.1568...0.06150.9120-0.0692-3.27313.7486-4.31944.7230-5.27815.1507-2.7235];u=[0.82510.09880.4628-0.91682.23250.07772.36540.34761.1473-1.9035...-0.92291.6400-0.84100.7599-0.4739-0.1784-1.7760-1.67221.2959-0.0591...-1.0576-1.00711.1342-0.07400.67590.52210.99540.5271-1.76560.4936...0.9591-3.1293-0.3604-0.42510.4185-0.6728-0.00272.11451.1157];h=zeros(40,6);h(1,:)=[-z(1)00u(1)00];h(2,:)=[-z(2)-z(1)0u(2)u(1)0];fori=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;fori=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.3026.-9.05528.1735-5.90043.9870-2.24860.9525-0.5325-1.52270.42001.0786...-1.55790.6640-1.42222.6444-2.95723.6340-3.12813.8334-3.25421.1568...0.06150.9120-0.0692-3.27313.7486-4.31944.7230-5.27815.1507-2.7235];u=[0.82510.09880.4628-0.91682.23250.07772.36540.34761.1473-1.9035...-0.92291.6400-0.84100.7599-0.4739-0.1784-1.7760-1.67221.2959-0.0591...-1.0576-1.00711.1342-0.07400.67590.52210.99540.5271-1.76560.4936...0.9591-3.1293-0.3604-0.42510.4185-0.6728-0.00272.11451.1157];h=zeros(40,6);h(1,:)=[-z(1)00u(1)00];h(2,:)=[-z(2)-z(1)0u(2)u(1)0];fori=3:1:40h(i,:)=[-z(i)-z(i-1)-z(i-2)u(i)u(i-1)u(i-2)];endp0=10"6;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];fori=2:1:40k(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)-h(i,:)*o(i-1,:)']]';end%辨識結果o(40,:)%誤差分析J=0;fori=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.3026.-9.05528.1735-5.90043.9870-2.24860.9525-0.5325-1.52270.42001.0786...-1.55790.6640-1.42222.6444-2.95723.6340-3.12813.8334-3.25421.1568...0.06150.9120-0.0692-3.27313.7486-4.31944.7230-5.27815.1507-2.7235];u=[0.82510.09880.4628-0.91682.23250.07772.36540.34761.1473-1.9035...-0.92291.6400-0.84100.7599-0.4739-0.1784-1.7760-1.67221.2959-0.0591...-1.0576-1.00711.1342-0.07400.67590.52210.99540.5271-1.76560.4936...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)00uf(1)00];hf(2,:)=[-zf(2)-zf(1)0uf(2)uf(1)0];o0=0.001;pf0=10"6;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,:)=[000];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];%k=2時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)=z(2)-hf(2,:)*of(2,:)';he(2,:)=[-e(1)00];ke(2,:)=[pe(:,:,2T)*he(2,:)'*inv(he(2,:)*pe(: