Dual variables of the lower level

This is a quick introduction to modeling and solving bilevel optimization with BilevelJuMP.

This modeling feature enables the implementation of workflows where one (or more) of the upper-level variables is the dual of a lower-level constraint. In particular, in the energy sector, it is common to model the energy prices as the dual variable associated with the energy demand equilibrium constraint. One example of an application that uses this feature is Fanzeres et al. (2019), which focuses on strategic bidding in auction-based energy markets. A small and simplified example of the modeled problem would be the model:

\[\begin{align} &\max_{\lambda, q_S} \quad \lambda \cdot g_S \\ &\textit{s.t.} \quad 0 \leq q_S \leq 100\\ &\hspace{28pt} (g_S, \lambda) \in \arg\min_{g_S, g_{1}, g_{2}, g_D} 50 g_{R1} + 100 g_{R2} + 1000 g_{D}\\ & \hspace{70pt} \textit{s.t.} \quad g_S \leq q_S \\ & \hspace{88pt} \quad 0 \leq g_S \leq 100 \\ & \hspace{88pt}\quad 0 \leq g_{1} \leq 40 \\ & \hspace{88pt}\quad 0 \leq g_{2} \leq 40 \\ & \hspace{88pt}\quad 0 \leq g_{D} \leq 100 \\ & \hspace{88pt}\quad g_S + g_{1} + g_{2} + g_D = 100 \quad : \quad \lambda \label{eq-dual-lambda} \end{align}\]

Where $\lambda$ is the dual of the load balance constraint (last constraint in the lower part), $g_S$, $g_{1}$, $g_2$ represent the generation of the strategic bidder and from two other (non-strategic) plants. $g_D$ represents the deficit in generation. Finally, $q_S$ is the quantity bid optimized by the strategic generator.

load packages

using BilevelJuMP
using Ipopt
using QuadraticToBinary
using HiGHS

BilevelJuMP.jl allows users to implement similar models using the function DualOf that binds a new variable in the upper level to an existing constraint in the lower level. The model can be written as:

model = BilevelModel()
@variable(Upper(model), 0 <= qS <= 100)
@variable(Lower(model), 0 <= gS <= 100)
@variable(Lower(model), 0 <= gR1 <= 40)
@variable(Lower(model), 0 <= gR2 <= 40)
@variable(Lower(model), 0 <= gD <= 100)
@objective(Lower(model), Min, 50gR1 + 100gR2 + 1000gD)
@constraint(Lower(model), gS <= qS)
@constraint(Lower(model), demand_equilibrium, gS + gR1 + gR2 + gD == 100)
@variable(Upper(model), lambda, DualOf(demand_equilibrium))
@objective(Upper(model), Max, lambda*gS)

\[ lambda\times gS \]

NLP solution

This model, can be solved by selecting a reformulation and a solver. Here we select Strong-Duality reformulation, the Ipopt solver and call optimizes to perform the reformulation and solve it.

BilevelJuMP.set_mode(model, BilevelJuMP.StrongDualityMode())
set_optimizer(model, Ipopt.Optimizer)
optimize!(model)
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:       17
Number of nonzeros in inequality constraint Jacobian.:       22
Number of nonzeros in Lagrangian Hessian.............:        3

Total number of variables............................:       15
                     variables with only lower bounds:        4
                variables with lower and upper bounds:        5
                     variables with only upper bounds:        5
Total number of equality constraints.................:        5
Total number of inequality constraints...............:        3
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        2

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.00e+03 1.28e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.4967374e-05 1.00e+03 1.27e+00  -1.0 3.44e+02    -  2.88e-05 2.90e-05h  1
   2  8.1838617e-04 1.00e+03 1.67e+00  -1.0 9.37e+02    -  2.13e-05 4.45e-04f  1
   3  2.6260004e-03 9.99e+02 1.64e+00  -1.0 9.16e+02    -  5.75e-04 5.48e-04h  1
   4  4.2939366e-01 9.72e+02 3.70e+00  -1.0 9.08e+02    -  7.13e-04 2.68e-02f  1
   5  8.9922740e-01 9.69e+02 3.69e+00  -1.0 7.99e+02    -  4.77e-03 3.50e-03f  1
   6  1.5988102e+00 9.68e+02 3.69e+00  -1.0 8.64e+02    -  4.71e-04 3.58e-04f  1
   7  8.1726295e+00 9.67e+02 3.65e+00  -1.0 1.79e+03    -  3.47e-04 1.12e-03f  1
   8  9.1199253e+00 9.67e+02 3.70e+00  -1.0 9.95e+02    -  8.44e-04 2.73e-04f  1
   9  1.1465721e+01 9.66e+02 3.60e+00  -1.0 2.01e+03    -  8.77e-05 6.96e-04f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.4960863e+01 9.64e+02 3.55e+00  -1.0 6.96e+02   0.0 1.22e-03 2.64e-03f  1
  11  1.6509553e+01 9.63e+02 5.48e+00  -1.0 6.42e+02  -0.5 9.07e-03 9.47e-04f  1
  12  2.0601057e+01 9.61e+02 5.71e+00  -1.0 5.08e+02  -1.0 2.37e-03 1.83e-03f  1
  13  2.8419542e+01 9.61e+02 5.09e+00  -1.0 1.48e+06    -  9.26e-07 1.35e-06f  1
  14  4.1802970e+01 9.61e+02 5.47e+00  -1.0 4.68e+06    -  8.13e-07 7.30e-07f  1
  15  7.4957309e+01 9.61e+02 4.56e+00  -1.0 2.03e+07    -  3.72e-07 4.18e-07f  1
  16  1.3881250e+02 9.61e+02 5.16e+00  -1.0 9.66e+06    -  1.75e-06 1.69e-06f  1
  17  1.5193761e+02 9.61e+02 5.21e+01  -1.0 1.00e+06    -  5.04e-05 3.34e-06f  1
  18  1.5247861e+02 9.61e+02 1.20e+02  -1.0 1.57e+03    -  4.47e-02 5.47e-05f  1
  19  8.6224578e+04 6.37e+04 7.78e+02  -1.0 9.64e+02    -  4.68e-04 9.24e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  8.6928452e+04 6.37e+04 7.83e+02  -1.0 1.39e+06    -  4.65e-05 1.36e-04f  1
  21  8.6931703e+04 6.37e+04 7.49e+02  -1.0 1.16e+04    -  7.85e-02 1.58e-04f  1
  22  8.6310221e+04 6.22e+04 7.75e+02  -1.0 4.34e+02    -  3.75e-01 2.32e-02h  1
  23  8.6303985e+04 6.22e+04 9.11e+02  -1.0 4.25e+02    -  3.24e-01 2.37e-04h  1
  24  8.6265395e+04 6.21e+04 7.98e+02  -1.0 7.31e+01    -  8.46e-06 7.03e-04h  1
  25  6.8704198e+04 4.52e+04 5.79e+02  -1.0 8.35e+01    -  6.37e-04 2.67e-01h  1
  26  5.4467029e+04 3.14e+04 4.06e+02  -1.0 5.85e+01    -  3.38e-01 2.99e-01h  1
  27  5.4323041e+04 3.13e+04 4.45e+02  -1.0 3.94e+01    -  2.91e-01 4.44e-03h  1
  28  5.4102419e+04 3.11e+04 4.01e+02  -1.0 4.12e+01    -  5.49e-05 6.40e-03h  1
  29  3.3729198e+04 1.25e+04 2.41e+02  -1.0 4.12e+01    -  2.15e-02 5.81e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  2.5938021e+04 5.43e+03 8.48e+01  -1.0 1.60e+01    -  9.96e-01 5.61e-01h  1
  31  2.0420516e+04 3.86e+02 6.64e+00  -1.0 6.82e+00    -  1.00e+00 9.22e-01h  1
  32  2.0336070e+04 3.09e+02 2.51e+02  -1.0 5.14e-01    -  1.00e+00 2.00e-01h  1
  33  2.0109529e+04 1.01e+02 1.91e+02  -1.0 4.11e-01    -  1.00e+00 6.72e-01h  1
  34  2.0046678e+04 4.36e+01 5.31e+02  -1.0 1.35e-01    -  1.00e+00 5.69e-01h  1
  35  2.0018605e+04 1.79e+01 1.20e+03  -1.0 5.80e-02    -  1.00e+00 5.90e-01h  1
  36  2.0007221e+04 7.43e+00 2.93e+03  -1.0 2.38e-02    -  1.00e+00 5.84e-01h  1
  37  2.0002490e+04 3.07e+00 7.02e+03  -1.0 9.86e-03    -  1.00e+00 5.86e-01h  1
  38  2.0000590e+04 1.27e+00 1.70e+04  -1.0 4.00e-03    -  1.00e+00 5.86e-01h  1
  39  1.9999939e+04 5.66e-01 4.39e+04  -1.0 1.37e-03    -  1.00e+00 5.56e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  1.9999916e+04 5.10e-01 2.05e+05  -1.0 2.97e-03    -  1.00e+00 9.97e-02f  3
  41  1.9999819e+04 1.55e-01 1.25e+05  -1.0 3.12e-03    -  1.00e+00 6.96e-01h  1
  42  1.9999819e+04 1.43e-01 7.31e+05  -1.0 1.22e-03    -  1.00e+00 7.51e-02f  4
  43  1.9999806e+04 2.90e-02 2.92e+05  -1.0 1.08e-03    -  1.00e+00 7.98e-01h  1
  44  1.9999805e+04 2.55e-02 3.26e+06  -1.0 1.66e-04    -  1.00e+00 1.21e-01f  3
  45  1.9999800e+04 6.64e-03 1.82e+06  -1.0 1.64e-04    -  1.00e+00 7.39e-01h  1
  46  1.9999800e+04 6.05e-03 1.31e+07  -1.0 5.25e-05    -  1.00e+00 8.04e-02f  4
  47  1.9999799e+04 1.96e-04 3.47e+06  -1.0 4.59e-05    -  1.00e+00 8.59e-01h  1
  48  1.9999799e+04 1.99e-04 2.80e+07  -1.0 9.79e-06    -  1.00e+00 3.05e-01f  2
  49  1.9999810e+04 9.97e-03 8.58e-03  -1.0 1.26e-05    -  1.00e+00 1.00e+00s 22
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  1.9999998e+04 3.36e-07 4.00e+00  -5.7 4.94e-03    -  1.00e+00 1.00e+00f  1
  51  1.9999999e+04 3.42e-11 4.43e+04  -5.7 3.90e-05    -  9.52e-01 9.98e-01f  1
  52  2.0000000e+04 5.17e-12 4.57e+02  -5.7 4.75e-06    -  1.00e+00 8.49e-01f  1
  53  2.0000000e+04 8.33e-14 1.03e-07  -5.7 5.91e-07    -  1.00e+00 1.00e+00f  1
  54  2.0000000e+04 1.14e-13 7.48e+01  -8.6 9.62e-08    -  9.59e-01 1.00e+00f  1
Cannot call restoration phase at almost feasible point,
but acceptable point from iteration 53 could be restored.

Number of Iterations....: 54

                                   (scaled)                 (unscaled)
Objective...............:  -2.0000000028637372e+04    2.0000000028637372e+04
Dual infeasibility......:   1.0265369750465189e-07    1.0265369750465189e-07
Constraint violation....:   8.3325260238143581e-14    8.3325260238143581e-14
Variable bound violation:   3.7706033850781751e-07    3.7706033850781751e-07
Complementarity.........:   2.1445435598259525e-06    2.1445435598259525e-06
Overall NLP error.......:   7.6462984588833153e-07    2.1445435598259525e-06


Number of objective function evaluations             = 91
Number of objective gradient evaluations             = 56
Number of equality constraint evaluations            = 91
Number of inequality constraint evaluations          = 91
Number of equality constraint Jacobian evaluations   = 56
Number of inequality constraint Jacobian evaluations = 56
Number of Lagrangian Hessian evaluations             = 55
Total seconds in IPOPT                               = 0.064

EXIT: Solved To Acceptable Level.

MIP solution

It is also possible to solve such problem by using a MIP formulation. The main issue is the product of variable in the upper level objective. However, this can be easily handled by using the package QuadraticToBinary.jl for automatic binary expansions. Because binary expansions require bounds on variables, we add the following lines:

set_lower_bound(lambda, 0.0)
set_upper_bound(lambda, 1000.0)

Then, as before, we set a solver (now HiGHS with the QuadraticToBinary.jl wrapper) and a solution method (now Fortuny-Amat and McCarl):

set_optimizer(model,
    ()->QuadraticToBinary.Optimizer{Float64}(HiGHS.Optimizer()))
BilevelJuMP.set_mode(model,
    BilevelJuMP.FortunyAmatMcCarlMode(dual_big_M = 100))
optimize!(model)
Running HiGHS 1.5.1 [date: 1970-01-01, git hash: 93f1876e4]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
71 rows, 69 cols, 193 nonzeros
Presolve: Infeasible

Solving report
  Status            Infeasible
  Primal bound      -inf
  Dual bound        inf
  Gap               0% (tolerance: 0.01%)
  Solution status   -
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             0
  LP iterations     0 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)
ERROR:   No invertible representation for getDualRay

More on DualOf usage

You might have a problem where you want duals of a vector of constraints like:

@constraint(Lower(model), reserves[i=1:3], (40 - gR1) + (40 - gR2) == 10 * i)
3-element Vector{ConstraintRef{BilevelModel, Int64, ScalarShape}}:
 reserves[1] : -gR1 - gR2 = -70
 reserves[2] : -gR1 - gR2 = -60
 reserves[3] : -gR1 - gR2 = -50

then you can do

@variable(Upper(model), reserve_dual[i=1:3], DualOf(reserves[i]))
3-element Vector{BilevelVariableRef}:
 reserve_dual[1]
 reserve_dual[2]
 reserve_dual[3]

or use anonymous variables

my_duals = []
for i in 1:3
    var = @variable(Upper(model), variable_type = DualOf(reserves[i]))
    push!(my_duals, var)
end
my_duals # a vector of anonimous variables
3-element Vector{Any}:
 anon
 anon
 anon

This page was generated using Literate.jl.