线性规划之单纯型算法

问题定义：

问题定义比较复杂，建议看《算法导论》里的线性规划一章。单纯型算法用于求解如下这类问题：

例：

求等式的最小值： -2X₁– 3X₂

且自变量满足如下约束：

X₁ + X₂ = 7

_X1 – 2X₂<= 4

_X1>= 0

将约束等式转换为标准型：

标准型的条件：

1. 求目标函数的最大值

2. 每个自变量都大于等于零（非负约束）

3.约束不等式，只有最小化约束

转换结果如下：

max 2X₁ – 3X₂ + 3X₃

并且满足：

X₁ + X₂- X₃ <= 7

-X₁ – X₂+ X₃ <= -7

X₁ – 2X₂+ 2X₃ <= 4

X₁, X₂, X₃ >= 0

将标准型转换成松弛型：

z = 2X₁– 3X₂ + 3X₃

X₄ = 7- X₁ - X₂ + X₃

X₅ = -7+ X₁ + X₂ - X₃

X₆ = 4- X₁ + 2X₂ - 2X₃

则基本变量为的下标集合 B = {4, 5, 6}

非基本变量的下标集合 N = {1, 2, 3}

C = [ 2 3 3 ]^T

^{b = [ 7 -7 4 ]T}

A	1	1	-1
-1	-1	1
1	-2	2

v = 0

具体实现时，将 A， b， C， v都存储在一个数组中

data	0v	2c	-3c	3c
7b	1A	1A	-1A
-7b	-1A	-1A	1A
4b	1A	-2A	2A

代码如下：

/*
*Copyright(c) Computer Science Department of XiaMen University 
*
*Authored by laimingxing on: 2012年 03月 03日 星期六 00:14:35 CST
*
* @desc:
*
* @history
*/
#include <iostream>
#include<algorithm>
#include <vector>
#include <string.h>
using namespace std;

const double unbounded = 10000000;
const int MAX = 10;
int n;
double data[MAX][MAX];//A,b,c,v都存储在data里面

void PIVOT( vector<int> &N, vector<int> &B, int l ,int e);
void simple( vector<int> &N, vector<int>&B);

int main(int argc, char* argv[])
{
	
//将标准行转换为松弛型
	//未知数个数
	const int nX = 3;
	//不等式个数 ,第一行为目标函数
	const int nEquation = 4;
	int i = 0, j = 0, k = 0;
	vector<int>B,N;
	double arr[ nEquation + 1 ][nX + 1 ] = {
	{ 2, -3, 3,0},{1, 1, -1, 7}, {-1, -1, 1, -7}, {1, -2, 2, 4}
				};
	
	//N
	for( k = 1; k <= nX; k++)
		N.push_back( k);
	//B
	for( ; k < nX + nEquation; k++)
		B.push_back( k);

	n = nX + nEquation ; 

	memset( data, 0, sizeof(int) * n * n);

	//c
	for( i = 1; i < n; i++)
		if( i <= nX) data[0][i] = arr[0][i - 1];
		else data[0][i] = 0;
	//A
	for( i = nX + 1; i < n; i++)	
	{
		for( j = 1; j <= nX ; j++)
		{
			data[i][j] = arr[ i - nX][j -1];
			//cout << data[i][j] << "\t";
		}
	}
	//b
	for( i = nX + 1; i < n; i++)
		data[i][0] = arr[ i - nX][nX];
 
	simple( N,B);
	return 0;
}
void simple( vector<int> &N, vector<int>&B)
{
	for( int j = 1; j <= n; j++)
	{
		if(data[0][j] <= 0) continue;

		double minBound = unbounded;
		int minIndex = 0;
		for( vector<int>::iterator i = B.begin(); i != B.end(); i++)	
		{
			if( data[*i][j] <= 0)continue;
			double temp = data[*i][0] / data[*i][j];
			if( temp < minBound)
			{
				minBound = temp;
				minIndex = *i;
			}
		}
		if( minBound > unbounded - 1)
			cout <<"Unbounded" << endl;
		else
			PIVOT( N, B,  minIndex, j);
	}	

	for( int i = 1 ; i <= n ; i++)
		if( find(B.begin(), B.end(),i) != B.end() )
			cout << "x"<<i << "  = " << data[i][0] << endl;
		else
			cout << "x"<<i << "  = " << 0 <<endl;
}

//N是非基本变量，B是基本变量
void PIVOT( vector<int> &N, vector<int> &B, int l ,int e)
{
	data[e][0] = data[l][0] / data[l][e];
	vector<int>::iterator  j;

	for( j = N.begin(); j != N.end(); j++)
	{
		if( *j == e)continue;
		data[e][*j] = data[l][*j] / data[l][e];
	}
	data[e][l] = 1/ data[l][e];

	for( vector<int>::iterator iter = B.begin(); iter != B.end(); iter++)
	{
		if( *iter == l) continue;
		data[*iter][0] -= data[*iter][e] * data[e][0];
		for( vector<int>::iterator it = N.begin(); it != N.end(); it++)
		{
			if( *it == e)continue;
			data[*iter][*it] -= data[*iter][e] * data[e][*it];
		}	
		data[*iter][l] = -1 * data[*iter][l] * data[e][l];
	}
	//Computer the objective function
	data[0][0] += data[0][e] * data[e][0];
	for( j = N.begin(); j != N.end(); j++)
	{
		if( *j == e) continue;
		data[0][*j] -= data[0][e] * data[e][*j];
	}
	data[0][l] = -1 * data[0][e] * data[e][l];
	//Compute new sets of basic and nonbasic variables.
	vector<int>::iterator temp_it = find( N.begin(), N.end(), e);
	if( temp_it != N.end() ) N.erase( temp_it);
	N.push_back(l);

	temp_it = find( B.begin(), B.end(), l);
	if( temp_it != B.end() ) B.erase( temp_it);
	B.push_back(e);

}