Getting started with BilevelJuMP
This is a quick introduction to modeling and solving bilevel optimization with BilevelJuMP.
If you are new to Julia, start with the Getting started with Julia from the JuMP documentation.
If you are new to JuMP, start with the Getting started with JuMP from the JuMP documentation.
Installation
BilevelJuMP is a JuMP extension that be installed by using the built-in package manager.
import Pkg
Pkg.add("BilevelJuMP")
That is all you need to model a bilevel optimization problem, but we want to also solve the problems. Therefore we need a solver, one such solver is HiGHS.Optimizer
, which is provided by the HiGHS.jl package.
import Pkg
Pkg.add("HiGHS")
A first example
We will solve the following bilevel optimization problem using BilevelJuMP and HiGHS. First we take a look in the entire code then we go through it step-by-step.
Here is the example from Dempe (2002), Chapter 3.2, Page 25:
\[\begin{align*} &\min_{x, y} && 3x + y \\ &\textit{s.t.} && x \leq 5 \\ & && y \leq 8 \\ & && y \geq 0 \\ & && x(y) \in \begin{aligned}[t] &\arg\min_{x} && -x\\ &\textit{s.t.} && x + y \leq 8\\ & && 4x + y \geq 8\\ & && 2x + y \leq 13\\ & && 2x - 7y \leq 0 \end{aligned} \end{align*}\]
Here is the complete code to model, solve and query results from the example:
julia> using BilevelJuMP
julia> using HiGHS
julia> model = BilevelModel(
HiGHS.Optimizer,
mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100))
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.BigMMode{Float64}(false, 100.0, 100.0, BilevelJuMP.ComplementBoundCache(Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.ConstraintIndex, BilevelJuMP.BilevelConstraintInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}()))
Solver name: HiGHS
julia> @variable(Lower(model), x)
x
julia> @variable(Upper(model), y)
y
julia> @objective(Upper(model), Min, 3x + y)
3 x + y
julia> @constraint(Upper(model), u1, x <= 5)
u1 : x ≤ 5
julia> @constraint(Upper(model), u2, y <= 8)
u2 : y ≤ 8
julia> @constraint(Upper(model), u3, y >= 0)
u3 : y ≥ 0
julia> @objective(Lower(model), Min, -x)
-x
julia> @constraint(Lower(model), l1, x + y <= 8)
l1 : x + y ≤ 8
julia> @constraint(Lower(model), l2, 4x + y >= 8)
l2 : 4 x + y ≥ 8
julia> @constraint(Lower(model), l3, 2x + y <= 13)
l3 : 2 x + y ≤ 13
julia> @constraint(Lower(model), l4, 2x - 7y <= 0)
l4 : 2 x - 7 y ≤ 0
julia> print(model)
Upper level:
Min 3 x + y
Subject to
u3 : y ≥ 0
u1 : x ≤ 5
u2 : y ≤ 8
Lower level:
Min -x
Subject to
l2 : 4 x + y ≥ 8
l1 : x + y ≤ 8
l3 : 2 x + y ≤ 13
l4 : 2 x - 7 y ≤ 0
julia> optimize!(model)
Running HiGHS 1.5.1 [date: 1970-01-01, git hash: 93f1876e4]
Copyright (c) 2023 HiGHS under MIT licence terms
WARNING: Row 0 has infeasibility of 1 from [lower, value, upper] = [ -1; 0; -1]
WARNING: Row 2 has infeasibility of 8 from [lower, value, upper] = [ 8; 0; inf]
Solution has num max sum
Col infeasibilities 0 0 0
Integer infeasibilities 0 0 0
Row infeasibilities 4 13 30
Row residuals 0 0 0
Attempting to find feasible solution of continuous variables for user-supplied values of discrete variables
Presolving model
13 rows, 6 cols, 24 nonzeros
Problem status detected on presolve: Infeasible
Model status : Infeasible
Objective value : 0.0000000000e+00
HiGHS run time : 0.00
Presolving model
13 rows, 10 cols, 32 nonzeros
13 rows, 10 cols, 32 nonzeros
Solving MIP model with:
13 rows
10 cols (4 binary, 0 integer, 0 implied int., 6 continuous)
32 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% 0 inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 6.133333333 6.133333333 0.00% 0 0 0 5 0.0s
Solving report
Status Optimal
Primal bound 6.13333333333
Dual bound 6.13333333333
Gap 0% (tolerance: 0.01%)
Solution status feasible
6.13333333333 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 6 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1
julia> primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
julia> dual_status(Lower(model))
FEASIBLE_POINT::ResultStatusCode = 1
julia> dual_status(Upper(model))
NO_SOLUTION::ResultStatusCode = 0
julia> objective_value(model)
6.133333333333333
julia> objective_value(Lower(model))
-1.8666666666666667
julia> objective_value(Upper(model))
6.133333333333333
julia> value(x)
1.8666666666666667
julia> value(y)
0.5333333333333332
julia> dual(l1)
1-element Vector{Float64}:
0.0
julia> dual(l2)
1-element Vector{Float64}:
49.75
Step-by-step
Once installed, BilevelJuMP can be loaded into julia:
julia> using BilevelJuMP
Note that JuMP comes inside BilevelJuMP, and does not need to be installed separately. Once loaded, all JuMP functions are exported along with the BilevelJuMP additional functions.
We include a solver, in this case HiGHS:
julia> using HiGHS
Just like regular JuMP has the Model
function to initialize an optimization problem, BilevelJuMP has the BilevelModel
function that takes a solver as a first positional argument, in this case HiGHS.Optimizer
and a mode
keyword argument that selects a bilevel solution method, in this case, BilevelJuMP.FortunyAmatMcCarlMode
. Note that BilevelJuMP.FortunyAmatMcCarlMode
takes two optional keyword arguments: primal_big_M
and dual_big_M
which have to be larger than the value of all primal and dual variavle sof the lower level respectively to guarantee that the solution is no eliminated.
julia> model = BilevelModel(
HiGHS.Optimizer,
mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100))
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.BigMMode{Float64}(false, 100.0, 100.0, BilevelJuMP.ComplementBoundCache(Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.ConstraintIndex, BilevelJuMP.BilevelConstraintInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}()))
Solver name: HiGHS
For more on mode
s and solutions methods, see XXX.
We can proceed, as usual in JuMP models and incrementally build our bilevel problem. We use the same macros as JuMP.
Variables are modeled using the @variable
macro, but in bilevel problems we must define which level the variable is decided, then we use the Upper
and Lower
constructors to direct variable to the proper levels:
julia> @variable(Lower(model), x)
x
julia> @variable(Upper(model), y)
y
The same goes for objective that are modeled with the @objective
macro:
julia> @objective(Upper(model), Min, 3x + y)
3 x + y
and constraints that are modeled with the @objective
macro:
julia> @constraint(Upper(model), u1, x <= 5)
u1 : x ≤ 5
julia> @constraint(Upper(model), u2, y <= 8)
u2 : y ≤ 8
julia> @constraint(Upper(model), u3, y >= 0)
u3 : y ≥ 0
repeat for the lower level:
julia> @objective(Lower(model), Min, -x)
-x
julia> @constraint(Lower(model), l1, x + y <= 8)
l1 : x + y ≤ 8
julia> @constraint(Lower(model), l2, 4x + y >= 8)
l2 : 4 x + y ≥ 8
julia> @constraint(Lower(model), l3, 2x + y <= 13)
l3 : 2 x + y ≤ 13
julia> @constraint(Lower(model), l4, 2x - 7y <= 0)
l4 : 2 x - 7 y ≤ 0
display the model
julia> print(model)
Upper level:
Min 3 x + y
Subject to
u3 : y ≥ 0
u1 : x ≤ 5
u2 : y ≤ 8
Lower level:
Min -x
Subject to
l2 : 4 x + y ≥ 8
l1 : x + y ≤ 8
l3 : 2 x + y ≤ 13
l4 : 2 x - 7 y ≤ 0
solve the bilevel problem, which will combine a mode
(in this case FortunyAmatMcCarlMode
) and a solver (in this case HiGHS
):
julia> optimize!(model)
Running HiGHS 1.5.1 [date: 1970-01-01, git hash: 93f1876e4]
Copyright (c) 2023 HiGHS under MIT licence terms
WARNING: Row 0 has infeasibility of 1 from [lower, value, upper] = [ -1; 0; -1]
WARNING: Row 2 has infeasibility of 8 from [lower, value, upper] = [ 8; 0; inf]
Solution has num max sum
Col infeasibilities 0 0 0
Integer infeasibilities 0 0 0
Row infeasibilities 4 13 30
Row residuals 0 0 0
Attempting to find feasible solution of continuous variables for user-supplied values of discrete variables
Presolving model
13 rows, 6 cols, 24 nonzeros
Problem status detected on presolve: Infeasible
Model status : Infeasible
Objective value : 0.0000000000e+00
HiGHS run time : 0.00
Presolving model
13 rows, 10 cols, 32 nonzeros
13 rows, 10 cols, 32 nonzeros
Solving MIP model with:
13 rows
10 cols (4 binary, 0 integer, 0 implied int., 6 continuous)
32 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% 0 inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 6.133333333 6.133333333 0.00% 0 0 0 5 0.0s
Solving report
Status Optimal
Primal bound 6.13333333333
Dual bound 6.13333333333
Gap 0% (tolerance: 0.01%)
Solution status feasible
6.13333333333 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 6 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
check the termination_status
to understand why the solver stopped:
julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1
check the primal_status
to check if there is a feasible solution available:
julia> primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
check the dual_status
to check if there is a dual solution available for the lower level:
julia> dual_status(Lower(model))
FEASIBLE_POINT::ResultStatusCode = 1
do the same for the upper level:
julia> dual_status(Upper(model))
NO_SOLUTION::ResultStatusCode = 0
Most method will not support upper level duals.
JuMP's dual_status
is not available to BilevelModel
's although you can query dual_status
of each level.
Query the objecive value of the bilevel model
julia> objective_value(model)
6.133333333333333
Query the objective value of the lower level and the upper level
julia> objective_value(Lower(model))
-1.8666666666666667
julia> objective_value(Upper(model))
6.133333333333333
Obtain primal solutions:
julia> value(x)
1.8666666666666667
julia> value(y)
0.5333333333333332
and dual solutions:
julia> dual(l1)
1-element Vector{Float64}:
0.0
julia> dual(l2)
1-element Vector{Float64}:
49.75
Model basics
We created a BilevelModel passing the optimizer and mode and initialization:
julia> model = BilevelModel(
HiGHS.Optimizer,
mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100))
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.BigMMode{Float64}(false, 100.0, 100.0, BilevelJuMP.ComplementBoundCache(Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.ConstraintIndex, BilevelJuMP.BilevelConstraintInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}()))
Solver name: HiGHS
We could do piece by piece
julia> model = BilevelModel()
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.NoMode{Float64}
No solver attached
julia> set_optimizer(model, HiGHS.Optimizer)
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.NoMode{Float64}
Solver name: HiGHS
julia> BilevelJuMP.set_mode(model,
BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100))
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.BigMMode{Float64}(false, 100.0, 100.0, BilevelJuMP.ComplementBoundCache(Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}(), Dict{MathOptInterface.ConstraintIndex, BilevelJuMP.BilevelConstraintInfo}(), Dict{MathOptInterface.VariableIndex, BilevelJuMP.BilevelVariableInfo}()))
Solver name: HiGHS
Both BilevelModel
and set_optimizer
take a optimizer constructor, in this case HiGHS.Optimizer
. Note that HiGHS.Optimizer()
returns an instance of the HiGHS.Optimizer
. Hence, and alternative way to pass this solver would be: set_optimizer(model, () -> HiGHS.Optimizer())
.
() -> HiGHS.Optimizer()
is a an anonymous function that returns an instance of the HiGHS.Optimizer
.
There is no equivalent of JuMP's direct_model
in BilevelJuMP.
Solver options
it is also possible to pass optimizers with attributes:
julia> model = BilevelModel(
optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false))
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.SOS1Mode{Float64}()
Solver name: HiGHS
or set such attributes separately:
julia> model = BilevelModel(HiGHS.Optimizer)
An Abstract JuMP Model
Feasibility problem with:
Variables: 0
Upper Constraints: 0
Lower Constraints: 0
Bilevel Model
Solution method: BilevelJuMP.SOS1Mode{Float64}()
Solver name: HiGHS
julia> set_attribute(model, "output_flag", false)
julia> get_attribute(model, "output_flag")
false
This page was generated using Literate.jl.