线性规划之单纯形算法

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

PIVOT

PIVOT\((N,B,A,b,c,v,l,e)\)
\ \ \ \ let \(\hat{A}\) be a new \(m\times n\) matrix
\ \ \ \ \(\hat{b}_{e}=b_l/a_{le}\)
\ \ \ \ for each \(j\in N-\{e\}\)
\ \ \ \ \ \ \ \ \(\hat{a}_{ej}=a_{lj}/a_{le}\)
\ \ \ \ \(\hat{a}_{el}=1/a_{le}\)
\ \ \ \ for each \(i\in B-\{l\}\)
\ \ \ \ \ \ \ \ \(\hat{b}_i=b_i-a_{ie}\hat{b}_e\)
\ \ \ \ \ \ \ \ for each \(j\in N-\{e\}\)
\ \ \ \ \ \ \ \ \ \ \ \ \(\hat{a}_{ij}=a_{ij}-a_{ie}\hat{a}_{ej}\)
\ \ \ \ \ \ \ \ \(\hat{a}_{il}=-a_{ie}\hat{a}_{el}\)
\ \ \ \ \(\hat{v}=v+c_e\hat{b}_e\)
\ \ \ \ for each \(j\in N-\{e\}\)
\ \ \ \ \ \ \ \ \(\hat{c}_j=c_j-c_e\hat{a}_{ej}\)
\ \ \ \ \(\hat{c}_l=-c_e\hat{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 \(j\in N\) has \(c_j>0\)
\ \ \ \ \ \ \ \ choose an index \(e\in N\) for which \(c_e>0\)
\ \ \ \ \ \ \ \ for each index \(i\in B\)
\ \ \ \ \ \ \ \ \ \ \ \ if \(a_{ie}>0\)
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(\Delta_i=b_i/a_{ie}\)
\ \ \ \ \ \ \ \ else
\ \ \ \ \ \ \ \ \ \ \ \ \(\Delta_i=\infty\)
\ \ \ \ \ \ \ \ choose an index \(l\in 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 \(i\in 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_k\geqslant 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; //原线性规划无解 
}