ここでは，数理最適化のモデリングについて学ぶ．
線形最適化
- 線形最適化
- 簡単な例題と双対問題
- 輸送問題，多品種輸送問題
- 栄養問題
AMPLによるモデル化
混合整数最適化問題
- 簡単な例題（パズル）
- 多制約ナップサック問題
モデリングのコツ
- 最大値の最小化
- 絶対値
- 離接制約
- if A then B 条件
機械学習と数理最適化の融合
数理最適化ソルバーとモデラー
以下の数理最適化ソルバー（モデラー，モデリング言語）を用いる．
数理最適化ソルバー Gurobi（商用，アカデミックフリー） https://www.gurobi.com/
- 独自のPythonインターフェイス（「あたらしい数理最適化」（近代科学社）で採用）
- 凸2次（制約）整数，2次錐最適化，非凸2次に対応している．
数理最適化モデラー PuLP （MITライセンス） https://coin-or.github.io/pulp/
- メインソルバーはCBC(COIN Branch & Cut; EPLライセンス）であり，混合整数最適化問題に対応している．
- ここでは，Gurobiと同じインターフェイスをもつラッパーであるmypulpパッケージを使う．
- CBCより高速なソルバーにSCIP https://www.scipopt.org/ がある．
最適化モデリング言語AMPL
- 無料版 (Community Edition) が https://ampl.com/ce/ でユーザー登録することによって使える．
- SCIPを含む無料のソルバーやGurobiを商用のソルバーのデモ版が同時にインストールされる．
- 非線形最適化ソルバーも呼び出すことができる．
練習問題を解くためには， 商用のGurobiのライセンスをとるか， mypulpを以下の方法でインストールしておく必要がある．
pip install mypulp線形最適化問題
線形最適化問題の一般形
線形最適化問題の一般形は以下のように書ける．
$$ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{R}^n \end{array} $$ここで，$x$ は変数を表す $n$次元実数ベクトル， $A$ は制約の左辺係数を表す $m \times n$ 行列，$c$ は目的関数の係数を表す $n$次元ベクトル， $b$ は制約式の右辺定数を表す $m$次元ベクトルである．
実際問題を解く際の注意
整数や非線形関数を含まない線形最適化問題は、比較的簡単に解くことができるが、大規模な実際問題を解決する際には、以下のような注意が必要である。
数値の桁数を不用意に大きくしない。例えば、数百億円規模の最適化問題において、小数点以下の桁数が非常に多い問題例においては、線形最適化でも非常に時間がかかることがある。必要最小限の桁数に丸め、目的関数値は10万（６桁）程度に設定すると、計算が安定する。
問題の疎性を考慮する。例えば、複数製品を工場から倉庫、倉庫から顧客まで運ぶサプライ・チェインの最適化を考える。固定費用を考えない場合には、線形最適化問題に帰着されるが、問題例の規模によっては解けないこともある。そういう場合には、工場で生産できない製品や、顧客需要のない製品に対しては経路からあらかじめ除外してネットワークを生成すると劇的に速度が改善されることがある。また、輸送費用が変わらない顧客を集約したり、製品のABC分析を行い、需要の小さい製品はとりあえず除外するか、集約して扱うことも推奨される。
実行不可能にならないように注意する。実際問題では、ユーザーは無理な注文をしがちだ。それをそのまま真に受けて定式化すると実行不可能になることがある。数理最適化ソルバーは、制約を満たす解が存在しないしないことを示してくれるが、どの制約のせいで解がないのかまでは示してくれない。Gurobiには、Irreducible Inconsistent Subsystem (IIS)を計算してくれるメソッドが準備されているが、定式化を行う際に、主要な制約の逸脱を許してソフト制約にしておくことを推奨する。また、全ての費用を合計した1つの目的関数を設定するのではなく、輸送費用、配送費用、倉庫費用など小分けにして計算して変数に保存しておき、最後にそれを合計したものを最適化すると、あとで個別の費用を計算する手間が省ける。
簡単な例題
簡単な例題で線形最適化問題を説明する．
トンコケ，コケトン，ミックスの丼を販売している． 各丼の販売数を変数とし， トンコケ丼 $x_1$，コケトン丼 $x_2$ ，ミックス丼 $x_3$ とする． 3種類の肉（順に豚肉，鶏肉，牛肉）の資源制約の下で，利益を最大化する問題は，以下の線形最適化問題になる．
$$ \begin{array}{l c c c c c} maximize & 15 x_1 & + 18 x_2 & +30 x_3 & & \\ s.t. & 2x_1 & + x_2 & + x_3 & \leq & 60 \\ & x_1 & + 2 x_2 & + x_3 &\leq & 60 \\ & & & x_3 &\leq & 30 \\ & & & x_1,x_2,x_3 & \geq & 0 \end{array} $$これをPythonを用いて解く．
from gurobipy import GRB, Model, quicksum, multidict, tuplelist
#from mypulp import GRB, Model, quicksum, multidict, tuplelist #mypulpを使う場合はここを生かす
model = Model() #モデルインスタンスの生成
x1 = model.addVar(name="x1") #トンコケ丼の変数（既定値は非負の実数）
x2 = model.addVar(name="x2") #コケトン丼の変数
x3 = model.addVar(ub=30.0, name="x3") #ミックス丼の変数（牛肉の上限をub引数で指定）
model.update() # Gurobiの怠惰な更新(lazy update)という仕様（忘れずに！）
model.addConstr(2 * x1 + x2 + x3 <= 60) #豚肉の制約を追加
model.addConstr(x1 + 2 * x2 + x3 <= 60) #鶏肉の制約
model.setObjective(15 * x1 + 18 * x2 + 30 * x3, GRB.MAXIMIZE) #目的関数（利益の最大化）
# SCIPを使用したい場合
# import pulp
# solver = pulp.SCIP()
# model.optimize(solver)
model.optimize() #最適化
if model.Status == GRB.Status.OPTIMAL: #Status属性で最適解が得られたことを確認
print("Opt. Value=", model.ObjVal) #最適値はObjVal属性
for v in model.getVars(): #getVars()メソッドで変数を得る
print(v.VarName, v.X) #変数名は VarName 属性， 最適解は X 属性
モデルファイルの出力
model.write("ファイル名.lp")で LPフォーマット(Linear Programming (LP) format)で保存
model.write("ファイル名.mps")でMPSフォーマット (Mathematical Programming System (MPS) format) で保存
- 可読性はないが，ほとんどの最適化ソルバーが対応している古典的な書式
PuLP (mypulp) だと print(model) でも画面にLPフォーマットを出力する．また，制約オブジェクトをprintすると，展開した式を確認することができる．
Gurobiでの制約オブジェクトの展開式（線形表現）は， model.getRow(制約オブジェクト) で得ることができる
model.write("lo1.lp")
model.write("lo1.mps")
print("LP file =========================")
!cat lo1.lp
print("MPS file =========================")
!cat lo1.mps
print(model)
#print(con) #mypulpの場合
print(model.getRow(con)) #Gurobiの場合
例題 双対問題
上の線形最適化問題において，資源（豚肉，鶏肉，牛肉）の価値を推定したい．
これは，先程の定式化を主問題としたときの双対問題を解けば良い．
主問題 $$ \begin{array}{l c c c c c} maximize & 15 x_1 & + 18 x_2 & +30 x_3 & & \\ s.t. & 2x_1 & + x_2 & + x_3 & \leq & 60 \\ & x_1 & + 2 x_2 & + x_3 &\leq & 60 \\ & & & x_3 &\leq & 30 \\ & & &x_1,x_2,x_3 & \geq & 0 \end{array} $$
双対問題
$$ \begin{array}{l c c c c c} minimize & 60 \pi_1 & + 60 \pi_2& +30 \pi_3 & & \\ s.t. & 2\pi_1 & + \pi_2 & & \geq & 15 \\ & \pi_1 & + 2\pi_2 & &\geq & 18 \\ & \pi_1 & +\pi_2 & +\pi_3 &\geq & 30 \\ & & & \pi_1,\pi_2,\pi_3 & \geq & 0 \end{array} $$実は，主問題を解けば双対問題の最適解も同時に得ることができる． 制約オブジェクトのPi属性に最適双対変数が入っている．
model = Model()
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
x3 = model.addVar(name="x3")
model.update()
c1 = model.addConstr(2 * x1 + x2 + x3 <= 60) #制約オブジェクトを保管
c2 = model.addConstr(x1 + 2 * x2 + x3 <= 60)
c3 = model.addConstr(x3 <=30)
model.setObjective(15 * x1 + 18 * x2 + 30 * x3, GRB.MAXIMIZE)
model.optimize()
print("Opt. Value=", model.ObjVal)
for v in model.getVars():
print(v.VarName, v.X)
print("豚肉=", c1.Pi, "鶏肉=", c2.Pi, "牛肉=", c3.Pi)
model = Model()
a = [5, 4, 2]
x = [model.addVar(name=f"x({i})") for i in range(3)]
L = quicksum(a[i] * x[i] for i in range(3))
model.update()
print(L)
name, height, weight = multidict(
{"Taro": [145, 30], "Hanako": [138, 34], "Simon": [150, 45]}
)
print(name)
print(height)
print(weight)
T = tuplelist(
[("Sara", "Apple"), ("Taro", "Pear"), ("Jiro", "Orange"), ("Simon", "Apple")]
)
print(T.select("*", "Apple"))
例題：多品種輸送問題
上で導入したquicksum, multidict, tuplelistの使用例として，以下の多品種輸送問題を考える．
製品（品種）の集合を $K$ とする． 工場 $j$ から顧客 $i$ に製品 $k$ が輸送される量を表す連続変数を $x_{ijk}$ とする． 顧客 $i$ における製品 $k$ の需要量を $d_{ik}$， 工場 $j$ の生産量上限（容量）を $M_j$， 顧客 $i$ と工場 $j$ 間に製品 $k$ の $1$ 単位の需要が移動するときにかかる輸送費用を $c_{ijk}$ とする．
上の記号を用いると，多品種輸送問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} \displaystyle\sum_{k \in K} c_{ijk} x_{ijk} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ijk} =d_{ik} & \forall i \in I; k \in K \\ & \displaystyle\sum_{i \in I} \displaystyle\sum_{k \in K} x_{ijk} \leq M_j & \forall j \in J \\ & x_{ijk} \geq 0 & \forall i \in I; j \in J; k \in K \end{array} $$最初の制約は，各製品ごとに需要が満たされることを表し， 2番目の制約は，工場で生産されるすべての製品の合計量が，工場の容量を超えないことを表す．
現実的には，すべての工場ですべての製品が製造可能ではない． 工場 $1$ では製品 $2,4$ を，工場 $2$ では製品 $1,2,3$ を，工場 $3$ では製品 $2,3,4$ を製造可能とする．
定式化は製品が輸送可能な経路だけを考えて，制約を作る必要がある． これはtuplelistを使うと綺麗に書ける．
d = {
(1, 1): 80,
(1, 2): 85,
(1, 3): 300,
(1, 4): 6,
(2, 1): 270,
(2, 2): 160,
(2, 3): 400,
(2, 4): 7,
(3, 1): 250,
(3, 2): 130,
(3, 3): 350,
(3, 4): 4,
(4, 1): 160,
(4, 2): 60,
(4, 3): 200,
(4, 4): 3,
(5, 1): 180,
(5, 2): 40,
(5, 3): 150,
(5, 4): 5,
}
I = set([i for (i, k) in d]) #顧客の集合
K = set([k for (i, k) in d]) #製品の集合
J, M = multidict({1: 3000, 2: 3000, 3: 3000}) #工場の容量
produce = {
1: [2, 4],
2: [1, 2, 3],
3: [2, 3, 4],
} # 工場1では製品2,4を，工場2では製品1,2,3を，工場3では製品2,3,4を製造可能
weight = {1: 5, 2: 2, 3: 3, 4: 4} # 製品の重量
cost = {
(1, 1): 4,
(1, 2): 6,
(1, 3): 9, # 地点間の輸送費用（これに製品の重量を乗じたものを輸送費用とする）
(2, 1): 5,
(2, 2): 4,
(2, 3): 7,
(3, 1): 6,
(3, 2): 3,
(3, 3): 4,
(4, 1): 8,
(4, 2): 5,
(4, 3): 3,
(5, 1): 10,
(5, 2): 8,
(5, 3): 4,
}
c = {}
for i in I:
for j in J:
for k in produce[j]:
c[i, j, k] = cost[i, j] * weight[k]
model = Model("multi-commodity transportation")
x = {}
for (i, j, k) in c:
x[i, j, k] = model.addVar(vtype="C", name=f"x({i},{j},{k})")
model.update()
arcs = tuplelist([(i, j, k) for (i, j, k) in x])
for i in I:
for k in K:
model.addConstr(
quicksum(x[i, j, k] for (i, j, k) in arcs.select(i, "*", k)) == d[i, k],
f"Demand({i},{k})",
)
for j in J:
model.addConstr(
quicksum(x[i, j, k] for (i, j, k) in arcs.select("*", j, "*")) <= M[j],
f"Capacity({j})",
)
model.setObjective(quicksum(c[i, j, k] * x[i, j, k] for (i, j, k) in x), GRB.MINIMIZE)
model.optimize()
print("Optimal value:", model.ObjVal)
EPS = 1.0e-6
for (i, j, k) in x:
if x[i, j, k].X > EPS:
print(f"{x[i,j,k].X} units of commodity {k} from {j} to {i}")
arcs.select(1, 1, "*")
例題：栄養問題
あなたは，某ハンバーガーショップの調査を命じられた健康オタクの諜報員だ．
あなたは任務のため，毎日ハンバーガーショップだけで食事をしなければならないが， 健康を守るため，なるべく政府の決めた栄養素の推奨値を遵守しようと考えている．
考慮する栄養素は，カロリー(Cal)，炭水化物(Carbo)，タンパク質(Protein)， ビタミンA(VitA)，ビタミンC(VitC)，カルシウム(Calc)，鉄分(Iron)であり， 1日に必要な量の上下限は，以下の表の通りとする．
現在，ハンバーガーショップで販売されている商品は， CQPounder, Big M, FFilet, Chicken, Fries, Milk, VegJuice の7種類だけであり， それぞれの価格と栄養素の含有量は，以下の表のようになっている．
さらに，調査費は限られているので，なるべく安い商品を購入するように命じられている． さて，どの商品を購入して食べれば，健康を維持できるだろうか？
| 栄養素 $N$ | Cal | Carbo | Protein | VitA | VitC | Calc | Iron | 価格 |
|---|---|---|---|---|---|---|---|---|
| 商品名 $F$ | $n_{ij}$ | $c_j$ | ||||||
| CQPounder | 556 | 39 | 30 | 147 | 10 | 221 | 2.4 | 360 |
| Big M | 556 | 46 | 26 | 97 | 9 | 142 | 2.4 | 320 |
| FFilet | 356 | 42 | 14 | 28 | 1 | 76 | 0.7 | 270 |
| Chicken | 431 | 45 | 20 | 9 | 2 | 37 | 0.9 | 290 |
| Fries | 249 | 30 | 3 | 0 | 5 | 7 | 0.6 | 190 |
| Milk | 138 | 10 | 7 | 80 | 2 | 227 | 0 | 170 |
| VegJuice | 69 | 17 | 1 | 750 | 2 | 18 | 0 | 100 |
| 上限 $a_i$ | 3000 | 375 | 60 | 750 | 100 | 900 | 7.5 | |
| 下限 $b_i$ | 2000 | 300 | 50 | 500 | 85 | 660 | 6.0 |
商品の集合を $F$ （Foodの略），栄養素の集合を $N$（Nutrientの略）とする． 栄養素 $i$ の1日の摂取量の下限を $a_i$， 上限を $b_i$ とし， 商品 $j$ の価格を $c_j$，含んでいる栄養素 $i$ の量を $n_{ij}$ とする． 商品 $j$ を購入する個数を非負の実数変数 $x_j$ で表すと，栄養問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j & \\ s.t. & a_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i & i \in N \\ & x_j \geq 0 & j \in F \end{array} $$さて．実際問題を解く際の注意で書いたように， ユーザーは無理な注文をしがちだ。 ここでは，制約の逸脱を表す以下の2つの変数を追加することによって，対処する．
- $d_i$: 不足変数
- $s_i$: 超過変数
制約は，以下のようにソフトな制約に変更する．
$$a_i - d_i \leq \displaystyle\sum_{j \in F} n_{ij} x_j \leq b_i +s_i \ \ i \in N$$
また，目的関数にも制約の逸脱ペナルティを追加する．ここで， $M$ は十分に大きな数とする．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j \in F} c_j x_j + \displaystyle\sum_{i \in N} M (d_i+s_i) & \end{array} $$F, c, n = multidict(
{
"CQPounder": [
360,
{
"Cal": 556,
"Carbo": 39,
"Protein": 30,
"VitA": 147,
"VitC": 10,
"Calc": 221,
"Iron": 2.4,
},
],
"Big M": [
320,
{
"Cal": 556,
"Carbo": 46,
"Protein": 26,
"VitA": 97,
"VitC": 9,
"Calc": 142,
"Iron": 2.4,
},
],
"FFilet": [
270,
{
"Cal": 356,
"Carbo": 42,
"Protein": 14,
"VitA": 28,
"VitC": 1,
"Calc": 76,
"Iron": 0.7,
},
],
"Chicken": [
290,
{
"Cal": 431,
"Carbo": 45,
"Protein": 20,
"VitA": 9,
"VitC": 2,
"Calc": 37,
"Iron": 0.9,
},
],
"Fries": [
190,
{
"Cal": 249,
"Carbo": 30,
"Protein": 3,
"VitA": 0,
"VitC": 5,
"Calc": 7,
"Iron": 0.6,
},
],
"Milk": [
170,
{
"Cal": 138,
"Carbo": 10,
"Protein": 7,
"VitA": 80,
"VitC": 2,
"Calc": 227,
"Iron": 0,
},
],
"VegJuice": [
100,
{
"Cal": 69,
"Carbo": 17,
"Protein": 1,
"VitA": 750,
"VitC": 2,
"Calc": 18,
"Iron": 0,
},
],
}
)
N, a, b = multidict(
{
"Cal": [2000, 3000],
"Carbo": [300, 375],
"Protein": [50, 60],
"VitA": [500, 750],
"VitC": [85, 100],
"Calc": [660, 900],
"Iron": [6.0, 7.5],
}
)
model = Model("modern_diet")
x, s, d = {}, {}, {}
for j in F:
x[j] = model.addVar(vtype="C", name=f"x({j})")
for i in N:
s[i] = model.addVar(vtype="C", name=f"surplus({i})")
d[i] = model.addVar(vtype="C", name=f"deficit({i})")
model.update()
for i in N:
model.addConstr(quicksum(n[j][i] * x[j] for j in F) >= a[i] - d[i], f"NutrLB({i})")
model.addConstr(quicksum(n[j][i] * x[j] for j in F) <= b[i] + s[i], f"NutrUB({i})")
model.setObjective(
quicksum(c[j] * x[j] for j in F) + quicksum(9999 * d[i] + 9999 * s[i] for i in N),
GRB.MINIMIZE,
)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
print("Optimal value:", model.ObjVal)
for j in F:
if x[j].X > 0:
print(j, x[j].X)
for i in N:
if d[i].X > 0:
print(f"deficit of {i} ={d[i].X}")
if s[i].X > 0:
print(f"surplus of {i} ={s[i].X}")
AMPLによるモデル化
最適化モデリング言語AMPLを用いて（小さな）栄養問題を解いてみよう．
AMPLの無料版が https://ampl.com/ce/ でユーザー登録することによって使える．
ログインすると License UUID: の後にIDが表示されるので，それをクリックしてクリップボードに記憶しておく．
Google Colabの場合には，以下のセルを実行して amplpy とAMPLをインストールする （途中でUUIDを入れるように促されるので，コピーする）．
!pip install -q amplpy
MODULES, LICENSE_UUID = ["highs"], None
from amplpy import tools
ampl = tools.ampl_notebook(modules=MODULES, license_uuid=LICENSE_UUID, globals_=globals())
インストール
Mac, Windows, Linuxの場合には，AMPLをインストールして，コマンドインターフェイス（amplもしくはamplide）を使う．
https://ampl.com/ce/ にログインし，OSにあった圧縮ファイルをダウンロードして展開する．
Macにインストールする際には，使用するコマンド（amplkey, amplなど）をFinderで右クリックで開いて，「開発元を検証できません。開いてもよろしいですか?」というメッセージが出たら「開く」を押す．
ライセンスファイルをダウンロードしたあとに，以下のコマンドでライセンスを登録する．
amplkey activate --uuid コピーしたUUIDamplを起動する．
統合開発環境 amplide を使いたい場合には，以下のコマンドを実行する．
xattr -rc amplide%%writefile diet.mod
set NUTR;
set FOOD;
param cost {FOOD} > 0;
param f_min {FOOD} >= 0;
param f_max {j in FOOD} >= f_min[j];
param n_min {NUTR} >= 0;
param n_max {i in NUTR} >= n_min[i];
param amt {NUTR,FOOD} >= 0;
var Buy {j in FOOD} >= f_min[j], <= f_max[j];
minimize Total_Cost: sum {j in FOOD} cost[j] * Buy[j];
subject to Diet {i in NUTR}:
n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];
%%writefile diet.dat
data;
set NUTR := A B1 B2 C ;
set FOOD := BEEF CHK FISH HAM MCH MTL SPG TUR ;
param: cost f_min f_max :=
BEEF 3.19 0 100
CHK 2.59 0 100
FISH 2.29 0 100
HAM 2.89 0 100
MCH 1.89 0 100
MTL 1.99 0 100
SPG 1.99 0 100
TUR 2.49 0 100 ;
param: n_min n_max :=
A 700 10000
C 700 10000
B1 700 10000
B2 700 10000 ;
param amt (tr):
A C B1 B2 :=
BEEF 60 20 10 15
CHK 8 0 20 20
FISH 8 10 15 10
HAM 40 40 35 10
MCH 15 35 15 15
MTL 70 30 15 15
SPG 25 50 25 15
TUR 60 20 15 10 ;
%%ampl_eval
reset;
model diet.mod;
data diet.dat;
option solver cbc;
solve;
display Buy;
問題 線形最適化(1）
あなたは業務用ジュースの販売会社の社長だ． いま，原料の ぶどう原液 $200$ $\ell$ とりんご原液 $100$ $\ell$を使って 2 種類のミックスジュース（商品名はA,B)を作ろうと思っている．
ジュースAを作るにはぶどう $3$ $\ell$とりんご $1$ $\ell$が必要で， ジュースBを作るにはぶどう $2$ $\ell$とりんご $2$ $\ell$が必要である． （なんとジュースは $4$ $\ell$入りの特大サイズなのだ！）
ジュースAは $3$ 千円，ジュースBは$4$ 千円の値段をつけた．
さて，ジュースAとジュースBをそれぞれ何本作れば，利益が最大になるだろうか？
問題 線形最適化(4）
あなたは丼チェーンの店長だ．店の主力製品は， トンコケ丼，コケトン丼，ミックス丼，ビーフ丼の4 種類で， トンコケ丼を1 杯作るには，$200$gの豚肉と $100$gの鶏肉， コケトン丼を1 杯作るには，$100$gの豚肉と $200$gの鶏肉， ミックス丼を1 杯作るには，豚肉，鶏肉，牛肉を $100$gずつ， 最後のビーフ丼は，牛肉だけを$300$g使う． ただし，ビーフ丼は限定商品のため1日 $10$ 杯しか作れない．
原料として使用できる豚，鶏，牛の肉は，最大1 日あたり $9$ キログラム，$9$kg， $6$kgで， 販売価格は，トンコケ丼1 杯 $1500$ 円，コケトン丼1 杯 $1800$ 円，ミックス丼1 杯 $2000$ 円， そしてビーフ丼は $5000$ 円だ．
さて，お店の利益を最大にするためには，あなたは丼を何杯ずつ作るように指示を出せばよいのだろうか？
また、上の問題において， 豚肉，鶏肉，牛肉の $100$gあたりの価値はいくらになるか計算せよ．
問題 輸送問題と双対性
あなたは，スポーツ用品販売チェインのオーナーだ． あなたは，店舗展開をしている5つの顧客 （需要地点）に対して， 3 つの自社工場で生産した1種類の製品を運ぶ必要がある． 工場の生産可能量（容量）と顧客の需要量は表のようになっている．
$$ \begin{array}{c c c | c c c c c } \hline 工場生産可能容量 & & & 顧需要客 & & & & \\ \hline 1 & 2 & 3 & 1 & 2 & 3 & 4 & 5 \\ \hline 500 & 500 & 500 & 80 & 270 & 250 & 160 & 180 \\ \hline \end{array} $$ただし，工場から顧客へ製品を輸送する際は必ず2箇所の倉庫のいずれかを経由しなければならない． 工場と倉庫間，倉庫と顧客間の輸送費用は， 以下の表のようになっているとしたとき， どのような輸送経路を選択すれば，総費用が最小になるであろうか？
$$ \begin{array}{c | c c c | c c c c c } \hline & 工場 & & & 顧客 & & & & \\ \hline 倉庫 & 1 & 2 & 3 & 1 & 2 & 3 & 4 & 5 \\ \hline 1 & 1 & 2 & 3 & 4 & 5 & 6 & 8 & 10 \\ 2 & 3 & 1 & 2 & 6 & 4 & 3 & 5 & 8 \\ \hline \end{array} $$上の輸送問題において，各顧客の追加注文は $1$ 単位あたりいくらの費用の増加をもたらすのか，双対性の概念を用いて計算せよ．
費用削減のためには，どの工場を拡張すれば良いかを，やはり双対性を用いて考えよ．
問題 ゼロ和ゲーム (1)
サッカーのPK（ペナルティキック）の最適戦略を考える．
いま，ゴールキーパーは（キッカーから見て）左に飛ぶか右に飛ぶかの2 つの戦略を持っており， キッカーは左に蹴るか右に蹴るかの2 つの戦略を持っているものとする． 得点が入る確率は，両選手の得意・不得意から以下のような確率になっているものとする（行がキーパーで，列がキッカーの戦略である）．
$$ \begin{array}{ c | c c} \hline & 左 & 右 \\ \hline 左 & 0.9 & 0.5 \\ 右 & 0.6 & 0.8 \\ \hline \end{array} $$キーパーは得点が入る確率を最小化したいし，キッカーは得点が入る確率を最大化したい． さて，両選手はどのような行動をとれば良いだろうか？
これは，ゼロ和ゲームの混合戦略を求める問題となり，確率的な行動をとることが最適戦略となる （つまりどちらに飛ぶかはサイコロを振って決めるのだ）．
キーパーが左に行く確率を変数 $L$，右に行く確率を変数 $R$ として線形最適化問題として定式化を行う． キッカーの戦略を（たとえば左に）固定した場合には，ゲームの値（お互いが最適戦略をとったときにゴールが決まる確率）$V$ は， $$ V \geq 0.9 L + 0.6 R $$ を満たす．キーパーはなるべく $V$ が小さくなるように変数を決め，キッカーは $V$ の値を最大化する．
混合整数最適化問題
（混合）整数最適化問題の一般形
一般の整数線形最適化問題は以下のように書ける．
$$ \begin{array}{ll} minimize & c^T x \\ s.t. & A x \leq b \\ & x \in \mathbf{Z}^n \end{array} $$ここで，$x$ は変数を表す $n$次元ベクトル，$A$ は制約の左辺係数を表す $m \times n$ 行列， $c$ は目的関数の係数を表す $n$次元ベクトル， $b$ は制約式の右辺定数を表す $m$次元ベクトルである．
変数の一部が整数でなく実数であることを許した問題を混合整数最適化問題とよぶ．
例題： 鶴亀算の拡張
簡単な例題で整数線形最適化問題を説明する．
宇宙ステーションに、地球人（頭1つで足２本）、犬星人（頭1つで足４本）、タコ星人（頭１つで足８本）が何人かいる。
頭の数の合計は32、足の数の合計は80である。人間以外の星人の数が最小の場合の、各星人の人数を求めよ。
model = Model("puzzle")
x = model.addVar(vtype="I", name="x")
y = model.addVar(vtype="I", name="y")
z = model.addVar(vtype="I", name="z")
model.update()
model.addConstr(x + y + z == 32, "Heads")
model.addConstr(2 * x + 4 * y + 8 * z == 80, "Legs")
model.setObjective(y + z, GRB.MINIMIZE)
model.optimize()
print("Opt. Val.=", model.ObjVal)
print("(x,y,z)=", (x.X, y.X, z.X))
例題：多制約ナップサック問題
あなたは，ぬいぐるみ専門の泥棒だ． ある晩，あなたは高級ぬいぐるみ店にこっそり忍び込んで，盗む物を選んでいる．
狙いはもちろん，マニアの間で高額で取り引きされているクマさん人形だ． クマさん人形は，現在4体販売されていて， それらの値段と重さと容積は， 以下のようになっている．
- 極小クマ: 16, 2, 3000
- 小クマ: 19, 3, 3500
- 中クマ: 23, 5100
- 大クマ: 28, 5, 7200
あなたは，転売価格の合計が最大になるようにクマさん人形を選んで逃げようと思っているが， あなたが逃走用に愛用しているナップサックはとても古く， $7$ kgより重い荷物を入れると，底がぬけてしまうし， 10000 $cm^3$ を超えた荷物を入れると破けてしまう．
さて，どのクマさん人形をもって逃げれば良いだろうか？
- $n$ 個のアイテム，$m$本の制約
- 各々のアイテム $j=1,2,\ldots,n$ の価値 $v_j (\geq 0)$
- アイテム $j$ の制約 $i=1,2,\ldots,m$ に対する重み $a_{ij} (\geq 0)$
制約 $i$ に対する制約の上限値 $b_i (\geq 0)$
アイテムの番号の集合 $J =\{1,2,\ldots,n \}$
- 制約の番号の集合 $I =\{1,2,\ldots,m\}$
- アイテム $j$ をナップサックに詰めるとき $1$， それ以外のとき $0$ になる $0$-$1$ 変数 $x_j$
上の記号を使うと，多制約ナップサック問題の定式化は以下のようになる．
$$ \begin{array}{l l l} maximize & \displaystyle\sum_{j \in J} v_j x_j & \\ s.t. & \displaystyle\sum_{j \in J} a_{ij} x_j \leq b_i & \forall i \in I \\ & x_j \in \{0,1\} & \forall j \in J \end{array} $$J, v = multidict({1: 16, 2: 19, 3: 23, 4: 28})
a = {
(1, 1): 2,
(1, 2): 3,
(1, 3): 4,
(1, 4): 5,
(2, 1): 3000,
(2, 2): 3500,
(2, 3): 5100,
(2, 4): 7200,
}
I, b = multidict({1: 7, 2: 10000})
model = Model("mkp")
x = {}
for j in J:
x[j] = model.addVar(vtype="B", name=f"x({j})")
model.update()
for i in I:
model.addConstr(quicksum(a[i, j] * x[j] for j in J) <= b[i], f"Capacity({i})")
model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)
model.optimize()
print("Optimal value=", model.ObjVal)
EPS = 1.0e-6
for v in model.getVars():
if v.X > EPS:
print(v.VarName, v.X)
J, v = multidict({1: 16, 2: 19, 3: 23, 4: 28})
a = {
(1, 1): 2,
(1, 2): 3,
(1, 3): 4,
(1, 4): 5,
(2, 1): 3000,
(2, 2): 3500,
(2, 3): 5100,
(2, 4): 7200,
}
I, b = multidict({1: 7, 2: 10000})
model = Model("mkp")
x = {}
for j in J:
x[j] = model.addVar(vtype="B", name=f"x({j})")
model.update()
for i in I:
model.addConstr(quicksum(a[i, j] * x[j] for j in J) <= b[i], f"Capacity({i})")
model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)
# いずれかのコメントを外す
# 1.
# model.addConstr( x[2] + x[3] <=1 )
# 2.
# model.addConstr( x[2] <= x[1])
# 3.
model.addConstr(x[1] + x[2] + x[4] >= 2)
model.optimize()
print("Optimal value=", model.ObjVal)
EPS = 1.0e-6
for v in model.getVars():
if v.X > EPS:
print(v.VarName, v.X)
最大値の最小化
実数変数 $x_1, x_2$ に対する2つの線形関数 $3x_1+4x_2$ と $2x_1+7x_2$ の大きい方を小さくしたい場合を考える． これは最大値の最小化問題であり，新しい実数変数 $z$ を導入し， $$ 3x_1 + 4x_2 \leq z, \ \ 2x_1 + 7x_2 \leq z $$ の制約を加えた後で， $z$ を最小化することによって，通常の線形最適化に帰着できる． これは，制約が何本あっても同じである．また，最大値の最小化も同様に定式化できる．
整数変数に対する線形式に対しては，「最大値の最小化」をしたいときに上のアプローチをとることは推奨されない．（小規模な問題例は別である．） これは，混合整数最適化ソルバーで用いている分枝限定法が，このタイプの制約に対して弱く，限界値の改善が難しいためである．
例題： 以下の制約の下で， $3x_1+4x_2$ と $2x_1+7x_2$ の大きい方を最小化せよ．
$$ x_1 + 2 x_2 \geq 12 \\ 2 x_2 + x_2 \geq 15 $$model = Model("min-max")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
z = model.addVar(name="z")
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(3 * x1 + 4 * x2 <= z)
model.addConstr(2 * x1 + 7 * x2 <= z)
model.setObjective(z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
最小値の最小化
最小値の最小化はさらにやっかいである． 最大値の最小化の場合と同様に，新しい実数変数 $z$ を導入し，今度は最小値と一致するように以下の制約を加える．
$$ 3x_1 + 4x_2 \geq z, \ \ 2x_1 + 7x_2 \geq z $$これらの制約だけで $z$ を最小化すると目的関数値は $-\infty$ になる． これを避けるためには，上の制約において，左辺の小さい方と $z$ が等しくなるという条件が必要になる． 大きな数 $M$ と $0$-$1$ 変数 $y$ を用いると，この条件は以下の制約で記述することができる．
$$ 3x_1 + 4x_2 \leq z+ My, \ \ 2x_1 + 7x_2 \leq z + M(1-y) $$$y=0$ のときには，左側の制約だけが意味をもち，$3x_1 + 4x_2 \leq z$ となる． これを元の制約 $3x_1 + 4x_2 \geq z$ と合わせることによって，$3x_1 + 4x_2 = z$ を得る． $y=1$ のときには，右側の制約だけが意味をもち， $2x_1 + 7x_2 \leq z$ となる． これを元の制約 $2x_1 + 7x_2 \geq z$ と合わせることによって $2x_1 + 7x_2 =z$ を得る．
例題
以下の制約の下で， $3x_1+4x_2$ と $2x_1+7x_2$ の小さい方を最小化せよ．
$$ x_1 + 2 x_2 \geq 12 \\ 2 x_2 + x_2 \geq 15 $$model = Model("min-min")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
z = model.addVar(name="z")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(3 * x1 + 4 * x2 >= z)
model.addConstr(2 * x1 + 7 * x2 >= z)
model.addConstr(3 * x1 + 4 * x2 <= z + M * y)
model.addConstr(2 * x1 + 7 * x2 <= z + M * (1 - y))
model.setObjective(z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
絶対値の最小化
実数変数 $x$ の絶対値 $|x|$ を「最小化」したいという要求は実務でしばしば現れる． 「最小化」なら，変数を1つ追加するだけで処理できる． まず，$|x|$ を表す新しい変数 $z$ を追加し， $z \geq x$ と $z \geq -x$ の2つの制約を追加する． $z$ を最小化すると，$x$ が非負のときには $z \geq x$ の制約が効いて，$x$ が負のときには $z \geq -x$ の制約が効いてくるので，$z$ は $|x|$ と一致する．
別の方法もある．この方法は，後述するように絶対値の「最大化」にも拡張できる． 2つの新しい非負の実数変数 $y$ と $z$ を導入する． まず，変数 $x$ を $x=y-z$ と2つの非負条件をもつ実数変数の差として記述する． $x$ の絶対値 $|x|$ を $y+z$ と記述すると， $y \geq 0, z \geq 0$ の制約が付加されているので，$x$ が負のときには $z$ が正（このとき $y$ は $0$ ）となり， $x$ が正のときには $y$ が正（このとき $z$ は $0$）になる． つまり，定式化内の $x$ をすべて $y-z$ で置き換え，最小化したい目的関数内の $|x|$ をすべて $y+z$ で置き換えれば良い．
例題
以下の線形最適化問題を，上で説明した2通りの定式化で解け．
$$ \begin{array}{ll} minimize & 2 x_1 + 3 x_2 + |x_1-x_2| \\ s.t. &x_1 + 2 x_2 \geq 12 \\ & 2 x_2 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} $$model = Model("absolute(1)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
A = model.addVar(name="A", lb=-1000000000.0) # 差を表す実数変数
z = model.addVar(name="z") # |A|を表す補助変数
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == A)
model.addConstr(z >= A)
model.addConstr(z >= -A)
model.setObjective(2 * x1 + 3 * x2 + z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, A.X)
model = Model("absolute(2)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y") # x_1-x_2 = y-z
z = model.addVar(name="z")
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == y - z)
model.setObjective(2 * x1 + 3 * x2 + y + z, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, y.X, z.X)
絶対値の最大化
絶対値を「最大化」したい場合には，$y$ と $z$ のいずれか一方だけを正にすることができるという制約が必要になる． そのためには，大きな数 $M$ と $0$-$1$ 変数 $\xi$ を使い，以下の制約を付加する． $$ y \leq M \xi $$ $$ z \leq M (1-\xi) $$ 変数 $\xi$ が $1$ のときには，$y$ が正になることができ，$0$ のときには $z$ が正になることができる．
例題
以下の絶対値が入った線形最適化問題を解け．
$$ \begin{array}{ll} maximize & -2 x_1 - 3 x_2 + |x_1-x_2| \\ s.t. &x_1 + 2 x_2 \geq 12 \\ & 2 x_2 + x_2 \geq 15 \\ & x_1, x_2 \geq 0 \end{array} $$model = Model("absolute(3)")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y") # |x_1-x_2| = y-z
z = model.addVar(name="z")
xi = model.addVar(name="xi", vtype="B")
M = 999999.0
model.update()
model.addConstr(x1 + 2 * x2 >= 12)
model.addConstr(2 * x1 + x2 >= 15)
model.addConstr(x1 - x2 == y - z)
model.addConstr(y <= M * xi)
model.addConstr(z <= M * (1 - xi))
model.setObjective(-2 * x1 - 3 * x2 + y + z, GRB.MAXIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X, y.X, z.X)
離接制約
either or の条件 とは，いずれか1つが満たされていなければならないという条件のことである．
$2x_1 + x_2 \leq 30, x_1 + 2x_2 \leq 30$ の2 本の制約のいずれかが成立するという条件
- 大きな数 $M$ と $0$-$1$変数 $y$ を用いて，
$$2x_1 + x_2 \leq 30 +My, \ \ \ x_1 + 2x_2 \leq 30+M(1-y)$$
$2x_1 + x_2 \leq 30, \ \ \ x_1 + 2x_2 \leq 30, \ \ \ 5x_1+x_2 \leq 50$ の3 本の制約のいずれかが成立するという条件
- 3 つの$0$-$1$変数 $y_1, y_2, y_3$ を用いて，
$$2x_1 + x_2 \leq 30 +M(1-y_1), \ \ \ x_1 + 2x_2 \leq 30+M(1-y_2), \ \ \ 5x_1+x_2 \leq 50+M(1-y_3)$$
$$y_1 +y_2 + y_3 \geq 1$$
例題
以下の「もしくは」の上限が入った線形最適化問題を解け．
$$ \begin{array}{ll} maximize & x_1 + x_2 \\ s.t. & 2x_1 + x_2 \leq 30 もしくは x_1 + 2x_2 \leq 40 \\ & x_1, x_2 \geq 0 \end{array} $$model = Model("disjunctive")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(2 * x1 + x2 <= 30 + M * y)
model.addConstr(x1 + 2 * x2 <= 40 + M * (1 - y))
model.setObjective(x1 + x2, GRB.MAXIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
if A then B 条件
実際の問題を考える際には、 「制約Aが成立している場合には，制約Bも成立しなければならない」という論理制約がしばしば出てくる。例えば、「倉庫を建設しなければそこに保管することができない」とか「この作業とこの作業は同時刻に処理することはできない」などが代表的だ。ここでは、その取扱いについて解説する。
事象が真か嘘かは0-1変数で表現できる。0-1変数 $x,y$ に対して、「$x$が真ならば$y$も真である」という論理制約は、$x<=y$ で表現できる。 いずれかが真であるという論理制約は $x+y=1$ となる。 もう少し難しいものとして、「「この条件が成立しているときには、この制約が必要だ」（if ... then ... ）という論理制約がある。 「制約Aが成立している場合には，制約Bも成立しなければならない」という条件は， 「Aが成立しない or Bが成立する」と同値である。 ０−１変数x,yのいずれかが1になるという制約は $x+y=1$ であるので、「NOT A or B」は0-1変数を用いて記述することができる。
「if A then B」は「NOT A or B」と同値
| A | B | if A then B | NOT A | NOT A or B |
|---|---|---|---|---|
| 偽 | 偽 | 真 | 真 | 真 |
| 偽 | 真 | 真 | 真 | 真 |
| 真 | 偽 | 偽 | 偽 | 偽 |
| 真 | 真 | 真 | 偽 | 真 |
例として，整数変数 $x_1,x_2$ に対する以下の論理制約を考えよう。
- if.. then の条件：
$$if \ \ \ x_1 +x_2 \leq 10 \ \ \ then \ \ \ x_1 \leq 5 $$
この条件は，
$$x_1+x_2 \geq 11 \ \ \ or \ \ \ x_1 \leq 5 $$
と同値であるので， 0-1変数 $y$ と大きな数 $M$ を用いて
- 離接制約で表現：
と書くことができる． $x_1, x_2$ が整数なので， 第1式は以下と同値
$$x_1 + x_2 \geq 11 - My$$- $x_1,x_2$ が整数でなく実数変数の場合の問題点：
数理最適化ソルバーでは，より大きい（$>$）と以上（$\geq$）の制約を区別しない．
つまり， $x_1 + x_2 > 10$ は $x_1 + x_2 \geq 10$ と同値である． 微少な値 $\epsilon >0$ を用いて， 以下のように書くことができる．
$$x_1 + x_2 \geq 10 +\epsilon - My$$例題
以下の論理条件を含む整数最適化問題を解け．
$$ \begin{array}{ll} minimize & x_1 + x_2 \\ s.t. & x_2 \geq 2 \\ & if \ \ \ x_1 \leq 2 \ \ \ then \ \ \ x_2 \geq 6 \\ & x_1, x_2 は非負の整数 \end{array} $$model = Model("if_then")
x1 = model.addVar(name="x1")
x2 = model.addVar(name="x2")
y = model.addVar(name="y", vtype="B")
M = 99999
model.update()
model.addConstr(x2 >= 2)
model.addConstr(x1 >= 3 - M * y)
model.addConstr(x2 >= 6 - M * (1 - y))
model.setObjective(x1 + x2, GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
print(x1.X, x2.X)
# 1. x=1 <=> x+y=0
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z = model.addVar(name="z", vtype="B")
# if x=1, then x+y=0
# x=0 or x+y =0
model.addConstr(x + (x + y) <= 1)
# if x+y=0, then x=1
# x+y>=1 or x>=1
model.addConstr(x + y >= 1 - 9999 * z)
model.addConstr(x >= 1 - 9999 * (1 - z))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X) # 旦那は嘘つき，奥さんは正直
# 2. x=1 <=> x+y<=1
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z1 = model.addVar(name="z1", vtype="B")
z2 = model.addVar(name="z2", vtype="B")
# if x=1, then x+y<=1
# x<=0 or x+y <=1
model.addConstr(x <= 9999 * z1)
model.addConstr(x + y <= 1 + 9999 * (1 - z1))
# if x+y<=1, then x=1
# x+y>=2 or x>=1
model.addConstr(x + y >= 2 - 9999 * z2)
model.addConstr(x >= 1 - 9999 * (1 - z2))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X) # 旦那は正直，奥さんは嘘つき
# 3. x=1 <=> (if x=1, then y=1) => x<=y
model = Model()
x = model.addVar(name="x", vtype="B")
y = model.addVar(name="y", vtype="B")
z = model.addVar(name="z", vtype="B")
# if x=1, then x<=y
model.addConstr(x <= y)
# if x<=y, then x=1
# x>=y+1 or x>=1
model.addConstr(x >= y + 1 - 9999 * z)
model.addConstr(x >= 1 - 9999 * (1 - z))
model.setObjective(x, GRB.MAXIMIZE)
model.optimize()
print(x.X, y.X) # 両方とも正直
問題 論理パズル (1)
ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる． 正直族は必ず本当のことを言い，嘘つき族は必ず嘘をつく．
またこの島には，夜になると狼に変身して人を襲う狼男が紛れ込んでいる． 狼男もこの島の住民なので，正直族か嘘つき族のいずれかに属する． あなたは，この島の人たちが狼男なのかの調査を依頼された．
3人のうち1 人が狼男であることが分かっている． A,B,C の3人組への証言は以下の通りである．
- A: 「わたしは狼男です．」
- B: 「わたしも狼男です．」
- C: 「わたしたちの中の高々1人が正直族です．」
さて，狼男は誰か？また誰が正直族で誰が嘘つき族か？
問題 絶対値
7 組の家族がみんなで食事会をしようと考えている． 4 人がけのテーブルが7 卓であり，親睦を深めるために，同じ家族の人たちは同じテーブルに座らないようにしたい．
家族の構成は以下の通りとしたとき，各テーブルの男女比をなるべく均等にする（女性と平均人数との差の絶対値の和を最小化する）座り方を求めよ． ただし，女性には (F) の記号を付してある．
- 磯野家: 波平，フネ (F)，カツオ，ワカメ(F)
- バカボン家: バカボンパパ，バカボンママ (F)，バカボン，ハジメ
- 野原家: ひろし，みさえ(F)，しんのすけ，ひまわり(F)
- のび家: のび助，玉子(F)，のび太
- 星家: 一徹，明子(F)，飛雄馬
- レイ家: テム，カマリア(F)，アムロ
ザビ家: デギン，ナルス(F)，ギレン，キシリア(F)，サスロ，ドズル，ガルマ
また、男女比をなるべく均等にしない（女性と平均人数との差の絶対値の和を最大化する）座り方を求めよ。
例題：多期間生産計画
あなたの会社では，2つの製品(P1, P2)を2台の機械(M1, M2)で製造している． 今後3ヶ月の需要が決まっており，製品P1に対しては 4000,8000,3000， 製品P2に対しては 1000,5000,5000である．
各月の機械を使用できる時間(h)も決まっており，機械M1に対しては700,300,1000，機械M2に対しては1500,400,300である．
製品は機械によって製造速度が異なり， 製品を1単位製造するのに必要な時間(h)は，以下の表のように与えられているものとする．
| M1 | M2 | |
|---|---|---|
| P1 | 0.15 | 0.16 |
| P2 | 0.12 | 0.14 |
製造の変動費用は，製品や機械によらず1時間あたり5ドル，製品1単位を1ヶ月持ち越したときの生じる在庫費用は0.1ドルとする． なお，初期在庫は0，最終在庫も0と仮定する．
最適な生産計画を立案せよ．
最終期に機械M1の使用可能な時間を100時間増やせる設備を100ドルで購入可能である． 購入すべきか？
機械のメンテナンスをするなら，どの月にどの機械に対してやるべきか？
from mypulp import Model, GRB, quicksum
T = 3
products = ["P1", "P2"]
machines = ["M1", "M2"]
demand = {"P1": [4000, 8000, 3000], "P2": [1000, 5000, 5000]}
capacity = {"M1": [700, 300, 1000], "M2": [1500, 400, 300]}
rate = {("P1", "M1"): 0.15, ("P1", "M2"): 0.16, ("P2", "M1"): 0.12, ("P2", "M2"): 0.14}
variable_cost = 5.0
inventory_cost = 0.1
model = Model()
x, I = {}, {}
for p in products:
for m in machines:
for t in range(T):
x[p, m, t] = model.addVar(name=f"x({p},{m},{t})")
for p in products:
I[p, -1] = I[p, T - 1] = 0.0
for t in range(T - 1):
I[p, t] = model.addVar(name=f"I({p},{t})")
model.update()
capacity_constraints = {}
for m in machines:
for t in range(T):
capacity_constraints[m, t] = model.addConstr(
quicksum(rate[p, m] * x[p, m, t] for p in products) <= capacity[m][t],
name=f"Capacity Constraint({m},{t})",
)
demand_constraints = {}
for p in products:
for t in range(T):
demand_constraints[p, t] = model.addConstr(
I[p, t - 1] + quicksum(x[p, m, t] for m in machines)
== demand[p][t] + I[p, t]
)
model.setObjective(
quicksum(
variable_cost * rate[p, m] * x[p, m, t]
for p in products
for m in machines
for t in range(T)
)
+ quicksum(inventory_cost * I[p, t] for p in products for t in range(T)),
GRB.MINIMIZE,
)
model.optimize()
print("Obj. Val.=", model.ObjVal)
for t in range(T):
for p in products:
for m in machines:
print(p, m, t, x[p, m, t].X)
余裕変数（slack)や双対変数(Pi)をみると，機械M1の1時間あたりの価値は0.33ドルなので，購入すべきでない．
余裕変数をみるとメインテナンスは，最初の月か最後の月にすべきである．
for c in capacity_constraints:
print(c, capacity_constraints[c].slack, capacity_constraints[c].Pi)
問題 2段階製造問題
ある会社が，3つの製品A,B,Cを製造している。 製造工程には、「切断」と「プレス」という2つの工程がある． すべての製品はこの2つの瞬間を通過しなければならない。
1日8時間使用可能な切削部門の生産能力は，製品Aは1時間に2000個、製品Bは1時間に1600個，製品Cは1時間に1100個である．
1日8時間使用可能なプレス部門の生産能力は，製品Aは1時間に1000個、製品Bは1時間に1500個、製品Cは1時間に2400個である．
製品切り替えの段取り時間（費用）はかからないものとする．
製品A,B,Cの1個あたりの利益をそれぞれ12,9,8（万円）としたとき，利益を最大化する製造計画を立案せよ．
問題 ナプキンのクリーニング（難）
あなたはホテルの宴会係だ．あなたは1週間に使用するナプキンを手配する必要がある． 各日の綺麗なナプキンの需要量は平日は $100$枚，土曜日と日曜日は $125$枚だ． 新しいナプキンを購入するには $100$円かかる． 使用したナプキンはクリーニング店で洗濯して綺麗なナプキンにすることができるが， 早いクリーニング店だと1日で1枚あたり $30$円かかり， 遅いクリーニング店だと2日で1枚あたり $10$円かかる． 月曜の朝のナプキンの在庫が $0$ としたとき，需要を満たす最適なナプキンの購入ならびにクリーニング計画をたてよ．
さらに，日曜末の在庫を月曜の朝に使うことができると仮定したときの最適なナプキンのクリーニング計画をたてよ．
機械学習と数理最適化の融合
Gurobi 10以降だと，機械（深層）学習モデルを数理最適化モデル内で記述できるようになった．
パッケージとしては，Python 3.9以上でインストールできる gurobi-machinelearning (https://github.com/Gurobi/gurobi-machinelearning) が必要となる．
簡単な例題で説明する．
$n$ 人の学生 $i$ に対して，奨学金 $x_i$ と成績 $SAT_i, GPA_i$ から入学する確率 $y_i$ をロジスティック回帰 $g$ で予測する．
予算の上限の下で，入学者数を最大にする最適化モデルを考える．
学生数の2割の総予算と，1人あたりの奨学金の上限を $2.5$ とする．
問題は以下のように定式化できる．
$$ \begin{array}{ ccc } maximize & \sum_{i=1}^n y_i \\ s.t. & \sum_{i=1}^n x_i \le 0.2 n \\ & y_i = g(x_i, SAT_i, GPA_i) & i = 1, \ldots, n \\ & 0 \le x_i \le 2.5 & i = 1, \ldots, n \end{array} $$import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
import gurobipy as gp
from gurobi_ml import add_predictor_constr
janos_data_url = "https://raw.githubusercontent.com/INFORMSJoC/2020.1023/master/data/"
historical_data = pd.read_csv(
janos_data_url + "college_student_enroll-s1-1.csv", index_col=0
)
features = ["merit", "SAT", "GPA"]
target = "enroll"
historical_data
regression = LogisticRegression(random_state=1)
pipe = make_pipeline(StandardScaler(), LogisticRegression(random_state=1))
pipe.fit(X=historical_data.loc[:, features], y=historical_data.loc[:, target])
studentsdata = pd.read_csv(janos_data_url + "college_applications6000.csv", index_col=0)
nstudents = 250
studentsdata = studentsdata.sample(nstudents)
studentsdata
feat_lb = studentsdata.copy()
feat_lb.loc[:, "merit"] = 0
feat_ub = studentsdata.copy()
feat_ub.loc[:, "merit"] = 2.5
feat_lb = feat_lb[features]
feat_ub = feat_ub[features]
feat_ub
m = gp.Model()
#行列形式で変数（特徴）を生成する．
feature_vars = m.addMVar(
feat_lb.shape, lb=feat_lb.to_numpy(), ub=feat_ub.to_numpy(), name="feats"
)
y = m.addMVar(nstudents, name="y")
#数理最適化の変数だけを切り出す
x = feature_vars[:, feat_lb.columns.get_indexer(["merit"])][:, 0]
m.setObjective(y.sum(), gp.GRB.MAXIMIZE)
m.addConstr(x.sum() <= 0.2 * nstudents)
m.update()
pred_constr = add_predictor_constr(
m, pipe, feature_vars, y, output_type="probability_1"
)
pred_constr.print_stats()
m.optimize()
print(
"Maximum error in approximating the regression {:.6}".format(
np.max(pred_constr.get_error())
)
)
print("Obj. Val.=", m.ObjVal)
for i,val in enumerate(x.X):
if val>0.001:
print(i, val, y[i].X)