🦀
Rustで線形計画問題を解いてみる
はじめに
最近数理最適化問題をいくつかPythonとJuliaで解いたのでRustでも解いてみようと思いました。今回は代表的な問題を2問解いていますができれば他の問題も解いていきたいです。
使用したクレート
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