# Two-type VRP

I am trying to formulate a VRP for two types of vehicles and a single depot. V1 represents, say EVs, and V2 represents gas-fueled vehicles. V1-type vehicles consume less fuel/energy in the first M time units. I have the following Gurobipy programming. The problem is that I get solutions with incomplete/disconnected tours, mainly because a vehicle traverses i to j and the reverse to satisfy the connectivity constraints. Although I added the MTZ sub-tour elimination, I still end up receiving such solutions. When I introduce the constraint: x[v,i,j] + x[v,j,i] == 1, I get an infeasible solution. What is wrong in this formulation?

Here are the set, parameter, and variable definitions. Set I represents all customer locations, Depot_loc is the depot location. All locations:= I \cup Depot_loc = J. V1 is the set of EV vehicles, V2 is the set of gas vehicles, and V1 \cup V2 = V. T[i,j] represents the travel time from i to j, and is asymmetric. U is the UB of travel time for each vehicle. Variables are x, y, t, t_xi, t_kappa. x[v,i,j] = 1 if v takes i,j. y[v,j] is the sub-tour elimination variable. Variable t[v] records the travel time for V2 type vehicle. t_xi[v] + t_kappa[v] records the travel time for V1 type vehicle, where t_xi[v] accounts for the cheaper first M time units and t_kappa[v] for the rest.

m = Model("Problem")

x = {}; y = {}

for v in V:

for i in J:

for j in J:

x[v,i,j] = m.addVar(vtype=GRB.BINARY, name="x%s"%str([v,i,j]))for v in V:

for j in I:

y[v,j] = m.addVar(vtype=GRB.CONTINUOUS, name = "y%s"%str([v,j]))

t = {v: m.addVar(vtype=GRB.CONTINUOUS, name="t%s"%str([v])) for v in V2}

t_xi = {v: m.addVar(vtype=GRB.CONTINUOUS, name="t_xi%s"%str([v])) for v in V1}

t_kappa = {v: m.addVar(vtype=GRB.CONTINUOUS, name="t_kappa%s"%str([v])) for v in V1}

m.setObjective(quicksum(First_cost*t_xi[v] +

Second_cost*t_kappa[v] for v in V1) +

quicksum(Third_cost*t[v] for v in V2), GRB.MINIMIZE)

m.addConstrs((quicksum(x[v,Depot_loc,j] for j in I) == 1 for v in V), name = "every_vehicle_starts_at_depot")

m.addConstrs((quicksum(x[v,i,Depot_loc] for i in I) == 1 for v in V), name = "every_vehicle_ends_at_depot")

m.addConstrs((quicksum(x[v,i,j] for v in V for j in J if i != j) >= 1 for i in I),name = "every_node_to_be_visited")

m.addConstrs((quicksum(T[i,j]*x[v,i,j] for i in J for j in J if i != j) <= U for v in V), name = "limited_travel_time")

m.addConstrs((t_xi[v]+t_kappa[v] == quicksum(T[i,j]*x[v,i,j]

for i in J for j in J if i != j) for v in V1), name = "travel_time_for_V1")

m.addConstrs((t[v] == quicksum(T[i,j]*x[v,i,j] for i in J for j in J if i != j) for v in V2), name = "travel_time_for_V2")

m.addConstrs((t_xi[v] <= M for v in V1), name = "limit_electric_travel")

m.addConstrs((quicksum(x[v,i,j] for i in I if i != j) == quicksum(x[v,j,j_prime] for j_prime in I if j_prime != j) for v in V for j in I), name = "connectivity")

m.addConstrs((y[v,j] >= y[v,i]-len(J)*(1-x[v,i,j])

for v in V for i in I for j in I if i != j), name = "subtour_elimination")

### BEGIN SOLVING ###

m.setParam('TimeLimit', 720000); m.setParam('MIPgap', 0.0005)

m.update(); m.optimize(); print("\n\n")

### END SOLVING ###

Please sign in to leave a comment.

## Comments

0 comments