Integer Programming with Simple Branch and Bound - Optimization for Decision Science

\begin{align*} \min \quad & Z = x + y_1 + 3 y_2 + 2 y_3 \\ \mathrm{s.t.} \quad & -x + 3 y_1 + 2 y_2 + y_3 \leq 0 \\ & -5 y_1 - 8 y_2 - 3 y_3 \leq -9 \\ & x \geq 0, y \in \{0,1\}^3 \end{align*}

(1)

import pyomo.environ as pyo
from pyomo.opt import SolverStatus, TerminationCondition


def create_relaxation():

    ## Create Pyomo Model
    m = pyo.ConcreteModel()

    ## Declare variables (all continuous for now)
    m.x = pyo.Var(domain=pyo.NonNegativeReals)
    m.y1 = pyo.Var(bounds=(0, 1))
    m.y2 = pyo.Var(bounds=(0, 1))
    m.y3 = pyo.Var(bounds=(0, 1))

    ## Declare constraints
    def con1_rule(m):
        return -1.0 * m.x + 3.0 * m.y1 + 2.0 * m.y2 + m.y3 <= 0

    m.con1 = pyo.Constraint(rule=con1_rule)

    def con2_rule(m):
        return -5.0 * m.y1 - 8.0 * m.y2 - 3.0 * m.y3 <= -9

    m.con2 = pyo.Constraint(rule=con2_rule)

    ## Declare objective
    def obj_rule(m):
        return m.x + m.y1 + 3 * m.y2 + 2 * m.y3

    m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)

    return m


def solve_print(m, verbose=False):

    # Solve the model.
    #
    # `load_solutions=False` matters here: several branch-and-bound nodes below
    # are infeasible, and appsi_highs raises rather than returning a results
    # object if you ask it to load a solution that does not exist. Deferring the
    # load lets us inspect the termination condition first, then load only when
    # there is something to load.
    opt = pyo.SolverFactory("appsi_highs")
    results = opt.solve(m, tee=verbose, load_solutions=False)

    if (results.solver.status == SolverStatus.ok) and (
        results.solver.termination_condition == TerminationCondition.optimal
    ):
        # Print solution when the solution is optimal
        m.solutions.load_from(results)
        print("x = ", pyo.value(m.x))
        print("y1 = ", pyo.value(m.y1))
        print("y2 = ", pyo.value(m.y2))
        print("y3 = ", pyo.value(m.y3))
        print("Z = ", pyo.value(m.obj))

    elif results.solver.termination_condition == TerminationCondition.infeasible:
        print("Infeasible.")

    else:
        # Something else is wrong
        print("Solver Status:", results.solver.status)

    # m.pprint()

Node 1 (Root)¶

Full LP relaxation.

m = create_relaxation()
solve_print(m, False)

x =  2.6
y1 =  0.2
y2 =  1.0
y3 =  0.0
Z =  5.8

Node 2¶

m = create_relaxation()
m.y1.fix(0.0)
solve_print(m, False)

x =  2.3333333333333335
y1 =  0.0
y2 =  1.0
y3 =  0.3333333333333333
Z =  6.0

Node 3¶

m = create_relaxation()
m.y1.fix(1.0)
solve_print(m, False)

x =  4.0
y1 =  1.0
y2 =  0.5
y3 =  0.0
Z =  6.5

Node 4¶

m = create_relaxation()
m.y1.fix(0.0)
m.y3.fix(0.0)
solve_print(m, False)

Infeasible.

Node 5¶

m = create_relaxation()
m.y1.fix(0.0)
m.y3.fix(1.0)
solve_print(m, False)

x =  2.5
y1 =  0.0
y2 =  0.75
y3 =  1.0
Z =  6.75

Node 6¶

m = create_relaxation()
m.y1.fix(1.0)
m.y2.fix(0.0)
solve_print(m, False)

Infeasible.

Node 7¶

m = create_relaxation()
m.y1.fix(1.0)
m.y2.fix(1.0)
solve_print(m, False)

x =  5.0
y1 =  1.0
y2 =  1.0
y3 =  0.0
Z =  9.0

Node 8¶

m = create_relaxation()
m.y1.fix(0.0)
m.y3.fix(1.0)
m.y2.fix(0.0)
solve_print(m, False)

Infeasible.

Node 9¶

m = create_relaxation()
m.y1.fix(0.0)
m.y3.fix(1.0)
m.y2.fix(1.0)
solve_print(m, False)

x =  3.0
y1 =  0.0
y2 =  1.0
y3 =  1.0
Z =  8.0