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.03 (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 modes 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.01 (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

Info

Most method will not support upper level duals.

Info

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

Warning

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.

Info

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

Info

View this file on Github.

This page was generated using Literate.jl.