python+gurobi 带时间窗带容量限制的路径规划 cvrptw

import numpy as np
rnd = np.random
rnd.seed(0)
import random
random.seed(0)
M=1000

# data

nbX = 5
source = 0  #出发点
puits = nbX  #回去点，共有nbx-1个顾客，2个车站

# sets
V = range(3)  #vehicles
X = list(range(nbX)) #clients
U = [(x1, x2) for x1 in X for x2 in X] #arcs

# 两点之间的距离
random.seed(0)
D = {u: int(random.uniform(1, 10)) for u in U}
# diagonale vide
for x in X:
    D[(x, x)] = 0
# inegalité triangulaire
for j in X:
    for i in X:
        for k in X:
            D[(i, k)] = min(D[(i, k)], D[(i, j)] + D[(j, k)])

# fenetres de temps
E  = [max(D[(source, x)], int(random.uniform(1, 5 * nbX))) for x in X]
L = [E[x] + int(random.uniform(20, 40)) for x in X]

E[source] = 0
Tmax = max(L) + max(D.values())
L[source] = Tmax

for x in X:
    U.append((x, puits))
    D[(x, puits)] = D[(source, x)]

X.append(puits)
E.append(E[source])
L.append(L[source])

# capacity
Q = 20
q = {i: rnd.randint(1, 10) for i in X}
q[source]=0
q[puits]=0

# service time
service_time = {i: rnd.randint(1, 5) for i in X}
service_time[source]=0
service_time[puits]=0

# optimization
from gurobipy import Model, GRB, quicksum
mdl = Model('CVRPTW')

# decision variables
x = mdl.addVars(U,V, vtype=GRB.BINARY)
T = {}
for i in X:
    # creat variablesT
    for k in V:
        T[i, k] = mdl.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='T_' + str(i) + '_' + str(k))

# #objective
mdl.modelSense = GRB.MINIMIZE
mdl.setObjective(quicksum(x[i, j,k]*D[i, j] for i,j in U for k in V))
#
# #constraints
mdl.addConstrs(quicksum(x[i,j,k] for k in V for j in range(1,nbX+1) if j != i) == 1 for i in range(1,nbX))
mdl.addConstrs(quicksum(x[0,j,k] for j in range(1,nbX)) ==1 for k in V )
mdl.addConstrs(quicksum(x[i,nbX,k] for i in range(1,nbX)) == 1 for k in V)
mdl.addConstrs((quicksum(x[i,j,k] for i in range(nbX) if j!=i)) ==
               (quicksum(x[j,i,k] for i in range(1,nbX+1) if j!=i))  for j in range(1,nbX) for k in V)
mdl.addConstrs(T[i,k]+D[i,j]+service_time[i]-T[j,k]<= M*(1-x[i,j,k]) for i,j in U for k in V)
mdl.addConstrs(T[i,k]<=L[i] for i in range(nbX+1) for k in V)
mdl.addConstrs(quicksum(q[i]*x[i,j,k] for i in range(nbX) for j in range(nbX+1) if i!=j)<=Q for k in V)
mdl.addConstrs(T[i,k]>=E[i] for i in range(nbX+1) for k in V)


#results
mdl.Params.MIPGap = 0.1
mdl.Params.TimeLimit = 30  # seconds
mdl.optimize()

A = [(x1,x2,x3) for x1,x2 in U for x3 in V]
active_arcs = [a for a in A if x[a].x > 0.99]

print("\n\n-----optimal value-----")
print(mdl.ObjVal)
print(active_arcs)
for key in T.keys():
    print(T[key],"=",T[key].x)

车辆可以早于时间窗左端到达，将超过时间窗左端的时间，作为惩罚项。车辆不能晚于时间窗右端。

带惩罚项
带容量限制
带开放的时间窗

import numpy as np
rnd = np.random
rnd.seed(0)
import random
random.seed(0)
M=1000

# data

nbX = 5
source = 0  #出发点
puits = nbX  #回去点，共有nbx-1个顾客，2个车站

# sets
V = range(3)  #vehicles
X = list(range(nbX)) #clients
U = [(x1, x2) for x1 in X for x2 in X] #arcs

# 两点之间的距离
random.seed(0)
D = {u: int(random.uniform(1, 10)) for u in U}
# diagonale vide
for x in X:
    D[(x, x)] = 0
# inegalité triangulaire
for j in X:
    for i in X:
        for k in X:
            D[(i, k)] = min(D[(i, k)], D[(i, j)] + D[(j, k)])

# fenetres de temps
E  = [max(D[(source, x)], int(random.uniform(1, 5 * nbX))) for x in X]
L = [E[x] + int(random.uniform(20, 40)) for x in X]

E[source] = 0
Tmax = max(L) + max(D.values())
L[source] = Tmax

for x in X:
    U.append((x, puits))
    D[(x, puits)] = D[(source, x)]

X.append(puits)
E.append(E[source])
L.append(L[source])

# capacity
Q = 20
q = {i: rnd.randint(1, 10) for i in X}
q[source]=0
q[puits]=0

# service time
service_time = {i: rnd.randint(1, 5) for i in X}
service_time[source]=0
service_time[puits]=0
from gurobipy import Model, GRB, quicksum, max_

mdl = Model('CVRPTW')
# decision variables
x = mdl.addVars(U,V, vtype=GRB.BINARY)
T = {}
for i in X:
    # creat variablesT
    for k in V:
        T[i, k] = mdl.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='T_' + str(i) + '_' + str(k))
p = mdl.addVars(range(1,nbX),V,vtype=GRB.CONTINUOUS)#1,2,3,4
# y = mdl.addVars(range(1,nbX),V,vtype=GRB.CONTINUOUS)
z = mdl.addVars(range(1,nbX),V,vtype=GRB.CONTINUOUS)

print(p)


#
# objective
mdl.modelSense = GRB.MINIMIZE
mdl.setObjective( quicksum(x[i, j,k]*D[i, j] for i,j in U for k in V)+quicksum(p[i,k] for i in range(1,nbX) for k in V))
#
# #constraints


mdl.addConstrs(quicksum(x[i,j,k] for k in V for j in range(1,nbX+1) if j != i) == 1 for i in range(1,nbX))
mdl.addConstrs(quicksum(x[0,j,k] for j in range(1,nbX)) ==1 for k in V )
mdl.addConstrs(quicksum(x[i,nbX,k] for i in range(1,nbX)) == 1 for k in V)
mdl.addConstrs((quicksum(x[i,j,k] for i in range(nbX) if j!=i)) ==
               (quicksum(x[j,i,k] for i in range(1,nbX+1) if j!=i))  for j in XX for k in V)
mdl.addConstrs(T[i,k]+D[i,j]+service_time[i]-T[j,k]<= M*(1-x[i,j,k]) for i,j in U for k in V)

mdl.addConstrs(T[i,k]<=L[i] for i in range(nbX+1) for k in V)
mdl.addConstrs(quicksum(q[i]*x[i,j,k] for i in range(nbX) for j in range(nbX+1) if i!=j)<=Q for k in V)

mdl.addConstrs(z[i,k] == quicksum(x[i,j,k] for j in range(1,nbX+1))*E[i]-T[i,k] for i in range(1,nbX) for k in V )

for i in range(1,nbX):
    for k in V:
        mdl.addConstr(p[i,k]==max_(0,z[i,k]))



mdl.Params.MIPGap = 0.1
mdl.Params.TimeLimit = 60  # seconds
mdl.optimize()
print(mdl.ObjVal)
A = [(x1,x2,x3) for x1,x2 in U for x3 in V]
print(active_arcs)
for key in T.keys():
    print(T[key],"=",T[key].x)