els_ro_aicdlg.cpp

实现一次完成增广最小二乘法,主要用于系统模型辨识

	{
		//model==TRUE,使用模型一
		z[0]=e[0];
		z[1]=1.5*z[0]+1.08*u[0]+e[1]-1.0*e[0];
		for(i=2;i<m_nwnLEN;i++)
		{
			z[i]=1.5*z[i-1]-0.7*z[i-2]+u[i-1]+0.5*u[i-2]+e[i]-1.0*e[i-1]+0.2*e[i-2];
		}
	}
	else
	{
		z[0]=e[0];
		z[1]=1.185*z[0]+1.08*u[0]+e[1]-1.0*e[0];
		z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+e[2]-1.0*e[1]+0.2*e[0];
		z[3]=1.185*z[2]-0.814*z[1]+0.518*z[0]+1.08*u[2]-0.745*u[1]+0.475*u[0]+e[3]-1.0*e[2]+0.2*e[1];
		z[4]=1.185*z[3]-0.814*z[2]+0.518*z[1]-0.349*z[0]+1.08*u[3]-0.745*u[2]+0.475*u[1]-0.253*u[0]+e[4]-1.0*e[3]+0.2*e[2];
		for(i=5;i<m_nwnLEN;i++)
		{
			z[i]=1.185*z[i-1]-0.814*z[i-2]+0.518*z[i-3]-0.349*z[i-4]+0.117*z[i-5]\
				+1.08*u[i-1]-0.745*u[i-2]+0.475*u[i-3]-0.253*u[i-4]+0.123*u[i-5]+e[i]-1.0*e[i-1]+0.2*e[i-2];
		}
	}

	WriteMatrix2File("u",u,m_nwnLEN,1);
	WriteMatrix2File("z",z,m_nwnLEN,1);
	WriteMatrix2File("e_CAR",e,m_nwnLEN,1);
	for(i=0;i<m_nwnLEN;i++)
	{
		ZZZ[i*3]=(int)(10*z[i]);
		if(!model)ZZZ[i*3]*=4;
		zzz[i*3]=z[i];
	}
	
	//1.首先辨识高阶CAR模型参数，辨识参数依阶次每次增加二个，并由此参数辨识结果计算新的白噪声
	aic(m_nwnLEN,2);
	PrintTextandGraph(2);//输出CAR模型参数到文件；

	//2.根据由CAR模型辨识逼近的估计结果计算误差，形成白噪声，重新进行CARMA模型辨识，辨识参数每次增加三个。
	if((fp2=fopen("e_CARMA.txt","r"))==NULL)return; 
    for(i=0;i<m_nwnLEN-10;i++)
    { 
		fscanf(fp2,"%lf ",&e[i]);//此时的e 为CARMA模型逼近的白噪声
	}
    fclose(fp2);
	for(i=m_nwnLEN-1;i>=10;i--)
	{
		e[i]=e[i-10];
	}
	for(i=0;i<10;i++)
	{
		e[i]=0;
		u[i]=0;
		z[i]=0;
	}
	WriteMatrix2File("e_CARMA",e,m_nwnLEN,1);
	aic(m_nwnLEN,3);
	PrintTextandGraph(3);//输出CARMA模型参数到文件，并输出到图形界面。

	free(z);free(u);free(e);	
}

void CELS_RO_AICDlg::aic(int num, int order_num)
{
	//结合AIC的模型阶次递推算法,本程序最大能算到10阶。
	//order_num=2，CAR模型辨识，辨识参数每次增加二个，并形成CARMA模型辨识用的白噪声
	//order_num=3，CARMA模型辨识，辨识参数每次增加三个，输出最终辨识结果
	int N=num-10;
    int i,j,nabc;//nabc为模型阶次。
	double sigama[1];
	double *X1T,*Y,*Ye,*e1,*X1,*X2T,*X2,*XTX_inv;
	double sita2[3],*sita_temp;
	double *X1TX1,*X1TY,X2TX2[9],*X2TX1,*X1TX2;
	double *B1,*AA,B2[9],BB[9],*AX2T,*BX2T;
	double *CC,*temp1,*temp2;
	nabc=1;	//na=nb=nc=1;
	X1T=(double*)calloc(order_num*10*N,sizeof(double));
	X2T=(double*)calloc(order_num*N,sizeof(double));
	X2=(double*)calloc(N*order_num,sizeof(double));
	Y=(double*)calloc(N,sizeof(double));
	Ye=(double*)calloc(N,sizeof(double));
	e1=(double*)calloc(N,sizeof(double));
	BX2T=(double*)calloc(order_num*N,sizeof(double));

	///////////////
	for(i=0;i<N;i++)
	{
		Y[i]=z[10+i];
	}
	//赋值给X1T矩阵
	for(j=0;j<N;j++)
	{
		X1T[0*N+j]=-z[9+j];
		X1T[1*N+j]= u[9+j];
		if(order_num==3)	X1T[2*N+j]= e[9+j];
	}

	//将X1T矩阵转置得X1矩阵。
	X1=(double*)calloc(N*order_num,sizeof(double));
	for(i=0;i<order_num;i++)
	{
		for(j=0;j<N;j++)
		{
			X1[j*order_num+i]=X1T[i*N+j];
		}
	}
	//计算X1TX1，并求逆。
	X1TX1=(double*)calloc(order_num*order_num,sizeof(double));
	brmul(X1T,X1,order_num,N,order_num,X1TX1);
	brinv(X1TX1,order_num);
	//求参数估计值sita
	X1TY=(double*)calloc(order_num,sizeof(double));
	brmul(X1T,Y,order_num,N,1,X1TY);
	brmul(X1TX1,X1TY,order_num,order_num,1,sita1);
	free(X1TY);
	//求预测误差矩阵e1。
	brmul(X1,sita1,N,order_num,1,Ye);
	for(i=0;i<N;i++)
	{
		e1[i]=Y[i]-Ye[i];
	}
	//求AIC值、并保存此时的阶次、估计参数sita值。
	brmul(e1,e1,1,N,1,sigama);
	AIC0=N*log10(sigama[0]/N)+order_num*(order_num*nabc);
	AIC[0]=AIC0;
	nabc0=1;
	memcpy(sita0,sita1,order_num*sizeof(double));
	if(order_num==2)
	{
		WriteMatrix2File("e_CARMA",e1,N,1);//生成CARMA模型用的白噪声e
	}
	for(i=10;i<N+10;i++)
	{
		ZZZ[i*3+order_num-1]=(int)(10*Ye[i-10]);
		if(!model)ZZZ[i*3+order_num-1]*=4;
		zzz[i*3+order_num-1]=Ye[i-10];
	}
	if(order_num==3)
		WriteMatrix2File("ZZZ",zzz,m_nwnLEN,3);

	//
	//na=nb=nc=n;   X=[X1，X2]
	for(int n=2;n<=10;n++)
	{	
		nabc=n;
		//模型阶次增加，重新赋值给X2T矩阵
		for(j=0;j<N;j++)
		{
			X2T[0*N+j]=-z[10-n+j];
			X2T[1*N+j]= u[10-n+j];
			if(order_num==3)	X2T[2*N+j]= e[10-n+j];
		}
		//X2T矩阵转置得X2矩阵
		for(i=0;i<order_num;i++)
		{
			for(j=0;j<N;j++)
			{
				X2[j*order_num+i]=X2T[i*N+j];
			}
		}
		X1TX2=(double*)calloc(order_num*(n-1)*order_num,sizeof(double));
		X2TX1=(double*)calloc(order_num*order_num*(n-1),sizeof(double));
		B1=(double*)calloc(order_num*(n-1)*order_num,sizeof(double));
		AA=(double*)calloc(order_num*(n-1)*order_num,sizeof(double));
		AX2T=(double*)calloc(order_num*(n-1)*N,sizeof(double));
		sita_temp=(double*)calloc(order_num*(n-1),sizeof(double));
		//阶次增加，重新计算增加个数后的所有参数估计值并保存。
		brmul(X2T,X1,order_num,N,order_num*(n-1),X2TX1);
		brmul(X1T,X2,order_num*(n-1),N,order_num,X1TX2);
		brmul(X2T,X2,order_num,N,order_num,X2TX2);
		brmul(X1TX1,X1TX2,order_num*(n-1),order_num*(n-1),order_num,B1);
		brmul(X2TX1,B1,order_num,order_num*(n-1),order_num,B2);
		for(i=0;i<order_num*order_num;i++)
		{
			BB[i]=X2TX2[i]-B2[i];
		}
		brinv(BB,order_num);
		brmul(B1,BB,order_num*(n-1),order_num,order_num,AA);
		brmul(AA,X2T,order_num*(n-1),order_num,N,AX2T);
		brmul(BB,X2T,order_num,order_num,N,BX2T);
		brmul(AX2T,e1,order_num*(n-1),N,1,sita_temp);
		for(i=0;i<order_num*(n-1);i++)
		{
			sita1[i]=sita1[i]-sita_temp[i];
		}
		brmul(BX2T,e1,order_num,N,1,sita2);
		memcpy(&sita1[order_num*(n-1)],sita2,order_num*sizeof(double));
		
		//把X2补充到X1里，形成新的X1矩阵，重新计算增加阶次后的参数估计误差
		memcpy(&X1T[order_num*(n-1)*N],X2T,order_num*N*sizeof(double));
		
		//将X1T矩阵转置得X1矩阵。
		free(X1);
		X1=(double*)calloc(N*order_num*n,sizeof(double));
		for(i=0;i<order_num*n;i++)
		{
			for(j=0;j<N;j++)
			{
				X1[j*order_num*n+i]=X1T[i*N+j];		//将X1T矩阵转置得X1矩阵。
			}
		}

		brmul(X1,sita1,N,order_num*n,1,Ye);
		for(i=0;i<N;i++)
		{
			e1[i]=Y[i]-Ye[i];
		}
		brmul(e1,e1,1,N,1,sigama);

		
		//利用AIC法得模型参数估计值及模型阶次。
		AIC[n-1]=N*log10(sigama[0]/N)+order_num*(order_num*nabc);

		if(AIC[n-1]<AIC0)
		{
			AIC0=AIC[n-1];
			nabc0=nabc;
			memcpy(sita0,sita1,order_num*n*sizeof(double));
			///////////
			if(order_num==2)
			{
				WriteMatrix2File("e_CARMA",e1,N,1);//生成CARMA模型用的白噪声e
			}
			for(i=10;i<N+10;i++)
			{
				ZZZ[i*3+order_num-1]=(int)(10*Ye[i-10]);
				if(!model)ZZZ[i*3+order_num-1]*=4;
				zzz[i*3+order_num-1]=Ye[i-10];
			}
			if(order_num==3)
				WriteMatrix2File("ZZZ",zzz,m_nwnLEN,3);
			/////////////////////////////
		}
		else
		{
			break;
		}
		
		//////////////////////////////////////////
		//重新构造XTX的逆阵，提供下一次递推用。
		XTX_inv=(double*)calloc(order_num*n*order_num*n,sizeof(double));
		CC=(double*)calloc(order_num*(n-1)*order_num*(n-1),sizeof(double));
		temp1=(double*)calloc(order_num*(n-1)*order_num*(n-1),sizeof(double));
		temp2=(double*)calloc(order_num*(n-1)*order_num*(n-1),sizeof(double));
		brmul(AA,X2TX1,order_num*(n-1),order_num,order_num*(n-1),temp1);
		brmul(temp1,X1TX1,order_num*(n-1),order_num*(n-1),order_num*(n-1),temp2);
		for(i=0;i<order_num*(n-1);i++)
		{
			for(j=0;j<order_num*(n-1);j++)
			{
				CC[i*order_num*(n-1)+j]=X1TX1[i*order_num*(n-1)+j]+temp2[i*order_num*(n-1)+j];
			}
		}
		for(i=0;i<order_num*(n-1);i++)
		{
			for(j=0;j<order_num*(n-1);j++)
			{
				XTX_inv[i*order_num*n+j]=CC[i*order_num*(n-1)+j];//C+AX2TX1C
			}
		}

		for(i=0;i<order_num;i++)
		{
			for(j=0;j<order_num*(n-1);j++)
			{
				XTX_inv[(order_num*(n-1)+i)*order_num*n+j]=-AA[order_num*j+i];//-A'
			}
		}
		for(i=0;i<order_num*(n-1);i++)
		{
			for(j=order_num*(n-1);j<order_num*n;j++)
			{
				XTX_inv[i*order_num*n+j]=-AA[i*order_num+j-order_num*(n-1)];//-A
			}
		}

		for(i=order_num*(n-1);i<order_num*n;i++)
		{
			for(j=order_num*(n-1);j<order_num*n;j++)
			{
				XTX_inv[i*order_num*n+j]=BB[(i-order_num*(n-1))*order_num+j-order_num*(n-1)];//B
			}
		}
		//
		free(temp1);
		free(temp2);
		free(CC);
		//
		free(X1TX1);
		X1TX1=(double*)calloc(order_num*n*order_num*n,sizeof(double));
		memcpy(X1TX1,XTX_inv,order_num*n*order_num*n*sizeof(double));
		////////////////////////////////////////////////
		free(X1TX2);
		free(X2TX1);
		free(B1);
		free(AA);
		free(AX2T);
		free(sita_temp);
		free(XTX_inv);
	}

	free(X1TX1);free(X1);
	free(X1T);free(X2T);free(X2);
	free(Y);free(e1);free(BX2T);
}

void CELS_RO_AICDlg::PrintTextandGraph(int order_num)
{
	FILE *fp4;
	CString filename;
	if(model)
		filename.Format("_模型一");
	else
		filename.Format("_模型二");
	if(order_num==2)
		filename.Format(filename+"_CAR");
	else
		filename.Format(filename+"_CARMA");
	filename.Format("sita"+filename+".txt");
    if((fp4=fopen(filename,"w"))==NULL) return;
	fprintf(fp4,"AIC最终辨识结果为：\n");
	fprintf(fp4,"\nna=%d",nabc0);
	fprintf(fp4,"\nnb=%d",nabc0);
	if(order_num==2)
	{
		fprintf(fp4,"\n\n");
		fprintf(fp4,"       a[i]              b[i]\n");
	}
	else
	{
		fprintf(fp4,"\nnc=%d\n\n",nabc0);
		fprintf(fp4,"       a[i]              b[i]              c[i]\n");
	}
	for(int i=0;i<nabc0;i++)
	{
		fprintf(fp4,"a[%d]=%10.6f   ",i+1,sita0[i*order_num]);
		fprintf(fp4,"b[%d]=%10.6f   ",i+1,sita0[i*order_num+1]);
		if(order_num==3)
		{
			fprintf(fp4,"c[%d]=%10.6f\n",i+1,sita0[i*order_num+2]);
		}
		else
		{
			fprintf(fp4,"\n");
		}
	}
	fprintf(fp4,"\n\nAIC值：\n"); 
	for(i=0;i<nabc0;i++)
	{
		fprintf(fp4,"AIC%d=   %10.6f\n",i+1,AIC[i]);
	}
	fclose(fp4);
	if(order_num==3)
	{
		int k;
		if(model)
		{
			//model==TRUE,使用模型一
			k=2;
			sita1[0]=-1.5;sita1[1]= 1.0;sita1[2]=-1.0;
			sita1[3]= 0.7;sita1[4]= 0.5;sita1[5]= 0.2;
			for(i=6;i<30;i++)
			{
				sita1[i]=0;
			}
		}
		else
		{
			//model==TRUE,使用模型一
			k=5;
			sita1[0]=-1.185;sita1[3]= 0.814;sita1[6]=-0.518;sita1[ 9]= 0.349;sita1[12]=-0.117;
			sita1[1]= 1.08 ;sita1[4]=-0.745;sita1[7]= 0.475;sita1[10]=-0.253;sita1[13]= 0.123;
			sita1[2]=-1.0  ;sita1[5]= 0.2  ;sita1[8]= 0.0  ;sita1[11]= 0.0  ;sita1[14]= 0.0  ;
		}
		CString cs;
		cs.Format("");
		for(int j=0;j<order_num;j++)
		{
			for(i=0;i<k;i++)
			{
				if(j==0)
					cs.Format(cs+"a%d=%.3f\n",i+1,sita1[j+3*i]);
				if(j==1)
					cs.Format(cs+"b%d=%.3f\n",i+1,sita1[j+3*i]);
				if(j==2)
				{
					if(i==2)break;
					cs.Format(cs+"c%d=%.3f\n",i+1,sita1[j+3*i]);
				}
			}
		}
		m_ctrStaticRefText.SetDrawText(cs);
		cs.Format("");
		for(j=0;j<order_num;j++)
		{
			for(i=0;i<nabc0;i++)
			{
				if(j==0)
					cs.Format(cs+"a%d=%.6f\n",i+1,sita0[j+order_num*i]);
				if(j==1)
					cs.Format(cs+"b%d=%.6f\n",i+1,sita0[j+order_num*i]);
				if(j==2)
					cs.Format(cs+"c%d=%.6f\n",i+1,sita0[j+order_num*i]);
			}
		}
		m_ctrStaticIDRefText.SetDrawText(cs);
		cs.Format("");
		for(j=0;j<order_num;j++)
		{
			for(i=0;i<nabc0;i++)
			{
				if(j==0)
					cs.Format(cs+"a%d: %.6f\n",i+1,fabs(sita0[j+order_num*i]-sita1[j+3*i]));
				if(j==1)
					cs.Format(cs+"b%d: %.6f\n",i+1,fabs(sita0[j+order_num*i]-sita1[j+3*i]));
				if(j==2)
					cs.Format(cs+"c%d: %.6f\n",i+1,fabs(sita0[j+order_num*i]-sita1[j+3*i]));
			}
		}
		m_ctrStaticErrText.SetDrawText(cs);
		delete m_tempInt;
		m_tempInt=new int[nabc0];
		for(i=0;i<nabc0;i++)
		{
			m_tempInt[i]=int((600+AIC[i]));
		}
		m_ctrStaticGraph.SetDrawData(m_nwnLEN,ZZZ,3,1);
		m_ctrStaticGraph.RedrawWindow();
		RedrawWindow();
	}

}



void CELS_RO_AICDlg::OnChangeEDITmVALn() 
{
	// TODO: If this is a RICHEDIT control, the control will not
	// send this notification unless you override the CDialog::OnInitDialog()
	// function and call CRichEditCtrl().SetEventMask()
	// with the ENM_CHANGE flag ORed into the mask.
	
	// TODO: Add your control notification handler code here
	
}

void CELS_RO_AICDlg::OnDestroy() 
{
	CDialog::OnDestroy();
	
	// TODO: Add your message handler code here
	free(ZZZ);	
	
}



void CELS_RO_AICDlg::OnRadioModel1() 
{
	// TODO: Add your control notification handler code here
	model=TRUE;
}

void CELS_RO_AICDlg::OnRadioModel2() 
{
	// TODO: Add your control notification handler code here
	model=FALSE;
}

void CELS_RO_AICDlg::OnLButtonDblClk(UINT nFlags, CPoint point) 
{
	// TODO: Add your message handler code here and/or call default
	m_nIndex++;
	int *tempZ;
	tempZ=(int*)calloc(3*m_nwnLEN,sizeof(int));
	for(int i=0;i<3;i++)
	{
		for(int j=0;j<m_nwnLEN;j++)
		{
			if(fmod(m_nIndex,4)==i)break;
			tempZ[i+j*3]=ZZZ[i+j*3];
		}
	}
	m_ctrStaticGraph.SetDrawData(m_nwnLEN,tempZ,3,1);
	m_ctrStaticGraph.RedrawWindow();
	RedrawWindow();
	if(m_nIndex==11&&model)
		MessageBox("  ");
	CDialog::OnLButtonDblClk(nFlags, point);
}