core 设计
linear-program 把变量、线性表达式、目标函数、约束和两阶段单纯形求解器放在 src 下的一个包中。本页给出求解器实现的数学内容，推导其停止判据，并说明代码如何围绕它们组织。
设计目标
让用户以教科书的记法编写线性规划：命名变量、带优化方向的目标函数，以及等式和不等式约束。
用教科书算法——在稠密单纯形表上的两阶段单纯形法——求解，并公开每一步（Lp::pivot、Lp::simplex_iteration、Lp::phase_1、Lp::phase_2），以便跟踪和讲授算法。
通过 luna-generic 的 trait 保持系数类型泛型，并把单纯形表存为 linear-algebra 的矩阵。
数学背景
变量 x ∈ R n x \in \mathbb{R}^n x ∈ R n 标准形
min x ∈ R n c ⊤ x subject to A x = b , x ≥ 0 , A ∈ R m × n , b ∈ R ≥ 0 m . \min_{x \in \mathbb{R}^n} \; c^\top x \quad \text{subject to} \quad A x = b, \quad x \ge 0, \qquad A \in \mathbb{R}^{m \times n},\; b \in \mathbb{R}^m_{\ge 0}. x ∈ R n min c ⊤ x subject to A x = b , x ≥ 0 , A ∈ R m × n , b ∈ R ≥ 0 m . 用本包编写的任何规划都可以在不改变最优解的前提下化为这种形式。Lp::to_standard 执行以下改写：
max c ⊤ x = − min ( − c ) ⊤ x , a ⊤ x = β , β < 0 ⟺ ( − a ) ⊤ x = − β , a ⊤ x ≤ β , β ≥ 0 ⟺ a ⊤ x + y = β , y ≥ 0 , a ⊤ x ≤ β , β < 0 ⟺ ( − a ) ⊤ x − y = − β , y ≥ 0 , a ⊤ x ≥ β , β > 0 ⟺ a ⊤ x − y = β , y ≥ 0 , a ⊤ x ≥ β , β ≤ 0 ⟺ ( − a ) ⊤ x + y = − β , y ≥ 0. \begin{aligned}
\max\; c^\top x &= -\min\; (-c)^\top x, \\
a^\top x = \beta,\ \beta < 0 \quad &\Longleftrightarrow \quad (-a)^\top x = -\beta, \\
a^\top x \le \beta,\ \beta \ge 0 \quad &\Longleftrightarrow \quad a^\top x + y = \beta,\ y \ge 0, \\
a^\top x \le \beta,\ \beta < 0 \quad &\Longleftrightarrow \quad (-a)^\top x - y = -\beta,\ y \ge 0, \\
a^\top x \ge \beta,\ \beta > 0 \quad &\Longleftrightarrow \quad a^\top x - y = \beta,\ y \ge 0, \\
a^\top x \ge \beta,\ \beta \le 0 \quad &\Longleftrightarrow \quad (-a)^\top x + y = -\beta,\ y \ge 0.
\end{aligned} max c ⊤ x a ⊤ x = β , β < 0 a ⊤ x ≤ β , β ≥ 0 a ⊤ x ≤ β , β < 0 a ⊤ x ≥ β , β > 0 a ⊤ x ≥ β , β ≤ 0 = − min ( − c ) ⊤ x , ⟺ ( − a ) ⊤ x = − β , ⟺ a ⊤ x + y = β , y ≥ 0 , ⟺ ( − a ) ⊤ x − y = − β , y ≥ 0 , ⟺ a ⊤ x − y = β , y ≥ 0 , ⟺ ( − a ) ⊤ x + y = − β , y ≥ 0. 每个等价关系成立，是因为 y y y a ⊤ x ≤ β a^\top x \le \beta a ⊤ x ≤ β y = β − a ⊤ x y = \beta - a^\top x y = β − a ⊤ x − 1 -1 − 1 − 1 -1 − 1 β < 0 \beta < 0 β < 0 β ≤ 0 \beta \le 0 β ≤ 0
标准形假定每个 变量都满足 x ≥ 0 x \ge 0 x ≥ 0 Variable 中保存的边界不会被使用。
基与基本解
设 A A A m m m 基 是由 m m m B B B A B A_B A B N N N 基本解
x N = 0 , x B = A B − 1 b , x_N = 0, \qquad x_B = A_B^{-1} b, x N = 0 , x B = A B − 1 b , 当 x B ≥ 0 x_B \ge 0 x B ≥ 0 可行 的。线性规划基本定理指出：若规划有最优解，则它有最优的基本可行解，因此单纯形法只需访问各个基。1 1 参见例如 Chvátal，Linear Programming （1983）第 3 章，或 Bertsimas 与 Tsitsiklis，Introduction to Linear Optimization （1997）定理 2.7。
单纯形表与检验数
求解器把规划存为一个 ( m + 1 ) × ( n + 1 ) (m + 1) \times (n + 1) ( m + 1 ) × ( n + 1 )
T = [ c ⊤ 0 A b ] . T = \begin{bmatrix} c^\top & 0 \\ A & b \end{bmatrix}. T = [ c ⊤ A 0 b ] . 第 0 行是目标行，第 1 , … , m 1, \dots, m 1 , … , m T r s ≠ 0 T_{rs} \ne 0 T r s = 0 ( r , s ) (r, s) ( r , s ) 转轴 是一步 Gauss–Jordan 消元：将第 r r r T r s T_{rs} T r s i i i r r r T i s T_{is} T i s s s s e r e_r e r Lp::pivot。
假设一系列转轴已经把基 B B B
T B = [ c ⊤ − y ⊤ A − y ⊤ b A B − 1 A A B − 1 b ] T_B = \begin{bmatrix} c^\top - y^\top A & -y^\top b \\ A_B^{-1} A & A_B^{-1} b \end{bmatrix} T B = [ c ⊤ − y ⊤ A A B − 1 A − y ⊤ b A B − 1 b ] 其中 y ∈ R m y \in \mathbb{R}^m y ∈ R m c B ⊤ − y ⊤ A B = 0 c_B^\top - y^\top A_B = 0 c B ⊤ − y ⊤ A B = 0
y ⊤ = c B ⊤ A B − 1 . y^\top = c_B^\top A_B^{-1}. y ⊤ = c B ⊤ A B − 1 . 代入 y y y
T B = [ c ˉ ⊤ − z B A ˉ b ˉ ] , c ˉ ⊤ = c ⊤ − c B ⊤ A B − 1 A , z B = c B ⊤ A B − 1 b , A ˉ = A B − 1 A , b ˉ = A B − 1 b . T_B = \begin{bmatrix} \bar c^\top & -z_B \\ \bar A & \bar b \end{bmatrix},
\qquad
\bar c^\top = c^\top - c_B^\top A_B^{-1} A, \quad
z_B = c_B^\top A_B^{-1} b, \quad
\bar A = A_B^{-1} A, \quad
\bar b = A_B^{-1} b. T B = [ c ˉ ⊤ A ˉ − z B b ˉ ] , c ˉ ⊤ = c ⊤ − c B ⊤ A B − 1 A , z B = c B ⊤ A B − 1 b , A ˉ = A B − 1 A , b ˉ = A B − 1 b . 其中 c ˉ j \bar c_j c ˉ j 检验数 （约简成本），b ˉ = x B \bar b = x_B b ˉ = x B 取负 的目标值 − z B = − c ⊤ x -z_B = -c^\top x − z B = − c ⊤ x Lp::phase_2 返回 − z -z − z Lp::two_stage 对最小化目标再次取负的原因。
检验数描述目标沿约束如何变化。对任意满足 A x = b Ax = b A x = b x x x x B = b ˉ − A B − 1 A N x N x_B = \bar b - A_B^{-1} A_N x_N x B = b ˉ − A B − 1 A N x N
c ⊤ x = c B ⊤ x B + c N ⊤ x N = c B ⊤ b ˉ − c B ⊤ A B − 1 A N x N + c N ⊤ x N = z B + c ˉ N ⊤ x N . \begin{aligned}
c^\top x &= c_B^\top x_B + c_N^\top x_N \\
&= c_B^\top \bar b - c_B^\top A_B^{-1} A_N x_N + c_N^\top x_N \\
&= z_B + \bar c_N^\top x_N .
\end{aligned} c ⊤ x = c B ⊤ x B + c N ⊤ x N = c B ⊤ b ˉ − c B ⊤ A B − 1 A N x N + c N ⊤ x N = z B + c ˉ N ⊤ x N . 最优性
若 c ˉ ≥ 0 \bar c \ge 0 c ˉ ≥ 0 每个可行的 x x x x N ≥ 0 x_N \ge 0 x N ≥ 0 c ⊤ x = z B + c ˉ N ⊤ x N ≥ z B c^\top x = z_B + \bar c_N^\top x_N \ge z_B c ⊤ x = z B + c ˉ N ⊤ x N ≥ z B z B z_B z B Lp::simplex_iteration 停止，这正是该判据（c ˉ B = 0 \bar c_B = 0 c ˉ B = 0
选择进基列
若某个 c ˉ s < 0 \bar c_s < 0 c ˉ s < 0 x s x_s x s 0 0 0 ∣ c ˉ s ∣ |\bar c_s| ∣ c ˉ s ∣
s = arg min j c ˉ j , c ˉ s < 0 , s = \operatorname*{arg\,min}_{j} \; \bar c_j, \qquad \bar c_s < 0, s = j arg min c ˉ j , c ˉ s < 0 , 即 Dantzig 规则，相同时取最左列。
比值检验与无界性
将 x s x_s x s t ≥ 0 t \ge 0 t ≥ 0 x B ( t ) = b ˉ − t a ˉ s x_B(t) = \bar b - t\, \bar a_s x B ( t ) = b ˉ − t a ˉ s a ˉ s \bar a_s a ˉ s A ˉ \bar A A ˉ s s s z B + t c ˉ s z_B + t\, \bar c_s z B + t c ˉ s
若 a ˉ s ≤ 0 \bar a_s \le 0 a ˉ s ≤ 0 x B ( t ) x_B(t) x B ( t ) t t t t ≥ 0 t \ge 0 t ≥ 0 x ( t ) x(t) x ( t ) z B + t c ˉ s → − ∞ z_B + t\,\bar c_s \to -\infty z B + t c ˉ s → − ∞ Lp::simplex_iteration 以 Problem is unbounded abort。
否则，可行性要求对每个 a ˉ i s > 0 \bar a_{is} > 0 a ˉ i s > 0 b ˉ i − t a ˉ i s ≥ 0 \bar b_i - t\, \bar a_{is} \ge 0 b ˉ i − t a ˉ i s ≥ 0
t ∗ = min i : a ˉ i s > 0 b ˉ i a ˉ i s , t^\ast = \min_{i \,:\, \bar a_{is} > 0} \frac{\bar b_i}{\bar a_{is}}, t ∗ = i : a ˉ i s > 0 min a ˉ i s b ˉ i , 即最小比值检验 。取到最小值的行 r r r t ∗ t^\ast t ∗ ( r , s ) (r, s) ( r , s ) B ′ = B ∖ { B r } ∪ { s } B' = B \setminus \{B_r\} \cup \{s\} B ′ = B ∖ { B r } ∪ { s } x ( t ∗ ) ≥ 0 x(t^\ast) \ge 0 x ( t ∗ ) ≥ 0
z B ′ = z B + t ∗ c ˉ s ≤ z B . z_{B'} = z_B + t^\ast\, \bar c_s \le z_B . z B ′ = z B + t ∗ c ˉ s ≤ z B . 相同时求解器取最上面的行。
终止性与退化
若每次转轴都有 t ∗ > 0 t^\ast > 0 t ∗ > 0 非退化 ），目标严格下降，基不会重复，方法至多经过 ( n m ) \binom{n}{m} ( m n ) b ˉ r = 0 \bar b_r = 0 b ˉ r = 0 t ∗ = 0 t^\ast = 0 t ∗ = 0 2 2 E. M. L. Beale，“Cycling in the dual simplex algorithm”，Naval Research Logistics Quarterly 2（1955），给出了一个含三个约束的循环例子。R. G. Bland，“New finite pivoting rules for the simplex method”，Mathematics of Operations Research 2（1977），证明了在候选者中总选择下标最小的进基变量和出基变量可以防止循环。 max_iterations 次转轴（默认 1000）后停止并打印消息，此时返回的单纯形表不是最优的。
设计决策
带人工变量的两阶段法
问题。 单纯形法需要从一个基本可行解出发，而标准形规划并不自带这样的解。
可选方案。 (a) 大 M 法：在原目标中加入代价为大数 M M M
选择。 (b)。大 M 法需要一个数值 M M M
第一阶段。 对每个没有可用单位列的约束行 i i i i i i 1 1 1 a i ≥ 0 a_i \ge 0 a i ≥ 0 R R R
min w = ∑ i ∈ R a i subject to A x + ∑ i ∈ R a i e i = b , x ≥ 0 , a ≥ 0. \min\; w = \sum_{i \in R} a_i \quad \text{subject to} \quad A x + \sum_{i \in R} a_i e_i = b, \quad x \ge 0,\ a \ge 0 . min w = i ∈ R ∑ a i subject to A x + i ∈ R ∑ a i e i = b , x ≥ 0 , a ≥ 0. 其初始基由人工列和已有的单位列组成，基本解为 x B = b ≥ 0 x_B = b \ge 0 x B = b ≥ 0 to_standard 已使 b b b
原规划可行当且仅当第一阶段的最优值 w ∗ = 0 w^\ast = 0 w ∗ = 0 若 x x x A x = b Ax = b A x = b ( x , a = 0 ) (x, a = 0) ( x , a = 0 ) w = 0 w = 0 w = 0 w ≥ 0 w \ge 0 w ≥ 0 w ∗ = 0 w^\ast = 0 w ∗ = 0 w ∗ = 0 w^\ast = 0 w ∗ = 0 a ≥ 0 a \ge 0 a ≥ 0 a = 0 a = 0 a = 0 x x x A x = b Ax = b A x = b
第一阶段单纯形表的典范形式。 第一阶段目标的第 0 行起初为 [ 0 ∣ 1 ⊤ ∣ 0 ] [\,0 \mid \mathbf 1^\top \mid 0\,] [ 0 ∣ 1 ⊤ ∣ 0 ] 1 1 1 A B = I A_B = I A B = I y ⊤ = c B ⊤ A B − 1 y^\top = c_B^\top A_B^{-1} y ⊤ = c B ⊤ A B − 1 y = 1 R y = \mathbf 1_R y = 1 R
[ − ∑ i ∈ R A i ⋅ ∣ 0 ∣ − ∑ i ∈ R b i ] , \Bigl[\; -\sum_{i \in R} A_{i\cdot} \;\Big|\; 0 \;\Big|\; -\sum_{i \in R} b_i \;\Bigr], [ − i ∈ R ∑ A i ⋅ 0 − i ∈ R ∑ b i ] , 即从第 0 行中减去每个人工行所得。Lp::phase_1 在调用 simplex_iteration 之前恰好执行这些减法。
第二阶段。 若 w ∗ ≠ 0 w^\ast \ne 0 w ∗ = 0 Lp::phase_2 会 abort：规划不可行。否则它删去人工列，把原代价 c c c j j j c j c_j c j simplex_iteration。
复用单位列作为初始基
对于 b ≥ 0 b \ge 0 b ≥ 0 <= 行，其松弛变量的列等于单位向量，因此无需人工变量即可作为初始基变量。构造第一阶段单纯形表的私有辅助函数在某列有一项等于 1 1 1 0 0 0 <= 且右端项非负的规划完全不需要人工变量，第一阶段实际上被跳过：其单纯形表保留原目标，所以第一阶段已经在优化它。
稠密单纯形表
单纯形表是稠密的 @mutable.Matrix。每次转轴需要 ( m + 1 ) ( n + 1 ) (m + 1)(n + 1) ( m + 1 ) ( n + 1 ) O ( m n ) O(mn) O ( mn )
以 abort 代替返回错误
不可行和无界的规划、零主元、未知的优化方向以及未知的关系字符串都会以 abort 终止程序。本包早于 Luna-Flow 返回带结构化错误的 Result 的约定。调用者无法从这些情况中恢复；参见下文的边界。
数值容差
单纯形迭代把各项与零做精确比较：由舍入产生的 − 10 − 17 -10^{-17} − 1 0 − 17 ApproximatelyZero：判定 w ∗ w^\ast w ∗ Double，阈值为 ∣ x ∣ < 10 − 15 |x| < 10^{-15} ∣ x ∣ < 1 0 − 15
转轴中的舍入误差与所涉及元素的大小成比例。浮点运算中的一次更新 t i j − t i s t r j t_{ij} - t_{is}\, t_{rj} t ij − t i s t r j
fl ( t i j − t i s t r j ) = ( t i j − t i s t r j ( 1 + δ 1 ) ) ( 1 + δ 2 ) , ∣ δ 1 ∣ , ∣ δ 2 ∣ ≤ u = 2 − 53 ≈ 1.11 × 10 − 16 , \operatorname{fl}\bigl(t_{ij} - t_{is}\, t_{rj}\bigr) = \bigl(t_{ij} - t_{is}\, t_{rj}(1 + \delta_1)\bigr)(1 + \delta_2), \qquad |\delta_1|, |\delta_2| \le u = 2^{-53} \approx 1.11 \times 10^{-16}, fl ( t ij − t i s t r j ) = ( t ij − t i s t r j ( 1 + δ 1 ) ) ( 1 + δ 2 ) , ∣ δ 1 ∣ , ∣ δ 2 ∣ ≤ u = 2 − 53 ≈ 1.11 × 1 0 − 16 , 因此其绝对误差大约以 u ( ∣ t i j ∣ + 2 ∣ t i s t r j ∣ ) u\,(|t_{ij}| + 2|t_{is}\, t_{rj}|) u ( ∣ t ij ∣ + 2∣ t i s t r j ∣ ) M M M k k k w ∗ w^\ast w ∗ k u M k\,u\,M k u M 10 − 15 ≈ 9 u 10^{-15} \approx 9u 1 0 − 15 ≈ 9 u 1 1 1 w ∗ = − 4.44 × 10 − 16 = − 4 u w^\ast = -4.44 \times 10^{-16} = -4u w ∗ = − 4.44 × 1 0 − 16 = − 4 u b b b
正确性与不变量
求解器在转轴之间维持以下不变量：
典范形式。 约束行中的每个基列都是单位向量，且第 0 行在每个基列上为零。Lp::pivot 由构造保持这一点，Lp::phase_1 为人工基建立它，Lp::phase_2 在替换第 0 行后重新建立它。原始可行性。 右端项 b ˉ \bar b b ˉ to_standard 使 b ≥ 0 b \ge 0 b ≥ 0 第 0 行中的目标。 第 0 行的最后一项是 − z B -z_B − z B
在这些不变量下，停止判据如上所推导是正确的：没有负检验数意味着最优，而负检验数配上非正列意味着无界。
实现中存在一些已知缺陷，会在特定情形下破坏这些保证。这里列出它们以便用户规避；代码本身未作修改。
目标行对齐。 Lp 用 Obj_func::to_vector 构造第 0 行，它按变量名顺序取已存储的系数并跳过零系数。因此只要某个目标系数为零，或变量未按名称顺序声明，第 0 行就与列顺序不一致，求解器优化的是被打乱的目标。例如，在 x 1 ≥ 1 x_1 \ge 1 x 1 ≥ 1 x 2 ≥ 2 x_2 \ge 2 x 2 ≥ 2 x 2 x_2 x 2 1 1 1 2 2 2 基的识别。 第二阶段把单位列识别为基变量，并在每行取第一个这样的列。当两列是相同的单位向量时，即使目标值正确，报告的解也可能指错变量。退化的人工变量。 若第一阶段结束后某个人工变量仍以零值留在基中，第二阶段会直接删去其列而不先将其转出，该行在第二阶段中就没有基列。人工变量索引。 对没有人工变量的行，Lp::phase_1 收到的人工列为 0，它无法与第 0 列区分。若此时第 0 列的第一阶段目标系数恰为 1 1 1
被否决的方案
大 M 法。 出于“带人工变量的两阶段法”中所述的数值原因而被否决。Bland 规则。 它保证终止，但在非退化规划上通常比 Dantzig 规则需要多得多的转轴。本包保留 Dantzig 规则，转而限制迭代次数。修正单纯形法与内点法。 它们在大型稀疏规划上扩展性更好，但会隐藏单纯形表，而本包有意公开单纯形表。通过代换实施变量边界。 边界 l ≤ x ≤ u l \le x \le u l ≤ x ≤ u x = l + x ′ x = l + x' x = l + x ′ x ′ ≤ u − l x' \le u - l x ′ ≤ u − l
边界
本包刻意不做以下事情：
以值的形式返回错误。不可行和无界的规划、零主元以及非法参数都会使程序 abort；
实施 x ≥ 0 x \ge 0 x ≥ 0
在每次单纯形运行 1000 次迭代的上限之外防止循环；
利用稀疏性、热启动或基分解；
计算对偶值或灵敏度区间，尽管最终基的向量 y = ( A B − 1 ) ⊤ c B y = (A_B^{-1})^\top c_B y = ( A B − 1 ) ⊤ c B
求解整数、混合整数或非线性规划。
