系統(tǒng)辨識(shí)實(shí)驗(yàn)報(bào)告1

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

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

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

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

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

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

7、in函數(shù)中給出3個(gè)選擇：求u(k)；求y(k)；退出程序。其相對(duì)應(yīng)的函數(shù)名稱(chēng)為gener，ouput，exit。3. 在gener子函數(shù)中產(chǎn)生M序列u(k)保存到DU.txt文本文件中。4. 在output子函數(shù)中，通過(guò)對(duì)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)用來(lái)達(dá)到生成y(k)序列并保存到DY.txt文本文件中。四程序設(shè)計(jì)框圖：gener子函數(shù)產(chǎn)生DUmain函數(shù)input子函數(shù)輸入DU，DVoutput子函數(shù)產(chǎn)生DYdeal子函數(shù)得DY退出程序五程序運(yùn)

8、行流程圖：1. main函數(shù)：2. gener函數(shù)：3. output函數(shù)：實(shí)驗(yàn)3一實(shí)驗(yàn)內(nèi)容及要求：編寫(xiě)并調(diào)試動(dòng)態(tài)離散時(shí)間模型LS成批算法程序。要求：（1）原始數(shù)據(jù)由DU和DY讀出。（2）調(diào)用求逆及相乘子程序。（3）顯示參數(shù)辨識(shí)結(jié)果。二實(shí)驗(yàn)原理：1.批次處理的方法就是把所有的數(shù)據(jù)采集到一次性進(jìn)行處理，但前提是白色噪聲、及M序列所共同作用而產(chǎn)生的輸出，才能使用最小二乘法。雖然這種方法的計(jì)算量龐大，但經(jīng)常用于處理時(shí)不變系統(tǒng)，方法簡(jiǎn)單。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)方陣四實(shí)驗(yàn)運(yùn)行結(jié)果：實(shí)驗(yàn)4一實(shí)驗(yàn)內(nèi)容及要求：編寫(xiě)并調(diào)試動(dòng)態(tài)離散時(shí)間模型LS遞推算法程序。要求：（1）原始數(shù)據(jù)由DU和DY讀出。（2）顯示參數(shù)辨識(shí)結(jié)果。（3）設(shè)置選擇變量決定是否輸出中間結(jié)果。二實(shí)驗(yàn)原理：1.基本思路：新的估計(jì)值sita'(k+1)=老的估計(jì)值sita'(k)+修正項(xiàng)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ù)，為無(wú)窮小標(biāo)量。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)系圖：四實(shí)驗(yàn)運(yùn)行結(jié)果：1.main函數(shù)界面示意圖：6. 只顯示結(jié)果界面：7. 顯示過(guò)程的參考界面：心得體會(huì)我很喜歡這個(gè)課程的期末考核方式，不用再拘泥于在題目當(dāng)中對(duì)該課程的了解，而是通過(guò)4個(gè)C語(yǔ)言設(shè)計(jì)的練習(xí)來(lái)達(dá)到學(xué)習(xí)的目的，而且對(duì)以后的學(xué)習(xí)還有很大的幫助。在編寫(xiě)C語(yǔ)言的過(guò)程中，也遇到了一些阻礙，特別是在編寫(xiě)第3,4個(gè)的時(shí)候。比如：用的數(shù)組太多，并且未將其初始化，運(yùn)行出來(lái)的結(jié)果經(jīng)常是很長(zhǎng)的一段隨機(jī)數(shù)；或者一模一樣的程序有時(shí)候就可以正常運(yùn)行，有時(shí)候就總是出現(xiàn)報(bào)錯(cuò)這些都是讓我心塞了兩周的問(wèn)題。在這些問(wèn)題解決之后，

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

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

14、產(chǎn)設(shè)備，從而來(lái)達(dá)到最大的經(jīng)濟(jì)效益。這學(xué)期主要就是介紹方法，雖然可以處理一些特殊情況的有噪聲和沒(méi)有噪聲的系統(tǒng)，但對(duì)于那些噪聲模型都不可預(yù)測(cè)的系統(tǒng)中，還是顯得無(wú)能為力。希望以后還有機(jī)會(huì)學(xué)習(xí)其他的系統(tǒng)辨識(shí)方法。附錄（實(shí)驗(yàn)代碼）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;/* 對(duì)增廣矩陣進(jìn)行運(yùn)算（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+;/* 交換第一個(gè)數(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;/* 得到前一個(gè)矩陣的對(duì)角線上元素的乘積 */for(row=0;row<n;row+)value*=crowrow;/* 判斷乘積！若為零，則原矩陣沒(méi)有逆矩陣；反之亦然 */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;/* 得到標(biāo)準(zhǔn)的（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);/* 對(duì)P（ 4 * 4 ）和sita進(jìn)行初始化 */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)開(kāi)始 */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