Conic Bilevel and Mixed Mode
Here we present a simple bilevel program with a conic lower level model described in example 3.3 from Chi, et al. 2014.
\[\begin{align} &\max_{x \in \mathbb{R}} \quad x + 3y_1 \\ &\textit{s.t.} \quad 2 \leq x \leq 6\\ & \hspace{28pt} y(x) \in \arg\min_{y\in {\mathbb{R}^3}} -y_1\\ & \hspace{58pt} \textit{s.t.} \quad x + y_1 \leq 8 \\ & \hspace{76pt} \quad x + 4y_1 \geq 8 \\ & \hspace{76pt} \quad x + 2y_1 \leq 12 \\ & \hspace{76pt} \quad y \in {SOC}_3 \label{eq-soc} \end{align}\]
In this problem most of the constraints are regular linear constraints while the last one is a second order cone constraint. Such constraint ensures that the vector $y$ belongs to a second order cone of dimension $3$, that is: $y_1 \geq \sqrt{y_2^2 + y_3^2}$.
This problem can be encoded using regular JuMP syntax for conic programs:
using BilevelJuMP
model = BilevelModel()
@variable(Upper(model), x)
@variable(Lower(model), y[i=1:3])
@objective(Upper(model), Min, x + 3y[1])
@constraint(Upper(model), x >= 2)
@constraint(Upper(model), x <= 6)
@objective(Lower(model), Min, - y[1])
@constraint(Lower(model), con1, x + y[1] <= 8)
@constraint(Lower(model), con2, x + 4y[1] >= 8)
@constraint(Lower(model), con3, x + 2y[1] <= 12)
@constraint(Lower(model), con4, y in SecondOrderCone())\[ [y_{1}, y_{2}, y_{3}] \in \text{MathOptInterface.SecondOrderCone(3)} \]
NLP solution and start values
We can set, for instance, the product reformulation and selected Ipopt as a solver. As Ipopt does not have native support for second order cones, we use the non-default MOI bridge SOCtoNonConvexQuad to convert second order cones into quadratic constraints.
using Ipopt
BilevelJuMP.set_mode(model, BilevelJuMP.ProductMode(1e-5))
set_optimizer(model,
()->MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(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...: 6
Number of nonzeros in inequality constraint Jacobian.: 43
Number of nonzeros in Lagrangian Hessian.............: 18
Total number of variables............................: 10
variables with only lower bounds: 1
variables with lower and upper bounds: 0
variables with only upper bounds: 2
Total number of equality constraints.................: 3
Total number of inequality constraints...............: 14
inequality constraints with only lower bounds: 5
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 9
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 8.00e+00 6.57e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.5275382e-01 7.82e+00 1.55e+00 -1.0 4.46e+00 - 1.13e-02 2.81e-02h 1
2 2.7495070e-01 7.67e+00 9.17e+01 -1.0 4.97e+00 - 6.55e-01 1.89e-02h 1
3 3.2925519e-01 7.60e+00 1.85e+02 -1.0 1.30e+01 - 6.90e-03 3.32e-03h 1
4 7.1724628e-01 7.09e+00 8.94e+01 -1.0 4.41e+02 - 7.39e-04 1.13e-03h 1
5 7.1999985e-01 7.09e+00 1.39e+04 -1.0 1.81e+02 2.0 5.57e-04 1.62e-05h 1
6 1.0099676e+00 6.77e+00 4.14e+03 -1.0 3.99e+02 1.5 2.63e-04 7.86e-04h 1
7 4.8649217e+00 2.56e+00 1.57e+05 -1.0 3.70e+02 1.0 3.07e-03 1.12e-02h 2
8 5.7218502e+00 1.70e+00 1.80e+05 -1.0 1.45e+02 1.5 1.88e-02 5.89e-03h 1
9 6.4378977e+00 1.16e+00 1.70e+05 -1.0 1.31e+02 1.0 2.54e-02 9.46e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 6.4476810e+00 1.15e+00 1.70e+05 -1.0 8.09e+01 0.5 1.13e-02 2.14e-04h 1
11 6.5116212e+00 1.09e+00 1.70e+05 -1.0 6.95e+01 0.9 1.42e-04 9.29e-04h 1
12 7.1915775e+00 9.31e-01 1.63e+05 -1.0 1.70e+03 0.5 4.97e-08 4.01e-04f 1
13 7.1805348e+00 9.23e-01 1.61e+05 -1.0 8.79e+00 - 1.05e-02 9.45e-03f 1
14 7.1809370e+00 9.22e-01 1.61e+05 -1.0 4.67e+00 1.8 3.66e-02 1.05e-04h 1
15 7.1931710e+00 9.22e-01 1.61e+05 -1.0 5.62e+01 - 1.32e-04 2.21e-04h 1
16r 7.1931710e+00 9.22e-01 1.00e+03 0.9 0.00e+00 - 0.00e+00 3.54e-07R 4
17r 7.1803638e+00 9.31e-01 9.96e+02 0.9 7.38e+03 - 1.72e-02 9.36e-04f 1
18 7.1802698e+00 9.27e-01 2.42e+01 -1.0 1.80e+00 1.3 1.05e-02 4.39e-03f 1
19 7.1824437e+00 9.24e-01 2.18e+03 -1.0 1.44e+00 0.8 1.74e-01 3.27e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 7.1841671e+00 9.24e-01 5.11e+05 -1.0 8.68e+00 0.4 8.64e-02 3.66e-04h 1
21 7.1782285e+00 9.24e-01 1.72e+08 -1.0 2.05e+02 - 1.33e-02 3.87e-05h 1
22r 7.1782285e+00 9.24e-01 1.00e+03 0.1 0.00e+00 - 0.00e+00 4.09e-07R 3
23r 7.2536425e+00 9.25e-01 1.21e+03 0.1 1.19e+03 - 9.06e-03 1.17e-03f 1
24 7.2536348e+00 9.25e-01 1.10e+05 -1.0 9.53e-01 -0.1 9.22e-01 1.70e-05h 1
25r 7.2536348e+00 9.25e-01 1.00e+03 -0.0 0.00e+00 - 0.00e+00 2.59e-07R 6
26r 7.5343972e+00 9.34e-01 1.22e+03 -0.0 9.54e+01 - 8.79e-03 3.97e-03f 1
27r 7.5327092e+00 9.46e-01 9.90e+02 -0.0 1.70e+00 2.0 1.33e-01 1.52e-02f 1
28r 7.5181606e+00 9.80e-01 9.71e+02 -0.0 7.28e-01 2.4 2.49e-03 1.46e-01f 1
29r 7.5135600e+00 9.88e-01 9.55e+02 -0.0 7.04e-01 2.9 6.85e-02 2.06e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 7.5135662e+00 9.88e-01 2.66e+03 -0.0 8.57e-02 3.3 1.00e+00 1.86e-02h 1
31r 7.5227783e+00 9.94e-01 9.97e+00 -0.0 1.58e-02 2.8 1.00e+00 1.00e+00f 1
32r 7.5237814e+00 1.00e+00 1.29e+00 -1.4 6.12e-03 2.3 1.00e+00 1.00e+00f 1
33r 7.5248565e+00 1.00e+00 1.37e-01 -2.1 1.95e-03 1.8 1.00e+00 1.00e+00f 1
34r 7.5265535e+00 1.00e+00 7.02e+01 -3.2 3.47e-03 - 1.00e+00 7.46e-01f 1
35r 9.2295707e+00 1.00e+00 6.44e+02 -3.2 2.52e+00 - 4.06e-01 9.10e-01f 1
36r 9.1585683e+00 1.00e+00 2.89e+02 -3.2 2.61e-01 1.4 6.81e-01 1.00e+00h 1
37r 9.1584002e+00 1.00e+00 8.13e+01 -3.2 1.15e-03 1.8 5.41e-01 8.53e-01h 1
38r 9.1611165e+00 1.00e+00 4.80e+01 -3.2 3.20e-03 1.3 1.57e-01 1.00e+00f 1
39r 9.1870694e+00 9.90e-01 5.54e+01 -3.2 3.68e-02 1.7 3.29e-02 8.23e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 9.1893066e+00 9.72e-01 4.56e+01 -3.2 1.82e-02 3.1 6.00e-01 1.00e+00f 1
41r 9.1979046e+00 9.52e-01 2.81e+01 -3.2 6.66e-02 2.6 3.62e-01 3.05e-01f 1
42r 9.3297133e+00 8.04e-01 7.97e+01 -3.2 1.08e+00 2.1 5.48e-02 1.62e-01f 1
43 9.3307037e+00 8.04e-01 6.64e+03 -1.0 2.93e+01 -0.6 6.97e-02 4.51e-05f 1
44 9.9071010e+00 7.70e-01 1.39e+04 -1.0 1.82e+01 -1.1 1.23e-01 4.22e-02f 2
45 1.1984910e+01 5.61e-01 2.05e+04 -1.0 1.22e+01 -0.6 3.78e-01 2.72e-01H 1
46 1.2009065e+01 3.55e-01 9.73e+04 -1.0 5.61e-01 0.7 1.00e+00 3.66e-01h 1
47 1.2007891e+01 2.63e-01 2.24e+05 -1.0 3.55e-01 1.1 1.00e+00 2.59e-01h 1
48 1.2015766e+01 1.20e-01 2.26e+05 -1.0 2.63e-01 0.6 1.00e+00 5.43e-01h 1
49 1.2021975e+01 5.78e-02 3.75e+05 -1.0 1.20e-01 1.1 1.00e+00 5.20e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 1.2032201e+01 2.04e-02 4.94e+05 -1.0 7.12e-02 0.6 1.00e+00 6.46e-01h 1
51 1.2053842e+01 6.82e-03 8.67e+05 -1.0 7.96e-02 0.1 1.00e+00 6.66e-01h 1
52 1.2068410e+01 8.41e-04 4.74e+05 -1.0 3.85e-02 0.5 1.00e+00 8.77e-01h 1
53 1.2084847e+01 0.00e+00 5.33e+00 -1.0 3.41e-02 0.1 1.00e+00 1.00e+00f 1
54 1.2023010e+01 2.22e-16 1.51e+05 -5.7 1.26e-01 -0.4 9.84e-01 1.00e+00f 1
55 1.1999979e+01 0.00e+00 2.91e+02 -5.7 5.58e-02 - 9.98e-01 8.16e-01h 1
56 1.1999973e+01 0.00e+00 4.25e-06 -5.7 3.02e-04 - 1.00e+00 1.00e+00h 1
57 1.1999970e+01 0.00e+00 5.91e+00 -8.6 5.01e-06 - 9.68e-01 1.00e+00h 1
58 1.1999970e+01 0.00e+00 2.61e-10 -8.6 1.18e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 58
(scaled) (unscaled)
Objective...............: 1.1999969855039765e+01 1.1999969855039765e+01
Dual infeasibility......: 2.6066157088495915e-10 2.6066157088495915e-10
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5705110524827392e-09 2.5705110524827392e-09
Overall NLP error.......: 2.5705110524827392e-09 2.5705110524827392e-09
Number of objective function evaluations = 77
Number of objective gradient evaluations = 43
Number of equality constraint evaluations = 77
Number of inequality constraint evaluations = 77
Number of equality constraint Jacobian evaluations = 62
Number of inequality constraint Jacobian evaluations = 62
Number of Lagrangian Hessian evaluations = 58
Total seconds in IPOPT = 0.031
EXIT: Optimal Solution Found.The user could also use the alternative JuMP syntax:
set_start_value(x, 6)
set_dual_start_value(con2, 0)0
MIP solution and mixed mode
Alternatively, we could have used a Mixed Integer Second Order Cone Program (MISOCP) solver together with binary expansions. Complementarity of conic constraints is more difficult to handle because they require a sum of products that cannot be reformulated with other methods. Therefore, we rely on product reformulation for conic constraints. However, we can use other reformulations like indicator constraints for the non-conic constraints. Mixing the two of them can be done with Mixed Mode.
The following code describes how to solve the problem with a MISOCP based solver.
using Xpress
using QuadraticToBinary
set_optimizer(model,
()->QuadraticToBinary.Optimizer{Float64}(Xpress.Optimizer(),lb=-10,ub=10))
BilevelJuMP.set_mode(model,
BilevelJuMP.MixedMode(default = BilevelJuMP.IndicatorMode()))
BilevelJuMP.set_mode(con4, BilevelJuMP.ProductMode(1e-5))
optimize!(model)This code was not executed because Xpress requires a commercial license. Other solvers supporting MISOCP could be used such as Gurobi and CPLEX.
We set the reformulation method as Mixed Mode and selected Indicator constraints to be the default for the case in which we do not explicitly specify the reformulation. Then we set product mode for the second order cone reformulation.
Binary expansions require bounded variables, hence the QuadraticToBinary meta-solver accepts fallback to upper and lower bounds (\texttt{ub} and \texttt{lb}), used for variables with no explicit bounds.
This page was generated using Literate.jl.