diff --git a/final_ls_cuts.py b/final_ls_cuts.py
new file mode 100644
index 000000000..cb58c1624
--- /dev/null
+++ b/final_ls_cuts.py
@@ -0,0 +1,173 @@
+##@file lotsizing_lazy.py
+# @brief solve the single-item lot-sizing problem.
+"""
+Approaches:
+    - sils: solve the problem using the standard formulation
+    - sils_cut: solve the problem using cutting planes
+
+Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
+"""
+from pyscipopt import Model, Conshdlr, quicksum, multidict, SCIP_RESULT, SCIP_PRESOLTIMING, SCIP_PROPTIMING
+
+
+class Conshdlr_sils(Conshdlr):
+
+    def addcut(self, checkonly, sol):
+        D, Ts = self.data
+        y, x, I = self.model.data
+        cutsadded = False
+
+        for ell in Ts:
+            lhs = 0
+            S, L = [], []
+            for t in range(1, ell + 1):
+                yt = self.model.getSolVal(sol, y[t])
+                xt = self.model.getSolVal(sol, 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]:
+                if checkonly:
+                    return True
+                else:
+                    # add cutting plane constraint
+                    self.model.addCons(quicksum([x[t] for t in L]) + \
+                                       quicksum(D[t, ell] * y[t] for t in S)
+                                       >= D[1, ell], removable=True)
+                cutsadded = True
+        return cutsadded
+
+    def conscheck(self, constraints, solution, checkintegrality, checklprows, printreason, completely):
+        if self.addcut(checkonly=True, sol=solution):
+            return {"result": SCIP_RESULT.INFEASIBLE}
+        else:
+            return {"result": SCIP_RESULT.FEASIBLE}
+
+    def consenfolp(self, constraints, nusefulconss, solinfeasible):
+        if self.addcut(checkonly=False, sol=None):
+            return {"result": SCIP_RESULT.CONSADDED}
+        else:
+            return {"result": SCIP_RESULT.FEASIBLE}
+
+    def conslock(self, constraint, locktype, nlockspos, nlocksneg):
+        pass
+
+
+def sils(T, f, c, d, h):
+    """sils -- LP lotsizing for the single item lot sizing problem
+    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 a model, ready to be solved.
+    """
+    model = Model("single item lotsizing")
+    Ts = range(1, T + 1)
+    M = sum(d[t] for t in Ts)
+    y, x, I = {}, {}, {}
+    for t in Ts:
+        y[t] = model.addVar(vtype="I", ub=1, name="y(%s)" % t)
+        x[t] = model.addVar(vtype="C", ub=M, name="x(%s)" % t)
+        I[t] = model.addVar(vtype="C", name="I(%s)" % t)
+    I[0] = 0
+
+    for t in Ts:
+        model.addCons(x[t] <= M * y[t], "ConstrUB(%s)" % t)
+        model.addCons(I[t - 1] + x[t] == I[t] + d[t], "FlowCons(%s)" % t)
+
+    model.setObjective( \
+        quicksum(f[t] * y[t] + c[t] * x[t] + h[t] * I[t] for t in Ts), \
+        "minimize")
+
+    model.data = y, x, I
+    return model
+
+
+def sils_cut(T, f, c, d, h, conshdlr):
+    """solve_sils -- solve the lot sizing problem with cutting planes
+       - start with a relaxed model
+       - used lazy constraints to elimitate fractional setup variables with cutting planes
+    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:
+        model.chgVarType(y[t], "C")
+    model.addVar(vtype="B", name="fake")  # for making the problem MIP
+
+    # 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
+
+    # include the lot sizing constraint handler
+    model.includeConshdlr(conshdlr, "SILS", "Constraint handler for single item lot sizing",
+                          sepapriority=0, enfopriority=-1, chckpriority=-1, sepafreq=-1, propfreq=-1,
+                          eagerfreq=-1, maxprerounds=0, delaysepa=False, delayprop=False, needscons=False,
+                          presoltiming=SCIP_PRESOLTIMING.FAST, proptiming=SCIP_PROPTIMING.BEFORELP)
+    conshdlr.data = D, Ts
+
+    model.data = y, x, I
+    return model
+
+
+def mk_example():
+    """mk_example: book example for the single item lot sizing"""
+    T = 5
+    _, f, c, d, h = multidict({
+        1: [3, 1, 5, 1],
+        2: [3, 1, 7, 1],
+        3: [3, 3, 3, 1],
+        4: [3, 3, 6, 1],
+        5: [3, 3, 4, 1],
+    })
+
+    return T, f, c, d, h
+
+
+if __name__ == "__main__":
+    T, f, c, d, h = mk_example()
+
+    model = sils(T, f, c, d, h)
+    y, x, I = model.data
+    model.optimize()
+    sils_obj = model.getObjVal()
+    print("\nOptimal value [standard]:", model.getObjVal())
+    print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t", "fix", "var", "h", "dem", "y", "x", "I"))
+    for t in range(1, T + 1):
+        print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % (
+        t, f[t], c[t], h[t], d[t], model.getVal(y[t]), model.getVal(x[t]), model.getVal(I[t])))
+
+    conshdlr = Conshdlr_sils()
+    model = sils_cut(T, f, c, d, h, conshdlr)
+    model.setBoolParam("misc/allowstrongdualreds", 0)
+    model.optimize()
+    sils_cut_obj = model.getObjVal()
+    y, x, I = model.data
+    print("\nOptimal value [cutting planes]:", model.getObjVal())
+    print("%8s%8s%8s%8s%8s%8s%12s%12s" % ("t", "fix", "var", "h", "dem", "y", "x", "I"))
+    for t in range(1, T + 1):
+        print("%8d%8d%8d%8d%8d%8.1f%12.1f%12.1f" % (
+        t, f[t], c[t], h[t], d[t], model.getVal(y[t]), model.getVal(x[t]), model.getVal(I[t])))
+
+    assert abs(sils_obj - sils_cut_obj) <= 1e-6