内点法.cpp

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

#include<iostream.h>
#include<math.h>
#include<iomanip.h>
const int n=2;
static double Rk=10;//Rk是惩罚因子
static int p=1;
double fun_original(double x[n]);//既是主函数，也是F(X,Rk)之一
double a_e_f(double x[n]);
double fun(double x[n]);//让单纯形法求解此函数最小值
void g_fun(double x[n],double g0[n]);//牛顿法子函数开始,求出与Rk对应的梯度向量.
void Hesse_fun(double x[n],double Hesse[n][n]);//求出与Rk对应的Hesse矩阵
int H(double g0[n],double c);
void inv(double a[n][n],double c[n][n]);//求Hesse矩阵的逆矩阵
void xiang_cheng(double a[n][n],double b[n],double c[n]);//n*n矩阵乘以n*1矩阵
void newton(double x[n]);//牛顿法
void main()
{
	int k=1;
	int i;
	double c_r=0.5;//c是惩罚因子的缩小系数
	double cc=0.000001;//cc是内点法求最终最优解的终止限
	double x[n]={2,2};//给定的初始点在某个增广目标函数的求解域之中
	do
	{
		newton(x);
		if(fabs(a_e_f(x))<cc)//小心它的值为负，要加上绝对值
			break;
		Rk=c_r*Rk;
	}while(1);
	cout<<"用内点法求解："<<endl;
	cout<<"minf(X)=(X1+1)^3/3+X2"<<endl;
	cout<<"     g1(X)=X1-1>=0"<<endl<<"s.t."<<endl;;
	cout<<"     g2(X)=X2>=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<<"Rk="<<Rk<<endl;
}
double fun_original(double x[n])//既是主函数，也是F(X,Rk)之一
{
	double y;
	y=pow(x[0]+1,3)/3+x[1];
	return y;
}
double a_e_f(double x[n])
{
	double y;
	y=Rk*(1/(x[0]-1)+1/x[1]);
	return y;
}
double fun(double x[n])//让单纯形法求解此函数最小值
{
	double y;
	y=fun_original(x)+a_e_f(x);
	return y;
}
void g_fun(double x[n],double g0[n])//牛顿法子函数开始,求出与Rk对应的梯度向量.
{
	g0[0]=pow(x[0]+1,2)-Rk/pow(x[1]-1,2);
	g0[1]=1-Rk/pow(x[1],2);
}
void Hesse_fun(double x[n],double Hesse[n][n])//求出与Rk对应的Hesse矩阵
{
	Hesse[0][0]=2*(x[0]+1)+2*Rk/pow(x[1]-1,3);
	Hesse[0][1]=Hesse[1][0]=0;
	Hesse[1][1]=2*Rk/pow(x[1],3);
}
int H(double g0[n],double c)
{
	double s=0;
	for(int i=0;i<n;i++)
	{
		s+=pow(g0[i],2);
	}
	if(sqrt(s)<c)
		return 1;
	else 
		return 0;
}
void inv(double a[n][n],double c[n][n])//求Hesse矩阵的逆矩阵
{
	double m[n];//辅助乘数
	double T;//存储换行临时变量
	double temp=0;
	double t;//最大列主元
	int tap1=0,tap2=0;//最大列主元下标
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
			{
				c[i][j]=1;
			}
			else
			{
				c[i][j]=0;
			}
		}
	}
	for(int s=0;s<n;s++)
	{
		t=a[s][s];//赋初值
		for(int p=s;p<n;p++)//选取最大列主元
		{
			if(fabs(a[p][s])>t)
			{
				t=a[p][s];
				tap1=p;
				tap2=s;
			}
		}
		if(t==0)
		{cout<<"列主元为零，失败！"<<endl;break;}
		if(tap1!=tap2)//换行
		{
			for(int i=0;i<n;i++)////////////////////////////////s
			{
				T=a[tap1][i];
				a[tap1][i]=a[s][i];
				a[s][i]=T;
				T=c[tap1][i];//逆矩阵换行
				c[tap1][i]=c[s][i];
				c[s][i]=T;
			}
		}
		tap1=tap2=0;//置零
		for(int j=0;j<n;j++)//单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响)
		{
			a[s][j]=a[s][j]/t;
			c[s][j]=c[s][j]/t;
		}		
		for(int i=0;i<n;i++)//消元
		{
			if(i!=s)
			{
				m[i]=a[i][s]/a[s][s];
			    for(int j=0;j<n;j++)//注意：逆矩阵每列都消元，而原矩阵则不受影响
				{
					a[i][j]=a[i][j]-m[i]*a[s][j];
					c[i][j]=c[i][j]-m[i]*c[s][j];
				}
			}
		}
	}
}
void xiang_cheng(double a[n][n],double b[n],double c[n])//n*n矩阵乘以n*1矩阵
{
	int i,j,k;
	if(n==n)
	{
		//cout<<"两矩阵阶数相符可以相乘！"<<endl;
		for(i=0;i<n;i++)
		{
			c[i]=0;
		}
		for(k=0;k<n;k++)//c矩阵的第k行
		{
			for(j=0;j<n;j++)
			{
				c[k]+=a[k][j]*b[j];
			}
		}
	}
	else
		cout<<"两个矩阵阶数不符，不能相乘"<<endl;
}
void newton(double x[n])
{
	int i;
	double g0[n];//梯度的值
	double c_H=0.000001;//H准则限值
	double Hesse[n][n];//Hesse矩阵
	double inv_Hesse[n][n];//Hesse矩阵的逆
	double temp[n];//保存矩阵相乘结果

	g_fun(x,g0);//返回梯度的初值		
	do
	{
		Hesse_fun(x,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变
		inv_Hesse[0][1]=inv_Hesse[1][0]=0;
		inv_Hesse[0][0]=1/Hesse[0][0];
		inv_Hesse[1][1]=1/Hesse[1][1];
		xiang_cheng(inv_Hesse,g0,temp);
	    for(i=0;i<n;i++)
		{
			x[i]=x[i]-temp[i];
		}
		g_fun(x,g0);
	}while(H(g0,c_H)==0);
	//cout<<"使用次"<<p<<"牛顿法"<<endl;
	p++;
}