# BilevelOptimization

This package is a Julia toolbox based on JuMP.jl for solving bilevel optimization problems. These are encountered in various applications, including power grids, security games, market equilibria or chemical reaction optimization.

## Generic bilevel linear problems (BLP)

The generic problem can be written:

```
min_{x} F(x,y)
such that
G * x + H * y ⩽ q
y ∈ arg min_y {d^T y + x^T F * y
such that
A * x + B * y ⩽ b
y ⩾ 0
}
x_j integer ∀ j ∈ Jx
```

`x`

represents the upper-level decision variable and `y`

the lower-level one.
`y`

is thus the solution to a parametric optimization sub-problem, depending
on the value of `x`

.
The required data describing this problem are
the feasibility domains of the upper and lower level and the coefficients
of the objective functions. All these are regrouped within the `BilevelLP`

type of this package.

The formulation is made as general as possible
for the problem to remain approachable with plain Mixed-Integer Solvers
(CBC, GLPK, SCIP, Gurobi, CPLEX). For a simple linear-linear problem,
the user can set `Jx = ∅`

and `F`

as a zero matrix of appropriate dimension.
The problem can be made as complex as wanted **at the upper level**,
as long as JuMP and the solver used support the constraints and objective.

### Usage:

The main function is `build_bilevel_lp`

, which will build the JuMP model
or modify it for the bilevel problem.

The signature:

`build_blp_model(bp::BilevelOptimization.BilevelLP, solver; comp_method)`

builds the model from scratch. It will return a `Tuple`

: `(m, x, y, λ, s)`

with:

`m`

the JuMP model`x`

the upper-level variable vector`y`

the lower-level variable vector`λ`

the dual of the lower-level problem`s`

the lower-level slack variable vector

The function can also be called with a model already built:

`build_blp_model(m::JuMP.Model, bp::BilevelOptimization.BilevelLP, x, y; comp_method)`

In which case it will add the lower-level optimality constraints, it returns the same tuple.

If the user is not willing to describe the whole problem using a `BilevelLP`

,
the following signature can be used:
`build_blp_model(m::JuMP.Model, B::M, d, s; comp_method)`

With `B`

the lower-level constraint matrix, d the lower-level objective and `s`

the lower-level
slack variable. Only the KKT conditions are added to the model in that case (not the lower-level
feasibility constraints).

## Installation

The package can be installed using Julia `Pkg`

tool:

```
julia> ]
(v1.0) pkg> add BilevelOptimization
```

You will also need an optimization solver up and running with JuMP.

## Testing

Tests can be performed using `Pkg`

:

```
julia> ]
(v1.0) pkg> test BilevelOptimization
```

## API documentation

From the Julia REPL, type `?`

to show the help prompt, then type the
identifier you want the documentation for.

```
julia> import BilevelOptimization
help?> BilevelOptimization.BilevelLP
A bilevel linear optimization problem of the form:
min cx^T * x + cy^T * y
s.t. G x + H y <= q
x_j ∈ [xl_j,xu_j]
x_j ∈ ℤ ∀ j ∈ Jx
y ∈ arg min {
d^T * y + x^T * F * y
s.t. A x + B y <= b
y_j ∈ [yl_j,yu_j]
}
Note that integer variables are allowed at the upper level.
```

## Resolution method

The "hard" part of the reduction of a bilevel problem is the set of complementarity constraints of the form:

`λ ⋅ (b - Ax - By) = 0`

These constraints cannot be handled directly, different methods have been
developed in the literature and implemented in this package.
The standard way is to give a different algorithm in `build_blp_model`

:

`build_blp_model(args..., comp_method::ComplementarityMethod = my_method)`

### Special ordered sets 1

Special-ordered Sets of type 1 or SOS1 are used by default for complementarity constraints. The option to pass is:

`build_blp_model(args..., comp_method = SOS1Complementarity())`

### Dual and primal bounds

The most common technique for these constraints is the linearization of the constraint with a formulation developed in Fortuny-Amat and McCarl, 1981, using so-called big-M constraints.

`build_blp_model(args..., comp_method = BoundComplementarity(MD, MP))`

`MD`

, `MP`

are primal and dual bounds, both can be either a scalar
for one bound per variable type or an abstract vector for one bound per
variable.

### Custom method

Users can create a custom method for the complementarity constraint,
by creating a type `T <: ComplementarityMethod`

(sub-typing is optional but
helps for clarity). They also need to implement a method:

`add_complementarity_constraint(m, cm::T, s, λ)`

Within which the complementarity constraints are added to the JuMP model `m`

.

## The toll-setting problem

As a special application of the above model, the module `BilevelFlowProblems`

offers the following problem:

- The upper-level, acting as a leader of the Stackelberg game, chooses taxes to set on some arcs of a directed graph.
- The lower-level, acting as the follower, makes a minimum-cost flow with a given minimum amount from the source to the sink.
- Each arc has an invariant base cost and a tax level decided upon by the leader.

This has been investigated in the literature as the "toll-setting problem". The required data include:

- the initial cost of each arc for all
`i,j`

- which edges can be taxed by the leader for all
`i,j`

- the tax options (at which level can each edge be taxed) for all
`i,j,k`

- flow capacities of each edge
- the minimum flow the follower has to pass from source to sink

### Example:

```
init_cost = [
0. 1. 1. 4.
0. 0. 0. 1.
0. 0. 0. 1.
0. 0. 0. 0.
]
taxable_edges = [
false true true false
false false false true
false false false true
false false false false
]
tax_options = zeros(4,4,5)
for i in 1:4, j in 1:4
if taxable_edges[i,j]
tax_options[i,j,:] .= (0.0, 0.5, 1.0, 1.5, 2.0)
end
end
capacities = [
0. 3. 2. 3.
0. 0. 0. 2.
0. 0. 0. 1.
0. 0. 0. 0.
]
minflow = 3.
BilevelFlowProblem(init_cost,taxable_edges,capacities,tax_options, minflow)
(m, r, y, f, λ) = build_blp_model(bfp, optimizer_with_attributes(Cbc.Optimizer, "LogLevel" => 0))
st = JuMP.optimize!(m)
# objective_value(m) ≈ 6.
# for j in 1:size(r)[2]
# for i in 1:size(r)[1]
# @test JuMP.value(r[i,j]) ≈ sum(JuMP.value.(y[i,j,:]) .* bfp.tax_options[i,j,:]) * JuMP.value(f[i,j])
# end
# end
```

## Questions, issues, contributions

Problems with the package and its usage can be explained through Github issues, ideally with a minimal working example showing the problem. Pull requests (PR) are welcome.

Please read detailed information in **CONTRIBUTING.md**.

## Related packages

- BilevelJuMP.jl is more flexible on the form of the problem input, and has been designed for MathOptInterface from the beginning, the two packages may be merged at some point
- Complementarity.jl solving a generic class including bilevel problems using non-linear techniques
- MibS for problems where the lower-level also includes integer variables. KKT conditions can therefore not be used and other branching and cutting plane techniques are leveraged.
- YALMIP includes a bilevel solver and offers roughly the same features (and a bit more) as BilevelOptimization.jl

## Citing

See *CITATION.bib*, prefer citing the paper published in the Journal of Open-Source Software.