Rとomprパッケージで数理最適化を実装する

# 1. 概要
この記事では、R言語の数理最適化パッケージ「ompr」を使って、混合整数計画問題の実装方法を紹介します。
omprは比較的新しいパッケージであり、日本語情報がまだ少ないため、この記事が皆さんの役に立てば幸いです。

pythonであればpulpが有名ですが、pulpと比較してomprの使い勝手や学習コストは同程度かと感じています。
むしろRならではのパイプ記法を使える分、Rユーザーには馴染みやすいかもしれません。
なお、本記事では最適化の数理面での解説は行わず、Rでの実装方法のみ紹介します。

# 2. 必要なパッケージ
```R
install.packages("dplyr")
install.packages("ompr")
install.packages("ompr.roi")
install.packages("ROI.plugin.glpk")
```
ompr.roiは、ROI(R Optimization Infrastructure)という最適化問題を解くためのフレームワーク上で様々なソルバーを利用できるようにする補助的なパッケージです。
今回はGLPKというソルバーを利用するので、ROI.plugin.glpkをインストールします。

# 3. omprを使って混合整数計画問題を解く方法
## 3.1. 問題例
下記の簡単なモデルを考えます。
xとyの2つの変数について、制約条件を満たすように目的関数を最大化するものです。
(なお、本章の問題例は https://dirkschumacher.github.io/ompr/ を引用・一部改変したものです)

- 変数
x >= 0 （整数変数）
1 <= y <= 5 （連続変数）
- 目的関数
maximize : x + y
- 制約条件
x + y <= 11.25

## 3.2. モデルの構築
```R
library(dplyr)
library(ompr)

model <- MIPModel() |>  # モデル作成
  add_variable(x, type = "integer", lb = 0) |>  # 変数xの作成
  add_variable(y, type = "continuous", lb = 1, ub = 5) |>  # 変数yの作成
  set_objective(x + y, "max") |>  # 目的変数の設定
  add_constraint(x + y <= 11.25)  # 制約条件の設定
```
- MIPModel : 混合整数計画問題のモデル作成
- add_variable : 変数を追加 / 整数or連続値orバイナリの指定 / 最小値と最大値をlb(Low Bound)とub(Upper Bound)で指定可
- set_objective : 目的関数の追加
- add_constraint : 制約条件の追加

## 3.3. ソルバーの選択と実行
```R
library(ompr.roi)
library(ROI.plugin.glpk)

result <- model |>  
  solve_model(with_ROI(solver = "glpk"))  # ソルバーの指定と実行
```

## 3.4. 結果の取得
resultから結果を取得できます。
get_solution()を使えば個別の結果も取得できます。
```R
result$status  # 最適解が得られているか確認
# "success"
result$objective_value  # 目的関数の結果
# 11.25
result$solution  # 変数の結果
#    x     y 
# 10.00  1.25

get_solution(result, x) # 変数xの解を取得
# x: 10.00
get_solution(result, y) # 変数yの解を取得
# y: 1.25
```

このくらい簡単な数式であれば直観的に記述できます。
見ての通りパイプを利用して、「モデル宣言→変数追加→目的関数追加→制約条件追加...」のように、Rユーザーなら馴染みのある形でモデル作成が出来ます。


# 4. データフレームを使った実装
一般的にデータはデータフレーム型で持つことが多いので、データフレームを使った実装方法を紹介します。

実装例としてナップザック問題を使います。
ナップザック問題とは、限られた容量のナップザックに、価値や重量・サイズの異なるアイテムを組み合わせて入れることを想定します。
この際、重量やサイズ容量を超えない範囲で、選んだアイテムの総価値を最大化する問題です。

## 4.1. 問題例
下記のモデルを考えます。
イメージとして、ナップザックの重量制限とサイズ制限を超えないように幾つかのお菓子を入れていき、満足度が最大となるようにします。
（重量制限とサイズ制限について、例えばお菓子のうち、綿飴は軽いけどサイズが大きいのでナップザックに入らない...みたいなケースを想定して制約条件としています）

- 変数
x[i] ∈ {0, 1, 2} (i = 1, 2, ..., n)
- 目的関数
最大化 Σ(items$value[i] * x[i]) (i = 1, 2, ..., n)
- 制約条件
Σ(items$weight[i] * x[i]) <= m (i = 1, 2, ..., n)
Σ(items$size[i] * x[i]) <= s (i = 1, 2, ..., n)

ここで、nはお菓子の種類、mはナップザックの重量制限、sはナップザックのサイズ制限を表します。

## 4.2. データ準備
```R:データの準備
items <- data.frame(
  okashi = c("お菓子A", "お菓子B", "お菓子C", "お菓子D", "お菓子E"),
  value = c(2, 3, 1, 3, 4),  # お菓子の価値
  weight = c(3, 4, 2, 5, 6),  # お菓子の重量
  size = c(7, 8, 2, 3, 10)  # お菓子のサイズ
)
n <- nrow(items) # お菓子の数:5個
m <- 20 # ナップザックの重量制限
s <- 30 # ナップザックのサイズ制限
```
## 4.3. モデルの構築
上記のデータテーブルを元に、モデルを作成します。
```R
model <- MIPModel() |>  
  add_variable(x[i], i = 1:n, type = "integer", lb = 0, ub = 2) |>  # n個分だけ変数作成
  set_objective(sum_expr(items$value[i] * x[i], i = 1:n), "max") |>  # 目的関数の設定
  add_constraint(sum_expr(items$weight[i] * x[i], i = 1:n) <= m) |>   # 重量の制約条件
  add_constraint(sum_expr(items$size[i] * x[i], i = 1:n) <= s)  # サイズの制約条件
```
- add_variable : "i = 1:n"と指定することで、nの数だけ変数を一括で作成可能。今回の例では変数x[1]~x[5]を作成し、それぞれお菓子A〜お菓子Eの数に対応します。ちなみに、add_variable(x[i, j], i = 1:n, j = 1:n)のようにiとjを用いて変数の行列を表すことができます。
- sum_expr : 線形和を計算。例えば目的関数で、「お菓子[i]の個数 × お菓子[i]の重量」の合計を算出しています。

## 4.4. 結果の取得
```R
result <- solve_model(model, with_ROI(solver = "glpk"))
result$status #最適解が得られたか確認
# "success"

result$solution
# x[1] x[2] x[3] x[4] x[5] 
#  0    2    0    1    1 
```
結果として、お菓子Bが2個・お菓子Dが1個・お菓子Eが1個の時に満足度が最大となります。
このようにデータテーブルを用いてもモデル作成ができます。


# 5. まとめ
この記事では、R言語の新しいパッケージ"ompr"を使って混合整数計画問題を解く方法を紹介しました。具体的な問題定式化やモデル構築の手順、ソルバーの選択と実行、そして結果の解釈まで、一通りの手順を説明しました。

"ompr"はシンプルで直感的な記述可能で、R言語を使って混合整数計画問題を解く際に便利なツールです。この記事が皆様の参考になれば幸いです。

# 参考文献
https://dirkschumacher.github.io/ompr/
https://www.r-orms.org/mixed-integer-linear-programming/packages/modelling-milp/