数理最適化モデリング言語 PySIMPLE の紹介

この記事は「数理最適化 Advent Calendar 2020」 5 日目の記事です．
4 日目は @kievdhia さんによる「確率計画法とその周辺の紹介」でした．

0. はじめに

はじめまして，NTT データ数理システム の数理計画部でシニアリサーチャーをしている池田です．普段は，じゃがりこを並べたり，数理最適化モデリング言語の開発をしています．

要するに中の人なのですが，皆さまは数理最適化の仕事と言えば何を思い浮かべるでしょうか．おそらく，ソルバー開発や，ソルバーを用いた問題解決といったことが真っ先に出てくるはずです．実際は，上記に関わるさまざま仕事が存在し，私がしている数理最適化モデリング言語の開発なども存在します．今回は，こちらの紹介をしていきます．

1. 数理最適化モデリング言語

数理最適化などの数学的なアルゴリズムを実現するソフトウェアに，具体的な数理モデルを与えるには，ユーザが数理モデルをプログラムとして記述する必要があります．

ソルバーの API は，係数行列を与えるなど，しばしば数式で記述された定式化と乖離が発生します．そこで，定式化された数理モデルをより自然に記述するインターフェースとして，数理最適化問題を記述する専用のモデリング言語があると便利です．

数理最適化モデリング言語には，PuLP，Pyomo，AMPL など様々な種類がありますが，私は PySIMPLE というものを開発しています．

例えば，1 日目の @___monta___ さんによる「感動した！ min-max max-min の最適化表現 手法」の定式化は PySIMPLE を用いると，次のように記述することができます．

from pysimple import *

###### 定数定義
# 格子の数
SIZE = 10
# ばらまくPointの数
N = 5

# 添字
i = Element(value=range(SIZE))
j = Element(value=range(SIZE))
k = Element(value=range(SIZE))
l = Element(value=range(SIZE))

# 2 点間の距離を算出する
distvalue = {(x1, x2, y1, y2): sqrt((x1-x2)**2+(y1-y2)**2)/SIZE for x1, x2, y1, y2 in product(*[range(SIZE)]*4)}
dist = Parameter(index=(i,k,j,l), value=distvalue)

###### 変数定義
# どこに Point を置くか
x = BinaryVariable(index=(i,j))
# 今回のキモ 多目的のための変数(どんなに良くても対角の sqrt(2) 以下)
z = Variable(lb=0, ub=1.5)

# 問題定義
p = Problem(name='SpreadingPointsProblem', type=max)

# z の最大化問題として定義
p += z

p += Sum(x[i,j], (i,j)) == N, 'N 個の Point が配置される'

# 確認した 2 点に同時に点が配置されているとしたばあい、その距離は z 以上
BIGM = z.ub
ikjl = (i!=k) | (j!=l)
p += z <= (1-x[ikjl(0,2)])*BIGM + (1-x[ikjl(1,3)])*BIGM + dist[ikjl], '2 点の距離は z 以上'

p.solve(silent=True)
assert p.status == NuoptStatus.OPTIMAL
print(p.objective.val)

# 結果出力
print(x[x[i,j].val==1].val)  # 1 が立っている部分だけ
Printf(','.join(['{:.0f}']*SIZE), *(x[i,j_].val for j_, in j.set))  # 表形式で

2. 輸送問題

LP の例として次の輸送問題を考えてみましょう．

構成要素	名称	説明
添字	$d \in \{d0, d1, .. \}$	倉庫
添字	$c \in \{c0, c1, .. \}$	顧客
定数	$cost_{d, c}$	輸送コスト
定数	$upper_d$	倉庫上限
定数	$lower_c$	顧客需要
変数	$z_{d,c}$	輸送量
目的関数	$\sum_{d,c} cost_{d,c}\cdot z_{d,c}$	輸送コスト
制約式	$\sum_c z_{d,c} \le upper_d$	倉庫上限
制約式	$\sum_d z_{d,c} \ge lower_c$	顧客需要

この定式化を PySIMPLE を使ってモデリングすると，次のようになります．

### データ用意
from random import seed, randint

DSIZE, CSIZE = 2, 3

class value:
    seed(0)
    upper = {f'd{d}': randint(1, 100000) for d in range(DSIZE)}
    lower = {f'c{c}': randint(1, 10000) for c in range(CSIZE)}
    cost = {(f'd{d}',f'c{c}'): randint(1, 100) for d in range(DSIZE) for c in range(CSIZE)}


### 最適化
from pysimple import Element, Parameter, Variable, Sum, Problem, NuoptStatus, Printf

d = Element(value=value.upper.keys()) # 倉庫
c = Element(value=value.lower.keys()) # 顧客

cost = Parameter(index=(d,c), value=value.cost)  # 倉庫から顧客への輸送コスト
upper = Parameter(index=d, value=value.upper)    # 倉庫取扱量上限
lower = Parameter(index=c, value=value.lower)    # 顧客需要量下限

z = Variable(index=(d,c), lb=0)  # 倉庫から顧客への輸送量

problem = Problem(name='輸送問題')
problem += Sum(cost[d,c]*z[d,c], (d,c))
problem += Sum(z[d,c], c) <= upper[d], '倉庫上限'
problem += Sum(z[d,c], d) >= lower[c], '顧客需要'

print(problem)
problem.solve(silent=True)
assert problem.status == NuoptStatus.OPTIMAL

dc = z[d,c].val>0.1  # 値が入った部分
print(z[dc].val)
Printf('倉庫{} -> 顧客{}: {:.0f}', dc, z[dc].val)

添字を使うことで，定式化と同じように記述できています．
この出力は次のようになります．

Problem(name='輸送問題', type=min):
[constraints]
倉庫上限:
-z[d0,c0]-z[d0,c1]-z[d0,c2]>=-50495
-z[d1,c0]-z[d1,c1]-z[d1,c2]>=-99347
顧客需要:
z[d0,c0]+z[d1,c0]>=6891
z[d0,c1]+z[d1,c1]>=664
z[d0,c2]+z[d1,c2]>=4243

[objective]
Sum((cost[d,c]*z[d,c])[d,c], (d,c)):
66*z[d0,c0]+63*z[d0,c1]+52*z[d0,c2]+39*z[d1,c0]+62*z[d1,c1]+46*z[d1,c2]

z[d1,c0].val=6890.999999521996
z[d1,c1].val=663.9999587327063
z[d1,c2].val=4242.999993921036
倉庫d1 -> 顧客c0: 6891
倉庫d1 -> 顧客c1: 664
倉庫d1 -> 顧客c2: 4243

3. 輸送問題(スパース ver.)

今度は輸送経路のうち，一部の経路のみが定義されるように変更してみましょう．定式化は以下のようになります．

構成要素	名称	説明
添字	$dc \in DC=\{(d0,c1), (d1,c0), .. \}$	倉庫から顧客への輸送経路(一部)
定数	$cost_{d,c}, (d,c) \in DC$	輸送コスト
定数	$upper_d, d \in DC(0)$	倉庫上限
定数	$lower_c, c \in DC(1)$	顧客需要
変数	$z_{d,c}, (d,c) \in DC$	輸送量
目的関数	$\sum_{d,c, (d,c) \in DC} cost_{d,c}\cdot z_{d,c}$	輸送コスト
制約式	$\sum_{c, (d,c) \in DC} z_{d,c} \le upper_d$	倉庫上限
制約式	$\sum_{d, (d,c) \in DC} z_{d,c} \ge lower_c$	顧客需要

ここで $DC(0)$, $DC(1)$ は，集合 $DC$ の 1 次元目，2 次元目を射影した集合です．
この定式化のモデリングは，次のようになります．

### データ用意
from random import seed, randint, random

DSIZE, CSIZE, RATE = 2, 3, 0.5

class value:
    seed(0)
    upper = {f'd{d}': randint(1, 100000) for d in range(DSIZE)}
    lower = {f'c{c}': randint(1, 10000) for c in range(CSIZE)}
    cost = {(f'd{d}',f'c{c}'): randint(1, 100) for d in range(DSIZE) for c in range(CSIZE) if random() < RATE}


### 最適化
from pysimple import Element, Parameter, Variable, Sum, Problem, NuoptStatus, Printf

dc = Element(value=value.cost.keys())  # 倉庫と顧客の疎な添字

cost = Parameter(index=dc, value=value.cost)       # 倉庫から顧客への輸送コスト
upper = Parameter(index=dc(0), value=value.upper)  # 倉庫取扱量上限
lower = Parameter(index=dc(1), value=value.lower)  # 顧客需要量下限

z = Variable(index=dc, lb=0)  # 倉庫から顧客への輸送量

problem = Problem(name='輸送問題')
problem += Sum(cost[dc]*z[dc])
problem += Sum(z[dc], dc(1)) <= upper[dc(0)], '倉庫上限'
problem += Sum(z[dc], dc(0)) >= lower[dc(1)], '顧客需要'

print(problem)
problem.solve(silent=True)
assert problem.status == NuoptStatus.OPTIMAL

dc0 = z[dc].val>0.1  # 値が入った部分
print(z[dc0].val)
Printf('倉庫{} -> 顧客{}: {:.0f}', dc0, z[dc0].val)

この出力は次のようになります．

Problem(name='輸送問題', type=min):
[constraints]
倉庫上限:
-z[d0,c1]>=-50495
-z[d1,c0]-z[d1,c2]>=-99347
顧客需要:
z[d0,c1]>=664
z[d1,c0]>=6891
z[d1,c2]>=4243

[objective]
Sum((cost[dc]*z[dc])[dc], dc):
39*z[d0,c1]+28*z[d1,c0]+97*z[d1,c2]

z[d0,c1].val=664.0000169446201
z[d1,c0].val=6891.000023601436
z[d1,c2].val=4243.0000068127865
倉庫d0 -> 顧客c1: 664
倉庫d1 -> 顧客c0: 6891
倉庫d1 -> 顧客c2: 4243

先ほどのモデリングと比べて，一部の経路しか定義されていないことが分かりますね．

4. おわりに

現在の PySIMPLE は弊社の数理最適化パッケージ Numerical Optimizer に付属するモデリング言語となっています．来年 3 月リリースのバージョンでは MIP(混合整数計画問題) に加え，QP(二次計画問題)，メタヒューリスティクスアルゴリズム も利用できるようになります．

明日は 梅谷先生(@umepon) による「資材切出し問題と列生成法」です．