系统辨识实验报告

1、系统辨识实验报告学院：信息科学与技术学院专业：自动化日期：2016/4/26目录实验14一实验内容及要求：4二实验原理：4三软件设计思想：4四程序结构框图：5五运行示意图：5实验28一实验内容及要求：8二实验原理：8三软件设计思想：9四程序设计框图：9五程序运行流程图：10实验312一实验内容及要求：12二实验原理：12三程序数据流程图：12四实验运行结果：13实验414一实验内容及要求：14二实验原理：14三数据递推关系图：14四实验运行结果：15心得体会16附录（实验代码）171.LabWork1172.LabWork2213.LabWork3234.LabWork426实验1一实验内容及

2、要求：1.编出矩阵A与B相乘得到的矩阵R的运算计算机程序要求：（1）A和B的维数及数值可通过键盘及数据文件输入（2）计算结果R可由屏幕及文件输出2.将1改写为子程序3查找有关的资料，读懂及调通矩阵求逆程序，并改写为子程序。二实验原理：1. 两个矩阵A、B相乘得到C矩阵，首先要满足的条件是A的列数与B行数相等，否则不能相乘。当满足条件后，根据C（i,k）=Ai,j*B(j,k)可以求得C矩阵。2. 当求矩阵的逆时，首先要判断其是否为方阵，若是则可以对其进行下一步的操作。本次实验中求逆主要是通过构造一个增广矩阵（FangZ | E）矩阵的初等行变换得到（E | FZNi）的这样的一个矩阵就可以求得

3、矩阵的逆。若矩阵FangZ不是满秩矩阵时，FangZ没有FZNi 。通过这样的求逆方式，避免了大方阵的求取行列式运算。三软件设计思想：1. 确定该软件的功能主要有：键盘输入两个矩阵然后相乘；文本data输入两个矩阵将结果放在文本result中；键盘输入一个方阵求得其逆矩阵。其中前两个的矩阵相乘运算部分设置为一个函数Mul。2. 在main函数中提供两个关于矩阵的选择：multiplitation；invertion。其相对应的子函数为MulOp(a,b,c)，Inv()。3. 在MulOp(a,b,c)子函数中，有两种输入矩阵的方式：way1，way2。相对应的功能为键盘输入，文本输入。并且两

4、者在处理矩阵时，都调用了Mul函数。4. 在Inv()子函数中，输入和显示原矩阵，和其相应的逆矩阵。调用qiuni（double FangZM,double FZNiM,int n）子函数，可以都得到原矩阵的逆矩阵。但当原矩阵不可逆时，系统输出为”The array is not invertible!“。四程序结构框图：way1键盘输入矩阵Mul矩阵乘积子函数MulOp子函数矩阵乘积way2文本输入矩阵main函数Inv子函数矩阵求逆qiuni子函数退出程序五运行示意图：1. main函数的主界面：2. MulOp子函数的界面：3. Inv子函数的界面：4. 通过键盘操作计算两个矩阵的乘积：

5、5. 求方阵的逆矩阵：实验2一实验内容及要求：编写并调试动态模型仿真程序：模型：y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)已知：白噪声v(k)数据文件为DV，数据长度为L=500要求：（1） 产生长度为L的M序列数据文件DU （2） 产生长度为L的模型输出数据文件DY 二实验原理：由于在现实中，白噪声序列很难求，所以寻找到M序列在一定程度上可以代替白色噪声序列。由L=500，所以n=9。根据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级未参加反馈。系统产生M序列的结构流程图：初始化寄存器a=1,1,0,0,0,1,0,1,0移位寄存器c=0,0,0,1,0,0,0,0,1M序列的长度为L=500i=1Y输出M序列iLNM序列DU(i)=a(i)寄存器向前移1位an=i=0n-1ciai(mod 2)三软件设计思想：1. 该软件的主要功能是：产生M序列赋给u(k)保存在DU.txt文件中；由u(k)和v(k)求得y(k)保存在DY.txt文件中。2. 在main函数

7、中给出3个选择：求u(k)；求y(k)；退出程序。其相对应的函数名称为gener，ouput，exit。3. 在gener子函数中产生M序列u(k)保存到DU.txt文本文件中。4. 在output子函数中，通过对input子函数（读入v(k),u(k)的数据）、deal子函数(由公式y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)求y(k)的调用来达到生成y(k)序列并保存到DY.txt文本文件中。四程序设计框图：gener子函数产生DUmain函数input子函数输入DU，DVoutput子函数产生DYdeal子函数得DY退出程序五程序运行流程图

8、：1. main函数：2. gener函数：3. output函数：实验3一实验内容及要求：编写并调试动态离散时间模型LS成批算法程序。要求：（1）原始数据由DU和DY读出。（2）调用求逆及相乘子程序。（3）显示参数辨识结果。二实验原理：1.批次处理的方法就是把所有的数据采集到一次性进行处理，但前提是白色噪声、及M序列所共同作用而产生的输出，才能使用最小二乘法。虽然这种方法的计算量庞大，但经常用于处理时不变系统，方法简单。2.构造模型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)-y(4

9、98)-y(497)u(2)u(1)u(3)u(2)u(498)u(497)三程序数据流程图：转置Y(m,1)列矩阵TX(4,m)矩阵X(m,4)矩阵xishu(4,1)列矩阵相乘相乘ji2(4,m)矩阵求逆相乘ni(4,4)方阵ji1(4,4)方阵四实验运行结果：实验4一实验内容及要求：编写并调试动态离散时间模型LS递推算法程序。要求：（1）原始数据由DU和DY读出。（2）显示参数辨识结果。（3）设置选择变量决定是否输出中间结果。二实验原理：1.基本思路：新的估计值sita(k+1)=老的估计值sita(k)+修正项1. 基本模型：y(k)=-a1y(k-1)-a2y(k-2)+b1u(k-

10、1)+b2u(k-2)+v(k)。2. 引入信息矩阵P(k)，维数为4*4，初始化为1051012*E的一单位阵。gama(k)为修正系数，为无穷小标量。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. gama(k+1)=1/1+x(k+1)*P(k)*x(k+1);三数据递推关系图：四

11、实验运行结果：1.main函数界面示意图：6. 只显示结果界面：7. 显示过程的参考界面：心得体会我很喜欢这个课程的期末考核方式，不用再拘泥于在题目当中对该课程的了解，而是通过4个C语言设计的练习来达到学习的目的，而且对以后的学习还有很大的帮助。在编写C语言的过程中，也遇到了一些阻碍，特别是在编写第3,4个的时候。比如：用的数组太多，并且未将其初始化，运行出来的结果经常是很长的一段随机数；或者一模一样的程序有时候就可以正常运行，有时候就总是出现报错这些都是让我心塞了两周的问题。在这些问题解决之后，运行出来的结果却与实际模型参数的出入有点大，于是又重新查找第2个实验是否是M序列产生的方式有问题。

12、通过对初始化寄存器赋给不同的值，可以让结果与真实模型参数之间的误差达到最小。通过这学期对系统辨识这个课程的学习，让我了解到系统辨识在建立数学模型的方面的重要性。对于不了解系统的工作机理的人来说，也有了一个可以知道系统模型的构建方法。但是也必须要我们对所构建模型有一个较为清醒的认识，比如要知道模型的类型，如果是动态模型，则要知道模型的参数或者阶次（但是权函数模型就不需要事前知道模型的参数或者阶次，不知道这种方法有没有什么我们不知道的缺陷）。在处理数据上，又一次让我了解到世上没有两全其美的事情，计算的复杂程度和精度就好像鱼和熊掌不能兼得。也正是因为这样，才会成为促使人们在这方面的不断寻找最优化算法

13、的动力。在对模型进行研究时，都是从最简单的模型开始研究，比如在没有噪声的情况下，所得到的模型参数就为系统模型真实参数。进一步，在白色噪声的情况下，所得到的参数就为系统模型真实参数的估计通过这样的推理，就可以得到LS在系统辨识问题上的普适性，不管系统受怎样的干扰，只要通过一系列的变换，最后还是通过LS来解决模型参数的确定工作。虽然我们对化学工艺不怎么了解，但通过学习系统辨识，让我们对已经投运后的工业生产设备的性能有一个跟踪性的了解，可以对系统做出一个较为完整的评估，从而来改善生产工艺或者是及时更换生产设备，从而来达到最大的经济效益。这学期主要就是介绍方法，虽然可以处理一些特殊情况的有噪声和没有噪

14、声的系统，但对于那些噪声模型都不可预测的系统中，还是显得无能为力。希望以后还有机会学习其他的系统辨识方法。附录（实验代码）1. LabWork1#define M 100double FangZMM=0,FZNiMM=0;/* 求矩阵的乘法函数 */void Mul(double aM,double bM,double cM,int n1,int n2,int n4)int i,j,k;for(i=0;in1;i+)for(k=0;kn4;k+)for(j=0;jn2;j+)cik+=aij*bjk;/* 求方阵的逆 */void qiuni(double FangZM,double FZNi

15、M,int n)double EMM,value=1.0,stemp,temp=0.0,san,cMM;int i,j,row,nextrow,flag=0;int col,switchtime=0;/* 构造原矩阵的增广矩阵 */for(i=0;in;i+)for(j=0;jn;j+)if(j=i) Eij=1;else Eij=0;for(i=0;in;i+)for(j=0;jn;j+)cij=FangZij;for(j=0;j（E|FZNi） */for(row=0;rown-1;row+)nextrow=row+1;if(cnextrowrow=0)while(cnextrowrow=

16、0)nextrow+;if(nextrow=n)flag=1;break;if(flag=1) continue;switchtime+;/* 交换第一个数字为0的行 */for(col=0;col2*n;col+)stemp=crowcol;crowcol=cnextrowcol;cnextrowcol=stemp;for(nextrow=row+1;nextrown;nextrow+)temp=cnextrowrow*1.0/crowrow;for(col=0;col2*n;col+)cnextrowcol+=-temp*crowcol;/* 得到前一个矩阵的对角线上元素的乘积 */for

17、(row=0;row=1;row-)for(i=row-1;i=0;i-) san=cirow*1.0/crowrow;for(col=row;col2*n;col+)cicol-=crowcol*san;/* 得到标准的（E|A） */for(row=0;rown;row+)for(col=0;coln;col+)crowcol+n/=crowrow;for(i=0;in;i+)for(j=0;jn;j+)FZNiij=cij+n;2. LabWork2/* 处理v(k),u(k)的数据得输出y(k) */void deal(double DVM,double DUM,double DYM)

18、int k;DY0=0;DY1=0;for(k=2;k500;k+)DYk=1.5*DYk-1-0.7*DYk-2+DUk-1+0.5*DUk-2+DVk;/* 产生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 open the file!n);exit(0);for(i=1;iM;i+)sum=0;for(j

19、=0;j9;j+)sum+=aj*cj;aa0=sum%2;for(j=1;j9;j+)aaj=aj-1;for(j=0;j9;j+)aj=aaj;DUi=a8;for(i=0;iM;i+)fprintf(fp,%d ,DUi);fclose(fp);printf(nnnttPlease refer to the file(DU.txt)!nn);system(pause);3. LabWork3#define M 500#define N 4void MUL1(double TXM,double XN,double ji1N)int i,j,k;for(i=0;iN;i+)for(k=0;k

20、N;k+)for(j=0;jM-3;j+)ji1ik+=TXij*Xjk;void MUL2(double niN,double TXM,double ji2M)int i,j,k;for(i=0;iN;i+)for(k=0;kM-3;k+)for(j=0;jN;j+)ji2ik+=niij*TXjk;void MUL3(double ji2M,double YM,double xishuN)int i,j;for(i=0;iN;i+)for(j=0;jM-3;j+)xishui+=ji2ij*Yj;void main()system(cls);system(color 75);int i,j

21、;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;iM-3;i+)for(j=0;j2;j+)Xij=-DYi+j+1;for(j=0;j2;j+)Xij+2=DUi+j+1;temp=0;temp=Xi0;Xi0=Xi1;Xi1=temp;temp=0;temp=Xi2;Xi2=Xi3;Xi3=temp;/* 得到Y（498 * 1）=DY(1) DY(498)的列矩阵 */for(i=0;iM-3;i+)Yi

22、=DYi+3;/* 得到X（498 * 4）的转置矩阵 TX（4 * 498） */for(i=0;iN;i+)for(j=0;jM-3;j+)TXij=Xji; MUL1(TX,X,ji1); /* 得到TX*X的乘积ji1 */qiuni(ji1,ni,N); /* 得到ji1的逆矩阵ni */MUL2(ni,TX,ji2); /* 得到ni*TX的乘积ji2 */MUL3(ji2,Y,xishu); /* 得到ji2*Y的乘积xishu */printf(tthe xishu are:nn);for(i=0;iN;i+)printf(nnt %.4f n,xishui);printf(n

23、n);4. LabWork4#define M 500#define N 4double sitaNN;void MUL(double aN,double bN,double cN,int n1,int n2,int n4)int i,j,k;for(i=0;in1;i+)for(k=0;kn4;k+)for(j=0;jn2;j+)cik+=aij*bjk;void trans(double dN,double eN,int n5,int n6)int i,j;for(i=0;in6;i+)for(j=0;jn5;j+)eij=dji;void consult()system(cls);int

24、 i,j,k;double DUM,DYM;double txNN=0,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;iN;i+)for(j=0;jN;j+)if(i=j) pij=10e9;else pij=0;for(i=0;iN;i+)sitai0=0.0;/* 大循环开始 */for(i=2;iM;i+)/* 求y(i+1) */y=DYi;/* 求x(i+1)的(4*1)的列阵 */for(k=0;k1;k+)for(j=0;j2;j+)xjk=-DYi+j-2;for(j=0;j2;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)