lp.cpp

工程优化和数值计算中常用的优化程序

#include<stdio.h>
#include<math.h>
#include<iostream.h>

/////////////////////////////////////////////////////////
/* 求换入基 */
int Rj(double **matrix, int s,int n)
{
	int i;

	for(i=0;i<s;i++)
		if(fabs(matrix[n][i])>=0.000001)
			if(matrix[n][i]<0) return 0;
	return 1;		
}

/////////////////////////////////////////////////////////
/* 求换入基 */
int Min(double **matrix, int s,int n)
{
	int i,temp=0;
	double min=matrix[n][0];

	for(i=1;i<s;i++)
		if(min>matrix[n][i])
		{
			min=matrix[n][i];
			temp=i;
		}
	return temp;
}


/////////////////////////////////////////////////////////
/* 判断无界情况 */
int Check(int in, double **matrix, int s,int n)
{
	int i;
	double max1=-1;

	for(i=0;i<n;i++)
		if(fabs(matrix[i][in])>=0.000001&&max1<matrix[i][s]/matrix[i][in])
			max1=matrix[i][s]/matrix[i][in];
	if(max1<0)
		return 1;
	return 0;
}

/////////////////////////////////////////////////////////
/* 求换出基 */
int SearchOut(int *temp,int in, double **matrix, int s, int *a,int n) 
{
	int i;
	double min=10000;

	for(i=0;i<n;i++)
		if(fabs(matrix[i][in])>=0.000001&&(matrix[i][s]/matrix[i][in]>=0)&&min>matrix[i][s]/matrix[i][in])
		{
			min=matrix[i][s]/matrix[i][in];
			*temp=i;
		}
	for(i=0;i<s;i++)
		if(a[i]==1&&matrix[*temp][i]==1) return i;
	return -1;	////!!!!!!!!
}

/////////////////////////////////////////////////////////
/* 主元化1 */
void Mto(int in,int temp, double **matrix, int s)
{
	int i;

	for(i=0;i<=s;i++)
		if(i!=in)
			matrix[temp][i]=matrix[temp][i]/matrix[temp][in];
	matrix[temp][in]=1;
}

/////////////////////////////////////////////////////////
/* 初等变换 */
void Be(int temp,int in, double **matrix, int s,int n)
{
	int i,j;

	double c;
	for(i=0;i<=n;i++)
	{
		c=matrix[i][in]/matrix[temp][in];
		if(i!=temp)
			for(j=0;j<=s;j++)
				matrix[i][j]=matrix[i][j]-matrix[temp][j]*c;
	}
}

/////////////////////////////////////////////////////////
/* 改变a[]的值 */
void Achange(int in,int out,int *a)
{
	int temp=a[in];

	a[in]=a[out];
	a[out]=temp;
}


/////////////////////////////////////////////////////////
void Process(double c[][100],int row,int vol,int n)
{
	int i;
	
	for(i=0;i<n;i++)
		if(i!=row) c[i][vol]=0;
}


////////////////////////////////////////////////////
//线性规划求解
double solve_lp(int m,				//变量数
			   int n,				//约束数		xi>=0不算,默认xi>=0
			   double **matrix,		//上面n行约束左边系数阵,最后一行是目标函数的系数
			   double *b,			//约束右边
			   int *code,			//约束不等式符号  0:<=		 1:=		2:>=
			   int type				//求最大/小值	  0:最小     1:最大
			   )
{
	int a[100];					/* 记录基础,非基础的解的情况,0:非基础,1:基础 */
	double x[100];				/* 记录解的数组 */
	int indexe,indexl,indexg;	/* 剩余变量,松弛变量,人工变量 */ 
	int s;						//标准型中变量个数

	int i,j,k;
	double nget[100][100],nlet[100][100],net[100][100]; /* 剩余变量数组,松弛变量数组,人工变量数组 */

	int in,out,temp=0;

	// 化标准形 ////////////////////////////////////////////////////////////
	// 加入变量
	indexe=indexl=indexg=0;
	for(i=0;i<n;i++)
	{
		if(code[i]==0)
		{
			nlet[i][indexl++]=1; 
			Process(nlet,i,indexl-1,n);
		}
		if(code[i]==1)
		{
			net[i][indexg++]=1; 
			Process(net,i,indexg-1,n); 
		}
		if(code[i]==2)
		{
			net[i][indexg++]=1;
			nget[i][indexe++]=-1;
			Process(net,i,indexg-1,n); Process(nget,i,indexe-1,n);
		}
	}
	s=indexe+indexl+indexg+m;

	//合并
	for(i=0;i<n;i++)
	{
		for(j=m;j<m+indexe;j++)
			if(nget[i][j-m]!=-1) matrix[i][j]=0;
			else matrix[i][j]=-1;
		for(j=m+indexe;j<m+indexe+indexl;j++)
			if(nlet[i][j-m-indexe]!=1) matrix[i][j]=0;
			else matrix[i][j]=1;
		for(j=m+indexe+indexl;j<s;j++)
			if(net[i][j-m-indexe-indexl]!=1) matrix[i][j]=0;
			else matrix[i][j]=1;
		matrix[i][s]=b[i];
	}
	for(i=m;i<m+indexe+indexl;i++)
		matrix[n][i]=0;
	for(i=m+indexe+indexl;i<s;i++)
		matrix[n][i]=100;
	matrix[n][s]=0;

	//初始化a[] 
	for(i=0;i<m+indexe;i++)		a[i]=0;
	for(i=m+indexe;i<s;i++)		a[i]=1;

	// 消去人工变量
	if(indexg!=0)
	{
		for(i=m+indexe+indexl;i<s;i++)
		{
			for(j=0;j<n;j++)
				if(matrix[j][i]==1)
				{
					for(k=0;k<=s;k++)
						matrix[n][k]=matrix[n][k]-matrix[j][k]*100;
					j=n;
				}
		}
	}

	// 单纯型算法 
	while(1)
	{
		// 基础可行解 
		for(i=0;i<n;i++)
			for(j=0;j<s;j++)
				if(matrix[i][j]==1&&a[j]==1)
				{
					x[j]=matrix[i][s];
					j=s;
				}
			for(i=0;i<s;i++)
				if(a[i]==0) x[i]=0;

		//Print(); /* 打印 */
		//Result(); /* 打印结果 */

		if(!Rj(matrix,s,n)) 
			in=Min(matrix, s,n); /* 求换入基 */
		else
		{
			if(indexg!=0)		 /* 判断人工变量 */
			{
				for(i=m+indexe+indexl;i<s;i++)
				if(fabs(x[i])>=0.000001)
				{
					//printf("No Answer\n");
					return -9999;
				}
			}

			//最后结果
			if(type==0) return(-matrix[n][s]);
			else return(matrix[n][s]);
		}

		if( Check(in,matrix,s,n) )
		{ /* 判断无界情况 */
			printf("No Delimition\n");
			return -9999;
		}
		out=SearchOut(&temp,in,matrix,s,a,n); /* 求换出基 */
		Mto(in,temp,matrix,s); /* 主元化1 */
		Be(temp,in,matrix,s,n); /* 初等变换 */
		Achange(in,out,a); /* 改变a[]的值 */
	}
}

void main()
{
	int i;

	int m=4,n=2,type=0;
	
	double **matrix=new double*[n+1];
	for (i=0; i<n+1; i++)		matrix[i] = new double[m];

	double *b=new double[n];
	int *code=new int[n];

	b[0]=6; b[1]=5;
	code[0]=1; code[1]=1;
	matrix[0][0]=2;matrix[0][1]=0;matrix[0][2]=-4;matrix[0][3]=1;
	matrix[1][0]=-1;matrix[1][1]=1;matrix[1][2]=3;matrix[1][3]=0;0
	matrix[2][0]=1;matrix[2][1]=-3;matrix[2][2]=2;matrix[2][3]=4;

	cout<<solve_lp(m,n,matrix,b,code,type);


}