Pythonでの列生成の実行

PuLPをインストールする -> http://www.pulp.org これで混合整数計画問題が解けるようになる.
pyファイルの冒頭に
from pulp import *
と書いて,PuLPを使えるようにする.
集合分割問題 \[ \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_r} x_{rv} + y_k =1, \forall k \in K \] \[ \sum_{r \in R_v} x_{rv} \leq 1 \mbox{ and } v \in V \] \[ x_{rv} \in \{0,1\}, y_k \in \{0,1\} \] 列生成の方法を知るための数式展開は,この文書を参照-> PDF

列生成の全体の手順
全体の擬似コードは,つぎのとおり.
-----------
1   initialize()
2   Horizon = DTime(tmin,tmax,tint)
3   Networks={} #辞書
4   for v in Vehicles:
5        Networks[v],costdict[v] = genNetwork(v,Orders,Vehicles[v],Horizon)
6
7   duals_o,duals_v=masterSolve(Orders,Vehicles,Routes,True)
8
9   existsBetterRoute=True
10  while existsBetterRoute==True:
11      Routes,existsBetterRoutes=subSolve(Routes,Orders,Vehicles,
             Networks,costdict,CargosIncluded,duals_o,duals_v)
12      duals_o,duals_v=masterSolve(Orders,Vehcles,Routes,True)
13
14  masterSolve(Orders,Vehicles,Routes,False)     
1行目: initialize()で,必要な初期化を行う.

2行目は,ネットワークを定義するのに必要な,計画期間を表すHorizonを定義する. Horizonは,クラスDTimeのインスタンスである.

Vehicle vに対して
3-5行目: genNetwork() で, 部分問題を定義するネットワークを生成する. 生成したものをNetworks[v]として覚えておく.

7行目: ルートがない状況で, masterSolve()によって限定主問題(線形緩和問題)を実行する.

10-12行目: while ループで, 効率のよいルートを生成する.
追加するべき効率的なルートが見つかる限り,ループを実行する
subSolve()でネットワーク上の最短路問題を解き, 効率的なルートを求める. 追加するルートが見つかったらRoutesに追加し, existGoodRoutesにTrueを返す.
12行目:新しく加えたルートを加えて限定主問題を解く.
ルートが見つからなくなったら,whileループを終了する.

14行目: 生成したルートを使って, 変数の0-1条件をつけて集合分割問題として解く.


主問題を定義する
--
1   def masterSolve(Orders,Vehicles,Routes,relax=True):
2      prob = LpProblem("Ship Scheduling Problem",LpMinimize)
3      if relax:
4              vartype = LpContinuous
5      else:
6          vartype = LpBinary
7      spotVars = LpVariable.dicts("Spot",Orderes,0,None,vartype)
8      routeVars = LpVariable.dicts("Route",Routes,0,None,vartype)
9      spotCost=1000
10     prob += lpSum([routeVars[i]*Routes[i].cost for i in Routes] 
11                            + [spotVars[i]*spotCost for i in Orders] )
12     for o in Orders:
13             prob+=lpSum([routeVars[r] for r in Routes for to in Routes[r].orders if 
14                            o==to.name])+spotVars[o]==1,o
15     for v in Vehicles:
16             prob+=lpSum([-routeVars[r] for r in Routes 
17                            if Routes[r].Vehicle.name==v])>==-1,v
18     prob.solve()
19     duals_o,duals_v={}
20     for o in Orders:
21          duals_o[o] = prob.constraints[o].pi
22    for v in Vehicles: 
23          duals_v[v] = prob.constraints[v].pi
24     return duals_o,duals_v
--
1行目: masterSolve()では, (限定)主問題を定義して解く.主問題の定義には, オーダーの集合,船舶の集合,ルートの集合が必要なので,それぞれ Orders,Vehicles,Routesを引数とする.relaxは,緩和問題として解くか,0-1整数計画問題(集合分割問題)として解くかを指定する.

2行目
   prob = LpProblem("Ship Scheduling Problem",LpMinimize)
で,問題を定義する.2番目の引数は,最適化する方向(最大化か最小化か)を指定する.
3-6行目:
   if relax:
           vartype = LpContinuous
   else:
       vartype = LpBinary
でprobの決定変数の種類を指定する. 決定変数の整数条件を緩和したいときにはrelaxをTrueと指定する. LpContinuousは変数が連続であることを表す.relaxをFalseなら,問題の変数は0-1(LpBinary)にする.

7-8行目:
spotVars = LpVariable.dicts("Spot",Orders,0,None,vartype)
routeVars = LpVariable.dicts("Route",Routes,0,None,vartype)
で変数を生成して,spotVars,routeVarsに格納する.spotVarsは,Ordersの要素数に等しい数の変数を生成する. spotVarsは,数式のy_kにあたる.

10-11行目: 目的関数を,次で生成する.
   prob += lpSum([routeVars[r]*Routes[r].cost for r in Routes] 
                     + [spotVars[k]*spotCost for k in Orders] )
これは,数式の \[ \mbox{min} \sum_{v \in V} \sum_{r \in R_v} C_{rv} x_{rv} + \sum_{k \in K}F_{k}y_k \] にあたる. [routeVars[r]*Routes[r].cost for r in Routes]は,ルートiに対応する変数routeVars[r]に, そのルートのコストRoutes[r].costを掛けて足したもの. Routesは,数式では \[ \cup_{v \in V} R_v \] にあたる

[spotVars[k]*spotCost for k in Orders]は,Ordersの各要素iに対応する変数spotVars[k]に,そのコストspotCostをかけて足す.ここでは,コストはすべて同じとするので,同じspotCostを使っている.

12-17行目: 制約式を定義する.最初の制約式
for k in Orders:
        prob += lpSum([routesVars[r] for r in Routes for 
                      to in Routes[r].orders if k==to.name])+spotVars[k]==1,k
は,数式の \[ \mbox{s.t.} \sum_{(r,v):r \in R_v \mbox{ and } k \in R_r} x_{rv} + y_k =1, \forall k \in K \] にあたる. Routes[r].ordersには,ルートrが処理するオーダーの名前が入っている.この中に,いま考えているオーダーkがあれば,そのルートに対応する変数routesVars[r]を足す. したがって,lpSum([routeVars[r] for r in Routes for to in Routes[r].orders if k==to.name])は, オーダーkを処理するルートの変数の和を表す. 二番目の制約式
for v in Vehicles:
     prob += lpSum([-routeVars[r] for r in Routes if Routes[r].Vehicle.name==v]) >= -1,v
は,数式の \[ \sum_{r \in R_v} x_{rv} \leq 1 \mbox{ and } v \in V \] にあたる. 各船舶に対して一つ制約式が定義される. 各船舶vに対して, 船舶v用のルートの変数の和をとっている. Routes[r].Vehicle.nameには, ルートrを処理する船舶Vehicleの名前が入っているので, lpSum([-routeVars[r] for r in Routes if Routes[r].Vehicle.name==v])>=-1によって船舶vに対する 制約式を表すことができる.

ここまでで,主問題が定義できた.
定義した限定主問題を解く
   prob.solve()
で限定主問題を解く.

19-23行目: 双対変数の値をとりだし,覚えておく. この双対変数の値は, 次の部分問題の定義(subSolve())で使う.

クラスの定義
船舶のルートを表すクラスRouteの定義は次のとおり
class Route:
    def __init__(self,name,OrdersToProcess=None,Vehicle=None,cost=0):
           self.name = name;
           self.cost = cost;
           self.orders = OrdersToProcess
           self.Vehicle = Vehicle
self.ordersは,このルートが処理するオーダーの集合を記憶しておく.
self.Vehicleには,このルートを処理する船舶を表すVehicleを記憶しておく.

つぎに,船舶を表すクラスVehicleはつぎのとおり.
class Vehicle:
        def __init__(self,name,empty_speed,full_speed,capacity,init_port)
                        self.name=name
                        self.empty_speed=empty_speed
                        self.full_speed=full_speed
                        self.capacity=capacity
                        self.init_port=init_port
nameには名前, empty_speedには空船時の速度を表すリストを格納する.full_speedは, 貨物を積んだときの速度を表すリストを格納する. empty_speedよりもfull_speedの方を小さくする.例えば,
empty_speed = [20.0,19.0,18.0],
full_speed  = [18.0,17.0,16.0]
とする.

つぎに, オーダーを表すクラスOrderを定義は次のとおり
class Order:
	def __init__(self,name,load_date,disc_date,load_pt,disc_pt,quant):	
			self.name=name
			self.load_pt=load_pt
			self.disc_pt=disc_pt
			self.load_date=load_date
			self.disc_date=disc_date
			self.quant=quant
load_ptは積み港, disc_ptは揚げ港, load_dateは積み日, disc_dateは揚げ日, quantは量を表している.

部分問題を定義して解く
1 def subSolve(Routes,Orders,Vehicles,Networks,costdict,CargosIncluded,duals_o,duals_v):
2  	 moreRoutes={}
3	 existsBetterRoute=False
4	 for v in Vehicles:
5		 for (f,t) in [(tf,tt) for tf in Networks[v].nodes_iter() for tt in Networks[v][tf]]:
6			 if f=="s":
7				 Networks[v][f][t]['weight']=costdict[v][f,t]+duals_v[v]
8			 else:
9				 orders_of_tail,type_of_tail=f[0].split("_")[0],f[0].split("_")[1]
10				 if type_of_tail=="disc":
11					 Networks[v][f][t]['weight']=costdict[v][f,t]-duals_o[order_of_tail]
12		pre,dis=nx.floyd_warshall_predecessor_and_distance(Networks[v])
13		source,sink='s','t'
14		cur_nodx=sink
15		shortestpath=[]
16		shortestpath+='t'
17		while cur_node<>source:
18			cur_node=pre[source][cur_node]
19			shortestpath+=[cur_node]
20		shortestpath.reverse()
21		shortestpath_cost=dis['s']['t']
22		if shortestpath_cost<-0.0000001:
23			existsBetterRoute=True
24			cargoIncluded=[]
25			OrdersToBeProcessed=[]
26			for i in [ti for ti in shortestpath if not (ti==('s') or ti==('t'))]:
27				cargo_name = i[0].split("_")[0]
28				cargoIncluded+=[cargo_num]
29				OrdersToBeProcessed+=[Orders[cargo_name]]
30			for e in sorted(set(cargoIncluded)):
31				shortestpath_cost+=duals_o[e]
32			shortestpath_cost-=duals_v[v]
33			newRouteName="R"+str(len(Routes)+1)
34			Routes[newRouteName]=Route(newRouteName,sorted(set(OrdersToBeProcessed)),Vehicles[v],shortestpath_cost)
35			CargosIncluded[newRouteName]=sorted(set(cargoIncluded))
36	return (Routes,existsBetterRoutes)
5-11行目:ネットワークの枝のコストを, masterSolve()で得られた双対変数duals_oとduals_vの値を使って定義する.
12行目: Floyd-Warshallアルゴリズムによってソースsからシンクtへの最短路を求める.その解として, 始点から各点への最小コストの値を 入れたdisと,最短路において各ノードの一つ前のノードを表すpreを得る.
13-20行目:preの内容から最短路を作って,shortestpathに格納する.
21行目: 最短路のコストをshortestpath_costに記憶する.
22行目で, 得られた最短路のコストが負かどうかを判定する. 負であれば, その最短路を追加するための処理をする. 具体的には, subSolve()の引数Routesに, この最短路を追加する.
23-33行目:class Routeのインスタンス生成の際に与えるのは, ルートの名前, ルートで処理するオーダーの集合, ルートを処理する船舶, 最短路のコスト(双対変数を除いたもの)であるので,これらの情報を作る.
34行目:最短路を,Routeクラスのインスタンスとして生成し, Routesに追加する
36行目:subSolve()の終了にあたり, RoutesとexistsBetterRoutesを返す.

部分問題を解くためのネットワークを定義する.
部分問題は, 船ごとに定義する. これは, 時間制約,容量制約などのさまざまな制約のつく最短路問題になる. これをそのまま解くのは難しいので, 時刻を離散化して, 時空間ネットワーク上の最短路問題として定式化しなおす.これは制約なしの最短路問題となり,簡単に解ける. また, 制約条件を満たすか否かが枝の定義の有無になり,取扱いが容易になる.

時空間ネットワークのノードは, オーダーと離散時刻の組として定義される. 時空間ネットワーク上の点(k,tau1)は, オーダーkの積み港に離散時刻tau1に到着することを表す. この状態からオーダーkを処理した直後,kの揚げ港からオーダーlの積み港に移動したとする.こうしてオーダーlの積み港に離散時刻tau2に到着するとき, ノード(k,tau1)からノード(l,tau2)に枝を定義する. ただし, このときに何からの制約条件を破っているときは枝は定義しない.

まず, 計画期間を表すクラスDTimeを定義する.このクラスにはメンバ関数を3つ定義する.
class DTime:
        def ___init__(self,tmin,tmax,tint)
        def toCont(self,dt)       
        def toDisc(self,cont)
コンストラクタ__init__()の引数として,考察対象とする最も早い時刻tmin, 最も遅い時刻tmax,離散時刻の刻み幅tint,を連続時間で与える.
toCont(dt)は,離散時刻dtに対応する連続時刻を返す. toDisc(cont)は,連続時刻contに対応する離散時刻を返す.

一隻の船舶に対するネットワークを定義する関数をgenNetwork(name,Orders,Horizon,Vehicle)として定義する.ネットワークを定義するには, オーダーの集合,計画期間,船舶を引数として与える. genNetwork()の具体的な処理は次のとおりである.
def genNetwork(name,Orders,Horizon,Vehicle):
        G=nx.DiGraph()
        dist, dtime, ctime = initialize_dist(), Cargotime2()
        tnodes=[('s'),('t')]
        for (o,t) in [(ss,tt) for ss in Orders for tt in dtime.idx]:
                tnodes+=[(Orders[o].name+"_load",t),(Orders[o].name+"_disc",t)]

        for o in [to for to in Orders if Orders[to].quant<=Vehicle.capacity]:
                for i in range(5):
                        vspeed = float(Vehicle.empty_speed[i])
                        mov_time = float(dist[Vehicle.init_port,Orders[o].load_pt])/vspeed
                        arr_time = 0.0 + mov_time
                        sta_time = max(arr_time,(Orders[o].load_date-1)*24+8.0)
                        G.add_weighted_edges_from([('s', (Orders[o].name+'_load',dtime.toDisc(arr_time)-1) ,mov_time*consumption(vspeed))])

        for o in [to for to in Orders if Orders[to].quant<=Vehicle.capacity]:
                for t in [tt for tt in dtime.idx if dtime.toCont(tt)<=(Orders[o].load_date-1)*24.0+15.0]:
                        from_node = (Orders[o].name+"_load",t)
                        for i in range(5):
                                vspeed = float(Vehicle.full_speed[i])
                                disc_arr_time,cost,in_time_flag=check_if_in_time_intracargo(dtime.toCont(t),Orders[o],dist,ctime,vspeed)
                                if in_time_flag==True:
                                        to_node = (Orders[o].name+"_disc",dtime.toDisc(disc_arr_time)-1)
                                        G.add_weighted_edges_from([(from_node,to_node,cost)])

        for o in [to for to in Orders if Orders[to].quant<=Vehicle.capacity]:
                for t in [tt for tt in dtime.idx if dtime.toCont(tt)<=(Orders[o].disc_date-1)*24.0+15.0]:
                        from_node,to_node = (Orders[o].name+"_disc",t),"t"
                        end_time,in_time_flag = check_if_in_time_discharging(dtime.toCont(t),Orders[o],ctime)
                        if in_time_flag==True:
                                G.add_weighted_edges_from([(from_node,to_node,0.0)])

        for e in [te for te in itertools.permutations(Orders,2) if Orders[te[0]].quant<=Vehicle.capacity and Orders[te[1]].quant<=Vehicle.capacity]:
                from_order,to_order = e[0],e[1]
                for t in [tt for tt in dtime.idx if dtime.toCont(tt)<=(Orders[from_order].disc_date-1)*24+15.0 and dtime.toCont(tt)>=(Orders[from_order].load_date-1)*24+8.0]:
                                from_node = (Orders[from_order].name+"_disc",t)
                                from_end_time,in_time_flag=check_if_in_time_discharging(dtime.toCont(t),Orders[from_order],ctime)
                                if in_time_flag==True:
                                        for i in range(5):
                                                vspeed = float(Vehicle.empty_speed[i])
                                                mov_time = float(dist[Orders[from_order].load_pt,Orders[to_order].disc_pt])/vspeed
                                                to_arr_time = from_end_time+mov_time
                                                if to_arr_time <= (Orders[to_order].load_date-1)*24.0+15.0:
                                                        to_node = (Orders[to_order].name+"_load",dtime.toDisc(to_arr_time)-1)
                                                        G.add_weighted_edges_from([(from_node,to_node,mov_time*consumption(vspeed))])
        costdict={}
        for (n,d) in [(tn,td) for tn in G.nodes_iter() for td in G[tn]]:
                costdict[n,d]=G[n][d]['weight']
        return (G,costdict)