系統(tǒng)辨識實驗報告1

1、系統(tǒng)辨識實驗報告學(xué)院：信息科學(xué)與技術(shù)學(xué)院專業(yè)：自動化日期：2016/4/26目錄實驗14一實驗內(nèi)容及要求：4二實驗原理：4三軟件設(shè)計思想：4四程序結(jié)構(gòu)框圖：5五運行示意圖：5實驗28一實驗內(nèi)容及要求：8二實驗原理：8三軟件設(shè)計思想：9四程序設(shè)計框圖：9五程序運行流程圖：10實驗312一實驗內(nèi)容及要求：12二實驗原理：12三程序數(shù)據(jù)流程圖：12四實驗運行結(jié)果：13實驗414一實驗內(nèi)容及要求：14二實驗原理：14三數(shù)據(jù)遞推關(guān)系圖：14四實驗運行結(jié)果：15心得體會16附錄（實驗代碼）171.LabWork1172.LabWork2213.LabWork3234.LabWork426實驗1一實驗內(nèi)容及

2、要求：1.編出矩陣A與B相乘得到的矩陣R的運算計算機程序要求：（1）A和B的維數(shù)及數(shù)值可通過鍵盤及數(shù)據(jù)文件輸入（2）計算結(jié)果R可由屏幕及文件輸出2.將1改寫為子程序3查找有關(guān)的資料，讀懂及調(diào)通矩陣求逆程序，并改寫為子程序。二實驗原理：1. 兩個矩陣A、B相乘得到C矩陣，首先要滿足的條件是A的列數(shù)與B行數(shù)相等，否則不能相乘。當滿足條件后，根據(jù)C（i,k）=Ai,j*B(j,k)可以求得C矩陣。2. 當求矩陣的逆時，首先要判斷其是否為方陣，若是則可以對其進行下一步的操作。本次實驗中求逆主要是通過構(gòu)造一個增廣矩陣（FangZ | E）矩陣的初等行變換得到（E | FZNi）的這樣的一個矩陣就可以求得

3、矩陣的逆。若矩陣FangZ不是滿秩矩陣時，F(xiàn)angZ沒有FZNi 。通過這樣的求逆方式，避免了大方陣的求取行列式運算。三軟件設(shè)計思想：1. 確定該軟件的功能主要有：鍵盤輸入兩個矩陣然后相乘；文本data輸入兩個矩陣將結(jié)果放在文本result中；鍵盤輸入一個方陣求得其逆矩陣。其中前兩個的矩陣相乘運算部分設(shè)置為一個函數(shù)Mul。2. 在main函數(shù)中提供兩個關(guān)于矩陣的選擇：multiplitation；invertion。其相對應(yīng)的子函數(shù)為MulOp(a,b,c)，Inv()。3. 在MulOp(a,b,c)子函數(shù)中，有兩種輸入矩陣的方式：way1，way2。相對應(yīng)的功能為鍵盤輸入，文本輸入。并且兩

4、者在處理矩陣時，都調(diào)用了Mul函數(shù)。4. 在Inv()子函數(shù)中，輸入和顯示原矩陣，和其相應(yīng)的逆矩陣。調(diào)用qiuni（double FangZM,double FZNiM,int n）子函數(shù)，可以都得到原矩陣的逆矩陣。但當原矩陣不可逆時，系統(tǒng)輸出為”The array is not invertible!“。四程序結(jié)構(gòu)框圖：way1鍵盤輸入矩陣Mul矩陣乘積子函數(shù)MulOp子函數(shù)矩陣乘積way2文本輸入矩陣main函數(shù)Inv子函數(shù)矩陣求逆qiuni子函數(shù)退出程序五運行示意圖：1. main函數(shù)的主界面：2. MulOp子函數(shù)的界面：3. Inv子函數(shù)的界面：4. 通過鍵盤操作計算兩個矩陣的乘積：

5、5. 求方陣的逆矩陣：實驗2一實驗內(nèi)容及要求：編寫并調(diào)試動態(tài)模型仿真程序：模型：y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)已知：白噪聲v(k)數(shù)據(jù)文件為DV，數(shù)據(jù)長度為L=500要求：（1） 產(chǎn)生長度為L的M序列數(shù)據(jù)文件DU （2） 產(chǎn)生長度為L的模型輸出數(shù)據(jù)文件DY 二實驗原理：由于在現(xiàn)實中，白噪聲序列很難求，所以尋找到M序列在一定程度上可以代替白色噪聲序列。由L=500，所以n=9。根據(jù)M序列的特征方程：fx=c0+c1*x+c2*x2+cnxn可知9階移位寄存器的多項式為fx=x9+x4+1，及可得c=0,0,0,1,0,0,0,0,19

6、級線性移位寄存器：+C8C1C0輸出a0a1a8a7圖中Ci表示反饋的兩種可能連接方式，Ci=1表示連線接通，第9-i級加入反饋中；Ci=0表示連線斷開，第9-i級未參加反饋。系統(tǒng)產(chǎn)生M序列的結(jié)構(gòu)流程圖：初始化寄存器a=1,1,0,0,0,1,0,1,0移位寄存器c=0,0,0,1,0,0,0,0,1M序列的長度為L=500i=1Y輸出M序列i<LNM序列DU(i)=a(i)寄存器向前移1位an=i=0n-1ciai(mod 2)三軟件設(shè)計思想：1. 該軟件的主要功能是：產(chǎn)生M序列賦給u(k)保存在DU.txt文件中；由u(k)和v(k)求得y(k)保存在DY.txt文件中。2. 在ma

7、in函數(shù)中給出3個選擇：求u(k)；求y(k)；退出程序。其相對應(yīng)的函數(shù)名稱為gener，ouput，exit。3. 在gener子函數(shù)中產(chǎn)生M序列u(k)保存到DU.txt文本文件中。4. 在output子函數(shù)中，通過對input子函數(shù)（讀入v(k),u(k)的數(shù)據(jù)）、deal子函數(shù)(由公式y(tǒng)(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)求y(k)的調(diào)用來達到生成y(k)序列并保存到DY.txt文本文件中。四程序設(shè)計框圖：gener子函數(shù)產(chǎn)生DUmain函數(shù)input子函數(shù)輸入DU，DVoutput子函數(shù)產(chǎn)生DYdeal子函數(shù)得DY退出程序五程序運

8、行流程圖：1. main函數(shù)：2. gener函數(shù)：3. output函數(shù)：實驗3一實驗內(nèi)容及要求：編寫并調(diào)試動態(tài)離散時間模型LS成批算法程序。要求：（1）原始數(shù)據(jù)由DU和DY讀出。（2）調(diào)用求逆及相乘子程序。（3）顯示參數(shù)辨識結(jié)果。二實驗原理：1.批次處理的方法就是把所有的數(shù)據(jù)采集到一次性進行處理，但前提是白色噪聲、及M序列所共同作用而產(chǎn)生的輸出，才能使用最小二乘法。雖然這種方法的計算量龐大，但經(jīng)常用于處理時不變系統(tǒng)，方法簡單。2.構(gòu)造模型Y=X*sita+V3Y=DY(3) DY(4) DY(499) 已知n=2,L=500,可知m=497所以有X=-y(2)-y(1)-y(3)-y(2)

9、-y(498)-y(497)u(2)u(1)u(3)u(2)u(498)u(497)三程序數(shù)據(jù)流程圖：轉(zhuǎn)置Y(m,1)列矩陣TX(4,m)矩陣X(m,4)矩陣xishu(4,1)列矩陣相乘相乘ji2(4,m)矩陣求逆相乘ni(4,4)方陣ji1(4,4)方陣四實驗運行結(jié)果：實驗4一實驗內(nèi)容及要求：編寫并調(diào)試動態(tài)離散時間模型LS遞推算法程序。要求：（1）原始數(shù)據(jù)由DU和DY讀出。（2）顯示參數(shù)辨識結(jié)果。（3）設(shè)置選擇變量決定是否輸出中間結(jié)果。二實驗原理：1.基本思路：新的估計值sita'(k+1)=老的估計值sita'(k)+修正項1. 基本模型：y(k)=-a1y(k-1)-a

10、2y(k-2)+b1u(k-1)+b2u(k-2)+v(k)。2. 引入信息矩陣P(k)，維數(shù)為4*4，初始化為1051012*E的一單位陣。gama(k)為修正系數(shù)，為無窮小標量。x(k+1)=-y(k) -y(k-1) u(k) u(k-1)。sita'=a1 a2 b1 b2。y(k+1)=DYk+3。3. sita'(k+1)=sita'(k)+gama(k+1)*P(k)*x(k+1)*y(k+1)-x'(K=1)*sita'(k)；4. P(k+1)=P(k)-gama(k+1)*P(k)*x(k+1)*x'(k+1)*P(k)；5.

11、 gama(k+1)=1/1+x'(k+1)*P(k)*x(k+1);三數(shù)據(jù)遞推關(guān)系圖：四實驗運行結(jié)果：1.main函數(shù)界面示意圖：6. 只顯示結(jié)果界面：7. 顯示過程的參考界面：心得體會我很喜歡這個課程的期末考核方式，不用再拘泥于在題目當中對該課程的了解，而是通過4個C語言設(shè)計的練習來達到學(xué)習的目的，而且對以后的學(xué)習還有很大的幫助。在編寫C語言的過程中，也遇到了一些阻礙，特別是在編寫第3,4個的時候。比如：用的數(shù)組太多，并且未將其初始化，運行出來的結(jié)果經(jīng)常是很長的一段隨機數(shù)；或者一模一樣的程序有時候就可以正常運行，有時候就總是出現(xiàn)報錯這些都是讓我心塞了兩周的問題。在這些問題解決之后，

12、運行出來的結(jié)果卻與實際模型參數(shù)的出入有點大，于是又重新查找第2個實驗是否是M序列產(chǎn)生的方式有問題。通過對初始化寄存器賦給不同的值，可以讓結(jié)果與真實模型參數(shù)之間的誤差達到最小。通過這學(xué)期對系統(tǒng)辨識這個課程的學(xué)習，讓我了解到系統(tǒng)辨識在建立數(shù)學(xué)模型的方面的重要性。對于不了解系統(tǒng)的工作機理的人來說，也有了一個可以知道系統(tǒng)模型的構(gòu)建方法。但是也必須要我們對所構(gòu)建模型有一個較為清醒的認識，比如要知道模型的類型，如果是動態(tài)模型，則要知道模型的參數(shù)或者階次（但是權(quán)函數(shù)模型就不需要事前知道模型的參數(shù)或者階次，不知道這種方法有沒有什么我們不知道的缺陷）。在處理數(shù)據(jù)上，又一次讓我了解到世上沒有兩全其美的事情，計算的

13、復(fù)雜程度和精度就好像魚和熊掌不能兼得。也正是因為這樣，才會成為促使人們在這方面的不斷尋找最優(yōu)化算法的動力。在對模型進行研究時，都是從最簡單的模型開始研究，比如在沒有噪聲的情況下，所得到的模型參數(shù)就為系統(tǒng)模型真實參數(shù)。進一步，在白色噪聲的情況下，所得到的參數(shù)就為系統(tǒng)模型真實參數(shù)的估計通過這樣的推理，就可以得到LS在系統(tǒng)辨識問題上的普適性，不管系統(tǒng)受怎樣的干擾，只要通過一系列的變換，最后還是通過LS來解決模型參數(shù)的確定工作。雖然我們對化學(xué)工藝不怎么了解，但通過學(xué)習系統(tǒng)辨識，讓我們對已經(jīng)投運后的工業(yè)生產(chǎn)設(shè)備的性能有一個跟蹤性的了解，可以對系統(tǒng)做出一個較為完整的評估，從而來改善生產(chǎn)工藝或者是及時更換生

14、產(chǎn)設(shè)備，從而來達到最大的經(jīng)濟效益。這學(xué)期主要就是介紹方法，雖然可以處理一些特殊情況的有噪聲和沒有噪聲的系統(tǒng)，但對于那些噪聲模型都不可預(yù)測的系統(tǒng)中，還是顯得無能為力。希望以后還有機會學(xué)習其他的系統(tǒng)辨識方法。附錄（實驗代碼）1. LabWork1#define M 100double FangZMM=0,FZNiMM=0;/* 求矩陣的乘法函數(shù) */void Mul(double aM,double bM,double cM,int n1,int n2,int n4)int i,j,k;for(i=0;i<n1;i+)for(k=0;k<n4;k+)for(j=0;j<n2;j+

15、)cik+=aij*bjk;/* 求方陣的逆 */void qiuni(double FangZM,double FZNiM,int n)double EMM,value=1.0,stemp,temp=0.0,san,cMM;int i,j,row,nextrow,flag=0;int col,switchtime=0;/* 構(gòu)造原矩陣的增廣矩陣 */for(i=0;i<n;i+)for(j=0;j<n;j+)if(j=i) Eij=1;else Eij=0;for(i=0;i<n;i+)for(j=0;j<n;j+)cij=FangZij;for(j=0;j<n

16、;j+)cij+n=Eij;/* 對增廣矩陣進行運算（FangZ|E）->（E|FZNi） */for(row=0;row<n-1;row+)nextrow=row+1;if(cnextrowrow=0)while(cnextrowrow=0)nextrow+;if(nextrow=n)flag=1;break;if(flag=1) continue;switchtime+;/* 交換第一個數(shù)字為0的行 */for(col=0;col<2*n;col+)stemp=crowcol;crowcol=cnextrowcol;cnextrowcol=stemp;for(nextro

17、w=row+1;nextrow<n;nextrow+)temp=cnextrowrow*1.0/crowrow;for(col=0;col<2*n;col+)cnextrowcol+=-temp*crowcol;/* 得到前一個矩陣的對角線上元素的乘積 */for(row=0;row<n;row+)value*=crowrow;/* 判斷乘積！若為零，則原矩陣沒有逆矩陣；反之亦然 */if(value=0) printf("nnThe array is not invertible!nn");elsefor(row=n-1;row>=1;row-)f

18、or(i=row-1;i>=0;i-) san=cirow*1.0/crowrow;for(col=row;col<2*n;col+)cicol-=crowcol*san;/* 得到標準的（E|A） */for(row=0;row<n;row+)for(col=0;col<n;col+)crowcol+n/=crowrow;for(i=0;i<n;i+)for(j=0;j<n;j+)FZNiij=cij+n;2. LabWork2/* 處理v(k),u(k)的數(shù)據(jù)得輸出y(k) */void deal(double DVM,double DUM,double

19、 DYM)int k;DY0=0;DY1=0;for(k=2;k<500;k+)DYk=1.5*DYk-1-0.7*DYk-2+DUk-1+0.5*DUk-2+DVk;/* 產(chǎn)生m序列u(k)輸出到文件 */void gener()int a9=1,1,0,0,0,1,0,1,0,c9=0,0,0,1,0,0,0,0,1,aa9;int DUM,temp,sum,i,j;DU0=a8;FILE *fp;system("cls");if(fp=fopen("DU.txt","w")=NULL)printf("cannot

20、 open the file!n");exit(0);for(i=1;i<M;i+)sum=0;for(j=0;j<9;j+)sum+=aj*cj;aa0=sum%2;for(j=1;j<9;j+)aaj=aj-1;for(j=0;j<9;j+)aj=aaj;DUi=a8;for(i=0;i<M;i+)fprintf(fp,"%d ",DUi);fclose(fp);printf("nnnttPlease refer to the file(DU.txt)!nn");system("pause"

21、);3. LabWork3#define M 500#define N 4void MUL1(double TXM,double XN,double ji1N)int i,j,k;for(i=0;i<N;i+)for(k=0;k<N;k+)for(j=0;j<M-3;j+)ji1ik+=TXij*Xjk;void MUL2(double niN,double TXM,double ji2M)int i,j,k;for(i=0;i<N;i+)for(k=0;k<M-3;k+)for(j=0;j<N;j+)ji2ik+=niij*TXjk;void MUL3(d

22、ouble ji2M,double YM,double xishuN)int i,j;for(i=0;i<N;i+)for(j=0;j<M-3;j+)xishui+=ji2ij*Yj;void main()system("cls");system("color 75");int i,j;double DUM,DYM;double XMN=0,temp,YM,TXNM,ji1NN=0,niNN=0,ji2NM=0,xishuN=0;read(DU,DY);/* 得到大矩陣X( 498 * 4 ) */for(i=0;i<M-3;i+)fo

23、r(j=0;j<2;j+)Xij=-DYi+j+1;for(j=0;j<2;j+)Xij+2=DUi+j+1;temp=0;temp=Xi0;Xi0=Xi1;Xi1=temp;temp=0;temp=Xi2;Xi2=Xi3;Xi3=temp;/* 得到Y(jié)（498 * 1）=DY(1) DY(498)的列矩陣 */for(i=0;i<M-3;i+)Yi=DYi+3;/* 得到X（498 * 4）的轉(zhuǎn)置矩陣 TX（4 * 498） */for(i=0;i<N;i+)for(j=0;j<M-3;j+)TXij=Xji; MUL1(TX,X,ji1); /* 得到TX*X

24、的乘積ji1 */qiuni(ji1,ni,N); /* 得到j(luò)i1的逆矩陣ni */MUL2(ni,TX,ji2); /* 得到ni*TX的乘積ji2 */MUL3(ji2,Y,xishu); /* 得到j(luò)i2*Y的乘積xishu */printf("tthe xishu are:nn");for(i=0;i<N;i+)printf("nnt %.4f n",xishui);printf("nn");4. LabWork4#define M 500#define N 4double sitaNN;void MUL(double

25、 aN,double bN,double cN,int n1,int n2,int n4)int i,j,k;for(i=0;i<n1;i+)for(k=0;k<n4;k+)for(j=0;j<n2;j+)cik+=aij*bjk;void trans(double dN,double eN,int n5,int n6)int i,j;for(i=0;i<n6;i+)for(j=0;j<n5;j+)eij=dji;void consult()system("cls");int i,j,k;double DUM,DYM;double txNN=0

26、,pNN,y,xNN,temp,gama;double ji1NN=0,ji2NN=0;double ji3NN=0,ji4NN=0;double ji5NN=0,ji6NN=0,ji7NN=0;read(DU,DY);/* 對P（ 4 * 4 ）和sita進行初始化 */for(i=0;i<N;i+)for(j=0;j<N;j+)if(i=j) pij=10e9;else pij=0;for(i=0;i<N;i+)sitai0=0.0;/* 大循環(huán)開始 */for(i=2;i<M;i+)/* 求y(i+1) */y=DYi;/* 求x(i+1)的(4*1)的列陣 */for(k=0;k<1;k+)for(j=0;j<2;j+)xjk=-DYi+j-2;for(j=0;j<2;j+)xj+2k=DUi+j-2;temp=x0k;x0k=x1k;x1k=temp;temp=x2k;x2k=x3k;x3k=temp;/* 求gama(i+1)的值gama=1/(1+tx(i+1)*p(i)*x(i+1) */trans(x,tx,4,1); /* x(4 1) 到 tx(1 4) */MUL(tx,p,ji1,1,4,4); /* tx(1 4)*p(4 4)=ji1(1 4) */MUL(ji