运筹学(一):线性规划(LP)

前言

参考清华大学经管类《运筹学》(第五版)

线性规划(LP)

基本概念

决策变量(Decision Variable):即规划问题中需要确定的能用数量表示的量

目标函数(Objective Function):关于决策变量的函数,也是决策者优化的目标

约束条件(Constraint):即决策变量需要满足的限制条件,通常表达为关于决策变量的等式或不等式

线性规划的标准型

标准型的形式为:
::: align-center
\begin{aligned}&\max z=2x_1+x_2\\&s.t.\left\{\begin{matrix}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\\\cdots\\a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n=b_m\\x_1,x_2,\cdots,x_n\ge 0\end{matrix}\right .\end{aligned}

:::
其中b_i\ge 0,上述标准型对应的矩阵形式为
::: align-center
\begin{aligned}&\max z=CX\\&s.t.\left\{\begin{matrix}AX=b\\X\ge 0\end{matrix}\right .\end{aligned}

:::
其中C为“目标函数系数”,b为“约束条件右边项”,矩阵A为约束条件的“系数矩阵”或“工艺矩阵”(秩为m,因为已经化为标准型,即假设m个约束条件均是线性无关的)
注意到约束条件均为等式,如果为不等式:

  • 若为“\le”型不等式,在等式左侧加一个松弛变量(Slack),可将其变为等式
  • 若为“\ge”型不等式,在等式左侧减一个剩余变量(Surplus),可将其变为等式

线性规划的图解法

图解法一般只用于解决两变量的线性规划问题,通用步骤:

  1. 画二维直角坐标系,非负约束构成坐标系的第一象限
  2. 画出每条约束所对应的区域(对不等式约束,首先画出等式先,再判明约束方向),并确定线性规划问题的可行域(即各条约束所对应区域的交集)
  3. 根据目标优化方向平移目标函数等值线,直到不能再平移为止,确定线性规划问题对应的最优点
  4. 根据最优点满足的等式构建联立方程组,从而求解出最优方案,并计算对应的目标函数值

解的情形:

  1. 存在多个解:等值线的斜率刚好和可行域的一个边界平行
  2. 可行域为空集:无解(Infeasible)
  3. 无有界最优解(Unbounded)

解的性质

几个概念

可行解与最优解

可行解(Feasible Soultion)和最优解(Optimal):满足线性规划模型约束条件的解称为可行解,所有可行解构成的集合称为可行域(Feasible Region),在可行域中,能使目标函数达到最大的解称为最优解(Optimium)

基、基向量与基变量

因为工艺矩阵A的秩为m,从列向量的角度讲至少存在m列彼此是线性独立的,故在A中存在一个m\times m阶子矩阵B,其秩为m(即B为满秩矩阵),称B为一个基阵,简称为基(Base),基阵中的每一个列向量P_{j}(j=1,2,\cdots, m)称为一个基向量,与之相对应的变量x_j称为基变量,其他列向量称为非基向量,除基变量之外的其他变量称为非基变量

基解、基可行解与可行基

将变量X按行分块为基变量和非基变量
::: align-center
X=\begin{bmatrix}X_B\\X_N\end{bmatrix}

:::
考虑约束条件AX=b,按照上述符号定义,方程组等价于
::: align-center
(B,N)\begin{bmatrix}X_B\\X_N\end{bmatrix}=b

:::
BX_B+NX_N=b
考虑到B是满秩矩阵它是可逆的,有
::: align-center
X_B=B^{-1}(b-NX_N)

:::
对应于基B,如果令X_N=0(即非基变量为全0),可以得到一个解
::: align-center
X=\begin{bmatrix}B^{-1}b\\0\end{bmatrix}

:::
我们称该解是对应于基B的基解,当然B^{-1}b不一定所有分量都为非负取值,只有当B^{-1}b\ge 0时,解X才是原线性规划问题的一个可行解,我们把满足非负约束的基解称为基可行解,对应于基可行解的基阵称为一个可行基

凸集与凸组合

Dn维欧氏空间的一点集,若D中任意两点X_1,X_2的连线上的所有点也属于D,则称D为一个凸集
因为两点之间的连线可以表示为
::: align-center
\alpha X_1+(1-\alpha)X_2,0<\alpha<1

:::
因此,凸集用数学语言描述为:对集合D中任意两点X_1,X_2,如果对任意\alpha\in (0,1),均有\alpha X_1+(1-\alpha)X_2\in D,则称集合D为一个凸集

X_i,i=1,2,\cdots,k,是n维欧式空间的k个点,若存在一组数\mu_i,i=1,2,\cdots,k,满足\mu_i\in[0,1],而且\mu_1+\mu_2+\cdots+\mu_k=1,那么
::: align-center
X=\mu_1X_1+\mu_2X_2+\cdots+\mu_kX_k

:::
是点X_1,X_2,\cdots,X_k凸组合

相应地,对于一个凸集D中的点X,如果不存在两个不同的点X_i\in DX_2\in D,使得X成为这两个连线上的一个点,则我们称点X为一个顶点

几个基本定理(不再证明)

  1. 如果一个线性规划问题的可行域非空,则其可行域是凸集
  2. 如果线性规划问题的可行域有界,则其最优值必可在某个顶点处有界
  3. 线性规划问题的基可行解X刚好对应于可行域上的某个顶点

求解线性规划问题的单纯性法

单纯型法(Simplex)搜索可行解的基本思路是:先找出一个初始的基可行解,然后判断该基可行解是否最优;如果否,则转换到相邻的能进一步改善目标函数值的基可行解,直到找到最优解为止

单纯性法的引入

依笔者理解,在一般情况(即约束条件的变量个数n严格大于m)下,我们有n-m个自由变量(线性代数概念上的自由变量),这些变量在运筹学上是通过松弛变量/剩余变量两个概念引入的,在线性代数中,我们通过基础解系中这些自由变量的组合可以获得代数方程的通解,而在运筹学上,我们通过调整这些自由变量的“松弛度”和“剩余度”则可能会让目标函数获得更好的取值结果。单纯型法就是这样工作的,我们为了能够更为方便的发现“如何调整才能让目标函数更好,且好的更快”考虑了如何确定最佳的换入和换出变量,通过这种“贪心”的方法,在有解的情况下我们最终会收敛到一个全局的Optimal。在几何上,当我们确定了某个基可行解,即图形上的某个顶点(见上方的定理)后,我们在做的是考虑“向此顶点周围的哪条边“移动”(或者说,转换,因为在后面接触到DP,也就是动态规划时会更好理解),到达另一个近邻顶点后能让目标函数更好,且好的更彻底”的问题。

现在我们用一个例子
::: align-center
\begin{aligned} &\max z=2x_1+3x_2\\ &s.t.\left\{\begin{matrix} x_1&+2x_2&+x_3&&&=8\\ 4x_1&&&+x_4&&=16\\ &4x_2&&&+x_5&=12\\ x_1,x_2\ge0 \end{matrix}\right . \end{aligned}

:::
引入单纯型法

正如在线性代数中获得基础解系一样,我们需要先获得一个基可行解,在求解步骤上与线性代数获得基础解系没有区别,在这个例子中,基阵的秩为3,由秩-零化度定理(Rank-Nullity Theorem)得到自由度为5-3=2,基变量应有3个,剩余2个为非基变量,获得基可行解的最简便方法为选择基变量为x_3,x_4,x_5,因为将线性方程组重排后,关于x_3,x_4,x_5的线性方程组系数矩阵已经是行最简形:
::: align-center
\left\{\begin{aligned} x_3&&&+x_1&+2x_2&=8\\ &x_4&&+4x_1&&=16\\ &&x_5&&+4x_2&=12 \end{aligned}\right .\Rightarrow\begin{pmatrix}1&0&0&1&2\\0&1&0&4&0\\0&0&1&0&4\end{pmatrix}\begin{pmatrix}x_3\\x_4\\x_5\\x_1\\x_2\end{pmatrix}=\begin{pmatrix}8\\16\\12\end{pmatrix}

:::
x_1=x_2=0,很容易得到基可行解为(x_1,x_2,x_3,x_4,x_5)^T=(0,0,8,16,12)^T
但在一般情况下,我们可能要考虑选择合适的可行解(在这个例子中,考虑任选3个)后通过高斯若尔当消元得到基可行解,或者后面的人工变量法

之后我们检查在此解的基础上,目标函数是否还可以最大化,目前的基变量固定为x_3,x_4,x_5,非基变量x_1,x_2有调整的空间(即可以增减),查看此时的目标函数z=2x_1+3x_2,发现不论x_1还是x_2增大都可以进一步改善目标函数,而显然增加x_2最大化的更快一些,此时我们确定x_2应被作为换入变量(Entering Baic Variable,又叫入基变量)(注意:换入是相对于基变量而言的,后面的换出也是),而对应的我们自然会想到随着x_2的变化目前的其他基变量会有什么变化,是否会有“临界情况”,再次查看约束条件我们发现最后一条的非负性约束还没有被用到(对于一般的线性规划问题我们总是假设有非负性约束条件,如果有特殊情况,按最后总结部分的表处理),于是我们考虑整理上方的式子为:
::: align-center
\left\{\begin{aligned} x_3&&&+x_1&+2x_2&=8\\ &x_4&&+4x_1&&=16\\ &&x_5&&+4x_2&=12 \end{aligned}\right .\Rightarrow\left\{\begin{aligned}&x_3=8-x_1-2x_2\\&x_4=16-4x_1\\&x_5=12-4x_2\end{aligned}\right .

:::
因为我们暂时不考虑x_1,保持原样,令x_1=0,之后套入对x_3,x_4,x_5非负性约束:
::: align-center
\left\{\begin{aligned}&x_3=8-2x_2\ge0\\&x_4=16\ge0\text{(恒成立)}\\&x_5=12-4x_2\ge0\end{aligned}\right .

:::
解方程组得到x_2可以增大,但不能超过3,x_2=3是临界情况,而此时x_5=0,这意味着x_5此时由基变量转换为了非基变量,那么我们自然就可以把x_5作为换出变量(Leaving Basic Variable,又称出基变量)

于是我们整理式子为:
::: align-center
\left\{\begin{aligned} x_3&&+2x_2&+x_1&&=8\\ &x_4&&+4x_1&&=16\\ &&+4x_2&&+x_5&=12 \end{aligned}\right .

:::
使用高斯若尔当消元,完毕后再写为式子:
::: align-center
\left\{\begin{aligned}&x_3=2-x_1+\frac{1}{2}x_5\\&x_4=16-4x_1\\&x_2=3-\frac{1}{4}x_5\end{aligned}\right .

:::
将其带入目标函数,将目标函数重写为非基变量(现在基变量已经变为x_3,x_4,x_2,目标函数中不应再出现x_2,直接带入最后一条等式即可)的组合,得到:
::: align-center
z=9+2x_1-\frac{3}{4}x_5

:::
同时令x_1=0,x_2=3,此时基可行解变为(x_1,x_2,x_3,x_4,x_5)^T=(0,3,2,16,0)^T,我们便完成了一次转换,在几何上这意味着我们发现由x_2=0的那个顶点转换到x_5=0的那个顶点可以进一步使目标函数最大化,于是我们向着x_1=0的那条“边”(实际上不是边而是超平面,因为维度已经达到5维)走到了新的顶点

再次观察目标函数,我们发现依然可以通过增加x_1使得目标函数进一步最大化,重复刚才的过程,容易得x_1是这次的换入变量,换出变量由非负性条件解得为x_3,迭代后的基可行解为(2,3,0,8,0)^T,又发现目标函数中x_5的系数变为正,于是考虑将x_5作为换入变量,换出变量解得为x_4,再经过一次迭代得到基可行解为(4,2,0,0,4)^T,目标函数变为:
::: align-center
z=14-1.5x_3-0.125x_4

:::
此时不论增加x_3还是x_4均不能再让目标函数再最大化(不要考虑减少x_3,x_4,因为我们得策略是始终“前进”,不允许出现“回头”的情况,一般情况下在标准形下我们只“前进”不“后退”,也不要让x_3,x_4为负,因为这会违反非负性条件),那么此时的基可行解就是最优解

单纯型表

读者可以发现使用上面的原型方法计算是十分啰嗦的,而且不便于抽象为计算机算法,针对这种问题,我们需要从纯代数的角度出发解决

单纯形表的布局方式如下表所示:
单纯形表.png

其中有几点需要注意:

  1. c_jC_Bc_j行填目标函数中所有决策变量的系数(部分教材中会把这些系数取相反数,对应的,下方的检验系数的意义会发生变化,这样做的好处是在检验系数对应行的b的对应列预先写0,在最后求得最优解时,这个entry的值就是目标函数的最值),而C_B列只填基变量
  2. 检验系数:对于基变量其检验系数总为0,对于非基变量x_{m+1}-x_{n},检验系数的计算公式(以x_j为例)为c_j-\sum_{i=1}^mc_ia_{ij}
  3. \theta用于最小比例检验(Minimum Ratio Test),计算公式为\min(\frac{b_i}{a_{ik}}|a_{ik}>0)=\frac{b_r}{a_{rk}}

计算步骤:

  1. 根据数学模型确定初始可行基和初始基可行解,建立初始单纯形表
  2. 计算各非基变量x_j的检验系数,如果所有检验系数为非正(某些教材中c_j行填相反数,则这里为非负,下同),则已得到规划问题的最优解,可终止计算;否则,转入下一步
  3. 在所有检验系数大于零的非基变量中,如果某个非基变量x_k对应的系数向量P_k\le 0,则此问题无有界最优解(因为违反了非负性约束),可终止计算;否则,转入下一步
  4. 在所有检验系数大于零的非基变量中,选择检验系数最大的非基变量x_k作为入基变量,按\theta规则计算,选最小的\theta对应X_Bl列的基变量作为出基变量,入基变量和出基变量确定系数矩阵A中的唯一entry:a_{rk}
  5. a_{rk}为中心进行换基迭代,把x_k所对应的列向量P_k=\begin{bmatrix}a_{1k}\\a_{2k}\\\vdots\\a_{rk}\\\vdots\\a_{mk}\end{bmatrix}变换为\begin{bmatrix}0\\0\\\vdots\\1\\\vdots\\0\end{bmatrix}
    X_B列中的x_r换为x_k,得到新的单纯形表,重复2-5直到终止

几种特殊情形

  1. 无穷多最优解的情形:当最终单纯形表中某个非基变量的检验系数为零有无穷多最优解,因为检验系数行代表着目标函数的增量,而如果这个决策变量又可以作为入基变量,即可以换入基可行性解,或几何上存在可以向该方向转换的“边”时,我们知道向这个方向转换不会给目标函数造成任何增量,则这条“边”上的所有点都是最优解
  2. 退化解的情形:如果通过两种基迭代得到的基变量不同,但是最优解完全相同,我们称此时的最优解为一个“退化解”,单纯形表中体现为某个基变量在b列的值为0,在几何上这代表沿着几条“边”走会走到同一个顶点,也即几个约束条件的方程交于同一个顶点,特别地,在二维(两个决策变量)的情况下如图所示
    退化解.png
  3. 极小化目标函数的问题:只需要改变对于检验变量的评判标准即可,在当前的定义(检验系数即是目标函数的增量)下,最优解在检验系数列全非负时取得,后同

人工变量法

此类方法适用于解决初始基可行解无法方便得出的情况

大M法

以问题:
::: align-center
\begin{aligned}&\min w=-3x_1+x_2+x_3\\ &s.t.\left\{\begin{aligned}&x_1-2x_2+x_3\le 11\\&-4x_1+x_2+2x_3\ge 3\\&-2x_1+x_3=1\\&x_1,x_2,x_3,x_4,x_5\ge 0\end{aligned}\right .\end{aligned}

:::
为例,将其变为标准形式:
::: align-center
\begin{aligned}&\max z=3x_1-x_2-x_3\\ &s.t.\left\{\begin{aligned}&x_1-2x_2+x_3+x_4=11\\&-4x_1+x_2+2x_3-x_5=3\\&-2x_1+x_3=1\\&x_1,x_2,x_3,x_4,x_5\ge 0\end{aligned}\right .\end{aligned}

:::
其中x_4是引入的松弛变量,x_5是引入的剩余变量,以上问题的稀疏矩阵很显然不包含任何三阶单位阵,无法直观地给出初始的基可行解。为了人为地“凑出”一个单位阵,我们可以考虑在约束条件中再引入两个非负的人工变量(Aritificial Dummy Variable)(记为x_6,x_7),即
::: align-center
\left\{\begin{aligned}&x_1-2x_2+x_3+x_4=11\\&-4x_1+x_2+2x_3-x_5+x_6=3\\&-2x_1+x_3+x_7=1\\&x_1,x_2,x_3,x_4,x_5\ge 0\end{aligned}\right .
:::
此时观察这三个方程组的最右侧的最后一个变量,容易发现此时(x_4,x_6,x_7)对应的列向量刚好构成一个单位阵(这里没有按列对齐写方程组,但是读者应能看出),然而上述加了人工变量的方程组的可行域显然和没加时是不同的,除非x_6,x_7取值0,也就是说,我们需要进行一定的等价变换,保证x_6,x_7取值为0,一种方法是对目标函数z=3x_1-x_2-x_3进行适当的调整,调整为:
::: align-center
\hat{z}=3x_1-x_2-x_3-Mx_6-Mx_7

:::
其中M是一个取值为无穷大的正数(注意到人工变量前的系数都是-M而不是-M_1-M_2,后者是没有必要的,因为效果相同,前者反而在后面的检验系数中更易判断大小),-Mx_6,-Mx_7可以理解为原目标函数的罚项,容易知道当原始问题的可行域非空时原始最大化问题的最优解一定有x_6=x_7=0(否则目标函数将无穷小)

接下来即可正常进行单纯形法,检验系数将变为含M的线性函数,使用极限的思想判断检验系数的大小即可(例如,3M-1很显然应比M-1更大,而M-1M-2极限相等,在标准的最大化问题单纯形法解题步骤中,很显然应该选前者),在迭代过程中系数带M的人工变量会很快变为出基变量,此时这些人工变量全部取0,它们的使命也就完成了

两阶段法

相信读者已经知道人工变量的使命是为方便地找到一个初始的基可行解,那也可以把前面的大M法分为两个阶段求解

  1. 构造一个辅助的线性规划模型,其目标是人工变量之和(在后续替换为其他非基变量的组合),优化方向是最小化,约束条件是引入人工变量后的等式约束,如果此阶段最终的单纯形表中所有人工变量均出基,也即取值为0,则第一阶段的最优解就是原始问题的一个基可行解,进入第二阶段计算,如果任意人工变量的取值为正,那么原始问题无解
  2. 在获得最小化人工变量(全0)的松弛解后,即可在第一阶段的单纯形表中删去人工变量,将目标函数替换为原始问题的目标函数,继续利用单纯形法迭代求得最优解

相比较而言,两阶段法中因为不需要引入M因而计算更为便利

额外的Notes和总结

在大M法中如果我们经过多次换基迭代已经达到了最优性条件,即所有非基变量检验系数都非正,但扔出现某个人工变量依然在基变量中且取值为正的情况,那么此时最优解对应的目标函数会是一个含M的线性函数(其系数为负),即目标函数的取值是负无穷大,这意味着原始线性规划问题的可行域是空集,因此大M法提供了一种判断线性规划问题可行域是否是空集的方法

附上约束条件的处理方法和单纯性法的解题步骤
约束的处理方法.png

单纯形法的解题步骤.png
游客

全部评论 (0)

暂无评论,快来抢沙发吧~