系统辨识实验报告

1、精选优质文档-倾情为你奉上系统辨识实验报告 C语言版专心-专注-专业目录实验一1. 请编出矩阵A与B相乘得到矩阵R 的运算计算机程序。要求：(1) A和B的维数及数值可通过键盘及数据文件输入；(2)计算结果R可由屏幕及文件输出。2. 将题1改写为子程序。3. 查找有关的资料，读懂及调通矩阵求逆子程序，并改写为子程序。实验分析如下：1. 两个矩阵相乘，无非就是三个维数：m，n，p。因为第一个矩阵的纵向维度必须要跟第二个矩阵的横向维度相同才可以相乘。通过键盘输入维数的话用scanf语句即可，通过文件输入维数的话可以先存在一个文件里三个数值，到时候直接用只读的方式读出来。2. 改写成子程序也很简单，

2、我一开始就是写的子程序，注意不要漏掉函数声明。3. 矩阵求逆的程序我查阅了相关书籍，对照原理大体读懂了程序的编写，因为网上和书本上的程序还存在一定的问题，我反复调试后找到了问题并解决了，最终可以实现矩阵的求逆。实验代码如下：/*矩阵相乘*/#include<stdio.h>void MatrixMultiply(int m,int n,int p,long Matrix1MAXMAX,long Matrix2MAXMAX,long MatrixResultMAXMAX);void main()long Matrix1MAXMAX,Matrix2MAXMAX;long MatrixR

3、esultMAXMAX,Temp;int i,j,m,n,p;/*输入两个矩阵的的行列数m,n,p*/printf("nPlease input m of Matrix1:n");scanf("%d",&m);printf("Please input n of Matrix1:n");scanf("%d",&n);printf("Please input p of Matrix2:n");scanf("%d",&p);/*输入第一个矩阵的每个元素*/p

4、rintf("nPlease input elements of Matrix1(%d*%d):n",m,n);for(i=0;i<m;i+)for(j=0;j<n;j+)scanf("%ld",&Temp);Matrix1ij=Temp;/*输入第二个矩阵的每个元素*/printf("nPlease input elements of Matrix2(%d*%d):n",n,p);for(i=0;i<n;i+)for(j=0;j<p;j+)scanf("%ld",&Temp

5、);Matrix2ij=Temp;/*调用函数进行乘法运算，结果放在MatrixResult 中*/MatrixMultiply(m,n,p,Matrix1,Matrix2,MatrixResult);/*打印输出结果矩阵*/printf("nResult matrix: n");for(i=0;i<m;i+)for(j=0;j<p;j+)printf("%ld ",MatrixResultij);printf("n");void MatrixMultiply(int m,int n,int p,long Matrix1M

6、AXMAX,long Matrix2MAXMAX,long MatrixResultMAXMAX)/* 矩阵维数：mn*np=mp */int i,j,k;long Sum;/*嵌套循环计算结果矩阵（m*p）的每个元素*/for(i=0;i<m;i+)for(j=0;j<p;j+)/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0;for(k=0;k<n;k+)Sum+=Matrix1ik*Matrix2kj;MatrixResultij=Sum;/*矩阵求逆*/#include <stdlib.h>#include <math.h>#inc

7、lude <stdio.h>/*输入：待求逆的数组首地址，矩阵的维数输出：若为1，表示求逆成功，结果存在原矩阵中否则，表示求逆失败，输出err*not inv*/int inv(double a, int n) int *is,*js,i,j,k,l,u,v; double d,p;is=(int*)malloc(n*sizeof(int);/记忆行交换信息空间js=(int*)malloc(n*sizeof(int);/记忆列交换信息空间for (k=0; k<=n-1; k+) /全选主元 d=0.0;for (i=k; i<=n-1; i+)for (j=k; j

8、<=n-1; j+) l=i*n+j; p=fabs(al); if (p>d) d=p; isk=i; jsk=j;if (d+1.0=1.0) free(is); free(js); printf("err*not invn");return(0);if (isk!=k) /行交换for (j=0; j<=n-1; j+)u=k*n+j; v=isk*n+j;p=au; au=av; av=p;if (jsk!=k) /列交换for (i=0; i<=n-1; i+)u=i*n+k; v=i*n+jsk;p=au; au=av; av=p;l=k

9、*n+k;al=1.0/al;for (j=0; j<=n-1; j+) /归一化计算if (j!=k) u=k*n+j; au=au*al;for (i=0; i<=n-1; i+) /消元计算if (i!=k)for (j=0; j<=n-1; j+)if (j!=k) u=i*n+j;au=au-ai*n+k*ak*n+j;for (i=0; i<=n-1; i+) /消元计算if (i!=k) u=i*n+k; au=-au*al;for (k=n-1; k>=0; k-) /恢复if (jsk!=k)for (j=0; j<=n-1; j+) /行

10、交换 u=k*n+j; v=jsk*n+j;p=au; au=av; av=p;if (isk!=k) /列交换for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+isk;p=au; au=av; av=p;free(is); free(js);return(1); /正常返回void main() /主函数测试 int i,j;double a44= 21860.39,20893.46,12.70,-487.26,/待求逆的矩阵20893.46,21860.39,15.35,12.70,12.70,15.35,499.00,-5.00,-487.26,12.70,

11、-5.00,500.00;double b44; /存储原矩阵for (i=0; i<=3; i+) /将待求逆的矩阵内容存入b中for (j=0; j<=3; j+)bij=aij; i=inv(&a00,4); /调用函数进行求逆，结果存入a中if (i!=0) printf("MAT A IS:n");/打印原矩阵for (i=0; i<=3; i+) for (j=0; j<=3; j+)printf("%13.4f",bij);printf("n");printf("n")

12、;printf("MAT A- IS:n");for (i=0; i<=3; i+) /打印矩阵求逆的结果for (j=0; j<=3; j+)printf("%13.4f",aij);printf("n"); 运行结果如下：实验二编写并调试动态模型仿真程序：模型：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实验分析如下：通过上式可以看出，（

13、-1.5，0.7，1，0.5）T就是我们要求的系统模型参数值。 根据公式 =（XTX）-1XTY 知，要想求得，必须先知道Y和X的值。已知本实验的n为2，m为500。由上式可知要求出X的值，必须先知道数据y(k)和v(k)。而v(k)是白噪声，已知的。系统的输入信号是M序列，根据老师给的程序以及查找相关资料，我了解到M序列的求法就是一个移位寄存器的原理，7级的移位寄存器满足的算法原则是：x1(k+1)=x3kx6k(代表异或逻辑关系)。将老师给的程序进行适当翻译后，u(k)也就求出来了。然后通过给出的动态模型就可以求出y(k)了。我首先将老师给的DV（白噪声）文件读出，用matlab做了测试，

14、得到的结果如下：这是很典型的白噪声的自相关函数图形。而接下来求出的DU（M序列），它的自相关函数结果如下：从中可以看出，如果将DV应用于本实验的话，是符合LS算法使用的前提要求的。而DU也基本上符合M序列的定义，是理想的输入信号。到现在我们就求出了DU和DY，DV也已知，基本工作就都准备好了，可以进行LS成批算法的编写。实验代码如下：/*求DU（M序列）*/#include<stdio.h>#define N 500#define n 7 /级数void main(void)FILE *fp=NULL;int xn,i,j,s,uN;for(i=0;i<n;i+)xi=1;f

15、or(i=0;i<N;i+)s=x3x6; /核心算法for(j=6;j>0;j-)xj=xj-1;x0=s;if(x6=0)ui=-1;elseui=1;for(i=0;i<N;i+)printf("u%d=%3dn",i+1,ui);fp=fopen("E:DU.txt","w+"); /将算出的M序列写到文件里for(i=0;i<N;i+)fprintf(fp,"%3d",ui);/*求DY*/#include<stdio.h>#define L 500void main(

16、)double yL=0.0,uL=0.0,vL=0.0;int i,k;FILE *fp1=NULL,*fp3=NULL,*fp2=NULL;/*读取文件DU*/fp1=fopen("E:DU.txt","r+");for(i=0;i<L;i+)fscanf(fp1,"%lf",(u+i);if(fgetc(fp1)=EOF)break;fclose(fp1);/*读取文件DV*/fp2=fopen("E:DV.txt","r+");for(i=0;i<L;i+)fscanf(fp

17、2,"%lf",(v+i);if(fgetc(fp2)=EOF)break;fclose(fp2);/*计算DY并写入文件里*/for(k=2;k<L;k+)yk=1.5*yk-1-0.7*yk-2+uk-1+0.5*uk-2+vi;fp3=fopen("E:DY.txt","w+");fprintf(fp3,"nn-输出序列DY-n");for(i=0;i<L;i+)fprintf(fp3,"%10.4f",yi);fclose(fp3);printf("-y-n&quo

18、t;);for(i=0;i<200;i+)printf("%10.4f",yi);if(i+1)%5=0)printf("n");printf("n");运行结果如下：其它详见附录。实验三编写并调试动态离散时间模型LS成批算法程序。要求：(1)原始数据由DU和DY读出；(2) 调用求逆及相乘子程序；(3)显示参数辨识结果。实验分析如下：由前面两个实验的准备，只要将它们联系起来，就可以实现LS成批算法了。LS成批算法流程如下：实验代码如下：#include<stdio.h>#include <stdlib.h&g

19、t;#include <math.h>#define L 500/*声明函数*/void MatrixMultiply1(int m,int n,int p,double Matrix14498, double Matrix24984,double MatrixResult44);void MatrixMultiply2(int m,int n,int p,double Matrix144, double Matrix24498,double MatrixResult4498);void MatrixMultiply3(int m,int n,int p,double Matrix

20、14498, double Matrix24981,double MatrixResult41);int inv(double a, int n);/*/void main()double yL=0.0,uL=0.0,vL=0.0,x4984;int i,k,j;FILE *fp1=NULL,*fp2=NULL,*fp3=NULL;/*读取文件u*/fp1=fopen("E:DU.txt","r+");for(i=0;i<L;i+)fscanf(fp1,"%lf",(u+i);if(fgetc(fp1)=EOF)break;fc

21、lose(fp1);/*读取文件v*/fp2=fopen("E:DV.txt","r+");for(i=0;i<L;i+)fscanf(fp2,"%lf",(v+i);if(fgetc(fp2)=EOF)break;fclose(fp2);/*计算DY*/for(k=2;k<L;k+)yk=1.5*yk-1-0.7*yk-2+uk-1+0.5*uk-2+vi;fp3=fopen("E:DY.txt","w+");for(i=0;i<L;i+)fprintf(fp3,"

22、%lfn",yi);fclose(fp3);/*计算X*/for(i=0;i<498;i+)xi0=-yi+1;xi1=-yi;xi2=ui+1;xi3=ui;/*计算最终结果*/double temp4498,temp24498,temp341;int m=0;double xtx44;double y14981=0;for(i=0;i<498;i+) for(j=0;j<4;j+) tempji=xij; /将x矩阵转置for(i=0;i<498;i+)y1i0=yi+2; /将前面算出的DY转化成498行，1列的矩阵/调用函数，求系数矩阵MatrixMu

23、ltiply1(4,498,4,temp,x,xtx);inv(&xtx00,4);MatrixMultiply2(4,4,498,xtx,temp,temp2);MatrixMultiply3(4,498,1,temp2,y1,temp3);/*最终结果输出*/ printf("a1=%6.4fna2=%6.4fnb1=%6.4fnb2=%6.4fn",temp300,temp310,temp320,temp330);/*函数部分*/void MatrixMultiply1(int m,int n,int p,double Matrix14498, double

24、Matrix24984,double MatrixResult44) /* 矩阵维数：mn*np=mp */int i,j,k;double Sum;/*嵌套循环计算结果矩阵（m*p）的每个元素*/for(i=0;i<m;i+)for(j=0;j<p;j+)/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0.0;for(k=0;k<n;k+)Sum+=Matrix1ik*Matrix2kj;MatrixResultij=Sum;void MatrixMultiply2(int m,int n,int p,double Matrix144, double Matrix

25、24498,double MatrixResult4498) /* 矩阵维数：mn*np=mp */int i,j,k;double Sum;/*嵌套循环计算结果矩阵（m*p）的每个元素*/for(i=0;i<m;i+)for(j=0;j<p;j+)/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0.0;for(k=0;k<n;k+)Sum+=Matrix1ik*Matrix2kj;MatrixResultij=Sum;void MatrixMultiply3(int m,int n,int p,double Matrix14498, double Matrix24

26、981,double MatrixResult41) /* 矩阵维数：mn*np=mp */int i,j,k;double Sum;/*嵌套循环计算结果矩阵（m*p）的每个元素*/for(i=0;i<m;i+)for(j=0;j<p;j+)/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0.0;for(k=0;k<n;k+)Sum+=Matrix1ik*Matrix2kj;MatrixResultij=Sum;int inv(double a, int n) int *is,*js,i,j,k,l,u,v; double d,p;is=(int*)malloc(n

27、*sizeof(int);/记忆行交换信息空间js=(int*)malloc(n*sizeof(int);/记忆列交换信息空间for (k=0; k<=n-1; k+) /全选主元 d=0.0;for (i=k; i<=n-1; i+)for (j=k; j<=n-1; j+) l=i*n+j; p=fabs(al); if (p>d) d=p; isk=i; jsk=j;if (d+1.0=1.0) free(is); free(js); printf("err*not invn");return(0);if (isk!=k) /行交换for (j

28、=0; j<=n-1; j+)u=k*n+j; v=isk*n+j;p=au; au=av; av=p;if (jsk!=k) /列交换for (i=0; i<=n-1; i+)u=i*n+k; v=i*n+jsk;p=au; au=av; av=p;l=k*n+k;al=1.0/al;for (j=0; j<=n-1; j+) /归一化计算if (j!=k) u=k*n+j; au=au*al;for (i=0; i<=n-1; i+) /消元计算if (i!=k)for (j=0; j<=n-1; j+)if (j!=k) u=i*n+j;au=au-ai*n

29、+k*ak*n+j;for (i=0; i<=n-1; i+) /消元计算if (i!=k) u=i*n+k; au=-au*al;for (k=n-1; k>=0; k-) /恢复 if (jsk!=k)for (j=0; j<=n-1; j+) /行交换 u=k*n+j; v=jsk*n+j;p=au; au=av; av=p;if (isk!=k) /列交换for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+isk;p=au; au=av; av=p;free(is); free(js);return(1); /正常返回运行结果如下：实验四编

30、写并调试动态离散时间模型LS递推算法程序。要求：(1)原始数据由DU和DY读出；(2)显示辨识结果；(3)设置选择变量决定是否输出中间结果。实验分析如下：根据老师给的流程图，一步步编写调试。LS递推算法流程如下：实验代码如下：/*u为输入变量，y为输出变量，theta为系统参数矩阵，p,gama,x为RLS算法中要用到的参数,xT为x的转置矩阵。z1,z2,z3为运算的中间变量。z变量用以判断是否需要输出中间结果，若为'y',则输出中间结果，否则，不输出中间结果。*/#include<stdio.h>#define L 500#define MAXSIZE 50vo

31、id main()void mulmat(double m1MAXSIZEMAXSIZE,double m2MAXSIZEMAXSIZE,double m3MAXSIZEMAXSIZE,int H1,int H2,int L1,int L2);double uL,yL,thetaMAXSIZEMAXSIZE=0,pMAXSIZEMAXSIZE,gama,xMAXSIZEMAXSIZE,xTMAXSIZEMAXSIZE;double z1MAXSIZEMAXSIZE=0,z2MAXSIZEMAXSIZE=0,z3MAXSIZEMAXSIZE=0;double alpha=1e8;int i,j,

32、k;char z;FILE *fp1=NULL,*fp2=NULL;/判断是否需要输出中间结果printf("是否需要输出中间结果？y/n:");scanf("%c",&z);printf("n");if(z='y')printf("a1a2b1b2n");/读入DU和DY数据fp1=fopen("E:DU.txt","r+");for(i=0;i<L;i+)fscanf(fp1,"%lf",(u+i);if(fgetc(fp

33、1)=EOF)break;fclose(fp1);fp2=fopen("E:DY.txt","r+");for(i=0;i<L;i+)fscanf(fp2,"%lf",(y+i);if(fgetc(fp2)=EOF)break;fclose(fp2);/p(0)和theta(0)赋初值for(i=0;i<4;i+)for(j=0;j<4;j+)if(i=j)pij=alpha;elsepij=0;thetai0=0;for(k=1;k<(L-1);k+)if(z='y') for(i=0;i&l

34、t;4;i+)printf("%4.2ftt",thetai0);printf("n");/测量yk+1,形成xk+1x00=-yk+1-1;x10=-yk+1-2;x20=uk+1-1;x30=uk+1-2;for(i=0;i<4;i+)xT0i=xi0;mulmat(xT,p,z1,1,4,4,4);mulmat(z1,x,z2,1,4,4,1);gama=1/(1+z200);mulmat(p,x,z1,4,4,4,1);mulmat(xT,theta,z2,1,4,4,1);z200=yk+1-z200;for(i=0;i<4;i+)

35、thetai0=thetai0+gama*z200*z1i0;mulmat(p,x,z1,4,4,4,1);mulmat(z1,xT,z2,4,1,1,4);mulmat(z2,p,z3,4,4,4,4);for(i=0;i<4;i+)for(j=0;j<4;j+)pij=pij-gama*z3ij;printf("a1=%4.2f,a2=%4.2f,b1=%4.2f,b2=%4.2fn",theta00,theta10,theta20,theta30);/矩阵相乘子函数void mulmat(double m1MAXSIZEMAXSIZE,double m2M

36、AXSIZEMAXSIZE,double m3MAXSIZEMAXSIZE,int H1,int H2,int L1,int L2)int i,j,k,H3,L3;H3=H1; L3=L2;for(i=0;i<H3;i+)for(j=0;j<L3;j+)m3ij=0;for(i=0;i<H3;i+)for(j=0;j<L3;j+)for(k=0;k<L1;k+)m3ij+=m1ik*m2kj;运行结果如下：实验总结：首先很开心程序成功运行出了结果，并且与准确值相比较偏差不大。在编写调试的过程中极大地锻炼了自己的编程能力，也在一次次错误中学会了冷静分析。因为是跟同学

37、一起做的，通过讨论及时地改正了自己理解的错误，也是程序成功运行地重要因素。实验过程中不仅锻炼我的C语言编程能力，还锻炼了我的matlab编程能力。当用C语言编写出问题，数值不太正确时，都会用matlab验证一下。尤其是求M序列时，验证了它的自相关函数。从而生成了一个比较理想的M序列信号。第三个实验中遇到了一些问题，我运用matlab求出的结果进行对照，最终发现是逆矩阵求出的结果有错误。经过仔细地研究发现是实参传递错了，改正后，运行出想要的代码。附录附录1：DU(M序列)-输出序列DU- 1 1 1 1 1 1 -1 -1 -1 -1 1 1 1 -1 1 1 1 1 -1 -1 1 -1 1

38、1 -1 -1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 1 1 -1 1 1 1 -1 -1 -1 1 1 1 1 1 1 1 -1 -1 -1 -1 1 1 1 -1 1 1

39、1 1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 1 1 -1 1 1 1 -1 -1 -1 1 1 1 1 1 1 1 -1 -1 -1

40、 -1 1 1 1 -1 1 1 1 1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 1 1 -1 1 1 1 -1 -1 -1 1 1 1

41、 1 1 1 1 -1 -1 -1 -1 1 1 1 -1 1 1 1 1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 1 1附录2：DV

42、3. 3. 2. -5.E-001 1.E-001 9.E-002 7.E-002 3.E-001 -5.E-001 9.E-002 -3.E-001 -4.E-002 -4.E-001 7.E-001 -4.E-001 -1.E-001 7.E-001 7.E-001 -7.E-002 1.E-001 1.E-002 -8.E-001 -2.E-001 1. 3.E-001 4.E-002 -1.E-001 2.E-001 5.E-002 -7.E-001 1. -4.E-001 6.E-001 -3.E-001 -1.E-001 4.E-001 4.E-001 4.E-001 1.E-00

43、1 -1.E-001 -1. 1.E-001 3.E-001 -5.E-001 1. 6.E-001 -9.E-001 6.E-001 -4.E-001 5.E-001 -9.E-001 -9.E-001 2.E-002 2.E-001 -8.E-001 1.E-001 7.E-001 -1.E-001 4.E-001 -5.E-001 -3.E-001 6.E-001 1.E-001 -1. -3.E-001 2.E-001 -6.E-001 -2.E-001 5.E-001 7.E-001 7.E-001 -3.E-002 -4.E-001 -7.E-001 -5.E-001 -4.E-0

44、01 1. 1.E-001 1. 1. 2.E-001 -8.E-001 -4.E-001 1. 7.E-001 -3.E-001 4.E-001 -7.E-002 9.E-001 6.E-001 4.E-001 8.E-001 9.E-002 -8.E-001 3.E-001 -1.E-001 6.E-001 8.E-001 3.E-001 9.E-001 -7.E-001 -1. -1.E-001 -1.E-001 -2.E-001 9.E-001 -3.E-001 6.E-001 1. -1.E-001 -5.E-001 2.E-001 1.E-002 2.E-001 2.E-001 3.E-001 -5.E-001 -4.E-001 -3.E-001 8.E-001 -3.E-001 1. 8.E-001 -2.E-001 -6.E-001 -3.E-001 3.E-001 -1. -5.E-001