# Solving mixed integer programming (MIP) problems using Julia+JuMP

Shuvomoy Das Gupta

July 10, 2021

In this blog, we will discuss how to solve a mixed integer programming (MIP) problem using Julia and JuMP.

### Problem

Let us try to write the JuMP code for the following standard form optimization problem:

\begin{aligned} & \text{minimize} && c^T x + d^T y\\ & \text{subject to} && A x + B y= f \\ & && x \succeq 0, y \succeq 0 \\ & && x \in \mathbb{R}^n, y \in \mathbb{\{0,1\}}^p, \end{aligned}

where $x,y$ are the decision variables, $A \in \mathbb{R}^{m \times n}, B \in \mathbb{R}^{m \times p}, c \in \mathbb{R}^n, d \in \mathbb{R}^p, f \in \mathbb{R}^m$. The data were randomly generated. The symbol $\succeq$ ($\preceq$) stands for element-wise greater (less) than or equal to.

### Data generation

Let us generate random data.

using Random

n = 15
p = 14
m = 13

A = randn(m,n)

B = randn(m,p)

c = abs.(randn(n,1))

d = abs.(randn(p,1))

x_rand_feas = abs.(randn(n,1))

y_rand_feas = bitrand(p,1)

f = A*x_rand_feas + B*y_rand_feas # to ensure that we have a feasible solution
using JuMP, Gurobi

### Function for solving MIP

Let us write the function, which will find an optimal solution. We have used the notions of warm-start and variable hint in the function. Information about them can be found at the following links:

function MIP_solver(A, B, c, d, f;
# if warm-start value is provided
warm_start = :off,
x_start_value = zeros(n, 1),
y_start_value = zeros(p,1),
# if variable hint is provided
variable_hint = :off,
x_hint =  zeros(n, 1),
y_hint = zeros(p,1)
)

sfMipModel = direct_model(Gurobi.Optimizer())
# using direct_model has results in smaller memory allocation

set_optimizer_attribute(sfMipModel, "MIPFocus", 3)
# If you are more interested in good quality feasible solutions, you can select MIPFocus=1.
# If you believe the solver is having no trouble finding the optimal solution, and wish to focus more
# attention on proving optimality, select MIPFocus=2.
# If the best objective bound is moving very slowly (or not at all), you may want to try MIPFocus=3 to focus on the bound.

set_optimizer_attribute(sfMipModel, "NonConvex", 2)
# "NonConvex" => 2 tells Gurobi to use its nonconvex algorithm

# some other Gurobi options relevant for MIP solving (commented here)
# -------------------------------------------------------------------

# set_optimizer_attribute(BnB_PEP_model, "MIPGap", 1e-4)
## https://www.gurobi.com/documentation/9.1/refman/mipgap2.html#parameter:MIPGap

# set_optimizer_attribute(BnB_PEP_model, "FeasibilityTol", 1e-4)
## https://www.gurobi.com/documentation/9.1/refman/feasibilitytol.html

# set_optimizer_attribute(BnB_PEP_model, "OptimalityTol", 1e-4)
## https://www.gurobi.com/documentation/9.1/refman/optimalitytol.html

@variable(sfMipModel, x[1:n] >=0)
@variable(sfMipModel, y[1:p] >= 0, Bin)

if warm_start == :on

for i in 1:n
set_start_value(x[i], x_start_value[i])
# alternatively the following code will do the same thing
# MOI.set(sfMipModel, Gurobi.VariableAttribute("Start"), x[i], x_start_value[i])
end

for i in 1:p
set_start_value(y[i], y_start_value[i])
# alternatively the following code will do the same thing
# MOI.set(sfMipModel, Gurobi.VariableAttribute("Start"), y[i], y_start_value[i])
end

end

if variable_hint == :on
for i in 1:n
MOI.set(sfMipModel, Gurobi.VariableAttribute("VarHintVal"), x[i], x_start_value[i])
MOI.set(sfMipModel, Gurobi.VariableAttribute("VarHintPri"), x[i], 10) # this number represents our confidence in the variable hint
end

for i in 1:p
MOI.set(sfMipModel, Gurobi.VariableAttribute("VarHintVal"), y[i], y_start_value[i])
MOI.set(sfMipModel, Gurobi.VariableAttribute("VarHintPri"), y[i], 10) # this number represents our confidence in the variable hint
end

end

@objective(sfMipModel, Min, sum(c[i] * x[i] for i in 1:n)+sum(d[i]*y[i] for i in 1:p))

for i in 1:m
@constraint(sfMipModel, sum(A[i,j]*x[j] for j in 1:n)+ sum(B[i,j]*y[j] for j in 1:p) == f[i])
end

optimize!(sfMipModel)

# store the important values and return them
status = JuMP.termination_status(sfMipModel)
x_sol = value.(x)
y_sol = value.(y)
obj_value = objective_value(sfMipModel)

return obj_value, x_sol, y_sol

end

### Solving from scratch

obj_value_1, x_sol_1, y_sol_1 = MIP_solver(A, B, c, d, f)
Solver output
-------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 13 rows, 29 columns and 377 nonzeros
Model fingerprint: 0xf6c02cdc
Variable types: 15 continuous, 14 integer (14 binary)
Coefficient statistics:
Matrix range     [2e-03, 3e+00]
Objective range  [4e-03, 2e+00]
Bounds range     [0e+00, 0e+00]
RHS range        [6e-01, 1e+01]
Presolve time: 0.00s
Presolved: 13 rows, 29 columns, 356 nonzeros
Variable types: 15 continuous, 14 integer (14 binary)
Presolved: 13 rows, 29 columns, 356 nonzeros

Root relaxation: objective 4.410336e+00, 16 iterations, 0.00 seconds

Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

0     0    4.41034    0    5          -    4.41034      -     -    0s
H    0     0                       7.0242168    4.41034  37.2%     -    0s
H    0     0                       6.5680788    4.41034  32.9%     -    0s
0     0    6.56808    0    4    6.56808    6.56808  0.00%     -    0s

Cutting planes:
Gomory: 5
MIR: 3
Flow cover: 6

Explored 1 nodes (24 simplex iterations) in 0.06 seconds
Thread count was 8 (of 8 available processors)

Solution count 2: 6.56808 7.02422

Optimal solution found (tolerance 1.00e-04)
Best objective 6.568078761353e+00, best bound 6.568078761353e+00, gap 0.0000%

User-callback calls 74, time in user-callback 0.00 sec

### Solving using warm-start

obj_value_2, x_sol_2, y_sol_2 = MIP_solver(A, B, c, d, f;
# if warm-start value is provided
warm_start = :on,
x_start_value = x_sol_1,
y_start_value = y_sol_1 )
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 13 rows, 29 columns and 377 nonzeros
Model fingerprint: 0xbeda5cf1
Variable types: 15 continuous, 14 integer (14 binary)
Coefficient statistics:
Matrix range     [2e-03, 3e+00]
Objective range  [4e-03, 2e+00]
Bounds range     [0e+00, 0e+00]
RHS range        [6e-01, 1e+01]

Loaded user MIP start with objective 6.56808 # this tells us that our provided warm-start values has been accepted by Gurobi as a feasible solution

Presolve time: 0.00s
Presolved: 13 rows, 29 columns, 356 nonzeros
Variable types: 15 continuous, 14 integer (14 binary)
Presolved: 13 rows, 29 columns, 356 nonzeros

Root relaxation: objective 4.410336e+00, 16 iterations, 0.00 seconds

Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

0     0    4.41034    0    5    6.56808    4.41034  32.9%     -    0s

Cutting planes:
Gomory: 5
MIR: 3
Flow cover: 6

Explored 1 nodes (16 simplex iterations) in 0.08 seconds
Thread count was 8 (of 8 available processors)

Solution count 1: 6.56808

Optimal solution found (tolerance 1.00e-04)
Best objective 6.568078761353e+00, best bound 6.568078761353e+00, gap 0.0000%

User-callback calls 71, time in user-callback 0.00 sec

### Solving using variable hint

obj_value_3, x_sol_3, y_sol_3 = MIP_solver(A, B, c, d, f;
# if variable hint is provided
variable_hint = :on,
x_hint =  x_sol_1,
y_hint = y_sol_1)
Solver output
-------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 13 rows, 29 columns and 377 nonzeros
Model fingerprint: 0x9ff923d5
Variable types: 15 continuous, 14 integer (14 binary)
Coefficient statistics:
Matrix range     [2e-03, 3e+00]
Objective range  [4e-03, 2e+00]
Bounds range     [0e+00, 0e+00]
RHS range        [6e-01, 1e+01]
Presolve time: 0.00s
Presolved: 13 rows, 29 columns, 356 nonzeros
Variable types: 15 continuous, 14 integer (14 binary)
Presolved: 13 rows, 29 columns, 356 nonzeros

Root relaxation: objective 4.410336e+00, 16 iterations, 0.00 seconds

Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

0     0    4.41034    0    5          -    4.41034      -     -    0s
H    0     0                       7.0242168    4.41034  37.2%     -    0s
H    0     0                       6.5680788    4.41034  32.9%     -    0s
0     0    6.56808    0    4    6.56808    6.56808  0.00%     -    0s

Cutting planes:
Gomory: 5
MIR: 3
Flow cover: 6

Explored 1 nodes (24 simplex iterations) in 0.05 seconds
Thread count was 8 (of 8 available processors)

Solution count 2: 6.56808 7.02422

Optimal solution found (tolerance 1.00e-04)
Best objective 6.568078761353e+00, best bound 6.568078761353e+00, gap 0.0000%

User-callback calls 71, time in user-callback 0.00 sec