算法问题之线性规划（单纯形算法）

与网上大多数单纯形算法不同，本文基本遵循《算法导论》这本书的内容。

基本概念

        在给定有限资源和竞争约束条件下，很多问题都可以表述为最大化或最小化某个目标。如果可以把目标指定为某些变量的一个线性函数，而且如果可以把资源的约束指定为这些变量的等式或不等式，则得到一个线性规划问题（linear-programming problem）。例如，加入你所在的公司生产两种产品：衬衫和手提袋。衬衫每件利润2美元，需要消耗1米布料和5粒扣子；手提袋每个利润3美元，需要消耗2米布料和2粒扣子。你有11米布料和20粒扣子，为了最大限度提高利润，该生产多少件衬衫、多少个手提袋呢？
        为了描述线性规划的性质和算法，使用规范的形式来表示。一般有标准的和松弛的两种形式。在标准型中的线性规划是约束为线性不等式的线性函数最大化。

而松弛型的线性规划是约束为线性等式的线性函数的最大化。

通常使用标准型来表示线性规划，但当描述单纯形算法的细节时，使用松弛型会比较方便。对于变量个数为2的问题，我们可以用图解法。

单纯形算法

        单纯形算法是求解线性规划的经典方法，虽然它的执行时间在最坏的情况下是非多项式的（指数时间复杂度），但是，在绝大部分情况下或者说实际运行过程中却是多项式时间。首先介绍单纯形（simplex）。如果有n个变量，每个约束定义了n维空间中的一个半空间，这些半空间的交集形成了可行区域称作为单纯形（simplex）。目标函数称为一个超平面，而且因为它的凸性，有一个最优解在单纯形的一个顶点上取得。
        单纯形算法以一个线性规划作为输入，输出它的一个最优解。它从单纯形的某个顶点开始，执行一系列的迭代。在每次迭代中，它沿着单纯形的一条边从当前定点移动一个目标值不小于（通常是大于）当前顶点的相邻顶点。当达到一个局部的最大值，即一个顶点的目标值大于其所有相邻顶点的目标值时，单纯形算法终止。因为可行区域是凸的而且目标函数是线性的额，所有局部最优解事实上就是全局最优解。考虑下列标准型的线性规划：

为了利用单纯形算法，必须将线性规划转成松弛型。

左边x4、x5和x6是基本变量，右边x1、x2和x3是非基本变量。基本解为（x1，x2，x3，x4，x5，x6）=（0，0，0，30，24，36），目标解z=0。在每次迭代中，目标是重新表达线性规划，来让基本解有一个更大的目标值。选择一个在目标函数中系数为正值的非基本变量Xe，而且尽可能增加Xe的数值而不违反任何约束。变量Xe变成基本变量，某个其他变量Xl变成非基本变量。其他基本变量和目标函数的值都可能改变。
        继续上面例子，让我们来考虑增加X1的值。当增加X1时，X4、X5、X6的值随之减小。因为对每个变量有一个非负的约束，所有不能允许它们之中的任何一个变成负值。如果X1增加到30以上，则X4称为负值，而当x分别增加到12和9时，X5和X6也成为负值。第三个约束是最紧的约束，因此我们互换X1和X6。

上面的操作为主元（pivot），一个主元选取一个非基本变量Xe和一个基本变量Xl，然后交换二者的角色。基本解为（x1，x2，x3，x4，x5，x6）=（9，0，0，21，6，0），目标解z=27。继续上面例子，将X3和X5互换。

基本解为（x1，x2，x3，x4，x5，x6）=（33/4，0，3/2，69/4，0，0），目标解z=111/4。继续上面例子，将X2和X3互换。

基本解为（x1，x2，x3，x4，x5，x6）=（8，4，0，18，0，0），目标解z=28。此时，目标函数中所有的系数都是负的，得到最优解。因此解是x1=8，x2=4，x3=0且目标值为28。

示例演示

        实际上实现单纯形算法还需要考虑无界、无解、退化等问题。在旋转过程中，可能会存在保持目标值不变的情况，这种现象称为退化。如何避免？一种方法是使用Bland规则：

• 换入变量Xe：目标条件中，第一个系数为正的作为换入变量
• 换出变量Xl：对所有的约束条件中，选择对Xe约束最紧的作为换出变量
  如何找到一个初始的基本可行解？如果所以bi>=0，说明初始的基本解就是基本可行解；若有bi<0，为了确定是否有可行解，可以指定一个辅助线性规划。如下图所示
  

详细内容请见《算法导论》这本书，下面按照书中的伪代码实现单纯形算法。

#include <iostream>
#include <vector>

/*
The Simplex algorithm aims to solve a linear program - optimizing a linear function subject
to linear constraints. As such it is useful for a very wide range of applications.

N.B. The linear program has to be given in *slack form*, which is as follows:
maximize
c_1 * x_1 + c_2 * x_2 + ... + c_n * x_n + v
subj. to
a_11 * x_1 + a_12 * x_2 + ... + a_1n * x_n + b_1 = s_1
a_21 * x_1 + a_22 * x_2 + ... + a_2n * x_n + b_2 = s_2
...
a_m1 * x_1 + a_m2 * x_2 + ... + a_mn * x_n + b_m = s_m
and
x_1, x_2, ..., x_n, s_1, s_2, ..., s_m >= 0

Every linear program can be translated into slack form; the parameters to specify are:
- the number of variables, n, and the number of constraints, m;
- the matrix A = [[A_11, A_12, ..., A_1n], ..., [A_m1, A_m2, ..., A_mn]];
- the vector b = [b_1, b_2, ..., b_m];
- the vector c = [c_1, c_2, ..., c_n];
- the variable v : free constant of objective function;                  

Complexity:    O(m^(n/2)) worst case
O(n + m) average case (common)
*/

typedef std::vector <std::vector<float>> Matrix;
typedef std::vector<float> CofVector;
typedef	std::vector<int> IndexVector;

class SimplexAlgorithm
{
public:
	enum ResultType
	{
		FEASIBLE,
		INFEASIBLE,
		UNBOUNDED
	};

	void simplex(const Matrix &A_input, const CofVector &b_input, const CofVector &c_input)
	{
		//check data valid
		if (A_input.empty() || b_input.empty() || c_input.empty())
			return;
		if (A_input.size() != b_input.size() || A_input[0].size() != c_input.size())
			return;
		PrintLinearProgramming(A_input, b_input, c_input);

		float v = .0f;
		Matrix A = A_input;
		CofVector b = b_input, c = c_input;
		IndexVector N(A[0].size(), 0), B(A.size(), 0);
		if (!initialzieSimplex(N, B, A, b, c, v))
		{
			_resulttype = INFEASIBLE;
			PrintResult();
			return;
		}
		if (!iterateSimplex(N, B, A, b, c, v))
		{
			_resulttype = UNBOUNDED;
			PrintResult();
			return;
		}
	
		_resulttype = FEASIBLE;
		_x.clear();
		for (int i = 0; i < N.size(); i++)
		{
			_x.push_back(.0f);
			for (int j = 0; j < B.size(); j++)
			{
				if (i == B[j])
					_x[i] = b[j];
			}
		}
		_v = v;
		PrintResult();
	}

private:
	void PrintLinearProgramming(const Matrix &A, const CofVector &b, const CofVector &c)
	{
		std::cout << "Linear Programming:" << std::endl;
		std::cout << "Maximize:  " << std::endl;
		for (int i = 0; i < c.size(); i++)
		{
			std::cout << " " << c[i] << "x" << i + 1;
			if (i != c.size() - 1)
				std::cout << " +";
		}
		std::cout << std::endl;

		std::cout << "Constraint: " << std::endl;
		for (int j = 0; j < b.size(); j++)
		{
			for (int i = 0; i < c.size(); i++)
			{
				std::cout << " " << A[j][i] << "x" << i + 1;
				if (i != c.size() - 1)
					std::cout << " +";
			}
			std::cout << " " << "=<" << " " << b[j] << std::endl;
		}
	}

	void PrintResult()
	{
		std::cout << "Result:  " << std::endl;
		switch (_resulttype)
		{
		case SimplexAlgorithm::FEASIBLE:
			std::cout << "z = " << _v << " ";
			for (int i = 0; i < _x.size(); i++)
				std::cout << "x" << i + 1 << " = " << _x[i] << " ";
			break;
		case SimplexAlgorithm::INFEASIBLE:
			std::cout << "infeasible";
			break;
		case SimplexAlgorithm::UNBOUNDED:
			std::cout << "unbounded";
			break;
		default:
			break;
		}
		std::cout << std::endl;
	}

	bool initialzieSimplex(IndexVector &N, IndexVector &B, Matrix &A, CofVector &b, CofVector &c, float &v)
	{
		//let k be the index of the minimum bi
		int k = 0;
		double min_b = b[0];
		for (int i = 1; i < B.size(); i++)
		{
			if (b[i] < min_b)
			{
				k = i;
				min_b = b[i];
			}
		}

		//is the initial basic solution feasible
		if (b[k] >= 0)
		{
			for (int i = 0; i < N.size(); i++)
				N[i] = i;
			for (int i = 0; i < B.size(); i++)
				B[i] = (int)N.size() + i;
			return true;
		}

		//form auxiliary LP by adding x0 to the left-hand side of each equation
		//and setting the objective function to -x0
		CofVector c_aux(N.size(), .0f); // -x0
		c_aux.push_back(-1.0f) ;
		CofVector b_aux = b;
		IndexVector N_aux;
		for (int i = 0; i <N.size(); i++)
			N_aux.push_back(i);
		IndexVector B_aux;
		for (int i = 0; i < B.size(); i++)
			B_aux.push_back((int)N.size() + i);
		//the index of x0
		const int kIndexX0 = B_aux.back() + 1;
		N_aux.push_back(kIndexX0);
		Matrix A_aux = A;
		for (int i = 0; i < B.size(); i++)
			A_aux[i].push_back(-1);

		int indexl_aux = k ;
		int indexe_aux = (int)N_aux.size() - 1;
		float v_aux = .0f;
		//the basic solution is now feasible for auxiliary LP
		pivot(N_aux, B_aux, A_aux, b_aux, c_aux, v_aux, indexl_aux, indexe_aux);
		//find an optimal solution
		iterateSimplex(N_aux, B_aux, A_aux, b_aux, c_aux, v_aux);

		//if the basic solution sets x0 = 0  
		if (v_aux == .0f) 
		{
			N.clear();
			for (auto &i : N_aux)
			{
				if (i != kIndexX0)
					N.push_back(i);
			}
			B = B_aux;
			for (int i = 0; i < A.size(); i++)
			{
				A[i].clear();
				for (int j = 0; j < N_aux.size(); j++)
				{
					if(N_aux[j] != kIndexX0)
						A[i].push_back(A_aux[i][j]);
				}
			}
			b = b_aux;
			
			CofVector c_temp(c.size(), .0f);
			//original objective function
			for(int i = 0; i < N.size(); i++) 
			{
				bool flag = false;
				for (int j = 0; j < N.size(); j++)
				{
					if (N[j] == i)
					{
						c_temp[j] = c[i];
						flag = true;
						break;
					}
				}
				if(flag)
					continue;

				for (int j = 0; j < B.size(); j++)
				{
					if (B[j] == i)
					{
						for (int m = 0; m < N.size(); m++)
							c_temp[m] += -c[i] * A[j][m];
						v += c[i] * b[j];
					}
				}	
			}
			c = c_temp;
			return true;
		}
		return false;
	}

	void pivot(IndexVector &N, IndexVector &B, Matrix &A, CofVector &b, CofVector &c, float &v, int indexl, int indexe)
	{
		IndexVector N_input = N, B_input = B;
		Matrix A_input = A;
		CofVector b_input = b, c_input = c;
		int indexl_input = indexl, indexe_input = indexe;
		indexe = indexl_input;  indexl = indexe_input;
		//compute the coefficients of the equation for new basic variable Xe
		b[indexe] = b_input[indexl_input] / A_input[indexl_input][indexe_input];
		for (int j = 0; j < N.size(); j++)
		{
			if (j != indexe_input)
				A[indexe][j] = A_input[indexl_input][j] / A_input[indexl_input][indexe_input];
		}
			
		A[indexe][indexl] = 1.0f / A_input[indexl_input][indexe_input];

		//Computer the coefficients of the remaining constraints
		for (int i = 0; i < B.size(); i++)
		{
			if (i != indexl_input) 
			{
				b[i] = b_input[i] - A_input[i][indexe_input] * b[indexe];
				for (int j = 0; j < N.size(); j++)
				{
					if (j != indexe_input)
						A[i][j] = A_input[i][j] - A_input[i][indexe_input] * A[indexe][j];
				}
				A[i][indexl] = -A_input[i][indexe_input] * A[indexe][indexl];
			}		
		}

		//Compute the objective function
		v = v + c_input[indexe_input] * b[indexe];
		for (int j = 0; j < N.size(); j++)
		{
			if (j != indexe_input)
				c[j] = c_input[j] - c_input[indexe_input] * A[indexe][j];
		}
		c[indexl] = -c_input[indexe_input] * A[indexe][indexl];

		//Compute new sets of basic and nonbasic variables
		B[indexl_input] = N_input[indexe_input];
		N[indexe_input] = B_input[indexl_input];
	}

	bool iterateSimplex(IndexVector &N, IndexVector &B, Matrix &A, CofVector &b, CofVector &c, float &v)
	{
		//find  index E
		int indexe = 0, indexl = 0;
		while (findIndexE(N, c, indexe))
		{
			if (findIndexL(B, A, b, indexl, indexe))
				pivot(N, B, A, b, c, v, indexl, indexe);
			else	
				return false;
		}
		return true;
	}

	bool findIndexE(const IndexVector &N, const CofVector &c, int &indexe)
	{
		indexe = -1;
		int min_e = INT_MAX;
		for (int i = 0; i < c.size(); i++)
		{
			if (c[i] > .0f && N[i] <  min_e)
			{
				indexe = i;
				min_e = N[i];
			}
		}
		if (indexe == -1)
			return false;
		return true;
	}

	bool findIndexL(const IndexVector &B, const Matrix &A, const CofVector &b, int &indexl, int indexe)
	{
		indexl = -1;
		float min_constr = FLT_MAX;
		for (int i = 0; i < B.size(); i++)
		{
			if (A[i][indexe] > .0f &&  b[i] / A[i][indexe] < min_constr)
			{
				indexl = i;
				min_constr = b[i] / A[i][indexe];
			}
		}
		if (indexl == -1)
			return false;
		return true;
	}

private:
	ResultType _resulttype = FEASIBLE;
	CofVector _x;
	float _v;

};

int main(int argc, char* argv[])
{
	/*	Simplex tests  */

	//Basic solution feasible
	//z=8, x1=8, x2=4, x3=0
	Matrix A = { { 1, 1, 3},
				 { 2, 2, 5},
				 { 4, 1, 2} };
	CofVector b = { 30, 24, 36 };
	CofVector c = { 3, 1, 2};

	unbounded;
	//Matrix A = { { 1, -1 },
	//			 { -1, 1 }};
	//CofVector b = { 1, 1};
	//CofVector c = { 1, 1};

	Basic solution infeasible but LP feasible;
	z=2, x1=14/9, x2=10/9
	//Matrix A = { { 2, -1 },
	//{ 1, -5 } };
	//CofVector b = { 2, -4 };
	//CofVector c = { 2, -1 };

	Basic solution infeasible but LP feasible;
	z=-1, x1=1, x2= 0
	//Matrix A = { { 1, 1 },
	//		     { -1, -1 } };
	//CofVector b = { 2, -1 };
	//CofVector c = { -1, -2 };

	Basic solution infeasible but LP infeasible;
	//Matrix A = { { -1, -1 },
	//		     { 1, 1 }};
	//CofVector b = { -2, 1 };
	//CofVector c = { -1, -2 };

	
	
	SimplexAlgorithm simplex;
	simplex.simplex(A, b, c);
	system("pause");
	return 0;
}

运行结果

参考资料

• 算法导论[M]
• 线性规划-单纯形算法详解