aicdlg.cpp

可用于同时辨识模型阶次和参数

// AICDlg.cpp : implementation file
//

#include "stdafx.h"
#include "AIC.h"
#include "AICDlg.h"

#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

#define A 179
#define PAI 3.1415926

/////////////////////////////////////////////////////////////////////////////
// CAboutDlg dialog used for App About

class CAboutDlg : public CDialog
{
public:
	CAboutDlg();

// Dialog Data
	//{{AFX_DATA(CAboutDlg)
	enum { IDD = IDD_ABOUTBOX };
	//}}AFX_DATA

	// ClassWizard generated virtual function overrides
	//{{AFX_VIRTUAL(CAboutDlg)
	protected:
	virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV support
	//}}AFX_VIRTUAL

// Implementation
protected:
	//{{AFX_MSG(CAboutDlg)
	//}}AFX_MSG
	DECLARE_MESSAGE_MAP()
};

CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)
{
	//{{AFX_DATA_INIT(CAboutDlg)
	//}}AFX_DATA_INIT
}

void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CAboutDlg)
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)
	//{{AFX_MSG_MAP(CAboutDlg)
		// No message handlers
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CAICDlg dialog

CAICDlg::CAICDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CAICDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CAICDlg)
	m_nN = 600;
	m_nMLen = 1023;
	m_nMLoc = 63;
	m_nMn = 10;
	m_dNmea = 0.0;
	m_dNvar = 0.1;
	m_nIDMethods = 0;
	//}}AFX_DATA_INIT
	// Note that LoadIcon does not require a subsequent DestroyIcon in Win32
	m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void CAICDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CAICDlg)
	DDX_Control(pDX, IDC_EDIT_Mn, m_editMn);
	DDX_Control(pDX, IDC_EDIT_MLoc, m_editMLoc);
	DDX_Control(pDX, IDC_EDIT_DataN, m_ctrN);
	DDX_Control(pDX, IDC_SPIN_Mn, m_ctrMn);
	DDX_Control(pDX, IDC_SPIN_Mloc, m_ctrMLoc);
	DDX_Text(pDX, IDC_EDIT_DataN, m_nN);
	DDX_Text(pDX, IDC_EDIT_Mlen, m_nMLen);
	DDX_Text(pDX, IDC_EDIT_MLoc, m_nMLoc);
	DDX_Text(pDX, IDC_EDIT_Mn, m_nMn);
	DDX_Text(pDX, IDC_EDIT_Nmea, m_dNmea);
	DDX_Text(pDX, IDC_EDIT_Nvar, m_dNvar);
	DDX_CBIndex(pDX, IDC_COMBO_IDMethods, m_nIDMethods);
	DDX_Control(pDX, IDC_STATIC_REFm, m_cStaticTextREFm);
	DDX_Control(pDX, IDC_STATIC_REFi, m_cStaticTextREFi);
	DDX_Control(pDX, IDC_STATIC_REFe, m_cStaticTextREFe);
	DDX_Control(pDX, IDC_STATIC_GRAPH, m_cStaticGraph);
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CAICDlg, CDialog)
	//{{AFX_MSG_MAP(CAICDlg)
	ON_WM_SYSCOMMAND()
	ON_WM_PAINT()
	ON_WM_QUERYDRAGICON()
	ON_BN_CLICKED(IDC_BUTTON_OK, OnButtonOk)
	ON_EN_CHANGE(IDC_EDIT_Mn, OnChangeEDITMn)
	ON_EN_CHANGE(IDC_EDIT_DataN, OnChangeEDITDataN)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CAICDlg message handlers

BOOL CAICDlg::OnInitDialog()
{
	CDialog::OnInitDialog();

	// Add "About..." menu item to system menu.

	// IDM_ABOUTBOX must be in the system command range.
	ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
	ASSERT(IDM_ABOUTBOX < 0xF000);

	CMenu* pSysMenu = GetSystemMenu(FALSE);
	if (pSysMenu != NULL)
	{
		CString strAboutMenu;
		strAboutMenu.LoadString(IDS_ABOUTBOX);
		if (!strAboutMenu.IsEmpty())
		{
			pSysMenu->AppendMenu(MF_SEPARATOR);
			pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
		}
	}

	// Set the icon for this dialog.  The framework does this automatically
	//  when the application's main window is not a dialog
	SetIcon(m_hIcon, TRUE);			// Set big icon
	SetIcon(m_hIcon, FALSE);		// Set small icon
	
	// TODO: Add extra initialization here
	m_ctrMLoc.SetBuddy(&m_editMLoc);
	m_ctrMLoc.SetRange(1,m_nMLen-m_nN);
	m_ctrMn.SetBuddy(&m_editMn);
	m_ctrMn.SetRange(1,13);
	m_tempInt=NULL;
	return TRUE;  // return TRUE  unless you set the focus to a control
}

void CAICDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
	if ((nID & 0xFFF0) == IDM_ABOUTBOX)
	{
		CAboutDlg dlgAbout;
		dlgAbout.DoModal();
	}
	else
	{
		CDialog::OnSysCommand(nID, lParam);
	}
}

// If you add a minimize button to your dialog, you will need the code below
//  to draw the icon.  For MFC applications using the document/view model,
//  this is automatically done for you by the framework.

void CAICDlg::OnPaint() 
{
	if (IsIconic())
	{
		CPaintDC dc(this); // device context for painting

		SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);

		// Center icon in client rectangle
		int cxIcon = GetSystemMetrics(SM_CXICON);
		int cyIcon = GetSystemMetrics(SM_CYICON);
		CRect rect;
		GetClientRect(&rect);
		int x = (rect.Width() - cxIcon + 1) / 2;
		int y = (rect.Height() - cyIcon + 1) / 2;

		// Draw the icon
		dc.DrawIcon(x, y, m_hIcon);
	}
	else
	{
		CDialog::OnPaint();
	}
}

// The system calls this to obtain the cursor to display while the user drags
//  the minimized window.
HCURSOR CAICDlg::OnQueryDragIcon()
{
	return (HCURSOR) m_hIcon;
}

void CAICDlg::aic(int n, int lu)
{
	//结合AIC的模型阶次递推算法
	//本程序最大能算到10阶。
	int N=n-10;
	FILE *fp1,*fp2,*fp4;
    int i,j,na,nb,na0,nb0;
	double AIC0,AIC[10],sigama[1];
	double *z,*u,*e,*X1T,*Y,*Ye,*e1,*X1,*X2T,*X2,*XTX_inv;
	double sita0[20],sita1[20],sita2[2],*sita_temp;
	double *X1TX1,*X1TY,X2TX2[4],*X2TX1,*X1TX2;
	double *B1,*AA,B2[4],BB[4],*AX2T,*BX2T;
	double *CC,*temp1,*temp2;
	na=1;nb=1;
	z=(double*)calloc(N+10,sizeof(double));
	u=(double*)calloc(N+10+lu,sizeof(double));
	e=(double*)calloc(N+10,sizeof(double));
	X1T=(double*)calloc(20*N,sizeof(double));
	X2T=(double*)calloc(2*N,sizeof(double));
	X2=(double*)calloc(N*2,sizeof(double));
	Y=(double*)calloc(N,sizeof(double));
	Ye=(double*)calloc(N,sizeof(double));
	e1=(double*)calloc(N,sizeof(double));
	BX2T=(double*)calloc(2*N,sizeof(double));

	if((fp1=fopen("MSerials.txt","r"))==NULL)return;
    if((fp2=fopen("WhiteNoise2.txt","r"))==NULL)return; 
 	//读入输入u[]和白噪声e[]
    for(i=0;i<N+10+lu;i++)
    { 
		fscanf(fp1,"%lf ",&u[i]);
	}
    for(i=0;i<N+10;i++)
    { 
		u[i]=u[i+lu];
	}
    for(i=0;i<N+10;i++)
    { 
		fscanf(fp2,"%lf ",&e[i]);
	}
	fclose(fp1);
    fclose(fp2);
	
	//用下式计算理论输出z[]
    z[0]=e[0];
	z[1]=1.185*z[0]+1.08*u[0]+e[1];
	z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+e[2];
	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];
	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];
    for(i=5;i<N+10;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];
	}
	
	WriteData2File("u_aic",u,N,1);
	WriteData2File("e_aic",e,N,1);
	WriteData2File("z_aic",z,N,1);
	///////////////
	for(i=0;i<N;i++)
	{
		Y[i]=z[10+i];
	}
	//na=nb=1;
	//赋值给X1T矩阵
	for(i=0;i<2;i++)
	{
		for(j=0;j<N;j++)
		{
			if(i==0)
			{
				X1T[i*N+j]=-z[9+j];
			}
			else
			{
				X1T[i*N+j]=u[9+j];
			}
		}
	}
	//将X1T矩阵转置得X1矩阵。
	X1=(double*)calloc(N*2,sizeof(double));
	for(i=0;i<2;i++)
	{
		for(j=0;j<N;j++)
		{
			X1[j*2+i]=X1T[i*N+j];
		}
	}
	//计算X1TX1，并求逆。
	X1TX1=(double*)calloc(2*2,sizeof(double));
	brmul(X1T,X1,2,N,2,X1TX1);
	brinv(X1TX1,2);
	//求参数估计值sita
	X1TY=(double*)calloc(2,sizeof(double));
	brmul(X1T,Y,2,N,1,X1TY);
	brmul(X1TX1,X1TY,2,2,1,sita1);
	free(X1TY);
	//求预测误差矩阵e1。
	brmul(X1,sita1,N,2,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)+2*(na+nb)+N;
	AIC[0]=AIC0;
	na0=1;nb0=1;
	memcpy(sita0,sita1,2*sizeof(double));
	//
	//na=nb=n;   X=[X1，X2]
	for(n=2;n<=10;n++)
	{	
		na=n;nb=n;
		//模型阶次增加，重新赋值给X2T矩阵
		for(i=0;i<2;i++)
		{
			for(j=0;j<N;j++)
			{
				if(i==0)
				{
					X2T[i*N+j]=-z[10-n+j];
				}
				else
				{
					X2T[i*N+j]=u[10-n+j];
				}
			}
		}
		//X2T矩阵转置得X2矩阵
		for(i=0;i<2;i++)
		{
			for(j=0;j<N;j++)
			{
				X2[j*2+i]=X2T[i*N+j];
			}
		}

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

		brmul(X1,sita1,N,2*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)+2*(na+nb)+N;
		if(AIC[n-1]<AIC0)
		{
			AIC0=AIC[n-1];
			na0=na;nb0=nb;
			memcpy(sita0,sita1,2*n*sizeof(double));
		}
		else
		{
			break;
		}
		//////////////////////////////////////////
		//重新构造XTX的逆阵，提供下一次递推用。
		XTX_inv=(double*)calloc(2*n*2*n,sizeof(double));
		CC=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
		temp1=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
		temp2=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
		brmul(AA,X2TX1,2*(n-1),2,2*(n-1),temp1);
		brmul(temp1,X1TX1,2*(n-1),2*(n-1),2*(n-1),temp2);
		for(i=0;i<2*(n-1);i++)
		{
			for(j=0;j<2*(n-1);j++)
			{
				CC[i*2*(n-1)+j]=X1TX1[i*2*(n-1)+j]+temp2[i*2*(n-1)+j];
			}
		}
		for(i=0;i<2*(n-1);i++)
		{
			for(j=0;j<2*(n-1);j++)
			{
				XTX_inv[i*2*n+j]=CC[i*2*(n-1)+j];
			}
		}
		for(j=0;j<2*(n-1);j++)
		{
			XTX_inv[2*(n-1)*2*n+j]=-AA[2*j];
			XTX_inv[(2*n-1)*2*n+j]=-AA[2*j+1];
		}
		for(i=0;i<2*(n-1);i++)
		{
			for(j=2*(n-1);j<2*n;j++)
			{
				XTX_inv[i*2*n+j]=-AA[i*2+j-2*(n-1)];
			}
		}
		for(i=2*(n-1);i<2*n;i++)
		{
			for(j=2*(n-1);j<2*n;j++)
			{
				XTX_inv[i*2*n+j]=BB[(i-2*(n-1))*2+j-2*(n-1)];
			}
		}
		//
		free(temp1);
		free(temp2);
		free(CC);
		//
		free(X1TX1);
		X1TX1=(double*)calloc(2*n*2*n,sizeof(double));
		memcpy(X1TX1,XTX_inv,2*n*2*n*sizeof(double));
		////////////////////////////////////////////////
		free(X1TX2);
		free(X2TX1);
		free(B1);
		free(AA);
		free(AX2T);
		free(sita_temp);
		free(XTX_inv);
	}

	free(X1TX1);free(X1);
	free(z);free(u);free(e);free(X1T);free(X2T);free(X2);
	free(Y);free(Ye);free(e1);free(BX2T);

    if((fp4=fopen("sita_aic.txt","w"))==NULL) return;
	fprintf(fp4,"AIC最终辨识结果sita[]为：\n");
	fprintf(fp4,"\nna=%d\n",na0);
	fprintf(fp4,"nb=%d\n\n",nb0);
	fprintf(fp4,"       a[i]              b[i]\n");
	for(i=0;i<na0;i++)
	{
		fprintf(fp4,"a[%d]=%10.6f   ",i+1,sita0[i*2]);
		fprintf(fp4,"b[%d]=%10.6f\n",i+1,sita0[i*2+1]);
	}
	fprintf(fp4,"\n"); 
	for(i=0;i<10;i++)
	{
		fprintf(fp4,"AIC%d=   %10.6f\n",i+1,AIC[i]);
	}
	fclose(fp4);
	sita1[0]=-1.185;sita1[2]= 0.814;sita1[4]=-0.518;sita1[6]= 0.349;sita1[8]=-0.117;
	sita1[1]= 1.08 ;sita1[3]=-0.745;sita1[5]= 0.475;sita1[7]=-0.253;sita1[9]= 0.123;
	CString cs;
	cs.Format("模型参数：");
	for(j=0;j<2;j++)
	{
		for(i=0;i<na0;i++)
		{
			if(j==0)
			{
				cs.Format(cs+"\na%d=%.3f",i+1,sita1[j+2*i]);