simplex.cpp

来自「改进的单纯形算法」· C++ 代码 · 共 973 行 · 第 1/2 页

CPP
973
字号
		m_prolongationvalue += pow(ppobjvalue[i],2.0);
	}
	
	if(m_prolongationvalue<m_lowvalue)
	{
	   	/* 继续延伸 */
		while((m_prolongationvalue<=m_lowvalue)&&(flag<4))
		{
	    	/* 去掉原来的最大值点,取而代之以原来的最小点,
		    	并用延伸点替换单纯形中目标值最小的点 */
			for(i=0; i<m_numvar; i++)
			{
				p_variable[m_highnumber][i] = p_variable[m_lownumber][i];
				p_variable[m_lownumber][i] = p_prolongation[i];
			}
			
			p_fitness[m_highnumber] = p_fitness[m_lownumber];
			m_highvalue = m_lowvalue;
			
			p_fitness[m_lownumber] = m_prolongationvalue;
			m_lowvalue = m_prolongationvalue;

			for(i=0;i<m_numpoint;i++)
			{
				p_bestobjvalue[i] = ppobjvalue[i];
			}

			for(i=0; i<m_numvar; i++)
			{
				if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
					p_prolongation[i] = p_variable[m_lownumber][i]-m_gamma*p_center[i];
				else
					p_prolongation[i] = p_variable[m_lownumber][i]+m_gamma*p_center[i];
			}
			
			/* 将参数存储到文件 */
			WriteToFem(p_prolongation);
			if(m_blErrorExit)return;
			/* 有限元计算 */
			ObjFunction();
			if(m_blErrorExit) return;
			/* 读入测点的计算值 */
			ReadFromFem(p_calculation);
			if(m_blErrorExit)return;
			
			m_prolongationvalue = 0.0;
			
			/* 计算测点的观测值与计算值之差 */
			for(i=0; i<m_numpoint; i++)
			{
				ppobjvalue[i] = p_observation[i]-p_calculation[i];
				m_prolongationvalue += pow(ppobjvalue[i],2.0);
			}
		}
		
		if(m_prolongationvalue<m_lowvalue)
		{
		    /* 去掉原来的最大值点,取而代之以原来的最小点,
			    并用延伸点替换单纯形中目标值最小的点 */
			for(i=0; i<m_numvar; i++)
			{
				p_variable[m_highnumber][i] = p_variable[m_lownumber][i];
				p_variable[m_lownumber][i] = p_prolongation[i];
			}
			
			p_fitness[m_highnumber] = p_fitness[m_lownumber];
			m_highvalue = m_lowvalue;
			
			p_fitness[m_lownumber] = m_prolongationvalue;
			m_lowvalue = m_prolongationvalue;

			for(i=0;i<m_numpoint;i++)
			{
				p_bestobjvalue[i] = ppobjvalue[i];
			}
		}
		else
		{
			/* 用延伸点替换单纯形中目标值最大的点 */
			for(i=0; i<m_numvar; i++)
			{
				p_variable[m_highnumber][i] = p_prolongation[i];
			}
			
			p_fitness[m_highnumber] = m_prolongationvalue;
			m_highvalue = m_prolongationvalue;
			
			/* 收缩 */
			Constriction();
			if(m_blErrorExit)return;
		}
	}
	else
	{
		if(m_prolongationvalue<m_highvalue)
		{
		    /* 用延伸点替换单纯形中目标值最大的点 */
			for(i=0; i<m_numvar; i++)
			{
				p_variable[m_highnumber][i] = p_prolongation[i];
			}
			
			p_fitness[m_highnumber] = m_prolongationvalue;
			m_highvalue = m_prolongationvalue;

			/* 收缩 */
			Constriction();
			if(m_blErrorExit)return;
		}
		else
		{
			/* 延伸失败,调整延伸系数,重新延伸 */
			flag = 0;
			while((m_prolongationvalue>m_highvalue)&&(flag<4))
			{
				/* 调整延伸系数 */
				m_gamma /= 2.0;
								
				for(i=0; i<m_numvar; i++)
				{
					if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
						p_prolongation[i] = p_variable[m_lownumber][i]-m_gamma*p_center[i];
					else
						p_prolongation[i] = p_variable[m_lownumber][i]+m_gamma*p_center[i];
				}
							
				/* 将参数存储到文件 */
				WriteToFem(p_prolongation);
				if(m_blErrorExit)return;
				/* 有限元计算 */
				ObjFunction();
				if(m_blErrorExit) return;
				/* 读入测点的计算值 */
				ReadFromFem(p_calculation);
				if(m_blErrorExit)return;
				
				m_prolongationvalue = 0.0;
				
				/* 计算测点的观测值与计算值之差 */
				for(i=0; i<m_numpoint; i++)
				{
					ppobjvalue[i] = p_observation[i]-p_calculation[i];
					m_prolongationvalue += pow(ppobjvalue[i],2.0);
				}
			}

			m_gamma = 1.0;
				
			if(m_prolongationvalue<m_lowvalue)
			{
			    /* 去掉原来的最大值点,取而代之以原来的最小点,
				   并用延伸点替换单纯形中目标值最小的点 */
				for(i=0; i<m_numvar; i++)
				{
					p_variable[m_highnumber][i] = p_variable[m_lownumber][i];
					p_variable[m_lownumber][i] = p_prolongation[i];
				}
				
				p_fitness[m_highnumber] = p_fitness[m_lownumber];
				m_highvalue = m_lowvalue;
				
				p_fitness[m_lownumber] = m_prolongationvalue;
				m_lowvalue = m_prolongationvalue;

				for(i=0;i<m_numpoint;i++)
				{
					p_bestobjvalue[i] = ppobjvalue[i];
				}
			}
			else
			{
				/* 用延伸点替换单纯形中目标值最大的点 */
				for(i=0; i<m_numvar; i++)
				{
					p_variable[m_highnumber][i] = p_prolongation[i];
				}
				
				p_fitness[m_highnumber] = m_prolongationvalue;
				m_highvalue = m_prolongationvalue;
				
				/* 收缩 */
				Constriction();
				if(m_blErrorExit)return;
			}
		}
	}

	delete [] ppobjvalue;
}

/* 收缩 */
void Simplex::Constriction()
{
	int i,j;
	double * pcobjvalue;
	pcobjvalue = new double[m_numpoint];
	
	/****************** 求收缩点及其目标值 ******************/

	for(i=0; i<m_numvar; i++)
	{
		if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
			m_constriction[i] = p_variable[m_highnumber][i]-m_beta*p_center[i];
		else
			m_constriction[i] = p_variable[m_highnumber][i]+m_beta*p_center[i];
	}
	
	/* 将参数存储到文件 */
	WriteToFem(m_constriction);
	if(m_blErrorExit)return;
	/* 有限元计算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 读入测点的计算值 */
	ReadFromFem(p_calculation);
	if(m_blErrorExit)return;
	
	m_constrictionvalue = 0.0;
	/* 计算测点的观测值与计算值之差 */
	for(i=0; i<m_numpoint; i++)
	{
		pcobjvalue[i] = p_observation[i]-p_calculation[i];
		m_constrictionvalue += pow(pcobjvalue[i],2.0);
	}
	/****************** 求收缩点及其目标值 ******************/
	
	if(m_constrictionvalue<m_lowvalue)
	{
	    /* 去掉原来的最大值点,取而代之以原来的最小点,
	    	并用收缩点替换单纯形中目标值最小的点 */
		for(i=0; i<m_numvar; i++)
		{
			p_variable[m_highnumber][i] = p_variable[m_lownumber][i];
			p_variable[m_lownumber][i] = m_constriction[i];
		}
		
		p_fitness[m_highnumber] = p_fitness[m_lownumber];
		m_highvalue = m_lowvalue;
		
		p_fitness[m_lownumber] = m_constrictionvalue;
		m_lowvalue = m_constrictionvalue;

		for(i=0;i<m_numpoint;i++)
		{
			p_bestobjvalue[i] = pcobjvalue[i];
		}
	}
	else
	{
		if(m_constrictionvalue<m_highvalue)
		{
			/* 用收缩点替换单纯形中目标值最大的点 */
			for(i=0; i<m_numvar; i++)
			{
				p_variable[m_highnumber][i] = p_prolongation[i];
			}
			
			p_fitness[m_highnumber] = m_prolongationvalue;
			m_highvalue = m_prolongationvalue;
			
			/* 缩小边长 */
			for(i=0; i<=m_numvar; i++)
				for(j=0; j<m_numvar; j++)
					p_variable[i][j] = 0.5*(p_variable[m_lownumber][j]+p_variable[i][j]);
				
			
			/******** 重新计算新的单纯形的目标值 ********/
			for(i=0; i<=m_numvar; i++)
			{
				p_fitness[i] = 0.0;
				
				/* 将参数存储到文件 */
				WriteToFem(p_variable[i]);
				if(m_blErrorExit)return;
				/* 有限元计算 */
				ObjFunction();
				if(m_blErrorExit) return;
				/* 读入测点的计算值 */
				ReadFromFem(p_calculation);
				if(m_blErrorExit)return;
				
				/* 计算测点的观测值与计算值之差的平方 */
				for(j=0; j<m_numpoint; j++)
				{
					pcobjvalue[j] = p_observation[j]-p_calculation[j];
					p_fitness[i] += pow(pcobjvalue[j],2.0);
				}
			}
		}
	}
}

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

/* 写参数值 */
void Simplex::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 Simplex::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();
}

/* 从BAA文件中读入待反演参数的性质标志以及初始值 */
void Simplex::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[0][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 Simplex::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 Simplex::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_lowvalue << endl;
	cResultFile << "反演参数值依次为:";

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

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

	cResultFile.close();
}

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?