jieci.cpp

系统辨识中

/*******************************************************************/
/*-----------递推最小二乘法估计及模型阶次检验实验------------------*/
/*******************************************************************/

#include <dos.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>

const int P=4,Np=(int)pow(2,P)-1,Len=600,realN=2; //P为M序列的阶次，Np循环周期,realN理想阶次
static int nL=2,L=150+nL*150,Nbeg=1,Nend=4,resultN;
static double *u=new double[Len];	//M序列
static double *z=new double[Len];	//有噪声的输出量
static double *v=new double[Len];	//白噪声
static double *e=new double[Len];	//有色噪声
static double *y=new double[Len];	//无噪声的输出量
static double *UU=new double[500*2*realN];

static double Sigma=0.5,a=1,miu=0.998;  //Sigma噪声标准差，a为M序列的幅度，miu遗忘因子
static double ab[4]={-1.5,0.7,1,0.5};   //真实参数
static double J[10],Lamda[10],Kj[10];   //损失函数、噪声标准差、模型静态增益估计
static double ZXB,deta1,deta2,detak;    // 噪信比、参数估计平方相对偏差、参数估计平方根偏差、静态增益估计相对偏差   
static double TN[10],Q[4];              //定阶法估计、二阶辨识结果、

const double ta[3]={3.09,3.03,3.01};    //风险水平0.05时L＝100，300，500的阀值

/*********************以下是与矩阵相关的一些运算*************************/
/*----------------------构造一个单位阵----------------------------------*/
void CreateIMatrix(int n,double *im)
{
	for(int i=0;i<n*n;i++)
	{
		if(int(i/n)==i%n)
			 *(im+i)=1;
		else *(im+i)=0;
	}
}
/*-----------------------矩阵加法-----------------------------------------*/
void AddMatrix(double *jia1,double *jia2,double *he,int row,int vol)
{
	int n=row*vol;
	for(int i=0;i<n;i++)
	{
		*(he+i)=*(jia1+i)+(*(jia2+i));
	}
}
/*-----------------------矩阵数乘-----------------------------------------*/
void MulNumMatrix(double *one,double num,double *ji,int row,int vol)
{
	int n=row*vol;
	for(int i=0;i<n;i++)
	{
		*(ji+i)=num*(*(one+i));
	}
}
/*------------------------矩阵乘法----------------------------------------*/
void MulMatrix(double *one,int row,int rv,int vol,double *two,double *ji)
{
	int i,j,k;
	double temp;
	for(i=0;i<row;i++)
	{
		for(j=0;j<vol;j++)
		{
			temp=0;
			for(k=0;k<rv;k++)
			{
				temp+=(*(one+rv*i+k))*(*(two+k*vol+j));
			}
			*(ji+i*vol+j)=temp;
		}
	}
}
/*********************以下是模型输出************************/
/*-------------------利用位运算生成M序列-------------------*/
void CreateU(int L)
{
	int i;
	unsigned m=2;		//给初始值u[0]=1
	for(i=0;i<L;i++)
	{
		m=m>>1;			//右移位
		if(m&1)
			u[i]=-a;	//逻辑为1，幅度－a
		else 
			u[i]=a;     //逻辑为0，幅度a
		if((m&9)==1||(m&9)==8)
			m=m|16;     //循环
	}
}
/*--------------------生成白噪声v-----------------------------*/
void CreateV(int L)
{
	double M=32768,A=179,x0=11,ksai,x11;//乘同余法参数
	int i,k;
	x11=x0;
	for(k=0;k<L;k++)
	{
		ksai=0;
		for(i=1;i<12;i++)
		{
			x11=A*x11;
			x11=fmod(x11,M);        //生成u[0,1]均匀分布随机数
			ksai=ksai+float(x11/M);
		}
		v[k]=Sigma*(ksai-6.0);      //生成正态分布的白噪声
	}
}
/*--------------------生成有色噪声e-----------------------------*/
void CreateE(int L)
{
	int k;
	for(k=0;k<L;k++)
	{
		if(k==0)e[k]=v[k];                //零阶系统
		else 
			if(k==1)
				e[k]=v[k]-ab[0]*e[k-1];	  //一阶系统
		    else 
				e[k]=v[k]-ab[0]*e[k-1]-ab[1]*e[k-2];//二阶系统
	}
}
/*----------------------传递函数的输出Z(k)-------------------*/
void CountZ(int L,double *ab)
{
	int k;
	y[0]=0;
	y[1]=u[0];
	for(k=2;k<L;k++)
		y[k]=-ab[0]*y[k-1]-ab[1]*y[k-2]+ab[2]*u[k-1]+ab[3]*u[k-2];//从方框图中得出的模型表示
	for(k=0;k<L;k++)
		z[k]=y[k]+e[k];      //含有有色噪声的输出量
}
/**************************以下是RLS同F－test子程序***********************/
/*-----------------------递推最小二乘法----------------------------------*/
void Fzhen(int n)
{
	double *thita=new double[2*n];//辨识出的参数
	double *hk=new double[2*n];   //已知数据
	double *P1=new double[4*n*n],*P2=new double[4*n*n];//P1、P2分别为前一时刻和当前时刻的P矩阵
	double *K=new double[2*n];    //增益矩阵K
	double te,tm,*ti=new double[4*n*n],*temp;
	int i,k;
	J[n]=0;                       //损失函数
	/*------------初始化thita和P-------------------*/
	CreateIMatrix(2*n,P1);
	MulNumMatrix(P1,1E+6,P1,2*n,2*n);   //为收敛快，P＝10(6)*I
  	for(i=0;i<2*n;i++)
	{
		thita[i]=0.001;                 //thita的初始值很小
	    UU[i]=thita[i];
	}
	/*----------------递推核心算法------------------*/
	for(k=n;k<L;k++)
	{
		for(i=0;i<n;i++)
		{
			hk[i]=-z[k-i-1];
			hk[i+n]=u[k-i-1];           //hk为已知的输入输出量
		}
		//计算K(k)
		MulMatrix(hk,1,2*n,2*n,P1,ti);
		MulMatrix(ti,1,2*n,1,hk,K);
		tm=1/(K[0]+miu);
		MulMatrix(P1,2*n,2*n,1,hk,ti);
		MulNumMatrix(ti,tm,K,2*n,1);
		//计算P(k)
		MulMatrix(K,2*n,1,2*n,hk,ti);
		MulNumMatrix(ti,-1,ti,2*n,2*n);
		for(i=0;i<2*n;i++)
			ti[i*2*n+i]+=1;
		MulNumMatrix(ti,1/miu,ti,2*n,2*n);
		MulMatrix(ti,2*n,2*n,2*n,P1,P2);
		//计算thita(k)
		MulMatrix(hk,1,2*n,1,thita,ti);
		te=z[k]-ti[0];
		MulNumMatrix(K,te,ti,2*n,1);
		AddMatrix(thita,ti,thita,2*n,1);
       	//计算J(k)
		J[n]=(J[n]+tm*te*te);
		temp=P1;
		P1=P2;
		P2=temp;
	}
	if(n==realN)
	{	      
		for(i=0;i<2*n;i++)
			Q[i]=thita[i];		//存储理想阶次下的辨识结果
	}
	Lamda[n]=sqrt(J[n]/(L-2*n));//估计噪声标准差
	/*--------------估计模型静态增益----------------*/
	te=0;tm=0;
	for(i=0;i<n;i++)
	{
		te+=thita[i];
		tm+=thita[i+n];
	}
	if(te==-1)
		Kj[n]=1E+10;
	else 
		Kj[n]=tm/(1+te);
/*------------保存递推过程中的辨识参数数据-----------*/
    FILE *file=fopen("最终结果.txt","w");
	if((file=fopen("最终结果.txt","w"))==NULL)
	{
		printf("cannot open the file.\n");
	    return;
	}
	fprintf(file,"%f\t\t%f\t\t%f\t\t%f\t\t\n",Q[0],Q[1],Q[2],Q[3]);
    fclose(file);
}
/*--------------------F-Test定阶--------------------------------*/
bool AcceptN(int n)
{
	if(n<1||n>8)return false;
	TN[n]=(J[n]/J[n+1]-1)*(L-2*n-2)/2;//统计量
	if(TN[n]<=ta[nL])return true;     //接受H：n+1>=n0；阶次估计值为n+1
	else return false;                //拒绝H：n+1>=n0
}
/*----------------------过程输出方差的计算------------------------------*/
void fangcha(double *a,double *b,double *at,double *bt,int n)
{
	int i,k=n-1;
	for(i=0;i<n;i++)               //具有周期性
	{
		at[k*n+i]=a[i];
		bt[k*n+i]=b[i];
	}
	for(k=n-2;k>=0;k--)
	{
		for(i=0;i<=k;i++)
		{
			at[k*n+i]=at[(k+1)*n+i]-(at[(k+1)*n+k+1]*at[(k+1)*n+k+1-i])/at[(k+1)*n];
			bt[k*n+i]=bt[(k+1)*n+i]-(bt[(k+1)*n+k+1]*at[(k+1)*n+k+1-i])/at[(k+1)*n];
		}
	}
}
/*--------------------------性能计算---------------------------*/
void CountG()
{
	/*------------计算噪声方差----------------*/
	double Deltae=0,e0=0;
	int i;
	for(i=0;i<L;i++)
		e0+=e[i];
	    e0/=L;                       //有色噪声均值
	for(i=0;i<L;i++)
		Deltae+=pow((e[i]-e0),2);   //有色噪声方差
	Deltae/=(L-1);
	/*-----------计算过程输出方差--------------*/
	double Deltay=0,a[3],b[3];
	a[2]=ab[1];
	a[1]=ab[0];
	a[0]=1;
	b[0]=0;
	b[1]=ab[2];
	b[2]=ab[3];

	int n=3;
	double *at=new double[n*n],*bt=new double[n*n];
	fangcha(a,b,at,bt,n);
	for(i=0;i<n;i++) 
		Deltay+=pow(bt[i*n+i],2)/at[i*n+i];
	    Deltay=Deltay/a[0];
	/*-----------计算噪信比---------------*/
    ZXB=sqrt(Deltae/Deltay);
	printf("系统输出噪信比:%f\n",ZXB);
	/*----------参数估计平方相对偏差------*/
	deta1=0;
	for(i=0;i<2*realN;i++)
		deta1+=pow((ab[i]-Q[i])/ab[i],2);
	    deta1=sqrt(deta1);
    printf("参数估计平方相对偏差:%f\n",deta1);
	/*---------参数估计平方根偏差-----------*/
	double temp=0;
	deta2=0;
	for(i=0;i<2*realN;i++)
	{
		deta2+=pow((ab[i]-Q[i]),2);
		temp+=pow(ab[i],2);
	}
	deta2=sqrt(deta2/temp);
    printf("参数估计平方根偏差:%f\n",deta2);
	/*---------静态增益估计相对偏差-----------*/
	temp=0,Kj[0]=0;
	for(i=0;i<realN;i++)
	{
		temp+=ab[i];        //ai之和
		Kj[0]+=ab[i+realN]; //bi之和
	}
	if(temp==-1)
		Kj[0]=1E+10;
	else 
		Kj[0]=Kj[0]/(1+temp);
	detak=sqrt(fabs((Kj[0]-Kj[realN])/Kj[0]));
    printf("静态增益估计相对偏差:%f\n",detak);
	return;
}
/*-------------------------主函数-----------------------------*/
void main()
{
loop:
    printf("----------递推最小二乘法及模型阶次检验实验-------------\n");
	printf("请输入标准差值(范围0～1):\n");
	scanf("%lf",&Sigma);
	getchar();
	if(Sigma<0||Sigma>1)  //标准差在[0,1]之间
	{
		printf("error.\n");
		goto loop;        //重新输入sigma
	}
	printf("遗忘因子：\n");
	scanf("%f",&miu);
	getchar();
	printf("数据长度：\n");
	scanf("%d",&L);
	getchar();
	printf("开始阶数：\n");
	scanf("%d",&Nbeg);
	getchar();
	printf("终止阶数：\n");
	scanf("%d",&Nend);
	getchar();

	if(Nbeg>=Nend)resultN=-2;
	CreateU(L);  //生成M序列
	CreateV(L);  //生成输出白噪声v
	CreateE(L);  //生成有色噪声e
	CountZ(L,ab);//传递函数的输出Z(k)
	int i=Nbeg;
	Fzhen(i);    
	for(i++;i<=Nend;i++)   //循环确定理想阶次
	{
		Fzhen(i);
		if(AcceptN(i-1))
			break;
	}
	printf("理想阶次：%d\n",realN);
  	if(Nbeg>realN||i<realN)
		Fzhen(realN);
	printf("辨识结果:\n");
	for(i=0;i<2*realN;i++)
	{
		printf("%f\t",Q[i]);
	}
	 printf("\n");
     printf("噪声标准差估计：%f\n",Lamda[realN]);
	 printf("静态增益估计：%f\n",Kj[realN]);
	   CountG();
 	if(i>Nend)resultN=-1;
	else resultN=i-1;
	printf("\n");
	goto loop;
  	return;
}