lsm.cpp

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

		/* 利用差分法计算偏导数 */
		for(int j=0; j<m_numpoint; j++)
		{
			daltevalue[j] = newcalculation[j] - p_calculation[j];
			p_matrix1[i][j] = daltevalue[j]/(p_variable[i]*m_delta);
		}
	}

	/* 释放临时内存空间 */
	delete[]newvariable;
	delete[]newcalculation;
    delete[]daltevalue; 
}

/* 计算矩阵的逆阵 */
void LSM::ReverseMatrix(double ** inmatrix,double ** outmatrix,int number)
{
	//function://高斯——约当法求逆阵*****
	//Parameter:inmatrix指向输入矩阵的地址，outmatrix是逆阵．number是矩阵的维数

	int i,j,k,l;
	double flag;  //用于判断矩阵是否正定的标志变量
	double temp=0.0;

    /* 改善矩阵的条件数，增加解的稳定性 */
	for(i=0; i<number; i++)
		for(j=0; j<number; j++)
			if(i==j)
				inmatrix[i][j] += p_dampcoefficient[i];

   	//*****建立一个单位阵*******//
	
	for(i=0;i<number;i++)
	{
		for(j=0;j<number;j++)
		{
			if(i==j)
				outmatrix[i][j]=1.0;
			else
				outmatrix[i][j]=0.0;
		}
	}

	/* 生成逆阵，存储在输出矩阵outmatrix中，输入矩阵inmatrix变成单位阵 */

	/* 判断是否正定，非正定报错，正定计算逆阵*/
	flag = inmatrix[0][0];

	if((fabs(flag)+1.0)==1.0)
	{
		m_blErrorExit = true;
		AfxMessageBox(GIDS_PROMPT_NOTASSURE_MATRIX);
	}
	else
	{
		for(i=0;i<number;i++)
		{
			for(j=0;j<number;j++)
			{
				if(i!=j)
				{
					temp=inmatrix[j][i]/inmatrix[i][i];			
					for(k=0; k<number; k++)
					{
						inmatrix[j][k]=inmatrix[j][k]-inmatrix[i][k]*temp;
						outmatrix[j][k]=outmatrix[j][k]-outmatrix[i][k]*temp;
					}
				}			
			}

			temp=inmatrix[i][i];

			for(l=0; l<number; l++)
			{
				inmatrix[i][l]=inmatrix[i][l]/temp;
				outmatrix[i][l]=outmatrix[i][l]/temp;
			}
		}
	}
}

/* 评价函数，用于判定是否达到精度要求 */
bool LSM::Criterion()
{
	bool flag;
	
	if((fabs(m_fitness)<=m_precision1)||(m_precisionnum>0))
	{
		m_errorflag = FALSE;
	}

	if((m_numlap<=m_maxnumlap)&&(m_errorflag))
		flag = TRUE;
	else
		flag = FALSE;

	return flag;
}

/* 计算 */
void LSM::Compute()
{
	Initialization();
	if(m_blErrorExit)return;

	while(Criterion()==TRUE)
	{
		m_oldfitness = m_fitness;
		Calculate_Delta();
		if(m_blErrorExit)return;
		AdjustVariable();
		if(m_blErrorExit) return;

		if(fabs((m_oldfitness-m_fitness)/m_oldfitness)<m_precision2)
			m_precisionnum++;
		/* 输出反分析优化结果 */
		WriteResult();
		if(m_blErrorExit) return;

		m_numlap++;
	}
}

/* 将变量存储到文件中 */
void LSM::WriteToFem(double * pvariable)
{
	CString str;

	CString csBackOutFileName = mGcsFileName+_T(".btof");//back analysis to fem
	ofstream cBackOutFile;
	cBackOutFile.open((LPCSTR)csBackOutFileName); 
	if (cBackOutFile.fail())
	{
		AfxMessageBox("Cannot open input file of back analysis!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	//write
    int index=1;

	cBackOutFile << m_numvar << endl;
    
	for(int i=0; i<m_numvar; i++)
	{
		str.Format("%d %d %d %d %.6f\n",index++, p_BAKindflag[i],
			p_iparameter[i], p_jparameter[i], pvariable[i]);
		
		cBackOutFile << (LPCSTR)str;
	}

	//边界荷载输出
	if (m_numbndload > 0)
	{
		index=1;
		cBackOutFile << m_numbndload << endl;
	
		for(int j=0; j<m_numbndload; j++)
		{
			str.Format("%d %d %d %d %d %d %.4f %.4f %.4f %.4f %.4f %.4f\n",
				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);
			
			cBackOutFile << (LPCSTR)str;
		}
	}

	cBackOutFile.close();
}

/* 有限元计算 */
void LSM::ObjFunction()
{
	if(EXCUTE_NUMERICAL_ANALYSIS()==false)
		m_blErrorExit=true;
}

/* 从文件中读取有限元计算结果 */
void LSM::ReadFromFem(double * calculatevalue)
{
	CString csFromFemFileName = mGcsFileName+_T(".ftob");//file from fem
	ifstream cFromFemFile;
	cFromFemFile.open((LPCSTR)csFromFemFileName); 
	if (cFromFemFile.fail())
	{
		AfxMessageBox("Cannot open input file of fem!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}

	int index;
	
	for(int i=0;i<m_numpoint;i++)
	{
		cFromFemFile >> index;
		cFromFemFile >> calculatevalue[i];
	}

	cFromFemFile.close();
}

void LSM::AdjustVariable()
{
	m_stepcoefficient = 1.0;

	int i,flag;
	flag=0;

	double dratio=0.0;
	double newfitness = 0.0;

	/* 临时存储新参数 */
	double *newvariable;
	/* 临时存储增量计算值 */
	double *newcalculation;
    /* 临时存储增量计算值的目标值 */
	double *newobjvalue; 

	newvariable = new double[m_numvar];
	newcalculation = new double[m_numpoint];
	newobjvalue = new double[m_numpoint];
	
	double oldvalue,newvalue;
	oldvalue=0.0;
	newvalue=0.0;

	for(i=0; i<m_numvar; i++)
		newvariable[i] = p_variable[i] + p_deltavariable[i]*m_stepcoefficient;
	
	/* 将新参数存储到文件 */
	WriteToFem(newvariable);
	if(m_blErrorExit) return;
	/* 有限元计算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 从文件中读取有限元计算结果 */
	ReadFromFem(newcalculation);
	if(m_blErrorExit) return;
	
	for(i=0; i<m_numpoint; i++)
		newfitness += pow((p_observation[i]-newcalculation[i]),2.0);

	if(newfitness<m_fitness)
	{
		while(flag<4)
		{
			if(newfitness<m_fitness)
			{
				/*******************保存结果**********************/
				
				m_fitness = newfitness;
				
				for(i=0; i<m_numvar; i++)
				{
					p_variable[i] = newvariable[i];
					p_calculation[i] = newcalculation[i];
				}
				
				/* 更新测点的观测值与计算值之差 ，目标矩阵*/
				for(i=0; i<m_numpoint; i++)
				{
					p_objvalue[i] = p_observation[i] - p_calculation[i];
				}
				
				/*******************保存结果**********************/
				
				for(i=0; i<m_numvar; i++)
					newvariable[i] = p_variable[i] + p_deltavariable[i]*m_stepcoefficient;
				
				/* 将新参数存储到文件 */
				WriteToFem(newvariable);
				if(m_blErrorExit) return;
				/* 有限元计算 */
				ObjFunction();
				if(m_blErrorExit) return;
				/* 从文件中读取有限元计算结果 */
				ReadFromFem(newcalculation);
				if(m_blErrorExit) return;
				
				newfitness=0.0;
				
				for(i=0; i<m_numpoint; i++)
					newfitness += pow((p_observation[i]-newcalculation[i]),2.0);
				
				flag++;
			}
			else
			{
					/* 计算新的增量 */
				Calculate_Delta();
				m_stepcoefficient = -1.0*m_stepcoefficient;

				for(i=0; i<m_numvar; i++)
					newvariable[i] = p_variable[i] + p_deltavariable[i]*m_stepcoefficient;
				
				/* 将新参数存储到文件 */
				WriteToFem(newvariable);
				if(m_blErrorExit) return;
				/* 有限元计算 */
				ObjFunction();
				if(m_blErrorExit) return;
				/* 从文件中读取有限元计算结果 */
				ReadFromFem(newcalculation);
				if(m_blErrorExit) return;
				
				newfitness=0.0;
				
				for(i=0; i<m_numpoint; i++)
					newfitness += pow((p_observation[i]-newcalculation[i]),2.0);
				
				flag++;
			}
		}
	}
	else
	{
		flag = 0;
		while(flag<2)
		{
			if(newfitness>=m_fitness)
			{
				/* 调整步长因子 */
				m_stepcoefficient = -1.0*m_stepcoefficient; 
				
				for(i=0; i<m_numvar; i++)
					newvariable[i] = p_variable[i] + p_deltavariable[i]*m_stepcoefficient;
				
				/* 将新参数存储到文件 */
				WriteToFem(newvariable);
				if(m_blErrorExit) return;
				/* 有限元计算 */
				ObjFunction();
				if(m_blErrorExit) return;
				/* 从文件中读取有限元计算结果 */
				ReadFromFem(newcalculation);
				if(m_blErrorExit) return;
				
				newfitness=0.0;
				
				for(i=0; i<m_numpoint; i++)
					newfitness += pow((p_observation[i]-newcalculation[i]),2.0);
				
				flag++;
			}
			else
			{
				/*******************保存结果**********************/
				m_fitness = newfitness;
				
				for(i=0; i<m_numvar; i++)
				{
					p_variable[i] = newvariable[i];
					p_calculation[i] = newcalculation[i];
				}
				
				/* 更新测点的观测值与计算值之差 ，目标矩阵*/
				for(i=0; i<m_numpoint; i++)
				{
					p_objvalue[i] = p_observation[i] - p_calculation[i];
				}
				/*******************保存结果**********************/
			}
		}
	}
	/* 释放临时内存 */
	delete[]newvariable;
	delete[]newcalculation;
    delete[]newobjvalue; 
}

/* 输出每一代演化的最佳结果及其适应值 */
void LSM::WriteResult()
{
	CString csToUserFileName = mGcsFileName+_T(".result");//file to user
	ofstream cResultFile;
	cResultFile.open((LPCSTR)csToUserFileName,ios::app);
	if (cResultFile.fail())
	{
		AfxMessageBox("Cannot open output result file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	cResultFile << "第" << m_numlap << "步优化结果:" << endl;
	cResultFile << "适应值:" << m_fitness << endl;
	cResultFile << "反演参数值依次为:";

	for(int i=0;i<m_numvar;i++)
		cResultFile << p_variable[i] <<"  ";
	cResultFile << endl;

	cResultFile << "各量测点量测值与计算值之差分别为:" ;
	for(i=0; i<m_numpoint; i++)
		cResultFile << p_objvalue[i] <<"  ";
	cResultFile << endl << endl;

	cResultFile.close();
}