系统辨识实验报告

1、实验一 矩阵相乘与求逆1、 实验要求1） 请编出矩阵A和B相乘得到矩阵R的运算计算机程序，要求A和B的维数及数值可以通过键盘及数据文件输入，计算结果R可由屏幕及文件输出。2） 将题一的程序改为子程序。3） 查找有关的资料，读懂及调通矩阵求逆子程序，并改写为子程序。2、 实验过程1） 键盘输入矩阵A和B，屏幕上显示计算结果矩阵R。#include#include int main() int i,j,k,m,n,p,q; printf(请输入A矩阵的行数：); scanf(%d,&m); printf(请输入A矩阵的列数：); scanf(%d,&n); printf(请输入B矩阵的行数：);

2、scanf(%d,&p); printf(请输入B矩阵的列数：); scanf(%d,&q); if(n!=p) printf(请检查A矩阵的列数与B矩阵的行数是否一致！); else double *matrixA,*matrixB,*matrixC; matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)malloc(sizeof(double)*m*n); matrixB=(double *)malloc(sizeof(double *)*p); matrixB0=(double *)malloc(sizeof(

3、double)*p*q); matrixC=(double *)malloc(sizeof(double *)*m); matrixC0=(double *)malloc(sizeof(double)*m*q); for(i=1;im;i+) matrixAi=matrixAi-1+n; for(i=1;ip;i+) matrixBi=matrixBi-1+q; for(i=1;im;i+) matrixCi=matrixCi-1+q; printf(请输入A矩阵的数：n); for(i=0;im;i+) for(j=0;jn;j+) scanf(%lf,&matrixAij); printf

4、(n输入的A矩阵是：n); for(i=0;im;i+) for(j=0;jn;j+) printf(%7.3lft,matrixAij); printf(n); printf(请输入B矩阵的数：n); for(i=0;ip;i+) for(j=0;jq;j+) scanf(%lf,&matrixBij); printf(n输入的B矩阵是：n); for(i=0;ip;i+) for(j=0;jq;j+) printf(%7.3lft,matrixBij); printf(n); for(i=0;im;i+) for(j=0;jq;j+) matrixCij=0; for(i=0;im;i+)

5、 for(j=0;jq;j+) for(k=0;kn;k+) matrixCij=matrixCij+matrixAik*matrixBkj; printf(n矩阵A*B的结果是：n); for(i=0;im;i+) for(j=0;jq;j+) printf(%7.3lft,matrixCij); printf(n); free (matrixA); free (matrixB); free (matrixC); getchar(); getchar();运行的结果如下：2） 实现文本的矩阵输入与计算结果的输出，并将相乘改为子程序。#includestdio.hint main() void

6、 matrix_mutipl(double *matrixA,double *matrixB,double *matrixC,int m,int n,int q); char ch;int m=0,n=0,p=0,q=0,i=0,j=0,k=0;char infileA80,infileB80,outfile80;FILE *fpin1,*fpin2,*fpin3;double *matrixA,*matrixB,*matrixC;printf(Enter the infileA name:n); scanf(%s,infileA); if(fpin1=fopen(infileA,r)=NUL

7、L) printf(Cant open file %sn,infileA); exit(0); do ch=fgetc(fpin1); if(ch=t) n+; if(ch=n) m+; while(!feof(fpin1); n=(n+m)/m; printf(输入的矩阵A行与列数分别为：%d,%dn,m,n); fclose(fpin1); fflush(stdin); matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)malloc(sizeof(double)*m*n); for(i=1;im;i+) mat

8、rixAi=matrixAi-1+n; for(i=0;im;i+) for(j=0;jn;j+) matrixAij=0; if(fpin1=fopen(infileA,r)=NULL) printf(Cant open file %sn,infileA); exit(0); do ch=fgetc(fpin1); while(!feof(fpin1)&(ch9); if (feof(fpin1) exit(0); ungetc(ch,fpin1); do i=k/n; j=k%n; fscanf(fpin1,%lf,&matrixAij); k+; while(!feof(fpin1)&k

9、m*n); fclose(fpin1); printf(输入的A矩阵为：n); for(i=0;im;i+) for(j=0;jn;j+) printf(%lf,matrixAij); printf(t); printf(n); fflush(stdin); printf(Enter the infileB name:n); scanf(%s,infileB); if(fpin2=fopen(infileB,r)=NULL) printf(Cant open file %sn,infileB); exit(0); do ch=fgetc(fpin2); if(ch=t) q+; if(ch=n

10、) p+; while(!feof(fpin2); q=(p+q)/p; printf(输入的矩阵B行与列数分别为：%d,%dn,p,q); fclose(fpin2); fflush(stdin); matrixB=(double *)malloc(sizeof(double *)*p); matrixB0=(double *)malloc(sizeof(double)*p*q); for(i=1;ip;i+) matrixBi=matrixBi-1+q; for(i=0;ip;i+) for(j=0;jq;j+) matrixBij=0; if(fpin2=fopen(infileB,r)

11、=NULL) printf(Cant open file %sn,infileB); exit(0); do ch=fgetc(fpin2); while(!feof(fpin2)&(ch9); if (feof(fpin2) exit(0); ungetc(ch,fpin2); k=0; do i=k/q; j=k%q; fscanf(fpin2,%lf,&matrixBij); k+; while(!feof(fpin2)&kp*q); fclose(fpin2); printf(输入的B矩阵为：n); for(i=0;ip;i+) for(j=0;jq;j+) printf(%lf,ma

12、trixBij); printf(t); printf(n); fflush(stdin); matrixC=(double *)malloc(sizeof(double *)*m); matrixC0=(double *)malloc(sizeof(double)*m*q); for(i=1;im;i+) matrixCi=matrixCi-1+q; for(i=0;im;i+) for(j=0;jq;j+) matrixCij=0; matrix_mutipl(double*)matrixA,(double*)matrixB,(double*)matrixC,m,n,q); printf(

13、矩阵A*B为：n); for(i=0;im;i+) for(j=0;jq;j+) printf(%lf,*(*(matrixC+i)+j); printf(t); printf(n); fflush(stdin); printf(Enter the outfile name:n); scanf(%s,outfile); if(fpin3=fopen(outfile,w)=NULL) printf(Cant open file %sn,outfile); exit(0); for(i=0;im;i+) for(j=0;jq;j+) fprintf(fpin3,%lf,*(*(matrixC+i)

14、+j); fprintf(fpin3,t); fprintf(fpin3,n); fflush(stdin); getchar(); void matrix_mutipl(double *matrixA,double *matrixB,double *matrixC,int m,int n,int q) int k=0,j=0,i=0; for(i=0;im;i+) for(j=0;jq;j+) for(k=0;kn;k+) *(*(matrixC+i)+j)=*(*(matrixC+i)+j)+*(*(matrixA+i)+k)*(*(*(matrixB+k)+j);运行结果如下： 文本文件

15、1.dat和2.dat需要放在程序的同一文件夹下，否则需要说明路径。3） 矩阵求逆子程序#includestdio.h#includemalloc.h#includemath.h void main() int inv(double *p,int n); double a44=1,2,0,0,1,2,0,0,0,0,3,0,0,0,0,1,*ab; ab=a0; int n=4,i=0,j; i=inv(ab,n); if(i!=0) for(i=0;in;i+) putchar(n); for(j=0;jn;j+) printf(%f ,aij); getchar();int inv(dou

16、ble *p,int n) void swap(double *a,double *b); int *is,*js,i,j,k,l; for(i=0;in;i+) putchar(n); for(j=0;jn;j+) printf(%f ,*(p+i*n+j); puts(nnnn); double temp,fmax; is=(int *)malloc(n*sizeof(int); js=(int *)malloc(n*sizeof(int); for(k=0;kn;k+) fmax=0.0; for(i=k;in;i+) for(j=k;jfmax) fmax=temp; isk=i;js

17、k=j; if(fmax+1.0)=1.0) free(is);free(js); printf(no inv); return(0); if(i=isk)!=k) for(j=0;jn;j+) swap(p+k*n+j,p+i*n+j); if(j=jsk)!=k) for(i=0;in;i+) swap(p+i*n+k,p+i*n+j); pk*n+k=1.0/pk*n+k; for(j=0;jn;j+) if(j!=k) pk*n+j*=pk*n+k; for(i=0;in;i+) if(i!=k) for(j=0;jn;j+) if(j!=k) pi*n+j=pi*n+j-pi*n+k

18、*pk*n+j; for(i=0;i=0;k-) if(j=jsk)!=k) for(i=0;in;i+) swap(p+j*n+i),(p+k*n+i); if(i=isk)!=k) for(j=0;jn;j+) swap(p+j*n+i),(p+j*n+k); free(is); free(js); return 1;void swap(double *a,double *b) double c; c=*a; *a=*b; *b=c;如上所示，本子程序先将输入的矩阵进行全对角占优的方法，避免了矩阵处理过程中的“病态”情况。3、 实验总结在整个实验过程中，出现了三个问题：1） 文本输入是实现

19、矩阵相乘的重要点，由于本人在实验的文本文件中加了一些提示语言，因此需要判断何时是数据，再进行矩阵数据的读入，因此增加了一个数据判断过程。2） 矩阵程序改为子程序时，需要实现地址之间的传递，我们不能简单地进行matrixAij处理，应该改为*(*(matrixA+i)+j)进行处理，而且还要实现二维数组的动态分配存储空间，如下：double *matrixA;matrixA=(double *)malloc(sizeof(double *)*m); matrixA0=(double *)malloc(sizeof(double)*m*n); for(i=1;im;i+) matrixAi=mat

20、rixAi-1+n; for(i=0;im;i+) for(j=0;jn;j+) matrixAij=0;最后的动态数组赋初值也是很重要的，否则在数据处理过程中会出现错误。3） 矩阵求逆的子程序中，需要先对矩阵进行全对角占优的处理，能很好避免数据处理过程中一些“病态”的处理，而且能避免复杂的计算量问题。实验二 M序列的生成与输出计算1、 实验要求编写并调试动态模型仿真程序模型：已知白噪声v（k）数据文件为DV，数据长度L=500要求：（1）产生长度为L的M序列数据文件（2） 产生长度为L的模型输出数据文件DY2、 实验过程1） 产生M序列 M序列是一种幅值为二值的伪随机序列，周期=（2n-1）

21、（其中n为系统的阶次，即辨识的参数的个数）。 有两种产生方法： 1.采用随机函数。即先随机一个周期的M序列的值（-a或a），再据此更新其他的数据即可。#include#include#include#includemain()int i,k,n=3,m=2,p;int *a,*mfuc;FILE *fp;p=pow(m,n)-1;a=(int *)malloc(sizeof(int *)*p);mfuc=(int *)malloc(sizeof(int *)*500);srand( (unsigned)time( NULL ) );for(i=0;ip;i+) ai=0;for(i=0;ip;

22、i+) ai=0.5+(float)rand()/(RAND_MAX+1); ai=2*ai-1; printf(%dn,ai);printf(nnnnn);for(i=0;i500;i+) k=i%p; mfuci=ak; printf(%dn,mfuci);if(fp=fopen(m.dat,w)=NULL)printf(Cannot open this file!n);exit(0);for(i=0;i500;i+)fprintf(fp,%d,mfuci);fputc(n,fp);fclose(fp);fflush(stdin);getchar();2. M序列用模二加和门电路生成。我们

23、以4阶的进行说明，首先给4个初始状态，如0001，再下一位的数据为该数据的前面第3、4位数据的异或值，一次循环下去，就可以得到一个周期为15的M序列，因为均值需要为0，我们这里命0为-1，实验程序如下：#includemain()int i=0;int a=1,b=1,c=-1,d=-1,e=0;int m=500;int *m1;m1=(int *)malloc(sizeof(int *)*m);FILE *fpin;printf(%dn,a);printf(%dn,b);printf(%dn,c);printf(%dn,d);m10=a;m11=b;m12=c;m13=d; for(i=4

24、;im;i+) e=a; a=b; b=c; c=d; if (e=a) d=-1; else d=1; printf(%dn,d); m1i=d; if(fpin=fopen(u.dat,w)=NULL) printf(Cannot open this file!n); exit(0); for(i=0;im;i+) fprintf(fpin,%d,m1i); fputc(n,fpin); fclose(fpin); free(m1); fflush(stdin); getchar(); 在上面程序运行过后，会生成一个u.dat的文件，即为以后的输入的数据。2） y值的输出。 输入的值选用前

25、面我们生成的M序列，这样可以开始计算对应的输出值，这里还有白噪声的干扰，具体的编程如下：#includemain()int i;int num=500;FILE *fpu,*fpy,*fpdv;double u500,dv500,y500;double u_k1=0,u_k2=0,y_k=0,y_k1=0,y_k2=0; for(i=0;inum;i+) ui=0;dvi=0; if(fpu=fopen(u.dat,r)=NULL) printf(Cannot open u.dat file!n); exit(0); i=0; do fscanf(fpu,%lf,&ui); i+; while

26、(!feof(fpu)&i500); fclose(fpu); if(fpdv=fopen(DV,r)=NULL) printf(Cannot open DV file!n); exit(0); i=0; do fscanf(fpdv,%lf,&dvi); i+; while(!feof(fpdv)&i500); fclose(fpdv); y0=0;y1= 0; for(i=2;inum;i+) u_k1=ui-1; u_k2=ui-2; y_k=1.5*y_k1-0.7*y_k2+u_k1+0.5*u_k2+dvi-2; yi=y_k; y_k2=y_k1; y_k1=yi; for(i=

27、0;inum;i+) printf(%lfn,yi); if(fpy=fopen(y.dat,w)=NULL) printf(Cannot open y.dat file!n);exit(0); for(i=0;inum;i+) fprintf(fpy,%lf,yi); fputc(n,fpy); fclose(fpy); fflush(stdin);getchar(); 通过上面的程序，我们可以得到一个y.dat的数据输出数据，方便我们以后调用。三、实验分析1） 采用M序列来替代白噪声，是因为它没有白噪声的缺点，如：要求无限长时间等，但它是接近白噪声的，是比较理想的辨识输入的激励信号。2）

28、数据的输出值，我们是根据离散函数进行计算的，是一种理论的数据输出，我们在此加上白噪声，以模仿实际的数据输出。实验三 成批最小二乘法（BLS）1、 实验要求 编写并调试动态离散时间模型LS成批算法程序。 要求： 1.原始数据由DU和DY读出； 2.调用求逆及相乘子程序； 3.显示参数辨识结果。二、实验过程1） 成批最小二乘法的计算流程系统模型：或写成：其中，u（k）我们选为M序列产生的数据，即u.dat里的数据，v（k）为白色随机噪声。即可写成以下形式：式中，对于k=n+1,n+2,.,n+m，以上方程可以构成一个线性方程组。可以写成：其中，参数估计其中u(k)应取：1）随机序列（白噪声）；2）

29、伪随即序列（M序列）。2） 具体的程序如下：#includemain() void matrix_mutipl(float *matrixA,float *matrixB,float *matrixC,int m,int n,int q); float inv(float *p,int n); void swap(double *a,double *b); int i,j; float u500; float *y,*x,*xt,*c,*d,*o,*yt;y=(float *)malloc(sizeof(float *)*500);y0=(float *)malloc(sizeof(float

30、)*500*1);for(i=1;i500;i+) yi=yi-1+1; x=(float *)malloc(sizeof(float *)*498);x0=(float *)malloc(sizeof(float)*498*4);for(i=1;i498;i+) xi=xi-1+4; xt=(float *)malloc(sizeof(float *)*4);xt0=(float *)malloc(sizeof(float)*4*498);for(i=1;i4;i+) xti=xti-1+498; c=(float *)malloc(sizeof(float *)*4);c0=(float

31、*)malloc(sizeof(float)*4*4);for(i=1;i4;i+) ci=ci-1+4; d=(float *)malloc(sizeof(float *)*4);d0=(float *)malloc(sizeof(float)*4*498);for(i=1;i4;i+) di=di-1+498; o=(float *)malloc(sizeof(float *)*4);o0=(float *)malloc(sizeof(float)*4*1);for(i=1;i4;i+) oi=oi-1+1; yt=(float *)malloc(sizeof(float *)*498);

32、yt0=(float *)malloc(sizeof(float)*498*1);for(i=1;i498;i+) yti=yti-1+1; for(i=0;i500;i+) ui=0;for(i=0;i500;i+)yi1=0; for(i=0;i498;i+)for(j=0;j4;j+)xij=0;for(i=0;i4;i+)for(j=0;j498;j+)xtij=0;for(i=0;i4;i+)for(j=0;j4;j+)cij=0; for(i=0;i4;i+)for(j=0;j498;j+)dij=0; for(i=0;i4;i+)for(j=0;j1;j+)oij=0; for(

33、i=0;i498;i+)for(j=0;j1;j+)ytij=0; FILE *fpdu,*fpdy;if(fpdu=fopen(u.dat,r)=NULL) printf(Cant open file u.datn); exit(0); i=0; do fscanf(fpdu,%f,&ui); i+; while(!feof(fpdu)&i500); fclose(fpdu); if(fpdy=fopen(y.dat,r)=NULL) printf(Cant open file y.datn); exit(0); i=0; do fscanf(fpdy,%g,&yi0); i+; while

34、(!feof(fpdy)&i500); fclose(fpdy);for(i=0;i498;i+)for(j=0;j1;j+)ytij=yi+2j; for(j=0;j2;j+) for(i=0;i498;i+)xij=(-1)*yi-j+10;for(j=2;j4;j+) for(i=0;i498;i+)xij=ui-j+3;for(i=0;i498;i+)for(j=0;j4;j+) xtji=xij; matrix_mutipl(float *)xt,(float *)x,(float *)c,4,498,4);i=inv(*c,4);if(i!=0) for(i=0;i4;i+) pu

35、tchar(n); for(j=0;j4;j+) printf(%f ,cij); matrix_mutipl(float *)c,(float *)xt,(float *)d,4,4,498);matrix_mutipl(float *)d,(float *)yt,(float *)o,4,498,1);printf(nn);for(i=0;i4;i+)printf(%gn,oi0);getchar();void matrix_mutipl(float *matrixA,float *matrixB,float *matrixC,int m,int n,int q) int k=0,j=0,

36、i=0; for(i=0;im;i+) for(j=0;jq;j+) for(k=0;kn;k+) *(*(matrixC+i)+j)=*(*(matrixC+i)+j)+*(*(matrixA+i)+k)*(*(*(matrixB+k)+j);float inv(float *p,int n) void swap(double *a,double *b); int *is,*js,i,j,k,l; for(i=0;in;i+) putchar(n); for(j=0;jn;j+) printf(%f ,*(p+i*n+j); puts(nnnn); double temp,fmax; is=

37、(int *)malloc(n*sizeof(int); js=(int *)malloc(n*sizeof(int); for(k=0;kn;k+) fmax=0.0; for(i=k;in;i+) for(j=k;jfmax) fmax=temp; isk=i;jsk=j; if(fmax+1.0)=1.0) free(is);free(js); printf(no inv); return(0); if(i=isk)!=k) for(j=0;jn;j+) swap(p+k*n+j,p+i*n+j);/交换指针 if(j=jsk)!=k) for(i=0;in;i+) swap(p+i*n

38、+k,p+i*n+j); /交换指针 pk*n+k=1.0/pk*n+k; for(j=0;jn;j+) if(j!=k) pk*n+j*=pk*n+k; for(i=0;in;i+) if(i!=k) for(j=0;jn;j+) if(j!=k) pi*n+j=pi*n+j-pi*n+k*pk*n+j; for(i=0;i=0;k-) if(j=jsk)!=k) for(i=0;in;i+) swap(p+j*n+i),(p+k*n+i); if(i=isk)!=k) for(j=0;jn;j+) swap(p+j*n+i),(p+j*n+k); free(is); free(js); return 1;void swap(double *a,double *b) double c; c=*a; *a=*b; *b=c;运行的实验结果如下：三、实验分析1） 在运算过程会出现求逆的过程，因此会不可避免地增加整个运算量，尤其是当整个系统的阶次较高时候，会占用较大内存。2）由于运算的结果相乘等处理过程中，会不可避免地出现偏差，如上，我们得到的最大的数据误差为4.8%，已经能较好地辨识系统参数。3）成批最小二乘法不适合系统的在线辨识。实验四 递推最小二乘法（RLS）1、 实验要求编写并调试动态离散时间模型LS递推算