Datawhale 202208 GitModel | 线性规划 上

0. 绪论

运筹优化的基本模型思路通常分为两部，需要解决"确定最佳方案""求出最佳的配比"此类的"最值"问题。对于优化类的问题,如何建模与如何算法求解也是最重要的部分。两个核心内容: 优化建模,下降算法。

我们涉及的优化模型由三个基本组成:

• 决策变量 $x_1,x_2,\cdots ,x_n$——影响目标的自变量$,$

• 目标函数 $f(x_1,x_2,\cdots ,x_n)$——确定最佳决策的指标依据$,$

• 约束条件 $c_i(x_1,x_2,\cdots ,x_n)⩽(=,⩾)0$——现实中的限制.

通常尽可能多的采用向量语言: 如利用向量 $\mathbf{x}=(x_1,x_2,\cdots ,x_n{)}^T$ 表示决策变量$,$ 或利用矩阵来表达如下的约束

$$A\mathbf{x}⩽\mathbf{b}\iff \{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n⩽b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n⩽b_2\\ ⋮\\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n⩽b_m\end{array},A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right].$$

此外$,$ 我们用 $R,R^n$ 代表实数域以及 $n$ 维实数 (向量) 空间$,f:R^n\to R$ 代表 $n$ 元实值函数.

用规范化的语言$,$ 我们现实中遇到的所有优化问题都可描述为如下的形式:

$$\begin{array}{rl}& \min f(\mathbf{x})\\ s.t.& \{\begin{array}{l}{c}_i(\mathbf{x})⩽0,i=1,2,\cdots ,m\\ c_j(\mathbf{x})=0,j=1,2,\cdots ,k\end{array}\end{array}\text{\hspace{1em}(0.1)(0.2)}$$

其中$,\mathbf{x}\in R^n,f:R^n\to R,$ 并且我们不要求目标函数 $f$ 以及约束函数 $c_i,c_j$ 有什么多余的性质如可微 (这一点对刚入门的初学者而言是比较反直觉的)$,$ 对于求极大值的问题$,$ 可以等价写为

$$\max f(\mathbf{x})\iff \min -f(\mathbf{x}).$$

优化问题 $0.1$ 中$,$ 约束 $0.2$ 构成的子集称为可行域 $D,$ $D$ 中的点是符合 $0.1$ 要求的点$,$ 称为可行解$,$ 令 $f(\mathbf{x})$ 最小的可行解 ${\mathbf{x}}_0$ 称为最优解$,$ 我们一般用

a r g m i n { f ( x ) ∣ x ∈ D } \mathrm{argmin} \{f(\boldsymbol{x})\mid \boldsymbol{x}\in D\} argmin{f(x)∣x∈D}

记为 $0.1$ 的最优解. 当 $D=R^n$ 时 $0.1$ 称为无约束优化问题$,$ 相应地 $D⊊R^n$ 时称之为约束优化问题.

我们涉及的下降算法基本思想是生成一列递减的点 { x n } n = 1 ∞ \{\boldsymbol{x}_n\}_{n=1}^{\infty} {xn​}n=1∞​ 来逼近 $0.1$ 的最优解$,$ 一般由以下几步组成:

1. 确定初始解 ${\mathbf{x}}_1$.

2. 确定迭代方式 ${\mathbf{x}}_{n+1}=g({\mathbf{x}}_n)$ 使得 $f({\mathbf{x}}_{n+1})⩽f({\mathbf{x}}_n),n=1,2,\cdots$.

3. 确定终止准则 (计算机只能迭代有限点$,$ 我们希望在有限步能尽可能接近最优解)——当 ${\mathbf{x}}_n$ 与 $f({\mathbf{x}}_n)$ 满足什么条件时$,$ 认为 ${\mathbf{x}}_n$ 最小值的近似点$,$ 终止算法.

迭代方式是下降算法最核心的一步$,$ 它决定算法生成的点列能否收敛得到最优解$,$ 迭代步长决定算法收敛的速度 (即算法运行的时间)$,$ 这在求解实际的优化问题过程中是非常关键的内容.

本课程介绍的是最简单的一类优化问题: 目标函数与约束条件都是线性的$,$ 称为线性规划问题(模型)$,$ 在日常生活中蕴含着非常多地线性规划问题$,$ 如生产策略$,$ 工作指派等等$,$ 所以建立线性规划模型事实上是在现实生活中非常常见的$,$ 而由于性质的优越$,$ 线性规划问题的下降算法设计也是非常简单的——即我们接下来要学习的单纯形法与分支定界法.

1. 线性规划（上）

🎯作业二：请你仿照上面代码，尝试求解以下模型？

$$\min W_{order}+W_{trans}$$
$$s.t.\{\begin{array}{l}{x}_3=z_{3,9+}+z_{3,10+}+z_{3,10-},\\ x_4=z_{4,10+}+z_{4,10-}+z_{4,11-},\\ x_i⩽s_i,i=3,4,\\ x_3+x_4=780,\\ z_{3,9+}⩾0,z_{4,11-}⩾0,\\ z_{i,10+}⩾0,z_{i,10-}⩾0,i=3,4,\end{array}$$
其中
$$W_{order}=155x_3+160x_4,$$
$$W_{trans}=73.3z_{3,9+}+101.1(z_{3,10+}+z_{3,10-})+93.1(z_{4,10+}+z_{4,10-})+78.9z_{4,11-}.$$

	0	1	2	3	4	5	6	7
模型变量	$x_3$	$x_4$	$z_{3,9^{+}}$	$z_{3,10^{-}}$	$z_{3,10^{+}}$	$z_{4,10^{-}}$	$z_{4,10^{+}}$	$z_{4,11^{-}}$
数组索引	0	1	2	3	4	5	6	7

# 代码如下

## 导入线性求解器，线性和整数求解的API都在其中
from ortools.linear_solver import pywraplp


## 声明求解器
solver = pywraplp.Solver.CreateSolver('GLOP')


## 创建变量 
infinity = solver.infinity()

x2 = solver.IntVar(0.0, infinity, 'z3_9ad')
x7 = solver.IntVar(0.0, infinity, 'z4_11re')

x4 = solver.IntVar(0.0, infinity, 'z3_10ad')
x6 = solver.IntVar(0.0, infinity, 'z4_10ad')

x3 = solver.IntVar(0.0, infinity, 'z3_10re')
x5 = solver.IntVar(0.0, infinity, 'z4_10re')

x0 = solver.IntVar(0.0, 1000, 'x3')
x1 = solver.IntVar(0.0, 2000, 'x4')

print('`NumVar`说明该变量是连续型变量，前面两个参数分别是该变量取值的上下界.如果没有上界的话设置为`solver.infinity()`. ')
print('变量个数 =', solver.NumVariables())


## 添加约束
solver.Add(x0 + x1 == 780)
solver.Add(x2 + x4 + x3 == x0)
solver.Add(x6 + x5 + x7 == x1)

print('约束个数 =', solver.NumConstraints())


## 定义目标函数
solver.Minimize( (155 * x0 +160 * x1) + (73.3 * x2 + 101.1 * (x3 + x4) + 93.1*(x5 + x6) + 78.9*x7) )


## 调用求解器
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('\n最优值 =', solver.Objective().Value())
    print('x3', x0.solution_value())
    print('x4', x1.solution_value())
    print('z3_9ad', x2.solution_value())
    print('z3_10re', x3.solution_value())
    print('z3_10ad', x4.solution_value())
    print('z4_10re', x5.solution_value())
    print('z4_10ad', x6.solution_value())
    print('z4_11re', x7.solution_value())
else:
    print('The problem does not have an optimal solution.')

没有上界的话设置为`solver.infinity()`. 
变量个数 = 8
约束个数 = 3

最优值 = 178074.0
x3 780.0
x4 0.0
z3_9ad 780.0
z3_10re 0.0
z3_10ad 0.0
z4_10re 0.0
z4_10ad 0.0
z4_11re 0.0

2. 线性规划（下）

2.3 线性规划模型的求解方式: 单纯形法

在讲解线性规划求解方法之前$,$ 我们引入必要的知识准备.

2.3.1 必要的知识准备一: 标准形式及化为标准形式的方法

对于一个线性规划$,$ 由于其约束是由线性等式以及线性不等式构成的$,$ 而在线性等式中寻找可行点 (比如求线性方程组) 比在线性不等式中寻找可行点的难度要小得多$,$ 对于任意一条线性不等式约束（不妨设为第$i$个约束）

$$a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n⩽b_i,$$

我们引入新的变量 $s_i⩾0$ 使得原约束变为等式约束
$$a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n+s_i=b_i,\text{\hspace{1em}(1.1)}$$
这样的变量 $s_i$ 我们称为松弛变量$,$ 同样的对于 $⩾$ 型的不等式约束的也可引入松弛变量解决问题
$$a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n⩾b_i\to a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n-s_i=b_i.$$

将一个线性规划问题中所有不等式约束等式化后我们得到线性规划的标准形式

$$\begin{array}{rl}& \min f(\mathbf{x})=c_1{x}_1+c_2{x}_2+\cdots +c_n{x}_n,\\ s.t.& \{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1,\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2,\\ ⋮\\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_m,\\ x_i⩾0,i=1,2,\cdots ,n\end{array}\end{array}\text{\hspace{1em}(1.2)}$$

其中$,$ $a_{ij}\in R,b_i⩾0,\mathbf{x}=(x_1,x_2,\cdots ,x_n)\in R^n$ 是 $n$ 个连续的决策变量 (含松弛变量) $,$ 具有 $m$ 条等式约束$,$ 由于松弛变量的存在一般假设 $n⩾m$ 且无法相互线性组合$,$ 我们利用矩阵也可以将其写成如下紧凑形式

$$minf(\mathbf{x})={\mathbf{c}}^T\mathbf{x},\text{\hspace{1em}(1.3)}$$
$$s.t.\{\begin{array}{l}A\mathbf{x}=\mathbf{b}\\ \mathbf{x}⩾0,\end{array}$$
其中 $A$ 为系数矩阵$,$ 从约束无法相互线性组合可知 $A$ 的秩必为 $m$.

除开引入松弛变量$,$ 对于一些常见的现实情况$,$ 我们需要一些小操作将其化为标准形式:

1. 当问题需要求解的最优解是最大值时$,$ 事实上等价于求解如下问题
   $$\max c_1{x}_1+c_2{x}_2+\cdots +c_n{x}_n\iff \min -(c_1{x}_1+c_2{x}_2+\cdots +c_n{x}_n).$$
2. 当问题的决策变量 $x_i$ 是非正的$,$ 或任意值时$,$ 我们做如下操作非负化决策变量
   $$x_i⩽0\to \text{引入}{x}_i^{\prime }=-x_i,$$
   $$x_i\in R^n\to \text{引入}{x}_i^{\prime },x_i^{\prime \prime }⩾0,x_i=x_i^{\prime }-x_i^{\prime \prime }.$$
3. 当约束的常数 $b_i$ 非正时$,$ 该约束等价于
   $$a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n=b_i\iff -(a_{i1}{x}_1+a_{i2}{x}_2+\cdots +a_{in}{x}_n)=-b_i.$$

集合 D = { a i 1 x 1 + a i 2 x 2 + ⋯ + a i n x n = b i ∣ x i ⩾ 0 , i = 1 , 2 , ⋯   , m } D=\{a_{i1}x_1+a_{i2}x_2+\cdots+a_{in}x_n= b_i\mid x_i\geqslant 0,i=1,2,\cdots,m\} D={ai1​x1​+ai2​x2​+⋯+ain​xn​=bi​∣xi​⩾0,i=1,2,⋯,m} 称为问题的可行域$,$ $D$ 中每个点都是问题的一个解$,$ 称为可行解$,$ 而使得问题达到最小值的可行解称为最优解$,$ 对于一个线性规划问题$,$ 其不一定有最优解$,$ 因为可行域可以无界$,$ 当然$,$ 若存在最优解$,$ 也不一定唯一$,$ 下述二图解释了线性规划无最优解或最优解不唯一的情况 (若我们构建方式不当$,$ 也有可能出现 $D=\varnothing$ 即不存在可行解的情况$,$ 为下文命题阐述严谨起见提出这一点).

2.3.2 从二维图解法出发引入单纯形法的思路

从优化理论出发$,$ 线性规划本质上是一个约束优化问题$,$ 无法简单地用无约束优化的一二阶条件求导数零点进行判断$,$ 幸运的是由于目标函数以及约束都是线性的$,$ 我们有更加简洁且实用的方法——单纯形法.

从高中学的线性规划知识我们知道对于变量数为 $2$ 的线性规划问题可以用图解法迅速解决$,$ 我们从图解法出发构建高维线性规划的解决方式. 事实上$,$ 由图解法我们可以观察得到两个基本事实:

1. 线性规划可行域是多边形.
2. 大多数存在唯一解的线性规划问题, 其解会在多边形可行域顶点取得.

那么对于高维的线性规划问题, 这两个特征会不会被保留下来呢$?$ 答案是会的, 我们需要借用凸集来推广高维的多边形.

定义 1.1.1 [凸集] 设 $C$ 是 $R^n (n\geqslant 2) $ 上的非空子集, 若对于任意 $x_1,x_2\in C$ 有 $\alpha x_1+(1-\alpha )x_2,\alpha \in [0,1]$ 则 $C$ 称为凸集. 若凸集 $C$ 上一点 $x$ 满足不存在 $x_1,x_2\in C$ 使得 $x=\alpha x_1+(1-\alpha) x_2(\alpha \in (0,1)) $ 则 $x$ 称为凸集的顶点.

从几何意义上来说$,$ 凸集势必会包含任意集合中两点的连线$,$ 而顶点只能是凸集中线段的端点 (等价于$,$ 只能是边界"线段"的交点) . 那么$,$ 可以仅由线段分解的多边形便是高维空间中最简单的非平凡凸集. 对于线性规划与凸集$,$ 显然从几何想象可以得到以下基本事实:

命题 1.1.2 线性规划若存在可行解$,$ 则问题可行域必定为凸多边形.

而目标函数是线性的其实保证了$,$ 可行域上的函数图像是一个倾斜的 $n$ 维凸多边形面$,$ 那么其最小值必定会在多边形某个顶点取到$!$

命题 1.1.3 线性规划若存在最优解$,$ 则最优解在凸多边形可行域的顶点集中取得.

那么$,$ 对于线性规划的最优解搜索范围$,$ 我们缩小到了有限的顶点集中$,$ 自然求解会方便的多$,$ 接下来的问题有二:

1. 如何计算线性约束可行域的顶点$?$

2. 枚举计算顶点的函数值可以找到最优值$,$ 但随着维数升高计算量太大了$,$ 有没有什么更好的最优顶点筛选方式$?$

2.3.3 求解顶点–基可行解

从本小节开始$,$ 我们正式开始展开单纯形法的讨论$,$ 由于数学要素较多$,$ 推荐在学习之前积累以下的前置知识:

1. 向量的线性相关性.
2. 向量组与矩阵的秩$,$ 极大线性无关组.
3. 线性方程组解的情况以及对应条件.

经过代数学的学习$,$ 我们可以敏感地感知到: 标准形式的线性规划可行域其实就是某个线性方程组的非负解集$,$ 那么自然会诞生这样一个问题: 可行域对应的线性方程组
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ ⋮\\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_m\\ x_i⩾0,i=1,2,\cdots ,n\end{array}\iff A\mathbf{x}=\mathbf{b}$$ 中什么样的解才是可行域对应的凸多边形顶点呢$?$ 利用我们高中以及大学的解析几何知识联想一下$,$ 二维空间中$,$ 二元一次方程组的唯一解对应着两条直线的交点$;$ 三维空间中$,$ 三元一次方程组的唯一解对应着三个平面的交点$,$ 那么若 $m$ 元一次方程有唯一解$,$ 则该解必定对应着 $n$ 维空间中 $m$ 个"超平面"的交点.

而代数知识告诉我们$,$ 什么样的线性方程组有唯一解$?$ 系数矩阵 $A$ 可逆的呗.可惜 $A$ 是一个长方阵 (列数 $n$ 比行数 $m$ 大) $,$ 仅仅是行满秩的$,$ 并无法保证其具有唯一解$,$ 那么我们将其表达成向量的线性组合形式$,$
$$A=\left[\begin{array}{cccc}{\mathbf{\alpha }}_1& {\mathbf{\alpha }}_2& \cdots & {\mathbf{\alpha }}_n\end{array}\right],{\mathbf{\alpha }}_i=\left[\begin{array}{c}{a}_{1i}\\ a_{2i}\\ ⋮\\ a_{ni}\end{array}\right],$$
$$A\mathbf{x}=\mathbf{b}\iff x_1{\mathbf{\alpha }}_1+x_2{\mathbf{\alpha }}_2+\cdots +x_n{\mathbf{\alpha }}_n=\mathbf{b},\text{\hspace{1em}(1.4)}$$
$A$ 的秩为 $m$ 代表着: 向量组 $\left[\begin{array}{cccc}{\mathbf{\alpha }}_1& {\mathbf{\alpha }}_2& \cdots & {\mathbf{\alpha }}_n\end{array}\right]$ 的秩为 $m,$ 即其中有 $m$ 个线性无关的向量$,$ 我们不妨将其挑出来: $B=\left[\begin{array}{cccc}{\mathbf{a}}_{i_1}& {\mathbf{a}}_{i_2}& \cdots & {\mathbf{a}}_{i_m}\end{array}\right],$ 这也代表着 $A$ 中对应的 $m\times m$ 子矩阵 $B$ 的秩是 $m$ (满秩) 从而是可逆的$,$ 事实上$,$ 我们只需要取 $x_k=0,k\ne i_j,j=1,2,\cdots ,m$ 即方程 $4$ 变为
$$\begin{array}{cccc}\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ ⋮\\ a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_m\end{array}& ⟶& \{\begin{array}{l}{a}_{11}0+\cdots +a_{1i_1}{x}_{i_1}+\cdots +a_{1i_1}{x}_{i_2}+\cdots +a_{1i_m}{x}_{i_m}+\cdots +a_{1n}0=b_1\\ a_{21}0+\cdots +a_{2i_1}{x}_{i_1}+\cdots +a_{2i_2}{x}_{i_2}+\cdots +a_{2i_m}{x}_{i_m}+\cdots +a_{2n}0=b_2\\ ⋮\\ a_{m1}0+\cdots +a_{m{i}_1}{x}_{i_1}+\cdots +a_{m{i}_2}{x}_{i_2}+\cdots +a_{m{i}_m}{x}_{i_m}+\cdots +a_{mn}0=b_m\end{array}& \\ \Updownarrow & & \Updownarrow & \\ x_1{\mathbf{\alpha }}_1+x_2{\mathbf{\alpha }}_2+\cdots +x_n{\mathbf{\alpha }}_n=\mathbf{b}& ⟶& 0{\mathbf{\alpha }}_1+\cdots +x_{i_1}{\mathbf{\alpha }}_{i_1}+\cdots +x_{i_2}{\mathbf{\alpha }}_{i_2}+\cdots +\cdots +x_{i_m}{\mathbf{\alpha }}_{i_m}+\cdots +0{\mathbf{\alpha }}_n=\mathbf{b}& \\ \Updownarrow & & \Updownarrow & \\ A\mathbf{x}=\mathbf{b}& ⟶& B{\mathbf{x}}^{\prime }=\mathbf{b}\Rightarrow {\mathbf{x}}^{\prime }=B^{-1}\mathbf{b}(B\text{行列满秩故可逆})& \end{array}$$
其中
$$B=\left[\begin{array}{cccc}{\mathbf{a}}_{i_1}& {\mathbf{a}}_{i_2}& \cdots & {\mathbf{a}}_{i_m}\end{array}\right]=\left[\begin{array}{cccc}{a}_{1i_1}& a_{1i_2}& \cdots & a_{1i_m}\\ a_{2i_1}& a_{2i_2}& \cdots & a_{2i_m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m{i}_1}& a_{m{i}_2}& \cdots & a_{m{i}_m}\end{array}\right],{\mathbf{x}}^{\prime }=\left[\begin{array}{c}{x}_{i1}\\ x_{i2}\\ ⋮\\ x_{im}\end{array}\right].$$

那么$(0,\cdots ,0,x_{i_1},\cdots ,x_{i_m},0,\cdots ,0)$ 便是 $A\mathbf{x}=\mathbf{b}$ 的一组解$,$ 并且也是子方程组 $B{\mathbf{x}}^{\prime }=\mathbf{b}$ 的唯一解$,$ 从而是约束 $A\mathbf{x}=\mathbf{b}$ 对应凸多边形的一个顶点$,$ 只要再保证 $x_{i_j}⩾0,j=1,2,\cdots ,m,$ 那么这样的解便是线性规划可行域凸多边形的一个顶点$!$

$A$ 中每组 $m$ 列线性无关的向量 $B=\left[\begin{array}{cccc}{\mathbf{a}}_{i_1}& {\mathbf{a}}_{i_2}& \cdots & {\mathbf{a}}_{i_m}\end{array}\right]$ (亦即$,$ 满秩的 $m\times m$ 子矩阵 $B$) 对应的子方程 $B{\mathbf{x}}^{\prime }=\mathbf{b}$ 唯一解$,$ 只要非负$,$ 便是可行域多边形的顶点之一$,$ 这是因为 $m$ 条等式线性约束在 $n$ 维代表了一个 $m$ 维的超"平面"$,$ 再加上变量非负 $x_j⩾0$ 的限制$,$ 这样的超平面又被无形的"坐标墙"框成了多边形$,$ 自然就必须有一些分量为 $0!$ 以三维空间为例子$,$ 顶点都是二维平面约束在坐标轴上的交点.

我们将 $A$ 中列满秩的$m\times m$ 子矩阵 $B$ 称为基$,$ 其对应的线性方程组解 $(0,\cdots ,0,x_{i_1},\cdots ,x_{i_m},0,\cdots ,0)$ 若有 $x_{i_j}⩾0,j=1,2,\cdots ,m$ 则称为基可行解. 基可行解最为重要的意义便是等价于可行域的顶点.

可以严谨地证明:

命题 1.1.3 线性规划问题 $2,3$ 中可行域中的点 $\mathbf{x}$ 是基可行解的必要条件是至少有 $n-m$ 个分量为 $0$.

命题 1.1.4 线性规划问题 $2,3$ 中凸多边形可行域的顶点集与线性规划的基可行解一一对应.

这边解决了我们在 1.3.3 提出的第一个问题:

凸多边形的顶点可直接计算基可行解求得.

参考内容

Github课程链接