core API

src 下的包就是 linear-program 的全部公开接口。其接口文件为 src/pkg.generated.mbti。包路径是 linear-program：模块名目前还没有 Luna-Flow/ 命名空间，因此其他包以 "linear-program" 导入它。

系数类型 V 是泛型的。每个函数都列出了它所需的 trait；Double 满足全部要求，也是求解器唯一经过测试的类型。矩阵是 linear-algebra 中 mutable 包的 @mutable.Matrix[V]。规划矩阵或单纯形表在第 0 行存放目标函数，其后每行存放一个等式约束；最后一列存放右端项。

本页的示例都是位于如下 moon.pkg 所在包中的测试：

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

并且共享以下声明：

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

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

变量

`Variable`

Variable 是带有下界和上界的具名决策变量。

type Variable
pub impl @luna-generic.One for Variable
pub impl Compare for Variable
pub impl Eq for Variable
pub impl Show for Variable

Variable 是抽象类型：只能用 Variable::new 创建值，不能直接读取其字段。两个变量名称相同即相等，与边界无关，并按名称排序。边界会被保存和打印，但求解器不会读取它们：规划中的每个变量都被视为非负。

`Variable::new`

Variable::new 由名称、下界和上界创建一个变量。

pub fn Variable::new(String, Double, Double) -> Self

若上界不应起任何限制作用，请使用 @double.max_value。Lp::to_standard 会把它新增的变量命名为 y1、y2……，因此不要把这些名称用于自己的变量。

`Variable::show_all`

show_all 打印变量名及其边界。

pub fn Variable::show_all(Self) -> String

test "variables" {
  let x = Variable::new("x1", 0.0, 100.0)
  inspect(x, content="x1")
  inspect(x.show_all(), content="x1 : [0, 100]")
  assert_true(x == Variable::new("x1", -5.0, 5.0))
  assert_true(x < Variable::new("x2", 0.0, 1.0))
}

`Variable::to_string`, `Variable::equal`, `Variable::compare`

它们是 Show、Eq 与 Compare 的方法，被提升为 Variable 的方法。

pub fn Variable::to_string(Self) -> String
pub fn Variable::equal(Self, Self) -> Bool
pub fn Variable::compare(Self, Self) -> Int

to_string 返回变量名。equal 与 compare 只看名称；compare 是名称的字典序，因此 "x10" 排在 "x2" 之前。新代码中请优先使用 "\{x}"、== 与 <。

`Variable::one`

Variable::one 是表示线性表达式常数项的无名变量。

pub fn Variable::one() -> Self

它是 luna-generic 的 One 方法提升而来的方法。其名称为空，边界为 @double.min_value 与 @double.max_value。Poly::from_var_array 会把没有对应变量的系数挂到它上面。

线性表达式

`Poly`

Poly[V] 是线性表达式：从变量到非零系数的映射，按变量名称排序保存。

type Poly[V] derive(Eq)
pub impl[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Neg for Poly[V]
pub impl[V : Eq + Show + @luna-generic.Semiring + @luna-generic.One] Show for Poly[V]

Poly 是可变的：add_term_inplace 会修改它，所有共享它的值都能看到变化。需要独立的表达式时请使用 copy。两个表达式当且仅当变量相同且系数相等时相等。

`Poly::new`, `Poly::from_var_array`

Poly::new 创建空表达式；Poly::from_var_array 由系数创建表达式。

pub fn[V] Poly::new() -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Poly::from_var_array(Array[V], Array[Variable]) -> Self[V]

Poly::from_var_array(coeffs, vars) 将 coeffs[i] 与 vars[i] 配对，并跳过零系数。位置超出最后一个变量的系数挂到常数项 Variable::one() 上；若有多个这样的系数，以最后一个为准。

`Poly::add_term_inplace`

add_term_inplace 原地向表达式加上某个变量的倍数。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Poly::add_term_inplace(Self[V], Variable, V) -> Unit

p.add_term_inplace(x, c) 把 c 加到 x 的系数上。加零不做任何事；系数变为零的项会被删除，因此永远不会存储零系数。

`Poly::copy`

copy 返回一个项相同、但不与原表达式共享存储的表达式。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Poly::copy(Self[V]) -> Self[V]

`Poly::to_array`, `Poly::to_vector`

to_array 列出已存储的系数；to_vector 为每个给定变量列出一个系数。

pub fn[V] Poly::to_array(Self[V]) -> Array[V]
pub fn[V : @luna-generic.Zero] Poly::to_vector(Self[V], Array[Variable]) -> Array[V]

to_array 按变量名顺序返回非零系数，不含变量本身。to_vector(vars) 按 vars 的顺序返回 vars.length() 个系数，未出现的变量取零；这正是规划矩阵中使用的稠密行。

`Poly::neg`, `Poly::equal`, `Poly::to_string`

它们是 Neg、Eq 与 Show 的方法，被提升为 Poly 的方法。

pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Poly::neg(Self[V]) -> Self[V]
pub fn[V : Eq] Poly::equal(Self[V], Self[V]) -> Bool
pub fn[V : Eq + Show + @luna-generic.Semiring + @luna-generic.One] Poly::to_string(Self[V]) -> String

-p 对每个系数取负并返回新表达式。Show 按名称顺序打印各项，以 " + " 连接，并省略等于一的系数，例如 3.5x1 + x2 + -2x3。

test "linear expressions" {
  let vars = nonneg(["x1", "x2", "x3"])
  let p : Poly[Double] = Poly::from_var_array([3.5, 0.0, -2.0], vars)
  inspect(p, content="3.5x1 + -2x3")
  p.add_term_inplace(vars[1], 1.0)
  p.add_term_inplace(vars[2], 2.0)
  inspect(p, content="3.5x1 + x2")
  inspect(-p, content="-3.5x1 + -1x2")
  debug_inspect(p.to_array(), content="[3.5, 1]")
  debug_inspect(p.to_vector(vars), content="[3.5, 1, 0]")
  assert_true(p.copy() == p)
}

目标函数

`Obj_func`

Obj_func[V] 是目标函数：一个线性表达式加上优化方向（最小化或最大化）。

type Obj_func[V]
pub impl[V : Eq + Show + @luna-generic.Semiring] Show for Obj_func[V]

Show 打印优化方向和表达式，例如 max 3x1 + 4x2。Obj_func::to_string 是提升而来的 Show 方法。

pub fn[V : Eq + Show + @luna-generic.Semiring] Obj_func::to_string(Self[V]) -> String

`Obj_func::new`, `Obj_func::from_array`

Obj_func::new 创建表达式为空的目标函数；Obj_func::from_array 还会同时设置表达式。

pub fn[V] Obj_func::new(String) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Obj_func::from_array(String, Array[V], Array[Variable]) -> Self[V]

优化方向为 "min"、"Min"、"MIN"、"minimize"、"Minimize" 之一时表示最小化，为 "max"、"Max"、"MAX"、"maximize"、"Maximize" 之一时表示最大化。其他字符串会导致 abort。Obj_func::from_array(sense, coeffs, vars) 用 Poly::from_var_array 构造表达式。

`Obj_func::set_poly`, `Obj_func::add_term_inplace`

set_poly 替换表达式；add_term_inplace 原地向其添加一项。

pub fn[V] Obj_func::set_poly(Self[V], Poly[V]) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero] Obj_func::add_term_inplace(Self[V], Variable, V) -> Unit

set_poly(p) 返回一个方向相同、表达式为 p 的新目标函数；它不会复制 p。

`Obj_func::judge_max`, `Obj_func::to_min`

judge_max 判断目标函数是否为最大化；to_min 返回等价的最小化目标函数。

pub fn[V] Obj_func::judge_max(Self[V]) -> Bool
pub fn[V : Eq + @luna-generic.Semiring + @luna-generic.Zero + Neg] Obj_func::to_min(Self[V]) -> Self[V]

to_min 对最小化目标原样返回，对最大化目标则将表达式取负，因为 $\max c^\top x = -\min\,(-c)^\top x$ 。

`Obj_func::to_vector`

to_vector 返回规划矩阵的目标行。

pub fn[V : @luna-generic.Zero] Obj_func::to_vector(Self[V], Int) -> Array[V]

to_vector(n) 取 to_array 的系数，并用零补齐到长度 n + 1，即单纯形表一行的宽度。由于 to_array 只按名称顺序列出非零系数，只有当每个变量的系数都非零且变量按名称顺序声明时，结果才与变量对齐。参见下文 Lp::two_stage 中的警告。

test "objectives" {
  let vars = nonneg(["x1", "x2"])
  let obj : Obj_func[Double] = Obj_func::from_array("max", [3.0, 4.0], vars)
  inspect(obj, content="max  3x1 + 4x2")
  assert_true(obj.judge_max())
  inspect(obj.to_min(), content="min  -3x1 + -4x2")
  debug_inspect(obj.to_vector(2), content="[3, 4, 0]")
}

约束

`Constraint`

Constraint[V] 是线性规划的约束集合：等式 poly = b 以及不等式 poly <= b 或 poly >= b。

type Constraint[V]
pub impl[V : Show + Eq + @luna-generic.Semiring] Show for Constraint[V]
pub fn[V : Show + Eq + @luna-generic.Semiring] Constraint::to_string(Self[V]) -> String

Show 在 s.t. 之后打印约束集，每行一个约束，等式在前。规划通过 Lp::cons_from_Array 获得约束；单独的 Constraint 值适合用来构造和打印约束集。

`Constraint::new`, `Constraint::from_array`

Constraint::new 创建空集合；Constraint::from_array 由系数行构造约束集。

pub fn[V] Constraint::new() -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Constraint::from_array(eq_array? : Array[(Array[V], V)], ineq_array? : Array[(Array[V], String, V)], Array[Variable]) -> Self[V]

等式行为 (coeffs, b)；不等式行为 (coeffs, relation, b)，其中 relation 必须恰好是 ">=" 或 "<="。其他关系会导致 abort。两个数组默认都为空。

`Constraint::add_eqpoly`, `Constraint::add_ineqpoly`

add_eqpoly 追加一个等式，add_ineqpoly 追加一个不等式。

pub fn[V] Constraint::add_eqpoly(Self[V], Poly[V], V) -> Unit
pub fn[V] Constraint::add_ineqpoly(Self[V], Poly[V], String, V) -> Unit

add_ineqpoly 不检查关系字符串；请传入 ">=" 或 "<="。

`Constraint::change_eqpoly`, `Constraint::change_ineqpoly`

它们替换从 1 开始计数的某个位置上的约束。

pub fn[V] Constraint::change_eqpoly(Self[V], Int, Poly[V], V) -> Unit
pub fn[V] Constraint::change_ineqpoly(Self[V], Int, Poly[V], String, V) -> Unit

change_eqpoly(i, poly, b) 替换第 i 个等式（从 1 开始计数）。当 i 不为正，或不小于等式数组的容量（数组的实现细节）时，两个函数都会 abort。

`Constraint::to_matrix`

to_matrix 为每个等式返回一个稠密行 [coefficients..., b]。

pub fn[V : @luna-generic.Zero] Constraint::to_matrix(Self[V], Array[Variable]) -> Array[Array[V]]

列的顺序与给定变量的顺序一致。不等式会被忽略，因此请对标准形规划的约束调用它。

test "constraints" {
  let vars = nonneg(["x1", "x2"])
  let c : Constraint[Double] = Constraint::from_array(
    eq_array=[([1.0, 1.0], 4.0)],
    ineq_array=[([2.0, 1.0], "<=", 6.0)],
    vars,
  )
  inspect(
    c,
    content=(
      #|s.t. x1 + x2 = 4
      #|     2x1 + x2 <= 6
      #|    
    ),
  )
  debug_inspect(c.to_matrix(vars), content="[[1, 1, 4]]")
}

线性规划

`Lp`

Lp[V] 是线性规划：包括变量、目标函数、约束，以及由它们构造的规划矩阵。

type Lp[V]
pub impl[V : Show + Eq + @luna-generic.Semiring] Show for Lp[V]
pub fn[V : Show + Eq + @luna-generic.Semiring] Lp::to_string(Self[V]) -> String

构造规划的函数都返回新的 Lp，不修改参数。Show 打印目标函数、约束以及每个变量及其边界。

`Lp::new`, `Lp::cons_from_Array`, `Lp::reset_obj_byarray`

Lp::new 创建一个没有约束的规划；另外两个函数分别替换其约束或目标系数。

pub fn[V : @luna-generic.Zero] Lp::new(Obj_func[V], Array[Variable]) -> Self[V]
pub fn[V : Eq + @luna-generic.Semiring] Lp::cons_from_Array(Self[V], eq_array? : Array[(Array[V], V)], ineq_array? : Array[(Array[V], String, V)]) -> Self[V]
pub fn[V : Compare + @luna-generic.Semiring + @luna-generic.Zero] Lp::reset_obj_byarray(Self[V], Array[V]) -> Self[V]

cons_from_Array 用规划的变量按 Constraint::from_array 的方式替换全部约束。reset_obj_byarray(coeffs) 保留优化方向并设置新的目标系数。

`Lp::from_matrix`

Lp::from_matrix 由规划矩阵创建只含等式约束的规划。

pub fn[V : Eq + @luna-generic.Semiring] Lp::from_matrix(@mutable.Matrix[V], Array[Variable], String) -> Self[V]

Lp::from_matrix(m, vars, sense) 从第 0 行（去掉最后一项）读取目标函数，并从其后每一行 $[a_i \mid b_i]$ 读取等式 $a_i \cdot x = b_i$ 。规划直接把 m 本身作为其矩阵，不做复制。

`Lp::to_standard`

to_standard 返回等价的标准形规划。

pub fn[V : Eq + @luna-generic.Semiring + Compare + Neg + @luna-generic.Zero] Lp::to_standard(Self[V]) -> Self[V]

结果为最小化问题，只含等式约束，且右端项非负：

最大化目标被取负（Obj_func::to_min）；
$b < 0$ 的等式乘以 $-1$ ；
每个不等式获得一个新变量 y1、y2……（按不等式顺序计数），在必要时先将该行乘以 $-1$ 使 $b \ge 0$ ，然后在 <= 行中以系数 $+1$ 、在 >= 行中以系数 $-1$ 加入。

新变量追加在规划变量之后。设计说明推导了每种情形。

`Lp::get_coeff_matrix`, `Lp::get_objfunc_vector`, `Lp::get_b_vector`

它们从规划矩阵中读取 $A$ 、 $c$ 与 $b$ 。

pub fn[V] Lp::get_coeff_matrix(Self[V]) -> @mutable.Matrix[V]
pub fn[V] Lp::get_objfunc_vector(Self[V]) -> Array[V]
pub fn[V] Lp::get_b_vector(Self[V]) -> Array[V]

规划矩阵只包含目标函数和等式约束，因此对于含不等式的规划，这些访问函数只有在 to_standard 之后才描述该规划。

test "standard form" {
  let vars = nonneg(["x1", "x2"])
  let lp : Lp[Double] = Lp::new(Obj_func::from_array("max", [3.0, 4.0], vars), vars)
    .cons_from_Array(ineq_array=[([2.0, 1.0], "<=", 12.0), ([1.0, 3.0], ">=", 2.0)])
  let s = lp.to_standard()
  inspect(
    s,
    content=(
      #|min  -3x1 + -4x2
      #|s.t. 2x1 + x2 + y1 = 12
      #|     x1 + 3x2 + -1y2 = 2
      #|      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, 2]")
  inspect(
    s.get_coeff_matrix(),
    content=(
      #||2, 1, 1, 0|
      #||1, 3, 0, -1|
    ),
  )
}

求解

`Lp::two_stage`

two_stage 以两阶段单纯形法求解标准形线性规划。

pub fn[V : @luna-generic.Zero + @luna-generic.One + Compare + Mul + Div + Sub + ApproximatelyZero + Show + Neg] Lp::two_stage(Self[V]) -> (Array[V], V)

请对 to_standard 的结果调用它。它按规划变量的顺序返回标准形所有变量（包括新增的 y 变量）的取值，以及标准形的最优值。标准形总是最小化，因此对于最大化规划，最大值是返回值的相反数。

当规划不可行（第一阶段结束时人工目标不为零）或无界时，two_stage 会 abort。每次单纯形运行在 1000 次迭代后停止；若尚未结束，会打印 Maximum iterations reached, may not have converged，随后结果照常返回，不报错。

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

`Lp::pivot`

pivot 原地对单纯形表执行一次 Gauss–Jordan 主元变换。

pub fn[V : @luna-generic.Zero + Compare + Div + Mul + Sub] Lp::pivot(@mutable.Matrix[V], Int, Int) -> Unit

Lp::pivot(t, r, s) 将第 r 行除以主元 $t_{rs}$ ，再从其他每一行 i（包括第 0 行）中减去新第 r 行的 $t_{is}$ 倍，使第 s 列成为单位向量 $e_r$ 。当位置超出矩阵或主元恰为零时会 abort。

`Lp::simplex_iteration`

simplex_iteration 原地在单纯形表上运行单纯形法。

pub fn[V : @luna-generic.Zero + Compare + Sub + Mul + Div + Show] Lp::simplex_iteration(@mutable.Matrix[V], max_iterations? : Int, debug? : Bool) -> Unit

单纯形表必须处于典范形式：第 0 行存放检验数（约简成本），最后一列为 $-z$ ，每个基列都是单位向量，且右端项非负。每次迭代选择第 0 行中最负的项所在列进基（相同时取最左列），并在 $\bar a_{is} > 0$ 的行中选择比值 $\bar b_i / \bar a_{is}$ 最小的行出基（相同时取最上面的行）。当第 0 行没有负项时停止；当进基列没有正项时以 Problem is unbounded abort；达到 max_iterations 次迭代（默认 1000）时打印消息并停止。设置 debug=true 时，会打印每次转轴及其后的单纯形表。

test "simplex by hand" {
  // max 3x1 + 4x2, 2x1 + x2 + y1 = 12, x1 + 3x2 + y2 = 10; 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],
  ])
  Lp::pivot(t, 2, 1) // x2 enters, y2 leaves
  inspect(t[0][4], content="13.333333333333334")
  Lp::simplex_iteration(t)
  inspect(t[0][4], content="22")
}

`Lp::phase_1`, `Lp::phase_2`

phase_1 与 phase_2 是 two_stage 的两个阶段。

pub fn[V : @luna-generic.Zero + @luna-generic.One + Compare + Sub + Compare + Div + Mul + Show] Lp::phase_1(@mutable.Matrix[V], Array[Int]) -> @mutable.Matrix[V]
pub fn[V : @luna-generic.Zero + @luna-generic.One + ApproximatelyZero + Div + Sub + Mul + Show + Compare] Lp::phase_2(Array[Variable], Array[V], @mutable.Matrix[V]) -> (Array[V], V)

Lp::phase_1(t, artificial_index) 接受一个扩展了人工列的单纯形表，其第 0 行在每个人工列中为 $1$ ，其余为 $0$ 。artificial_index[i] 是约束 i 的人工变量所在列，若该约束没有人工变量则为 0。它从第 0 行中减去每个带人工变量的约束行，使单纯形表成为典范形式，然后运行 simplex_iteration 并返回同一个矩阵。

Lp::phase_2(vars, c, t) 检查第一阶段是否达到零，把前 vars.length() 个变量的列和右端项复制到新的单纯形表中，把目标 c 写入第 0 行，从第 0 行中消去基列（在 ApproximatelyZero 意义下等于单位向量的列），运行 simplex_iteration，并返回基本解以及第 0 行的最后一项，即 $-z$ 。若第一阶段以非零目标结束，它以 W* from Phase1 isn't zero, Lp doesn't have solution abort。

由规划构造第一阶段单纯形表和 artificial_index 的辅助函数是私有的；除非自己构造单纯形表，否则请使用 two_stage。

数值容差

`ApproximatelyZero`

ApproximatelyZero 判断某个值在第二阶段中是否算作零。

pub trait ApproximatelyZero {
  fn is_zero_eps(Self) -> Bool
}
pub impl ApproximatelyZero for Int
pub impl ApproximatelyZero for Double

Int 只有等于 0 时才为零。Double 在 $|x| < 10^{-15}$ 时为零，这是一个绝对阈值。为自己的系数类型实现该 trait，即可将其用于 two_stage。

`double_equal_to_zero`

double_equal_to_zero 调用 is_zero_eps。

pub fn[V : ApproximatelyZero] double_equal_to_zero(V) -> Bool

test "tolerance" {
  assert_true(@lp.double_equal_to_zero(1.5e-17))
  assert_false(@lp.double_equal_to_zero(1.0e-12))
  assert_true(@lp.double_equal_to_zero(0))
}

已弃用

以下提升而来的方法为保持源码兼容而保留。它们在接口文件中隐藏，从其他包调用时会发出警告。

方法	替代
`Variable::not_equal`, `Poly::not_equal`	`a != b`
`Variable::op_lt`, `op_le`, `op_gt`, `op_ge`	`<`, `<=`, `>`, `>=`
`Variable::output`, `Poly::output`, `Obj_func::output`, `Constraint::output`, `Lp::output`	`to_string()` 或 `"\{x}"`

core API

变量

Variable

Variable::new

Variable::show_all

Variable::to_string, Variable::equal, Variable::compare

Variable::one

线性表达式

Poly

Poly::new, Poly::from_var_array

Poly::add_term_inplace

Poly::copy

Poly::to_array, Poly::to_vector

Poly::neg, Poly::equal, Poly::to_string

目标函数

Obj_func

Obj_func::new, Obj_func::from_array

Obj_func::set_poly, Obj_func::add_term_inplace

Obj_func::judge_max, Obj_func::to_min

Obj_func::to_vector

约束

Constraint

Constraint::new, Constraint::from_array

Constraint::add_eqpoly, Constraint::add_ineqpoly

Constraint::change_eqpoly, Constraint::change_ineqpoly

Constraint::to_matrix