code: def sils_cut(T,f,c,d,h):
      
      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
docs: 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.