MDVRPTW问题，但是是technician routing scheduling problem，设置有6个技术人员、6个客户以及两个中转站，包含时间窗约束，以下是表格数据和完整代码。
'''Multi-Depot Vehicle Route Problem With Time Window 带时间窗口的多网点车辆路由问题 2024/5/15 by sadbooooi'''
import sys
from collections import defaultdict
import xlrd
import gurobipy as gp
from gurobipy import GRB
# 创建技术人员类
class Technician():
def __init__(self, name, cap, depot):
self.name = name
self.cap = cap
self.depot = depot
def __str__(self):
return f"Technician: {self.name}\n Capacity: {self.cap}\n Depot: {self.depot}"
# 创建工作类
class Job():
def __init__(self, name, priority, duration, coveredBy):
self.name = name
self.priority = priority
self.duration = duration
self.coveredBy = coveredBy
def __str__(self):
about = f"Job: {self.name}\n Priority: {self.priority}\n Duration: {self.duration}\n Covered by: "
about += ", ".join([t.name for t in self.coveredBy])
return about
# 创建客户类
class Customer():
def __init__(self, name, loc, job, tStart, tEnd, tDue):
self.name = name
self.loc = loc
self.job = job
self.tStart = tStart
self.tEnd = tEnd
self.tDue = tDue
def __str__(self):
coveredBy = ", ".join([t.name for t in self.job.coveredBy])
return (f"Customer: {self.name}\n Location: {self.loc}\n Job: {self.job.name}\n Priority: {self.job.priority}\n "
f" Duration: {self.job.duration}\n Covered by: {coveredBy}\n Start time: {self.tStart}\n "
f"End time: {self.tEnd}\n Due time: {self.tDue}")
# 打开Excle工作表格
wb = xlrd.open_workbook(r'C:\Users\Administrator\Desktop\modeling-examples-master\technician_routing_scheduling\data-Sce0.xls')
# 打开名为data-Sce0.xlsx的Excel工作簿
# 读取技术人员数据
ws = wb.sheet_by_name('Technicians') # 从工作簿wb中获取名为“Technicians”的工作表
technicians = []
for i,t in enumerate(ws.col_values(0)[3:]):
# 创建技术人员对象
thisTech = Technician(*ws.row_values(3+i)[:3])
technicians.append(thisTech)
# 读取工作数据
jobs = []
for j, b in enumerate(ws.row_values(0)[3:]):
coveredBy = [t for i, t in enumerate(technicians) if ws.cell_value(3 + i, 3 + j) == 1]
# 遍历之前创建的技术人员列表technicians，并检查在技术人员行和工作列交叉的单元格（即第i+3行和第j+3列）的值是否为1，如果为1，则表示该技术人员负责这项工作，将其添加到coveredBy列表中
# 创建工作对象
thisJob = Job(*ws.col_values(3 + j)[:3], coveredBy)
jobs.append(thisJob)
# 读取位置数据
ws = wb.sheet_by_name('Locations')
locations = ws.col_values(0)[1:]
dist = {(l, l) : 0 for l in locations} # 对于每个位置l，字典中都会有一个键值对(l, l): 0，表示位置到自身的距离是0
for i,l1 in enumerate(locations):
for j,l2 in enumerate(locations):
if i < j: # 因为距离是对称的，所以只需要计算上三角或下三角的部分
dist[l1,l2] = ws.cell_value(1+i, 1+j)
dist[l2,l1] = dist[l1,l2]
# 读取客户数据
ws = wb.sheet_by_name('Customers')
customers = []
for i,c in enumerate(ws.col_values(0)[1:]):
for b in jobs:
if b.name == ws.cell_value(1+i, 2):
# 使用相应的工作对象创建客户对象
rowVals = ws.row_values(1+i)
# 打印行数据
thisCustomer = Customer(*rowVals[:2], b, *rowVals[3:])
customers.append(thisCustomer)
break
def solve_trs0(technicians, customers, dist):
# 建立有用的数据结构
K = [k.name for k in technicians]
C = [j.name for j in customers]
J = [j.loc for j in customers]
L = list(set([l[0] for l in dist.keys()]))
D = list(set([t.depot for t in technicians]))
cap = {k.name : k.cap for k in technicians}
loc = {j.name : j.loc for j in customers}
depot = {k.name : k.depot for k in technicians}
canCover = {j.name : [k.name for k in j.job.coveredBy] for j in customers}
dur = {j.name : j.job.duration for j in customers}
tStart = {j.name : j.tStart for j in customers}
tEnd = {j.name : j.tEnd for j in customers}
tDue = {j.name : j.tDue for j in customers}
priority = {j.name : j.job.priority for j in customers}
### 建立模型
m = gp.Model("trs0")
### 决策变量
# j是否由技术人员k指派完成
x = m.addVars(C, K, vtype=GRB.BINARY, name="x")
# 技术人员k是否执行一项工作
u = m.addVars(K, vtype=GRB.BINARY, name="u")
# 技术人员k是否从i前往j执行指派
y = m.addVars(L, L, K, vtype=GRB.BINARY, name="y")
# 技术员只能从同一个仓库往返
for k in technicians:
for d in D:
if k.depot != d:
for i in L:
y[i, d, k.name].ub = 0
y[d, i, k.name].ub = 0
# 服务开始时间
t = m.addVars(L, ub=600, name="t")
# 服务的延迟
z = m.addVars(C, name="z")
# 修正时间窗口上下限的人工变量
xa = m.addVars(C, name="xa")
xb = m.addVars(C, name="xb")
# 未完成的工作
g = m.addVars(C, vtype=GRB.BINARY, name="g")
### 添加约束
# 对于每一项工作，有且仅有一名技术人员完成 （1 and 2）
m.addConstrs((gp.quicksum(x[j, k] for k in canCover[j]) + g[j] == 1 for j in C), name="assignToJob")
m.addConstrs((x.sum(j, '*') <= 1 for j in C), name="assignOne")
# 确保不超过每一个技术人员的可用力，即时间容量 （3）
capLHS = {k: gp.quicksum(dur[j] * x[j, k] for j in C) + \
gp.quicksum(dist[i, j] * y[i, j, k] for i in L for j in L) for k in K} # 建立以技术人员k为键的字典capLHS存储其作业时间和通行时间
m.addConstrs((capLHS[k] <= cap[k] * u[k] for k in K), name="techCapacity")
# 确保技术人员分配到工作后前往所在工作地(4 and 5)
m.addConstrs((y.sum('*', loc[j], k) == x[j, k] for k in K for j in C), \
name="techTour1")
m.addConstrs((y.sum(loc[j], '*', k) == x[j, k] for k in K for j in C), \
name="techTour2")
# 对于每一个技术人员和仓库，确保其从一个仓库出发并返回(6 and 7)
m.addConstrs((gp.quicksum(y[j, depot[k], k] for j in J) == u[k] for k in K), \
name="sameDepot1")
m.addConstrs((gp.quicksum(y[depot[k], j, k] for j in J) == u[k] for k in K), \
name="sameDepot2")
# 时间约束，技术人员k从工作i到j，那么j的服务开始时间必须大于i的完成时间加上i到j的通行时间 （8） ?????????
M = {(i, j): 600 + dur[i] + dist[loc[i], loc[j]] for i in C for j in C}
m.addConstrs((t[loc[j]] >= t[loc[i]] + dur[i] + dist[loc[i], loc[j]] \
- M[i, j] * (1 - gp.quicksum(y[loc[i], loc[j], k] for k in K)) \
for i in C for j in C), name="tempoCustomer")
# t[loc[j]]表示客户j所在位置的服务开始时间，这个时间变量t[loc[j]]是通过模型中的约束条件动态确定的，而不需要显式地定义
M = {(i, j): 600 + dist[i, loc[j]] for i in D for j in C}
m.addConstrs((t[loc[j]] >= t[i] + dist[i, loc[j]] \
- M[i, j] * (1 - y.sum(i, loc[j], '*')) for i in D for j in C), \
name="tempoDepot") # 技术人员从仓库到客户位置的行程
# 时间窗口约束 （9 and 10）
m.addConstrs((t[loc[j]] + xa[j] >= tStart[j] for j in C), name="timeWinA")
m.addConstrs((t[loc[j]] - xb[j] <= tEnd[j] for j in C), name="timeWinB")
# 延迟约束 (11)
m.addConstrs((z[j] >= t[loc[j]] + dur[j] - tDue[j] for j in C), \
name="lateness")
### 目标函数
M = 6100
m.setObjective(z.prod(priority) + gp.quicksum(0.01 * M * priority[j] * (xa[j] + xb[j]) for j in C) +
gp.quicksum(M * priority[j] * g[j] for j in C), GRB.MINIMIZE)
m.write("TRS0.lp") # 通过m.write("TRS0.lp")将模型写入LP文件，以便查看模型的具体细节
m.optimize()
status = m.Status
if status in [GRB.INF_OR_UNBD, GRB.INFEASIBLE, GRB.UNBOUNDED]:
print("Model is either infeasible or unbounded.")
sys.exit(0) # 终止当前的 Python 脚本并返回一个退出代码
elif status != GRB.OPTIMAL:
print("Optimization terminated with status {}".format(status))
sys.exit(0)
### 打印结果
# 指派
print("")
for j in customers:
if g[j.name].X > 0.5: # .X用于返回决策变量中的最优解
jobStr = "Nobody assigned to {} ({}) in {}".format(j.name, j.job.name, j.loc)
else:
for k in K:
if x[j.name, k].X > 0.5:
jobStr = "{} assigned to {} ({}) in {}. Start at t={:.2f}.".format(k, j.name, j.job.name, j.loc,
t[j.loc].X)
if z[j.name].X > 1e-6:
jobStr += " {:.2f} minutes late.".format(z[j.name].X)
if xa[j.name].X > 1e-6:
jobStr += " Start time corrected by {:.2f} minutes.".format(xa[j.name].X)
if xb[j.name].X > 1e-6:
jobStr += " End time corrected by {:.2f} minutes.".format(xb[j.name].X)
print(jobStr)
# 技术人员
print("")
for k in technicians:
if u[k.name].X > 0.5:
cur = k.depot
route = k.depot
while True:
for j in customers:
if y[cur, j.loc, k.name].X > 0.5:
route += " -> {} (dist={}, t={:.2f}, proc={})".format(j.loc, dist[cur, j.loc], t[j.loc].X,
j.job.duration)
cur = j.loc
for i in D:
if y[cur, i, k.name].X > 0.5:
route += " -> {} (dist={})".format(i, dist[cur, i])
cur = i
break
if cur == k.depot:
break
print("{}'s route: {}".format(k.name, route))
else:
print("{} is not used".format(k.name))
# 利用率
print("")
for k in K:
used = capLHS[k].getValue()
total = cap[k]
util = used / cap[k] if cap[k] > 0 else 0
print("{}'s utilization is {:.2%} ({:.2f}/{:.2f})".format(k, util, \
used, cap[k]))
totUsed = sum(capLHS[k].getValue() for k in K)
totCap = sum(cap[k] for k in K)
totUtil = totUsed / totCap if totCap > 0 else 0
print("Total technician utilization is {:.2%} ({:.2f}/{:.2f})".format(totUtil, totUsed, totCap))
def printScen(scenStr):
sLen = len(scenStr)
print("\n" + "*"*sLen + "\n" + scenStr + "\n" + "*"*sLen + "\n")
if __name__ == "__main__":
# Base model
printScen("Solving base scenario model")
solve_trs0(technicians, customers, dist)
673
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
