PDPTW model is infeasible using python gurobi
回答済みinstance from Li & Lim's PDPTW benchmark 100tasks lc101.txt
The pdptw model is referenced to Parragh-2008-A survey on pickup and delivery problems Part II Transportation between pickup and delivery locations
Below is the complete code：
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 4 10:07:39 2020
@author: AtticusYuan
"""
from gurobipy import *
import numpy as np
import pandas as pd
import math
def ReadtxtData(DataPath,Vehicles,Locations,Demand,TimeWindow,ServiceTime,Request):
# 读取车辆信息
VehiclesInfo=pd.read_table(DataPath,nrows=1,names=['K','C','S'])
#print('VehiclesInfo:')
#print(VehiclesInfo)
for i in range(VehiclesInfo.iloc[0,0]):
Vehicles[i]=[VehiclesInfo.iloc[0,1],VehiclesInfo.iloc[0,2]] # 键i=车辆的序号，值依次为车辆容量、速度
#print(Vehicles)
# 读取Depot和任务信息
ColumnNames=['TaskNo','X','Y','Demand','ET','LT','ST','PI','DI']
TaskData=pd.read_table(DataPath,skiprows=[0],names=ColumnNames,index_col='TaskNo')
#print('TaskData:')
#print(TaskData)
# 提取Depot和取送货点的位置坐标Locations
nrows=TaskData.shape[0]
count = 0
for i in range(nrows):
if i not in Locations.keys():
Locations[i]=[count,TaskData.iloc[i,0],TaskData.iloc[i,1]]
count += 1
#print('Locations:',Locations)
# 提取Depot和取送货点的Demand
for i in range(nrows):
if i not in Demand.keys():
Demand[i]=TaskData.iloc[i,2]
#print('Demand',Demand)
# 提取Depot和取送货点的最早和最晚取送货时间及时间窗
EarliestTime=TaskData.sort_values(by='ET').iloc[0,3]
LatestTime=TaskData.sort_values(by='LT',ascending=False).iloc[0,4]
for i in range(nrows):
if i not in TimeWindow.keys():
TimeWindow[i]=[TaskData.iloc[i,3],TaskData.iloc[i,4]]
# print('TimeWindow',TimeWindow)
print('EarliestTime:',EarliestTime)
print('LatestTime:',LatestTime)
# 提取Depot和取送货点的服务时间ServiceTime
for i in range(nrows):
if i not in ServiceTime.keys():
ServiceTime[i]=TaskData.iloc[i,5]
#print('ServiceTime',ServiceTime)
# 提取运输Request
count = 0
for i in range(1,nrows):
if TaskData.iloc[i,6] == 0:
Request[count]=[i,TaskData.iloc[i,7]] # 将取送货点组合在一起
count += 1
#print('Request:')
#print(Request)
return [nrows,LatestTime]
def CalculateEuclidDistance(x1, y1, x2, y2): # 计算两点(x1, y1), (x2, y2)之间的Euclid距离
arr1=np.array(x1,y1)
arr2=np.array(x2,y2)
EuclidDistance=np.sqrt(np.sum(np.square(arr1-arr2)))
return EuclidDistance
def ConstructDistanceMatrix(Locations, DistanceMatrix): # 建立距离矩阵
for i in Locations.keys():
for j in Locations.keys():
DistanceMatrix[i, j] = CalculateEuclidDistance(Locations[i][1], Locations[i][2],
Locations[j][1], Locations[j][2])
# print('DistanceMatrix:')
# print(DistanceMatrix)
# 尝试获取字典DistanceMatrix中最大的值
key=max(DistanceMatrix, key=DistanceMatrix.get)
LongestDistance=DistanceMatrix[key]
# 以下尝试获取字典DistanceMatrix中最大的值
# SortedDistanceMatrix=sorted(DistanceMatrix.items(), key=lambda kv:(kv[1], kv[0]),reverse=True)
# i=0
# for key,value in SortedDistanceMatrix.items():
# while i<1:
# LongestDistance=value
# i+=1
print('LongestDistance:',LongestDistance)
return LongestDistance
def BuildGurobiModel(Vehicles,Locations,Demand,TimeWindow,ServiceTime,Request,DistanceMatrix,xsol, bsol, qsol):
qindex = {} # 存储变量qik的下标ik，表示车辆k即将离开i时的载货量
bindex = {} # 存储变量bik的下标ik，表示车辆k开始服务i的时间
TimeMatrix = {} # 节点之间车辆的行驶时间
xobjcof= {} # 变量xijk在目标函数中的系数,也就是距离
for k in Vehicles.keys():
for i in range(nrows):
qindex[i,k]=0
bindex[i,k]=0
for j in range(nrows):
TimeMatrix[i,j,k]=DistanceMatrix[i,j]/Vehicles[k][1]
xobjcof[i,j,k]=DistanceMatrix[i,j]
# 创建模型
model = Model()
# 创建变量
x = model.addVars(xobjcof.keys(), obj=xobjcof, vtype=GRB.BINARY, name='x') # 变量x_{ijk}
q = model.addVars(qindex.keys(), vtype=GRB.CONTINUOUS, name='q') # 变量q_{ik}
b = model.addVars(bindex.keys(), vtype=GRB.CONTINUOUS, name='b') # 变量b_{ik}
# (1) every vertex has to be served exactly once，除了depot外，每个节点只被服务一次
for i in range(1, nrows):
model.addConstr(x.sum(i,'*','*') == 1) # 星号*的用法
# (2-3) guarantee that every vehicle starts at the depot and returns to the depot
# at the end of its route. 每辆车必须从depot出发，最后回到depot
# (4) Flow conservation 流平衡约束
for k in Vehicles.keys():
model.addConstr(x.sum(0,'*',k) == 1) # 车辆k的起始停靠点
model.addConstr(x.sum('*',0,k) == 1) # 车辆k的结束停靠点
for i in range(1, nrows): # 车辆k在P和D的流平衡约束
model.addConstr(x.sum(i,'*',k) == x.sum('*',i,k))
# (5) Time variables are used to eliminate subtours,用时间变量来消除子回路约束
# 需要用大M法来线性化该约束，M定义为 2*(LatestTime+LongestDistance)
for i in xobjcof.keys():
model.addConstr(b[i[1],i[2]]>=b[i[0],i[2]]+ServiceTime[i[0]]+TimeMatrix[i]+
2*(LatestTime+LongestDistance)*x[i]-2*(LatestTime+LongestDistance))
# (6-7) guarantee that a vehicle’s capacity is not exceeded throughout its tour，载货量平衡与车辆载量约束
# (6) 载货量平衡约束，需要用大M法来线性化该约束，M定义为 100*车辆最大载量
for i in xobjcof.keys():
model.addConstr(q[i[1],i[2]]>=q[i[0],i[2]]+Demand[i[1]]+
x[i]*(100*Vehicles[i[2]][0])-(100*Vehicles[i[2]][0]))
# （7）车辆载量约束
for i in qindex.keys():
model.addConstr(q[i]>=0)
model.addConstr(q[i]>=Demand[i[0]])
model.addConstr(q[i]<=Vehicles[i[1]][0])
model.addConstr(q[i]<=Vehicles[i[1]][0]+Demand[i[0]])
# (8) both origin and destination of a request must be served by the same vehicle
# 保证取货后要有对应的送货，取货和送货由同一辆车完成
for i in range(len(Request)):
for j in range(len(Vehicles)):
model.addConstr(x.sum(Request[i][0],'*',j)==x.sum('*',Request[i][1],j))
# (9) delivery can only occur after pickup,先取后送货约束
for i in range(len(Request)):
for j in range(len(Vehicles)):
model.addConstr(b[Request[i][0],j]<=b[Request[i][1],j])
# (10-11) 节点时间窗约束，左时间窗，右时间窗
for i in range(len(Request)):
for j in range(len(Vehicles)):
# 运输请求i两个节点的左时间窗
model.addConstr(b[Request[i][0],j]>=TimeWindow[Request[i][0]][0])
model.addConstr(b[Request[i][1],j]>=TimeWindow[Request[i][1]][0])
# 节点右时间窗
# 运输请求i两个节点的右时间窗
model.addConstr(b[Request[i][0],j]<=TimeWindow[Request[i][0]][1])
model.addConstr(b[Request[i][1],j]<=TimeWindow[Request[i][1]][1])
# 以下约束可能导致模型不可行，因为bik是车辆k开始服务节点i的时间
#model.addConstr(b[Request[i][1],j]+ServiceTime[Request[i][1]]<=TimeWindow[Request[i][0]][1])
model.params.TimeLimit = 600 # 设定模型求解时间600s
model.write('PDPTWTestDataSolvedbyGurobi.lp')
model.write('PDPTWTestDataSolvedbyGurobi.mps')
model.setParam(GRB.Param.LogFile, 'PDPTWTestDataSolvedbyGurobi.log')
model.optimize()
for i in xobjcof.keys():
if x[i].x -0.5>0:
xsol[i] = 1
bsol[i[0], i[2]]=b[i[0],i[2]].x
bsol[i[1], i[2]]=b[i[1],i[2]].x
qsol[i[0], i[2]]=q[i[0],i[2]].x
qsol[i[1], i[2]]=q[i[1],i[2]].x
print('xsol:')
print(xsol)
print('bsol:')
print(bsol)
print('qsol:')
print(qsol)
try:
#读取数据
DataPath = 'lc101.txt'
Vehicles={}
Locations={}
Demand={}
TimeWindow={}
ServiceTime={}
Request={}
Val1=ReadtxtData(DataPath,Vehicles,Locations,Demand,TimeWindow,ServiceTime,Request)
nrows=Val1[0]
LatestTime=Val1[1]
DistanceMatrix = {}
LongestDistance=ConstructDistanceMatrix(Locations, DistanceMatrix)
xsol = {}
bsol = {}
qsol = {}
#建立模型
BuildGurobiModel(Vehicles,Locations,Demand,TimeWindow,ServiceTime,Request,DistanceMatrix,xsol, bsol, qsol)
# OutputResult(Locations, Vehicles, Orders, EarliestDate, xsol, asol, qsol, wsol)
except GurobiError as exception:
print('Error code ' + str(exception.errno) + ": " + str(exception))
except AttributeError:
print('Encountered an attribute error')
-
正式なコメント
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
-
-
Hi Xin,
I would try computing an IIS to determine which constraints are contributing to the infeasibility. This should help you identify if there is a problem with the constraints, or even the data itself. In Python:
model.computeIIS()
model.write('model.ilp')The resulting ILP file contains a subset of constraints/bounds that together are infeasible. However, removing any one of the constraints/bounds from this set results in a feasible subsystem.
Thanks,
Eli
0 -
-
Hi Eli,
Thank you for your help! I'll try to follow your advice now. May I ask you again if I have any questions?
Thanks,
Xin
0 -
-
-
Hi Xin,
Sure, I try to answer questions here when I have some spare time.
Eli
0 -
-
Would it be possible to reshare the code completely in English?
0 -
投稿コメントは受け付けていません。
コメント
5件のコメント