mils_echelon: echelon formulation for the multi-item, multi-stage lot-sizing problem

    Parameters:
        - T: number of periods
        - K: set of resources
        - P: set of items
        - f[t,p]: set-up costs (on period t, for product p)
        - g[t,p]: set-up times
        - c[t,p]: variable costs
        - d[t,p]: demand values
        - h[t,p]: holding costs
        - a[t,k,p]: amount of resource k for producing p in period t
        - M[t,k]: resource k upper bound on period t
        - UB[t,p]: upper bound of production time of product p in period t
        - phi[(i,j)]: units of i required to produce a unit of j (j parent of i)

def mils_echelon(T,K,P,f,g,c,d,h,a,M,UB,phi):
    """
    mils_echelon: echelon formulation for the multi-item, multi-stage lot-sizing problem

    Parameters:
        - T: number of periods
        - K: set of resources
        - P: set of items
        - f[t,p]: set-up costs (on period t, for product p)
        - g[t,p]: set-up times
        - c[t,p]: variable costs
        - d[t,p]: demand values
        - h[t,p]: holding costs
        - a[t,k,p]: amount of resource k for producing p in period t
        - M[t,k]: resource k upper bound on period t
        - UB[t,p]: upper bound of production time of product p in period t
        - phi[(i,j)]: units of i required to produce a unit of j (j parent of i)
    """
    rho = calc_rho(phi) # rho[(i,j)]: units of i required to produce a unit of j (j ancestor of i)

    model = Model("multi-stage lotsizing -- echelon formulation")

    y,x,E,H = {},{},{},{}
    Ts = range(1,T+1)
    for p in P:
        for t in Ts:
            y[t,p] = model.addVar(vtype="B", name="y(%s,%s)"%(t,p))
            x[t,p] = model.addVar(vtype="C", name="x(%s,%s)"%(t,p))
            H[t,p] = h[t,p] - sum([h[t,q]*phi[q,p] for (q,p2) in phi if p2 == p])
            E[t,p] = model.addVar(vtype="C", name="E(%s,%s)"%(t,p))        # echelon inventory
        E[0,p] = model.addVar(vtype="C", name="E(%s,%s)"%(0,p))    # echelon inventory

    for t in Ts:
        for p in P:
            # flow conservation constraints
            dsum = d[t,p] + sum([rho[p,q]*d[t,q] for (p2,q) in rho if p2 == p])
            model.addCons(E[t-1,p] + x[t,p] == E[t,p] + dsum, "FlowCons(%s,%s)"%(t,p))

            # capacity connection constraints
            model.addCons(x[t,p] <= UB[t,p]*y[t,p], "ConstrUB(%s,%s)"%(t,p))

        # time capacity constraints
        for k in K:
            model.addCons(quicksum(a[t,k,p]*x[t,p] + g[t,p]*y[t,p] for p in P) <= M[t,k],
                            "TimeUB(%s,%s)"%(t,k))


    # calculate echelon quantities
    for p in P:
        model.addCons(E[0,p] == 0, "EchelonInit(%s)"%(p))
        for t in Ts:
            model.addCons(E[t,p] >= quicksum(phi[p,q]*E[t,q] for (p2,q) in phi if p2 == p),
                            "EchelonLB(%s,%s)"%(t,p))

    model.setObjective(\
        quicksum(f[t,p]*y[t,p] + c[t,p]*x[t,p] + H[t,p]*E[t,p] for t in Ts for p in P), \
        "minimize")

    model.data = y,x,E
    return model