code: def vrp(V, c, m, q, Q):
      
  
      model = Model("vrp")
      vrp_conshdlr = VRPconshdlr()
  
      x = {}
      for i in V:
          for j in V:
              if j > i and i == V[0]:       # depot
                  x[i,j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)"%(i,j))
              elif j > i:
                  x[i,j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)"%(i,j))
  
      model.addCons(quicksum(x[V[0],j] for j in V[1:]) == 2*m, "DegreeDepot")
      for i in V[1:]:
          model.addCons(quicksum(x[j,i] for j in V if j < i) +
                          quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)
  
      model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j>i), "minimize")
      model.data = x
  
      return model, vrp_conshdlr
docs: solve_vrp -- solve the vehicle routing problem.
         - start with assignment model (depot has a special status)
         - add cuts until all components of the graph are connected
      Parameters:
          - V: set/list of nodes in the graph
          - c[i,j]: cost for traversing edge (i,j)
          - m: number of vehicles available
          - q[i]: demand for customer i
          - Q: vehicle capacity
      Returns the optimum objective value and the list of edges used.