estimview.cpp

单纯形法是一种不错的随机搜索算法

// EstimView.cpp : implementation of the CEstimView class
//

#include "stdafx.h"
#include "Estim.h"

#include "EstimDoc.h"
#include "EstimView.h"
#include "math.h"
#include "string.h"

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

/////////////////////////////////////////////////////////////////////////////
// CEstimView

IMPLEMENT_DYNCREATE(CEstimView, CView)

BEGIN_MESSAGE_MAP(CEstimView, CView)
	//{{AFX_MSG_MAP(CEstimView)
	ON_COMMAND(ID_Simp, OnSimp)
	ON_COMMAND(ID_Simp_FireFall, OnSimpFireFall)
	ON_COMMAND(ID_Simp_Heredity, OnSimpHeredity)
	ON_COMMAND(ID_SimpCalc, OnSimpCalc)
	ON_COMMAND(ID_SimulateFire, OnSimulateFire)
	ON_COMMAND(ID_Inherit, OnInherit)
	ON_COMMAND(ID_VariSimp, OnVariSimp)
	ON_COMMAND(ID_ChangeL_Vari, OnChangeLVari)
	//}}AFX_MSG_MAP
	// Standard printing commands
	ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CEstimView construction/destruction

CEstimView::CEstimView()
{
	XResult[0] =10;
	XResult[1] =10;
}

CEstimView::~CEstimView()
{
}

BOOL CEstimView::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CEstimView drawing

void CEstimView::OnDraw(CDC* pDC)
{
	CEstimDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	// TODO: add draw code for native data here
}

/////////////////////////////////////////////////////////////////////////////
// CEstimView printing

BOOL CEstimView::OnPreparePrinting(CPrintInfo* pInfo)
{
	// default preparation
	return DoPreparePrinting(pInfo);
}

void CEstimView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add extra initialization before printing
}

void CEstimView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CEstimView diagnostics

#ifdef _DEBUG
void CEstimView::AssertValid() const
{
	CView::AssertValid();
}

void CEstimView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CEstimDoc* CEstimView::GetDocument() // non-debug version is inline
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CEstimDoc)));
	return (CEstimDoc*)m_pDocument;
}
#endif //_DEBUG

///////////////////////////////////////////////////////////////////////////////

//单纯形法
void CEstimView::OnSimpCalc() 
{
	int l=5;
	iCycleNum=0;
	if(!SimpCalc(XResult, l))
		MessageBox("1计算单纯形过程中出错!");
	else
	{
		CString string;
		string.Format("单纯形法共循环了%d次，结果为：%lf,%lf,R(X)=%f,Fx=%lf",iCycleNum,XResult[0],XResult[1],RResult,FxMin);
		MessageBox(string);
	}

}
//改进的单纯形法
void CEstimView::OnSimp() 
{

	int l=5,NoEvolve=1;
	int iCalcNum=0;     //进行单纯形计算的次数
	bool bEnd= false;
	iCycleNum=0;
	if(!SimpCalc(XResult, l))
		MessageBox("2计算单纯形过程中出错!");

	while(!bEnd)
	{
		//改变步长，重新计算单纯形
		l=l+1;
		//顶点数据临时存储在Temp中
		long double XTemp[2],RTemp,dR;
		XTemp[0] = XResult[0];
		XTemp[1] = XResult[1];
		RTemp = RResult;
	
        if(!SimpCalc(XResult, l))
			MessageBox("3计算单纯形过程中出错!");
		else
			iCalcNum++;
		dR = RResult - RTemp;
		
		//Rx 连续变大20次了，结束循环，取Rx小的为结果
		if(dR>=0.0)
		{
			NoEvolve++;
			if(NoEvolve>20)   //连续20次都没有改进
			{
				XResult[0] = XTemp[0];
				XResult[1] = XTemp[1];
				RResult = RTemp;
				bEnd = true;
			}
			
		}
		//Rx变小了，结果再作顶点，重新循环，循环次数最多1000
		if(iCalcNum>1000)
			bEnd = true;
	}
	CString string;
	string.Format("改进的单纯形法共循环了%d次，结果为：步长：%d\n%lf,%lf,R(X)=%f,Fx=%lf",iCycleNum,l,XResult[0],XResult[1],RResult,FxMin);
	MessageBox(string);
}
//均匀变异的单纯形法
//优点：全局
//缺点：参数多的时候，复杂
void CEstimView::OnVariSimp() 
{
	int l=5;                           //步长
	int j;
	double dMin=-10.0,dMax=10.0;       //设定变异的范围
	long double XMin[2],RMin,dR;
	RMin=999999999.0;
	int iCalcNum=int(20/l);   //分区运用单纯形法
	for(iCycleNum=0;iCycleNum<iCalcNum;iCycleNum++)
	{
		XResult[0] =dMin+iCycleNum*(dMax-dMin)/iCalcNum;
		for(j=0;j<iCalcNum;j++)
		{
			XResult[1] =dMin+j*(dMax-dMin)/iCalcNum;
			if(!SimpCalc(XResult, l))
			{	MessageBox("2计算单纯形过程中出错!");
			}
			dR = RResult - RMin;
			
			//取Rx小的为结果
			if(dR<=0.0)
			{
				XMin[0] = XResult[0];
				XMin[1] = XResult[1];
				RMin = RResult;
			}
		}
	}	
	CString str;
	str.Format("均匀变异的单纯形法结果为：步长：%d\n%lf,%lf,R(X)=%f",l,XMin[0],XMin[1],RMin);
	MessageBox(str);
	
}
//变步长，随机变异
void CEstimView::OnChangeLVari() 
{	
	int l,i;                           //步长
	double dMin=-10.0,dMax=10.0;       //设定变异的范围
	long double XMin[2],RMin,dR;
	RMin=999999999.0;
	bool bEnd= false;
	int iChangeLNum=10,iVariNum=10;   //变步长、变异次数
	int iNoImproveNum=0;    //连续10次没有改进，结束循环
	int iCycleNum=0;
	while(!bEnd&&iCycleNum<1000)
	{
		l=5;
		iCycleNum++;
		if(!SimpCalc(XResult, l))
		MessageBox("2计算单纯形过程中出错!");
		dR = RResult - RMin;			
		//取Rx小的为结果
		if(dR<=0.0)
		{
			XMin[0] = XResult[0];
			XMin[1] = XResult[1];
			RMin = RResult;
			iNoImproveNum=0;
		}
		else 
			iNoImproveNum++;
		//变步长
		for(i=0;i<iChangeLNum;i++)   //10次
		{
			l=l+1;
			//顶点数据临时存储在Temp中
			if(!SimpCalc(XResult, l))
				MessageBox("3计算单纯形过程中出错!");
			dR = RResult - RMin;			
			//取Rx小的为结果
			if(dR<=0.0)
			{
				XMin[0] = XResult[0];
				XMin[1] = XResult[1];
				RMin = RResult;
				iNoImproveNum=0;
			}
			else 
				iNoImproveNum++;
		}
		//随机变异
		l=5;
		for(i=0;i<iVariNum;i++)    //次数
		{
			XResult[0]=dMin+rand()%1000/1000.0*20;
			XResult[1]=dMin+rand()%1000/1000.0*20;
			if(!SimpCalc(XResult, l))
			MessageBox("2计算单纯形过程中出错!");
			dR = RResult - RMin;			
			//取Rx小的为结果
			if(dR<=0.0)
			{
				XMin[0] = XResult[0];
				XMin[1] = XResult[1];
				RMin = RResult;
				iNoImproveNum=0;
			}
			else 
				iNoImproveNum++;
		}
		if (iNoImproveNum>20)
			bEnd = true;
	}
	
	CString str;
	str.Format("均匀变异的单纯形法结果为：步长：%d\n%lf,%lf,R(X)=%f",l,XMin[0],XMin[1],RMin);
	MessageBox(str);
}
//单纯形－模拟退火法
void CEstimView::OnSimpFireFall() 
{	
	int l=20;
	iCycleNum=0;
	SimpCalc(XResult,l);

	long double XTemp[2],RTemp;
	bool bEnd= false;
	
	while(!bEnd)
	{
		//顶点数据临时存储在Temp中
		XTemp[0] = XResult[0];
		XTemp[1] = XResult[1];
		RTemp = RResult;

		FireFallCalc();
		if( (RResult - RTemp)>0.000000001 ||(RResult - RTemp)<-0.00000001)
		{
			SimpCalc(XResult,l);
			if(iCycleNum>10000)
				bEnd = true;
		}
		else 
			bEnd = true;
	}
	CString string;
	string.Format("单纯形——模拟退火共循环了%d次，结果为：%lf,%lf,R(X)=%f,Fx=%lf",iCycleNum,XResult[0],XResult[1],RResult,FxMin);
	MessageBox(string);
}
//单纯形－遗传算法
void CEstimView::OnSimpHeredity() 
{
}
//模拟退火法
void CEstimView::OnSimulateFire() 
{
	XResult[0] =-0.8;
	XResult[1] =1.4;
	iCycleNum=0;
	RResult = R(XResult);
	
	FireFallCalc();
		
	CString string;
		string.Format("模拟退火法共循环了%d次，结果为：%lf,%lf,R(X)=%f,Fx=%lf",iCycleNum,XResult[0],XResult[1],RResult,FxMin);
		MessageBox(string);
}
//遗传算法
void CEstimView::OnInherit() 
{	
}

//求平方值
//参数：long double 
//返回：long double
long double CEstimView::Square(long double x)
{
	return x*x;
}
//单纯形法终止条件判断
//参数：三个顶点的目标函数值
//返回：满足终止条件（<0.007）时为true
bool CEstimView::SimpEndOrNot(long double r[3], int imin)
{
	
	int i;
	long double delta = 0.0000000001,error=0.0;
	for(i=0;i<3;i++)
	{
		error = error+Square((r[i]-r[imin]));
	}

	error = sqrt(error/3.0);
	if (error<= delta)
		return true;
	else 
		return false;
}
//求目标函数的大小
//参数：顶点数组（x1,x2）
//返回：目标函数大小Rx
long double CEstimView::R(long double x[2])
{
	long double l[5] ={-186.7340,-186.736,-186.735,-186.742,-186.739};
	long double Fx,Rx=0.0,temp1=0.0,temp2=0.0;
	for(int i=1;i<6;i++)
	{
		temp1=temp1+i*cos( (i+1)*x[0]+i );
		temp2=temp2+i*cos( (i+1)*x[1]+i );
	}
    Fx=temp1*temp2+0.5*( Square(x[0]+1.42513) + Square(x[1]+0.80032) );
	for(i=0;i<5;i++)
	{
		Rx=Rx+Square(Fx-l[i]);
	}
	FxMin =Fx;
	return Rx;
}
//运用单纯形法求新顶点
//参数：顶点数组（x1，x2）、步长l
//返回：新的顶点值X0（x1，x2）赋值给全局变量XResult[2]
//      计算成功，返回 true
bool CEstimView::SimpCalc(long double x[2], int l)
{
	//生成一个单纯形结果查看文件，便于查错
	/*
	CStdioFile FileTxt;
	CString strWrite,TextFileName = "D:\\单纯形法结果.txt";
	FileTxt.Open(TextFileName, CFile::modeCreate | CFile::modeWrite |CFile::typeText);
	*/
    long double X[3][2],c,d;
	X[0][0] = x[0];
	X[0][1] = x[1];
	//t为参数个数，t+1为顶点数
	int t=2;
	//求另外两个顶点值
	c=(sqrt(t+1.0) + t-1.0)*l/sqrt(2.0)/t;
	d=(sqrt(t+1.0) -1.0)*l/sqrt(2.0)/t;
	X[1][0] = X[0][0] +c;
	X[1][1] = X[0][1] +d;
	X[2][0] = X[0][0] +d;
	X[2][1] = X[0][1] +c;
	
	//Rx[3]为三个顶点的目标函数值,Xe,Re为新顶点的顶点值、目标函数值
	//Xc,Rc,Xr,Rr为中心点和反射点的
	//iMax,iMedium,iMin为三个顶点目标函数值的排序
	long double Rx[3],Xe[2],Re,Xc[2],Xr[2],Rr;
	int iMax,iMedium,iMin;
	//求目标函数值
	for(int i =0;i<3;i++)
	{
		Rx[i]= R(X[i]);
	}	
	//给目标函数值排序
	long double Ru=-1.7e308,Rl=1.7e308;
	for(i=0;i<3;i++)
	{
		if(Ru<Rx[i])
		{
			Ru = Rx[i];
			iMax = i;
		}
		if(Rl>Rx[i])
		{
			Rl = Rx[i];
			iMin = i;
		}
	}
	iMedium = 3-iMin-iMax;
	while(!SimpEndOrNot(Rx,iMin) )
	{
		for(i=0;i<2;i++)
		{
			//计算中心点Xc[2],Rc
			Xc[i] = ( X[iMin][i]+X[iMedium][i] )/2.0;
			//计算反射点Xr[2],Rr
			Xr[i]= X[iMin][i]+X[iMedium][i]-X[iMax][i];
		}
		Rr=R(Xr);
		//(1)伸长
		if(Rr<=Rx[iMin])
		{
			for(i=0;i<2;i++)
			{
				Xe[i] = 3.0*Xr[i]-2.0*Xc[i];
			}
			Re=R(Xe);
			if(Re<Rx[iMin])
			{
				X[iMax][0] = Xe[0];
				X[iMax][1] = Xe[1];
			}
			else
			{
				X[iMax][0] = Xr[0];
				X[iMax][1] = Xr[1];
			}
		}
		//(2)缩短
		else if(Rr>=Rx[iMax])
		{
			for(i=0;i<2;i++)
			{
				Xe[i] = 0.5*(Xc[i] +Xr[i]);
			}
			Re=R(Xe);
			if(Re<Rx[iMin])
			{
				X[iMax][0] = Xe[0];
				X[iMax][1] = Xe[1];
			}
			//缩小
			else
			{
				for(i=0;i<2;i++)
				{
					X[iMax][i] = (X[iMax][i]+X[iMin][i])/2.0;
					X[iMedium][i] = (X[iMedium][i]+X[iMin][i])/2.0;
				}
			}
		}
		else
		{
			X[iMax][0] = Xr[0];
			X[iMax][1] = Xr[1];
		}
		//重求目标函数值
		for(i=0;i<3;i++)
		{
			Rx[i] = R(X[i]);
		}

		//给目标函数值重新排序
		Ru=-1.7e308,Rl=1.7e308;
		for(i=0;i<3;i++)
		{
			if(Ru<Rx[i])
			{
				Ru = Rx[i];
				iMax = i;
			}
			if(Rl>Rx[i])
			{
				Rl = Rx[i];
				iMin = i;
			}
		}
		iMedium = 3-iMin-iMax;
		iCycleNum++;

		//写单纯形法结果文件
		/*
		strWrite.Format("%d,%lf,%lf,%lf\n",iCycleNum,X[iMin][0],X[iMin][1],Rx[iMin]);
		FileTxt.WriteString(strWrite);
		*/
	}
	XResult[0] = X[iMin][0];
	XResult[1] = X[iMin][1];
	RResult = Rx[iMin];

	//FileTxt.Close();
	return true;
}
//模拟退火法计算函数
//参数：
//返回：新的顶点数组
//      计算成功为 True
bool CEstimView::FireFallCalc()
{
	//
	long double T,XNew[2],RNew,DeltaR,Random,m,r=5.0;
	
	//k为循环次数
	int k=1;
	T = 100.0;
	//温度T虽然不断下降，但是不会为0，所以while(T)....必死循环
	while(T && iCycleNum<10000)
	{
		//求(－1，1)之间的随机数
		//m是为了求Random的符号
		m = RandomNumber(&r);
		if(m>0.5)
			m =1.0;
		else
			m = -1.0;
		Random = rand()%1000/1000.0;
		Random = Random*m;
		//X1=X0+Random*V0,V0为步长,取5.0,Random为(-1,1)间的随机数
		XNew[0] = XResult[0] + Random*4.0;
		XNew[1] = XResult[1] + Random*4.0;
		RNew = R(XNew);

		DeltaR = RNew - RResult;
		if(DeltaR<0)
		{
			//(1)王老师选择的降温方式：
			 T = T/( log10(10.0+sqrt(k)) );

			//(2)快速降温方式：见《快速模拟退火算法及应用》
			//T(k) = T0 * exp(-C*pow(k,1/N)
			//式中：C为给定常数，N为待反演参数的个数。可改写为：
			//T(k) = T0 *pow( a,pow(k,1/N) )
			//通常a选择[0.7,1]，1/N选择0.5或者1.0

		//	T = T*pow(0.8,pow(k,0.5) );
			XResult[0] = XNew[0];
			XResult[1] = XNew[1];
			RResult = RNew;
			k++;
		}
		else
		{
			double A,P;
			A = rand()%1000/1000.0;
			P = exp(-DeltaR/T);
			if(P>A)
			{
				XResult[0] = XNew[0];
				XResult[1] = XNew[1];
				RResult = RNew;
			}
		}
		++iCycleNum;
	}
	return true;
}
//产生一个随机数(0，1)
//参数：一个指向long double的指针变量
//返回：long double 型(0，1)间的随机变量
long double CEstimView::RandomNumber(long double *r)

{
	int m;
    long double s,u,v,p;
    s=65536.0; u=2053.0; v=13849.0;
    m=(int)(*r/s); 
	*r=*r-m*s;
    *r=u*(*r)+v; m=(int)(*r/s);
    *r=*r-m*s; 
	p=*r/s;
    return(p);
}