linear.cpp

组合数学中单纯形算法解决一般线性规划问题

#include  <conio.h>   
#include  <math.h>  

#include  <iostream>  
using namespace std;


//A-----------等式约束的系数矩阵   
//B-----------等式约束右端常数   
//C-----------所求最小值表达式的各项系数   
//m-----------系数矩阵行数   
//n-----------系数矩阵列数
//Result------返回最小值   
bool Line_Optimize1(double *A,double *c,int *X_nonbase,int *X_base,double *B,int m,int n,double &result)
{
	//double *B = new double[m];
	double *a = new double[m * m];

	double *tmp_a = new double[m * m];
	double *tmp_B = new double[m];
	double *tmp_c = new double[m];
	
	for(int i=0; i < m; i++)
	{
		for( int j=0; j < m; j++)
		{
			a[i * m + j] = A[i * n + X_base[j]] / A[ i * n + X_nonbase[i]];
		}
		B[i] = B[i] / A[i*n + X_nonbase[i]];
	}

	bool bSuccess = false;
	while(!bSuccess)
	{
		int k_cIndex = -1;
		int e_bIndex = 0;

		//  确定换入变量
		double c_max = 0;
		int i;
		for( i = 0; i< m; i++)
		{
			if(c[i] > c_max)
			{
				c_max = c[i];
				k_cIndex = i;
			}
		}

		//不存在换入变量c > 0则结束
		if(k_cIndex < 0)
		{
			bSuccess = true;
			break;
		}

		// 确定换出变量
		e_bIndex = -1;
		double min = 1E10;
		for( i = 0; i < m; i++)
		{
			if(a[i * m + k_cIndex] > 0)
			{
				double div_tmp = B[i] / a[i * m + k_cIndex] ;
				if( div_tmp < min)
				{
					e_bIndex = i;
					min = div_tmp;
				}
			}
		}

		//不存在换出变量a[i][e] > 0则结束
		if(e_bIndex < 0)
		{
			bSuccess = true;
			break;
		}

		int ch_temp = X_base[k_cIndex];
		X_base[k_cIndex] = X_nonbase[e_bIndex];
		X_nonbase[e_bIndex] = ch_temp;

		//转轴变换
		for( i = 0; i < m; i++)
		{
			for(int j = 0; j < m; j++)
			{
				tmp_a[i * m +j] = a[i * m +j];
			}

			tmp_B[i] = B[i];
			tmp_c[i] = c[i];
		}
		// b:
		for( i = 0; i < m; i++)
		{
			// b
			if(i != e_bIndex)
				B[i] = tmp_B[i] - a[ i * m + e_bIndex] * tmp_B[k_cIndex] / tmp_a[k_cIndex * m + e_bIndex];
			else
				B[i] = tmp_B[k_cIndex] / tmp_a[k_cIndex * m + e_bIndex];

			// c
			if( i != k_cIndex)
				c[i] = tmp_c[i] - tmp_c[e_bIndex] * tmp_a[ k_cIndex * m + i] / tmp_a[k_cIndex * m + e_bIndex];
			else
				c[k_cIndex] = - tmp_c[e_bIndex] / tmp_a[k_cIndex * m + e_bIndex];
		}

		// a
		for( i = 0; i < m; i++)
		{
			for(int j = 0; j < m; j++)
			{
				if( i == e_bIndex && j == k_cIndex)
					a[ i * m + j] = 1 / tmp_a[ k_cIndex * m + e_bIndex];
				else if( j == k_cIndex)
					a[ i * m + k_cIndex] =  - tmp_a[i * m + e_bIndex] / tmp_a[ k_cIndex * m + e_bIndex];
				else if( i == e_bIndex)
					a[ k_cIndex * m + j] =  tmp_a[k_cIndex * m + j] / tmp_a[ k_cIndex * m + e_bIndex];
				else
					a[ i * m + j] = tmp_a[ i * m + j] - tmp_a[i *m + e_bIndex] * tmp_a[k_cIndex * m + j] / tmp_a[k_cIndex * m + e_bIndex];
			}
		}

		//计算 z
		result = 0;
		for( i=0; i < m; i++)
		{
			if(X_nonbase[i] == 1)
				result += -B[i];

			if(X_nonbase[i] == 2)
				result += 3 * B[i];

			if(X_nonbase[i] == 4)
				result += -2 * B[i];
		}
	}

	delete []a ;//= new double[m * m];

	delete []tmp_a ;//= new double[m * m];
	delete []tmp_B ;//= new double[m];
	delete []tmp_c ;//= new double[m];

	return true;
}

// maxz = -X2 + 3X3 - 2X5
// s.t. X1 + 3X2 - X3 + 2X5 = 7
//		-2X2 + 4X3 + X4 = 12;
//		-4X2 + 3X3 + 8X5 + X6 = 10;
//
int main()
{
	double  A[]={   
		1,3,-1,0,2,0,   
		0,-2,4,1,0,0,   
		0,-4,3,0,8,1   
	};

	int X_nonbase[] = {0,3,5};	//初始离基变量索引
	int X_base[] = {1,2,4};		//初始入基变量索引
	double C[] = {-1,3,-2};
	double B[] = {7,12,10};
	double result = 0;			//目标函数值

	Line_Optimize1(A,C,X_nonbase,X_base,B,3,6,result);

	for  (int  i=0;  i<3;  i++)
	{   
		cout<<"The  No. X"<<X_nonbase[i]+1<<" = "<<B[i]<<endl;   
	}   

	cout  <<  "The  other  root  is  Zero."<<endl;   
	cout  <<  "The  Minmize  value  is:  "<<result<<endl;

	return 0;
}