mils_fl: facility location formulation for the multi-item lot-sizing problem

    Requires more variables, but gives a better solution because LB is
    better than the standard formulation.  It can be used as a
    heuristic method that is sometimes better than relax-and-fix.

    Parameters:
        - T: number of periods
        - P: set of products
        - 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
        - M[t]:   resource upper bound on period t
    Returns a model, ready to be solved.

def mils_fl(T,P,f,g,c,d,h,M):
    """
    mils_fl: facility location formulation for the multi-item lot-sizing problem

    Requires more variables, but gives a better solution because LB is
    better than the standard formulation.  It can be used as a
    heuristic method that is sometimes better than relax-and-fix.

    Parameters:
        - T: number of periods
        - P: set of products
        - 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
        - M[t]:   resource upper bound on period t
    Returns a model, ready to be solved.
    """
    Ts = range(1,T+1)

    model = Model("multi-item lotsizing -- facility location formulation")

    y,X = {},{}
    for p in P:
        for t in Ts:
            y[t,p] = model.addVar(vtype="B", name="y(%s,%s)"%(t,p))
            for s in range(1,t+1):
                X[s,t,p] = model.addVar(name="X(%s,%s,%s)"%(s,t,p))


    for t in Ts:
        # capacity constraints
        model.addCons(quicksum(X[t,s,p] for s in range(t,T+1) for p in P) + \
                        quicksum(g[t,p]*y[t,p] for p in P) <= M[t],
                        "Capacity(%s)"%(t))

        for p in P:
            # demand satisfaction constraints
            model.addCons(quicksum(X[s,t,p] for s in range(1,t+1)) == d[t,p], "Demand(%s,%s)"%(t,p))

            # connection constraints
            for s in range(1,t+1):
                model.addCons(X[s,t,p] <= d[t,p] * y[s,p], "Connect(%s,%s,%s)"%(s,t,p))

    C = {} # variable costs plus holding costs
    for p in P:
        for s in Ts:
            sumC = 0
            for t in range(s,T+1):
                C[s,t,p] = (c[s,p] + sumC)
                sumC += h[t,p]

    model.setObjective(quicksum(f[t,p]*y[t,p] for t in Ts for p in P) + \
                       quicksum(C[s,t,p]*X[s,t,p] for t in Ts for p in P for s in range(1,t+1)),
                       "minimize")


    model.data = y,X
    model.write("tmp.lp")
    return model