Dual variables of the lower level
BilevelJuMP supports use of duals of lower-problem constraints as variables in the upper problem via the DualOf() function.
DualOf() takes a named constraint from the lower problem as an argument. It is used in the expr argument of the @variable JuMP macro:
@constraint(Lower(model), some_constraint, x <= y)
@variable(Upper(model), lambda, DualOf(some_constraint))Example: strategic bidding in an energy market
In the energy sector, it is common to model market prices as the dual of the demand equilibrium constraint in the lower-level economic dispatch problem. One example of this approach is found in Fanzeres et al. (2019), which focuses on strategic energy offers in auction-based energy markets. A simplified example of the model is:
\[\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_{R1}, g_{R2}, 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_{R1} \leq 40 \\ & \hspace{88pt}\quad 0 \leq g_{R2} \leq 40 \\ & \hspace{88pt}\quad 0 \leq g_{D} \leq 100 \\ & \hspace{88pt}\quad g_S + g_{R1} + g_{R2} + g_D = 100 \quad : \quad \lambda \label{eq-dual-lambda} \end{align}\]
where:
\[S\]
is the strategically-bidding asset controlled by the upper-problem agent, where:\[q_S\]
is the quantity of generation offered into the market, to be optimized by the upper-problem agent to maximize revenue, and\[g_S\]
is the quantity of generation dispatched by the lower-problem system operator, which is no greater than $q_S$
\[g_1\]
and $g_2$ are the generation of two other non-strategic, price-taking generators;\[g_D\]
is the deficit in generation; and\[\lambda\]
is the dual of the load balance constraint
To implement this model in BilevelJuMP, first load the necessary packages:
using BilevelJuMP
using Ipopt
using QuadraticToBinary
using HiGHSInstantiate the model and fully describe the lower problem:
model = BilevelModel()
@variable(Upper(model), 0 <= qS <= 100)
@variable(Lower(model), 0 <= gS <= 100)
@variable(Lower(model), 0 <= g1 <= 40)
@variable(Lower(model), 0 <= g2 <= 40)
@variable(Lower(model), 0 <= gD <= 100)
@objective(Lower(model), Min, 50g1 + 100g2 + 1000gD)
@constraint(Lower(model), gS <= qS)
@constraint(Lower(model), demand_equilibrium, gS + g1 + g2 + gD == 100)\[ gS + g1 + g2 + gD = 100 \]
The BilevelJuMP.jl function DualOf() binds a new variable in the upper level to an existing constraint in the lower level:
@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 and the Ipopt solver, and call optimize!() to perform the reformulation and solve it.
BilevelJuMP.set_mode(model, BilevelJuMP.StrongDualityMode())
set_optimizer(model, Ipopt.Optimizer)
optimize!(model)┌ Warning: primal_var_dual_quad_slack field is deprecated, use primal_var_in_quad_obj_to_dual_slack_var instead
└ @ Dualization ~/.julia/packages/Dualization/lt5ES/src/structures.jl:251
┌ Warning: primal_parameter field is deprecated, use primal_parameter_to_dual_parameter instead
└ @ Dualization ~/.julia/packages/Dualization/lt5ES/src/structures.jl:248
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.10e-03 1.31e+07 -1.0 5.25e-05 - 1.00e+00 8.04e-02f 4
47 1.9999800e+04 8.68e-04 3.47e+06 -1.0 4.59e-05 - 1.00e+00 8.59e-01h 1
48 1.9999799e+04 1.53e-04 2.79e+07 -1.0 9.72e-06 - 1.00e+00 3.06e-01f 2
49 1.9999800e+04 2.45e-05 2.36e+07 -1.0 1.28e-05 - 1.00e+00 5.00e-01h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 1.9999800e+04 7.60e-09 3.29e-03 -1.0 9.82e-06 - 1.00e+00 1.00e+00h 1
51 1.9999998e+04 4.82e-07 5.32e+02 -8.6 4.93e-03 - 1.00e+00 1.00e+00f 1
52 1.9999999e+04 2.59e-10 9.47e+03 -8.6 4.66e-05 - 9.08e-01 9.96e-01f 1
53 2.0000000e+04 1.28e-13 4.68e+00 -8.6 4.69e-06 - 9.95e-01 1.00e+00f 1
54 2.0000000e+04 1.18e-13 9.70e-12 -8.6 2.11e-09 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 54
(scaled) (unscaled)
Objective...............: -2.0000000046432971e+04 2.0000000046432971e+04
Dual infeasibility......: 9.7010151132088956e-12 9.7010151132088956e-12
Constraint violation....: 1.1812332651272291e-13 1.1812332651272291e-13
Variable bound violation: 3.9995490652700028e-07 3.9995490652700028e-07
Complementarity.........: 4.9649033286658919e-09 4.9649033286658919e-09
Overall NLP error.......: 4.3798493233700494e-09 4.9649033286658919e-09
Number of objective function evaluations = 68
Number of objective gradient evaluations = 55
Number of equality constraint evaluations = 68
Number of inequality constraint evaluations = 68
Number of equality constraint Jacobian evaluations = 55
Number of inequality constraint Jacobian evaluations = 55
Number of Lagrangian Hessian evaluations = 54
Total seconds in IPOPT = 0.021
EXIT: Optimal Solution Found.MIP solution
It is also possible to solve such problem by using a MIP formulation. The main issue is the product of variables 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)┌ Warning: primal_var_dual_quad_slack field is deprecated, use primal_var_in_quad_obj_to_dual_slack_var instead
└ @ Dualization ~/.julia/packages/Dualization/lt5ES/src/structures.jl:251
┌ Warning: primal_parameter field is deprecated, use primal_parameter_to_dual_parameter instead
└ @ Dualization ~/.julia/packages/Dualization/lt5ES/src/structures.jl:248
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 getDualRayUsing DualOf() with vectors of constraints
DualOf() can also be used to create a vector of variables in the upper problem based on a vector of named constraints in the lower problem.
DualOf() requires a named constraint as an input, for example:
@constraint(Lower(model), reserves[i=1:3], (40 - g1) + (40 - g2) == 10 * i)3-element Vector{ConstraintRef{BilevelModel, Int64, ScalarShape}}:
reserves[1] : -g1 - g2 = -70
reserves[2] : -g1 - g2 = -60
reserves[3] : -g1 - g2 = -50You can use DualOf() with the built-in JuMP method for creating vectors of named variables:
@variable(Upper(model), reserve_dual[i=1:3], DualOf(reserves[i]))3-element Vector{BilevelVariableRef}:
reserve_dual[1]
reserve_dual[2]
reserve_dual[3]You can also use DualOf() to create a vector of 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 anonymous variables3-element Vector{Any}:
anon
anon
anonThis page was generated using Literate.jl.