提示：本文参考了scipy的linprog源码，对源码感兴趣的小伙伴可以直接去读源码，注释真的是非常详尽了，比代码都长。

1. 补充问题

上一节中的代码在运行时还有很多细节没有处理，这里补充两个比较重要的情况：

存在等式约束
如果有等式约束，那么就没法通过添加松弛变量直接给出初始可行解，需要用大M法或者两阶段法求解。计算机求解一般使用二阶段法，即首先给等式约束条件添加人工变量，使得问题有一个初始可行解。二阶段法是第一阶段最小化人工变量的和使其等于0，而大M法则是给人工变量添加一个非常大的系数M使其等于0。
注意人工变量和松弛变量/剩余变量不一样。松弛变量/剩余变量是用于将不等式变为等式（标准形），而人工变量则是在等式约束中额外添加一个变量，这个变量最后一定要等于0才行，否则就是没有初始可行解，不能进入第二阶段。
在第一阶段，每个不等式约束对应一个松弛变量，每个等式约束对应一个人工变量。另外，b<0的不等式约束也需要添加人工变量，因为对应的松弛变量等于b不能作为基变量。第一阶段需要添加一个额外的目标函数：min⁡Σxav\min \Sigma x_{av}minΣxav，xavx_{av}xav指的是所有的人工变量。
RHS <0
如果right hand side（也就是b）小于0，则这一行也无法给出初始可行解。此时的一般做法是左右都乘以-1，添加松弛变量变为等式约束，然后添加人工变量得到初始可行解。
没有有限最优解（目标函数无界）
对应的情况是：目标函数中有某个系数大于0，对应变量在约束条件中的系数全都≤0，问题结束，没有最优解。

首先假设问题除了不等式约束，还有等式约束：
求解问题为：
min cxcxcx
s.t.
Aubx≤bubA_{ub} x ≤ b_{ub}Aubx≤bub
Aeqx=beqA_{eq} x = b_{eq}Aeqx=beq

2. 基于scipy的python算法实现

下面代码的核心步骤是构造两阶段的标准单纯形矩阵，具体求解函数可参考运筹系列1，这里我们直接调用scipy的_solve_simplex函数。

首先定义一些变量：
T：用于求解的矩阵。第二阶段比第一阶段要多一行。
basis：基向量的下标集合，通过solve_simplex函数改变
maxiter：最大迭代次数
tol：求解精度
OptimizeResult：格式输出函数
solve_simplex：单纯形法求解函数。

完整的求解步骤如下：

from scipy.optimize._linprog_simplex import _solve_simplex # 求解标准单纯形矩阵的函数
from scipy.optimize import OptimizeResult # 格式化输出结果的函数
def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, maxiter=1000, tol=1.0E-12):cc = np.asarray(c)f0 = 0n = len(c)mub = len(b_ub)meq = len(b_eq)m = mub + meq  # 变量总数n_slack = mub  # 不等号添加松弛变量n_artificial = meq + np.count_nonzero(b_ub < 0)  # 等号/没有初始可行解的不等号添加人工变量Aub_rows, Aub_cols = A_ub.shapeAeq_rows, Aeq_cols = A_eq.shapeslcount = 0avcount = 0status = 0messages = {0: "Optimization terminated successfully.",1: "Iteration limit reached.",2: "Optimization failed. Unable to find a feasible"" starting point.",3: "Optimization failed. The problem appears to be unbounded.",4: "Optimization failed. Singular matrix encountered."}# 建立第一阶段单纯型表T。T = np.zeros([m + 2, n + n_slack + n_artificial + 1])  # 添加两行和一列T[-2, :n] = cc  # T倒数第二行用来记系数T[-2, -1] = f0b = T[:-2, -1]  # T倒数第一列用来记bif meq > 0:  # 给T赋值T[:meq, :n] = A_eqb[:meq] = b_eqif mub > 0:T[meq:meq + mub, :n] = A_ubb[meq:meq + mub] = b_ubnp.fill_diagonal(T[meq:m, n:n + n_slack], 1) # 剩余变量对应系数1basis = np.zeros(m, dtype=int)r_artificial = np.zeros(n_artificial, dtype=int)for i in range(m):if i < meq or b[i] < 0:basis[i] = n + n_slack + avcountr_artificial[avcount] = iavcount += 1if b[i] < 0:b[i] *= -1T[i, :-1] *= -1T[i, basis[i]] = 1T[-1, basis[i]] = 1 # T的最后一行是新增加的约束条件else:basis[i] = n + slcountslcount += 1for r in r_artificial:T[-1, :] = T[-1, :] - T[r, :] # 将人工变量的1全部变为0，T转为标准形。至此T构造完毕。# 第一阶段求解nit1, status = _solve_simplex(T, n, basis)# 如果第一阶段的目标函数为0，则删除人工变量，进入第二阶段if abs(T[-1, -1]) <= tol:T = T[:-1, :]T = np.delete(T, np.s_[n + n_slack:n + n_slack + n_artificial], 1)else:status = 2if status != 0:message = messages[status]return OptimizeResult(x=np.nan, fun=-T[-1, -1], nit=nit1, status=status, message=message, success=False)nit2, status = _solve_simplex(T, n, basis, maxiter=maxiter - nit1, phase=2, tol=tol, nit0=nit1)if status != 0:message = messages[status]return OptimizeResult(x=np.nan, fun=-T[-1, -1], nit=nit2, status=status, message=message, success=False)solution = np.zeros(n + n_slack + n_artificial)solution[basis[:m]] = T[:m, -1]x = solution[:n]slack = solution[n:n + n_slack]obj = -T[-1, -1]return OptimizeResult(x=x, fun=obj, nit=int(nit2), status=status, slack=slack, message=messages[status],success=(status == 0))

3. 例子说明

下面是例子：

import numpy as np
c=np.array([2,3,-5])
A_ub=np.array([[-2,5,-1]])
B_ub=np.array([-10])
A_eq=np.array([[1,1,1]])
B_eq=np.array([7])
res=linprog(-c,A_ub,B_ub,A_eq,B_eq)
print(res)

对应 max⁡2x1+3x2−5x3\max 2x_1+3x_2-5x_3max2x1+3x2−5x3
s.t.
−2x1+5x2−x3≤−10-2x_1+5x_2-x_3 ≤ -10−2x1+5x2−x3≤−10
x1+x2+x3=7x_1+x_2+x_3 =7x1+x2+x3=7

输出为

     fun: -14.571428571428571message: 'Optimization terminated successfully.'nit: 2slack: array([0.])status: 0success: Truex: array([6.42857143, 0.57142857, 0.        ])

下面具体说明一下中间过程：
首先，将max转为min问题，需要将c取负数，最后的结果相应的也取负数即可。
其次，转变为标准形，添加了2个slack variable，添加了1个artificial variable（-10小于0，因此需要一个人工变量）。
第一阶段的问题需要添加一个额外的约束，如下图：

此时矩阵T为：

调用标准单纯形法进行求解，最后的T如下：

其中T[-1,-1]代表第一阶段的最优值为0，可以继续进行第二阶段的求解。删除最后一行，以及标号为4、5的两列，求解后如下表：

因为只有2个约束条件，因此基变量也只有2个，剩下的两个变量=0。

4. 代码跟踪

我们如果用运筹系列1的代码进行跟踪，第一阶段：
初始化：

做一个简单的转换，将最后添加的两个人工变量变成基变量：

第一次pivoting：

第二次pivoting：

接下来进入第二阶段，首先剔除人工变量所在的列，然后将最后一行替换为原来的c:

进行消元，生成初始可行解：

然后将d，s送入求解器，直接就得到最优解：

参考文档

[1] Dantzig, George B., Linear programming and extensions. Rand Corporation Research Study Princeton Univ. Press, Princeton, NJ, 1963
[2] Hillier, S.H. and Lieberman, G.J. (1995), “Introduction to Mathematical Programming”, McGraw-Hill, Chapter 4.
[3] Bland, Robert G. New finite pivoting rules for the simplex method. Mathematics of Operations Research (2), 1977: pp. 103-107.
[4] Andersen, Erling D., and Knud D. Andersen. “The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm.” High performance optimization. Springer US, 2000. 197-232.
[5] Andersen, Erling D. “Finding all linearly dependent rows in large-scale linear programming.” Optimization Methods and Software 6.3 (1995): 219-227.
[6] Freund, Robert M. “Primal-Dual Interior-Point Methods for Linear Programming based on Newton’s Method.” Unpublished Course Notes, March 2004. Available 2/25/2017 at https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf
[7] Fourer, Robert. “Solving Linear Programs by Interior-Point Methods.” Unpublished Course Notes, August 26, 2005. Available 2/25/2017 at http://www.4er.org/CourseNotes/Book%20B/B-III.pdf
[8] Andersen, Erling D., and Knud D. Andersen. “Presolving in linear programming.” Mathematical Programming 71.2 (1995): 221-245.
[9] Bertsimas, Dimitris, and J. Tsitsiklis. “Introduction to linear programming.” Athena Scientific 1 (1997): 997.
[10] Andersen, Erling D., et al. Implementation of interior point methods for large scale linear programming. HEC/Universite de Geneve, 1996.

运筹系列2：线性规划两阶段法python代码相关推荐

两阶段法-Python实现
Python单纯形法-两阶段法单纯形法简介 Python代码 1.主函数 2.定义Simplex()大类 3. 将初值全部放入大矩阵T中 4.最优性检验函数 5.迭代函数: 6. 去人工变量 7.两 ...
线性规划问题及单纯形法-两阶段法
两阶段法两阶段法:用计算机处理数据时,只能用很大的数代替M,可能造成计算机上的错误,这个M无法确定,故采用两阶段法,和大M法是一致的. 第一阶段:在原线性规划问题中加入人工变量,使其目标函数值为人工 ...
利用两阶段法通过寻找基可行解求线性规划问题的最优解
算法介绍: java代码实现: package sy1; //标准化系数矩阵并加上人工变量 public class BzhAndJrg {public double A[][]; //原矩阵的系数 ...
两阶段法求解线性规划求解
用两阶段法求解 min f=2x1-x2+x3 s.t x1+2x2- x3=1 2x1+ x2+ x3=5 x1- x2+2x3=4 xi>=0,i=1,2,3 ...
运筹学两阶段法编程c语言,运筹学上机实验 - 单纯形方法的两阶段法
理论部分不解释了, 就是粘个实验课的代码,按照书上的算法写的,仅仅是把课本上的样例过了,有bug可能难免,欢迎指出. Sample 1. $$ \left\{ \begin{aligned} min ...
运筹学上机实验 - 单纯形方法的两阶段法
理论部分不解释了, 就是粘个实验课的代码,按照书上的算法写的,仅仅是把课本上的样例过了,有bug可能难免,欢迎指出. Sample 1. $$ \left\{ \begin{aligned} min ...
单纯形法；大M法；两阶段法
目录线性规划问题的标准形式 1.单纯形法 1.1定义 1.2思路 1.3计算步骤 1.4Python求解 1.5调用scipy包求解 2.大M法 2.1思路 3.两阶段法思路线性规划问题的标准形 ...
线性规划两阶段求解方法
想快速掌握线性规划两阶段求解方法,强烈推荐博文<三言两语讲清楚线性规划单纯形方法>. 百度百科给了下面一个例子,感觉其解法不容易看明白原理,换一种解释方法,应该很容易看明白两阶段法的原理. ...
凸优化笔记4（两阶段法）
前言下面先简要介绍两阶段法,在通过例题说明具体流程.重点看例题,有些问题在本刊其他文章中不予赘述. 一.两阶段法介绍大M法与两阶段法都是在原问题缺少初始可行基的情况下利用引人人工变量构造人工基,以 ...

运筹系列2：线性规划两阶段法python代码