Skip to main content

Constraint indexes are not read by the model

Answered

Comments

3 comments

  • Maliheh Aramon
    Gurobi Staff Gurobi Staff
    Hi Benay, 
     
    The real issue is that the constraint below does not guarantee that the consecutive visits to node \(i\) are at least 20 days apart. The constraint below just ensures that if node \(i\) is visited at time \(t\), it will be visited at time \(t+20\) as well. It does not guarantee anything for the intermediate times between time \(t\) and time \(t+20\).
    mdl.addConstr(v[i, t+20] >= v[i, t])

    Given your explanation, if node \(i\) is visited at time \(t\), i.e., \(v_{it} =1\), it should not be visited in the next 20 days, i.e., \(\sum_{j = t+1}^{t+20} v_{ij} = 0\).

    You can use the Gurobi indicator constraint API to implement:

    \[\mbox{if  } v_{it} =1 \Rightarrow \sum_{j = t+1}^{t+20} v_{ij} = 0\]

    Specifically, you can replace the above constraint with the constraint below in your code snippet:

    start_time, end_time = T[0], T[-1] + 1
    mdl.addConstrs(
    (v[i, t] == 1)
    >> (
    quicksum(
    v[i, j] for j in range(max(t + 1, start_time), min(t + 20, end_time))
    ) == 0
    )
    for i in I2
    for t in T
    )

    I tested this with your code snippet and I got:

    <gurobi.Var v[1,1] (value 1.0)>
    <gurobi.Var v[2,1] (value 1.0)>
    <gurobi.Var v[2,39] (value 1.0)>
    <gurobi.Var v[3,1] (value 1.0)>
    <gurobi.Var v[4,1] (value 1.0)>
    <gurobi.Var v[5,1] (value 1.0)>
    <gurobi.Var v[6,1] (value 1.0)>
    <gurobi.Var v[7,1] (value 1.0)>
    <gurobi.Var v[8,1] (value 1.0)>
    <gurobi.Var v[9,1] (value 1.0)>
    <gurobi.Var v[9,39] (value 1.0)>
    <gurobi.Var v[10,1] (value 1.0)>

    Best regards,

    Maliheh

    0
  • Benay Uslu
    First Question

    Thanks a lot, Maliheh, you really saved my day.

    Now that I handled the visiting period constraint, I got errors from the basic VRP constraints of my model. The main logic of my model is like this I have 10 nodes and a depot, 3 vehicles, and 40 days in my planning horizon. I try to construct a periodic vehicle routing schedule in which the visiting times of each node (regarding their frequency of visit) is satisfied, and a vehicle completes its tour (i.e, depart from depot (node 0), go to its assigned nodes for that period, and return to the depot at the end of the same period.). But each vehicle does not suppose to depart from depot at a given period, neither they should do that in every period. Everything should supposed to happen just when it is needed (and the visiting frequencies should determine this condition), and in the remaining times vehicles can remain idle in the depot, as this is what we essentially want when we try to minimize distance. similarly, not all the vehicles are supposed to be utilized if I can handle the tours with the same vehicle. I plan to solve this problem in a very large network consisting of 175 nodes, and I know when I do that I'm gonna need 3 vehicles or more. But since I cannot run my code in a simple pilot network, I cannot move forward. Here is my current code:

    # -*- coding: utf-8 -*-
    """
    Created on Thu Jan 6 21:58:28 2022

    @author: User
    """

    import random
    import pandas as pd
    import gurobipy as gp
    from gurobipy import*
    import matplotlib.pyplot as plt
    mdl = Model("periodic vrp")
    '''Sets and parameters'''

    N = []
    for i in range(11):
    N.append(i)

    I = []
    for i in range(1,11):
    I.append(i)

    I2 = [3,9]

    K = [1, 2, 3]

    T = []
    for i in range(1,41):
    T.append(i)

    T2 = []
    for i in range(1,21):
    T2.append(i)

    C=5

    random.seed(1)

    def distance(p1,p2):
    x_dist=(p1[0]-p2[0])
    y_dist=(p1[1]-p2[1])
    euclidian=(x_dist**2+y_dist**2)**0.5
    return euclidian

    D = [] # set of main distribution centers
    D_loc = {}
    for i in range(11):
    D.append(i)
    D_loc[i] = ((random.randint(0,100)),random.randint(0,100))

    DIST = {}
    for i in D:
    for j in D:
    DIST[(i,j)] = distance(D_loc[i],D_loc[j])

    q = {1:0.72,2:0.58,3:0.93,4:0.68,5:0.52,6:0.38,7:0.55,8:0.28,9:1.25,10:0.2}

    FREQ = {1:1,2:1,3:2,4:1,5:1,6:1,7:1,8:1,9:2,10:1}


    '''Decision Variables'''


    pairs= [(n,m,k,t) for n in N for m in N for k in K for t in T]
    x = mdl.addVars(pairs, vtype=GRB.BINARY, name ="x")

    bins=[(i,k,t) for i in I for k in K for t in T]
    y = mdl.addVars(bins, vtype=GRB.BINARY, name ="y")

    visits=[(i,t) for i in I for t in T]
    v = mdl.addVars(visits, vtype=GRB.BINARY, name ="v")

    mtzs=[(i,k,t) for i in I for k in K for t in T]
    u = mdl.addVars(mtzs, name ="u")

    '''Constraints'''

    #----------flow balance constraints and connection of decision variables to each other
    # for i in I:
    # for k in K:
    # for t in T:
    # mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == quicksum(x[i,n,k,t] for n in N if n !=i))

    for i in I:
    for k in K:
    for t in T:
    mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == y[i,k,t])
    for i in I:
    for k in K:
    for t in T:
    if t !=40:
    mdl.addConstr(quicksum(x[i,n,k,t+1] for n in N if n !=i) == y[i,k,t])


    #----------connection of decision variables to each other
    mdl.addConstrs(quicksum(y[i,k,t] for k in K) == v[i,t] for i in I for t in T)

    for n in I:
    for i in I:
    for t in T:
    if n !=i:
    mdl.addConstr(quicksum(x[n,i,k,t] for k in K) <= (v[n,t]+v[i,t])/2)

    for t in T:
    for j in I:
    mdl.addConstr(quicksum(x[n,i,k,t] for n in N for k in K) == v[i,t])
    #------------


    mdl.addConstrs(quicksum(y[i,k,t] for k in K for t in T) == FREQ[i] for i in I) #every node should be visited as much as its required frequency



    mdl.addConstrs(quicksum(x[0,i,k,40] for i in I) == 0 for k in K)
    mdl.addConstrs(quicksum(x[0,i,k,t] for i in I) <= 1 for t in T for k in K)
    #mdl.addConstrs(quicksum(x[i,0,k,t] for i in I) == 1 for t in T for k in K)

    mdl.addConstrs(u[j,k,t] - u[i,k,t] >= q[j]/FREQ[j]-C*(1-x[i,j,k,t]) for i in I for j in I for k in K for t in T if i!=j)
    mdl.addConstrs(q[i] <= u[i,k,t] for i in I for k in K for t in T)
    mdl.addConstrs(u[i,k,t] <= C for i in I for k in K for t in T)


    start_time, end_time = T[0], T[-1] + 1
    mdl.addConstrs(
    (v[i, t] == 1)
    >> (
    quicksum(
    v[i, j] for j in range(max(t + 1, start_time), min(t + 20, end_time))
    ) == 0
    )
    for i in I2
    for t in T
    )


    # mdl.addConstrs(v[i,t+20]+v[i,t+30] >= v[i,t]+v[i,t+10] for i in I4 for t in T4)

    '''Objective Function'''
    obj=quicksum(x[n,m,k,t]*DIST[(n,m)] for n in N for m in N for k in K for t in T)
    mdl.setObjective(obj, GRB.MINIMIZE)

    '''Solve'''

    mdl.optimize()

    print('Obj: %g' % mdl.objVal)


    for v in mdl.getVars():
    if v.x==1:
    print(v)

    and when I enforce the currently commented constraint of returning to the depot at the end of each period, which is:

    mdl.addConstrs(quicksum(x[i,0,k,t] for i in I) == 1 for t in T for k in K)

    i get infeasibility I guess because gurobi says that 

    AttributeError: Unable to retrieve attribute 'objVal'

    I really tried my best to work it out but could not find any way for a month. If you have any suggestions, I would be more than happy to hear that.

    PS: I already tried using gurobi's cbLazy constraints but it did not change anything for my case, or I could not implement it successfully, I am not sure to be honest. I directly retrieve the code 

    def subtourelim(model, where):
    if where == GRB.Callback.MIPSOL:
    # make a list of edges selected in the solution
    vals = model.cbGetSolution(model._vars)
    selected = gp.tuplelist((i, j) for i, j in model._vars.keys()
    if vals[i, j] > 0.5)
    # find the shortest cycle in the selected edge list
    tour = subtour(selected)
    if len(tour) < 40:
    # add subtour elimination constr. for every pair of cities in tour
    model.cbLazy(gp.quicksum(model._vars[i, j]
    for i, j in combinations(tour, 2))
    <= len(tour)-1)

    from the tsp example in here. It gives the same results with my subtour elimination constraints written with u_ijk variables.

    Thank you in advance for your help and cooperation.

    Best,

    0
  • Maliheh Aramon
    Gurobi Staff Gurobi Staff
    Hi Benay, 
     
    To debug your code and narrow down the infeasibility in your model, you can use the Gurobi Model.computeIIS() method to find an irreducible inconsistent subsystem (IIS) and write it into a file as below:
    mdl.computeIIS()
    mdl.write("model.ilp")
    To make things easier, it would be best to use the smallest problem size that results in infeasibility.
     
    Since your problem is an extension of the classical VRP problem, you might want to consider first implementing a classical VRP problem using one of the formulations available in the literature. You can then build on top of the base VRP, the different constraints that you have in your model.
     
    Best regards,
    Maliheh
    0

Please sign in to leave a comment.