code: def mctransp(I,J,K,c,d,M):
      
  
      model = Model("multi-commodity transportation")
  
      # Create variables
      x = {}
      for (i,j,k) in c:
          x[i,j,k] = model.addVar(vtype="C", name="x(%s,%s,%s)" % (i,j,k), obj=c[i,j,k])
  
      # tuplelist is a Gurobi data structure to manage lists of equal sized tuples - try itertools as alternative
      arcs = tuplelist([(i,j,k) for (i,j,k) in x])
  
      # Demand constraints
      for i in I:
          for k in K:
              model.addCons(sum(x[i,j,k] for (i,j,k) in arcs.select(i,"*",k)) == d[i,k], "Demand(%s,%s)" % (i,k))
  
      # Capacity constraints
      for j in J:
          model.addCons(sum(x[i,j,k] for (i,j,k) in arcs.select("*",j,"*")) <= M[j], "Capacity(%s)" % j)
  
      model.data = x
      return model
docs: mctransp -- model for solving the Multi-commodity Transportation Problem
      Parameters:
          - I: set of customers
          - J: set of facilities
          - K: set of commodities
          - c[i,j,k]: unit transportation cost on arc (i,j) for commodity k
          - d[i][k]: demand for commodity k at node i
          - M[j]: capacity
      Returns a model, ready to be solved.