系统辨识实验报告.docx

系统辨识实验报告1

系统辨识实验报告

学院：信息科学及技术学院
专业：自动化
日期：2016/4/26

目录

实验1 4
一．实验内容及要求： 4
二．实验原理： 4
三．软件设计思想： 4
四．程序结构框图： 5
五．运行示意图： 5
实验2 8
一．实验内容及要求： 8
二．实验原理： 8
三．软件设计思想： 9
四．程序设计框图： 9
五．程序运行流程图： 10
实验3 12
一．实验内容及要求： 12
二．实验原理： 12
三．程序数据流程图： 12
四．实验运行结果： 13
实验4 14
一．实验内容及要求： 14
二．实验原理： 14
三．数据递推关系图： 14
四．实验运行结果： 15
心得体会 16
附录（实验代码） 17
1. LabWork1 17
2. LabWork2 21
3. LabWork3 23
4. LabWork4 26

实验1
一．实验内容及要求：
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）的这样的一个矩阵就可以求得矩阵的逆。若矩阵FangZ不是满秩矩阵时，FangZ没有FZNi 。通过这样的求逆方式，避免了大方阵的求取行列式运算。
三．软件设计思想：
1. 确定该软件的功能主要有：键盘输入两个矩阵然后相乘；文本data输入两个矩阵将结果放在文本result中；键盘输入一个方阵求得其逆矩阵。其中前两个的矩阵相乘运算部分设置为一个函数Mul。
2. 在main函数中提供两个关于矩阵的选择：multiplitation；invertion。其相对应的子函数为MulOp(a,b,c)，Inv()。
3. 在MulOp(a,b,c)子函数中，有两种输入矩阵的方式：way1，way2。相对应的功能为键盘输入，文本输入。并且两者在处理矩阵时，都调用了Mul函数。
4. 在Inv()子函数中，输入和显示原矩阵，和其相应的逆矩阵。调用qiuni（double FangZ[][M],double FZNi[][M],int n）子函数，可以都得到原矩阵的逆矩阵。但当原矩阵不可逆时，系统输出为”The array is not invertible!“。
四．程序结构框图：

way1键盘输入矩阵

Mul矩阵乘积子函数

MulOp子函数
矩阵乘积

way2文本输入
矩阵

main函数

Inv子函数
矩阵求逆

qiuni子函数

退出程序

五．运行示意图：
1. main函数的主界面：

2. MulOp子函数的界面：

3. Inv子函数的界面：

4. 通过键盘操作计算两个矩阵的乘积：

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+…+cnx^n
可知9阶移位寄存器的多项式为fx=x9+x4+1，及可得c=[0,0,0,1,0,0,0,0,1]
9级线性移位寄存器：
+
+
+

C8
C1
C0

输出
a0
a1
a8
a7

图中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,1]

M序列的长度为L=500

i=1

Y
输出M序列
i<L

N

M序列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函数中给出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子函数
产生DU
main函数

input子函数
输入DU，DV

output子函数
产生DY

deal子函数
得DY

退出程序

五．程序运行流程图：
1. main函数：

2. gener函数：

3. output函数：

实验3
一．实验内容及要求：
编写并调试动态离散时间模型LS成批算法程序。
要求：（1）原始数据由DU和DY读出。
（2）调用求逆及相乘子程序。
（3）显示参数辨识结果。
二．实验原理：
1.批次处理的方法就是把所有的数据采集到一次性进行处理，但前提是白色噪声、及M序列所共同作用而产生的输出，才能使用最小二乘法。虽然这种方法的计算量庞大，但经常用于处理时不变系统，方法简单。
2.构造模型Y=X*sita+V
3．Y=[DY(3) DY(4) DY(499) ]已知n=2,L=500,可知m=497
所以有X=-y(2)-y(1)-y(3)-y(2)-y(498)-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-1)+b2u(k-2)+v(k)。
2. 引入信息矩阵P(k)，维数为4*4，初始化为10^5~10^12*E的一单位阵。gama(k)为修正系数，为无穷小标量。x(k+1)=[-y(k) -y(k-1) u(k) u(k-1)]。sita'=[a1 a2 b1 b2]。y(k+1)=DY[k+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)];
三．数据递推关系图：

四．实验运行结果：
1.main函数界面示意图：

6. 只显示结果界面：

7. 显示过程的参考界面：

心得体会
我很喜欢这个课程的期末考核方式，不用再拘泥于在题目当中对该课程的了解，而是通过4个C语言设计的练习来达到学习的目的，而且对以后的学习还有很大的帮助。
在编写C语言的过程中，也遇到了一些阻碍，特别是在编写第3,4个的时候。比如：用的数组太多，并且未将其初始化，运行出来的结果经常是很长的一段随机数；或者一模一样的程序有时候就可以正常运行，有时候就总是出现报错……这些都是让我心塞了两周的问题。在这些问题解决之后，运行出来的结果却及实际模型参数的出入有点大，于是又重新查找第2个实验是否是M序列产生的方式有问题。通过对初始化寄存器赋给不同的值，可以让结果及真实模型参数之间的误差达到最小。
通过这学期对系统辨识这个课程的学习，让我了解到系统辨识在建立数学模型的方面的重要性。对于不了解系统的工作机理的人来说，也有了一个可以知道系统模型的构建方法。但是也必须要我们对所构建模型有一个较为清醒的认识，比如要知道模型的类型，如果是动态模型，则要知道模型的参数或者阶次（但是权函数模型就不需要事前知道模型的参数或者阶次，不知道这种方法有没有什么我们不知道的缺陷）。在处理数据上，又一次让我了解到世上没有两全其美的事情，计算的复杂程度和精度就好像鱼和熊掌不能兼得。也正是因为这样，才会成为促使人们在这方面的不断寻找最优化算法的动力。在对模型进行研究时，都是从最简单的模型开始研究，比如在没有噪声的情况下，所得到的模型参数就为系统模型真实参数。进一步，在白色噪声的情况下，所得到的参数就为系统模型真实参数的估计……通过这样的推理，就可以得到LS在系统辨识问题上的普适性，不管系统受怎样的干扰，只要通过一系列的变换，最后还是通过LS来解决模型参数的确定工作。
虽然我们对化学工艺不怎么了解，但通过学习系统辨识，让我们对已经投运后的工业生产设备的性能有一个跟踪性的了解，可以对系统做出一个较为完整的评估，从而来改善生产工艺或者是及时更换生产设备，从而来达到最大的经济效益。
这学期主要就是介绍ＬＳ方法，虽然ＬＳ可以处理一些特殊情况的有噪声和没有噪声的系统，但对于那些噪声模型都不可预测的系统中，ＬＳ还是显得无能为力。希望以后还有机会学习其他的系统辨识方法。

附录（实验代码）
1. LabWork1
#define M 100
double FangZ[M][M]={0},FZNi[M][M]={0};

/* 求矩阵的乘法函数 */
void Mul(double a[][M],double b[][M],double c[][M],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++)
c[i][k]+=a[i][j]*b[j][k];
}
/* 求方阵的逆 */
void qiuni(double FangZ[][M],double FZNi[][M],int n)
{
double E[M][M],value=1.0,stemp,temp=0.0,san,c[M][M];
int i,j,row,nextrow,flag=0;
int col,switchtime=0;
/* 构造原矩阵的增广矩阵 */
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(j==i) E[i][j]=1;
else E[i][j]=0;
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
c[i][j]=FangZ[i][j];
for(j=0;j<n;j++)
c[i][j+n]=E[i][j];
}
/* 对增广矩阵进行运算（FangZ|E）->（E|FZNi） */
for(row=0;row<n-1;row++)
{
nextrow=row+1;
if(c[nextrow][row]==0)
{
while(c[nextrow][row]==0)
{
nextrow++;
if(nextrow==n)
{
flag=1;
break;
}
}
if(flag==1) continue;
switchtime++;
/* 交换第一个数字为0的行 */
for(col=0;col<2*n;col++)
{
stemp=c[row][col];
c[row][col]=c[nextrow][col];
c[nextrow][col]=stemp;
}
}
for(nextrow=row+1;nextrow<n;nextrow++)
{
temp=c[nextrow][row]*1.0/c[row][row];
for(col=0;col<2*n;col++)
c[nextrow][col]+=-temp*c[row][col];
}
}
/* 得到前一个矩阵的对角线上元素的乘积 */
for(row=0;row<n;row++)
value*=c[row][row];
/* 判断乘积！若为零，则原矩阵没有逆矩阵；反之亦然 */
if(value==0) printf("\n\nThe array is not invertible!\n\n");
else
{
for(row=n-1;row>=1;row--)
for(i=row-1;i>=0;i--)
{ san=c[i][row]*1.0/c[row][row];
for(col=row;col<2*n;col++)
c[i][col]-=c[row][col]*san;
}
/* 得到标准的（E|A） */
for(row=0;row<n;row++)
for(col=0;col<n;col++)
c[row][col+n]/=c[row][row];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
FZNi[i][j]=c[i][j+n];
}
}
2. LabWork2
/* 处理v(k),u(k)的数据得输出y(k) */
void deal(double DV[M],double DU[M],double DY[M])
{
int k;
DY[0]=0;
DY[1]=0;
for(k=2;k<500;k++)
DY[k]=1.5*DY[k-1]-0.7*DY[k-2]+DU[k-1]+0.5*DU[k-2]+DV[k];
}
/* 产生m序列u(k)输出到文件 */
void gener()
{
int a[9]={1,1,0,0,0,1,0,1,0},c[9]={0,0,0,1,0,0,0,0,1},aa[9];
int DU[M],temp,sum,i,j;
DU[0]=a[8];
FILE *fp;
system("cls");
if((fp=fopen("DU.txt","w"))==NULL)
{
printf("cannot open the file!\n");
exit(0);
}
for(i=1;i<M;i++)
{
sum=0;
for(j=0;j<9;j++)
sum+=a[j]*c[j];
aa[0]=sum%2;
for(j=1;j<9;j++)
aa[j]=a[j-1];
for(j=0;j<9;j++)
a[j]=aa[j];
DU[i]=a[8];
}
for(i=0;i<M;i++)
fprintf(fp,"%d ",DU[i]);
fclose(fp);
printf("\n\n\n\t\tPlease refer to the )!\n\n");
system("pause");
}
3. LabWork3
#define M 500
#define N 4
void MUL1(double TX[][M],double X[][N],double ji1[][N])
{
int i,j,k;
for(i=0;i<N;i++)
for(k=0;k<N;k++)
for(j=0;j<M-3;j++)
ji1[i][k]+=TX[i][j]*X[j][k];
}

void MUL2(double ni[][N],double TX[][M],double ji2[][M])
{
int i,j,k;
for(i=0;i<N;i++)
for(k=0;k<M-3;k++)
for(j=0;j<N;j++)
ji2[i][k]+=ni[i][j]*TX[j][k];
}

void MUL3(double ji2[][M],double Y[M],double xishu[N])
{
int i,j;
for(i=0;i<N;i++)
for(j=0;j<M-3;j++)
xishu[i]+=ji2[i][j]*Y[j];
}
void main()
{
system("cls");system("color 75");
int i,j;
double DU[M],DY[M];
double X[M][N]={0},temp,Y[M],TX[N][M],ji1[N][N]={0},ni[N][N]={0},ji2[N][M]={0},xishu[N]={0};
read(DU,DY);

/* 得到大矩阵X( 498 * 4 ) */
for(i=0;i<M-3;i++)
{
for(j=0;j<2;j++)
X[i][j]=-DY[i+j+1];
for(j=0;j<2;j++)
X[i][j+2]=DU[i+j+1];
temp=0;temp=X[i][0];X[i][0]=X[i][1];X[i][1]=temp;
temp=0;temp=X[i][2];X[i][2]=X[i][3];X[i][3]=temp;
}
/* 得到Y（498 * 1）=[DY(1) DY(498)]的列矩阵 */
for(i=0;i<M-3;i++)
Y[i]=DY[i+3];
/* 得到X（498 * 4）的转置矩阵 TX（4 * 498） */
for(i=0;i<N;i++)
for(j=0;j<M-3;j++)
TX[i][j]=X[j][i];
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:\n\n");
for(i=0;i<N;i++)
printf("\n\n\t %.4f \n",xishu[i]);
printf("\n\n");
}
4. LabWork4
#define M 500
#define N 4
double sita[N][N];
void MUL(double a[][N],double b[][N],double c[][N],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++)
c[i][k]+=a[i][j]*b[j][k];
}
void trans(double d[][N],double e[][N],int n5,int n6)
{
int i,j;
for(i=0;i<n6;i++)
for(j=0;j<n5;j++)
e[i][j]=d[j][i];
}
void consult()
{
system("cls");
int i,j,k;
double DU[M],DY[M];
double tx[N][N]={0},p[N][N],y,x[N][N],temp,gama;
double ji1[N][N]={0},ji2[N][N]={0};
double ji3[N][N]={0},ji4[N][N]={0};
double ji5[N][N]={0},ji6[N][N]={0},ji7[N][N]={0};
read(DU,DY);
/* 对P（ 4 * 4 ）和sita进行初始化 */
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(i==j) p[i][j]=10e9;
else p[i][j]=0;
}
for(i=0;i<N;i++)
sita[i][0]=0.0;
/* 大循环开始 */
for(i=2;i<M;i++)
{
/* 求y(i+1) */
y=DY[i];
/* 求x(i+1)的(4*1)的列阵 */
for(k=0;k<1;k++)
{
for(j=0;j<2;j++)
x[j][k]=-DY[i+j-2];
for(j=0;j<2;j++)
x[j+2][k]=DU[i+j-2];
temp=x[0][k];x[0][k]=x[1][k];x[1][k]=temp;
temp=x[2][k];x[2][k]=x[3][k];x[3][k]=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(ji1,x,ji2,1,4,1); /* ji1(1 4)*x(4 1)=ji2(1 1) */
gama=1.0/(1+ji2[0][0]);
/* 求sita(i+1)的（4*1）的列阵sita(i+1)=sita(i)+gama(i+1)*p(i)*x(i+1)*(y(i+1)-tx(i+1)*sita(i)) */
MUL(tx,sita,ji3,1,4,1); /* tx(1 4)*sita(4 1)=ji3(1 1) */
MUL(p,x,ji4,4,4,1); /* p(4 4)*x(4 1)=ji4(4 1) */
for(j=0;j<N;j++)
ji4[j][0]=ji4[j][0]*gama*(y-ji3[0][0]);
for(j=0;j<N;j++)
sita[j][0]=sita[j][0]+ji4[j][0];
/*求p(i+1)的（4*4）的矩阵 p(i+1)=p(i)-gama(i+1)*p(i)*x(i+1)*tx(i+1)*p(i) */
MUL(p,x,ji5,4,4,1); /* p(4 4)*x(4 1)=ji5(4 1) */
MUL(ji5,tx,ji6,4,1,4); /* ji5(4 1)*tx(1 4)=ji6(4 4) */
MUL(ji6,p,ji7,4,4,4); /* ji6(4 4)*p(4 4)=ji7(4 4) */
for(j=0;j<N;j++)
for(k=0;k<N;k++)
p[j][k]=p[j][k]-gama*ji7[j][k];
/* 使得计算过程当中使用的数组都清零 */
for(j=0;j<N;j++)
for(k=0;k<N;k++)
{
ji1[j][k]=0;ji2[j][k]=0;ji3[j][k]=0;ji4[j][k]=0;
ji5[j][k]=0;ji6[j][k]=0;ji7[j][k]=0;x[j][k]=0;tx[j][k]=0;
}
}
/* 结果显示 */
printf("\tthe xishu are:");
for(i=0;i<N;i++)
printf("\n\n\n\t\t %.4f \n",sita[i][0]);
printf("\n\n");

system("pause");
}
32 / 33