🪬
Rustで数理最適化に触れてみる
はじめに
本記事では、Rust製の混合整数線形計画モデリング用ライブラリ(crate)であるgood_lpを使って、いくつかの簡単な数理最適化問題に触れてみたいと思います。good_lpが対応しているソルバーの一覧はこちらです。有名どころだとcoin_cbcやscipなどが対応しているようです。
本記事ではデフォルトに指定されているcoin_cbcを使うことを前提にしたコードを書きます。
また、本記事で前提としてる環境やバージョンは下記の通りになります。
- macOS Ventura 13.3.1
- Rust: 1.67.1
- good_lp: 1.4.0
準備
ソルバーのダウンロード
brew install cbc
good_lpの追加（参考）
Cargo.toml
[dependencies]
good_lp = "1.4.0"
例1
まずは実数
この線形計画問題の実行可能解は下図の灰色に塗られた領域で表されます。図に示している通り、この例では丁度
ではそれを踏まえて、実際にこの問題good_lpを使って解いてみます。
コードは以下の通りで、実行するとちゃんと目的関数の値として4が得られることが確認出来ます。
main.rs
use std::error::Error;
use good_lp::{constraint, variable, variables, default_solver, Solution, SolverModel};
fn main() -> Result<(), Box<dyn Error>> {
let mut vars = variables!();
// 0以上の値をとる決定変数xとyを定義
let x = vars.add(variable().min(0));
let y = vars.add(variable().min(0));
let solution = vars.maximise(x + y) // 目的関数の設定
.using(default_solver) // ソルバーとしてcoin_cbcが選択される
.with(constraint!(3 * x - y <= 4)) // 制約の追加
.with(constraint!(x + 2 * y <= 6))
.with(constraint!(x + 6 * y <= 16))
.solve()?;
// 解の表示（最適解のコストは4）
println!("x + y = {}", solution.model().obj_value());
println!("x={}, y={}", solution.value(x), solution.value(y));
Ok(())
}
※ちなみに変数はvariablesマクロを使って以下のように定義することも可能です。
/*
以下3行のコードと同じく0以上の連続変数x,yが定義される
let mut vars = variables!();
let x = vars.add(variable().min(0));
let y = vars.add(variable().min(0));
*/
variables! {
vars:
x >= 0;
y >= 0;
}
例2: 割当問題
次に割当問題を解いてみたいと思います。割当問題はタスク
-
c_{i,j} i j - 制約1: 各エージェントは一つのタスクのみを担当する
- 制約2: 各タスクは一つのエージェントのみに割り当てられる
この問題はgood_lpを使って以下のようなコードで表現することが出来ます。
main.rs
use std::{error::Error, collections::HashMap};
use good_lp::{constraint, variable, variables, default_solver, SolverModel, Expression};
const NUM_TASKS: usize = 3;
const NUM_AGENTS: usize = 3;
fn main() -> Result<(), Box<dyn Error>> {
let costs = [ // c_ij: タスクiをエージェントjに割り当てるコスト
[4, 5, 6],
[6, 4, 3],
[9, 7, 4]
];
// x_ij: タスクiをエージェントjに割り当てる0-1変数
let mut variables = variables!();
let mut x = HashMap::new();
for i in 0..NUM_TASKS {
for j in 0..NUM_AGENTS {
x.insert((i, j), variables.add(variable().binary()));
}
}
// 目的関数の定義
let objective = x
.iter()
.map(|(&(i, j), &var)| costs[i][j] * var)
.sum::<Expression>();
// 問題の定義
let mut problem = variables.minimise(objective).using(default_solver);
// 制約追加:各エージェントは一つのタスクのみを担当する
for j in 0..NUM_AGENTS {
let agent_constraint = (0..NUM_TASKS)
.map(|i| x[&(i, j)])
.sum::<Expression>();
problem = problem.with(constraint!(agent_constraint == 1));
}
// 制約追加:各タスクは一つのエージェントのみに割り当てられる
for i in 0..NUM_TASKS {
let job_constraint = (0..NUM_AGENTS)
.map(|j| x[&(i, j)])
.sum::<Expression>();
problem = problem.with(constraint!(job_constraint == 1));
}
let solution = problem.solve()?;
// 解の表示（最適解のコストは12）
println!("objective value: {}", solution.model().obj_value());
Ok(())
}
例3: 輸送問題
最後に割当問題の拡張である輸送問題を解いてみます。輸送問題は工場
-
c_{i,j} i j - 制約1: 各顧客の需要を満たす
- 制約2: 各工場の生産能力を超えない
※全ての
main.rs
use std::{error::Error, collections::HashMap};
use good_lp::{constraint, variable, variables, default_solver, SolverModel, Expression};
const NUM_FACTORIES: usize = 3;
const NUM_CUSTOMERS: usize = 5;
fn main() -> Result<(), Box<dyn Error>> {
let costs = [ // c_ij: 工場iから顧客jへの輸送コスト
[3, 4, 5, 6, 7],
[4, 5, 6, 9, 8],
[6, 4, 3, 4, 6],
[9, 7, 4, 6, 4]
];
let capacities = [520, 450, 500]; // C_i: 工場iの生産能力
let demands = [100, 300, 210, 160, 200]; // d_j: 顧客jの需要
// x_ij: 工場iから顧客jへの輸送量を表す連続変数
let mut variables = variables!();
let mut x = HashMap::new();
for i in 0..NUM_FACTORIES {
for j in 0..NUM_CUSTOMERS {
x.insert((i, j), variables.add(variable().min(0)));
}
}
// 目的関数の定義
let objective = x
.iter()
.map(|(&(i, j), &var)| costs[i][j] * var)
.sum::<Expression>();
// 問題の定義
let mut problem = variables.minimise(objective).using(default_solver);
// 制約追加:各顧客の需要を満たす
for j in 0..NUM_CUSTOMERS {
let demand_constraint = (0..NUM_FACTORIES)
.map(|i| x[&(i, j)])
.sum::<Expression>();
problem = problem.with(constraint!(demand_constraint == demands[j]));
}
// 制約追加:各工場の生産能力を超えない
for i in 0..NUM_FACTORIES {
let capacity_constraint = (0..NUM_CUSTOMERS)
.map(|j| x[&(i, j)])
.sum::<Expression>();
problem = problem.with(constraint!(capacity_constraint <= capacities[i]));
}
let solution = problem.solve()?;
// 解の表示(最適化のコストは3890)
println!("objective value: {}", solution.model().obj_value());
Ok(())
}
参考文献
本記事を書くにあたって、以下の記事や書籍を参考にさせて頂きました。
Discussion