Formulating CVRP with different classes of vehicles, but unlimited vehicles in each class
Awaiting user inputHello,
I'm trying to formulate a CVRP like this.
#s is the matrix of the distances between any two nodes on the graph. First row & column (index 0) are measured from the DC
s = np.array([[0.0,60.12,62.48,64.81,106.91,74.37,86.4,69.21],[60.12,0.0,13.91,4.76,52.14,15.62,29.3,9.34],[62.48,13.91,0.0,15.25,60.72,24.75,37.97,18.19],[64.81,4.76,15.25,0.0,47.87,11.14,24.93,4.6],[106.91,52.14,60.72,47.87,0.0,36.77,22.98,43.56],[74.37,15.62,24.75,11.14,36.77,0.0,13.8,6.82],[86.4,29.3,37.97,24.93,22.98,13.8,0.0,20.57],[69.21,9.34,18.19,4.6,43.56,6.82,20.57,0.0]])
#DMax is the carrying capacities of the 3 vehicle types in cbm
DMax = np.array([3.0,6.0,15.0])
#c is the matrix of converting each distance to currency (e.g. PHP) based on the vehicle type deployed between those two nodes.
c = np.array([10.0,15.0,35.0])
#cnew is the cost of deploying any new vehicle
cnew = 2000.0
#d is the demand per node of the graph
d = np.array([10,10,10,10,10,10,10])
#M is the max number of nodes allowed on one ride of one vehicle
M = 2.0
#L is the max distance which any vehicle can cover in a single ride (i.e. in one day) due to urban traffic congestion
L = 120.0
# x[k][i][j] defines whether a vehicle of the kth type traversed between nodes i & j.
#The first index/ axis (sheets before rows & columns) is the index for vehicle type.
x = model.addMVar((3,8,8),vtype = GRB.BINARY,name='x')
#v is the number of vehicles of each type used, here there are 3 types of vehicles
v = model.addMVar((3),vtype=GRB.BINARY,name='v')
X0 = model.addMVar((3,8))
#Number of vehicles of a particular type OUT from the DC = Number of vehicles of that type IN the graph
model.addConstr(np.sum(x[k],axis(0)) == X0[k] for k in range(3))
model.addConstr(X0[k][0] == v[k] for k in range (3))
#No ride from the same node to same node
model.addConstr(x[k][i][i] == 0 for i in range(8) for k in range(3))
#For all nodes which are not 0 (i.e. the DC), the total indexes incoming are 1, so that 1 vehicle of a certain type from 1 node visits this node
#Receiving is defined with columns, i.e. the 3rd axis for any node.
model.addConstr(np.sum(x,axis(1,0)) == np.array([0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]))
#The total distance travelled by any vehicle is less than the max permissible distance
model.addConstr(np.sum(np.multiply(x[k],s)) <= v[k]*L for k in range(3))
#Sum of columns in any sheet (i.e. vehicle type) gives which nodes it's servicing. This multiplied by the demand weight of that node, gives the total weight that vehicle is carrying.
#The total weight is summed per vehicle type and compared to the vehicle's carrying capacity
model.addConstr(np.sum(np.multiply(np.sum(x[k],axis(0)),d)) <= v[k]*DMax[k] for k in range(3))
#Constrain the number of drops executed by a vehicle type
model.addConstr(np.sum(x[k],axis(0)) <= v[k]*M for k in range(3))
model.setObjective(np.sum(np.multiply(x[k],s))*c[k] + cnew*v[k] for k in range(3))
model.addConstr(np.sum(x[k],axis(0))[0] == v[k] for k in range(3))

Hi Tanay,
One way you can debug this is to use an interactive Python session (such as a Jupyter notebook) and set a value for k, and see what results:
k = 0
np.sum(x[k],axis(0)) == X0[k]Then you can drill down, inspecting smaller components, such as np.sum(x[k]), and see what they result in.
Note the axis(0) is not valid, as there is no object axis defined. Perhaps you meant
np.sum(x[k], axis=0)
or
np.sum(x[k], 0)
However note that in general you can't pass Gurobi objects (like the MVar x[k]) to numpy functions  numpy does not know what to do with these things. Instead we have a number of functions defined on the MVar class, like MVar.sum(), as well as some defined in the gurobipy namespace: MatrixFriendly API Functions.
I expect you will want to use the following:
x[k].sum(axis=0) == X0[k]
Lastly, it looks like you are trying to add multiple constraints at once (for k in range(3)), in which case you need to use Model.addConstrs not Model.addConstr, or alternatively loop:
for k in range(3):
model.addConstr(x[k].sum(axis=0) == X0[k]) Riley
0 
Hello Riley and Gurobi team,
Thank you for your suggestions. I implemented a version of the MillerTuckerZemlin constraints because in the multiple classes of vehicles, the solutions had subtours. However, now the solutions I'm getting are in clear violations of the constraints.
I'm getting a solution like x[0][0][i] = 1 for all i belonging to (1,2,3...7)
By the constraint that if j!=i (lets keep j as 1 and i as 0), then
demTillNode for the 1st node should be >= demTillNode for [0] + x[k][0][1] * dem of 1st node (which is 1.2)
So demTillNode for 1st node should've been >= 1.2 at least, but here you can see that it's 0. Could someone help me in what's going wrong? Thank you!
!pip install gurobipy
import gurobipy as gp
params = {
"WLSACCESSID":"xxxxxxxxxxxxx",
"WLSSECRET":"xxxxxxxxxxxx",
"LICENSEID":xxxxx,
}
env = gp.Env(params=params)
model= gp.Model(env=env)
from gurobipy import GRB
import math
#K is the count of Vehicle Types
K= 3
#N is the universe of Nodes with [0] being the DC
N = 8
#s is the matrix of the distances between any two nodes on the graph. First row & column (index 0) are measured from the DC
s = [[0.0,60.12,62.48,64.81,106.91,74.37,86.4,69.21],[60.12,0.0,13.91,4.76,52.14,15.62,29.3,9.34],[62.48,13.91,0.0,15.25,60.72,24.75,37.97,18.19],[64.81,4.76,15.25,0.0,47.87,11.14,24.93,4.6],[106.91,52.14,60.72,47.87,0.0,36.77,22.98,43.56],[74.37,15.62,24.75,11.14,36.77,0.0,13.8,6.82],[86.4,29.3,37.97,24.93,22.98,13.8,0.0,20.57],[69.21,9.34,18.19,4.6,43.56,6.82,20.57,0.0]]
#DMax is the carrying capacities of the 3 vehicle types in cbm
DMax = [3.0,6.0,15.0]
#c is the matrix of converting each distance to currency (e.g. PHP) based on the vehicle type deployed between those two nodes.
c = [10.0,15.0,35.0]
#cnew is the cost of deploying any new vehicle
cnew = 2000.0
#d is the demand per node of the graph
d = [0,1.2,1.5,4.9,0.5,1.7,1.2,1.0]
#M is the max number of nodes allowed on one ride of one vehicle
M = 2.0
#L is the max distance which any vehicle can cover in a single ride (i.e. in one day) due to urban traffic congestion
L = 150
# x[k][i][j] defines whether a vehicle of the kth type traversed from node i to node j
x = [[[ [] for _ in range(N)] for _ in range(N)] for _ in range(K)]
for k in range(K):
for i inrange(N):
for j inrange(N):
x[k][i][j] = model.addVar(vtype=GRB.BINARY, name= f'x[{k}][{i}][{j}]')
#Need z to multiply y with cost per distance as per each vehicle type
z = [[[ [] for _ in range(N)] for _ in range(N)] for _ in range(K)]
for k in range(K):
for i inrange(N):
for j inrange(N):
z[k][i][j] = model.addVar(name= f'z[{k}][{i}][{j}]')
lenTillNode = [[[] for _ in range(N)] for _ in range(K)]
distTillNode = [[[] for _ in range(N)] for _ in range(K)]
demTillNode = [[[] for _ in range(N)] for _ in range(K)]
for k in range(K):
for i inrange(N):
lenTillNode[k][i] = model.addVar(name=f'lenTillNode[{k}][{i}]')
distTillNode[k][i] = model.addVar(name=f'distTillNode[{k}][{i}]')
demTillNode[k][i] = model.addVar(name=f'demTillNode[{k}][{i}]')
#No ride from the same node to same node and no ride to the DC
for k in range(K):
for i inrange(N):
for j inrange(N):
model.addConstr(x[k][i][i] == 0, name='Intranode travel set as zero')
model.addConstr(x[k][i][0] == 0, name='The DC has no receiving')
model.addConstr(x[k][i][j] + x[k][j][i] <= 1, name='No backhaul or looping allowed')
#For all nodes which are not 0 (i.e. the DC), the total indexes incoming are 1, so that 1 vehicle of a certain type from 1 node visits this node
for j in range(1,N):
model.addConstr(sum(x[k][i][j] for i in range(N) for k in range(K)) == 1, name='Each node touched just once')
#Define vehicletype wise variable costs based on distance travelled
for k in range(K):
for i inrange(N):
for j inrange(N):
model.addConstr(x[k][i][j] * s[i][j] * c[k] == z[k][i][j], name='Binary edges assigned with costs')
#For any edge, a vehicle of typek can traverse iff a vehicle of the same type has reached. We assume that vehicles can't change their type at nodes.
#for k in range(K):
#for i in range(1,N):
#model.addConstr(sum(x[k][i][j] for j in range (1,N)) <= sum(x[k][jdash0][i] for jdash0 in range(N)),name='Conservation of vehicles at node')
#model.addConstr(sum(x[k][i][j] for j in range(1,N)) <= 1, name='Max one vehicle exits node')
for k in range(K):
model.addConstr(lenTillNode[k][0] == 0)
model.addConstr(distTillNode[k][0] == 0)
model.addConstr(demTillNode[k][0] == 0)
for i inrange(1,N):
model.addConstr(lenTillNode[k][i] >= 1)
model.addConstr(lenTillNode[k][i] <= M)
model.addConstr(distTillNode[k][i] <= L)
model.addConstr(demTillNode[k][i] <= DMax[k])
for j inrange(1,N):
if j != i:
model.addConstr(lenTillNode[k][j] >= lenTillNode[k][i] + x[k][i][j])
model.addConstr(distTillNode[k][j] >= distTillNode[k][i] + x[k][i][j] * s[i][j])
model.addConstr(demTillNode[k][j] >= demTillNode[k][i] + x[k][i][j] * d[j])
model.setObjective(sum(z[k][i][j] for i in range(N) for j in range(N) for k in range(K)) + cnew * sum(x[k][0][i] for i in range(N)), GRB.MINIMIZE)
model.optimize()
sysNames = model.getAttr('VarName',model.getVars())
sysValues = model.getAttr('X',model.getVars())
for idx in range(len(sysNames)):
print(sysNames[idx],' ',sysValues[idx])0 
Hi Tanay,
My advice is to write your model to a LP file, and open the file in a text editor. You will want to verify that the constraint you believe to be violated is present in the model. You can then manually calculate the values of the LHS and RHS using the solution values and verify the constraint is violated. I expect one of these two steps will fail, but the learning will help fix the model.
 Riley
0
Please sign in to leave a comment.
Comments
3 comments