I'm trying to solve belows CVRPTW problem using LP by using gurobi:
The arcs represent the travel time and at the nodes we see the earliest possible arrival time, the latest possible arrival time, the service time and the demand of the customer. And important constraint is that a vehicle can arrive at a customer before its service time but then only it will have to wait until the time window opens. In that case the travel time is not determined by the arcs but by the waiting time.
From this image it becomes clear that the optimal routes are:
- Vehicle 1: [1, 3, 2]
- Vehicle 2: [4] Customer 4 can't be added to the first vehicle as the time is at: MAX(1 (0 -> 1) + 1 (service 1) + 1 (1 -> 3) + 1 (service 3) + 1 (3 -> 2), 8 (earliest 2)) = 8 Next 8 + 1 (service 2) + 1 (2 -> 4) is already 10 so at the latest service time of customer 4.
This is the class I'm using to solve the problem:
class LP_CVRPTW():
def __init__(self, df_orders, duration_matrix, capacity):
self.df_orders = df_orders
self.duration_matrix = duration_matrix
self.et = list(df_orders["earliest"])
self.lt = list(df_orders["latest"])
self.st = list(df_orders["service"])
self.depot_open = df_orders["earliest"].values[0]
self.depot_close = df_orders["latest"].values[0]
self.Demand = list(df_orders["demand"])
self.capacity = capacity
self.n = len(df_orders)
# depot, customers = coordinates[0, :], coordinates[1:, :]
self.M = 100**100 # big number
self.X = None
self.Z = None
self.Y = None
def run_model(self):
m = Model("MVRP")
x, y, z = {}, {}, {} #intialize the decision variables
'''distance matrix (32*32 array)'''
for (i, j), value in self.duration_matrix.items():
'''variable_1: X[i,j] =(0,1), i,j = Nodes'''
x[i, j] = m.addVar(vtype=GRB.BINARY, name="x%d,%d" % (i, j))
if i == j:
self.duration_matrix[(i, j)] = self.M ## big 'M'
self.duration_matrix[(i, j)] = self.M
continue
m.update()
'''variable_2: y[j] = cumulative demand covered'''
for j in range(self.n):
y[j] = m.addVar(vtype=GRB.INTEGER, name="y%d" % (j)) # cumulative demand satisfied variable
z[j] = m.addVar(vtype=GRB.INTEGER, name="z%d" % (j)) # cumulative time variable
m.update()
'''constraint_1: sum(x[i,j]) = 1, for i = 1,2,...,32'''
# vehicles leaving each customer node is equal to one
for i in range(self.n - 1):
m.addConstr(quicksum(x[(i + 1, j)] for j in range(self.n)) == 1)
m.update()
''' constraint_2: sum(x[i,j] =1 for j = 1,2,.....,32)'''
# vehicles arriving to each customer node is equal to one
for j in range(self.n - 1):
m.addConstr(quicksum(x[(i, j + 1)] for i in range(self.n)) == 1)
m.update()
# '''constraint_3: sum(x[0,j]) = 5'''
# # vehicles leaving depot is equal to no_of_vehicles
# m.addConstr(quicksum(x[(0, j)] for j in range(n)) == no_of_vehicles)
# m.update()
'''constraint_4: sum(x[i,0]) = 5'''
# vehicles arriving depot is smaller than no_of_vehicles
# m.addConstr(quicksum(x[(i, 0)] for i in range(n)) <= num_vehicles)
# m.update()
''' Either of the constraint_5 or the constrain_6 can eliminate sub-tours independently'''
'''constraint_5: capacity of vehicle and also eliminating sub-tours'''
for j in range(self.n - 1):
m.addConstr(y[j + 1] <= self.capacity)
m.addConstr(y[j + 1] >= self.Demand[j + 1])
for i in range(self.n - 1):
m.addConstr(y[j + 1] >= y[i + 1] + self.Demand[j + 1] * (x[i + 1, j + 1]) - Q * (1 - (x[i + 1, j + 1])))
m.update()
'''constraint_6: time-windows and also eliminating sub-tours'''
for i in range(self.n - 1):
# assumption: service starts at depot, depot == depot_open
# cumulative travel time should start after the earliest service start time
m.addConstr(z[i + 1] >= (self.et[i + 1] - self.depot_open))
# cumulative travel time can't be started after the latest service start time
m.addConstr(z[i + 1] <= (self.lt[i + 1] - self.depot_open))
# cumulative travel time should be between time windows of the depot
m.addConstr(z[i + 1] <= (self.depot_close - self.depot_open))
for j in range(self.n - 1):
# taking the linear distance from one node to other as travelling time in seconds between those nodes
travel_time_seconds = (self.st[j + 1] + self.duration_matrix[j + 1, i + 1]) # convert to seconds
m.addConstr(z[i + 1] >= z[j + 1] + travel_time_seconds * x[j + 1, i + 1] - self.M * (1 - x[j + 1, i + 1]))
m.update()
'''objective function'''
m.setObjective(quicksum(quicksum(x[(i, j)]*self.duration_matrix[(i, j)] for j in range(self.n)) for i in range(self.n)),GRB.MINIMIZE)
m.update()
'''optimize'''
m.optimize()
'''retrieve the solution'''
sol_y, sol_x, sol_z = m.getAttr('x', y), m.getAttr('x', x), m.getAttr('x', z)
X, Y, Z = np.empty([self.n, self.n]), np.empty([self.n]), np.empty([self.n])
for i in range(self.n):
Y[i] = sol_y[i]
Z[i] = sol_z[i]
for j in range(self.n):
X[i, j] = int(sol_x[i, j])
self.X = X
self.Y = Y
self.Z = Z
def _best_routes(self):
tours = [[i, j] for i in range(self.X.shape[0]) for j in range(self.X.shape[1]) if self.X[i, j] ==1]
routes = []
tours_0 = [t for t in tours if t[0] == 0]
for t in tours_0:
end = t[-1]
route = [end]
while end != 0:
f = [r for r in tours if r[0] == end][0]
end = f[-1]
if end != 0:
route.append(end)
routes.append(route)
return routes
def plot_tours(self):
routes = self._best_routes
tours = [[i, j] for i in range(self.X.shape[0]) for j in range(self.X.shape[1]) if self.X[i, j] ==1]
for idx, row in self.df_orders.iterrows():
plt.annotate(idx, (row["lon"], row["lat"]))
for t, tour in enumerate(tours):
plt.plot(self.df_orders["lon"][tour], self.df_orders["lat"][tour], color = "black", linewidth=0.5)
plt.scatter(self.df_orders["lon"][1:], self.df_orders["lat"][1:], marker = 'x', color = 'g', label = "customers")
plt.scatter(self.df_orders["lon"][0], self.df_orders["lat"][0], marker = "o", color = 'b', label = "depot")
plt.xlabel("lon"), plt.ylabel("lat"), plt.title("Tours"), plt.legend(loc = 4)
plt.show()
return routes
@property
def best_routes(self):
return self._best_routes()
And this is my input:
df_orders =
order_id demand date earliest latest service lat lon
0 0 0 0 2023-01-01 0 10 0 0 0
1 1 1 4 2023-01-01 0 10 1 0 0
2 2 2 2 2023-01-01 8 10 1 0 0
3 3 3 5 2023-01-01 0 10 1 0 0
4 4 4 2 2023-01-01 9 10 0 0 0
duration_matrix = {
(0, 1): 1, (1, 0): 1,
(0, 2): 1, (2, 0): 1,
(0, 3): 3, (3, 0): 3,
(0, 4): 1, (4, 0): 1,
(1, 2): 0, (2, 1): 0,
(1, 3): 1, (3, 1): 1,
(1, 4): 1, (4, 1): 1,
(2, 3): 1, (3, 2): 1,
(2, 4): 1, (4, 2): 1,
(3, 4): 4, (4, 3): 1,
(0, 0): 0, (1, 1): 0,
(2, 2): 0, (3, 3):0,
(4, 4): 0
}
# Assume same to duration_matrix
distance_matrix = duration_matrix
vrp_test6 = LP_CVRPTW(
df_orders= df_orders,
duration_matrix= duration_matrix,
capacity = 100
)
vrp_test6.run_model()
The problem is at constraint 6 but I can't seem to figure out what the problem is.
Kind regards,
Darius