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.

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

    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