mils: standard formulation for the multi-item lot-sizing problem
    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(T,P,f,g,c,d,h,M):
    """
    mils: standard formulation for the multi-item lot-sizing problem
    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_callback(model,where):
        # remember to set     model.params.DualReductions = 0     before using!
        if where != GRB.Callback.MIPSOL and where != GRB.Callback.MIPNODE:
            return
        for p in P:
            for ell in Ts:
                lhs = 0
                S,L = [],[]
                for t in range(1,ell+1):
                    yt = model.cbGetSolution(y[t,p])
                    xt = model.cbGetSolution(x[t,p])
                    if D[t,ell,p]*yt < xt:
                        S.append(t)
                        lhs += D[t,ell,p]*yt
                    else:
                        L.append(t)
                        lhs += xt
                if lhs < D[1,ell,p]:
                    # add cutting plane constraint
                    model.cbLazy(quicksum(x[t,p] for t in L) +\
                                 quicksum(D[t,ell,p] * y[t,p] for t in S)
                                 >= D[1,ell,p])
        return

    model = Model("standard multi-item lotsizing")

    y,x,I = {},{},{}
    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))
            I[t,p] = model.addVar(vtype="C", name="I(%s,%s)"%(t,p))
        I[0,p] = 0

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

        for p in P:
            # flow conservation constraints
            model.addCons(I[t-1,p] + x[t,p] == I[t,p] + d[t,p], "FlowCons(%s,%s)"%(t,p))

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

            # tighten constraints
            model.addCons(x[t,p] <= d[t,p]*y[t,p] + I[t,p], "Tighten(%s,%s)"%(t,p))

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

    # compute D[i,j,p] = sum_{t=i}^j d[t,p]
    D = {}
    for p in P:
        for t in Ts:
            s = 0
            for j in range(t,T+1):
                s += d[j,p]
                D[t,j,p] = s

    model.data = y,x,I
    return model,mils_callback