外点法.cpp

用C++编的一些最优化作业中的程序,有Newton法,DFP法,共轭梯度法,单纯形法,内点法,外点法,内外点法,都能使用,我已经全部调试过了

#include<iostream.h>
#include<math.h>
int const n=2;
static double Mk=1;
double fun_original(double x[n]);//既是主函数，也是F(X,Mk)之一
double g_fun(double x[n]);//约束函数
double h_fun(double x[n]);;//等式函数
double a_e_f(double x[n]);
double fun(double x[n]);//让单纯形法求解此函数最小值
void compare(double x[n+1][n],double LGH[3][n],double LGH_tap[3]);//在调用之前LGH先给x[n+1]赋值
void zx(double x[n+1][n],double LGH[3][n],double LGH_tap[3],double xn_1[n],double xn_2[n]);
double H(double x[n+1][n],double LGH[3][n],double c);//判别准则
void danchun(double original_point[n]);
void main()
{
	int k=1;
	int i;
	double r=10;
	double c=0.000001;
	double x[n]={0,0.5};//给定的初始点在某个增广目标函数的求解域之中
	do
	{
		danchun(x);
		//ln(x2)>=0时F(X,Mk)=-X1+X2;无极小点，不用管这个增广目标函数
		if(Mk*a_e_f(x)<c)
			break;
		Mk=r*Mk;
	}while(1);
	cout<<"用外点法求解："<<endl;
	cout<<"minf(X)=-X1+X2"<<endl;
	cout<<"     g1(X)=lnX2>=0"<<endl<<"s.t."<<endl;;
	cout<<"     g2(X)=X1+X2-1=0"<<endl;
	cout<<"结果是："<<endl;
	cout<<"Xk="<<endl;
	for(i=0;i<n;i++)
		cout<<x[i]<<" ";
	cout<<endl;
	cout<<"f(Xk)="<<fun_original(x)<<endl;
	//cout<<"Mk="<<Mk<<endl;
}
double fun_original(double x[n])//既是主函数，也是F(X,Mk)之一
{
	double y;
	y=-x[0]+x[1];
	return y;
}
double g_fun(double x[n])//约束函数
{
	double y;
	y=pow(log(x[1]),2);
	return y;
}
double h_fun(double x[n])//等式函数
{
	double y;
	y=x[0]+x[1]-1;
	return y;
}
double a_e_f(double x[n])
{
	double y;
	y=Mk*(pow(h_fun(x),2)+pow(g_fun(x),2));
	return y;
}					
double fun(double x[n])//让单纯形法求解此函数最小值
{
	double y;
	y=fun_original(x)+a_e_f(x);
	return y;
}
void compare(double x[n+1][n],double LGH[3][n],double LGH_tap[3])//在调用之前LGH先给x[n+1]赋值
{
	int i;
	double y[n+1];
	for(i=0;i<n+1;i++)
		y[i]=fun(x[i]);//y[i]=fun(x[i][n]);//////////////////////////////////////////////////?关于传值疑问

	double temp[2];
	temp[0]=y[0];
	for(i=0;i<n+1;i++)//temp[0]存储最小的值XL,temp[1]存储对应x的下标
	{
		if(y[i]<=temp[0])
		{
			temp[0]=y[i];
			temp[1]=i;
		}
	}
	LGH_tap[0]=temp[1];//记录标记
	for(i=0;i<n;i++)
		LGH[0][i]=x[int(temp[1])][i];
	
	for(i=0;i<n+1;i++)//temp[0]存储最大的值XH,temp[1]存储对应x的下标
	{
		if(y[i]>=temp[0])
		{
			temp[0]=y[i];
			temp[1]=i;
		}
	}
	LGH_tap[2]=temp[1];//
	for(i=0;i<n;i++)
		LGH[2][i]=x[int(temp[1])][i];

	temp[0]=y[0];
	if(LGH_tap[2]==0)
	{
		temp[0]=y[1];
	}
	for(i=0;i<n+1;i++)//寻找次大值
	{		
		if(i==LGH_tap[2])//if(i==temp[1])
			continue;
		if(y[i]>=temp[0])
		{
			temp[0]=y[i];
			temp[1]=i;
		}
	}
	LGH_tap[1]=temp[1];
	for(i=0;i<n;i++)
		LGH[1][i]=x[int(temp[1])][i];
}
void zx(double x[n+1][n],double LGH[3][n],double LGH_tap[3],double xn_1[n],double xn_2[n])
{
	int i,j;
	for(i=0;i<n;i++)
		xn_1[i]=xn_2[i]=0;

	for(i=0;i<n+1;i++)//返回重心
	{
		if(i==LGH_tap[2])//轮流加n次（去掉一次）
			continue;
		for(j=0;j<n;j++)
		{		
			xn_1[j]+=x[i][j];
		}	
	}
	for(i=0;i<n;i++)
		xn_1[i]/=n;
	for(i=0;i<n;i++)//返回反射点
		xn_2[i]=2*xn_1[i]-LGH[2][i];
}
double H(double x[n+1][n],double LGH[3][n],double c)//判别准则
{
	int i;
	double s=0;
	for(i=0;i<n+1;i++)//小心i从0到n即(i=0;i<n+1;i++)
	{
		s+=pow(fun(x[i])-fun(LGH[0]),2);//s+=pow(fun(x[i][n])-fun(LGH[0][n]),2);//////////////////////////////////////////////////?关于传值疑问
	}
	if(s<=c)
		return 1;
	else 
		return 0;
}
void danchun(double original_point[n])
{
	double x[n+1][n]={0};
	double c=0.000001;
	double a=2;
	double b=0.5;
	double LGH[3][n]={0};
	double LGH_tap[3]={0};
	double xn_1[n]={0};
	double xn_2[n]={0};
	double xn_3[n]={0};
	double xn_4[n]={0};
	int i,j,tap=0;
	for(i=0;i<n;i++)
	x[0][i]=original_point[i];
	for(i=1;i<n+1;i++)//假设只给出一个初始点X0(1,1),求另外两个初始点
	{
		for(j=0;j<n;j++)
		{
			if(j+1==i)
				x[i][j]=x[0][j]+1;
			else
				x[i][j]=x[0][j];
		}
	}
	compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它们在x[n+1][n]中的下标
	zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射点X(n+2)
	while(H(x,LGH,c)==0)
	{
		if(fun(xn_2)<fun(LGH[0]))//if(fun(xn_2)<fun(LGH[0][n]))//////////////////////////////////////////////////?关于传值疑问
		{
			for(i=0;i<n;i++)
			{
				xn_3[i]=xn_1[i]+a*(xn_2[i]-xn_1[i]);//扩张
			}
			if(fun(xn_3)<fun(LGH[0]))
			{
				for(i=0;i<n;i++)
				{
					x[int(LGH_tap[2])][i]=LGH[2][i]=xn_3[i];
				}
				//goto step 2;
			}
			else
			{
				for(i=0;i<n;i++)
				{
					x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];
				}
				//goto step 2;
			}
		}
		else 
		{
			if(fun(xn_2)<fun(LGH[1]))
			{
				for(i=0;i<n;i++)
				{
					x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];
				}
				//goto step 2;
			}
			else 
			{
				if(fun(xn_2)<fun(LGH[2]))//else if(fun(xn_2)<fun(LGH[2][n]))
				{
					for(i=0;i<n;i++)
					{
						x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];
					}
				}
				else//means (fun(xn_2)>=fun(LGH[2][n]))
				{
					for(i=0;i<n;i++)
					{
						xn_4[i]=xn_1[i]+b*(xn_2[i]-xn_1[i]);//压缩
					}
					//但另出来出来
					if(fun(xn_4)<fun(LGH[2]))
					{
						for(i=0;i<n;i++)
						{				
							x[int(LGH_tap[2])][i]=LGH[2][i]=xn_4[i];
						}
						//goto step 2;
					}
					else
					{
						for(i=0;i<n+1;i++)
						{
							if(i==LGH_tap[0])
								continue;
							for(j=0;j<n;j++)
							{
								x[i][j]=(LGH[0][j]+x[i][j])/2;
							}
						}
					}//goto step 2
				}			
			}
		}

		//That's it!!!
		compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它们在x[n+1][n]中的下标
     	zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射点X(n+2)
		tap++;				
	}
	//cout<<"迭代次数:"<<tap<<endl;
	for(i=0;i<n;i++)
		original_point[i]=LGH[0][i];//把在某个Mk值下增广目标函数的最优值返回
}
//单纯小形法结束