core チュートリアル

このチュートリアルでは、線形計画を記述し、標準形にし、linear-program の 2 段階シンプレックス法で解く方法を示します。2 変数の生産計画問題から始め、最後にシンプレックス表の上でシンプレックス法を手で実行します。各ステップの背後にある数学は設計ノートにあります。

クイックスタート

モジュール名は Luna-Flow/ 名前空間のない linear-program で、まだ 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$ です。解が 4 成分あるのは to_standard が不等式ごとにスラック変数を 1 つ加えるためで、値の符号が反転しているのは標準形が $-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$ は $y_1 \ge 0$ のもとで $2x_1 + x_2 - y_1 = 12$ になりました。

最小化問題を解く

最小化問題では符号を変える必要はありません。返り値が最小値です。

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")
}

ここではどのスラック列も単位列ではないので、第 1 段階は 2 つの人工変数を加えてゼロまで下げ、続いて第 2 段階が第 1 段階の見つけた実行可能点から本来の目的関数を最小化します。

等式制約の計画を行列で書く

すべての制約が等式なら、計画全体が 1 つの行列に収まります。第 0 行に目的係数を、以降の各行に制約 $[a_i \mid b_i]$ を 1 つずつ入れます。次の計画は $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$ です。第 0 行の最後の要素は from_matrix によって無視されます。

式を項ごとに組み立てる

Poly と Obj_func は 1 項ずつ組み立てることもでき、係数がデータから来る場合に便利です。

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$ の被約費用がまだ負であることは、もう 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 でピボットすると誤った結果になります。
許容誤差は絶対的です。 第 2 段階は $|x| < 10^{-15}$ をゼロとみなします。係数が大きい計画では丸め誤差がこのしきい値を超え、実行可能な計画が実行不可能と報告されることがあります。係数が適度な大きさになるよう行をスケールしてください。
変数に y1、y2、… と名付けないでください。 to_standard はこれらの名前をスラック変数に使い、変数は名前で識別されます。

次のステップ

API リファレンスは、Lp::phase_1 と Lp::phase_2 を含むすべての型と関数を説明しています。
設計ノートは、標準形、最適性と非有界性の判定、および 2 つの段階を導出しています。
linear-algebra はシンプレックス表に使う Matrix 型を、luna-generic は係数型を制約する trait を提供します。