编程手养成计划 Day16

学习线性规划，我参考了几位大佬的文章，这里先把链接放出来方便查阅：

一、什么是线性规划

（一）线性规划解决什么问题？

线性规划（Linear Programming 简记LP）研究的是在线性不等式或等式的限制条件下，使得某一线性目标取得最大（或最小）的问题。

（二）线性规划的特点

1. 优点

有统一算法，任何线性规划问题都能求解，解决多变量最优决策的方法。

2. 缺点

对于数据的准确性要求高，只能对线性的问题进行规划约束，而且计算量大，有由线性规划演变的非线性规划法等等后续的方法弥补，但是计算量增加许多。

（三）线性规划的数学表达式

$$max{x}_1+2x_2(1.1)$$
$$s.t.\begin{array}{r}{x}_1+x_2\le 3(1.2)\\ x_2\le 2(1.3)\\ x_1\ge 0(1.4)\\ x_2\ge 0(1.5)\end{array}$$
其中，max表示求解最大值，s.t.(subject to)表示受到约束。我们不难发现上面的所有式子均是线性的，而线性规划的目的即在满足s.t.的前提下寻找max或min。

（四）线性规划的几何意义

如图：根据约束条件，我们可以得到图中蓝色区域，即可行域。红色线表示所求式子的一种取值，即在红色线上取$(x_1,x_2)$使得$(x_1+2x_2)$的值不变，而我们的目标是求得什么时候$(x_1+2x_2)$的值最大

二、线性规划算法

（一）图解法

1.算法原理

我们可以令$x_1+2x_2=b$，那么就有$x_2=-\frac{1}{2}{x}_1+\frac{1}{2}b$，根据我们初中就学过的$y=kx+b$可以知道，红色线与$x_2$轴交点就是$\frac{1}{2}b$的取值，所以当我们向上平移红色线到超出可行域时，我们就得到了$(x_1+2x_2)$的最大值。

2.算法优点

图解法可以直观、简便的得到最优解。

3.算法缺点

一旦变量超过两个，很难画出图像。

4.启示

最优解总是出现在可行域构成的多面体的顶点处。
(关于这一点的证明，请大家自行移步文章最上面的链接)

（二）枚举法

根据图解法，我们得到了一条结论：最优解总是出现在可行域构成的多面体的顶点处。
那么我们只需要算出所有顶点，带入到式子中并找到最优解即可。

1.顶点的计算

（1）计算顶点的条件：
在计算顶点之前，我希望先明确标准形式下约束矩阵需要满足的条件：$A\in m\times n$，需满足$m\le n$，并且矩阵A行满秩。以保证约束条件不会冗余，同时方程有解。

（2）计算顶点的步骤：

• 从 $n$个变量中任意抽取其中$m$ 个，然后将剩余的 $n-m$个变量赋值为0。
• 抽取得到的$m$个变量构成$m\times m$的方程组，求解这个方程组即可获得顶点。
• 判断这个顶点是否满足约束$x\ge 0$，若不满足则舍弃掉，若满足则保存该顶点。

2.多胞体表示定理

（1）重要性质：
一旦多胞体的顶点确定下来了，则该多胞体就被确定下来了。
（意思就是顶点确定后，这个图形就会确定）
（2）多胞体表示定理：
多胞体表示定理：多胞体$P$中的任意一点可以表示为所有$P$的顶点的凸组合。 设多胞体 P = { x ∈ R n ∣ A x = b , x ≥ 0 } P=\{{x∈R^n|Ax=b, x≥0\}} P={x∈Rn∣Ax=b,x≥0}是非空有界的，设$P$的顶点集合为 V = { v 1 , v 2 , . . . , v K } V=\{{v_1,v_2,...,v_K\}} V={v1​,v2​,...,vK​}，上述定理可以表示为 $\forall x\in P$ 有如下表达式成立：
$$x=\sum _{i=1}^K{\theta }_i{v}_i$$
其中${\sum }_{i=1}^K{\theta }_i=1,{\theta }_i\ge 0,\forall i=1,...,K$。
（3）线性规划基本定理：
对于标准形式的线性规划问题，如果该问题存在有界的最优解，那么至少有一个最优解在顶点上。 对于 $\forall x\in P$ 有如下表达式成立：
$$\begin{array}{rlr}{c}^{T}x& =c^T\sum _{i=1}^K{\theta }_i{\upsilon }_i& (1)\\ & =\sum _{i=1}^K{\theta }_i(c^T{\upsilon }_i)& (2)\\ & \ge \underset{i\in V}{\min }\left\{c^T{\upsilon }_i\right\}\left(\sum _{i=1}^K{\theta }_i\right)& (3)\\ & =c^T{\upsilon }_{\min }& (4)\end{array}$$
表达式(1)成立是因为多面体表示定理，表达式(2)成立就是简单的一个顺序变化，表达式(3)成立是因为是最小值，表达式(4)成立是因为 ${\sum }_{i=1}^K{\theta }_i=1$。

3. 无解与解无界

（1）无解：
即约束条件中有冲突。
（2）解无界：
不说那么多理论，直接上图：

我们发现，将红色线向下平移后，仍然有点在可行域内，所以此时对红色线求最小值则为负无穷。同理向上平移后则得到最大值为正无穷。
此时对于 $\forall x\in P$ 有如下表达式成立：
$$\begin{array}{rlr}{c}^{T}x& =c^T(\sum _{i=1}^K{\theta }_i{\upsilon }_i+d)& (1)\\ & \ge \underset{i\in V}{\min }\left\{c^T{\upsilon }_i\right\}\left(\sum _{i=1}^K{\theta }_i\right)& (3)\\ & =c^T{\upsilon }_{\min }& (4)\end{array}$$
我们发现，这个式子与上面的公式极其相似，区别在于：当 $c^{T}d\le 0$ 时，会导致无界，因此该情况不是我们线性规划基本定理讨论的范围。
即，我们默认 $c^{T}d\ge 0$。

4. 枚举法的缺点：

算法复杂度高，上述步骤需要循环$C_m^n$次才能求得到最优解。

（三）单纯型法

1. 大致框架

• Step 1：从一个初始的基可行解出发；
• Step 2：检查是否是最优解（最优解条件），若是最优解则停止，否则进入下一步；
• Step 3：从当前顶点移动到更好的顶点，然后跳转到 Step 2；

2. 详细算法及证明

看了一下，因为计算过程很是复杂，不好写出来，因为这篇主要还是备赛笔记嘛，以实现功能为主，所以，在这里放出来学习链接吧，如果以后有机会再把详细过程补上：

详细算法请看：
【十分钟弄懂单纯形法】
具体的理论证明请看：
【线性规划(三)】单纯型法（上）
另外这位博主也详细介绍了求解线性规划的其他方法，大家可以去拜读。

三、代码实现

接下来我们将围绕一道例题，学习使用编程语言解决DP问题的方法。题目如下：
$$\begin{gathered} z=\min -x_0+4x_1 \\ s.t.=\{\begin{array}{rl}-3x_0+x_1& \le 6,\\ -x_0-2x_1& \ge -4,\\ x_1& \ge -3\end{array} \end{gathered}$$

（一）Python

1. scipy.optimize.linprog()函数

（1）函数输入

下面是scipy.optimize.linprog()的完整参数列表：

def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,
            bounds=None, method='interior-point', callback=None,
            options=None, x0=None)

具体参数信息：

参数	意义	数据类型
c	目标函数的决策变量对应的系数向量(行列向量都可以，下同)	一维数组
A_ub	不等式约束组成的决策变量系数矩阵	二维数组
b_ub	由A_ub对应不等式顺序的阈值向量	一维数组
A_eq	等式约束组成的决策变量系数矩阵	二维数组
b_eq	由A_ub对应等式顺序的阈值向量	一维数组
bounds	表示决策变量x连续的定义域的n×2维矩阵，None表示无穷	(min, max)序列对
method	调用的求解方法	字符串
callback	选择的回调函数	scipy.optimize.OptimizeResult
options	求解器选择的字典	字典
x0	初始假设的决策变量向量，若可行linprog会对方法优化	一维数组

注：

• bounds：

1. None：使用None表示没有界限，默认情况下，界限为（0，None）（所有决策变量均为非负数）
2. 如果提供一个元组(min, max)，则最小值和最大值将用作所有决策变量的界限。

• method

1. 算法，{'simplex', 'revised simplex', 'interior-point', 'highs', 'highs-ds', 'highs-ipm'}以上算法可选。
   在 osgeo.cn 的scipy文档中对method的方法给出以下解释：

方法 'highs-ds' 是C++高性能双重修订单工实现(HSOL)的包装器 [13], [14]. 
方法 'highs-ipm' 是C++实现的 i n前面-p oint m 方法 [13]; 它的特点是交叉例程，因此它与单纯形求解器一样精确。
方法 'highs' 自动在两者之间进行选择。
有关涉及以下内容的新代码 linprog ，我们建议显式选择这三个方法值中的一个，而不是 'interior-point' (默认)(内点法)， 'revised simplex' ，以及 'simplex' (旧版)。

• options

1. maxiter：整数，要执行的最大迭代次数

2. disp：布尔值，设置为True以打印收敛消息，默认值：False

3. autoscale：布尔值，设置为True以自动执行平衡，如果约束中的数值分开几个数量级，请考虑使用此选项，默认值：False

4. presolve：布尔值，设置为False可禁用自动预解析，默认值：True

5. rr：布尔值，设置为False可禁用自动移除冗余，默认值：True

（2）函数输出

参数	意义
x	在满足约束的情况下将目标函数最小化的决策变量的值
fun	目标函数的最佳值$(c^{T}x)$
slack	不等式约束的松弛值（名义上为正值）$b_{ub}-A_{ub}x$
con	等式约束的残差（名义上为零）$b_{eq}-Aeqx$
success	当算法成功找到最佳解决方案时为 True
status	表示算法退出状态的整数
nit	在所有阶段中执行的迭代总数
message	算法退出状态的字符串描述符

关于status的值:

• 0 : 优化成功终止
• 1 : 达到了迭代限制
• 2 : 问题似乎不可行
• 3 : 问题似乎是不收敛
• 4 : 遇到数值困难

（3）程序实例

"""
这个问题没有以'linprog'所接受的形式呈现。这一点可以
通过将 "大于 "不等式约束转换为 "小于 "不等式约束，可以很容易地解决这个问题。
不等式约束转换成 "小于 "不等式约束，就可以轻松解决。
将两边都乘以:math:`-1`。还请注意，最后一个
约束实际上是简单的约束：数学：`-3 \leq x_1 \leq \infty`。
最后，由于在 :math:`x_0`上没有约束，我们必须明确地
指定边界 :math:`-\infty \leq x_0 \leq \infty`，因为
默认情况下，变量是非负的。在收集系数后
成数组和图元后，这个问题的输入是
"""

>>> c = [-1, 4]
>>> A = [[-3, 1], [1, 2]]
>>> b = [6, 4]
>>> x0_bounds = (None, None)
>>> x1_bounds = (-3, None)
>>> from scipy.optimize import linprog
>>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])

# 请注意，`linprog'的默认方法是'interior-point'，其本质是的性质是近似的。

>>> print(res)
     con: array([], dtype=float64)
     fun: -21.99999984082494 # may vary
 message: 'Optimization terminated successfully.'
     nit: 6 # may vary
   slack: array([3.89999997e+01, 8.46872439e-08] # may vary
  status: 0
 success: True
       x: array([ 9.99999989, -2.99999999]) # may vary

# 如果你需要更高的精确度，请尝试 'revised simplex'.

>>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds], method='revised simplex')
>>> print(res)

2. pulp库

（1）首先，我们要定义一个规划问题：

prob = LpProblem("Problem", LpMinimize)

LpProblem 是定义问题的构造函数。
　　"Problem"是用户定义的问题名（用于输出信息）。
　　参数 sense 用来指定求最小值/最大值问题，可选参数值：LpMinimize、LpMaximize 。
（2）定义决策变量

x0 = LpVariable("x0")
x1 = LpVariable("x1", -3)

LpVariable 是定义决策变量的函数，下面是完整参数：
LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None)
　　‘x1’ 是用户定义的变量名。
　　参数 lowBound、upBound 用来设定决策变量的下界、上界；可以不定义下界/上界，默认的下界/上界是负无穷/正无穷。
　　参数 cat 用来设定变量类型，可选参数值：‘Continuous’表示连续变量（默认值）、’ Integer ' 表示离散变量（用于整数规划问题）、’ Binary ’ 表示0/1变量（用于0/1规划问题）。
（3）添加目标函数

prob += 4 * x1 - x0

添加目标函数使用 “问题名 += 目标函数式” 格式。
（4）添加约束条件

prob += (-3 * x0 + x1 <= 6)
prob += (-3 * x0 - 2 * x1 >= -4)

添加约束条件使用 “问题名 += 约束条件表达式” 格式。
　　约束条件可以是等式约束或不等式约束，不等式约束可以是 小于等于 或 大于等于，分别使用关键字">=“、”<=“和”=="。
（5）求解

prob.solve()		# 求解
print("Status:", LpStatus[prob.status])	# 打印求解状态
for a in prob.variables():	# 输出每个变量的值
  print(a.name, "=", a.varValue)
print("F(x) = ", pulp.value(prob.objective))	# 输出最优解的目标函数值

求解状态输出 Optimal表示最优解。
solve() 是求解函数。PuLP默认采用 CBC 求解器来求解优化问题，也可以调用其它的优化器来求解，如：GLPK，COIN CLP/CBC，CPLEX，和GUROBI，但需要另外安装。
（6）完整代码

from pulp import *
# 定义规划问题
prob = LpProblem("Problem", LpMinimize)
# 添加决策变量
x0 = LpVariable("x0")
x1 = LpVariable("x1", -3)
# 添加目标函数
prob += -x0 + 4 * x1
# 添加约束条件
prob += (-3 * x0 + x1 <= 6)
prob += (-x0 - 2 * x1 >= -4)
# 求解
prob.solve()
print("Status:", LpStatus[prob.status])
for a in prob.variables():	
   print(a.name, "=", a.varValue)
print("F(x) = ", pulp.value(prob.objective))

关于更多的Pulp库求解线性规划问题的方法，请看：
Python数模笔记-PuLP库（1）线性规划入门
Python数模笔记-PuLP库（2）线性规划进阶
Python数模笔记-PuLP库（3）线性规划实例

（二）Matlab

先放出官方文档：
linprog
Matlab中求解线性规划的命令为：
[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub)
式中：

• x返回决策向量的取值;
• fval返回目标函数的最优值;
• f为价值向量;
• A和b对应线性不等式约束;
• Aeq和beq对应线性等式约束；
• lb和ub分别对应决策向量的下界向量和上界向量；

PS:
(1)若约束条件中间不齐全，则可以用空向量[]来替代缺失部分，如：
[x,fval] = linprog(f,A,b,[],beq,[],ub)
(2)若约束条件只有前面条件，则可以直接省略参数，如：
[x,fval] = linprog(f,A,b)

完整代码如下：