lsm.cpp

改进的最小二乘法实现代码

//最小二乘法类cpp文件
//Edit by Wang Shimin
//2004.8.1～～2004.8.24
//本文件为最新文件，修改于2006.3.12

// LSM.cpp: implementation of the LSM class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "LSM.h"
#include "math.h"

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

LSM::LSM()
{
	m_stepcoefficient = 1.0; /* 步长因子初始化为1.0 */ 
	m_oldfitness = 0.0;   /* 上一步适应值初始化为零 */
	m_precisionnum = 0;   /* 连续不优化次数初始化为零 */
	m_fitness = 0.0;      /* 本步适应值初始化为零 */
	p_Bndload = NULL;     /* 边界荷载 */
	m_errorflag = TRUE;   /* 收敛标志变量 */
	m_numlap = 1;     /* 循环次数初始化为零 */
	p_matrix4 =NULL;   /* 存储系数阵的转置矩阵与目标矩阵的乘积矩阵 */
	p_matrix3 =NULL;   /* 存储系数阵与其转置矩阵的乘积矩阵的逆阵 */
	p_matrix2 =NULL;   /* 存储系数阵与其转置矩阵的乘积矩阵 */
	p_matrix1 = NULL;  /* 存储系数阵 */
	p_objvalue = NULL;  /* 存储目标值，计算值与观测值之间的差值 */
	p_calculation = NULL;  /* 存储计算值 */
	p_observation = NULL;  /* 存储观测值 */
	p_variable = NULL;   /* 存储反演参数 */
	p_deltavariable = NULL; /* 存储反演参数增量 */
	p_BAKindflag = NULL;/* 反演种类项标志 */
	p_iparameter = NULL;/* 反演材料种类项标志 */
	p_jparameter = NULL;/* 反演材料性质项标志 */
	p_dampcoefficient = NULL;  /* 阻尼因子 */

	m_blErrorExit = false;
}

LSM::~LSM()
{
	/* 释放动态内存 */
	delete[]p_calculation;
	delete[]p_observation;
	delete[]p_objvalue;
	delete[]p_variable;
	delete[]p_deltavariable;
	delete[]p_dampcoefficient;
	
	for(int i=0; i<m_numvar; i++)
	{
		delete[]p_matrix1[i];
		delete[]p_matrix2[i];
		delete[]p_matrix3[i];
	}

	delete[]p_matrix4;
	delete[]p_matrix3;
	delete[]p_matrix2;
	delete[]p_matrix1;

	delete[]p_BAKindflag;
	delete[]p_iparameter;
	delete[]p_jparameter;
	delete[]p_Bndload;
}

/* 初始化程序中的参数 */
void LSM::Initialization()
{
	CString csBaaFileName = mGcsFileName+_T(".BAA");
	ifstream cBaaFile;
	cBaaFile.open((LPCSTR)csBaaFileName); 
	if (cBaaFile.fail())
	{
		AfxMessageBox("Cannot open input BAA file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	
	int itemp;
	cBaaFile >> itemp;
	cBaaFile >> m_numvar;   /* 待反演参数的个数 */
	cBaaFile.close();


	CString csMEAFileName = mGcsFileName+_T(".MEA");
	ifstream cMEAFile;
	cMEAFile.open((LPCSTR)csMEAFileName); 
	if (cMEAFile.fail())
	{
		AfxMessageBox("Cannot open input MEA file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	
	cMEAFile >> m_numpoint;  /* 从MEA文件中读取观测点的个数 */
	cMEAFile.close();

	if(p_variable)
		delete[]p_variable;
	p_variable = new double[m_numvar];   /* 存储反演参数 */

	if(p_deltavariable)
		delete[]p_deltavariable;
	p_deltavariable = new double[m_numvar];  /* 存储反演参数增量 */

	if(p_calculation)
		delete[]p_calculation;
	p_calculation = new double[m_numpoint];  /* 存储计算值 */

	if(p_observation)
		delete[]p_observation;
	p_observation = new double[m_numpoint];  /* 存储观测值 */

	if(p_objvalue)
		delete[]p_objvalue;
	p_objvalue = new double[m_numpoint];  /* 存储目标值，即观测值与计算值之差 */

	/* 为存储识别待反演参数性质的数组分配内存空间 */

	if(p_BAKindflag)
		delete[]p_BAKindflag;

	if(p_iparameter)
		delete[]p_iparameter;

	if(p_jparameter)
		delete[]p_jparameter;
		
	if(p_dampcoefficient)
		delete[]p_dampcoefficient;

	p_BAKindflag = new int[m_numvar];
	p_iparameter = new int[m_numvar];
	p_jparameter = new int[m_numvar];
	p_dampcoefficient = new double[m_numvar];

	int i;
	for(i=0; i<m_numvar; i++)
		p_dampcoefficient[i] = 0.00000000000001; /* 阻尼因子 */

	/* 读入待反演参数 */
	ReadBAA();
	if(m_blErrorExit) return;
	/* 读入算法参数 */
	ReadINP();
	if(m_blErrorExit) return;

	/* 读入测点的观测值 */
	ReadMEA();
	if(m_blErrorExit) return;

	m_delta = 0.01; /* 差分精度 */
	
	/* 将待反演参数的增量值初始化为零 */
	for(i=0; i<m_numvar; i++)
		p_deltavariable[i] = 0.0;

	/* 将参数存储到文件 */
	WriteToFem(p_variable);
	if(m_blErrorExit) return;
	/* 有限元计算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 读入测点的计算值 */
	ReadFromFem(p_calculation);
	if(m_blErrorExit) return;

	/* 计算测点的观测值与计算值之差 ，目标矩阵*/
	for(i=0; i<m_numpoint; i++)
		p_objvalue[i] = p_observation[i] - p_calculation[i];
	
	for(i=0; i<m_numpoint; i++)
	{
		m_fitness += pow((p_objvalue[i]),2.0);
	}

	/* 为矩阵分配动态内存空间 */
	if(p_matrix1)
		delete[]p_matrix1;
	p_matrix1 = new double *[m_numvar];
	
	for(i=0; i<m_numvar; i++)
		p_matrix1[i]=new double[m_numpoint];
	
	if(p_matrix2)
		delete[]p_matrix2;
	p_matrix2 = new double *[m_numvar];
	
	for(i=0; i<m_numvar; i++)
		p_matrix2[i]=new double[m_numvar];

	if(p_matrix3)
		delete[]p_matrix3;
	p_matrix3 = new double *[m_numvar];
	
	for(i=0; i<m_numvar; i++)
		p_matrix3[i]=new double[m_numvar];

	if(p_matrix4)
		delete[]p_matrix4;
	p_matrix4 = new double[m_numvar];
}

/* 从.INP文件中读取最大迭代次数和两个控制迭代的精度参数 */
void LSM::ReadINP()
{
	CString csINPFileName = mGcsFileName+_T(".INP");
	ifstream cINPFile;
	cINPFile.open((LPCSTR)csINPFileName); 
	if (cINPFile.fail())
	{
		AfxMessageBox("Cannot open input INP file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}

	char tempch[256];

	cINPFile.getline(tempch,256);
	cINPFile >> m_maxnumlap >> m_precision1 >> m_precision2;

	//从用户输入读取每个待反演参数的阻尼系数
//    for(int i=0; i<m_numvar; i++)
//		cINPFile >> p_dampcoefficient[i];

	cINPFile.close();
}

/* 从BAA文件中读入待反演参数的性质标志以及初始值 */
void LSM::ReadBAA()
{
	char tempchar[256];
	int  index;

	CString csBaaFileName = mGcsFileName+_T(".BAA");
	ifstream cBaaFile;
	cBaaFile.open((LPCSTR)csBaaFileName); 
	if (cBaaFile.fail())
	{
		AfxMessageBox("Cannot open input BAA file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	
	cBaaFile.getline(tempchar,256);

	//读入各个参数的上下限，并计算出各个参数的精度
	for(int i=0; i<m_numvar; i++)
	{
		cBaaFile >> index;
		cBaaFile >> p_BAKindflag[i];
		cBaaFile >> p_iparameter[i];
		cBaaFile >> p_jparameter[i];
		cBaaFile >> p_variable[i];
		cBaaFile.getline(tempchar,256);
	}

	cBaaFile >> m_numbndload;

	if(m_numbndload>0)
	{
		if(p_Bndload)
		{
			delete[]p_Bndload ;
			p_Bndload = NULL ;
		}

		p_Bndload = new bndload[m_numbndload];

		for(int j=0; j<m_numbndload; j++)
		{
			cBaaFile >> p_Bndload[j].index >> p_Bndload[j].wAddStage
				>> p_Bndload[j].wAddStep >> p_Bndload[j].iNode1
				>> p_Bndload[j].iNode2 >> p_Bndload[j].iCoordFlag
				>> p_Bndload[j].BeginLoad1 >> p_Bndload[j].BeginLoad2
				>> p_Bndload[j].BeginLoad3 >> p_Bndload[j].EndLoad1
				>> p_Bndload[j].EndLoad2 >> p_Bndload[j].EndLoad3;
		
			cBaaFile.getline(tempchar,256);
		}
	}
	
	cBaaFile.close();
}

//从MEA文件中读取监测数据
void LSM::ReadMEA()
{
	CString csFromMeaFileName = mGcsFileName+_T(".MEA");//file from fem
	ifstream cFromMeaFile;
	cFromMeaFile.open((LPCSTR)csFromMeaFileName); 
	if (cFromMeaFile.fail())
	{
		AfxMessageBox("Cannot open input file of mea!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}

	int iTemp;
	double dTemp;
	char tempch[256];

	cFromMeaFile.getline(tempch,256);
	for(int i=0;i<m_numpoint;i++)
	{
		cFromMeaFile >> iTemp >> dTemp >> dTemp >> dTemp  >> iTemp  >> iTemp
			 >> iTemp  >> iTemp  >> iTemp  >> iTemp  >> iTemp;
		cFromMeaFile >> p_observation[i];
		cFromMeaFile.getline(tempch,256);
	}

	cFromMeaFile.close();
}

/* 计算反演参数的增量值 */
void LSM::Calculate_Delta()
{
	int i,j,k;

	/* 计算系数矩阵p_matrix1 */
	CoefficientMatrix();

    /* 计算系数矩阵与其转置矩阵的乘积矩阵，即p_matrix2 */
	for(i=0; i<m_numvar; i++)
	{
		for(j=0; j<m_numvar; j++)
		{
			p_matrix2[i][j] = 0.0;

			for(k=0; k<m_numpoint; k++)
				p_matrix2[i][j] += p_matrix1[i][k]*p_matrix1[j][k];
		}
	}

	/* 计算矩阵p_matrix2的转置矩阵，存储在p_matrix3中，在此计算过程中,p_matrix2变成单位阵 */
	ReverseMatrix(p_matrix2,p_matrix3,m_numvar);
	if(m_blErrorExit)return;
	/* 计算系数阵与目标值矩阵的乘积矩阵p_matrix4 */
	for(i=0; i<m_numvar; i++)
	{
		p_matrix4[i] = 0.0;

		for(j=0; j<m_numpoint; j++)
			p_matrix4[i] += p_matrix1[i][j]*p_objvalue[j];
	}
	
	/* 计算待反演参数增量值 */
	for(i=0; i<m_numvar; i++)
	{
		p_deltavariable[i] = 0.0;

		for(j=0; j<m_numvar; j++)
			p_deltavariable[i] += p_matrix3[i][j]*p_matrix4[j];
	}
}

/* 计算系数矩阵 */
void LSM::CoefficientMatrix()
{
	/* 临时存储新参数 */
	double *newvariable;
	/* 临时存储增量的有限元计算值 */
	double *newcalculation;
    /* 求偏导时临时存储两次有限元计算值的差值 */
	double *daltevalue; 
	
	/* 为指针变量开辟动态内存空间 */
	newvariable = new double[m_numvar];
	newcalculation = new double[m_numpoint];
	daltevalue = new double[m_numpoint];
	
	for(int i=0; i<m_numvar; i++)
		newvariable[i] = p_variable[i];
	
	for(i=0; i<m_numvar; i++)
	{
		/* 计算新参数 */
		newvariable[i] = p_variable[i]*(1.0+m_delta);
		/* 将新参数存储到文件 */
		WriteToFem(newvariable);
		if(m_blErrorExit) return;
		/* 有限元计算 */
		ObjFunction();
		if(m_blErrorExit) return;
		/* 从文件中读取有限元计算结果 */
		ReadFromFem(newcalculation);
		if(m_blErrorExit) return;

		newvariable[i] = p_variable[i];