simplex.cpp

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

//单纯形法类cpp文件
//Edit by Wang Shimin
//2004.8.25～2004.9.10
//本文件为最新文件，修改于2006.3.12

// Simplex.cpp: implementation of the Simplex class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Simplex.h"
#include "math.h"
#include "fstream"

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

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

Simplex::Simplex()
{
	/*	初始化参数值 */
	m_precisionnum = 0; /* 连续多少步不进化则终止 */
	m_maxnumlap = 10;   /* 最大循环次数 */
	m_numlap = 1;       /* 循环次数 */
	m_precision = 0.000001;  /* 收敛精度 */

//	m_sidelength = 1.0;  /* 单纯形边长 */
	m_oldfitness = 0.0; /* 上一步最优目标值 */
//	m_parametereq = (sqrt((double)m_numvar+1.0)+(double)m_numvar-1.0)/((double)m_numvar*sqrt(2.0))*m_sidelength;   /* 参数Q */
//	m_parameterep = m_parametereq-m_sidelength/sqrt(2.0);   /* 参数P */
	m_parameterep = 0.1; /* 参数P */
	m_parametereq = 0.2;/* 参数Q */
	m_alfa = 1.0;   /* 反射系数，大于0 */
	m_beta = 0.5;   /* 收缩系数，大于0，小于1 */
	m_gamma = 1.0;  /* 延伸系数,大于0 */

	/*将指针置零*/
	p_Bndload = NULL;
	p_BAKindflag = NULL;
	p_iparameter = NULL;
	p_jparameter = NULL;
	m_constriction = NULL;
	p_prolongation = NULL;
	p_center = NULL;
	p_echo = NULL;
	p_variable = NULL;
	p_bestobjvalue = NULL;
	p_fitness = NULL;
	p_observation = NULL;
	p_calculation = NULL;

	m_blErrorExit = false;
}

Simplex::~Simplex()
{
	/*释放内存空间*/
	for(int i=0; i<=m_numvar; i++)
	{
		delete[]p_variable[i];
	}
	delete[]p_variable;

	delete[]m_constriction;
	delete[]p_prolongation;
	delete[]p_center;
	delete[]p_echo;
	delete[]p_bestobjvalue;
	delete[]p_fitness;
	delete[]p_observation;
	delete[]p_calculation;
	delete[]p_BAKindflag;
	delete[]p_iparameter;
	delete[]p_jparameter;
 	delete[]p_Bndload;
}

/* 计算 */
void Simplex::Compute()
{
	/* 初始化函数 */
	Initialization();
	if(m_blErrorExit)return;
	
	while(Criterion())
	{
		m_oldfitness = m_lowvalue;
		/* 反射 */
		Echo();
		if(m_blErrorExit)return;

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

		m_numlap++;
	}
}

/* 初始化函数 */
void Simplex::Initialization()
{
	double ** pobjvalue;     /* 存储目标值，计算值与观测值之间的差值  */
	pobjvalue = NULL;

	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();

	int i,j;

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

	if(p_BAKindflag)
		delete[]p_BAKindflag;

	if(p_iparameter)
		delete[]p_iparameter;

	if(p_jparameter)
		delete[]p_jparameter;

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

	/* 存储量测值与计算值之差 */
	if(pobjvalue)
		delete[]pobjvalue;
	pobjvalue = new double *[m_numvar+1]; 
	for(i=0; i<=m_numvar; i++)
		pobjvalue[i]=new double[m_numpoint];
	/* 存储量测值与计算值之差 */

	/* 为单纯形分配动态内存空间 */
	if(p_variable)
		delete[]p_variable;
	p_variable = new double *[m_numvar+1];

	for(i=0; i<=m_numvar; i++)
		p_variable[i]=new double[m_numvar];

	/* *********** 建立初始单纯形 ************* */
	
	/* 从BAA文件中读入初始点 */
	ReadBAA();
	if(m_blErrorExit)return;

	/* 形成单纯形 */
	for(i=1; i<=m_numvar; i++)
	{
		for(j=0; j<m_numvar; j++)
		{
			if(j==(i-1))
				p_variable[i][j] = p_variable[0][j]*(1.0+m_parameterep);
			else
				p_variable[i][j] = p_variable[0][j]*(1.0+m_parametereq);
		}
	}

	if(m_constriction)
		delete[]m_constriction;
	m_constriction = new double[m_numvar];  /* 为收缩点变量分配动态内存空间 */

	if(p_prolongation)
		delete[]p_prolongation;
	p_prolongation = new double[m_numvar];  /* 为延伸点变量分配动态内存空间 */

	if(p_center)
		delete[]p_center;
	p_center = new double[m_numvar];  /* 为中心点变量分配动态内存空间 */

	if(p_echo)
		delete[]p_echo;
	p_echo = new double[m_numvar];   /* 为反射点变量分配动态内存空间 */ 

	if(p_bestobjvalue)
		delete[]p_bestobjvalue;
	p_bestobjvalue = new double[m_numpoint]; /* 存储量测值与计算值之差 */
	
	if(p_fitness)
		delete[]p_fitness;
	p_fitness = new double[m_numvar+1]; /* 存储单纯形中的目标值 */

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

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

	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++)
		{
			pobjvalue[i][j] = p_observation[j]-p_calculation[j];
			p_fitness[i] += pow(pobjvalue[i][j],2.0);
		}
	}
	
	/* 求最大值和最小值 */
	m_highvalue = p_fitness[0];
	m_lowvalue = p_fitness[0];
	m_highnumber = 0;
	m_lownumber = 0;
	
	for(i=1; i<=m_numvar; i++)
	{
		if(p_fitness[i]>m_highvalue)
		{
			m_highvalue = p_fitness[i];
			m_highnumber = i;
		}
		
		if(p_fitness[i]<m_lowvalue)
		{
			m_lowvalue = p_fitness[i];
			m_lownumber = i;
		}
	}
	
	m_oldfitness = m_lowvalue;
	
	for(i=0;i<m_numpoint;i++)
	{
		p_bestobjvalue[i] = pobjvalue[m_lownumber][i];
	}
	
	delete [] pobjvalue;
}

/* 判别函数 */
bool Simplex::Criterion()
{
	bool flag;
	flag = TRUE;

	if((m_numlap>m_maxnumlap)||(m_precisionnum>0))
		flag = FALSE;

	return flag;
}

/* 反射 */
void Simplex::Echo()
{
	int i,j;
	double sum;
	double * peobjvalue;
	peobjvalue = new double[m_numpoint];
	bool flag;
	flag=TRUE;

	/* 求最大值和最小值 */
	m_highvalue = p_fitness[0];
	m_lowvalue = p_fitness[0];
	m_highnumber = 0;
	m_lownumber = 0;

	for(i=1; i<=m_numvar; i++)
	{
		if(p_fitness[i]>m_highvalue)
		{
			m_highvalue = p_fitness[i];
			m_highnumber = i;
		}
		
		if(p_fitness[i]<m_lowvalue)
		{
			m_lowvalue = p_fitness[i];
			m_lownumber = i;
		}
	}
	
	for(i=0; i<m_numvar; i++)
	{
		sum = 0.0;
		
		for(j=0; j<=m_numvar; j++)
			sum += p_variable[j][i]-p_variable[m_highnumber][i];
	
		p_center[i] = fabs(sum/double(m_numvar));

		if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
			p_echo[i] = p_variable[m_lownumber][i]-m_alfa*p_center[i];
		else
			p_echo[i] = p_variable[m_lownumber][i]+m_alfa*p_center[i];
	}

	/****************** 求中心点的目标值 *******************/

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

	/* 计算测点的观测值与计算值之差 */
	for(i=0; i<m_numpoint; i++)
		m_centervalue += pow((p_observation[i]-p_calculation[i]),2.0);

	/****************** 求中心点的目标值 *******************/

	/****************** 求反射点的目标值 *******************/

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

	/* 计算测点的观测值与计算值之差 */
	for(i=0; i<m_numpoint; i++)
	{
		peobjvalue[i] = p_observation[i]-p_calculation[i];
		m_echovalue += pow(peobjvalue[i],2.0);
	}

	/****************** 求反射点的目标值 *******************/
	
	if(m_echovalue<=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_echo[i];
		}

		p_fitness[m_highnumber] = p_fitness[m_lownumber];
		m_highvalue = m_lowvalue;

		p_fitness[m_lownumber] = m_echovalue;
		m_lowvalue = m_echovalue;
		
		for(i=0;i<m_numpoint;i++)
		{
			p_bestobjvalue[i] = peobjvalue[i];
		}
		
		/* 延伸 */
		Prolongation();
		if(m_blErrorExit)return;
	}
	else
	{
		for(i=0; i<m_numvar; i++)
		{
			if((i!=m_highnumber)&&(p_fitness[i]>m_echovalue))
				flag=FALSE;
		}
		
		if(flag=TRUE)
		{
			if(m_echovalue>=m_highvalue)
			{
				/* 收缩 */
				Constriction();
				if(m_blErrorExit)return;
			}
			else
			{
				/* 用反射点替换单纯形中目标值最大的点 */
				for(i=0; i<m_numvar; i++)
					p_variable[m_highnumber][i] = p_echo[i];
				
				p_fitness[m_highnumber] = m_echovalue;
				m_highvalue = m_echovalue;
				
				/* 收缩 */
				Constriction();
				if(m_blErrorExit)return;
			}
		}
		else
		{
			/* 用反射点替换单纯形中目标值最大的点 */
			for(i=0; i<m_numvar; i++)
				p_variable[m_highnumber][i] = p_echo[i];
			
			p_fitness[m_highnumber] = m_echovalue;
			m_highvalue = m_echovalue;
		}
	}

	delete [] peobjvalue;
}

/* 延伸 */
void Simplex::Prolongation()
{
	int i,flag;

	flag = 0;

	double * ppobjvalue;
	ppobjvalue = new double[m_numpoint];
	
    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];