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.

def mctransp(I,J,K,c,d,M):
    """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.
    """

    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