Pythonで整数計画問題・線形計画問題を解く(PuLP編)

実装方法その2

PuLPのDocumentを全面的に参考にして,集合分割問題(vehicle routing problemのための)の異なる実装方法を試してみた.こららの方が使いやすいことがあるかもしれない.定式化はその1とおなじなので,その1のところの記述を参考にしてください.
(以前書いた実装方法その1へ飛ぶ)

これから出てくるコードは,表示のために適宜改行しているので,コピーペーストすると動かないので,適宜改行を削除してください.

まず,集合分割問題の定義は次のとおり


1:  def DefMasterProblem(vehicles, vehicle_feasible_routes,cargos,route_cost, vartype):
2:     feasible_routes={}
3:     for v in vehicles:
4:        feasible_routes+=vehicle_feasible_routes[v]
5:     penalty=100000
6:     x=pulp.LpVariable.dicts('R',feasible_routes,lowBound=0,upBound=1, cat=vartype)
7:     y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=vartype)
8:     routing_model = pulp.LpProblem("Routing model",pulp.LpMinimize)
9:     routing_model += sum([route_cost[route]*x[route] for route in feasible_routes])+
                               sum([penalty*y[c] for c in cargos])
10:    for c in cargos:
11:        routing_model += sum([x[route] for route in feasible_routes if 
                               c in route[1]])+y[c]==1,"Must_assign_%s"%c
12:    for v in vehicles:
13:        routing_model += sum([x[route] for route in 
                                vehicle_feasible_routes[v]])<=1,"Vehicle_assign_%s",v
14:    return routing_model,x,y
5つの引数をとります.まず,vehiclesは,運搬車の名前を要素とするtupleです.

vehicles=('v1','v2')
などとします.つぎに,vehicle_feasible_routesは,各運搬車が実行可能なルートを表す辞書です.キーは,vehiclesの要素です.まず空で初期化するとすると,次のようになります.

vehicle_feasible_routes={}
vehicle_feasible_routes['v1']=[]
vehicle_feasible_routes['v2']=[]
そして,v1の実行可能ルート'A'->'B'が見つかったとすると,

vehicle_feasible_routes['v1']+=[('v1',('A','B'))]
としてvehicle_feasible_routes['v1']に追加します.cargosには,処理しなければならない貨物の集合をリストとして与えます.例えば

cargos='A B C'.split()
route_costは,vehicle_feasible_routes[v]内の要素,すなわちルートのコストを,ルートをキーとした辞書として与えます.たとえば,vehicle_feasible_routes['v1']の要素('v1',('A','B'))のコストを指定するには

route_cost[('v1',('A','B'))]=1000
などと指定します(1000はルートのコスト).

まず,各運搬車の実行可能ルートをすべて含んだリストpossible_routesを作ります.このpossible_routesを使ってルートrに対して変数x[r]を生成します.また,cargosを使って各貨物cに対して変数y[c]を生成します


6:  x=pulp.LpVariable.dicts('R',possible_routes,lowBound=0,upBound=1,cat=vartype)
7:  y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=vartype)
次に,LPを生成します.

8:  routing_model = pulp.LpProblem('Routing model',pulp.LpMinimize)
このLPに目的関数を設定します.

9:  routing_model += sum([route_cost[route]*x[route] for route in feasible_routes]) +
                      sum([penalty*y[c] for c in cargos])
ここで,右辺は線形式であることに注意してください.線形であれば目的関数として設定されます.次に,制約式を設定します.

10:  for c in cargos:
11:       routing_model += sum([x[route] for route in feasible routes if c 
                               in route[1]])+y[c]==1,"Must_assign_%s"%c
12:  for v in vehicles: 
13:       routing_model += sum([x[route] for route in vehicle_feasible_routes[v]])<=1,"Vehicle_assign_%s"%v
10,11行目は,各貨物はいずれかのルートで処理されなければならない,という条件です.貨物cが処理できないときは,y[c]=1とするとこの制約式は満たされます.rの制約式の名前は"Must_assign_A"などとなり,あとで双対変数の値を取り出すときに使います.12,13行目は,各運搬車は1つのルートで運航するか,運航しない(割り当てられるルートはない)という条件です.この制約式の名前は"Vehicle_assign_v1"などとなります.

こうして生成したLP(routing_model)をreturnで返します.

DefMasterProblem(...)で生成したLPを解くのが,次のSolveMasterProblem(model,x,y,vehicles_feasible_routes,vehicles,cargos)です.


1:  def SolveMasterProblem(model,x,y,vehicle_feasible_routes,vehicles,cargos):
2:      model.solve(COIN(msg=0))
3:      for v in vehicles:
4:           for route in vehicle_feasible_routes[v]:
5:               if x[route].value()>0.0001:
6:                   print route,x[route].value(),
7:      print ""
8:     for c in cargos:
9:           if y[c].value()>0.0001:
10:               print c
11:     print ""
12:     dual_cargo,dual_vehicle={},{}
13:     for c in cargos:
14:         dual_cargo[c]=model.constraints["Must_assign_%s"%c].pi
15:     for v in vehicles:
16:         if "Vehicle_assign_%s"%v in model.constraints.keys():
17:             dual_vehicle[v]=model.constraints["Vehicle_assign_%s"%v].pi
18:     return dual_cargo,dual_vehicle
まず,2行目でLPを解きます.3行目から11行目は,得られた解を表示するためのものです.12行目では,双対変数の値を記憶する辞書dual_cargo,dual_vehicleを初期化します.13,14行目では,各荷物cに関する制約式"Must_assign_c"に対する双対変数の値を,dual_cargo[c]に代入します.双対変数を取り出すには,model.constraints["Must_assign_c"].piとします.ここで注意するのは,model.solve(COIN(msg=0))のかわりにmodel.solve()としてGLPKを用いたときは,model.constraints["Must_assign_c"].piは正しい値を返しません(仕様のようです).15,16行目では,各運搬車vに対する 制約式"Vehicle_assign_v"に対する双対変数の値を,dual_vehicle[v]に代入します. 18行目では,双対変数を格納した辞書を返します.


実装方法その1
PuLPをインストールする. PuLP HP
PuLPを使うと,整数計画問題・線形計画問題がPython上で解ける.

from pulp import *
とすると,空を飛べるようになる
集合分割問題 \[ \mbox{min} \sum_{v \in V} \sum_{r \in R_v} C_{rv} x_{rv} + \sum_{k \in K}F_{k}y_k \] \[ \mbox{s.t.} \sum_{(r,v):r \in R_v \mbox{ and } k \in R_v} x_{rv} + y_k =1, \forall k \in K \] \[ \sum_{r \in R_v} x_{rv} \leq 1, \forall v \in V \] \[ x_{rv} \in \{0,1\}, y_k \in \{0,1\} \]

これをPuLPを使って定式化する.

集合Kは,処理したい貨物の集合を表す.Vは運搬車の集合を表す.R_vは,運搬車vが実行可能なルートの集合を表す. x_rvは,運搬車vがルートrを採用するとき1,それ以外のとき0をとる変数である.そのルートのコストをC_rvとする.だから目的関数の第一項は,採用するルートのコストの和を表している.y_kは,貨物kをどの採用ルートでも処理しないときに1,それ以外のとき0をとる変数である.だから目的関数の第二項は,積み残しの荷物に対する罰金を表す.

最初の制約式は,「それぞれの荷物kは,ある運搬車によって処理されるか,それともどの運搬車にも処理されないかのいずれか」であることを表している.第一項は,「貨物kを処理するルートに対する変数の和」を表している.したがって,採用されるルートのいずれかで貨物kが処理されるなら,この第一項が1になる.

二番目の制約式は,各運搬車が採用するルートは高々1つであることを表している.

この集合分割問題をPythonとPuLPで書くと,下のようになる.
コード中,tOrdersは上の数式のK,tVehiclesは上の数式のV, tRoutesは上の数式のR_vをvに関してunionをとったものに対応している.

まず,Routesをクラスとして定義している(3〜8行目).これは,ルートのコスト(self.cost, 数式のC_rvに対応),処理する貨物のリスト(self.orders),ルートを処理する運搬車(self.vehicle)を持つ.

貨物の集合,実行可能ルートの集合、運搬車の集合はいずれも辞書で表す(9行目).貨物の集合tOrdersに,5つの要素[1,2,3,4,5]を入れる.実行可能ルートの集合tRoutesには,4つのルートを入れる.各ルートは,クラスRoutesのオブジェクトとする(13〜16行目).運搬車の集合tVehiclesに,2つの要素を入れる(17,18行目).

最小化問題として,問題を定義する(20行目)
変数の種類を,relaxがTrueかFalseかで示すことにする.relaxがTrueなら,変数を連続変数にする(23行目).relaxがFalseなら,変数を0-1変数にする(25行目).

x,yに対応する変数を定義する(27,28行目).1番目の引数は,変数の名前(の先頭部分)を指定する.2番目の引数には辞書を与える.この辞書の各要素に対して、一つずつ変数が定義される.辞書のキーを添字に与えると,そのキーに対応する変数が得られる.たとえば,辞書tOrdersにoというキーがあったとする.このキーに対して定義された変数は,yVars[o]によって得られる.5番目の引数で変数のタイプを指定する.いま,vartypeをLpBinaryとしているので.この変数は0-1変数として定義される.

目的関数を定義するには,

prob += LpSum(...)
とする(30行目).右辺が線形式なら,目的関数として扱われる.右辺が不等式・等式なら,制約式として扱われる.

次に,最初の制約式を定義している(31,32行目).貨物一つに対して制約式を一つ定義している(31行目のfor文).ここでは,

prob+=lpSum(...)+yVars[k]==1, "k"+str(k)
と,右辺が等式になっている.したがって,(目的関数ではなく)制約式として扱われる.....==1,のあとの"k"+str(k)で,この制約式の名前を指定している. lpSum()の中は,xVars[r]の和をとっているのだが,和をとる際に条件を付けている.それは,実行可能ルートrが処理する貨物の中に,いま考えている貨物kが入っているときだけ,和の中に入れる,という条件である.実行可能ルートrが処理する貨物はリストtRoutes[r].ordersに入っている.したがって,"if k in tRoutes[r].orders"を入れればよい.

第二の制約式も,最初の制約式と同様に定義できる(34行目).運搬車1つに対して制約式を一つ定義している.ここでは,xVars[r]の和をとるときに"if tRoutes[r].vehicle==v"という条件を付けている.これは,rは運搬車vの実行可能ルートであるときだけ和の中に入れる,ということだから,上の数式における?sum_{r ?in R_v}を表していることがわかる.

ここまでで,集合分割問題が定義できた.あとは解けばよい.解くには,

prob.solve()

とする(38行目).ソルバとしてCBCを使いたいときは,
prob.solve(COIN(msg=1))
と明示すればよいようだ("msg=1"はメッセージの出力レベル).
得られた最適値を得るには,value(prob.objective)とする(40行目).また,各変数の値を得るには,41〜43行目のようにすればよい.

双対変数の値を得るには,prob.constraints[AAA].piのAAAの部分に,制約式の名前を指定する.制約式の名前は,制約式を生成するとき32,34行目で指定したものである.ソルバとしてGLPKを使った場合は,双対変数変数の値を得るこの命令を実行しようとすると,エラーがでる.prob.solve(COIN(msg=1))とソルバを明示すると,この命令はエラーなく実行できるようである.

  1 from pulp import *
  2 
  3 class Routes:
  4         def __init__(self,_name,_orders=None,_vehicle=None,_cost=0):
  5                 self.name=_name;
  6                 self.orders=_orders
  7                 self.vehicle=_vehicle
  8                 self.cost=_cost;
  9 tOrders,tRoutes,tVehicles={},{},{}
 10 for i in range(5):
 11         tOrders[i]=i
 12 print "tOrders:",tOrders
 13 tRoutes[0]=Routes("r0",[0,1,2],0,100)
 14 tRoutes[1]=Routes("r1",[1,2,4],0,200)
 15 tRoutes[2]=Routes("r2",[1,3,4],1,300)
 16 tRoutes[3]=Routes("r3",[3,4],1,100)
 17 tVehicles[0]=0
 18 tVehicles[1]=1
 19 
 20 prob = LpProblem("Set Partitioning Problem",LpMinimize)
 21 relax=False
 22 if relax:
 23         vartype = LpContinous
 24 else:
 25         vartype = LpBinary
 26 
 27 yVars = LpVariable.dicts("y",tOrders,0,None,vartype)
 28 xVars = LpVariable.dicts("x",tRoutes,0,None,vartype)
 29 Fcost=10000
 30 prob += lpSum([xVars[r]*tRoutes[r].cost for r in tRoutes] + [yVars[k]*Fcost for k in tOrders])
 31 for k in tOrders:
 32         prob += lpSum([xVars[r] for 
               r in tRoutes if k in tRoutes[r].orders]) + yVars[k]==1,"k"+str(k)
 33 for v in tVehicles:
 34         prob += lpSum([-xVars[r] for 
               r in tRoutes if tRoutes[r].vehicle==v])>=-1,"v"+str(v)
 35     
 36 
 37 print prob
 38 #prob.solve()
 39 prob.solve(COIN(msg=1))
 40 print "cost=",value(prob.objective)
 41 for v in prob.variables():
 42         if v.varValue>=0.01:
 43                 print v.name,"=",v.varValue
 44 
 45 for k in tOrders:
 46         print prob.constraints["k"+str(k)].pi
 47 for v in tVehicles:
 48         print prob.constraints["v"+str(v)].pi

LpVariableの使い方

LpVariableは4つのパラメータを持つ.一番目は名前,二番目は下限,三番目は上限,四番目はデータの種類(連続か離散か).下限と上限を指定しなければ,下限と上限はないという設定になる.Noneと明示しても,下限と上限はないという設定になる. データの種類はLpContinuous,LpIntegerかLpBinaryで,設定しなければLpContinuousになる.


x=LpVariable("x",0,1,LpContinous)
指定する引数を明示してもよい.たとえば次のようにしてもよい.

x=LpVariable("x",lowBound=0,upBound=1,cat=LpContinuous)

変数LpVariablesは,辞書として定義することもできる. (A Blending Problemの例). 例えば,Ingredientsというリストが


Ingredients=['CHICKEN','BEEF','MUTTON']
と定義されていたら,Ingredientsの各要素を参照キーとする変数は,次の命令で生成できる.

ingredient_vars = LpVariable.dicts("Ingr",Ingredients,0)
LpVariable.dicts()の引数も,一番目が変数に付けられる名前,二番目が下限,三番目が上限,四番目がデータの種類. このとき,変数の名前は,"Ingr_[参照キー]"となる.例えば,MUTTONを参照キーとする変数の名前は,Ingr_MUTTONとなる.
Ingredientsの各要素に対するコストを辞書として用意する.

costs={'CHICKEN':0.013,'MUTTON':0.2,'BEEF':0.1}
この時,Ingredientsの各要素に対する変数にコストを掛けた線形関数は,次の命令で得られる.

lpSum([costs[i]*ingredients_vars[i] for i in Ingredients])

LpVariable.dicts()は辞書をつかえば,複数の添え字をもつ変数も簡単に扱える.例えば,1,2,3という3つのノードからなる完全グラフの各枝に,変数を定義するとする.これは,次の命令で実現できる.


Nodes = [1,2,3]
Arcs = [(f,t) for f in Nodes for t in Node if f!=t]
Arcs_vars = LpVariable.dicts("Arc",Arcs,lowBound=0,upBound=2,cat=LpContinuous) 
Arcsの各要素に,次のようにコストを定義することにする.

costs={(1,2):1,(1,3):2,(2,1):21,(2,3):23,(3,1):31,(3,2):32};
すると,アークに対応する変数にそのコストを掛けて和をとった線形式は,次の命令で得られる.

lpSum([Arcs_vars[i]*costs[i] for i in Arcs])
これを使えば,最短路を求めるコードは簡単に書ける(<=まちがってました.下のコードでは,部分巡回路を生成してしまう可能性があります(2014/03/19)).例えば次のように.

rom pulp import *
Nodes=[1,2,3]
Arcs=[(f,t) for f in Nodes for t in Nodes if f!=t]
Arcs_vars=LpVariable.dicts("Arc",Arcs,cat=LpBinary)
costs={(1,2):1,(1,3):2,(2,1):21,(2,3):23,(3,1):31,(3,2):32}
prob=LpProblem("Network Problem",LpMinimize)
prob+=lpSum([Arcs_vars[e]*costs[e] for e in Arcs])
source,sink=1,3;
for i in [n for n in Nodes if n!=source and n!=sink]:
        prob+=lpSum([Arcs_vars[e] for e in Arcs if e[0]==i]+
        [-1.0*Arcs_vars[e] for e in Arcs if 
        e[1]==i])==0,"flow"+str(i)
prob+=lpSum([Arcs_vars[e] for e in Arcs if 
        e[0]==source])==1,"flowsource"
prob+=lpSum([Arcs_vars[e] for e in Arcs if 
        e[1]==sink])==1,"flowsink"
print prob
prob.solve()
for v in prob.variables():
        if v.varValue>=0.01:
                print v.name,"=",v.varValue
print "cost=",value(prob.objective)

コストの他に,移動時間を考えたければ,各アークの移動時間を表す辞書を用意する. 始点sourceから終点sinkまでの総移動時間に上限Tを課したければ,timesの値を用いた制約式を追加すればよい. 今度は6つのノードからなる例を用いる.


from pulp import *
Nodes=[1,2,3,4,5,6]
Arcs=[(1,2),(1,4),(2,3),(2,4),(2,5),(3,5),(3,6),(4,5),(5,6)]
Arcs_vars=LpVariable.dicts("Arc",Arcs,cat=LpBinary)
costs={(1,2):1,(1,4):2,(2,3):3,(2,4):3,(2,5):2,(3,5):1,(3,6):3,(4,5):4,(5,6):2}
times={(1,2):2,(1,4):1,(2,3):1,(2,4):2,(2,5):2,(3,5):2,(3,6):2,(4,5):1,(5,6):5}
prob=LpProblem("Network Problem",LpMinimize)
prob+=lpSum([Arcs_vars[e]*costs[e] for e in Arcs])
source,sink=1,6;
for i in [n for n in Nodes if n!=source and n!=sink]:
        prob+=lpSum([Arcs_vars[e] for e in Arcs if e[0]==i]+
        [-1.0*Arcs_vars[e] for e in Arcs if 
        e[1]==i])==0,"flow"+str(i)
prob+=lpSum([Arcs_vars[e] for e in Arcs if 
        e[0]==source])==1,"flowsource"
prob+=lpSum([Arcs_vars[e] for e in Arcs if 
		e[1]==sink])==1,"flowsink"
prob+=lpSum([Arcs_vars[e]*times[e] for e in Arcs])<=maxT
print prob
prob.solve()
for v in prob.variables():
        if v.varValue>=0.01:
                print v.name,"=",v.varValue
print "cost=",value(prob.objective)
	
	

コメントする