solve_sils -- solve the lot sizing problem with cutting planes
       - start with a relaxed model
       - add cuts until there are no fractional setup variables
    Parameters:
        - T: number of periods
        - P: set of products
        - f[t]: set-up costs (on period t)
        - c[t]: variable costs
        - d[t]: demand values
        - h[t]: holding costs
    Returns the final model solved, with all necessary cuts added.

def sils_cut(T,f,c,d,h):
    """solve_sils -- solve the lot sizing problem with cutting planes
       - start with a relaxed model
       - add cuts until there are no fractional setup variables
    Parameters:
        - T: number of periods
        - P: set of products
        - f[t]: set-up costs (on period t)
        - c[t]: variable costs
        - d[t]: demand values
        - h[t]: holding costs
    Returns the final model solved, with all necessary cuts added.
    """
    Ts = range(1,T+1)

    model = sils(T,f,c,d,h)
    y,x,I = model.data

    # relax integer variables
    for t in Ts:
        y[t].vtype = "C"

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

    EPS = 1.e-6
    cuts = True
    while cuts:
        model.optimize()
        cuts = False
        for ell in Ts:
            lhs = 0
            S,L = [],[]
            for t in range(1,ell+1):
                yt = model.getVal(y[t])
                xt = model.getVal(x[t])
                if D[t,ell]*yt < xt:
                    S.append(t)
                    lhs += D[t,ell]*yt
                else:
                    L.append(t)
                    lhs += xt
            if lhs < D[1,ell]:
                # add cutting plane constraint
                model.addCons(quicksum([x[t] for t in L]) +\
                                quicksum(D[t,ell] * y[t] for t in S)
                                >= D[1,ell])
                cuts = True

    model.data = y,x,I
    return model