Skip to main content

Formulating CVRP with different classes of vehicles, but unlimited vehicles in each class

Awaiting user input
Tanay Deshpande
  • Edited
First Comment
First Question

Hello,

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))
 
The syntax error is in the bolded rows, which don't allow me to access a specific element of the array. For example, without resorting to defining X0, even this gives an error:
model.addConstr(np.sum(x[k],axis(0))[0] == v[k] for k in range(3))
Maybe I'm having trouble in understanding how to create a temporary variable (the 1st element) based on one of the model variables (v), which then affects another model variable.
 
I found this post: https://support.gurobi.com/hc/en-us/community/posts/13087603579025-Accessing-optimization-variables-during-optimization in which the user Daniel discusses accessing values of the variables during optimization, and I understand that's not possible.
 
Could someone help on this please? Thank you!

Related to

0

Comments

3 comments

Date Votes
  • Riley Clement
    Gurobi Staff Gurobi Staff

    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: Matrix-Friendly 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
  • Tanay Deshpande
    First Comment
    First Question

    Hello Riley and Gurobi team,

     

    Thank you for your suggestions. I implemented a version of the Miller-Tucker-Zemlin 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 vehicle-type 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 type-k 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
  • Riley Clement
    Gurobi Staff Gurobi Staff

    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.

Didn't find what you were looking for?

New post