PuLPとPICOSによる数理最適化

PuLPとは

機能 :  LPのモデリングツール

守備範囲 : LPとIP 整数変数も扱うことができます. ドキュメントはHPから入手可能.このドキュメントからの引用を多く使用します.

ソルバは別途用意する必要があります.今回はCBCを使用します.私が調べた範囲では、GLPKを使うと、LPの双対変数の値を取得することができません。ここでは後ほど双対変数の値を用いる場面がありますが、その部分はGLPKでは実行することができないことに留意ください。

線形計画問題の定義

次の線形計画問題をPuLPで定義します.

(画像表示準備)

In [67]:
from IPython.display import Image
fig = Image(url='http://mathopt.sakura.ne.jp/LP.png')
fig 
Out[67]:

使う準備 : importする

In [68]:
import pulp

問題の定義 :

In [69]:
prob = pulp.LpProblem("myProblem",pulp.LpMinimize)

最大化のときは,LpMaximize.生成した問題に,probと名前を付けて覚えておきます。

変数の定義 :

In [70]:
z1 = pulp.LpVariable("z1",0,3)

三つの引数はそれぞれ、変数の名前、下限値、上限値です。

In [71]:
z2 = pulp.LpVariable("z2",0,1)

制約式の定義 :

In [72]:
prob += z1 + z2 <= 2

+=の右側が(線形)不等式であれば,それが制約式としてprobに追加されます.結果をprintしてみると,

In [73]:
print prob
myProblem:
MINIMIZE
None
SUBJECT TO
_C1: z1 + z2 <= 2

VARIABLES
z1 <= 3 Continuous
z2 <= 1 Continuous

目的関数の定義 :

In [74]:
prob += -1*z1 - 2*z2 

制約式の定義のところで述べましたが、+=の右側が線形不等式ではなく線形の関数であれば、それは目的関数として定義されます。

In [75]:
print prob
myProblem:
MINIMIZE
-1*z1 + -2*z2 + 0
SUBJECT TO
_C1: z1 + z2 <= 2

VARIABLES
z1 <= 3 Continuous
z2 <= 1 Continuous

では、目的関数を二度定義するとどうなるでしょうか?重複して定義すると、最後に定義したものが目的関数として設定されます。

In [76]:
prob+=20*z1+30*z2
In [77]:
print prob
myProblem:
MINIMIZE
20*z1 + 30*z2 + 0
SUBJECT TO
_C1: z1 + z2 <= 2

VARIABLES
z1 <= 3 Continuous
z2 <= 1 Continuous

これで問題が定義できました

In [78]:
prob += -z1-2*z2

解く

こうして定義した問題probを解くには,solve()を実行します.複数のソルバが実行可能な環境であれば,引数を指定することで,使用するソルバを指定できます.

In [79]:
status = prob.solve()

実行結果の状態は、LpStatus[status]で取得します。

In [80]:
print pulp.LpStatus[status]
Optimal

最適値を得るには、pulp.value(prob.objective)を実行します.

In [81]:
print pulp.value(prob.objective)
-3.0

最適解を得るには、次のように実行します。

In [82]:
print pulp.value(z1), pulp.value(z2)
1.0 1.0

双対変数の取得

LPを解いたときには,各制約式に対する双対変数の値を取得することができます.ただし,ソルバとしてGLPKを用いたときは取得できないようです。私の環境が原因か、仕様なのかはいまだ不明です。

双対変数を取得するためには,LPの制約式に名前を付けておくと便利です。先ほどの問題を一から定義し直し、今度は制約式に名前を付けましょう。

In [83]:
prob_d = pulp.LpProblem("myProblem",pulp.LpMinimize)
x1=pulp.LpVariable("x1",0,3)
x2=pulp.LpVariable("x2",0,1)
prob_d+=-x1-2*x2
prob_d+=x1+x2<=2,"First_Constraint"
print prob_d
myProblem:
MINIMIZE
-1*x1 + -2*x2 + 0
SUBJECT TO
First_Constraint: x1 + x2 <= 2

VARIABLES
x1 <= 3 Continuous
x2 <= 1 Continuous

さて、この線形計画問題を解いて、制約式に対する双対変数の値を取得します。まずは解いて、最適解を得ます。

In [84]:
status = prob_d.solve()
print pulp.LpStatus[status]
Optimal

つぎに、制約式の名前を指定して、それに対応する双対変数の値を得ます。

In [85]:
print prob_d.constraints.keys()
['First_Constraint']
In [86]:
print prob_d.constraints["First_Constraint"].pi
-1.0

集合分割問題定式化 : Vehicle Routing Problem

配送計画問題の集合分割問題定式化は,次のように書くことができます. 記法は若干複雑ですが,集合分割問題です.

In [87]:
fig = Image(url='http://mathopt.sakura.ne.jp/VRP.png')
fig 
Out[87]:

Vehicle Routing Problem(配送計画問題)の集合分割問題定式化をPuLPを用いて記述してみます.まずは,最小化問題として問題を定義します.

In [88]:
routing_model=pulp.LpProblem("VRP",pulp.LpMinimize)

集合分割問題 : 変数の定義

リストから変数を生成する

配送計画問題の集合分割問題定式化では、貨物の集合を、その部分集合に分割します。そのために、貨物の部分集合一つにつき、0-1変数を一つ定義します。部分集合の数が膨大になると、それらの一つ一つに名前を付けて変数を定義する必要があります。それを簡単に実現することのできる命令がpulpには用意されています。

いま、貨物の部分集合の名前を要素とするリストpow_sがあるとします。pow_sは、pow_s=['s1','s2','s3']と、三つの要素からなるとします.PuLPには、リストを与えると、その各要素に対して変数を定義することのできる命令があります。それはつぎのものです.

In [89]:
pow_s=['s1','s2','s3']
pow_s_vars=pulp.LpVariable.dicts('x', pow_s, lowBound=0, upBound=1, cat=pulp.LpBinary)
print pow_s_vars
{'s3': x_s3, 's2': x_s2, 's1': x_s1}

この命令により、3つの変数が生成され、それらは辞書pow_s_varsにおさめられました。この辞書のキーは、生成に用いたリストの各要素です。キーとして要素を指定すると、その要素に対して生成された変数が値として得られます。例えば、s1に対して定義された変数を得るには、pow_s_vars['s1']とします.

In [90]:
print pow_s_vars['s1']
x_s1

命令pulp.LpVariable.dicts()の第一引数には、生成する変数につける名前の最初の部分を指定します。今の例ではxを指定したので、生成される変数の名前は、リストの要素の先頭に’x_’を付加した’x_s1'になります。

第二引数には、リストを指定します。第三引数には、変数の下限値を指定します。第四引数には上限値を指定します。第五引数には変数のタイプを指定します。ここでは0-1変数を定義したいので、LpBinaryを 指定します。連続変数を定義したい場合は、LpContinuousを指定します。

なお、"lowBound=","upBound=","cat="は省略することができます。

運搬車ごとに実行可能ルートを定義する

ここでは,性能の違う2台の運搬車'v0'と'v1'が使用可能とします.v0で実行可能なルートを要素とするリストを,v_feas_rt['v0']=[]と定義しておきます.同様に,v1が実行可能なルートを要素とするリストを,v_feas_rt['v1']=[]と定義しておきます.v0の実行可能なルートはs0,s1,s2の3つ,v1の実行可能なルートはs3,s4,s5の3つとします.さらに処理するべき貨物はc0,c1,c2,c3の4つとします.これらの貨物を要素とするリストcargosを定義します.

In [91]:
cargos=["c0","c1","c2","c3"];  print cargos
['c0', 'c1', 'c2', 'c3']

実行可能ルートs0において処理する貨物は,c0とします.このとき,v_feas_rt['v0']の要素として,('v0',('c0'))を追加します.すなわち,最初の要素に運搬車の名前を,2番目の要素としてそのルートで処理する貨物のタプルを指定します.同様に,ルートs1において処理する貨物はc0とc1として,このルートをv_feas_rt['v0']に追加します.さらに,ルートs2において処理する貨物はc1とc3とします.

In [92]:
v_feas_rt={}; v_feas_rt['v0']=[]
v_feas_rt['v0']+=[('v0',('c0'))]
v_feas_rt['v0']+=[('v0',('c0','c1'))]
v_feas_rt['v0']+=[('v0',('c1','c3'))]
print v_feas_rt['v0']
[('v0', 'c0'), ('v0', ('c0', 'c1')), ('v0', ('c1', 'c3'))]

各ルートのコストを値としてもつ辞書route_costもあわせて定義します.

In [93]:
route_cost={}
In [94]:
route_cost[('v0',('c0'))]=100
route_cost[('v0',('c0','c1'))]=200
route_cost[('v0',('c1','c3'))]=300

運搬車v1の実行可能ルートも同様にv_feas_rt['v1']に設定します.

In [95]:
v_feas_rt['v1']=[]
v_feas_rt['v1']+=[('v1',('c3'))]
v_feas_rt['v1']+=[('v1',('c0','c2'))]
v_feas_rt['v1']+=[('v1',('c0','c3'))]
print v_feas_rt['v1']
[('v1', 'c3'), ('v1', ('c0', 'c2')), ('v1', ('c0', 'c3'))]
In [96]:
route_cost[('v1',('c3'))]=150
route_cost[('v1',('c0','c2'))]=220
route_cost[('v1',('c0','c3'))]=180
print route_cost
{('v1', ('c0', 'c3')): 180, ('v1', ('c0', 'c2')): 220, ('v0', ('c0', 'c1')): 200, ('v0', ('c1', 'c3')): 300, ('v0', 'c0'): 100, ('v1', 'c3'): 150}

こうして定義したルートを,集合分割問題の制約行列として表わすと,次のようになります.

x[0,0] x[0,1] x[0,2] x[1,0] x[1,1] x[1,2] y[c0] y[c1] y[c2] y[c3]
c0 1 1 1 1 1
c1 1 1 1
c2 1 1
c3 1 1 1 1

集合分割問題 : 実行可能ルートごとに変数を定義する

先ほど運搬車ごとの実行可能ルートを表すv_feas_rtを定義しましたが,これらの要素を総てあわせたfeas_rtを用意しておきます.

In [97]:
vehicles,feas_rt=['v0','v1'],[]
for v in vehicles:
        feas_rt+=v_feas_rt[v]
print feas_rt
[('v0', 'c0'), ('v0', ('c0', 'c1')), ('v0', ('c1', 'c3')), ('v1', 'c3'), ('v1', ('c0', 'c2')), ('v1', ('c0', 'c3'))]

全ての実行可能ルートは,feas_rtの要素として得られます.そこで,LpVariable.dicts()を用いてfeas_rtの各要素に対して変数を定義します.

In [98]:
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
print x
{('v1', ('c0', 'c3')): x_('v1',_('c0',_'c3')), ('v1', ('c0', 'c2')): x_('v1',_('c0',_'c2')), ('v0', ('c0', 'c1')): x_('v0',_('c0',_'c1')), ('v0', ('c1', 'c3')): x_('v0',_('c1',_'c3')), ('v0', 'c0'): x_('v0',_'c0'), ('v1', 'c3'): x_('v1',_'c3')}

こうして,各実行可能ルートに対して変数xが定義できました.変数としては,xの他にyも定義します.yは,各貨物に対して1つ定義します

.これは,貨物がどのルートでも処理されないときに1となる0-1変数です.この変数があることにより,集合分割問題は常に実行可能になります.

In [99]:
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
print y
{'c3': y_c3, 'c2': y_c2, 'c1': y_c1, 'c0': y_c0}

ここで,xとyは,0以上1以下の連続変数として定義しました.本来は0-1変数ですが,列生成ではLP緩和問題として解いて双対変数の値を取り出して使うので,最初から連続変数として緩和してあります.

集合分割問題:制約式の定義

まず,貨物を要素とするタプル cargos=('c0','c1','c2','c3')を定義します

.制約式は各貨物に対して1つずつ定義します.そのためにこのタプルをつかいます.例として,貨物c0に対する制約式を取り上げます.c0に対する制約式は,c0を含む実行可能ルートに対する変数を足したものに,y[c0]を足したもの(線形式)が,1に等しい,というものとして定義します.ここで,c0を含む実行可能ルートに対する変数は,次の命令で得られます.

In [100]:
[x[route] for route in feas_rt if 'c0' in route[1] ]
Out[100]:
[x_('v0',_'c0'),
 x_('v0',_('c0',_'c1')),
 x_('v1',_('c0',_'c2')),
 x_('v1',_('c0',_'c3'))]

これらの変数の和を得るには,このリストをsum()の引数に与えます.結局,問題routing_modelに,各貨物に対する制約式を追加するにはつぎの命令を実行します.

In [101]:
for c in cargos:
    routing_model+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
In [102]:
print routing_model
VRP:
MINIMIZE
None
SUBJECT TO
Must_assign_c0: x_('v0',_'c0') + x_('v0',_('c0',_'c1'))
 + x_('v1',_('c0',_'c2')) + x_('v1',_('c0',_'c3')) + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1')) + x_('v0',_('c1',_'c3')) + y_c1 = 1

Must_assign_c2: x_('v1',_('c0',_'c2')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c1',_'c3')) + x_('v1',_'c3')
 + x_('v1',_('c0',_'c3')) + y_c3 = 1

VARIABLES
x_('v0',_'c0') <= 1 Continuous
x_('v0',_('c0',_'c1')) <= 1 Continuous
x_('v0',_('c1',_'c3')) <= 1 Continuous
x_('v1',_'c3') <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c0',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

これで,各貨物を処理するという制約式を定義できました.制約式としてもう一種類のものを定義します.それは,各運搬車が実行するルートの数は1以下,という制約です.運搬車vが処理するルートに対応する変数は[x[route] for route in v_feas_rt[v]]として得られるので,各船に対してこの制約式を定義するにはつぎの命令を実行します.

In [103]:
for v in vehicles:
    routing_model += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
print routing_model
VRP:
MINIMIZE
None
SUBJECT TO
Must_assign_c0: x_('v0',_'c0') + x_('v0',_('c0',_'c1'))
 + x_('v1',_('c0',_'c2')) + x_('v1',_('c0',_'c3')) + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1')) + x_('v0',_('c1',_'c3')) + y_c1 = 1

Must_assign_c2: x_('v1',_('c0',_'c2')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c1',_'c3')) + x_('v1',_'c3')
 + x_('v1',_('c0',_'c3')) + y_c3 = 1

Vehicle_assign_v0: x_('v0',_'c0') + x_('v0',_('c0',_'c1'))
 + x_('v0',_('c1',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_'c3') + x_('v1',_('c0',_'c2'))
 + x_('v1',_('c0',_'c3')) <= 1

VARIABLES
x_('v0',_'c0') <= 1 Continuous
x_('v0',_('c0',_'c1')) <= 1 Continuous
x_('v0',_('c1',_'c3')) <= 1 Continuous
x_('v1',_'c3') <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c0',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

集合分割問題:目的関数の定義

次に,目的関数を定義します.目的関数は,採用するルートのコストの和と,変数yにペナルティを示す大きな係数をかけたものの和として定義します.

In [104]:
routing_model+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print routing_model
VRP:
MINIMIZE
100*x_('v0',_'c0') + 200*x_('v0',_('c0',_'c1')) + 300*x_('v0',_('c1',_'c3')) + 150*x_('v1',_'c3') + 220*x_('v1',_('c0',_'c2')) + 180*x_('v1',_('c0',_'c3')) + 1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0
SUBJECT TO
Must_assign_c0: x_('v0',_'c0') + x_('v0',_('c0',_'c1'))
 + x_('v1',_('c0',_'c2')) + x_('v1',_('c0',_'c3')) + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1')) + x_('v0',_('c1',_'c3')) + y_c1 = 1

Must_assign_c2: x_('v1',_('c0',_'c2')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c1',_'c3')) + x_('v1',_'c3')
 + x_('v1',_('c0',_'c3')) + y_c3 = 1

Vehicle_assign_v0: x_('v0',_'c0') + x_('v0',_('c0',_'c1'))
 + x_('v0',_('c1',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_'c3') + x_('v1',_('c0',_'c2'))
 + x_('v1',_('c0',_'c3')) <= 1

VARIABLES
x_('v0',_'c0') <= 1 Continuous
x_('v0',_('c0',_'c1')) <= 1 Continuous
x_('v0',_('c1',_'c3')) <= 1 Continuous
x_('v1',_'c3') <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c0',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

集合分割問題:双対変数を取り出す

これで,集合分割問題の線形緩和問題が定義できました.さて,列生成法では,この線形緩和問題を解いて得られる双対変数の値を用いて,効率的な実行可能ルートを新たに生成・追加していきます.双対変数は制約式毎に得られます.貨物の割当制約式に対する双対変数の値を辞書dual_cargoに,運搬車の使用回数の制約式に対する双対変数の値を辞書dual_shipに保存することにします.

In [105]:
dual_cargo,dual_vehicle={},{}

まず,定義した線形緩和問題を解きます.

In [106]:
status = routing_model.solve()
print pulp.LpStatus[status]
Optimal

さて,routing_modelを定義するとき,制約式に名前を付けました.この名前を用いて,制約式に対する双対変数の値を取得します.貨物の割当制約式に対する双対変数の値は

In [107]:
for c in cargos:
    dual_cargo[c]=routing_model.constraints["Must_assign_%s"%c].pi

運搬車の使用回数の制約式に対する双対変数の値は

In [108]:
for v in vehicles:
    dual_vehicle[v]=routing_model.constraints["Vehicle_assign_%s"%v].pi
In [109]:
print dual_cargo,dual_vehicle
{'c3': 150.0, 'c2': 1000.0, 'c1': 150.0, 'c0': 0.0} {'v0': 0.0, 'v1': 0.0}

集合分割問題:列を生成する

networkxで部分問題を定義する

上の例では,各運搬車の実行可能ルートが既知であるとしていました.しかし,実際には実行可能ルートはすべてはわからないときもあります.このようなときは,実行可能ルートがs-tパスとして表わされるようなネットワークを定義して,その上のパスとして求めるのが便利です.ここでは,ノード集合{s,c0,c1,c2,c3,t},アーク集合{(s,c0),(s,c1),(c0,c1),(c0,c2),(c0,c3),(c1,c3),(c2,t),(c3,c2),(c3,t)}から定義されるネットワークを用います.各アークには適当にコストを定義しておきます。

In [110]:
import networkx as nx
G={}
G['v0']=nx.DiGraph()
G['v0'].add_weighted_edges_from([('s','c0',100),('s','c1',50),('c0','c1',20),('c0','c2',40),('c0','c3',30),('c1','c3',35),('c2','t',70),('c3','c2',15),('c3','t',45)])
G['v1']=nx.DiGraph()
G['v1'].add_weighted_edges_from([('s','c0',150),('s','c1',10),('c0','c2',40),('c0','c3',30),('c1','c3',35),('c2','t',70),('c3','t',45)])
In [111]:
import networkx as nx
pos=nx.shell_layout(G['v0']);
nx.draw_shell(G['v0'],node_size=1000,node_color='w');
nx.draw_networkx_edge_labels(G['v0'],pos);
In [112]:
import networkx as nx
pos=nx.shell_layout(G['v1']);
nx.draw_shell(G['v1'],node_size=1000,node_color='w');
nx.draw_networkx_edge_labels(G['v1'],pos);

このくらいの規模のネットワークだと,すべてのs-tパスを列挙してもたいした数にはなりません.しかし,実務規模の問題例だと,s-tパスの数は数十万以上になることは普通にあります.したがって,すべてのパスを列挙することは避けたほうがよいです.そこで,実行可能ルートの部分集合で定義した集合分割問題の線形緩和問題の双対変数の値を用いて,反復的にルート(s-tパス)を生成する方法が有効です.

列が空の集合分割問題(線形緩和問題) を定義する

まず,実行可能ルートが空の状態で,集合分割問題を定義します.

In [113]:
import pulp
cargos=["c0","c1","c2","c3"]
vehicles,feas_rt=['v0','v1'],[]
v_feas_rt={}
v_feas_rt['v0'],v_feas_rt['v1']=[],[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]

master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
#for v in vehicles:
#    master_p += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p
MasterProblem:
MINIMIZE
1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0.0
SUBJECT TO
Must_assign_c0: y_c0 = 1

Must_assign_c1: y_c1 = 1

Must_assign_c2: y_c2 = 1

Must_assign_c3: y_c3 = 1

VARIABLES
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

In [114]:
status=master_p.solve()
print pulp.LpStatus[status]
Optimal

双対変数の値を取り出し,部分問題のコストに反映する

次に,各貨物の割り当て制約式に対する双対変数の値を取り出します.

In [115]:
dual_cargo={}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
print dual_cargo
{'c3': 1000.0, 'c2': 1000.0, 'c1': 1000.0, 'c0': 1000.0}

これらの値を使って,ネットワークの各アークのコストを定義し直します.具体的には,アーク(i,j)の始点iが貨物を表わす(=始点sでない)ときには,コストを,[アークの(i,j)の重み]-dual_cargo[i]に更新します.

実際の配送計画問題では,こうしてアークのコストを定義し直したネットワーク上で(負のコストを持つ)最短路問題を解くことによって,集合分割問題に追加するべき実行ルートを生成します.ここでは簡単のため,最短路問題を解くことはせず,全列挙したs-tバスの中から最も費用がちいさくなるものを選ぶことにします.

先ほど定義したネットワーク上でのs-tパスを列挙するには,networkxのall_simple_paths()という命令を使うことができます.

さて,all_simple_paths(G['v0'])で得られる各s-tパスに対して,含まれるアークのコストの和を求めます.

In [116]:
print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
	
s-t paths for v0
['s', 'c1', 'c3', 'c2', 't'] -2830.0
['s', 'c1', 'c3', 't'] -1870.0
['s', 'c0', 'c3', 'c2', 't'] -2785.0
['s', 'c0', 'c3', 't'] -1825.0
['s', 'c0', 'c2', 't'] -1790.0
['s', 'c0', 'c1', 'c3', 'c2', 't'] -3760.0
['s', 'c0', 'c1', 'c3', 't'] -2800.0

これらから,最もコストのちいさいs-tパスは,['s', 'c0', 'c1', 'c3', 'c2', 't'] であることが分かります.同様に,運搬車v1に対するネットワーク上のs-tパスのコストを求めます.

In [117]:
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
	
['s', 'c1', 'c3', 't'] -1870.0
['s', 'c0', 'c3', 't'] -1825.0
['s', 'c0', 'c2', 't'] -1790.0

このうち最もコストの小さいs-tパスは,['s', 'c1', 'c3', 't'] であることが分かります.

s-tパスを集合分割問題の列に追加する

ということで,v0に対する実行可能ルートとして,['s', 'c0', 'c1', 'c3', 'c2', 't']を,v1に対する実行可能ルートとして,['s', 'c1', 'c3', 't']を追加することにします.具体的には,v_feas_rt['v0'],v_feas_rt['v1']に,そのルートで処理される貨物を要素とする要素を加えます.同時に,route_costにもコストを表わす要素を追加しますが,このコストとしては,双対変数の値を取り除いた正味のアークのコストの和を指定します.

In [118]:
route_cost={}
v_feas_rt['v0']+=[('v0',('c0','c1','c2','c3'))]
v_feas_rt['v1']+=[('v1',('c1','c3'))]
route_cost[('v0',('c0','c1','c2','c3'))]=240
route_cost[('v1',('c1','c3'))]=90

これで,一度の反復が終了です.次に,二番の反復を開始します.

集合分割問題の線形緩和を再度解く

In [119]:
import pulp
cargos=["c0","c1","c2","c3"]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p
MasterProblem:
MINIMIZE
240*x_('v0',_('c0',_'c1',_'c2',_'c3')) + 90*x_('v1',_('c1',_'c3')) + 1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0
SUBJECT TO
Must_assign_c0: x_('v0',_('c0',_'c1',_'c2',_'c3')) + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3'))
 + y_c1 = 1

Must_assign_c2: x_('v0',_('c0',_'c1',_'c2',_'c3')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3'))
 + y_c3 = 1

Vehicle_assign_v0: x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_('c1',_'c3')) <= 1

VARIABLES
x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1 Continuous
x_('v1',_('c1',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

これで新しい集合分割問題(の線形緩和問題)が定義されました.これを解きます.今度は,船舶の使用回数に関する制約式が空ではないので,これらの制約式に対する双対変数も取り出します.

In [120]:
status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle
Optimal
{'c3': 0.0, 'c2': 1000.0, 'c1': 90.0, 'c0': 1000.0} {'v0': 0.0, 'v1': 0.0}

最も効率の良いs-tパスをみつける

これらの値をアークのコストに反映した状態で,最もコストの小さいs-tパスを求めます.

In [121]:
print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	
s-t paths for v0
['s', 'c1', 'c3', 'c2', 't'] -920.0
['s', 'c1', 'c3', 't'] 40.0
['s', 'c0', 'c3', 'c2', 't'] -1785.0
['s', 'c0', 'c3', 't'] -825.0
['s', 'c0', 'c2', 't'] -1790.0
['s', 'c0', 'c1', 'c3', 'c2', 't'] -1850.0
['s', 'c0', 'c1', 'c3', 't'] -890.0
s-t paths for v1
['s', 'c1', 'c3', 't'] 0.0
['s', 'c0', 'c3', 't'] -775.0
['s', 'c0', 'c2', 't'] -1740.0

v0に対するネットワーク上で最小のコストをもつs-tパスは,['s', 'c0', 'c1', 'c3', 'c2', 't'],v1に対するネットワーク上では['s', 'c0', 'c2', 't']であることがわかりました.v0に対する最小コストのs-tパスは先ほど追加したのと同じですので,今回は,v1に対するルート['s', 'c0', 'c2', 't']のみ追加します.

In [122]:
v_feas_rt['v1']+=[('v1',('c0','c2'))]
route_cost[('v1',('c0','c2'))]=260
In [123]:
import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
print feas_rt
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p
[('v0', ('c0', 'c1', 'c2', 'c3')), ('v1', ('c1', 'c3')), ('v1', ('c0', 'c2'))]
MasterProblem:
MINIMIZE
240*x_('v0',_('c0',_'c1',_'c2',_'c3')) + 260*x_('v1',_('c0',_'c2')) + 90*x_('v1',_('c1',_'c3')) + 1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0
SUBJECT TO
Must_assign_c0: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c0',_'c2'))
 + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3'))
 + y_c1 = 1

Must_assign_c2: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c0',_'c2'))
 + y_c2 = 1

Must_assign_c3: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3'))
 + y_c3 = 1

Vehicle_assign_v0: x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_('c0',_'c2')) + x_('v1',_('c1',_'c3')) <= 1

VARIABLES
x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c1',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

In [124]:
status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle
Optimal
{'c3': 0.0, 'c2': 240.0, 'c1': 0.0, 'c0': 0.0} {'v0': 0.0, 'v1': 0.0}
In [125]:
print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	
s-t paths for v0
['s', 'c1', 'c3', 'c2', 't'] -70.0
['s', 'c1', 'c3', 't'] 130.0
['s', 'c0', 'c3', 'c2', 't'] -25.0
['s', 'c0', 'c3', 't'] 175.0
['s', 'c0', 'c2', 't'] -30.0
['s', 'c0', 'c1', 'c3', 'c2', 't'] 0.0
['s', 'c0', 'c1', 'c3', 't'] 200.0
s-t paths for v1
['s', 'c1', 'c3', 't'] 90.0
['s', 'c0', 'c3', 't'] 225.0
['s', 'c0', 'c2', 't'] 20.0

v0に対する最小コストのs-t パスは,['s', 'c1', 'c3', 'c2', 't'] ,v1に対するs-tパスは['s', 'c0', 'c2', 't'] ですが,v1に対する パスはコストが正なので,追加する必要はありません.v0に対するs-tパスのみ追加することにします。

In [126]:
v_feas_rt['v0']+=[('v0',('c1','c2','c3'))]
route_cost[('v0',('c1','c2','c3'))]=170
In [127]:
import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p
MasterProblem:
MINIMIZE
240*x_('v0',_('c0',_'c1',_'c2',_'c3')) + 170*x_('v0',_('c1',_'c2',_'c3')) + 260*x_('v1',_('c0',_'c2')) + 90*x_('v1',_('c1',_'c3')) + 1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0
SUBJECT TO
Must_assign_c0: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v1',_('c0',_'c2'))
 + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3')) + y_c1 = 1

Must_assign_c2: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c1',_'c2',_'c3')) + x_('v1',_('c0',_'c2')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3')) + y_c3 = 1

Vehicle_assign_v0: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c1',_'c2',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_('c0',_'c2')) + x_('v1',_('c1',_'c3')) <= 1

VARIABLES
x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1 Continuous
x_('v0',_('c1',_'c2',_'c3')) <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c1',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

この問題を解くと,つぎのようになります.

In [128]:
status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle
Optimal
{'c3': 0.0, 'c2': 170.0, 'c1': 0.0, 'c0': 70.0} {'v0': 0.0, 'v1': 0.0}
In [129]:
print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	
s-t paths for v0
['s', 'c1', 'c3', 'c2', 't'] 0.0
['s', 'c1', 'c3', 't'] 130.0
['s', 'c0', 'c3', 'c2', 't'] -25.0
['s', 'c0', 'c3', 't'] 105.0
['s', 'c0', 'c2', 't'] -30.0
['s', 'c0', 'c1', 'c3', 'c2', 't'] 0.0
['s', 'c0', 'c1', 'c3', 't'] 130.0
s-t paths for v1
['s', 'c1', 'c3', 't'] 90.0
['s', 'c0', 'c3', 't'] 155.0
['s', 'c0', 'c2', 't'] 20.0

v0のs-tパスで,コストが負になるものが2つありますので,両方入れてしまいましょう。 

In [130]:
v_feas_rt['v0']+=[('v0',('c0','c2','c3'))]
v_feas_rt['v0']+=[('v0',('c0','c2'))]
route_cost[('v0',('c0','c2','c3'))]=215
route_cost[('v0',('c0','c2'))]=210
In [131]:
import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p
MasterProblem:
MINIMIZE
240*x_('v0',_('c0',_'c1',_'c2',_'c3')) + 210*x_('v0',_('c0',_'c2')) + 215*x_('v0',_('c0',_'c2',_'c3')) + 170*x_('v0',_('c1',_'c2',_'c3')) + 260*x_('v1',_('c0',_'c2')) + 90*x_('v1',_('c1',_'c3')) + 1000*y_c0 + 1000*y_c1 + 1000*y_c2 + 1000*y_c3 + 0
SUBJECT TO
Must_assign_c0: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v0',_('c0',_'c2'))
 + x_('v0',_('c0',_'c2',_'c3')) + x_('v1',_('c0',_'c2')) + y_c0 = 1

Must_assign_c1: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c1',_'c2',_'c3')) + x_('v1',_('c1',_'c3')) + y_c1 = 1

Must_assign_c2: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v0',_('c0',_'c2'))
 + x_('v0',_('c0',_'c2',_'c3')) + x_('v0',_('c1',_'c2',_'c3'))
 + x_('v1',_('c0',_'c2')) + y_c2 = 1

Must_assign_c3: x_('v0',_('c0',_'c1',_'c2',_'c3'))
 + x_('v0',_('c0',_'c2',_'c3')) + x_('v0',_('c1',_'c2',_'c3'))
 + x_('v1',_('c1',_'c3')) + y_c3 = 1

Vehicle_assign_v0: x_('v0',_('c0',_'c1',_'c2',_'c3')) + x_('v0',_('c0',_'c2'))
 + x_('v0',_('c0',_'c2',_'c3')) + x_('v0',_('c1',_'c2',_'c3')) <= 1

Vehicle_assign_v1: x_('v1',_('c0',_'c2')) + x_('v1',_('c1',_'c3')) <= 1

VARIABLES
x_('v0',_('c0',_'c1',_'c2',_'c3')) <= 1 Continuous
x_('v0',_('c0',_'c2')) <= 1 Continuous
x_('v0',_('c0',_'c2',_'c3')) <= 1 Continuous
x_('v0',_('c1',_'c2',_'c3')) <= 1 Continuous
x_('v1',_('c0',_'c2')) <= 1 Continuous
x_('v1',_('c1',_'c3')) <= 1 Continuous
y_c0 <= 1 Continuous
y_c1 <= 1 Continuous
y_c2 <= 1 Continuous
y_c3 <= 1 Continuous

In [132]:
status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle
Optimal
{'c3': 0.0, 'c2': 140.0, 'c1': 30.0, 'c0': 70.0} {'v0': 0.0, 'v1': 0.0}
In [133]:
print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	
s-t paths for v0
['s', 'c1', 'c3', 'c2', 't'] 0.0
['s', 'c1', 'c3', 't'] 100.0
['s', 'c0', 'c3', 'c2', 't'] 5.0
['s', 'c0', 'c3', 't'] 105.0
['s', 'c0', 'c2', 't'] 0.0
['s', 'c0', 'c1', 'c3', 'c2', 't'] 0.0
['s', 'c0', 'c1', 'c3', 't'] 100.0
s-t paths for v1
['s', 'c1', 'c3', 't'] 60.0
['s', 'c0', 'c3', 't'] 155.0
['s', 'c0', 'c2', 't'] 50.0

これで,すべてのコストが非負になりました.これは,もう追加すべきs-t パスがないことを示しています.この状態が,線形緩和問題が最適に解けていることを表わしています.

集合分割問題 : 0-1変数に戻して解く

さて,ここから,もとの集合分割問題の近似解を得ます.具体的には,0-1の連続変数として定義していたものを0-1変数として集合分割問題を定義します.そして,この集合分割 問題の解を,もとのVRPの近似解とします.

In [134]:
import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
spp=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpBinary)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpBinary)
for c in cargos:
    spp+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    spp+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
spp+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
status=spp.solve()
print pulp.LpStatus[status]

for v in spp.variables():
    if v.varValue>0.1:
        print v.name,v.varValue
Optimal
x_('v0',_('c0',_'c1',_'c2',_'c3')) 1.0