🦀

Rustで線形計画問題を解いてみる

7 min read

はじめに

最近数理最適化問題をいくつかPythonとJuliaで解いたのでRustでも解いてみようと思いました。今回は代表的な問題を2問解いていますができれば他の問題も解いていきたいです。

この記事のコード(GitHub)

使用したクレート

good_lpとは

  • good_lpクレートとは簡単で使いやすい線形計画問題用のモデラ。

  • 特徴と制限(リポジトリより引用)

Features and limitations
- Linear programming. This crate currently supports only the definition of linear programs. You cannot use it with quadratic functions. For instance: you can maximise 3 * x + y, but not 3 * x * y.
- Continuous and integer variables. good_lp itself supports mixed integer-linear programming (MILP), but not all underlying solvers support integer variables.
- Not a solver. This crate uses other rust crates to provide the solvers. There is no solving algorithm in good_lp itself. If you have an issue with a solver, report it to the solver directly. See below for the list of supported solvers.

0-1ナップサック問題

参考: Numerical Optimizer SIMPLE例題集 V19 - 2.5 ナップサック問題

規定の容量以内で荷物の総価値を最大化する問題。今回は0-1ナップサック問題なので各荷物は1つまでしか入れられない。

例題(上記参考サイトから)

ナップサックの容量: 65

品物 1個あたりの価値 1個あたりのサイズ
缶コーヒー 120 10
水入りペットボトル 130 12
バナナ 80 7
りんご 100 9
おにぎり 250 21
パン 185 16

方針

  • 変数: 各荷物を入れる個数(binary)
  • 目的関数: ナップサックに入れた荷物の価値の合計
  • 制約
    • 各品物は0個以上
    • ナップサックに入っている荷物の総サイズが65以下

ソースコード

src/01knapsack.rs
use std::error::Error;

use good_lp::{constraint, default_solver, variables, Solution, SolverModel};

fn main() -> Result<(), Box<dyn Error>> {
    let yen: Vec<i32> = vec![100, 130, 70, 40, 90, 30, 20, 20];
    let weight: Vec<i32> = vec![80, 70, 60, 50, 40, 30, 20, 10];
    // 変数の定義
    variables!(vars:x[yen.len()] (binary));
    let solution = vars
        // 目的関数(最大化)の定義
        .maximise(
            x[0] * yen[0]
                + x[1] * yen[1]
                + x[2] * yen[2]
                + x[3] * yen[3]
                + x[4] * yen[4]
                + x[5] * yen[5]
                + x[6] * yen[6]
                + x[7] * yen[7],
        )
        // defaul_solverはCargo.tpmlで選択済み(今回はcoin_cbc)
        .using(default_solver)
        // 制約の定義
        .with(constraint!(
            x[0] * weight[0]
                + x[1] * weight[1]
                + x[2] * weight[2]
                + x[3] * weight[3]
                + x[4] * weight[4]
                + x[5] * weight[5]
                + x[6] * weight[6]
                + x[7] * weight[7]
                <= 100
        ))
        .solve()?;
    // 結果の表示
    for i in 0..yen.len() {
        print!("x{} : {}\n", i, solution.value(x[i]));
    }
    println!(
        "Total yen: {}",
        solution.eval(
            solution.value(x[0]) as i32 * yen[0]
                + solution.value(x[1]) as i32 * yen[1]
                + solution.value(x[2]) as i32 * yen[2]
                + solution.value(x[3]) as i32 * yen[3]
                + solution.value(x[4]) as i32 * yen[4]
                + solution.value(x[5]) as i32 * yen[5]
                + solution.value(x[6]) as i32 * yen[6]
                + solution.value(x[7]) as i32 * yen[7]
        )
    );
    Ok(())
}

実行結果

Math-Optim-Rust via 🦀 v1.55.0
❯ cargo run --bin 01knapsack

~ 省略 ~

Objective value:                170.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.02

Total time (CPU seconds):       0.01   (Wallclock seconds):       0.03

x0 : 0
x1 : 1
x2 : 0
x3 : 0
x4 : 0
x5 : 0
x6 : 1
x7 : 1
Total yen: 170

多期間計画問題

参考: Numerical Optimizer SIMPLE例題集 V19 - 2.3 多期間計画問題

他機関にわたって意思を決定する問題。各期間ごとに変数を定義するのが一般的。

例題(上記参考サイトから)

  • 原料使用料
I II
A 2 7
B 5 3
  • 生産/在庫コスト
I II
生産 75 50
在庫 8 7
  • 月毎の各製品の出荷量
I II
1ヶ月目 30 20
2ヶ月目 60 50
3ヶ月目 80 90
  • 月毎の原材料の利用可能量
A B
1ヶ月目 920 790
2ヶ月目 750 600
3ヶ月目 500 480

方針

  • 変数: 各製品の各期間における生産量及び在庫量
  • 目的関数: すべての期間を通しての総コスト(最小)
  • 制約
    • 変数は非負の数
    • 各期間においてそれぞれの原料の利用可能料を超過しない
    • 各期間においてそれぞれの製品の出荷量を満たす

ソースコード

src/multi_period_planning.rs
use std::error::Error;

use good_lp::{constraint, default_solver, variables, Solution, SolverModel};

fn main() -> Result<(), Box<dyn Error>> {
    variables! {
        vars:
            0 <= xi1;
            0 <= xii1;
            0 <= xi2;
            0 <= xii2;
            0 <= xi3;
            0 <= xii3;
            0 <= yi1;
            0 <= yii1;
            0 <= yi2;
            0 <= yii2;
    }
    let solution = vars
        .minimise(
            75 * xi1
                + 50 * xii1
                + 8 * yi1
                + 7 * yii1
                + 75 * xi2
                + 50 * xii2
                + 8 * yi2
                + 7 * yii2
                + 75 * xi3
                + 50 * xii3,
        )
        .using(default_solver)
        .with(constraint!(2 * xi1 + 7 * xii1 <= 920))
        .with(constraint!(5 * xi1 + 3 * xii1 <= 790))
        .with(constraint!(2 * xi2 + 7 * xii2 <= 750))
        .with(constraint!(5 * xi2 + 3 * xii2 <= 600))
        .with(constraint!(2 * xi3 + 7 * xii3 <= 500))
        .with(constraint!(5 * xi3 + 3 * xii3 <= 480))
        .with(constraint!(xi1 - yi1 == 30))
        .with(constraint!(xii1 - yii1 == 20))
        .with(constraint!(xi2 + yi1 - yi2 == 60))
        .with(constraint!(xii2 + yii1 - yii2 == 50))
        .with(constraint!(xi3 + yi2 == 80))
        .with(constraint!(xii3 + yii2 == 90))
        .solve()?;

    println!("製品Iの1ヶ月目の生産量: {}", solution.value(xi1));
    println!("製品IIの1ヶ月目の生産量: {}", solution.value(xii1));
    println!("製品Iの2ヶ月目の生産量: {}", solution.value(xi2));
    println!("製品IIの2ヶ月目の生産量: {}", solution.value(xii2));
    println!("製品Iの3ヶ月目の生産量: {}", solution.value(xii3));
    println!("製品IIの3ヶ月目の生産量: {}", solution.value(xii3));
    println!("製品Iの1ヶ月目の在庫量: {}", solution.value(yi1));
    println!("製品IIの1ヶ月目の在庫量: {}", solution.value(yii1));
    println!("製品Iの2ヶ月目の在庫量: {}", solution.value(yi2));
    println!("製品IIの2ヶ月目の在庫量: {}", solution.value(yii2));
    Ok(())
}

実行結果

Math-Optim-Rust via 🦀 v1.55.0
❯ cargo run --bin mpp

~ 省略 ~

Optimal - objective value 21199.172
After Postsolve, objective 21199.172, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 21199.17241 - 5 iterations time 0.002, Presolve 0.00
製品Iの1ヶ月目の生産量: 38
製品IIの1ヶ月目の生産量: 20
製品Iの2ヶ月目の生産量: 67.86206896551724
製品IIの2ヶ月目の生産量: 86.89655172413794
製品Iの3ヶ月目の生産量: 53.10344827586207
製品IIの3ヶ月目の生産量: 53.10344827586207
製品Iの1ヶ月目の在庫量: 8
製品IIの1ヶ月目の在庫量: 0
製品Iの2ヶ月目の在庫量: 15.862068965517238
製品IIの2ヶ月目の在庫量: 36.89655172413793

終わりに

今回は自分のRustにおけるスマートな書き方がわからなかったため、内積計算もすべてハードコーディングしてある。将来的には.sum()などを用いた可搬性をもたせたコードを書きたい。もしこの記事の誤り、もっとこういうコードを書いたほうがいいなどのご意見がありましたらissue/PRをよろしくおねがいします。

自分のブログに書いた同記事

Discussion

ログインするとコメントできます