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.6066070303533018e-10    2.6066070303533018e-10
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5705110524849051e-09    2.5705110524849051e-09
Overall NLP error.......:   2.5705110524849051e-09    2.5705110524849051e-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.044

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)
Info

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.