core 教程

本教程介绍如何编写线性规划、将其化为标准形，并用 linear-program 的两阶段单纯形法求解。它从一个两变量的生产问题开始，最后在单纯形表上手动执行单纯形法。每一步背后的数学见设计说明。

快速入门

该模块名为 linear-program，没有 Luna-Flow/ 命名空间，目前也尚未以 Luna-Flow 组织的名义发布到 mooncakes。请从本地检出使用：把该检出及其依赖 luna-generic 与 linear-algebra 的检出加入 moon.work 的 members，并在 moon.mod 的 import 块中列出它：

import {
  "linear-program@0.1.0",
  "Luna-Flow/linear-algebra@0.4.7",
}

然后在 moon.pkg 中导入该包：

import {
  "linear-program" @lp,
  "Luna-Flow/linear-algebra/mutable" @la,
  "moonbitlang/core/double",
}

最小的实用程序在 $2x_1 + x_2 \le 12$、$x_1 + 3x_2 \le 10$ 和 $x_1, x_2 \ge 0$ 的约束下最大化 $3x_1 + 4x_2$：

using @lp {type Variable, type Poly, type Obj_func, type Lp}

test "quick start" {
  let x = [Variable::new("x1", 0.0, 100.0), Variable::new("x2", 0.0, 100.0)]
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("max", [3.0, 4.0], x), x)
    .cons_from_Array(ineq_array=[([2.0, 1.0], "<=", 12.0), ([1.0, 3.0], "<=", 10.0)])
  let (solution, value) = lp.to_standard().two_stage()
  debug_inspect(solution, content="[5.199999999999999, 1.6000000000000005, 0, 0]")
  inspect(-value, content="22")
}

最优解为 $x_1 = 5.2$、$x_2 = 1.6$，最优值为 $22$（忽略舍入误差）。解有四个分量，是因为 to_standard 为每个不等式添加了一个松弛变量；值取了负号，是因为标准形最小化的是 $-3x_1 - 4x_2$。

日常任务

编写规划并读回

变量带有名称和边界；目标函数带有优化方向和系数；约束是系数行。Show 以常见记法打印规划：

fn nonneg(names : Array[String]) -> Array[Variable] {
  names.map(n => Variable::new(n, 0.0, @double.max_value))
}

test "write a program" {
  let x = nonneg(["x1", "x2", "x3"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("max", [1.0, 2.0, 3.0], x), x)
    .cons_from_Array(
      eq_array=[([2.0, 2.0, 3.0], 100.0)],
      ineq_array=[([2.5, 0.0, 4.7], ">=", 0.0), ([-1.7, 2.5, 12.0], "<=", 0.0)],
    )
  inspect(
    lp,
    content=(
      #|max  x1 + 2x2 + 3x3
      #|s.t. 2x1 + 2x2 + 3x3 = 100
      #|     2.5x1 + 4.7x3 >= 0
      #|    -1.7x1 + 2.5x2 + 12x3 <= 0
      #|     x1 : [0, 1.7976931348623157e+308] x2 : [0, 1.7976931348623157e+308] x3 : [0, 1.7976931348623157e+308]
      #|
    ),
  )
}

等式行为 (coefficients, b)；不等式行为 (coefficients, relation, b)，关系必须恰好写作 ">=" 或 "<="。后续示例都会用到辅助函数 nonneg。

将规划化为标准形

求解器处理的是在 $Ax = b$、$b \ge 0$、$x \ge 0$ 约束下的 $\min c^\top x$。to_standard 生成这种形式，访问函数返回 $c$、$A$ 和 $b$：

test "standard form" {
  let x = nonneg(["x1", "x2"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("min", [3.0, 4.0], x), x)
    .cons_from_Array(ineq_array=[([2.0, 1.0], ">=", 12.0), ([1.0, 3.0], ">=", 10.0)])
  let s = lp.to_standard()
  inspect(
    s,
    content=(
      #|min  3x1 + 4x2
      #|s.t. 2x1 + x2 + -1y1 = 12
      #|     x1 + 3x2 + -1y2 = 10
      #|      x1 : [0, 1.7976931348623157e+308] x2 : [0, 1.7976931348623157e+308] y1 : [0, 1.7976931348623157e+308] y2 : [0, 1.7976931348623157e+308]
      #|
    ),
  )
  debug_inspect(s.get_objfunc_vector(), content="[3, 4, 0, 0]")
  debug_inspect(s.get_b_vector(), content="[12, 10]")
}

每个 >= 行都获得了一个系数为 $-1$ 的剩余变量，因此 $2x_1 + x_2 \ge 12$ 变成了 $2x_1 + x_2 - y_1 = 12$，且 $y_1 \ge 0$。

求解最小化问题

最小化问题不需要变号：返回值就是最小值。

test "minimise" {
  let x = nonneg(["x1", "x2"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("min", [3.0, 4.0], x), x)
    .cons_from_Array(ineq_array=[([2.0, 1.0], ">=", 12.0), ([1.0, 3.0], ">=", 10.0)])
  let (solution, value) = lp.to_standard().two_stage()
  debug_inspect(solution, content="[5.199999999999999, 1.6000000000000005, 0, 0]")
  inspect(value, content="22")
}

这里没有任何松弛列是单位列，因此第一阶段添加两个人工变量并将其驱动到零，随后第二阶段从第一阶段找到的可行点出发最小化真正的目标。

用矩阵编写含等式约束的规划

当所有约束都是等式时，整个规划可以放进一个矩阵：第 0 行存放目标系数，其后每行存放一个约束 $[a_i \mid b_i]$。下面的规划在 $x_1 + 2x_2 + 3x_3 = 10$ 与 $2x_1 + x_2 + x_3 = 8$ 的约束下最小化 $-x_1 - 2x_2 - 3x_3$：

test "from a matrix" {
  let x = nonneg(["x1", "x2", "x3"])
  let m : @la.Matrix[Double] = @la.Matrix::from_2d_array([
    [-1.0, -2.0, -3.0, 0.0],
    [1.0, 2.0, 3.0, 10.0],
    [2.0, 1.0, 1.0, 8.0],
  ])
  let lp = Lp::from_matrix(m, x, "min")
  inspect(
    lp,
    content=(
      #|min  -1x1 + -2x2 + -3x3
      #|s.t. x1 + 2x2 + 3x3 = 10
      #|     2x1 + x2 + x3 = 8
      #|      x1 : [0, 1.7976931348623157e+308] x2 : [0, 1.7976931348623157e+308] x3 : [0, 1.7976931348623157e+308]
      #|
    ),
  )
  let (solution, value) = lp.to_standard().two_stage()
  debug_inspect(solution, content="[1.9999999999999991, 4.000000000000001, 0]")
  inspect(value, content="-10")
}

最优解为 $x = (2, 4, 0)$，最优值为 $-10$。from_matrix 会忽略第 0 行的最后一项。

逐项构造表达式

Poly 与 Obj_func 也可以逐项构造，这在系数来自数据时很方便：

test "term by term" {
  let x = nonneg(["x1", "x2"])
  let p : Poly[Double] = Poly::new()
  p.add_term_inplace(x[0], 3.0)
  p.add_term_inplace(x[1], 4.0)
  p.add_term_inplace(x[0], 1.0)
  let obj = Obj_func::new("max").set_poly(p)
  inspect(obj, content="max  4x1 + 4x2")
  let lp = Lp::new(obj, x).cons_from_Array(ineq_array=[([1.0, 1.0], "<=", 5.0)])
  let (_, value) = lp.to_standard().two_stage()
  inspect(-value, content="20")
}

add_term_inplace 会累加到已有系数上，因此 x1 最终为 $3 + 1 = 4$。

进阶

手动执行单纯形法

two_stage 的各个步骤都是公开的。当规划已经处于典范形式时（例如所有约束都是 $b \ge 0$ 的 <=，且松弛列构成单位矩阵），你可以自己写出单纯形表并在其上转轴。第 0 行存放检验数和 $-z$：

test "pivot by hand" {
  // max 3x1 + 4x2  ==  min -3x1 - 4x2, slack variables y1, y2 basic
  let t : @la.Matrix[Double] = @la.Matrix::from_2d_array([
    [-3.0, -4.0, 0.0, 0.0, 0.0],
    [2.0, 1.0, 1.0, 0.0, 12.0],
    [1.0, 3.0, 0.0, 1.0, 10.0],
  ])
  // x2 has the most negative reduced cost; the ratios are 12/1 and 10/3, so row 2 leaves
  Lp::pivot(t, 2, 1)
  inspect(
    t,
    content=(
      #||-1.6666666666666667, 0, 0, 1.3333333333333333, 13.333333333333334|
      #||1.6666666666666667, 0, 1, -0.3333333333333333, 8.666666666666666|
      #||0.3333333333333333, 1, 0, 0.3333333333333333, 3.3333333333333335|
    ),
  )
  // finish with Dantzig's rule
  Lp::simplex_iteration(t)
  inspect(t[0][4], content="22")
}

第一次转轴之后，目标从 $0$ 改进到 $40/3$；$x_1$ 仍有负的检验数，说明再转轴一次可以进一步改进到 $22$。向 simplex_iteration 传入 debug=true 可以打印每次转轴。

使用其他系数类型

每个函数都在系数类型 V 上泛型。求解器需要 Zero、One、Compare、算术运算符、Show 和 ApproximatelyZero。要使用自己的精确数类型，请为其实现这些 trait，包括容差 trait；对精确类型而言，它就是与零的普通相等比较：

impl @lp.ApproximatelyZero for Rational with fn is_zero_eps(x) {
  x == Rational::zero()
}

使用精确算术时，求解器返回精确结果，容差也就无关紧要。

自行检查结果

two_stage 遇到不可行和无界的规划时会 abort 而不是返回错误，因此当规划来自用户数据时，请在调用前验证输入。求解之后，把解代回标准形的约束以确认可行性：

test "check feasibility" {
  let x = nonneg(["x1", "x2"])
  let s : Lp[Double] = Lp::new(Obj_func::from_array("max", [3.0, 4.0], x), x)
    .cons_from_Array(ineq_array=[([2.0, 1.0], "<=", 12.0), ([1.0, 3.0], "<=", 10.0)])
    .to_standard()
  let (solution, _) = s.two_stage()
  let a = s.get_coeff_matrix()
  let b = s.get_b_vector()
  for i in 0..<b.length() {
    let mut lhs = 0.0
    for j in 0..<solution.length() {
      lhs = lhs + a[i][j] * solution[j]
    }
    assert_true((lhs - b[i]).abs() < 1.0e-9)
  }
}

常见陷阱

• 返回值属于标准形。 to_standard 把每个规划都变成最小化问题，因此对于最大化问题，请对 two_stage 的返回值取负。
• 解包含松弛变量。 其长度是标准形的变量个数；你自己的变量按声明顺序排在最前面。
• 目标系数必须对齐。 由于一个已知缺陷，目标行是按变量名顺序由非零系数填充的。请给每个变量一个非零目标系数，并按名称顺序声明变量（x1、x2……，不要把 x2 放在 x10 之前），否则求解器会优化另一个目标。
• 边界不会被实施。 Variable::new("x", 0.0, 100.0) 保存了上界 $100$，但求解器只假定 $x \ge 0$。如需上界，请显式添加约束 ([1.0, ...], "<=", 100.0)。
• 不可行和无界的规划会 abort。 没有可供检查的 Result。
• 使用 Double 系数。 Int 满足 trait 约束，但整数除法会截断，因此在 Int 上转轴会得到错误结果。
• 容差是绝对的。 第二阶段把 $|x| < 10^{-15}$ 视为零。对于系数较大的规划，舍入误差可能超过该阈值，可行的规划可能被报告为不可行；请缩放各行，使系数大小适中。
• 不要把变量命名为 y1、y2…… to_standard 用这些名称表示松弛变量，而变量是按名称识别的。

下一步

• API 参考记录了每个类型和函数，包括 Lp::phase_1 与 Lp::phase_2。
• 设计说明推导了标准形、最优性与无界性判据以及两个阶段。
• linear-algebra 提供了用于单纯形表的 Matrix 类型，luna-generic 提供了约束系数类型的 trait。