线性规划之单纯形算法

首先给出转动操作的伪代码（摘自算法导论）：

PIVOT

PIVOT((N,B,A,b,c,v,l,e))
let (hat{A}) be a new (m imes n) matrix
(hat{b}_{e}=b_l/a_{le})
for each (jin N-{e})
(hat{a}_{ej}=a_{lj}/a_{le})
(hat{a}_{el}=1/a_{le})
for each (iin B-{l})
(hat{b}_i=b_i-a_{ie}hat{b}_e)
for each (jin N-{e})
(hat{a}_{ij}=a_{ij}-a_{ie}hat{a}_{ej})
(hat{a}_{il}=-a_{ie}hat{a}_{el})
(hat{v}=v+c_ehat{b}_e)
for each (jin N-{e})
(hat{c}_j=c_j-c_ehat{a}_{ej})
(hat{c}_l=-c_ehat{a}_{el})
(hat{N}=N-{e}cup {l})
(hat{B}=B-{l}cup {e})
return((hat{N},hat{B},hat{A},hat{b},hat{c},hat{v}))
C++代码：

//输入为一个松弛型
void PIVOT(int l, int e, int flag) { //l为替出变量下标，e为替入变量下标
	b[e] = b[l]/A[l][e]; //一系列更换约束系数的操作 
	for(int j = flag; j <= n+m; ++j) A[e][j] = A[l][j]/A[l][e];
	for(int i = flag; i <= n+m; ++i) {
		if(i == e) continue;
		b[i] = b[i]-A[i][e]*b[e];
		for(int j = flag; j <= n+m; ++j) {
			if(j == e) continue;
			A[i][j] = A[i][j]-A[i][e]*A[e][j];
		}
		A[i][e] = 0;
	}
	v = v+c[e]*b[e]; //更改目标函数 
	for(int i = flag; i <= n+m; ++i) {
		if(i == e) continue;
		c[i] = c[i]-c[e]*A[e][i];
	}
	c[e] = 0;
	B[e] = 1, B[l] = 0; //因为N是B的补集，不用记录
}

然后是单纯形主体：

SIMPLEX

SIMPLEX((A,b,c))
((N,B,A,b,c,v)=)INITIALIZE-SIMPLEX((A,b,c))
let (Delta) be a new vector of length (m)
while som index (jin N) has (c_j>0)
choose an index (ein N) for which (c_e>0)
for each index (iin B)
if (a_{ie}>0)
(Delta_i=b_i/a_{ie})
else
(Delta_i=infty)
choose an index (lin B) that minimizes (Delta_{l})
if (Delta_{l}=infty)
return "unbounded"
else ((N,B,A,b,c,v)=)PIVOT((N,B,A,b,c,v,l,e))
for (i=1) to (n)
if (iin B)
(overline{x}_i=b_i)
else
(overline{x}_i=0)
return ((overline{x}_1,overline{x}_2,dots,overline{x}_n))

int SIMPLEX() { //-1无解，0无界，1有限最优解
	if(!INITIALIZE_SIMPLEX()) return -1;
	while(1) {
		int e = -1, l = -1;
		for(int i = n+m; i >= 1; --i)
			if(!dcmp(c[i], 0) && c[i] > 0) e = i; //Bland规则
		if(e == -1) break;
		double det = INF;
		for(int i = n+m; i >= 1; --i) {
			if(!B[i] || A[i][e] < 1 || dcmp(A[i][e], 0)) continue;
			if(b[i]/A[i][e] < det || dcmp(det, b[i]/A[i][e])) l = i, det = b[i]/A[i][e];
		}
		if(dcmp(det, INF)) return 0; //无界 
		else PIVOT(l, e, 1);
	}
	return 1; //有限的最优解 
}

最后是初始化可行解：

INITIALIZE-SIMPLE

INITIALIZE-SIMPLEX((A,b,c))
let (k) be the index of the minimum (b_i)
if (b_kgeqslant 0)
return (({1,2,dots,n},{n+1,n+2,dots,n+m},A,b,c,0))
form (L_{aux}) by adding (-x_0) to the left-hand side of each constraint and setting the objective function to (-x_0)
let ((N,B,A,b,c,v)) be the resulting slack form for (L_{aux})
(l=n+k)
((N,B,A,b,c,v)=)PIVOT((N,B,A,b,c,v,l,0))
iterate the while loop of lines 3-12 of SIMPLEX until an optimal solution to (L_{aux}) is found
if the optimal solution to (L_{aux}) sets (overline{x}_0) to 0
if (overline{x}_0) is basic
({color{Red} {perform one (degenerate) pivot to make it nonbasic}})
from the final slack form of (L_{aux}),remove (x_0) from the constraints and restore the original objective function of (L),but replace each basic variable in this objective function by the right-hand side of its associated constraint
return the modified final slack form
else return "infeasible"
P.S.标红色的那句话不懂什么意思啊， 确实是原文，但我看网上的伪代码里都没有这个判断，这种情况真的会出现吗

//传入一个标准型
//要么返回无解，要么传回一个松弛型的基本可行解 
bool INITIALIZE_SIMPLEX() {
	for(int i = n+1; i <= n+m; ++i) B[i] = A[i][i] = 1;
	int k = n+1;
	for(int i = n+1; i <= n+m; ++i) if(b[i] < b[k]) k = i;
	if(b[k] >= 0) { //基本解可行，直接返回 
		v = 0;
		return 1;
	}
	for(int i = n+1; i <= n+m; ++i) A[i][0] = -1;
	static double c0[N+M+5]; //先拷贝A,c并使z=-x0
	memcpy(c0, c, sizeof c);
	memset(c, 0, sizeof c);
	c[0] = -1;
	PIVOT(k, 0, 0); //需要一组使辅助线性规划成立的基本可行解
	while(1) { //寻找辅助线性规划的最优解 
		int e = -1, l = -1;
		for(int i = n+m; i >= 0; --i)
			if(!dcmp(c[i], 0) && c[i] > 0) e = i; //Bland规则
		if(e == -1) break;
		double det = INF;
		for(int i = n+m; i >= 0; --i) {
			if(!B[i] || A[i][e] < 0 || dcmp(A[i][e], 0)) continue;
			if(b[i]/A[i][e] < det || dcmp(det, b[i]/A[i][e])) l = i, det = b[i]/A[i][e];
		}
		PIVOT(l, e, 0);
	}
	if(dcmp(b[0], 0)) { //找到一个辅助线性规划的最优解，也就是原问题的一个基本可行解
		memcpy(c, c0, sizeof c);
		for(int i = 0; i <= n+m; ++i)
			if(B[i] && !dcmp(c[i], 0)) {
				v = v+c[i]*b[i];
				for(int j = 0; j <= n+m; ++j) {
					if(j == i) continue;
					c[j] = c[j]-c[i]*A[i][j];
				}
				c[i] = 0;
			}
		c[0] = 0;
		for(int i = 0; i <= n+m; ++i)
			if(B[i]) A[i][0] = 0;
		return 1;
	}
	else return 0; //原线性规划无解 
}