Skip to main content

VRP - Model and Code

Answered
Truc Vinh Do
Gurobi-versary
Conversationalist
First Question

Hi all,

I have a model, and I don't know how write code in Gurobi.

Could you please help me? 

Thank you very much.

 

 

0

Comments

12 comments

Date Votes
  • Mario Ruthmair
    Gurobi Staff Gurobi Staff

    Hi,

    Here are some resources that you could use to formulate your model:

    Do you have a more specific question regarding formulating this VRP in Python?

    Best regards,
    Mario

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    Yes Mario,

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    # -*- coding: utf-8 -*-
    """
    Created on Thu May 25 09:04:38 2023

    @author: dvtruc
    """
    import math
    import random
    import gurobipy as gp
    from gurobipy import *
    from gurobipy import GRB
    md1=Model("VRP")

    nbN=10
    nbK=4
    nbW=2
    demandMin=2
    demandMax=5

    random.seed(0)
    depot=[0,1]
    N=[*range(len(depot), len(depot)+ nbN)]
    node=depot +N
    K=[*range(1, nbK+1)]
    W=[* range(1, nbW+1)]
    arcs1=[(i,j) for i in N for j in N if i!=j]
    arcs2=[(i,j) for i in depot for j in N]
    #ei={i: random.randint(0,200) for i in N}
    ei=[]
    for i in range(nbN):
        ei.append(random.randint(0,200))

    #li={i: random.randint(500,800) for i in N}
    li=[]
    for i in range(nbN):
        li.append(random.randint(500,800))

    #qi={i: random.randint(demandMin,demandMax) for i in N}
    qi=[]
    for i in range(nbN):
        ei.append(random.randint(demandMin,demandMax))
    tiR=[]
    for i in range(nbN):
        tiR.append(random.randint(500,800))
    tiS=[]
    for i in range(nbN):
        tiS.append(random.randint(5,10))
        
    points = [(random.randint(0, 999), random.randint(0, 999)) for i in N]

    tijN = {(i, j): math.ceil(
            math.sqrt(sum((points[i-len(depot)][k] - points[j-len(depot)][k]) ** 2 for k in range(2)))
        )
        for i in N
        for j in N
        if i != j
    }

    tkD=[]
    for i in range(nbN):
        tkD.append(random.randint(5,10))

    points2 = [(random.randint(0, 999), random.randint(0, 999)) for i in depot]
    tki = {(i, j): math.ceil(
            math.sqrt(sum((points2[i][k] - points[j-len(depot)][k]) ** 2 for k in range(2)))
        )
        for i in depot
        for j in N
    }
    tH=800
    Q=[]
    for i in range(nbK):
        Q.append(random.randint(5,10))
    M=1000000

    pair_kw = [(k,w) for k in K for w in W]
    tkwL=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='tkwL')
    tkwB=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='tkwB')
    pair_kiw=[(k,i,w) for k in K for i in N for w in W]
    tkiw=md1.addVars(pair_kiw,vtype=GRB.INTEGER, name='tkiw')
    deltakw=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='deltakw')
    deltakiw=md1.addVars(pair_kiw,vtype=GRB.INTEGER, name="deltakiw")
    gamakwL=md1.addVars(pair_kw,vtype=GRB.BINARY, name="gamakwL")
    muykiw=md1.addVars(pair_kiw,vtype=GRB.BINARY, name="muykiw")
    nukiw=md1.addVars(pair_kiw,vtype=GRB.BINARY, name="nukiw")

    pair_ijkw=[(i,j,k,w) for i in N for j in N for k in K for w in W]
    epsilonijkw=md1.addVars(pair_ijkw,vtype=GRB.BINARY, name="epsilonijkw")

    md1.setObjective(quicksum(epsilonijkw[i,j,k,w] for i in N for j in N for k in K for w in W), GRB.MINIMIZE)
        

    I write the code above and I have some errors on:

    Line 9: import gurobipy as gp # unused
    Line 10:from gurobipy import * # unused

    Line 12:md1=Model("VRP") # Model undefined

    Line 90: md1.setObjective(quicksum(epsilonijkw[i,j,k,w] for i in N for j in N for k in K for w in W), GRB.MINIMIZE)

    quicksum undefined.

    Please help me

     

     

    0
  • Elisabeth Rodríguez Heck
    Gurobi Staff Gurobi Staff

    Hi,

    The Model() and quicksum() methods are part of the gurobipy package. Since you imported gurobipy as gp, you have to call these models using the gp prefix, for example, instead of writing

    md1=Model("VRP")

    you have to write

    md1=gp.Model("VRP")

    Similarly for quicksum.

    Best regards,

    Elisabeth

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    Thank you very much Dr.Elisabeth.

    For

    please help how to write the constraint.

    Tks.

     

    0
  • Mario Ruthmair
    Gurobi Staff Gurobi Staff

    Hi,

    Could you please post the code for the balance constraints that you already have and that does not work? Then we can see if we can spot the mistake.

    Best regards,
    Mario

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    md1.addConstrs(muykiw[k,l,w] +  gp.quicksum(epsilonijkw[i,l,k,w] for i in N) == nuykiw[k,l,w] + gp.quicksum(epsilonijkw[l,j,k,w] for j in N) for l in N for k in K for w in W )
    md1.addConstrs(nuykiw[k,l,w] + gp.quicksum(epsilonijkw[l,j,k,w] for j in N) == lamdakiw[k,i,w]  for l in N for k in K for w in W )

    Best regards,

    Trúc

    0
  • Mario Ruthmair
    Gurobi Staff Gurobi Staff

    It looks ok in terms of syntax, what is the error you get?

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    At first, I don't know how to write. While I am waiting for the help, I try to write.

    Thank you very much for your help.

    Best regards,

    Trúc

    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    Here is my code with "Model is infeasible".

    If you can, please help me.


    import math
    import random
    import gurobipy as gp
    #from gurobipy import *
    from gurobipy import GRB
    md1=gp.Model("VRP")

    nbN=10
    nbK=4
    nbW=2
    demandMin=2
    demandMax=5

    random.seed(0)
    depot=[0,1]
    N=[*range(len(depot), len(depot)+ nbN)]
    node=depot +N
    K=[*range(1, nbK+1)]
    W=[* range(1, nbW+1)]
    arcs1=[(i,j) for i in N for j in N if i!=j]
    arcs2=[(i,j) for i in depot for j in N]
    ei=[]
    for i in range(nbN):
        ei.append(random.randint(0,200))

    #li={i: random.randint(500,800) for i in N}
    li=[]
    for i in range(nbN):
        li.append(random.randint(500,800))

    qi=[]
    for i in range(nbN):
        qi.append(random.randint(demandMin,demandMax))
    tiR=[]
    for i in range(nbN):
        tiR.append(random.randint(500,800))
    tiS=[]
    for i in range(nbN):
        tiS.append(random.randint(5,10))
        
    points = [(random.randint(0, 999), random.randint(0, 999)) for i in N]

    tijN = {(i, j): math.ceil(
            math.sqrt(sum((points[i-len(depot)][k] - points[j-len(depot)][k]) ** 2 for k in range(2)))
        )
        for i in N
        for j in N
        if i != j
    }

    tkD=[]
    for i in range(nbN):
        tkD.append(random.randint(5,10))

    points2 = [(random.randint(0, 999), random.randint(0, 999)) for i in depot]
    tki = {(i, j): math.ceil(
            math.sqrt(sum((points2[i][k] - points[j-len(depot)][k]) ** 2 for k in range(2)))
        )
        for i in depot
        for j in N
    }
    VK=[]
    VK={i:random.randint(0, len(depot)-1) for i in K}
    tki2 = {(i, j): math.ceil(
            math.sqrt(sum((points2[VK[i]][k] - points[j-len(K)][k]) ** 2 for k in range(2)))
        )
        for i in K
        for j in N
    }
    tH=800
    Q=[]
    for i in range(nbK):
        Q.append(random.randint(5,10))
    M=1000000

    pair_kw = [(k,w) for k in K for w in W]
    tkwL=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='tkwL')
    tkwB=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='tkwB')
    pair_kiw=[(k,i,w) for k in K for i in N for w in W]
    tkiw=md1.addVars(pair_kiw,vtype=GRB.INTEGER, name='tkiw')
    deltakw=md1.addVars(pair_kw,vtype=GRB.INTEGER, name='deltakw')
    deltakiw=md1.addVars(pair_kiw,vtype=GRB.INTEGER, name="deltakiw")
    gamakwL=md1.addVars(pair_kw,vtype=GRB.BINARY, name="gamakwL")
    muykiw=md1.addVars(pair_kiw,vtype=GRB.BINARY, name="muykiw")
    nuykiw=md1.addVars(pair_kiw,vtype=GRB.BINARY, name="nuykiw")

    pair_ijkw=[(i,j,k,w) for i in N for j in N for k in K for w in W]
    epsilonijkw=md1.addVars(pair_ijkw,vtype=GRB.BINARY, name="epsilonijkw")
    lamdakiw=md1.addVars(pair_kiw,vtype=GRB.BINARY, name="lamdakiw")

    md1.setObjective(gp.quicksum(tijN[i,j]*epsilonijkw[i,j,k,w] for i in N for j in N for k in K for w in W if i!=j)\
        +gp.quicksum(tki2[k,i]*muykiw[k,i,w] for k in K for i in N for w in W) \
        +gp.quicksum(tki2[k,i]*nuykiw[k,i,w] for k in K for i in N for w in W) \
        , GRB.MINIMIZE)
    md1.addConstrs(gp.quicksum(lamdakiw[k,i,w] for k in K for w in W) ==1 for i in N)
    md1.addConstrs(muykiw[k,l,w] +  gp.quicksum(epsilonijkw[i,l,k,w] for i in N) == nuykiw[k,l,w] + gp.quicksum(epsilonijkw[l,j,k,w] for j in N) for l in N for k in K for w in W )
    md1.addConstrs(nuykiw[k,l,w] + gp.quicksum(epsilonijkw[l,j,k,w] for j in N) == lamdakiw[k,i,w]  for l in N for k in K for w in W )
    md1.addConstrs(gp.quicksum(muykiw[k,i,w] for i in N) == gp.quicksum(nuykiw[k,i,w] for i in N) for k in K for w in W)
    md1.addConstrs(gp.quicksum(nuykiw[k,i,w] for i in N) == gamakwL[k,w] for k in K for w in W)
    md1.addConstrs(gamakwL[k,w+1]<=gamakwL[k,w] for k in K for w in W if w<=nbW-1)
    md1.addConstrs(gp.quicksum(lamdakiw[k,i,w] for i in N) <= gamakwL[k,w]*nbN for k in K for w in W)
    md1.addConstrs(deltakw[k,w]==gp.quicksum(lamdakiw[k,i,w]*qi[i-len(depot)] for i in N) for k in K for w in W)
    md1.addConstrs(deltakw[k,w] <= gamakwL[k,w]*Q[k-1] for k in K for w in W)
    md1.addConstrs(tkiw[k,i,w] <= lamdakiw[k,i,w]*li[i-len(depot)] for i in N for k in K for w in W)
    md1.addConstrs(tkwL[k,w] >= lamdakiw[k,i,w]*(tiR[i-len(depot)]+tkD[k]) - M*(1-gamakwL[k,w]) for i in N for k in K for w  in W)
    md1.addConstrs(tkwL[k,w+1] >= tkwB[k,w] + tkD[k] - M*(1-gamakwL[k,w+1]) for k in K for w in W if w<=nbW-1)
    md1.addConstrs(tkiw[k,i,w] >= tkwL[k,w] + tki2[k,i] - M*(1 - muykiw[k,i,w]) for i in N for k in K for w in W)
    md1.addConstrs(tkwB[k,w] >= tkiw[k,i,w] + tiS[i-len(depot)] + tki2[k,i] - M*(1 - nuykiw[k,i,w]) for i in N for k in K for w in W)
    md1.addConstrs(tkwB[k,w] >= ei[i-len(depot)] + tiS[i-len(depot)] + tki2[k,i] - M*(1 - nuykiw[k,i,w]) for i in N for k in K for w in W)
    md1.addConstrs(tkwB[k,w] >= tkiw[k,i,w] + tiS[i-len(depot)] + tijN[i,j] - M*(1 - epsilonijkw[i,j,k,w]) for i in N for j in N for k in K for w in W if i!=j)
    md1.addConstrs(tkwB[k,w] >= ei[i-len(depot)] + tiS[i-len(depot)] + tijN[i,j] - M*(1 - epsilonijkw[i,j,k,w]) for i in N for j in N for k in K for w in W if i!=j)
    md1.addConstrs(tkwB[k,w]<=tH for k in K for w in W )
    md1.optimize()

     

     

    0
  • Mario Ruthmair
    Gurobi Staff Gurobi Staff

    Here is an article on how to deal with infeasible models.
    In general, there could be several reasons why your model is infeasible:

    • The VRP instance itself could be infeasible.
    • There are issues in your constraints but this is not easy to debug just by looking at the code. You could iteratively remove sets of constraints and see if it gets feasible. Then you might have an idea where the bug is located.
    0
  • Truc Vinh Do
    Gurobi-versary
    Conversationalist
    First Question

    Thank you very much.

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post