OI-wiki 单纯形法全解析从标准形式到两阶段实现与竞赛实战【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki导读本文是 OI-wiki 中 docs/math/simplex.md 的深度展开系统讲解线性规划的标准解法——单纯形法simplex method。文章以 OI / ICPC 算法竞赛为背景从标准形式与基本可行解出发逐步推导转轴操作、终止条件、压缩单纯形表的更新流程与几何直觉并深入讲解松弛形式、两阶段法、大 M 法等工程实现细节最后结合仓库内可运行的参考实现simplex_0.cpp与例题NOI2008 志愿者招募给出可直接落地的完整方案。读完本文你将能够理解单纯形法的每一步数学依据掌握基于压缩单纯形表的二阶段单纯形法 C 实现并能在竞赛中正确地建模、求解与判读最优 / 无界 / 不可行线性规划问题。前置知识线性规划基础标准形式、可行域、对偶问题、互补松弛条件。引入算法竞赛中的单纯形法算法竞赛中经常使用单纯形法解决线性规划问题。不过需要注意的是OI 中遇到的线性规划问题大多具有更特殊的结构如全幺模矩阵对应的网络流模型常常可以转化为网络流问题求解因此单纯形法在竞赛中并不常用效率也不如专门为网络流设计的算法。它真正的价值在于当问题无法归约到网络流等特殊结构时单纯形法是算法竞赛中最通用、最常应用的线性规划求解方法。更完整的算法对比单纯形法、椭球法、内点法参见 docs/math/linear-programming.md。基本概念标准形式与一个完整例子假设要求解如下一个有 $n$ 个决策变量和 $mn$ 个约束的标准形式线性规划问题$$ \begin{aligned} \min_{x}; z c^Tx \ \text{subject to } Ax b, \ x \ge 0. \end{aligned} $$不妨假设这 $m$ 个等式约束确定的线性方程组有解且 $A$ 满秩则 $\operatorname{rank}A m \le n$。一个例子在严格叙述单纯形法的步骤之前先考察一个具体的例子方便理解整个迭代过程。例子考虑线性规划问题$$ \begin{aligned} \max; 10 x_1 12 x_2 12 x_3 \ \text{subject to } x_1 2 x_2 2x_3 \le 20, \ 2x_1 x_2 2x_3 \le 20, \ 2x_1 2x_2 x_3 \le 20,\ x_1,x_2,x_3 \ge 0. \end{aligned} $$通过添加松弛变量就得到它的标准形式$$ \begin{aligned} \min; -10 x_1 - 12 x_2 - 12 x_3 \ \text{subject to } x_1 2 x_2 2x_3 x_4 20, \ 2x_1 x_2 2x_3 x_5 20, \ 2x_1 2x_2 x_3 x_6 20,\ x_1,x_2,x_3,x_4,x_5,x_6 \ge 0. \end{aligned} $$观察该问题的等式约束它们其实相当于将变量 $x_4,x_5,x_6$ 由变量 $x_1,x_2,x_3$ 表示。将原问题稍微整理一下就有$$ \begin{array}{rrrrrr} \min_{x_i\ge 0} z 0 -10x_1 -12x_2 -12x_3;\ \text{subject to} x_4 20 -x_1 -2x_2 -2x_3, \ x_5 20 -2x_1 -x_2 -2x_3, \ x_6 20 -2x_1 -2x_2 -x_3. \ \end{array} $$从这个形式中可以清楚地看到如果令 $x_1x_2x_30$就得到原问题的一组可行解 $x (0,0,0,20,20,20)^T$对应的价值为 $z0$。为方便叙述称那些设为零的变量 $x_1,x_2,x_3$ 为非基变量剩下的变量 $x_4,x_5,x_6$ 为基变量。这组可行解显然不是最优解。只要适当地增加 $x_1,x_2,x_3$ 的值使得 $x_4,x_5,x_6$ 仍然是非负数就可以保持解仍然可行。而且因为目标函数中 $x_1,x_2,x_3$ 的系数都是严格的负数增加它们的值一定会降低目标函数的值。比如选择增加 $x_1$为了尽可能多地降低目标函数需要尽可能多地增加 $x_1$但为了保证解仍然可行需要保证 $x_4,x_5,x_6\ge 0$。因此 $x_1$ 最多可以增加到$$ \min\left{\dfrac{20}{1},\dfrac{20}{2},\dfrac{20}{2}\right} 10. $$此时可行解变为 $x (10,0,0,10,0,0)^T$。因为 $x_1$ 成为了基变量为了回到最初的情形三个基变量由三个非基变量表示需要选择一个新的非基变量。因为 $x_5,x_6$ 都是零可以选择其中任意一个作为非基变量不妨选择 $x_5$。将$$ x_1 10 - 0.5x_5 - 0.5x_2 - x_3 $$代入原来的问题就可以将原问题改写为$$ \begin{array}{rrrrrr} \min_{x_i\ge 0} z -100 5x_5 -7x_2 -2x_3;\ \text{subject to} x_4 10 0.5x_5 -1.5x_2 -x_3, \ x_1 10 -0.5x_5 -0.5x_2 -x_3, \ x_6 0 x_5 -x_2 x_3. \ \end{array} $$继续观察当前的目标函数非基变量 $x_3$ 的系数仍然是负数可以考虑增加 $x_3$ 的值。为了保证 $x_4,x_1,x_6\ge 0$$x_3$ 最多只能增加到$$ \min\left{\dfrac{10}{1},\dfrac{10}{1}\right} 10. $$注意到 $x_6$ 的表达式中 $x_3$ 的系数是正数所以无论怎么增加 $x_3$都不会使 $x_6$ 变为负数——这正是这次大括号中只有两项的原因。因为 $x_3$ 增加到 $10$ 时 $x_1,x_4$ 都变为零可以任选其中一个作为新的非基变量不妨选择 $x_4$。将$$ x_3 10 0.5x_5 - 1.5x_2 - x_4 $$代入上述问题问题变形为$$ \begin{array}{rrrrrr} \min_{x_i\ge 0} z -120 4x_5 -4x_2 2x_4;\ \text{subject to} x_3 10 0.5x_5 -1.5x_2 -x_4, \ x_1 0 -x_5 x_2 x_4, \ x_6 10 1.5x_5 -2.5x_2 -x_4. \ \end{array} $$代入 $x_5x_2x_40$就能从这个形式中读出当前可行解是 $x (0,0,10,0,0,10)^T$对应价值为 $z-120$。重复之前的操作因为 $x_2$ 的系数是负数可以增加它的值但为了保持 $x_3,x_6$ 非负只能增加到$$ \min\left{\dfrac{10}{1.5},\dfrac{10}{2.5}\right} 4. $$最小值出现在变量 $x_6$ 的表达式中所以它在 $x_24$ 时变为零。将$$ x_2 4 0.6x_5 - 0.4x_6 - 0.4x_4 $$代入上述问题原问题改写为$$ \begin{array}{rrrrrr} \min_{x_i\ge 0} z -136 1.6x_5 1.6x_6 3.6x_4;\ \text{subject to} x_3 4 -0.4x_5 0.6x_6 -0.4x_4, \ x_1 4 -0.4x_5 -0.4x_6 0.6x_4, \ x_2 4 1.5x_5 -2.5x_6 -x_4. \ \end{array} $$仍然令非基变量 $x_5,x_6,x_4$ 为零得到当前可行解为 $x (4,4,4,0,0,0)^T$对应价值为 $z-136$。因为目标函数中所有非基变量的系数都是正数无法继续改进目标函数所以当前可行解就是最优解算法终止。这个例子中算法从一组可行解出发不断地改进目标函数直到无法继续改进——这就是单纯形法的基本思想。基本可行解由于 $A$ 是满秩的总是可以选取大小为 $m$ 的子集 $B\subseteq{1,2,\cdots,n}$ 使得 $A_B$ 是可逆方阵。由此可以将 $x_B$ 由剩下的变量 $x_N$ 表示$$ x_B A_B^{-1}b - A_B^{-1}A_Nx_N. $$其中 $N{1,2,\cdots,n}\setminus B$矩阵 $A_B,A_N$ 分别为矩阵 $A$ 中标号 $i\in B$ 和标号 $i\in N$ 的列组成的子矩阵向量 $x_B,x_N$ 分别为向量 $x$ 中相应分量的子向量。如果 $i\in B$称 $x_i$ 为基变量basic variable否则称 $x_i$ 为非基变量non-basic variable。基变量的全体称为一组基basis本文用对应的标号集合 $B$ 表示一组基。提示「基」「基」这个名称可以从线性代数的角度理解设 $A$ 的全体列向量张成的线性空间为 $V$那么基 $B$ 对应的列向量就是空间 $V$ 的一组基。在基变量 $x_B$ 的表达式中令 $x_N0$就得到全体等式约束的一组解$$ x (x_B,x_N) (A_B^{-1}b,0). $$这样得到的解称为线性规划问题的一个基本解basic solution。如果它还满足所有非负约束即 $x\ge 0$那么它也是原问题的一个可行解称为基本可行解basic feasible solution, BFS。在单纯形法的迭代过程中需要始终保持当前的解为一组基本可行解。转轴单纯形法的每次迭代称为一次转轴pivoting。从结果上看每次转轴总是移除一个旧的基变量、再添加一个新的基变量进而改进目标函数的值。提示「转轴」「转轴」这个名称同样可以从线性代数的角度理解基 $B$ 对应的列向量是空间 $V$ 的一组基它们对应着在该基的表示下空间 $V$ 的一组坐标轴。因此转轴的过程就是将某条坐标轴旋转到新位置的过程。为了确定需要添加的基变量将目标函数利用非基变量表示为$$ \begin{aligned} c^Tx c^T_Bx_B c^T_Nx_N \ c_B^TA_B^{-1}b (c_N^T - c_B^TA_B^{-1}A_N)x_N. \end{aligned} $$令 $x_N0$就得到目标函数在当前基本可行解处的价值 $zc_B^TA_B^{-1}b$。表达式第二项的系数表示 $x_N$ 改变时目标函数的变化$$ \tilde c_N \dfrac{\partial z}{\partial x_N} c_N - A_N^T(A_B^{-1})^Tc_B. $$注意到 $c_B - A_B^T(A_B^{-1})^Tc_B 0$所以可以记向量$$ \tilde c (\tilde c_B^T,\tilde c_N^T)^T c - A^T(A_B^{-1})^Tc_B $$为线性规划问题在基本可行解 $x$ 处的约化成本reduced cost。分量 $\tilde c_i0$ 说明增加变量 $x_i$ 的值可以改进原问题的目标函数。这样的变量只能是一个非基变量它称为本次转轴的入基变量entering variable——因为在转轴后 $x_i$ 将变为基变量不再恒设为零但依然有可能等于零。选择完入基变量后还需要选择要移除的旧基变量。为此只需确定在增加 $x_i$ 的过程中哪个现有的基变量最先变为零。将 $x_N(x_i,x_{N\setminus{i}})(x_i,0)$ 代入 $x_B$ 的表达式就有$$ x_B A_B^{-1}b - A_B^{-1}A_ix_i. $$因此 $x_i$ 可以增加的最大量等于$$ \theta \min\left{\dfrac{(A_B^{-1}b)_j}{(A_B^{-1}A_i)_j}:(A_B^{-1}A_i)_j0\right}. $$最先变为零的变量就是使得该表达式取得最小值的下标 $j$ 对应的基变量 $x_{B_j}$。它是增加 $x_i$ 过程中的「瓶颈」——继续增加 $x_i$ 将使得 $x_{B_j}$ 成为负值。这一变量就是本次转轴的出基变量leaving variable确定出基变量的方法称为最小比值检验minimum ratio test。设入基变量为 $x_i$出基变量为 $x_{i}$转轴之后基变量就是 $x_{B\setminus{i}\cup{i}}$非基变量就是 $x_{N\setminus{i}\cup{i}}$。终止条件单纯形法就是从一组基本可行解出发不断转轴的过程。上一节的讨论并不完整它忽略了一些特殊情形有些对应着算法终止有些则需要额外处理。入基变量不存在即 $\tilde c\ge 0$。此时当前基本可行解就是最优解算法终止。要严格证明这一点需要用到互补松弛条件。令 $y(A_B^{-1})^Tc_B$。注意到在整个算法过程中始终保持 $x$ 是可行解且互补松弛条件成立即$$ x^T(c-A^Ty) \tilde c^Tx \tilde c_B^Tx_B \tilde c_N^Tx_N 0. $$因此只要 $y$ 是对偶问题的可行解即 $A^Ty\le c$就能得到 $x$ 和 $y$ 分别是原问题和对偶问题的最优解。这个条件就是 $\tilde c\ge 0$即不存在入基变量。提示「影子价格」向量 $y(A_B^{-1})^Tc_B$ 常称为对偶向量dual vector。当 $B$ 对应的基本可行解是原问题的最优解时$y$ 是对偶问题的最优解。因此利用单纯形法求解原问题时也会顺带获得对偶问题的最优解。因为 $y$ 是当前价值关于约束常量的偏导数即 $\dfrac{\partial(c^Tx)}{\partial b} (A_B^{-1})^Tc_B y$所以它也称为影子价格shadow price。出基变量不存在即 $A_B^{-1}A_i\le 0$。此时转轴不存在任何「瓶颈」可以无限增加 $x_i$ 来改进目标函数直到 $-\infty$说明原问题无界算法终止。入基变量和出基变量的选择可能不唯一。不恰当的选择方式可能导致过多转轴甚至使算法陷入循环而无法终止。这类情形的处理比较复杂需要应用转轴规则防止循环并减少转轴次数。单纯形表具体实现转轴过程时只需要维护每次转轴之后线性规划问题的系数矩阵$$ \tilde T_B \begin{pmatrix} -z_B \tilde c^T_N \ x A_B^{-1}A_N \end{pmatrix}\begin{pmatrix} -c_B^TA_B^{-1}b c^T - c_B^TA_B^{-1}A_N \ A_B^{-1}b A_B^{-1}A_N \end{pmatrix}. $$它对应着线性规划问题$$ \begin{array}{rrrr} \min_{x\ge 0} c_B^TA_B^{-1}b \tilde c_N^Tx_N; \ \text{subject to} x_B A_B^{-1}b - A_B^{-1}A_Nx_N. \end{array} $$矩阵 $\tilde T_B$ 称为该线性规划问题相对于基 $B$ 的压缩单纯形表condensed simplex tableau表的左上角 $(\tilde T_B)_{00}$ 为当前解的价值的相反数第 $0$ 行、第 $i$ 列的量 $(\tilde T_B){0i}$ 为第 $i$ 个非基变量 $x{N_i}$ 的约化成本第 $j$ 行、第 $0$ 列的量 $(\tilde T_B){j0}$ 为第 $j$ 个基变量 $x{B_j}$ 的值$A_B^{-1}A_N$ 就是用非基变量 $x_N$ 表示基变量 $x_B$ 的表达式中的系数。容易看出转轴需要的所有信息都可以从压缩单纯形表中直接获得。利用压缩单纯形表单次转轴包括如下操作选取列 $i1,\cdots,n-m$使得 $(\tilde T_B){0i}0$。如果不存在这样的 $i$那么当前解就是最优解量 $-(\tilde T_B){00}$ 就是最优价值。选取行 $j1,\cdots,m$使得 $(\tilde T_B){ji}0$ 且 $(\tilde T_B){j0}/(\tilde T_B)_{ji}$ 最小。如果不存在这样的 $j$那么原问题无界。令变量 $x_{N_i}$ 入基、变量 $x_{B_j}$ 出基并更新单纯形表。更新单纯形表在更新单纯形表之前第 $j$ 行表示等式$$ x_{B_j} (\tilde T_B){j0} - \sum{i1}^{n-m}(\tilde T_B){ji}x{N_i}. $$为了更新单纯形表需要用 $x_{N\setminus{N_i}\cup{B_j}}$ 表示 $x_{N_i}$也就是$$ x_{N_i} \dfrac{(\tilde T_B){j0}}{(\tilde T_B){ji}} - \dfrac{1}{(\tilde T_B){ji}}x{B_j} - \sum_{i\neq i}\dfrac{(\tilde T_B){ji}}{(\tilde T_B){ji}}x_{N_{i}}. $$将它代入其余的式子中就得到$$ x_{B_{j}} \left((\tilde T_B){j0} - (\tilde T_B){ji}\dfrac{(\tilde T_B){j0}}{(\tilde T_B){ji}}\right) \dfrac{(\tilde T_B){ji}}{(\tilde T_B){ji}}x_{B_j} - \sum_{i\neq i}\left((\tilde T_B){ji}-(\tilde T_B){ji}\dfrac{(\tilde T_B){ji}}{(\tilde T_B){ji}}\right)x_{N_i}. $$第 $0$ 行类似只是等式左侧变为 $-z$。虽然式子看起来复杂但实现时只需要分两步更新第 $j$ 行令 $\alpha(\tilde T_B)_{ji}$将第 $i$ 列数字置为 $1$然后将整行所有数字同除以 $\alpha$更新第 $j\neq j$ 行令 $\beta(\tilde T_B)_{ji}$将第 $j$ 列数字置为 $0$然后将整行数字同时减去 $\beta$ 倍的第 $j$ 行数字。提示「单纯形表」单纯形表simplex tableau指矩阵$$ T_B \begin{pmatrix} -z \tilde c^T \ x A_B^{-1}A \end{pmatrix}\begin{pmatrix} -c_B^TA_B^{-1}b c^T - c_B^TA_B^{-1}A \ A_B^{-1}b A_B^{-1}A \end{pmatrix}. $$相较于压缩单纯形表它多了 $m$ 列分别对应 $m$ 个基变量且对应第 $j$ 个基变量的列一定是 $e_j$该向量只在第 $j$ 行处取值为 $1$其余行均为 $0$。因为这些列并没有提供多余的信息实现单纯形法时常省略这些列就得到了压缩单纯形表。利用单纯形表可以更方便地理解更新步骤因为所有单纯形表 $T_B$ 都可以由同一个矩阵 $T_0$ 左乘一个与基有关的可逆矩阵 $L_B$ 得到即$$ T_B \begin{pmatrix} -c_B^TA_B^{-1}b c^T - c_B^TA_B^{-1}A \ A_B^{-1}b A_B^{-1}A \end{pmatrix}\begin{pmatrix} 1 -c_B^TA_B^{-1} \ O A_B^{-1} \end{pmatrix} \begin{pmatrix} 0 c^T \ b A \end{pmatrix}L_BT_0, $$所以这些单纯形表和 $T_0$ 之间可以通过若干次初等行变换相互转化。更新单纯形表时只需要施行初等行变换使得入基变量对应列变为 $e_j$ 即可迁移到压缩单纯形表上就是上文给出的两步操作。更新压缩单纯形表的参考实现位于 docs/math/code/simplex/simplex_0.cpppivot函数// Pivot on (N[x], B[y]). void pivot(int x, int y) { std::swap(N[x], B[y]); long double v -1 / tab[x][y]; for (int j 0; j m 1; j) { tab[x][j] j y ? -v : v * tab[x][j]; } for (int i 0; i n; i) { if (i x) continue; v tab[i][y]; tab[i][y] 0; for (int j 0; j m 1; j) { tab[i][j] v * tab[x][j]; } } }由该实现可知单次更新单纯形表的时间复杂度为 $O(mn)$。稍后讨论转轴规则时会说明确定出基变量和入基变量的复杂度同样不会超过 $O(mn)$因此单次转轴的时间复杂度就是 $O(mn)$。例子的完整单纯形表推演为方便理解此处列出前文例子利用压缩单纯形表计算的详细步骤。初始时压缩单纯形表如下$$ \begin{array}{|l|c|ccc|} \hline x_1 x_2 x_3 \ \hline 0 -10 -12 -12 \ \hline x_4 20 1 2 2 \ x_5 20 2 1 2 \ x_6 20 2 2 1 \ \hline \end{array} $$根据第 $0$ 行的约化成本可以选择 $x_1,x_2,x_3$ 入基令 $x_1$ 入基。再根据最小比值检验可以选择 $x_5,x_6$ 出基令 $x_5$ 出基。更新压缩单纯形表$$ \begin{array}{|l|c|ccc|} \hline x_5 x_2 x_3 \ \hline 100 5 -7 -2 \ \hline x_4 10 -0.5 1.5 1 \ x_1 10 0.5 0.5 1 \ x_6 0 -1 1 -1 \ \hline \end{array} $$根据第 $0$ 行的约化成本可以选择 $x_2,x_3$ 入基令 $x_3$ 入基。再根据最小比值检验可以选择 $x_4,x_1$ 出基令 $x_4$ 出基。更新压缩单纯形表$$ \begin{array}{|l|c|ccc|} \hline x_5 x_2 x_4 \ \hline 120 4 -4 2 \ \hline x_3 10 -0.5 1.5 1 \ x_1 0 1 -1 -1 \ x_6 10 -1.5 2.5 1 \ \hline \end{array} $$根据第 $0$ 行的约化成本只能选择 $x_2$ 入基令 $x_2$ 入基。再根据最小比值检验只能选择 $x_6$ 出基令 $x_6$ 出基。利用前述初等行变换更新单纯形表$$ \begin{array}{|l|c|ccc|} \hline x_5 x_4 x_6 \ \hline 136 1.6 3.6 1.6 \ \hline x_3 4 0.4 0.4 -0.6 \ x_1 4 0.4 -0.6 0.4 \ x_2 4 -0.6 0.4 0.4 \ \hline \end{array} $$根据第 $0$ 行的约化成本不存在入基变量。因此当前解 $x(4,4,4,0,0,0)^T$ 就是最优解最小化问题的最优价值为 $-136$。除了利用单纯形表实现单纯形法之外还可以使用修正单纯形法revised simplex method它进一步改进了时空复杂度将单次更新的复杂度降低到 $O(m^2)$在 $m\ll n$ 或 $A$ 是稀疏矩阵的情形下尤为高效。几何背景顶点、边与转轴的直观理解对线性规划问题的可行域$$ \mathcal D {x\in\mathbf R^n : Ax b,~ x\ge 0} $$的分析指出线性规划问题的最优解如果存在必然可以选取为可行域 $\mathcal D$ 的顶点。求解线性规划问题就转化为在所有顶点解里找到价值函数最优的那个。每一个顶点的坐标都可以通过 $n$ 个紧约束联立得到的方程组求解得到。对于标准形式的约束所有 $m$ 个等式约束一定是紧的剩下的 $n-m$ 个约束只能从非负约束中选取。选取这些非负约束作为紧约束就相当于将相应的决策变量 $x_N$ 设为 $0$相应地方程组 $Ax b$ 退化为关于剩余 $m$ 个决策变量 $x_B$ 的线性方程组 $A_Bx_B b$。只要 $A_B$ 可逆就可以解得 $x_B A_B^{-1}b$从而得到一个解 $(x_B,x_N)(A_B^{-1}b,0)$如果 $x_B\ge 0$这就是 $\mathcal D$ 的一个顶点坐标。容易看出顶点解的概念和前文定义的基本可行解是一致的。因此只要在所有基本可行解内找到最优的那个就能获得原问题的最优解。虽然这大幅简化了问题但可行域的顶点个数是指数级的穷举并不现实。为了解决这一困难可以考虑沿着可行域的边移动从一个顶点移动到与之相邻的顶点。因为相邻的顶点必定位于同一条边上所以它们至少满足 $n-1$ 条相同的紧约束也就是说相邻顶点对应的紧约束能且仅能相差一个。因此对于一个基本可行解 $x$只要将它的一个基变量换成一个非基变量就能得到一个相邻的adjacent基本可行解 $x$——这正是转轴操作。所以单纯形法从一个基本可行解出发、不断转轴改进目标函数的过程其实就是在相应的可行域上从一个顶点出发、不断向相邻顶点移动进而改进目标函数的过程。以下图为例来自 docs/math/images/simplex-geo.svg本文讨论的例子中可行域是一个有五个顶点的三维多面体。前述求解过程从几何直观上看就对应着多面体的顶点间的如下路径$$ (0,0,0) \rightarrow (0,0,10) \rightarrow (10,0,0) \rightarrow (4,4,4). $$实现细节从松弛形式到两阶段法利用单纯形表已经能够求解许多线性规划问题但对于最一般的情形仍有许多细节值得深入探讨。松弛形式保证矩阵满秩将一般形式的线性规划问题转化为标准形式的方法已经讨论过了。但为了方便利用单纯形法求解还需要保证系数矩阵 $A$ 满秩。虽然先转换为标准形式再消去线性相关的约束是可行的但为了求解简便通常采用如下策略将线性规划问题转化为不等式形式inequality form即 $\min{c^Tx : Ax \le b,~ x \ge 0}$ 的形式通过添加松弛变量 $s$将问题转化为标准形式$\min{c^Tx : Ax s b,~ x\ge 0,~ s \ge 0}$。这样做的好处是最终得到的标准形式的系数矩阵 $(A,I)$总是满秩的且总是存在未必可行的基本解 $(x,s)(0,b)$。这种特殊的标准形式也称为松弛形式slack form。初始基本可行解两阶段法与人工变量前文描述的单纯形法总是假定已知一组基本可行解。有些时候这很容易例如如果上述松弛形式中 $b\ge 0$那么 $(x,s)(0,b)$ 就是一组基本可行解前文的数值例子就是这种情形。对于一般的情形可以采用两阶段法two-phase method第一阶段求解一个可行性线性规划问题获得原问题的一个基本可行解第二阶段从这个基本可行解出发应用单纯形法求解原问题。假设有标准形式的问题 $\min{c^Tx : Ax b \ge 0,~ x\ge 0}$。在第一阶段中需要求解问题$$ \min{1^Tx_a : Ax x_a b,~ x\ge 0,~ s\ge 0}. $$这本质上是可行性线性规划问题其中新添加的变量 $x_a$ 称为人工变量artificial variable。它一定有基本可行解 $(x,x_a)(0,b)$所以可以直接用单纯形法求解如果该问题的最优价值严格大于 $0$那么不存在 $x\ge 0$ 使得 $Axb$即原问题不可行如果该问题的最优价值等于 $0$那么它的最优解中人工变量只能是零。如果仍有一些人工变量是基变量可以通过若干次转轴将它们出基。当所有人工变量都是非基变量时第一阶段得到的基本解就可以用作第二阶段的初始基本可行解。不显式引入人工变量的第一阶段实现实现第一阶段时没有必要显式地引入人工变量。对于任意选取的初始基 $B$有$$ x_B A_B^{-1}A_Nx_N A_B^{-1}b. $$如果 $(A_B^{-1}b)j\ge 0$那么无需引入人工变量否则需要额外引入人工变量 $x^{-}{B_j}$即$$ x_{B_j} - x^-{B_j} (A_B^{-1}A_N){(j)}x_N (A_B^{-1}b)_j. $$记下标集合 $L:{j:(A_B^{-1}b)_j0}$那么一阶段的单纯形表都是如下单纯形表经过若干初等行变换得到的$$ \begin{array}{|r|c|cccc|} \hline x_N x_{B_{\sim L}} x_{B_L} x_{B_{L}}^- \ \hline 0 0^T 0^T 0^T 1^T \ \hline x_{B_{\sim L}} (A_B^{-1}b){\sim L} (A_B^{-1}A_N){(\sim L)} I O O \ x_{B_{L}}^- (A_B^{-1}b)L (A_B^{-1}A_N){(L)} O I -I \ \hline \end{array} $$类似于将单纯形表简化为压缩单纯形表可以将其适当简化将 $x_{B_{L}}^-$ 行左乘以 $1^T$ 然后加到第 $0$ 行上然后略去后三列$$ \begin{array}{|r|c|cccc|} \hline x_N \ \hline 1^Tb_L 1^T(A_B^{-1}A_N)L \ \hline x{B_{\sim L}} (A_B^{-1}b){\sim L} (A_B^{-1}A_N){(\sim L)} \ -x_{B_{L}} (A_B^{-1}b)_L (A_B^{-1}A_N)_L \ \hline \end{array} $$这与正常的压缩单纯形表大体一致只是最后一行的变量上标记了负号表示该行仍然含有人工变量即原来的松弛变量仍然不可行。利用该表转轴过程如下如果 $L\varnothing$算法终止。否则根据第 $0$ 行的约化成本为负这一条件选取入基变量 $x_{N_i}$。如果不存在则原问题不可行算法终止。再根据第 $i$ 列选取出基变量 $x_{B_j}$。仍然利用最小比值检验但要求同时保证当前的可行变量可行、当前的不可行变量不可行即选取$$ \arg\min_{j}\left{\dfrac{(\tilde T_B){j0}}{(\tilde T_B){ji}}:(j\notin L\land(\tilde T_B){ji}0)\lor(j\in L\land(\tilde T_B){ji}0)\right} $$作为出基变量所在行。如果存在多个这样的出基变量优先选取不可行的出基变量。令 $x_{N_i}$ 入基、$x_{B_j}$ 出基并更新单纯形表。如果 $j\in L$那么将 $j$ 移出 $L$即取消该行的负号标记并且将 $(\tilde T_B)_{0i}$ 加一。之所以可以略去人工变量所在列是因为如果它们仍然是基变量那么它们对应的列就是 $e_j$无需记录如果它们不再是基变量那么它们就不会再次入基也无需记录。人工变量出基时需要替换成相应的非人工变量这正是最后一步的目的。参考实现位于 docs/math/code/simplex/simplex_0.cppinitialize函数// First phase: find an initial BFS. // Return false if no feasible solution. bool initialize() { int neg_count 0; for (int j 0; j m; j) { if (tab[n][j] -eps) { for (int i 0; i n; i) { tab[i][m 1] tab[i][j]; } B[j] ~B[j]; neg_count; } } while (neg_count) { int x -1; long double mi -eps; for (int i 0; i n; i) { if (tab[i][m 1] mi) { x i; mi tab[i][m 1]; } } if (x -1) return false; int y -1; mi INFINITY; for (int j 0; j m; j) { if ((B[j] 0 tab[x][j] -eps) || (B[j] 0 tab[x][j] eps)) { auto tmp tab[n][j] / tab[x][j] (B[j] 0 ? -eps : eps); if (tmp mi) { y j; mi tmp; } } } if (B[y] 0) { --neg_count; B[y] ~B[y]; pivot(x, y); tab[x][m 1] 1; } else { pivot(x, y); } } return true; }注意在一阶段开始前额外添加一行用于记录一阶段的目标函数。转轴时对整个表进行转轴包括第二阶段的目标函数。这样在第一阶段完成时第二阶段的目标函数也一并相应地更新可以直接开始二阶段的单纯形法。大 M 法两阶段法也可以通过一次单纯形法实现只需要取充分大的正数 $M$直接求解问题$$ \min{c^Tx M1^Tx_a : Ax x_a b,~ x\ge 0,~ s\ge 0} $$就可以得到原问题的最优解。实现时并不会赋予 $M$ 具体的数值而是将它视为一个未知的充分大的正数进行运算。这种方法称为大 $M$ 法big $M$ method。警告朴素算法的实际效率是指数级的因为松弛形式总是存在初始基本解只是未必可行一种简单的寻找初始基本可行解的想法是从一个不可行的基本解出发反复利用转轴操作将不可行的基变量出基并选择对应行中同样为负的数字对应列的非基变量入基直到所有基变量都是非负数为止。参考实现见 docs/math/code/simplex/simplex_2.cpp 的initialize函数。这样做虽然简单但相较于二阶段法它没有一个描述当前基的不可行程度的目标函数因此缺乏明确的改进方向。实际测试可以发现相较于通常只需要 $O(m)$ 转轴次数的二阶段法或大 $M$ 法这一朴素算法通常需要 $O(2^m)$ 的转轴次数且当 $n,m$ 较大时容易陷入循环。虽然朴素算法转轴次数中的常数很小但仅仅适用于 $n,m50$ 的情形。转轴规则转轴时如果出现多个可选的入基变量或出基变量就需要用到转轴规则pivot rule来决定选择哪个变量入基或出基。利用单纯形表本节讨论的所有规则都能在 $O(mn)$ 时间内找到入基和出基变量因此单次转轴的时间复杂度仍然是 $O(mn)$。入基变量的选择往往决定了算法终止前转轴的次数常见规则如下选择最先找到的入基变量选择标号最小的入基变量Bland 规则的一部分选择约化成本绝对值即 $|c_i|$最大的入基变量Dantzig 规则选择单次转轴价值函数改进即 $|c_i|\theta_i$最大的入基变量选择对应着最陡的边的入基变量即沿着边移动单位长度引起的价值函数改进即 $|c_i|/|A_B^{-1}A_i|$最大随机选择一个入基变量。实践中最陡边规则的效率最高。通常认为适当的转轴规则可以在大致 $2m$ 次转轴内得到大多数问题的最优解。但是对于目前所有已知的转轴规则都存在特殊构造的例子最经典的是 Klee–Minty 立方体反例能将转轴次数卡到指数级。这也正是单纯形法实践中运行效率相当优秀、但理论最差复杂度是指数级的原因。出基变量的选择往往决定了算法是否会陷入循环。如果存在多个最优价值相同的基本可行解算法就有可能一直在这些基本可行解之间循环。这些情形并不常见因此很多单纯形法的实现并不会指定出基变量的选择规则。常见的避免循环的规则有两种Bland 规则总是选择标号最小的入基变量和出基变量。字典序规则总是选择 $(A_B^{-1}A_i)_j0$ 且$$ \left(\dfrac{(A_B^{-1}b)_j}{(A_B^{-1}A_i)j},\dfrac{(A_B^{-1}){j1}}{(A_B^{-1}A_i)j},\cdots,\dfrac{(A_B^{-1}){jm}}{(A_B^{-1}A_i)_j}\right) $$字典序最小的行号 $j$ 对应的出基变量 $x_{B_j}$。此时入基变量的选择不重要。注意如果线性规划问题是松弛形式的那么这些量都可以从前文所述形式的单纯形表 $T_B$ 中直接找到否则可以在找到一组初始基本解未必可行后利用这些初始基中的基变量对应列顺序保持固定的系数作为 $A_B^{-1}$ 的系数。Bland 规则效率很低因为规则本身对入基和出基变量的选择方法是一样的很容易造成同一个变量反复入基再出基。相对来说字典序规则更为实用——它相当于对线性规划问题中的参数进行微扰使得不存在最优价值相同的基本可行解也就不存在循环的可能性。参考实现基于压缩单纯形表的二阶段单纯形法仓库在 docs/math/code/simplex/simplex_0.cpp 提供了一份完整的、可直接编译运行的参考实现对应 Luogu P13337【模板】线性规划。其整体结构如下int m, n; // Number of constraints and variables. std::vectorstd::vectorlong double tab; // Compressed tableau (transposed) with // first-phase objective attached. std::vectorint B, N; // Basic and nonbasic variables. constexpr long double eps 1e-12l; // Precision.关键设计点使用转置存储的压缩单纯形表tab尺寸为 $(n1)\times(m2)$前 $n$ 行对应决策变量第 $n$ 行是目标函数右侧常量列额外的一列用于附着第一阶段的辅助目标B、N分别维护基变量与非基变量的标号基变量初始化为 $n,n1,\dots,nm-1$松弛变量的标号非基变量初始化为 $0,\dots,n-1$使用long double与阈值eps 1e-12处理浮点误差solve()串联三个阶段initialize()返回false表示不可行输出Infeasiblesimplex()返回false表示无界输出Unbounded否则输出最优价值与最优解。其中第二阶段的主循环如下// Second phase: find an optimal BFS. // Return false if the problem is unbounded. bool simplex() { while (true) { int x -1; long double mi -eps; for (int i 0; i n; i) { if (tab[i][m] mi) { x i; mi tab[i][m]; } } // No column with a negative reduced cost Optimal. if (x -1) break; int y -1; mi INFINITY; for (int j 0; j m; j) { if (tab[x][j] eps) continue; if (tab[n][j] / tab[x][j] mi) { y j; mi tab[n][j] / tab[x][j]; } } // No row with a positive ratio Unbounded. if (y -1) return false; pivot(x, y); } return true; }输入格式与 docs/math/examples/simplex/simplex_0.in 一致第一行为变量数 $n$ 与约束数 $m$第二行给出 $n$ 个目标函数系数随后是 $m$ 个约束每个约束为一列数据先读入 $0\sim n-1$ 行的变量系数再读入第 $n$ 行的约束右侧常量 $b$。仓库同时提供了对应的标准答案文件 docs/math/examples/simplex/simplex_0.ans最优价值与每个决策变量的取值可用于自测。输出判定-1输出Infeasible1输出Unbounded0输出最优价值tab[n][m]以及各决策变量取值只有当前是基变量时才取tab[n][j]的值否则为 0。例题NOI2008 志愿者招募仓库在 docs/math/code/simplex/simplex_1.cpp 提供了「NOI2008」志愿者招募」的完整实现。题目如下总共 $n$ 天活动需要招募志愿者第 $i$ 天至少需要 $b_i$ 位志愿者共有 $m$ 类志愿者第 $j$ 类志愿者可以服务的日期为连续区间 $[l_j,r_j]$单位招募成本为 $c_i$。求最优招募方案使招募志愿者的成本最低。建模设第 $j$ 类志愿者招募 $x_j$ 位可以列出线性规划问题$$ \begin{align*} \max_{x}; \sum_{j1}^mc_jx_j \ \text{subject to } \sum_{i1}^n a_{ij}x_j \ge b_i,~i1,\cdots,n,\ x_j\ge 0,~j1,\cdots,m. \end{align*} $$其中系数$$ a_{ij} \begin{cases} 1, l_j\le i\le r_j,\ 0, \text{otherwise.} \end{cases} $$原问题没有显然的初始可行解。因此不妨考虑其对偶问题$$ \begin{align*} \min_{y}; \sum_{i1}^n b_iy_i \ \text{subject to } \sum_{j1}^na_{ij}y_i \le c_j,~j1,\cdots,m,\ y_i\ge 0,~i1,\cdots,n. \end{align*} $$通过添加松弛变量容易得到一组初始可行解可以直接略过一阶段通过单纯形法求解。根据对偶原理得到的解就是原问题的解。这种「原问题没有显然的初始可行解就转而对偶化以获取初始可行解」的技巧是竞赛中使用单纯形法的重要实战手段。simplex_1.cpp的输入为 $n,m$、$n$ 个 $c_i$随后 $m$ 行每行给出 $l_j,r_j,b_j$利用区间系数矩阵的稀疏结构for (int j s; j t; j) a[i][j] 1;建出对偶问题后直接运行单纯形法最终以(int)(simplex() 0.5)输出整数答案。复杂度与适用性小结| 环节 | 复杂度 | 说明 | | :-: | :-: | :- | | 单次转轴单纯形表更新 | $O(mn)$ | 来自 simplex_0.cpp 的pivot实现 | | 单次转轴修正单纯形法 | $O(m^2)$ | 在 $m\ll n$ 或 $A$ 稀疏时尤为高效 | | 二阶段法 / 大 $M$ 法转轴次数 | 通常 $O(m)$ | 实践表现优秀 | | 朴素初始化转轴次数 | 通常 $O(2^m)$ | 仅适用于 $n,m50$见 simplex_2.cpp |需要注意的是单纯形法最差情形复杂度是指数级的Klee–Minty 反例可以卡掉所有已知转轴规则但它仍是算法竞赛中最常应用且实践效率相当优秀的线性规划求解方法。习题以下习题可用于巩固本文内容均为线性规划建模与单纯形法实战Luogu P13337【模板】线性规划配套参考实现为 docs/math/code/simplex/simplex_0.cpp测试数据见 docs/math/examples/simplex/UOJ#179. 线性规划配套朴素初始化版本见 simplex_2.cpp其头部注释注明已在 UOJ 上验证通过Luogu P4232 无意识之外的捉迷藏Codeforces 1430 G. Yet Another DAG ProblemAtCoder Beginner Contest 231 H - Minimum Coloring参考资料线性规划之单纯形法【超详解 图解】2016 国家集训队论文算法导论Matoušek, Jiří, and Bernd Gärtner. Understanding and using linear programming. Vol. 1. Berlin: Springer, 2007.Inayatullah, Syed, Nasir Touheed, and Muhammad Imtiaz. A streamlined artificial variable free version of simplex method. PloS one 10, no. 3 (2015): e0116156.Floudas, Christodoulos A., and Panos M. Pardalos, eds. Encyclopedia of optimization. Springer Science Business Media, 2008.【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考