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 HiGHSBilevelJuMP.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 getDualRayMore 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 = -50then 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 variables3-element Vector{Any}:
anon
anon
anonThis page was generated using Literate.jl.