最適化ズブの素人がGEKKOを使って輸送問題を解く

はじめに

皆さんはGEKKOを知っているだろうか？はじめに言っておくとGEKKOは最適化問題を計算機で解くための"ソルバー"と呼ばれるツールを簡単に動かすためのツールである（こういうものをモデラーと呼ぶらしい）。pythonで動くのでインタープリタの強みを活かしてサクサク最適化問題を解くことができる。
僕とGEKKOの出会いは突然だった。それは、会社での出来事。業務の中で最適化問題をかじる事になった時。隣のデスクの最適化おじさん（※最適化問題を愛しているおじさんの事）が呟いた。
「GEKKOってのがあるらしいよ」
なにそれ、かっこいい。。。
僕の脳裏には某ロボアニメのテロリスト集団「ゲッコー○テイト」のメンバーの顔が走馬灯のように駆け抜けた。
実はソルバーやモデラーは他にもいっぱい種類があるのだが、僕はGEKKOを使うことに決めた。理由は名前がかっこいいから。

ここには触ってみた感想とかを適当に綴っておく。

制約付き最適化問題と輸送問題

GEKKOの話の前に一般的な制約付き最適化問題やその中の輸送問題について書く。実は素人なので、間違っている内容もあるかもしれない。

制約付き最適化問題とは

簡単にいうと、いくつかの制約条件を満たしつつ、目的関数を最小（あるいは最大）化する変数を探しだす問題である。
線形・非線形、連続・不連続などで多少勝手が変わるが、基本的には以下のように定式化される。

minimize~~f(x) \\
subject~to~~x \in X

$x \in X$で表現してあるのが制約式にあたる。ここでは、等号や不等号で変数を縛る条件が付与される。
さらに詳しくは、ネットに資料や記事が溢れているので割愛する。

最適輸送問題

読んで字の如く、輸送を最適化する。

例えば、上の図のように、東京と大阪にそれぞれ倉庫があり、同種の在庫がストックしてあるとする。福岡のAさんが商品を頼んだ時、単純に考えるなら距離が近い大阪の倉庫から発送されるだろう。これは、輸送コストに応じて輸送を最適化している事になる。
では、輸送元も先も複数ある場合はどうなってくるだろうか。考えられる輸送ルートは輸送元数×輸送先数存在する事になる。そして、各ルートそれぞれが異なる輸送コストを持つ。その上でトータルの輸送コストを最小に抑える輸送量の組み合わせを探さなくてはならない。
仮に、ある輸送先だけに注目するなら話は簡単である。その輸送先に延びるルートのうちコストが一番低いものを選ぶだけでよい。じかし実際には、輸送量には幾つかの制約があるのでこれはできない。例えば、倉庫は無限の在庫を格納できるわけではない。在庫の兼ね合いによっては、少しコストの高いルートも選択されうることになる。最小化したいのはあくまで全体のコストなのだ。
ある条件の基、何かを最小化する。というわけで、輸送問題は条件付き最適化問題に落とし込むことができる。

GEKKO

概要

最適化モデラーである。名前がかっこいい。詳しくは公式ドキュメントを参照されたし。

インストール

インストールは簡単で、下の１行で済んだ。
pip install gekko
実はここがGEKKOのいいところ。ソルバーを別でビルドしたりせずに、この１行だけで準備完了となる。

実際にやってみた

適当に考えた練習問題をGEKKOで解いてみた。

問題１

倉庫から消費者への米の輸送を考える。

今回は以下のような設定をおく。

• 倉庫はそれぞれキャパが異なる
• 消費者はそれぞれ消費量が異なる
• 倉庫-消費者間の輸送経路は決まっており、輸送コストが与えられる
• 供給の総量=需要の総量とする

その上で、倉庫-消費者間の各ルートに何キロの米を流せばいいかを求めたい。

定式化

この問題は連続変数の線形最適化問題に落とし込める。
まず、各倉庫にある米の量(供給)をベクトル$ \boldsymbol{s} $、消費者が求めている米の量(需要)をベクトル$\boldsymbol{d}$で表す。
そして、$s_i$から$d_j$への輸送量を$\gamma_{ij}$、輸送コストを$c_{ij}$とする。
この問題では、需要と供給の送料は同じとしているので、テーブルにまとめると以下のようになる。ちなみに、輸送量は常に正であることに注意。

ここで、最適化関数をモデリングする。
トータルのコストは単純に輸送量と輸送コストの線形和だと考えて
$$\sum_{i,j} \gamma_{ij} \cdot c_{ij}$$
とする。
以上より、以下の様な最適化問題に定式化される。

minimize~~\sum_{i,j} \gamma_{ij} \cdot c_{ij} \\
subject~~to \sum_j \gamma_{ij} = s_i \\
~~~~~~~~~~~\sum_i \gamma_{ij} = d_j \\
~~~~~~~~~~~\gamma_{ij} \geq 0

実装

GEKKOで実装していく。コード全文はコチラ。
まずは必要なライブラリのインポート。

from gekko import GEKKO
import numpy as np

今回使うのはたったこれだけである

続いて、前提として与えられる定数を定義しておく。数は適当に決めている。

#### 問題設定 ####
ware_house_num = 3 # 倉庫の数
shop_num = 2 # 客の数

supry = np.array([10,5,2]) # 倉庫の方の供給
demand = np.array([7,10]) # 客の方の需要
c_npa = np.array([[1.,2.,3.], [4.,5.,6.]]) # 輸送コスト

さらに、この手順が本当に必要か定かではないが、需要と供給のベクトルをGEKKOのm.Paramという型に変換しておく

supry_pa = m.Array(m.Param,(supry.shape))
for i in range(len(supry)):
    supry_pa[i].value = supry[i]

demand_pa = m.Array(m.Param,(demand.shape))
for j in range(len(demand)):
    demand_pa[j].value = demand[j]

また、コスト行列は定数なので,m.Constという型に変換する

c = m.Array(m.Const,(c_npa.shape))
for i in range(shop_num):
  for j in range(ware_house_num):
    c[i,j].value = c_npa[i, j]

そして、変数を初期化する。この値も適当。※ところで、上の定式化で使ってる文字とプログラムの変数名が揃ってないのはご容赦を…

# 変数の初期化
s = m.Array(m.Var,(shop_num, ware_house_num,))
for i in range(shop_num):
  for j in range(ware_house_num):
    s[i,j].value = 2.0
    s[i,j].lower = 0.0

ここで、変数の値と下限を設定している。この他にも、上限や整数縛り、好きな名前を付けたりすることもできる。

一番重要な目的関数と制約式は次のように書く

# 最適化関数
m.Obj(m.sum([m.sum([s[i,j] * c[i,j] for i in range(shop_num)]) for j in range(ware_house_num)]))
# 制約式
for i in range(shop_num):
    m.Equation(m.sum([s[i,j] for j in range(ware_house_num)]) == demand_pa[i])
for j in range(ware_house_num):
    m.Equation(m.sum([s[i,j] for i in range(shop_num)]) == supry_pa[j])

ここまで書けたらあとは、ワンポチでソルバーが面倒な最適化アルゴを自動で実行してくれる。具体的には次のメソッドを呼ぶ。

m.solve()

エラー無くだらだらと計算過程のログなどが出力されていたらおそらく成功である。
ログの下の方のObjective :と書かれた欄には目的関数の値が示されているが、今回注目したいのはどちらかというと、その値を取る時の変数の値なので、それを確認したかったら上で定義した変数を出力するといい。この場合はprint(s)などとすれば確認できる。

問題２

さらに少し複雑な問題を考えてみる。米は農家から倉庫、倉庫から消費者へ渡ると想定する。また、倉庫では何故か一定量中抜きが行われているとする。この時の輸送量を最適化する。

定式化

農家の供給量を$\boldsymbol{s}$、消費者の需要量を$\boldsymbol{d}$倉庫での中抜き量を$\boldsymbol{a}$とする。また農家から倉庫への輸送量と輸送コストを$\bar{\gamma_{kj}}, \bar{c_{kj}}$、倉庫から消費者へのそれを$\gamma_{ji}, c_{ji}$とすると以下のように定式化される。

minimize~~\sum_{k,j} \bar{\gamma_{kj}} \cdot \bar{c_{kj}} + \sum_{j,i} \gamma_{ji} \cdot c_{ji} \\
subject~~to \sum_j \bar{\gamma_{kj}} = s_k \\
~~~~~~~~~~~\sum_i \gamma_{ji} = \sum_k \bar{\gamma_{kj}} - a_j \\
~~~~~~~~~~~\sum_j \gamma_{ji} = d_i \\
~~~~~~~~~~~\bar{\gamma_{kj}},\gamma_{ji} \geq 0

実装

問題１とそんなに変わらないので説明は省略。コードはコチラを参照されたし

感想

定式化した問題を簡単にコードに落とすことができた。最初は「あんまり使ってる人いなさそうだし、大丈夫か？」と思ったが、なかなかいいおもちゃを手に入れたかもしれない。
さながら『ねだるな、勝ち取れ、さすれば与えられん』といったところか（言いたいだけ）。

おわりに

久しぶりに記事を書いた。どうやら前記事を書いた時から３年以上も経過しているらしい、、、３年てやばいな

そして、定式化や実装面で上述の最適化おじさんに大変有益なアドバイスをいくつも頂きました。ありがとう最適化おじさん。。。

追記

その後、少し使ってみるとGEKKOの遅さに気づくことになる。実は、GEKKOの裏で走っているIPOPTというソルバーは非線形問題用で、上の様な線形問題を解くときは線形問題用のソルバーを使った方が速く解けるらしい。というわけでもうGEKKOは使わない。フォーエバーGEKKO...