code: def transp(I,J,c,d,M):
      
  
      model = Model("transportation")
  
      # Create variables
      x = {}
  
      for i in I:
          for j in J:
              x[i,j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j))
  
      # Demand constraints
      for i in I:
          model.addCons(quicksum(x[i,j] for j in J if (i,j) in x) == d[i], name="Demand(%s)" % i)
  
      # Capacity constraints
      for j in J:
          model.addCons(quicksum(x[i,j] for i in I if (i,j) in x) <= M[j], name="Capacity(%s)" % j)
  
      # Objective
      model.setObjective(quicksum(c[i,j]*x[i,j]  for (i,j) in x), "minimize")
  
      model.optimize()
  
      model.data = x
      return model
docs: transp -- model for solving the transportation problem
      Parameters:
          I - set of customers
          J - set of facilities
          c[i,j] - unit transportation cost on arc (i,j)
          d[i] - demand at node i
          M[j] - capacity
      Returns a model, ready to be solved.