与网上大多数单纯形算法不同,本文基本遵循《算法导论》这本书的内容。
基本概念
在给定有限资源和竞争约束条件下,很多问题都可以表述为最大化或最小化某个目标。如果可以把目标指定为某些变量的一个线性函数,而且如果可以把资源的约束指定为这些变量的等式或不等式,则得到一个线性规划问题(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]
- 线性规划-单纯形算法详解