3.3 Integer cuts:

Often, it can be important to find not only the “best” solution, but a number of solutions that are equally optimal, or close to optimal. For discrete optimization problems, this can be done using something known as an integer cut. Consider again the knapsack problem where the choice of which items to select is a discrete variable \(x_{i} \forall i \in A\). Let \(x_{i}^{*}\) be a particular set of x values we want to remove from the feasible solution space. We define an integer cut using two sets. The first set \(S_{0}\) contains the indices for those variables whose current solution is 0, and the second set \(S_{1}\) consists of indices for those variables whose current solution is 1. Given these two sets, an integer cut constraint that would prevent such a solution from appearing again is defined by,

\[\sum _{i \in S_{0}}x[i] + \sum _{i \in \in S_{1}}(1 - x[i]) \geq 1\]

Below, we write a loop that solves the problem 5 times, adding an integer cut to remove the previous solution, and printing the value of the objective function and the solution at each iteration of the loop.

import pyomo.environ as pyo

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

def obj_rule(m):
    return sum( b[i]*m.x[i] for i in A )
model.obj = pyo.Objective(rule=obj_rule, sense = pyo.maximize )

def weight_con_rule(m):
    return sum( w[i]*m.x[i] for i in A ) <= W_max
model.weight_con = pyo.Constraint(rule=weight_con_rule)

opt = pyo.SolverFactory('glpk')


# create the ConstraintList to hold the integer cuts
model.int_cuts = pyo.ConstraintList()

# loop 5 times
for l in range(5):
    # solve the problem
    result_obj = opt.solve(model)

    # print the solution
    output_str = 'Obj: ' + str(pyo.value(model.obj))
    for i in A:
        output_str += "  x[%s]: %f" % (str(i), pyo.value(model.x[i]))
    print(output_str)           
    
    # add the integer cut based on the current solution
    cut_expr = 0
    for i in A:
        if pyo.value(model.x[i]) < 0.5:
            cut_expr += model.x[i]
        else:
            cut_expr += (1.0 - model.x[i])
    model.int_cuts.add(cut_expr >= 1)

Obj: 25.0  x[hammer]: 1.000000  x[wrench]: 0.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 20.0  x[hammer]: 0.000000  x[wrench]: 1.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 19.0  x[hammer]: 1.000000  x[wrench]: 0.000000  x[screwdriver]: 0.000000  x[towel]: 1.000000

Obj: 17.0  x[hammer]: 0.000000  x[wrench]: 0.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 14.0  x[hammer]: 0.000000  x[wrench]: 1.000000  x[screwdriver]: 0.000000  x[towel]: 1.000000