Merge pull request #196 from joaquimg/jg/docs

add tutorials and reorganize package

.gitignore

@@ -5,6 +5,7 @@ deps/deps.jl
 Manifest.toml
 docs/build
 docs/src/examples/*.md
+docs/src/tutorials/*.md
 docs/src/.DS_Store
 .DS_Store
 gurobi.log
@@ -13,3 +14,5 @@ test.jl
 
 *.mps
 *.aux
+/docs/src/tutorials/conic_lower.md
+/docs/src/tutorials/lower_duals.md

Project.toml

@@ -8,9 +8,11 @@ Dualization = "191a621a-6537-11e9-281d-650236a99e60"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
 Dualization = "0.5.6"
 JuMP = "1.7.0"
 MathOptInterface = "1.2"
+Reexport = "1"
 julia = "1.6"

README.md

@@ -23,7 +23,7 @@ The currently available methods are based on re-writing the problem using the KK
 ## Example
 
 ```julia
-using JuMP, BilevelJuMP, SCIP
+using BilevelJuMP, SCIP
 
 model = BilevelModel(SCIP.Optimizer, mode = BilevelJuMP.SOS1Mode())
 
@@ -91,81 +91,3 @@ reformulation which add the constraint enforcing primal dual equality. The optio
 is `BilevelJuMP.StrongDualityMode(eps)` where `eps` is the tolerance on the enforcing
 constraint.
 
-### Note on QuadraticToBinary
-
-[QuadraticToBinary.jl](https://github.com/joaquimg/QuadraticToBinary.jl) is a
-package that converts quadratic terms in constraints and objective. To do so
-the pack acts like a solver on top of the real solver and most data is forwarded
-directly to the solver itself. For many solvers it is enough to use:
-
-```julia
-SOLVER = Xpress.Optimizer()
-Q_SOLVER = QuadraticToBinary.Optimizer{Float64}(SOLVER)
-BilevelModel(()->Q_SOLVER, mode = BilevelJuMP.ProductMode(1e-5))
-```
-
-However, this might lead to some solver not supporting certain functionality like SCIP.
-In this case we need to:
-
-```julia
-SOLVER = SCIP.Optimizer()
-CACHED_SOLVER = MOI.Utilities.CachingOptimizer(
-    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), SOLVER)
-Q_SOLVER = QuadraticToBinary.Optimizer{Float64}(CACHED_SOLVER)
-BilevelModel(()->Q_SOLVER, mode = BilevelJuMP.ProductMode(1e-5))
-```
-Note that we used `()->Q_SOLVER` instead of just `Q_SOLVER` because `BilevelModel`
-requires as constructor and not an instance of an object.
-
-## Advanced Features
-
-### Lower level dual variables
-
-Suppose you have a constraint `b` in the lower level:
-
-```julia
-@constraint(Lower(model), b, ...)
-```
-
-It is possible to access the dual variable of `b` to use it in the upper level:
-
-```julia
-@variable(Upper(model), lambda, DualOf(b))
-```
-
-### Conic lower level
-
-BilevelJuMP allows users to write conic models in the lower level. However,
-solving this kind of problems is much harder and requires complex solution
-methods. Mosek's Conic MIP can be used with the aid of
-[QuadraticToBinary.jl](https://github.com/joaquimg/QuadraticToBinary.jl).
-
-It is also possible to solve Second Order Cone constrained models with Ipopt.
-In this case we need to add a special, non-standard bridge, to Ipopt as follows:
-
-```julia
-IPO_OPT = Ipopt.Optimizer(print_level=0)
-IPO = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(IPO_OPT)
-BilevelModel(()->IPO, mode = BilevelJuMP.ProductMode(1e-5))
-```
-
-## Troubleshooting
-
-* Cbc has known bugs in its SOS1 constraints, so `BilevelJuMP.SOS1Mode` might
-not work properly with Cbc.
-
-* For anonymous variables with `DualOf` use:
-```julia
-@variable(Upper(model, variable_type = DualOf(my_lower_constraint))
-```
-
-* Nonconvex/nonconcave/nonpsd objective/constraint error in a MIP solver.
-See section on
-[`QuadraticToBinary.jl`](#note-on-quadratictobinary)
-or if you are using
-[`Gurobi`](https://github.com/jump-dev/Gurobi.jl#using-gurobi-v90-and-you-got-an-error-like-q-not-psd)
-use:
-```julia
-model = BilevelModel(Gurobi.Optimizer, mode = BilevelJuMP.SOS1Mode()) #or other mode
-set_optimizer_attribute(model, "NonConvex", 2)
-```

benchmarks/forecast.jl

@@ -1,5 +1,5 @@
 using Random
-using JuMP, BilevelJuMP
+using BilevelJuMP
 
 function bench_forecast(prods, samples, optimizer, mode, seed = 1234)
     rng = Random.MersenneTwister(seed)

benchmarks/rand.jl

@@ -1,5 +1,5 @@
 using SparseArrays, LinearAlgebra
-using JuMP, BilevelJuMP
+using BilevelJuMP
 using Random
 
 function bench_rand(rows, cols, density, optimizer, mode, seed = 1234)

benchmarks/runbench.jl

@@ -1,6 +1,6 @@
 using Random
 using Dates
-using JuMP, BilevelJuMP
+using BilevelJuMP
 
 using MathOptInterface
 MOI = MathOptInterface
@@ -149,28 +149,28 @@ SOLVERS = [
     #=
         PrimalDual BIN 10
     =#
-    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "glpk_sd10"),
-    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "mosek_sd10"),
-    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "gurobi_sd10"),
-    # (()->QB(cpx(),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "cplex_sd10"),
-    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "xpress_sd10"),
-    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "cbc_sd10"), #TODO
-    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityEqualityMode(), "scip_sd10"),
+    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "glpk_sd10"),
+    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "mosek_sd10"),
+    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "gurobi_sd10"),
+    # (()->QB(cpx(),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "cplex_sd10"),
+    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "xpress_sd10"),
+    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "cbc_sd10"), #TODO
+    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-10,ub=10), BilevelJuMP.StrongDualityMode(), "scip_sd10"),
     #=
         PrimalDual BIN 100
     =#
-    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "glpk_sd100"),
-    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "mosek_sd100"),
-    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "gurobi_sd100"),
-    # (()->QB(cpx(),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "cplex_sd100"),
-    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "xpress_sd100"),
-    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "cbc_sd100"),
-    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityEqualityMode(), "scip_sd100"),
+    # (()->QB(GLPK.Optimizer(tm_lim=MAX_TIME*1_000),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "glpk_sd100"),
+    # (()->QB(Mosek.Optimizer(MIO_MAX_TIME=MAX_TIME*1.0,OPTIMIZER_MAX_TIME=MAX_TIME*1.0),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "mosek_sd100"),
+    # (()->QB(Gurobi.Optimizer(TimeLimit=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "gurobi_sd100"),
+    # (()->QB(cpx(),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "cplex_sd100"),
+    # (()->QB(Xpress.Optimizer(MAXTIME=-MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "xpress_sd100"),
+    # (()->QB(cache(Cbc.Optimizer(seconds=MAX_TIME*1.0)),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "cbc_sd100"),
+    # (()->QB(SCIP.Optimizer(limits_time=MAX_TIME*1),lb=-100,ub=100), BilevelJuMP.StrongDualityMode(), "scip_sd100"),
     #=
         PrimalDual NLP (DONE)
     =#
-    # (with_att(Ipopt.Optimizer, "max_cpu_time" => MAX_TIME*1.0), BilevelJuMP.StrongDualityEqualityMode(), "ipopt_sd"),
-    # (new_knitro, BilevelJuMP.StrongDualityEqualityMode(), "knitro_sd"),
+    # (with_att(Ipopt.Optimizer, "max_cpu_time" => MAX_TIME*1.0), BilevelJuMP.StrongDualityMode(), "ipopt_sd"),
+    # (new_knitro, BilevelJuMP.StrongDualityMode(), "knitro_sd"),
     #=
         Product NLP (DONE)
     =#
@@ -185,8 +185,8 @@ SOLVERS = [
     =#
     # (() -> AmplNLWriter.Optimizer("bonmin"), BilevelJuMP.ProductMode(1e-7), "bonmin_prod"),
     # (() -> AmplNLWriter.Optimizer("couenne"), BilevelJuMP.ProductMode(1e-7), "couenne_prod"),
-    # (() -> AmplNLWriter.Optimizer("bonmin"), BilevelJuMP.StrongDualityEqualityMode(), "bonmin_sd"),
-    # (() -> AmplNLWriter.Optimizer("couenne"), BilevelJuMP.StrongDualityEqualityMode(), "couenne_sd"),
+    # (() -> AmplNLWriter.Optimizer("bonmin"), BilevelJuMP.StrongDualityMode(), "bonmin_sd"),
+    # (() -> AmplNLWriter.Optimizer("couenne"), BilevelJuMP.StrongDualityMode(), "couenne_sd"),
 ]
 
 PROBLEMS = [:SVR, :TOLL, :FORECAST, :RAND]

benchmarks/svr.jl

@@ -1,4 +1,4 @@
-using JuMP, BilevelJuMP
+using BilevelJuMP
 using Random
 
 function bench_svr(dim, sample, optimizer, mode, seed = 1234)

benchmarks/toll.jl

@@ -1,4 +1,4 @@
-using JuMP, BilevelJuMP
+using BilevelJuMP
 using Random
 
 function bench_toll(nodes, optimizer, mode, seed = 1234)

docs/Project.toml

@@ -2,13 +2,20 @@
 BilevelJuMP = "485130c0-026e-11ea-0f1a-6992cd14145c"
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MibS_jll = "118347d2-e127-56b9-9933-6abf9cc1adc1"
+QuadraticToBinary = "014a38d5-7acb-4e20-b6c0-4fe5c2344fd1"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Documenter = "~0.26"
+HiGHS = "=1.5.1"
 Ipopt = "=1.0.2"
 MibS_jll = "=1.1.3"
+QuadraticToBinary = "=0.4.0"
+SCIP = "=0.11.12"

docs/make.jl

@@ -4,6 +4,7 @@ using Literate
 using Test
 
 const _EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")
+const _TUTORIAL_DIR = joinpath(@__DIR__, "src", "tutorials")
 
 """
     _include_sandbox(filename)
@@ -53,6 +54,15 @@ function literate_directory(dir)
 end
 
 literate_directory(_EXAMPLE_DIR)
+literate_directory(_TUTORIAL_DIR)
+
+_REPL_FILES = ["getting_started.md",]
+for file in _REPL_FILES
+    filename = joinpath(@__DIR__, "src", "tutorials", file)
+    content = read(filename, String)
+    content = replace(content, "@example" => "@repl")
+    write(filename, content)
+end
 
 makedocs(;
     modules = [BilevelJuMP],
@@ -61,16 +71,35 @@ makedocs(;
     format = Documenter.HTML(;
         mathengine = Documenter.MathJax2(),
         prettyurls = get(ENV, "CI", nothing) == "true",
+        collapselevel = 1,
     ),
     sitename = "BilevelJuMP.jl",
     authors = "Joaquim Garcia",
     pages = [
         "Home" => "index.md",
-        "Manual" => "manual.md",
-        "Examples" => [
-            joinpath("examples", f) for
-            f in readdir(_EXAMPLE_DIR) if endswith(f, ".md")
-        ],
+        # "Manual" => "manual.md",
+        "Tutorials" => joinpath.("tutorials", [
+            "getting_started.md",
+            "modes.md",
+            "lower_duals.md",
+            "conic_lower.md",
+            "quad_to_bin.md",
+        ]),
+        "Tutorials" => joinpath.("examples", [
+            "FOBP_example2.md",
+            "FOBP_example3.md",
+            "FOBP_example4.md",
+            "FOBP_example5.md",
+            "DTMP_example1.md",
+            "PHTP_example1.md",
+            "PHTP_example2.md",
+            "SOCBLP_example1.md",
+            "MibS_example1.md",
+            "MibS_example2.md",
+        ]),
+        "Background Information" => "background.md",
+        "API Reference" => "reference.md",
+        "Troubleshooting" => "troubleshooting.md",
     ],
 )
 

docs/src/background.md

@@ -0,0 +1,58 @@
+# Bilevel Optimization
+
+Bilevel optimization is a vast discipline with a long (50+ years) history.
+We will not attempt to present a bilevel optimization introduction here.
+Instead, we point the reader to the excelent text:
+
+[A Gentle and Incomplete Introduction to Bilevel Optimization](https://optimization-online.org/wp-content/uploads/2021/06/8450-1.pdf)
+
+by Yasmine Beck and Martin Schmidt.
+
+The interested reader can find more information in books:
+
+* [Dempe 2002](https://doi.org/10.1007/b101970)
+* [Bard 2013](https://doi.org/10.1007/978-1-4757-2836-1)
+* [Dempe et al. 2015](https://doi.org/10.1007/978-3-662-45827-3)
+
+and in other reviews:
+
+* [Vicente and Calamai 1994](https://doi.org/10.1007/BF01096458)
+* [Colson et al. 2007](https://doi.org/10.1007/s10479-007-0176-2)
+* [Kalashnikov et al. 2015](https://doi.org/10.1155/2015/310301)
+* [Dempe (2018)](https://optimization-online.org/wp-content/uploads/2018/08/6773.pdf)
+
+## Bilevel Optimization in BilevelJuMP
+
+In BileveJuMP focus on the following bilevel problem form:
+
+```math
+\begin{aligned}
+    &\min_{x, \textbf{y}, z} && f_0(x, \textbf{y}, z) \\
+    &\textit{s.t.} && f_i(x, \textbf{y}, z) \in \mathcal{S}_i, \quad i = 1 \ldots k, \\
+        &&& x(z), \textbf{y}(z) \in
+     \begin{aligned}[t]
+        &\arg\min_{x, \textbf{y}} && \frac{1}{2}{[x, z]}^{\top} Q {[x, z]} + {a_0}^{\top} x+ {d_0}^{\top} z + b_0\\
+            &\textit{s.t.} && A_i x + D_i z + b_i \in \mathcal{C}_i \quad : \ y_i \quad, \quad i = 1 \ldots m,
+     \end{aligned}
+\end{aligned}
+```
+
+where ``z \in {\mathbb{R}^l}`` and ``x \in {\mathbb{R}^n}`` are, respectively, from upper and lower-level primal decision variables. Vectors of lower level dual decision variables are represented individually, by ``y_1 \in \mathbb{R}^{p_1}, \ldots, y_m \in \mathbb{R}^{p_m}``, or jointly, by the ``m``-tuple ``\textbf{y} = (y_1, \ldots, y_m)``.
+The square brackets ``[\cdot, \cdot]`` represent the stacking of two vectors or scalars. Thus,
+``[x,z]`` is a ``(n+l)``--vector with the elements of ``x`` and ``z`` stacked. The numbers of constraints in the upper and lower problems are given by ``k`` and ``m``, respectively. ``Q``, ``a_i``, ``d_i``, ``b_i``, ``A_i``, ``D_i`` are matrices (upper case) and vectors (lower case) of constants of the lower level problem and ``\mathcal{C}_i \subset \mathbb{R}^{p_i}`` are convex conic sets. The functions ``f_i`` can be linear, quadratic, or non-linear. The sets ``\mathcal{S}_i\in \mathbb{R}^{q_i}`` can be convex cones, as in the lower level, but can also represent other sets, such as the sets of integers or binary variables.
+We use the "function-in-set" notation following the MOI definition of mathematical optimization problems.
+As in traditional bilevel programming, ``z`` is decided in the upper level and passed to the lower level as a parameter and ``x`` might be seen as an upper-level variable constrained to be an optimal solution of the lower level.
+Also, we only consider optimistic bilevel problems. In short, the solution of the lower level will be the one that optimizes the upper level in case of degeneracy.
+
+For more information see our paper:
+
+```
+@article{diasgarcia2022bileveljump,
+    title={{BilevelJuMP. jl}: {M}odeling and solving bilevel optimization in {J}ulia},
+    author={{Dias Garcia}, Joaquim and Bodin, Guilherme and Street, Alexandre},
+    journal={arXiv preprint arXiv:2205.02307},
+    year={2022}
+}
+```
+
+Here is the [pdf](https://arxiv.org/pdf/2205.02307.pdf).