dualprog.c

对偶单纯形法求线性规划最优解

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

//  对偶单纯形法通用处理函数库
//  Author: lubangjian
//	Date:   2004/11/20

#define SUCCESS	1
#define FAIL	0

/*分配内存*/
void *_malloc(char *position,int memsize)
{
	void *mem;
	if ((mem=malloc(memsize))==NULL)
	{
		printf("%s:内存分配失败\n",position);
		return NULL;
	}
	return mem;
}

//功能说明:取离基变量
//入参说明:p--限制条件参数矩阵
//		   b--基解
//		   row--矩阵列数
//		   col--矩阵行数
//返 回 值:>0--进基变量列号,-1--找不到进基变量,-2--当前已为最优解
int get_out_para(double *b,double *p,int row,int col)
{
	int colcount,rowcount,rowmark;
	
	//取最小判别数值
	for (rowcount=1,rowmark=0;rowcount<row;rowcount++)
		if (b[rowmark]>b[rowcount])
			rowmark=rowcount;	
	if (b[rowmark]>=0.0)//当前找到最优解
		return -2;
	//检查是否进基变量系数矩阵所在列是否全小于0
	for (colcount=0;colcount<col;colcount++)
	{
		if (p[rowmark*col+colcount]<0.0)
			break;
	}
	if (colcount>=col)//当前不可能有最优解
		return -1;
	return rowmark;//返回离基变量所在行号
}


//功能说明:计算进基变量值,并由此确定离基变量
//入参说明:entp--进基变量参数
//		   p--限制条件参数矩阵
//		   b--基解
//		   sigma--进基变量判别数存放区
//		   row--矩阵列数
//		   col--矩阵行数
int get_ent_para(int outp,double *p,double *b,double *sigma,int row,int col)
{
	int colcount,find=0;
	double maxsigma=0;
	int colmark;

	for (colcount=0;colcount<col;colcount++)
	{
		if (p[outp*col+colcount]<0)
		{
			find++;
			if ((find==1)||(maxsigma<sigma[colcount]/p[outp*col+colcount]))//为第一次计算或找到更大值
			{
				maxsigma=sigma[colcount]/p[outp*col+colcount];
				colmark=colcount;
			}
		}
	}

	return colmark;//返回进基变量所在列号
}

//功能说明:输出最优化结果
int output(int *B,double *b,double *cb,int row,int col,double *solve,int printres)
{
	int rowcount,colcount;
	double funcval;
	char *position="output";

	for (colcount=0;colcount<col;colcount++)
		solve[colcount]=0.0;
	
	for (rowcount=0,funcval=0;rowcount<row;rowcount++)
	{
		solve[B[rowcount]]=b[rowcount];//最优解
		funcval+=cb[rowcount]*b[rowcount];
	}

	if (printres)//需要输出结果
	{
		printf("最优解:(");
		for (colcount=0;colcount<col;colcount++)
		{
			if (colcount>0)//此处默认为第一个变量不可能为大M
				printf(",");
			printf("%f",solve[colcount]);
		}
		printf(")\n函数值:%f\n",funcval);
	}
	return SUCCESS;
}

//功能说明:输出中间结果
int output_mid(double *cj,
			   double *b,
			   int *B,
			   double *p,
			   int row,
			   int col,
			   double *sigma,
			   double *cb,
			   double *cjo)
{
	int rowcount,colcount;

	printf("\tcj   \t|");
	for (colcount=0;colcount<col;colcount++)
		printf(" %f",cj[colcount]);
	printf("\n");
	for (rowcount=0;rowcount<row;rowcount++)
	{
		printf("%f|P%d|%f|",cb[rowcount],B[rowcount]+1,b[rowcount]);
		for (colcount=0;colcount<col;colcount++)
			printf(" %f",p[rowcount*col+colcount]);
		printf("\n");	
	}
	printf("\tsigma  \t|");
	for (colcount=0;colcount<col;colcount++)
		printf(" %f",sigma[colcount]);
	printf("\n\n");
	return 0;
}

#define Free_Alocated_mem	if (cb) free(cb);\
							if (cjo) free(cjo);\
							if (sigma) free(sigma);\
							if (b) free(b);\
							if (B) free(B);\
							if (p) free(p)
//功能说明:大M单纯形法
//入参说明:cj--原形参数列表
//		   ent_b--初始基解
//		   ent_B--基变量参数
//		   ent_p--限制条件参数矩阵
//		   row--矩阵列数
//		   col--矩阵行数
//		   printres--是否打印结果
//		   needmid--是否需要中间结果
//		   solve--解
//返 回 值:SUCCESS--成功,FAIL--失败
int dualprog(double *cj,double *ent_b,int *ent_B,double *ent_p,int row,int col,int printres,int needmid,double *solve)
{
	double *cb=NULL;
	double *cjo=NULL;
	int entp,outp;//进基变量,离基变量
	double *sigma=NULL;
	double tempo,tempb;
	int res;
	int colcount,rowcount;
	double *b=NULL,*p=NULL;
	int *B=NULL;
	char *position="dualprog";

	if (!(b=(double *)_malloc(position,sizeof(double)*row)))
		return FAIL;
	if (!(B=(int *)_malloc(position,sizeof(int)*row)))
	{
		Free_Alocated_mem;
		return FAIL;
	}
	if (!(p=(double *)_malloc(position,sizeof(double)*row*col)))
	{
		Free_Alocated_mem;
		return FAIL;
	}
	if (!(sigma=(double *)_malloc(position,sizeof(double)*col)))//为判别数分配空间
	{
		Free_Alocated_mem;
		return FAIL;
	}
	if (!(cb=(double *)_malloc(position,sizeof(double)*row)))//为CB分配空间
	{
		Free_Alocated_mem;
		return FAIL;
	}
	memcpy(b,ent_b,sizeof(double)*row);
	memcpy(B,ent_B,sizeof(int)*row);
	memcpy(p,ent_p,sizeof(double)*row*col);

	for (rowcount=0;rowcount<row;rowcount++)
		cb[rowcount]=cj[B[rowcount]];
	//计算初始判别数
	for (colcount=0;colcount<col;colcount++)
	{
		tempo=0.0;
		for (rowcount=0;rowcount<row;rowcount++)//计算一列的判别数
			tempo+=p[rowcount*col+colcount]*cj[B[rowcount]];
		sigma[colcount]=cj[colcount]-tempo;
	}
	while(1)//逐次迭代取最优解
	{
		res=get_out_para(b,p,row,col);
		//检查基解
		
		if (res==-2)//b>0,当前已找到最优解
		{
			if (needmid)
			{
				//printf("entp=%d,outp=%d\n",entp+1,outp+1);
				output_mid(cj,b,B,p,row,col,sigma,cb,cjo);
			}
			output(B,b,cb,row,col,solve,printres);//输出最优结果
			return SUCCESS;
		}
		else if (res==-1)//当前没有最优解
		{
			return FAIL;
			printf("无法找到最优解\n");
		}
		else
			outp=res;
		entp=get_ent_para(outp,p,b,sigma,row,col);//计算进基变判别数,取得进基变量
		if (needmid)//需要输出中间结果
			output_mid(cj,b,B,p,row,col,sigma,cb,cjo);
		B[outp]=entp;//更新B值
		cb[outp]=cj[entp];//更新cb值
		//计算非离基行参数矩阵值
		tempb=b[outp]/p[outp*col+entp];
		for (rowcount=0;rowcount<row;rowcount++)
		{
			if (rowcount==outp)//当前为离基变量所在行
			{
				b[rowcount]/=p[outp*col+entp];/*计算与离基变量在同行的b值*/
				continue;
			}
			b[rowcount]-=tempb*p[rowcount*col+entp];/*计算与离基变量不同行b值*/
			tempo=p[rowcount*col+entp]/p[outp*col+entp];
			for (colcount=0;colcount<col;colcount++)
				p[rowcount*col+colcount]-=p[outp*col+colcount]*tempo;
		}
		//计算系数矩阵离基量所在行的值和新的判别数
		tempo=sigma[entp]/p[outp*col+entp];
		for (colcount=0;colcount<col;colcount++)
		{
			sigma[colcount]-=p[outp*col+colcount]*tempo;//新的判别数值
			if (colcount!=entp)
				p[outp*col+colcount]/=p[outp*col+entp];//与进基变量不在同列的系数矩阵值
		}
		p[outp*col+entp]=1.0;//
	}

	Free_Alocated_mem;
	return FAIL;
}

//调用例子
int main(int argc,char **argv)
{
	int needmid=1;//默认为不需中间结果
	double cj[6]={12,8,16,12,0,0},solve[6];
	double b[2]={-2,-3};
	int B[2]={4,5};
	double p[2][6]={{-2,-1,-4,0,1,0},{-2,-2,0,-4,0,1}};

	if (argc>1)
		needmid=atoi(argv[1]);

	dualprog(cj,b,B,*p,2,6,1,needmid,solve);
}