Pythonによる最適化　第3回

前回はPythonのモジュールnetworkXを紹介した．今回は，科学技術計算用モジュールSciPyを用いた問題解決について解説する．例とするのは，確率的在庫モデル（新聞売り子モデル）と呼ばれる古典的な問題である．問題は以下のように記述できる．

新聞の売り子が，1種類の新聞を販売している．新聞が売れる量（需要量）は，経験からある程度推測できると仮定し，確率変数として与えられているものとする．いま，売れ残りのときの在庫費用hと，品切れのときの品切れ費用bの和が最小になるように仕入れ量xを決めるものとする．どれだけの量を仕入れれば良いだろうか?

新聞の売れ行きが正規分布（負の場合には0）と仮定してGurobiで解いてみよう．まず，ランダムな需要を発生させる．通常のrandomモジュールのgaussメソッドを用いて正規分布にしたがう需要を発生させるのも手であるが，ここではもう少し厳密に需要iが発生する確率prob[i]をSciPyモジュールを用いてか発生させよう．まず，SciPyの正規分布に関するクラスnormを読み込む．

from scipy.stats import norm 

次いで，平均mu，標準偏差sigmaの正規分布のインスタンスを生成する．

mu=100   #需要の平均
sigma=10 #需要の標準偏差
ND=norm(mu,sigma)

正規分布オブジェクトのcdfメソッドは累積分布関数であり，引数以下になる確率を返す．したがって，需要dが0以下になる確率は以下のように設定すれば良い．

prob[0]=ND.cdf(0)
d[0]=0

同様にpdfメソッドは確率密度関数であり，引数になる確率を返す．したがって需要がi=1..Sになる確率を以下のように設定する．ちなみにSはシナリオの数でここでは1000に設定してある．

for i in range(1,S):
    d[i]=i
    prob[i]=ND.pdf(i)

このように設定したパラメータを用いると新聞売り子問題は，以下のようにGurobiで定式化できる．ここで，xが発注量，B[s],I[s]はシナリオsにおける品切れ量と在庫量，b,hはそれぞれ品切れ，在庫費用である．

x=model.addVar(vtype="C",name="x”)
for s in range(S):
    B[s]=model.addVar(vtype="C",name="B(%s)"%(s))
    I[s]=model.addVar(vtype="C",name="I(%s)"%(s))
model.update()
model.addConstr(x+B[s]==d[s]+I[s])

model.setObjective(quicksum(prob[s]*(h*I[s]+b*B[s]) for s in range(S) ),
GRB.MINIMIZE)

実は新聞売り子問題に対しては解析的な解を計算することができる．しかし，他の分布やリスクを考慮したモデル，CVaRを評価関数としたモデル，多期間や多段階の在庫モデルなど，様々なバリエーションに対しても上の定式化は容易に拡張できるのが利点である．実際に，筆者の近著（サプライ・チェインにおけるリスク管理と人道支援ロジスティクス）では，上のモデルの拡張に対して実験的解析を行い，リスクを考慮した在庫モデルに対する様々な洞察を得ている．