10.4 线性规划

一、线性规划

线性规划（Linear programming）是在线性代数的基础上再结合两个新的思想：不等式（inequalities） 和最小化（minimization）。这个问题仍然是从矩阵方程 $A\mathbf{x}=\mathbf{b}$ 开始，但是只考虑非负解向量，我们要求 $\mathbf{x}\ge 0$（意味着 $\mathbf{x}$ 的分量都不能为负）。一般来说，未知数要比方程的个数多，即矩阵的 $n>m$，因此，如果方程 $A\mathbf{x}=\mathbf{b}$ 存在有 $\mathbf{x}\ge 0$ 的解，那么可能会有很多个这样的解。线性规划是研究如何寻找解 ${\mathbf{x}}^{\ast }\ge 0$ 使得下面的成本最小：

成本为 $c_1{x}_1+c_2{x}_2+\cdots c_n{x}_n$，寻找线性方程组 $A\mathbf{x}=\mathbf{b}$ 的非负解向量 ${\mathbf{x}}^{\ast }$，使得成本最小。

因此，线性规划问题是由一个矩阵 $A$ 和两个向量 $\mathbf{b}$ 和 $\mathbf{c}$ 来定义的：

1. $m\times n$ 的矩阵 $A$ 有 $n>m$：如 $A=\left[\begin{array}{ccc}1& 1& 2\end{array}\right]$（一个方程，三个未知数）
2. $m$ 个方程组 $A\mathbf{x}=\mathbf{b}$ 中，$\mathbf{b}$ 有 $m$ 个分量：如 $\mathbf{b}=\left[\begin{array}{c}4\end{array}\right]$
3. 成本系数向量（cost vector） 有 $n$ 个分量：如 $\mathbf{c}=\left[\begin{array}{ccc}5& 3& 8\end{array}\right]$.

那么这个问题就是最小化成本 $\mathbf{c}\cdot \mathbf{x}$，其中约束条件为 $A\mathbf{x}=\mathbf{b}$ 和 $\mathbf{x}\ge 0：$$$最小化5x_1+3x_2+8x_3，约束条件是x_1+x_2+2x_3=4和x_1,x_2,x_3\ge 0$$我们这里直接跳到了问题部分，没有解释问题的来源。线性规划实际上是数学在管理学中最重要的应用，开发相关快速求解的算法和程序代码是非常有意义的。一般来说，寻找 ${\mathbf{x}}^{\ast }$ 要比直接求解 $A\mathbf{x}=\mathbf{b}$ 更难，因为对于解向量要满足额外的约束条件：${\mathbf{x}}^{\ast }\ge 0$ 和成本 ${\mathbf{c}}^T{\mathbf{x}}^{\ast }$ 最小。求解该问题后，会解释问题的背景，并介绍著名的单纯形法和内点法。
首先观察约束条件：$A\mathbf{x}=\mathbf{b}$ 和 $\mathbf{x}\ge 0$. 方程 $x_1+x_2+2x_3=4$ 表示的是三维空间中的一个平面，加上非负性的条件 $x_1\ge 0,x_2\ge 0,x_3\ge 0$ 会将这个平面限制到了一个三角形区域，所以解向量 ${\mathbf{x}}^{\ast }$ 一定位于 Figure 10.5 中的三角形 $PQR$ 上。

在这个三角形 $PQR$ 的内部，$\mathbf{x}$ 所有的分量都是正的；在这个三角形的边上，有一个分量为零；而在它三个顶点 $P,Q$ 和 $R$ 处，有两个分量为零。最优解（optimal solution）${\mathbf{x}}^{\ast }$ 一定是这些顶点中的一个！ 下面说明为什么。
这个三角形区域包含所有满足约束条件 $A\mathbf{x}=\mathbf{b}$ 和 $\mathbf{x}\ge 0$ 的向量 $\mathbf{x}$，这些 $\mathbf{x}$ 称为可行点（feasible points）， 而这个三角形区域称为可行域（feasible set）. 然后，我们要在可行域中寻找能够使得成本 $\mathbf{c}\cdot \mathbf{x}$ 最小化的可行点：

$$在三角形PQR中寻找使得成本5x_1+3x_2+8x_3最小的点{\mathbf{x}}^{\ast }$$

使成本为零的向量位于平面 $5x_1+3x_2+8x_3=0$ 上，但是这个平面三角形并不相交，所以我们无法在 $\mathbf{x}$ 满足约束条件的情况下，达到零成本。因此需要逐渐增加成本 $C$，直至平面 $5x_1+3x_2+8x_3=C$ 与三角形相交。随着 $C$ 的增加，我们可以得到一族向着三角形区域移动的互相平行的平面。
首先与三角形区域相交的平面 $5x_1+3x_2+8x_3=C$ 有最小的 $C$. 交点就是最优解 ${\mathbf{x}}^{\ast }$. 而这个交点一定是顶点 $P,Q$ 或 $R$ 中的一个，一个平移的平面在经过三角形顶点之前不可能会先经过其内点！因此我们需要比较在每个顶点处的成本 $5x_1+3x_2+8x_3$：

$$P=(4,0,0)对应成本20Q=(0,4,0)对应成本12R=(0,0,2)对应成本16$$

顶点 $Q$ 对应的成本最小，所以 ${\mathbf{x}}^{\ast }=(0,4,0)$ 就是这个线性规划问题的最优解。
如果成本系数向量 $\mathbf{c}$ 发生改变，那么这一族互相平行的平面就会发生倾斜。如果改变比较小，$Q$ 仍然是最优解。对于成本 $\mathbf{c}\cdot \mathbf{x}=5x_1+4x_2+7x_3$，最优解 ${\mathbf{x}}^{\ast }$ 就会移动到 $R=(0,0,2)$，此时最小成本是 $7\cdot 2=14$.

注1： 有些线性规划问题是寻求利润最大化而不是成本最小化，这在数学原理上几乎是相同的。此时选择的是一个较大的 $C$，然后随着 $C$ 的减小，这一族平行的平面向着原点移动（不再是远离了）。第一个交点仍然是顶点。

注2： 约束条件 $A\mathbf{x}=\mathbf{b}$ 和 $\mathbf{x}\ge 0$ 可能无法同时满足。如方程 $x_1+x_2+x_3=-1$ 就没有满足 $\mathbf{x}\ge 0$ 的解，此时可行域是空集。

注3： 可行域也有可能的无解的。如果约束条件是 $x_1+x_2-2x_3=4$，此时一个较大的正向量 $(100,100,98)$ 就是一个可行点，而更大的向量 $(1000,1000,998)$ 也是一个可行点。此时平面 $A\mathbf{x}=\mathbf{b}$ 将不会再被截断为一个三角形了。其中两个顶点 $P$ 和 $Q$ 仍然可能是最优解 ${\mathbf{x}}^{\ast }$，但是顶点 $R$ 移动到了无穷远处。

注4： 对于一个无解的可行域，最小成本可能是 $-\infty$（负无穷大）。假设成本函数是 $-x_1-x_2+x_3$，则向量 $(100,100,98)$ 对应的成本 $C=-102$；向量 $(1000,1000,98)$ 对应的成本 $C=-1002$. 此时我们会获得报酬，而不是支付成本。在实际应用中这种情况不可能发生，但是理论上是可能的，因为 $A,\mathbf{b}$ 和 $\mathbf{c}$ 可能会产生无法预料的三角形区域和成本函数。

二、原问题和对偶问题

下面举一个例子来解释 $A,\mathbf{b},\mathbf{c}$. 未知数 $x_1,x_2,x_3$ 分别表示一个博士、一个学生和一台机器的工作时间，每小时的成本分别为 $$5$，$$3$ 和 $$8$，工作的时间不能为负：$x_1\ge 0,x_2\ge 0,x_3\ge 0$. 博士和学生都是每小时完成一道题目，而机器每小时可以完成两道题。现在有四道题目需要完成，原则上它们可以分工合作完成：$x_1+x_2+2x_3=4$.
现在的问题是用最小的成本 ${\mathbf{c}}^T\mathbf{x}$ 完成这四道题目。
如果这三者同时工作，那么这项工作需要一个小时完成：$x_1=x_2=x_3=1$，成本是 $5+3+8=16$ 美元。但是我们可以发现博士和学生的解题速度是一样的，而学生的成本要更低，所以博士会被学生替代（这个问题就变得比较实际）。当学生工作两个小时，机器工作一小时也可以完成全部的四道题目，此时的成本是 $6+8=14$ 美元。由于博士并没有参与解题，所以 $x_1=0$，此时这个解在图中三角形区域的边 $QR$ 上，因此最优解应该是这四道题目全部都由学生完成（在 $Q$ 点）或全部由机器完成（在 $R$ 点）。而在本例题中，学生用四个小时来解这四道题的成本是 $$12$ —— 这是最小的成本。
由于约束条件 $A\mathbf{x}=\mathbf{b}$ 只有一个方程，所以顶点 $(0,4,0)$ 只有一个非零分量。当约束条件 $A\mathbf{x}=\mathbf{b}$ 有 $m$ 个方程时，顶点有 $m$ 个非零分量。尽管可以令其中 $n-m$ 个自由变量为零，然后求解 $A\mathbf{x}=\mathbf{b}$ 剩下的 $m$ 个变量，但是，我们并不知道应该选择求解哪 $m$ 个变量。
可行域顶点的个数就是从 $n$ 个分量中选择 $m$ 个分量方法的总数，即 $C_n^m$. 从 $n$ 个数中选择 $m$ 个数，这与概率密切相关。如有 $20$ 个未知数 $m=8$ 个方程（这些数仍然很小），可行域有 $C_{20}^8=\frac{20!}{12!8!}=125970$ 个顶点。
我们为了找到最小的成本检验 $3$ 个顶点没有什么问题，但是要检验这 $125970$ 个顶点，我们不能使用这个方法。而使用单纯形法要快得多。

对偶问题（The Daul Problem）： 在线性规划中，问题通常是成对出现的，有最大值问题和最小值问题 —— 原问题和它的 “对偶” 问题。原问题使用矩阵 $A$ 和两个向量 $\mathbf{b}$ 和 $\mathbf{c}$ 来表示，对偶问题是使用 $A$ 的转置，然后交换一下向量 $\mathbf{b}$ 和 $\mathbf{c}$ 来表示：目标是最大化 $\mathbf{b}\cdot \mathbf{y}$. 下面是原问题的对偶问题：
一个作弊者通过出售问题的答案来解决问题。 每个问题的费用是 $y$ 美元，总共是 $4y$ 美元（此时 $\mathbf{b}=\left[\begin{array}{c}4\end{array}\right]$ 进入到收入中）。很明显，这个作弊者的收费不能超过博士、学生或者机器的成本，即：$y\le 5,y\le 3$ 且 $2y\le 8$（此时 $\mathbf{c}=(5,3,8)$ 进入到了不等式约束条件中），这个作弊者需要最大化他的收入 $4y$.

对偶问题$在约束条件A^T\mathbf{y}\le \mathbf{c}下，最大化收入\mathbf{b}\cdot \mathbf{y}\le \mathbf{c}$.

当 $y=3$ 时，收入 $4y=12$ 美元达到最大，对偶问题的最大值（$$12$）等于原问题的最小值（$$12$），这就是对偶性。

对偶问题中，如果其中一个问题有最优解向量（${\mathbf{x}}^{\ast }$ 或 ${\mathbf{y}}^{\ast }$），那么另一也有最优解向量。
$$最小化成本\mathbf{c}\cdot {\mathbf{x}}^{\ast }等于最大化收入\mathbf{b}\cdot {\mathbf{y}}^{\ast }$$

本课程前面有介绍行图和列图，第一个 “对偶定理” 是关于秩的：矩阵的线性无关的最大行数等于线性无关的最大列数。这个定理和上面的对偶问题一样，在矩阵规模较小时很容易进行验证。下面证明该定理中比较容易的部分 —— 作弊者的收入 ${\mathbf{b}}^T\mathbf{y}$ 不会超过实际的解题成本：$$如果A\mathbf{x}=\mathbf{b},\mathbf{x}\ge 0,A^T\mathbf{y}\le \mathbf{c}则{\mathbf{b}}^T\mathbf{y}=(A\mathbf{x}{)}^T\mathbf{y}={\mathbf{x}}^T(A^T\mathbf{y})\le {\mathbf{x}}^T\mathbf{c}(\mathrm{10.4.1})$$完整的对偶定理说的是，当 ${\mathbf{b}}^T\mathbf{y}$ 达到最大值时，${\mathbf{x}}^T\mathbf{c}$ 达到最小值，它们相等：$\mathbf{b}\cdot {\mathbf{y}}^{\ast }=\mathbf{c}\cdot {\mathbf{x}}^{\ast }$. 注意（10.4.1）中最后一步有不等式符号 $\le$：
$$\mathbf{x}\ge 0和\mathbf{s}=\mathbf{c}-A^T\mathbf{y}\ge 0的点积{\mathbf{x}}^T\mathbf{s}\ge 0，这个就是{\mathbf{x}}^T{A}^T\mathbf{y}\le {\mathbf{x}}^T\mathbf{c}$$

$两式相等需要{\mathbf{x}}^T\mathbf{s}=0$，因此最优解对于每个 $j$ 都有 ${\mathbf{x}}_j^{\ast }=0或{\mathbf{s}}_j^{\ast }=0$.

三、单纯形法

消元法是求解线性方程组的主要方法，单纯形法（simplex method）是求解线性不等式的主要方法，这里只阐述单纯形法的原理。单纯形法是从可行域一个顶点出发搜索使成本更低的相邻顶点，最终（实际应用中会很快）到达成本最低的那个顶点。
顶点对应一个 $n$ 维的 $\mathbf{x}\ge 0$ 的向量，且满足 $m$ 个方程组 $A\mathbf{x}=\mathbf{b}$，它最多有 $m$ 个正分量，另外 $n-m$ 个分量均为零（这些对应着自由变量，带回方程组可以求得 $m$ 个基本变量的值，$\mathbf{x}$ 的所有分量一定能够都是非负的，否则就不是顶点）。对于相邻的顶点，$x$ 的一个零分量会变为正值，而一个正分量会变为零。
单纯形法的每一步都必须要决定哪个分量应该通过将其变为正值来 “换入（enter）”，哪个分量应该把它变为零来 “换出（leaves）”，决定的依据是降低总成本。这样单纯形法会向最优解 ${\mathbf{x}}^{\ast }$ 移动一步。
下面说明整体思想：观察当前顶点的每个零分量，如果将其从 $0$ 变为 $1$，而为了保持 $A\mathbf{x}=\mathbf{b}$ 成立，其它的非零分量将不得不调整，我们通过回代来求新的 $\mathbf{x}$，并计算总成本 $\mathbf{c}\cdot \mathbf{x}$ 的改变量。这个改变量就是该分量 “降低的成本（reduced cost）” $r$，换入的分量是就是取得最小负值（most negative） $r$ 的那一个，$r$ 是该分量最大的单位成本减少量。

【例1】上述的原问题，假设当前的顶点是 $P=(4,0,0)$，表示所有的工作全部由博士完成（总成本是 $$20$）. 如果学生工作一小时，$\mathbf{x}=(3,1,0)$ 对应的成本将下降到 $$18$，减少的成本是 $r=-2$ 美元. 如果机器工作一小时，则 $\mathbf{x}=(2,0,1)$ 对应的成本也是 $$18$，减少到成本也为 $r=-2$ 美元。这种情况下单纯形法可以选择学生或机器的工作时间作为换入分量。
即使在这样一个小的例子中，第一步也可能无法立即得到最优的 ${\mathbf{x}}^{\ast }$。单纯性法首先选择一个换入分量，然后才可以确定换入分量的数值。我们将换入分量从 $0$ 变为 $1$ 后再计算 $r$，但是一单位的分量通常太大或太小。本例中，单纯形法选择换出分量（博士），最优解将会向着顶点 $Q$ 或 $R$ 移动。
换入的分量的值越大，成本就越低。当其中一个正分量（我们调整其余的分量以满足 $A\mathbf{x}=\mathbf{b}$）变为零时，我们就要停止该这个换入过程。换出分量是第一个取到零的那个正分量 $x_i$，此时我们就找到了一个相邻的顶点，然后从新的顶点开始搜索下一个要换入和换出的分量。
当所有成本的改变量都为正时，当前的顶点就是最优解 ${\mathbf{x}}^{\ast }$. 此时没有零分量可以在成本 $\mathbf{c}\cdot \mathbf{x}$ 不增加的情况下变为正值，没有新的变量要换入了，这个问题就解决了（同样也求得的对偶问题的 ${\mathbf{y}}^{\ast }$）。

注： 通常需要 $\alpha n$ 步就可以得到 ${\mathbf{x}}^{\ast }$，其中 $\alpha$ 不会太大。但是也有需要经过指数级单纯性法步数才能找到最优解的例子，进而人们发展出了另一种求解方法，这个方法保证了使用更少（但更困难）的步骤求得 ${\mathbf{x}}^{\ast }$，但是需要经过可行域的内部。

【例2】最小化成本 $\mathbf{c}\cdot \mathbf{x}=3x_1+x_2+9x_3+x_4$. 约束条件是 $\mathbf{x}\ge 0$ 和下面两个方程 $A\mathbf{x}=\mathbf{b}$：$$\{\begin{array}{r}{x}_1+2x_3+x_4=4\\ x_2+x_3-x_4=2\end{array}m=2个方程，n=4个未知数$$起始顶点是 $\mathbf{x}=(4,2,0,0)$，对应的成本为 $\mathbf{c}\cdot \mathbf{x}=14$. $\mathbf{x}$ 有 $m=2$ 个非零分量和 $n-m=2$ 个零分量，该顶点零分量是 $x_3$ 和 $x_4$. 现在需要判断要换入的分量是 $x_3$ 还是 $x_4$（变为非零分量），我们先将这两个分量分别增加一个单位：$$\begin{gathered} 如果x_3=1而x_4=0，则\mathbf{x}=(2,1,1,0)对应的成本为16. \\ 如果x_4=1而x_3=0，则\mathbf{x}=(3,3,0,1)对应的成本为13. \end{gathered}$$将上述两个成本与成本 $14$ 进行比较，换入 $x_3$ 的成本改变量为 $r=2$，此时为正值，是无效换入。换入 $x_4$ 的成本改变量为 $r=-1$，此时是负值是有效换入。所以，我们换入的分量是 $x_4$.
那么 $x_4$ 可以换入多少呢？一单位的 $x_4$ 会使得 $x_1$ 从 $4$ 降到 $3$，四单位就会使得 $x_1$ 从 $4$ 降到零（同时 $x_2$ 会增加到 $6$），因此换出的分量是 $x_1$，新的顶点为 $\mathbf{x}=(0,6,0,4)$，对应的成本变为了 $\mathbf{x}\cdot \mathbf{c}=10$. 这个就是最优的 ${\mathbf{x}}^{\ast }$，但是为了证明我们还需要尝试进一步的单纯性步骤，此时我们从顶点 $(0,6,0,4)$ 出发，假设换入 $x_1$ 和 $x_3$：$$从顶点(0,6,0,4)出发\begin{array}{l}如果x_1=1而x_3=0，则\mathbf{x}=(1,5,0,3)对应的成本为11.\\ 如果x_3=1而x_1=0，则\mathbf{x}=(0,3,1,2)对应的成本为14.\end{array}$$这两个成本都比 $10$ 要高，它们的 $r$ 都是正值 —— 此时不需要移动。因此，当前的顶点 $(0,6,0,4)$ 就是最优解 ${\mathbf{x}}^{\ast }$.
这些计算可以程序化，每个单纯性步骤需要求解三个有相同系数矩阵 $B$ 的线性方程组，矩阵 $B$ 是由矩阵 $A$ 的 $m$ 个基本列向量组成的 $m\times m$ 的矩阵。当一个列换入，而旧的列换出时，有快速更新 $B^{-1}$ 的方法。这就是大多数代码处理单纯性法的方式。
最优解 ${\mathbf{x}}^{\ast }$ 中有 $m$ 个非零分量，最优解 ${\mathbf{y}}^{\ast }$ 是 $m$ 个方程 $A^T{\mathbf{y}}^{\ast }=\mathbf{c}$ 的解向量，然后得到的最优关系 ${\mathbf{x}}^T\mathbf{s}=0$ 就是对偶性：要么 $x_j^{\ast }=0$，要么 ${\mathbf{s}}^{\ast }=\mathbf{c}-A^T{\mathbf{y}}^{\ast }$ 中的 “松弛量（slack）” $s_j^{\ast }=0$.
当 ${\mathbf{x}}^{\ast }=(0,4,0)$ 是最优顶点 $Q$ 时，作弊者设定的价格为 $y^{\ast }=3$.

四、内点法

单纯形法是沿着可行域的边移动，最终到达最优顶点 ${\mathbf{x}}^{\ast }$. 内点法（Interior point methods）是在可行域内部移动（其中 $\mathbf{x}>0$），这个方法是期望能直接找到最优解 ${\mathbf{x}}^{\ast }$，并且效果很好。
固定在可行域内部搜索的一种方法是在边界处设置障碍，我们添加一个额外的对数项成本，因此当任何的 $x_j$ 趋近于零时，这个成本都会趋于无穷大，最优解向量满足 $\mathbf{x}>0$，参数 $\theta$ 的值较小，并且会趋于零。

障碍问题 $约束条件为A\mathbf{x}=\mathbf{b}时，最小化$$${\mathbf{c}}^T\mathbf{x}-\theta (\log x_1+\log x_2+\cdots +\log x_n)(\mathrm{10.4.2})$$

这个成本是非线性的（线性规划中的不等式已经体现了非线性），由于当 $x_j=0$ 时，$\log x_j$ 为无穷大，所以不再需要约束条件 $x_j\ge 0$.
对于每个 $\theta$，障碍问题都是原问题的一个近似问题，我们对这 $m$ 个约束方程 $A\mathbf{x}=\mathbf{b}$ 引入拉格朗日乘子（Lagrange multipliers）$y_1,y_2,\cdots ,y_m$，这是处理约束条件一个很好的方法：$$\mathbf{y}来自拉格朗日乘子L(x,y,\theta )={\mathbf{c}}^T\mathbf{x}-\theta \left(\sum _{i=1}^n\log x_i\right)-{\mathbf{y}}^T(A\mathbf{x}-\mathbf{b})(\mathrm{10.4.3})$$由 $\frac{\partial L}{\partial y}=0$ 我们可以反过来得到 $A\mathbf{x}=\mathbf{b}$. 偏导数 $\frac{\partial L}{\partial x_j}$ 有重大的意义！

障碍问题中的最优解 $\frac{\partial L}{\partial x_j}=c_j-\frac{\theta }{x_j}-(A^T\mathbf{y}{)}_j=0$ 即${\mathbf{x}}_j{\mathbf{s}}_j=\theta (\mathrm{10.4.4})$

原问题中有 $x_j{s}_j=0$，而障碍问题是 $x_j{s}_j=\theta$，最优解 ${\mathbf{x}}^{\ast }(\theta )$ 位于通往 ${\mathbf{x}}^{\ast }(0)$ 的中心路径（central path） 上。这 $n$ 个最优方程 $x_j{s}_j=\theta$ 是非线性的，我们可以用牛顿迭代法来求解。
当前的障碍问题中的 $\mathbf{x},\mathbf{y},\mathbf{s}$ 满足 $A\mathbf{x}=\mathbf{b},\mathbf{x}\ge 0$ 和 $A^T\mathbf{y}+\mathbf{s}=\mathbf{c}$，但是并不满足 $x_j{s}_j=\theta$，因为现在是在可行域内，还不是最优解。牛顿迭代法选取步长 $\Delta \mathbf{x},\Delta \mathbf{y},\Delta \mathbf{s}$，并且忽略掉 $(x+\Delta x)(s+\Delta s)=\theta$ 中的二次项 $\Delta \mathbf{x}\Delta \mathbf{s}$，然后使用下面的线性方程组来修正 $\mathbf{x},\mathbf{y},\mathbf{s}$：$$牛顿迭代法\begin{array}{r}A\Delta \mathbf{x}=0\\ A^T\Delta \mathbf{y}+\Delta \mathbf{s}=0\\ s_j\Delta x_j+x_j\Delta s_j=\theta -x_j{s}_j\end{array}(\mathrm{10.4.5})$$牛顿迭代法对每个 $\theta$ 都能达到二阶收敛，然后令 $\theta$ 趋近于零。对偶间隙（duality gap）${\mathbf{x}}^T\mathbf{s}$ 通常在迭代 $20$ 至 $60$ 步后会小于 $10^{-8}$.
内点法在商业软件中几乎是 “按原样” 使用，它用于解决一大类线性和非线性的问题。