6. Solving Mixed Integer Programs with Branch & Cut#

6.1. Mixed Integer Programs (MIPs)#

Mixed integer programs (MIPs) are those where some of the decision variables are discrete while others are continuous. These problems may come in a number of varieties, such as mixed integer quadratically constrainted programs (MIQCP), which contain at least one constraint with a quadratic term. Another variety is mixed integer quadratic programs (MIQPs), whose constraints are linear but the objective function is quadratic.

Here, we will focus on mixed integer linear programs (MILPs), where the objective and all the constraints are linear. In MILPs (as well as in other MIPs), the discreteness of the problem is often imposed by modeling a choice among a fixed number of options (sometimes as simple as a “yes or no” decision) or having to work with a whole number of resource. In comparison to purely linear programs of comparable sizes, MIPs are often more difficult to solve. This is due to the fact that in going from continuous variables to integer ones, the geometry of the feasible region gets more complicated.

6.2. Branch and cut#

How can we solve these problems? If all variables were continuous, this would reduce to a linear program, which we know how to solve with the simplex method (or with other methods like the interior point method). Doing so is called a relaxation, i.e. approximating a difficult problem by one that is easier to solve. The relaxed LP solution gives an upper bound on the true optimum (assuming one exists), since the feasible set is larger. If the LP solution happens to have integers for the correct variables, then we got lucky and ended up finding the optimum.

Suppose however that the relaxed solution contains non-integer values for integer variables. What then? One idea is to add a cuts, i.e. linear constraints that exclude the relaxed LP solution but don’t exclude any integer points. We then end up with more LPs but with stricker relaxations, which we can continue solving and cutting similarly.

The branching part of this algorithm comes from the fact that as we add cuts and end up with more problems, some of them can be eliminated from consideration without having to solve them. For example, imagine that we solved the LP relaxation of an MILP and the relaxed solution contains \(x_1=1.5\) where \(x_1\) is supposed to be an integer. Assuming that these won’t eliminate any integer solutions, we may create two LPs by adding constraints \(x_1\leq 1\) and \(x_1\geq 2\). In both these LPs, the original relaxed solution is infeasible. If the former has an integer \(x_1\) with a better objective function value when the latter has a non-integer \(x_1\) with a worse objective value, then one can prune the generation of further LPs from the latter. This is because any further LP will add additional constaints, therefore limiting the feasible region which will make the objective function even worse. Since the former already yielded an integer solution with a better objective value, any further LPs will fail to outperform it.

6.3. Example#

6.3.1. Problem#

This example is from [Wil13] Chapter 12.15.

A number of power stations are committed to meeting the following electricity load demands over a day:

12 p.m. to 6 a.m.	15 000 MW
6 a.m. to 9 a.m.	30 000 MW
9 a.m. to 3 p.m.	25 000 MW
3 p.m. to 6 p.m.	40 000 MW
6 p.m. to 12 p.m.	27 000 MW

There are three types of generating unit available: 12 of type 1, 10 of type 2 and ﬁve of type 3. Each generator has to work between a minimum and a maximum level. There is an hourly cost of running each generator at minimum level. In addition, there is an extra hourly cost for each megawatt at which a unit is operated above the minimum level. Starting up a generator also involves a cost. All this information is given in the table below (with costs in £).

In addition to meeting the estimated load demands there must be sufﬁcient generators working at any time to make it possible to meet an increase in load of up to 15%. This increase would have to be accomplished by adjusting the output of generators already operating within their permitted limits.

	Minimum level	Maximum level	Cost per hour at minimum	Cost per hour per megawatt above minimum	Cost
Type 1	850 MW	2000 MW	1000	2	2000
Type 2	1250 MW	1750 MW	2600	1.30	1000
Type 3	1500 MW	4000 MW	3000	3	500

Which generators should be working in which periods of the day to minimise total cost?

6.3.2. Modeling#

This model has a total of 55 constraints and 30 simple upper bounds. There are 45 variables, of which 30 are general integer variables.

6.3.2.1. Variables#

\[ \begin{align}\begin{aligned}n_{ij} = \text{number of generating units of type } i \text{ working in period } j\\s_{ij} = \text{number of generators of type } i \text{ started up in period } j\\x_{ij} = \text{total output rate from generators of type } i \text{ in period } j\end{aligned}\end{align} \]

\(x_{ij}\) are continuous, the rest are integer variables.

6.3.2.2. Constraints#

Demand must be met in each period:

\[\sum_ i x_{ij} \geq D_j \text{ for all j},\]

where \(D_j\) is demand given in period \(j\).

Output must lie within the limits of the generators working:

\[ \begin{align}\begin{aligned}x_{ij}\geq m_in_{ij} \text{ for all} i \text{ and }j\\x_{ij}\leq M_in_{ij} \text{ for all} i \text{ and }j\end{aligned}\end{align} \]

where \(m_i\) and \(M_i\) are the given minimum and maximum output levels for generators of type \(i\).

The extra guaranteed load requirement must be able to be met witout starting up any more generators:

\[\sum_iM_in_{ij}\geq \frac{115}{100}D_j \text{ for all }j\]

The number of generators started in period \(j\) must equal the increase in number

\[s_{ij}\geq n_{ij}-n_{ij-1} \text{for all } i \text{ and } j\]

where \(n_{ij}\) is the number of generators started in period \(j\) (when \(j=1\), period \(j-1\) is taken as 5).

In addition, all the integer variables have simple upper bounds corresponding to the total number of generators of each type.

6.3.2.3. Objective function (to be minimized)#

\[\text{Cost} = \sum_{i,j}C_{ij}(x_{ij}-m_in_{ij})+\sum_{i,j}E_{ij}n_{ij}+\sum_{i,j}F_is_{ij}\]

where \(C_{ij}\) are costs per hour per megawatt above minimum level multipled by the number of hours in the period; \(E_{ij}\) are costs per hour for operating at minimum level multipled by the number of hours in the period and \(F_i\) are start-up costs.

6.3.3. Solution#

First we import packages and define the relevant constants

using JuMP
import HiGHS

model = Model(HiGHS.Optimizer)
n_types = [12, 10, 5]
min_level = [850, 1250, 1500]
max_level = [2000, 1750, 4000]
cost_per_hour_at_min = [1000, 2600, 3000]
cost_per_hour_per_mw = [2, 1.30, 3]
cost_per_startup = [2000, 1000, 500]
demand = [15000, 30000, 25000, 40000, 27000]
hours_per_period = [6, 3, 6, 3, 6]

Create the model.

@variable(model, 0 <= n[i = 1:3, j = 1:5] <= n_types[i], Int)
@variable(model, 0 <= s[i = 1:3, j=1:5] <= n_types[i], Int)
@variable(model, 0 <= x[i = 1:3, j=1:5])

# Demand is met
@constraint(model, demand_met[j in 1:5], sum(x[i,j] for i in 1:3) >= demand[j])
# Power generated is within the limits
@constraint(model, output_lb[i in 1:3, j in 1:5], x[i,j] >= min_level[i]*n[i,j])
@constraint(model, output_ub[i in 1:3, j in 1:5], x[i,j] <= max_level[i]*n[i,j])
# increased load margin is satisfied
@constraint(model, extra_load[j in 1:5], sum(max_level[i]*n[i,j] for i in 1:3) >= 1.15*demand[j] )
# The number of started generators match
@constraint(model, continuity[i in 1:3, j in 2:5], s[i,j] >= n[i,j] - n[i,j-1])
@constraint(model, continuity2[i in 1:3], s[i,1] >= n[i,1] - n[i,5])

@objective(model, Min, 
    sum(cost_per_hour_per_mw[i]*hours_per_period[j]*(x[i,j]-min_level[i]*n[i,j]) for i in 1:3, j in 1:5) 
    + sum(cost_per_hour_at_min[i]*hours_per_period[j]*n[i,j] for i in 1:3, j in 1:5) 
    + sum(cost_per_startup[i]*s[i,j] for i in 1:3, j in 1:5))

print(model)

\[\begin{split} \begin{aligned} \min\quad & 12 x_{1,1} - 4200 n_{1,1} + 6 x_{1,2} - 2100 n_{1,2} + 12 x_{1,3} - 4200 n_{1,3} + 6 x_{1,4} - 2100 n_{1,4} + 12 x_{1,5} - 4200 n_{1,5} + 7.800000000000001 x_{2,1} + 5850 n_{2,1} + 3.9000000000000004 x_{2,2} + 2925 n_{2,2} + 7.800000000000001 x_{2,3} + 5850 n_{2,3} + 3.9000000000000004 x_{2,4} + 2925 n_{2,4} + 7.800000000000001 x_{2,5} + 5850 n_{2,5} + 18 x_{3,1} - 9000 n_{3,1} + 9 x_{3,2} - 4500 n_{3,2} + 18 x_{3,3} - 9000 n_{3,3} + 9 x_{3,4} - 4500 n_{3,4} + 18 x_{3,5} - 9000 n_{3,5} + 2000 s_{1,1} + 2000 s_{1,2} + 2000 s_{1,3} + 2000 s_{1,4} + 2000 s_{1,5} + 1000 s_{2,1} + 1000 s_{2,2} + 1000 s_{2,3} + 1000 s_{2,4} + 1000 s_{2,5} + 500 s_{3,1} + 500 s_{3,2} + 500 s_{3,3} + 500 s_{3,4} + 500 s_{3,5}\\ \text{Subject to} \quad & x_{1,1} + x_{2,1} + x_{3,1} \geq 15000\\ & x_{1,2} + x_{2,2} + x_{3,2} \geq 30000\\ & x_{1,3} + x_{2,3} + x_{3,3} \geq 25000\\ & x_{1,4} + x_{2,4} + x_{3,4} \geq 40000\\ & x_{1,5} + x_{2,5} + x_{3,5} \geq 27000\\ & -850 n_{1,1} + x_{1,1} \geq 0\\ & -1250 n_{2,1} + x_{2,1} \geq 0\\ & -1500 n_{3,1} + x_{3,1} \geq 0\\ & -850 n_{1,2} + x_{1,2} \geq 0\\ & -1250 n_{2,2} + x_{2,2} \geq 0\\ & -1500 n_{3,2} + x_{3,2} \geq 0\\ & -850 n_{1,3} + x_{1,3} \geq 0\\ & -1250 n_{2,3} + x_{2,3} \geq 0\\ & -1500 n_{3,3} + x_{3,3} \geq 0\\ & -850 n_{1,4} + x_{1,4} \geq 0\\ & -1250 n_{2,4} + x_{2,4} \geq 0\\ & -1500 n_{3,4} + x_{3,4} \geq 0\\ & -850 n_{1,5} + x_{1,5} \geq 0\\ & -1250 n_{2,5} + x_{2,5} \geq 0\\ & -1500 n_{3,5} + x_{3,5} \geq 0\\ & 2000 n_{1,1} + 1750 n_{2,1} + 4000 n_{3,1} \geq 17250\\ & 2000 n_{1,2} + 1750 n_{2,2} + 4000 n_{3,2} \geq 34500\\ & 2000 n_{1,3} + 1750 n_{2,3} + 4000 n_{3,3} \geq 28749.999999999996\\ & 2000 n_{1,4} + 1750 n_{2,4} + 4000 n_{3,4} \geq 46000\\ & 2000 n_{1,5} + 1750 n_{2,5} + 4000 n_{3,5} \geq 31049.999999999996\\ & n_{1,1} - n_{1,2} + s_{1,2} \geq 0\\ & n_{2,1} - n_{2,2} + s_{2,2} \geq 0\\ & n_{3,1} - n_{3,2} + s_{3,2} \geq 0\\ & n_{1,2} - n_{1,3} + s_{1,3} \geq 0\\ & n_{2,2} - n_{2,3} + s_{2,3} \geq 0\\ & n_{3,2} - n_{3,3} + s_{3,3} \geq 0\\ & n_{1,3} - n_{1,4} + s_{1,4} \geq 0\\ & n_{2,3} - n_{2,4} + s_{2,4} \geq 0\\ & n_{3,3} - n_{3,4} + s_{3,4} \geq 0\\ & n_{1,4} - n_{1,5} + s_{1,5} \geq 0\\ & n_{2,4} - n_{2,5} + s_{2,5} \geq 0\\ & n_{3,4} - n_{3,5} + s_{3,5} \geq 0\\ & -n_{1,1} + n_{1,5} + s_{1,1} \geq 0\\ & -n_{2,1} + n_{2,5} + s_{2,1} \geq 0\\ & -n_{3,1} + n_{3,5} + s_{3,1} \geq 0\\ & -2000 n_{1,1} + x_{1,1} \leq 0\\ & -1750 n_{2,1} + x_{2,1} \leq 0\\ & -4000 n_{3,1} + x_{3,1} \leq 0\\ & -2000 n_{1,2} + x_{1,2} \leq 0\\ & -1750 n_{2,2} + x_{2,2} \leq 0\\ & -4000 n_{3,2} + x_{3,2} \leq 0\\ & -2000 n_{1,3} + x_{1,3} \leq 0\\ & -1750 n_{2,3} + x_{2,3} \leq 0\\ & -4000 n_{3,3} + x_{3,3} \leq 0\\ & -2000 n_{1,4} + x_{1,4} \leq 0\\ & [[\ldots\text{60 constraints skipped}\ldots]] \\ & n_{2,4} \leq 10\\ & n_{3,4} \leq 5\\ & n_{1,5} \leq 12\\ & n_{2,5} \leq 10\\ & n_{3,5} \leq 5\\ & s_{1,1} \leq 12\\ & s_{2,1} \leq 10\\ & s_{3,1} \leq 5\\ & s_{1,2} \leq 12\\ & s_{2,2} \leq 10\\ & s_{3,2} \leq 5\\ & s_{1,3} \leq 12\\ & s_{2,3} \leq 10\\ & s_{3,3} \leq 5\\ & s_{1,4} \leq 12\\ & s_{2,4} \leq 10\\ & s_{3,4} \leq 5\\ & s_{1,5} \leq 12\\ & s_{2,5} \leq 10\\ & s_{3,5} \leq 5\\ & n_{1,1} \in \mathbb{Z}\\ & n_{2,1} \in \mathbb{Z}\\ & n_{3,1} \in \mathbb{Z}\\ & n_{1,2} \in \mathbb{Z}\\ & n_{2,2} \in \mathbb{Z}\\ & n_{3,2} \in \mathbb{Z}\\ & n_{1,3} \in \mathbb{Z}\\ & n_{2,3} \in \mathbb{Z}\\ & n_{3,3} \in \mathbb{Z}\\ & n_{1,4} \in \mathbb{Z}\\ & n_{2,4} \in \mathbb{Z}\\ & n_{3,4} \in \mathbb{Z}\\ & n_{1,5} \in \mathbb{Z}\\ & n_{2,5} \in \mathbb{Z}\\ & n_{3,5} \in \mathbb{Z}\\ & s_{1,1} \in \mathbb{Z}\\ & s_{2,1} \in \mathbb{Z}\\ & s_{3,1} \in \mathbb{Z}\\ & s_{1,2} \in \mathbb{Z}\\ & s_{2,2} \in \mathbb{Z}\\ & s_{3,2} \in \mathbb{Z}\\ & s_{1,3} \in \mathbb{Z}\\ & s_{2,3} \in \mathbb{Z}\\ & s_{3,3} \in \mathbb{Z}\\ & s_{1,4} \in \mathbb{Z}\\ & s_{2,4} \in \mathbb{Z}\\ & s_{3,4} \in \mathbb{Z}\\ & s_{1,5} \in \mathbb{Z}\\ & s_{2,5} \in \mathbb{Z}\\ & s_{3,5} \in \mathbb{Z}\\ \end{aligned} \end{split}\]

optimize!(model)

Running HiGHS 1.7.1 (git hash: 43329e528): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 4e+03]
  Cost   [4e+00, 9e+03]
  Bound  [5e+00, 1e+01]
  RHS    [2e+04, 5e+04]
Presolving model
55 rows, 45 cols, 135 nonzeros  0s
55 rows, 45 cols, 135 nonzeros  0s

Solving MIP model with:
   55 rows
   45 cols (0 binary, 30 integer, 0 implied int., 15 continuous)
   135 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   -313500         inf                  inf        0      0      0         0     0.0s
 R       0       0         0   0.00%   985514.285714   991230             0.58%        0      0      0        24     0.0s
 T       0       0         0   0.00%   986542.857143   988540             0.20%       20      4      0        30     0.0s

Solving report
  Status            Optimal
  Primal bound      988540
  Dual bound        988540
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    988540 (objective)
                    0 (bound viol.)
                    5.3290705182e-15 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     30 (total)
                    0 (strong br.)
                    6 (separation)
                    0 (heuristics)

is_solved_and_feasible(model)

true

print(objective_value(model))
print(value.(n), "\n")
print(value.(s), "\n")
print(value.(x))

988540.0

[12.0 12.0 12.0 12.0 12.0; 3.0000000000000004 8.0 8.0 8.999999999999995 8.999999999999995; -0.0 -0.0 1.2126596023639043e-15 2.0 -0.0]
[0.0 0.0 0.0 0.0 0.0; 0.0 5.0 0.0 0.9999999999999947 0.0; 0.0 0.0 1.2126596023639043e-15 1.999999999999999 0.0]
[10200.0 16000.0 10999.999999999998 21250.000000000007 11250.00000000001; 4800.0 14000.0 14000.0 15749.99999999999 15749.99999999999; 0.0 -0.0 1.8189894035458565e-12 3000.0 0.0]