このブログis何
Googleが開発したOR-Toolsを使って数理最適化を学んでいきます。
内容は殆どがこの本を参考にしたものです。
Practical Python AI Projects: Mathematical Models of Optimization Problems with Google OR-Tools (English Edition)
著者によるコード一覧はこちら
本ブログに記載するコードは著者オリジナルの物と殆ど変わらず、一部日本語に直したりしている程度です。
数学的な厳密さや、理論の解説は優先度低めになってます。ご勘弁を。
今回は : 線形計画～製油編～
線形計画のモデリングをします。線形計画ってなんやねんとかはサボります、マジごめんなさい。
今回は「原料となる原油から欲しい石油製品を良い感じに作って利益を最大化しようぜ
」って例題です。
原油の使用量なら変数が連続値でも許されそうだからこんな例題なんやと思います、多分。
場面設定
あなたは夏休みの宿題である自由研究に取り組む小学生です。社会科見学で製油所に行ったあなたは、「効率よく石油製品を作って儲けてえ…」と思い、これを自由研究の題材に選びました。
作りたい石油製品にはそれぞれオクタン価の指定があり、原料となる原油にもそれぞれオクタン価があります。(オクタン価はノッキングの起こりにくさだそうです。調べてもよく分からなかったので勝手にアルコール度数みたいなもんだと考えてます)
オクタン価の異なる原油を良い感じに混ぜて、特定のオクタン価をもつ石油製品を作らねばなりません。原油は種類によって1バレルあたりの費用が異なり、石油製品は種類によって1バレルあたりの売価が異なります。
文章にすると意味わからんくなってくるので表をご覧ください。
| 原油 | オクタン価 | 所有量 | 費用/バレル |
|---|---|---|---|
| R0 | 99 | 782 | 55.34 |
| R1 | 94 | 894 | 54.12 |
| R2 | 84 | 631 | 53.68 |
| R3 | 92 | 648 | 57.03 |
| R4 | 87 | 956 | 54.81 |
| R5 | 97 | 647 | 56.25 |
| R6 | 81 | 689 | 57.55 |
| R7 | 96 | 609 | 58.21 |
| 製品 | オクタン価 | 需要下限 | 需要上限 | 売価 |
|---|---|---|---|---|
| F0 | 88 | 415 | 11707 | 61.97 |
| F1 | 94 | 199 | 7761 | 62.04 |
| F2 | 90 | 479 | 12596 | 61.99 |
ここで重要な仮定として、原油を混ぜたもののオクタン価は線形関数で決まるとします。
つまり、オクタン価が80の原油と90の原油を4:6の比率で混ぜたときには、
80\times\frac{4}{10}+90\times\frac{6}{10}=86
と、オクタン価86の製品が生まれます。食塩水の濃度の問題とかでやった気がしますね。
さて、あなたは利益を最大化したいわけですから、(作った石油製品による売上)-(使った原油による費用)を最大にする必要があります。
定式化
「利益を最大化するような〇〇」を知りたいというのが今回の問題ですが、〇〇はなんでしょうか？
「それぞれの原油の使用量」だと少し足りません。ある石油製品Xに注目したとき、Xを作るのにどの原油をどれだけ使ったかの内訳が分からないと困るからです。
そこで、
G_{ij}:原油iを、製品j のために使う量
という変数を用意します。したがって（原油の種類）×（製品の種類）の行列として決定変数が定義されます。この表の中身が分かればハッピー。
これにより、原油の総使用量Rと製品の総生成量Fについて、
R_i = \sum_{j} G_{ij}:原油i の総使用量。行方向に和をとる\\
F_j = \sum_{i} G_{ij}:製品j の総生成量。列方向に和をとる
が成り立ちます。
さあこれで目的関数を定義できます。原油iの単価を$c_i$、製品jの売価を$p_j$とすると、
max \sum_{j}F_jp_j-\sum_{i}R_ic_i
が目的関数になります。
最後に制約条件ですが、$R_i,F_j$の取りうる範囲に関しての制約は面倒くさいので省略します。まあ所有量以上は使えないとか、需要の範囲で作らないといけないとかです。
今回大事なのはオクタン価の話です。原油はただ混ぜるのではなく、生成物が特定のオクタン価を持つように混ぜねばなりません。
原油iのオクタン価を$O_i$、製品jのオクタン価を$o_j$とすると、（既にややこしい）
\sum_{i}O_iG_{ij}=F_jo_j\\
が成り立つ必要があります。これが満たされなければ、欲しい製品を作れてないことになります。
とりあえずこれで定式化は出来ました。それでは実装しましょう。
実装
今回は必要な関数を記述したコードのみ紹介します。著者GitHubに必要なコードは全て載っています。
実行するのにtest_gas_blend.py、OR-Toolsを使いやすくするためにmy_or_tools.pyを用います。
from my_or_tools import SolVal,ObjVal,newSolver
def solve_gas(C, D): #原油のリスト,製品のリストを引数にとる
s = newSolver('Gas blending problem') #ソルバーを定義
nR,nF = len(C),len(D) #原油の数,製品の数を取得
Roc,Rmax,Rcost = 0,1,2 #原油のオクタン価,所有量,費用の列番号
Foc,Fmin,Fmax,Fprice = 0,1,2,3 #製品のオクタン価,需要下限,需要上限,売価の列番号
G = [[s.NumVar(0.0,10000,'')
for j in range(nF)] for i in range(nR)] #変数Gijの行列を用意
R = [s.NumVar(0,C[i][Rmax],'') for i in range(nR)] #Gijの行方向の和である変数Rを用意
F = [s.NumVar(D[j][Fmin],D[j][Fmax],'') for j in range(nF)] #Gijの列方向の和である変数Fを用意
for i in range(nR):
s.Add(R[i] == sum(G[i][j] for j in range(nF))) #Rが満たすべき等式
for j in range(nF):
s.Add(F[j] == sum(G[i][j] for i in range(nR))) #Fが満たすべき等式
for j in range(nF):
s.Add(F[j]*D[j][Foc] == #オクタン価に関する制約
s.Sum([G[i][j]*C[i][Roc] for i in range(nR)]))
Cost = s.Sum(R[i]*C[i][Rcost] for i in range(nR)) #原油の使用による費用
Price = s.Sum(F[j]*D[j][Fprice] for j in range(nF)) #製品の生成による売上
s.Maximize(Price - Cost) #利益を最大化する
rc = s.Solve()
return rc,ObjVal(s),SolVal(G)
my_or_tools.pyを使うことでソルバーの定義などが楽になっています。
結果
右下の値がこの最適解による利益となっています。
| 製品0 | 製品1 | 製品2 | 使用量 | 費用 | |
|---|---|---|---|---|---|
| 原油0 | 542.5 | 239.5 | 782.0 | 43275.88 | |
| 原油1 | 894.0 | 894.0 | 48383.28 | ||
| 原油2 | 631.0 | 631.0 | 33872.08 | ||
| 原油3 | 648.0 | 648.0 | 36955.44 | ||
| 原油4 | 704.41 | 251.59 | 956.0 | 52398.36 | |
| 原油5 | 647.0 | 647.0 | 36393.75 | ||
| 原油6 | 449.5 | 239.5 | 689.0 | 39651.95 | |
| 原油7 | 50.93 | 558.07 | 609.0 | 35499.89 | |
| 生成量 | 2378.83 | 2988.67 | 479.0 | ||
| 売上 | 147385.32 | 186037.28 | 29693.21 | 36375.18 |
たかし君は無事夏休みの自由研究を提出することができました。めでたしめでたし。
Qiitaで数式や表を書くのは大変やな。慣れてないだけかもしれんけど。
自分の復習用なんでアルゴリズムの解説はほとんどなく、定式化とOR-Toolsによる実装だけになりそうです、すんません。