系統(tǒng)辨識作業(yè)及答案

1、一.問答題1 .介紹系統(tǒng)辨識的步驟。答：(1)先驗知識和建模目的的依據(jù)；(2)實驗設(shè)計；(3)結(jié)構(gòu)辨識；(4)參數(shù)估計；(5)模型適用性檢驗。2 .考慮單輸入單輸出隨機系統(tǒng)，狀態(tài)空間模型101x(k1)=20x(k)1u(k)y(k)=11X(k)v(k)轉(zhuǎn)換成ARMA模型。答：ARMA模型的特點是u(k)=0,10x(k1)x(k)20y(k)=11x(k)v(k)3.設(shè)有一個五級移位寄存器，反饋取自第2級和第3級輸出的模2加法和。試說明:(1) 其輸出序列是什么？(2) 是否是M序列？(3) 它與反饋取自第4級與第3級輸出模2加法和所得的序列有何不同？(4) 其逆M序列是什么？答：(1)設(shè)

2、設(shè)輸入序列11111CD1111190111017001112510011(2)01111100011118100112601001(3)00111111001119010012710100(4)1001112010012010100281101001001(13)1010021110102900111(6)10100141101022111013001110110101511101(23)01110310011111101160111024001113210011其輸出序列為：1111100101不是M序列第4級與第3級模2相加結(jié)果CD1111191100117011112501100011

3、11100110018001112610110(300111111011019000112701011(4)0001112010112010001281010110001(13)1010121010002911010(6)01000141101022001003011101001001511101(23)10010311111010010161111024110013201111不同點：第2級和第3級模二相加產(chǎn)生的序列，是從第4時刻開始，每隔7個時刻重復(fù)一次;15個時刻重復(fù)一次。第4級與第3級模2相加產(chǎn)生的，序列，是從第2時刻開始每隔第5級與第4級模2相加結(jié)果如下：已知其為M序歹U。11111

4、9000100111110100011701010(30011111110001810101(4)000111201100191101000001(13)001102011101(6)10000141001121111100100015010012211111001001610100M序列：方波信號：逆重復(fù)M:1111101001011000101000100100010101111014 .畫出廣義最小二乘法的離線迭代算法的簡單計算框圖。答：廣義最小二乘法的離線迭代算法的簡單計算框圖如下:5 .考慮如下數(shù)學(xué)模型y=asinx+bcosx,試用x(k),y(k),k=1,2,3,.,N輸入輸出

5、數(shù)據(jù)估計系統(tǒng)參數(shù)a,bo答：y1=asinx1bcosx1y2=asinx2bcosx2ayN=asinxNbcosxN二Y-二E二.I-(7-)-1-TY6 .利用最小二乘算法辨識如下模型參數(shù)z(k)-1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)其中，v(k)是零均值白噪聲。當(dāng)模型階次為2時，可以獲得準(zhǔn)確的辨識結(jié)果，而模型階次取3時，只能得到如下一組模型參數(shù)辨識結(jié)果(括號內(nèi)為模型參數(shù)真值)：a1=-1.08884(-1.5)a2=0.08326(0.7)a3=0.28781(0.0)b1=1.00000(1.0)b2=0.91116(0.5)b3=0.2

6、0558(0.0)顯然，辨識結(jié)果已經(jīng)遠遠偏離了模型參數(shù)真值，試從理論上解釋為什么會出現(xiàn)這種現(xiàn)象。答：又于n階系統(tǒng)與n+1階系統(tǒng)參數(shù)估計之間有如下的關(guān)系：對于n+1階系統(tǒng)Az1y(k)=Bzuke(k)設(shè)其待估參數(shù)為Tn1：-|boaibi.an4%dbn.1-丁(2)則u(1)=u(n)-A/Y-中Mn)由題目知n=2時系統(tǒng)參數(shù)為準(zhǔn)確值，則n=3時按照上式去計算，估算出的系數(shù)必遠遠偏離系統(tǒng)模型參數(shù)值。7 .請說明閉環(huán)系統(tǒng)不可辨識的原因。答：閉環(huán)系統(tǒng)不可辨識的原因：反饋使得一個閉環(huán)系統(tǒng)對不同的輸入常產(chǎn)生差不多相同的輸出，觀測的輸入輸出數(shù)據(jù)所包含的信息比開環(huán)辨識少的多；輸入信號與噪聲因反饋而相關(guān)：

7、有偏估計，非一致性估計；在閉環(huán)條件下，用開環(huán)辨識方法系統(tǒng)的參數(shù)有時也是不可唯一辨識的。8 .設(shè)閉環(huán)系統(tǒng)前向通道模型為y(k)-1.4y(k-1)-0.45y(k-2)u(k-1)0.7u(k-2)；(k)反饋調(diào)節(jié)器為u(k)=y(k)0.2y(k-1)試畫出其閉環(huán)系統(tǒng)框圖，并判斷系統(tǒng)是否可辨識？答：系統(tǒng)是可以辨識的，由于為非奇異，故在u(k)=y(k)+0.2y(k-1)條件下，參數(shù)是可以辨識的。閉環(huán)系統(tǒng)框圖如下圖所示：t(n,n1)二J(n)-J(n1)J(n1)9.對系統(tǒng)模型階次進行辨識，得到1階-4階的參數(shù)估計，性能指標(biāo)與系統(tǒng)模型階次的關(guān)系如下表所示，利用F檢驗法判斷系統(tǒng)模型的階次。n=

8、1n=2n=3n=4n=5J51.814.6312.4612.4112.40解：由F檢驗法原理知若t(n,n三3.09則可以接受系統(tǒng)階數(shù)。由計算得，t(1,2)=4.13,t(2,3)=0.49,t(3,4)=0.0034,t(4,5)=0所以系統(tǒng)的階數(shù)為3。二.編程題1. (1)編程產(chǎn)生一組正態(tài)分布的白噪聲信號，它的均值和方差以及長度可隨意調(diào)整，將產(chǎn)生的白噪聲信號存入數(shù)據(jù)文件data1.txt(2)編程產(chǎn)生一組M序列信號，它的幅值和長度可隨意調(diào)整，將產(chǎn)生的M序列存入數(shù)據(jù)文件data2.txt(3)編程產(chǎn)生一組逆重復(fù)M序列信號，它的幅值和長度可隨意調(diào)整，將產(chǎn)生的逆重復(fù)M序列存入數(shù)據(jù)文件data

9、3.txt解：(1)functiony=WNoise(N,E,VAR)%N為長度E為均值VAR為方差y=randn(1,N);y=y-mean(y);y=y/std(y);y=E+sqrt(VAR)*y;plot(y)title('嚴(yán)曉龍實驗：產(chǎn)生一組正態(tài)分布的白噪聲信號)savedata1.txty-ascii調(diào)用函數(shù)實驗：WNoise(400,0,1),得到數(shù)據(jù)見data1.txt,如圖所示:UFigure1".FileEditViewnsertToolsDesktopWindowHelpdlHIQlW裝國口國IQ嚴(yán)曉龍實瞼，產(chǎn)生一組正態(tài)分布的白噪聲信號：1iij.t.J

10、I."0501C0150200度。打。350400(2) functionseq=mseq(a,L,N)%a為M序列幅值N為長度L為移位單位數(shù)register=randint(1,L)%寄存器初始化p=zeros(1,L);%特征向量p(L-1:L)=1;%默認最后兩個寄存器相加temp=0;fori=1:Nseq(i)=a*register(L);temp=sum(register.*p);register(2:L)=register(1:L-1);%移位register(1)=mod(temp,2);endx=0:1/5:1.2;stairs(seq);grid;set(gca,

11、'ylim',-0.2,1.2);ylabel('M序列')title('嚴(yán)曉龍實驗：移位寄存器產(chǎn)生的M序列')savedata2.txtseq-ascii調(diào)用函數(shù)實驗：mseq(2,40,15),得到數(shù)據(jù)見data2.txt,和下圖:RieEditViewInsertToolsDesktopWindowHelpdd會昵Xb3I0嚴(yán)曉龍實驗：移位寄存器產(chǎn)生的M序列251oI0.武當(dāng)IrB1.i'Ih-.1-J051015(3) functionseq=invM(a,L,N)%a為M序列幅值N為長度L為移位單位數(shù)register=randi

12、nt(1,L)%寄存器初始化p=zeros(1,L);%特征向量p(L-1:L)=1;%默認最后兩個寄存器相加temp=0;fori=1:2AL-1seq(i)=register(L);temp=sum(register.*p);register(2:L)=register(1:L-1);%移位register(1)=mod(temp,2);endseq=seqseq;fori=1:2*(2AL-1)ifmod(i,2)=1invm(i)=1;elseinvm(i)=0;endseq(i)=a*xor(seq(i),invm(i);endfori=1:Nifmod(i,2*(2AL-1)=0m

13、seq(i)=seq(2*(2AL-1);elsemseq(i)=seq(mod(i,2*(2AL-1);endendseq=mseq;stairs(seq);grid;set(gca,'ylim',-0.2,1.2);title('嚴(yán)曉龍實驗：產(chǎn)生一組逆重復(fù)M序列信號,)savedata3.txtseq-ascii調(diào)用函數(shù)實驗：invM(1,10,40),得到數(shù)據(jù)見data3.txt,和下圖2.12.mat中的數(shù)據(jù)是單輸入單輸出系統(tǒng)進行采樣后100對輸入輸出數(shù)據(jù)，其中input表示系統(tǒng)的輸入數(shù)據(jù)，output表示受到噪聲污染后的系統(tǒng)的輸出數(shù)據(jù)。在辨識過程中，可以認為噪

14、聲具有正態(tài)分布，其均值為0。(1) 判斷該系統(tǒng)的階次(方法不限)(2) 利用遞推最小二乘法進行參數(shù)估計。解:J是模型階數(shù)的辨識，一般說來低階模型描述粗糙，高階模型精度高。殘差平方總和NT9模型階數(shù)的函數(shù)J(n)八.(y(k)-二三)2kX-在不同的模型階數(shù)的假設(shè)下，參數(shù)估計得到的J(n)值亦不同。討論如下當(dāng)n=1時程序如下：啟動matlab,打開12.mat;運行下面程序u=zeros(100,1);%構(gòu)造輸入矩陣z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;u(i,1)=input(i);z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;z(i,1)=outp

15、ut(i);r=100;forp=1:(r-2)%利用循環(huán)生成觀測矩陣h(p,:)=-z(p+1)u(p+1);%endhl=h;forb=1:(r-2)%生成輸出矩陣zl(b,:)=z(b+2);zl'endzl'%根據(jù)最小二乘法公式進行參數(shù)辯識c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;a1=1)a2=c(2)j=0;fork=4:100;hl=-z(k-1);u(k-1)'x=hl*c;y=z(k)-x;s=y*y;j=j+s;end仿真結(jié)果如下a1=-0.2576a2=0.6985j=0.8556al-0,257

16、6a2=0.69851 =0.855K(2)當(dāng)n=2時程序如下(輸入輸出數(shù)據(jù)同上，只給出不同于一階系統(tǒng)的程序不同之處)其中U、Z分別是作業(yè)要求給出得的輸入輸出，數(shù)據(jù)輸入同上。啟動matlab,打開12.mat;運行下面程序u=zeros(100,1);%構(gòu)造輸入矩陣z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;u(i,1)=input(i);z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;z(i,1)=output(i);r=100;%利用循環(huán)生成觀測矩陣。forp=1:(r-2)h(p,:)=-z(p+1)-z(p)u(p+1)u(p);endhl=h;%生

17、成輸出矩陣。forb=1:(r-2)zl(b,:)=z(b+2);zl'endzl'%根據(jù)最小二乘法公式進行參數(shù)辯識c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;%輸出辯識參數(shù)a1=31)a2=c(2)b1=c(3)b2=c(4)j=0；%求(n)fork=4:100;%開始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2)'x=h1'*c;y=z(k)-x;s=y*y;j=j+s;end仿真結(jié)果如下：a1=-0.4362a2=0.2407b1=1.8844b2=-1.1313j=0.5977當(dāng)n

18、=3時程序如下啟動matlab,打開12.mat;運行下面程序u=zeros(100,1);%構(gòu)造輸入矩陣z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;u(i,1)=input(i);z=zeros(100,1);%構(gòu)造輸出矩陣i=1:1:100;z(i,1)=output(i);r=100;forp=2:(r-1)h(p,:)=-z(p+1)-z(p)-z(p-1)u(p+1)u(p)u(p-1);endhl=h;forb=2:(r-1)zl(b,:)=z(b+1);zl'endzl'c1=hl'*hl;c2=inv(c1);c3=hl'*zl;c=c2*c3;a1=c(1)a2=32)a3=c(3)b1=c(4)b2=c(5)b3=c(6)j=0；