Add documentation

docs/Project.toml

@@ -0,0 +1,7 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+MathOptComplements = "3144fb79-1414-4d28-b1d1-63dadd798e24"
+PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"

docs/make.jl

@@ -0,0 +1,33 @@
+using Documenter
+using DocumenterInterLinks
+using MathOptComplements
+
+links = InterLinks(
+    "MathOptInterface" => "https://jump.dev/MathOptInterface.jl/stable/",
+)
+
+makedocs(
+    sitename = "MathOptComplements.jl",
+    format = Documenter.HTML(
+        assets = ["assets/favicon.ico"],
+        prettyurls = Base.get(ENV, "CI", nothing) == "true",
+        mathengine = Documenter.KaTeX()
+    ),
+    modules = [MathOptComplements],
+    repo = "https://github.com/blegat/MathOptComplements.jl/blob/{commit}{path}#{line}",
+    checkdocs = :none,
+    clean=true,
+    pages = [
+        "Home" => "index.md",
+        "Quickstart" => "quickstart.md",
+        "Tutorials" => [
+            "Equilibrium problem" => "equilibrium.md",
+        ],
+        "API reference" => [
+            "Reformulations" => "api/reformulation.md",
+            "Relaxations" => "api/relaxation.md",
+            "Bridges" => "api/bridges.md",
+        ]
+    ],
+    plugins = [links],
+)

docs/src/api/bridges.md

@@ -0,0 +1,33 @@
+```@meta
+CurrentModule = MathOptComplements
+```
+
+# Complementarity bridges
+
+```@docs
+
+Bridges.SpecifySetTypeBridge
+Bridges.VerticalBridge
+Bridges.ComplementsVectorizeBridge
+Bridges.SplitIntervalBridge
+Bridges.FlipSignBridge
+```
+
+
+## SOS1 reformulation
+
+```@docs
+Bridges.ToSOS1Bridge
+```
+
+## Nonlinear reformulation
+
+```@docs
+Bridges.NonlinearBridge
+```
+
+## Utilities
+
+```@docs
+Bridges.add_all_bridges
+```

docs/src/api/reformulation.md

@@ -0,0 +1,13 @@
+```@meta
+CurrentModule = MathOptComplements
+```
+
+# Reformulations
+
+```@docs
+
+ComplementsWithSetType
+DefaultComplementarityReformulation
+ComplementarityReformulation
+
+```

docs/src/api/relaxation.md

@@ -0,0 +1,14 @@
+```@meta
+CurrentModule = MathOptComplements
+```
+
+# Nonlinear relaxations
+
+```@docs
+AbstractComplementarityRelaxation
+Bridges.ScholtesRelaxation
+Bridges.FischerBurmeisterRelaxation
+Bridges.LiuFukushimaRelaxation
+Bridges.KanzowSchwarzRelaxation
+
+```

docs/src/equilibrium.md

@@ -0,0 +1,78 @@
+```@meta
+CurrentModule = MathOptComplements
+```
+## Solving an equilibrium problem with MathOptComplements
+
+In this tutorial, we show how to solve an equilbrium problem
+using the different methods implemented in MathOptComplements.
+
+We take an example from the [JuMP documentation](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/complementarity/#Electricity-consumption),
+implementing a risk neutral competitive equilibrium between a producer and a consumer of electricity.
+
+```@example tutorial_equilibrium
+using JuMP
+
+function equilibrium_model()
+    I = 90_000                     # Annualized capital cost
+    C = 60                         # Operation cost per MWh
+    τ = 8_760                      # Hours per year
+    θ = [0.2, 0.2, 0.2, 0.2, 0.2]  # Scenario probabilities
+    A = [300, 350, 400, 450, 500]  # Utility function coefficients
+    B = 1                          # Utility function coefficients
+    model = Model()
+    @variable(model, x >= 0, start = 1)           # Installed capacity
+    @variable(model, Q[ω=1:5] >= 0, start = 1)  # Consumption
+    @variable(model, Y[ω=1:5] >= 0, start = 1)  # Production
+    @variable(model, P[ω=1:5], start = 1)       # Electricity price
+    @variable(model, μ[ω=1:5] >= 0, start = 1)  # Capital scarcity margin
+    # Unit investment cost equals annualized scarcity profit or investment is 0
+    @constraint(model, I - τ * θ' * μ ⟂ x)
+    # Difference between price and scarcity margin is equal to operation cost
+    @constraint(model, [ω = 1:5], C - (P[ω] - μ[ω]) ⟂ Y[ω])
+    # Price is equal to consumer's marginal utility
+    @constraint(model, [ω = 1:5], P[ω] - (A[ω] - B * Q[ω]) ⟂ Q[ω])
+    # Production is equal to consumption
+    @constraint(model, [ω = 1:5], Y[ω] - Q[ω] ⟂ P[ω])
+    # Production does not exceed capacity
+    @constraint(model, [ω = 1:5], x - Y[ω] ⟂ μ[ω])
+    return model
+end
+```
+
+This instance is featuring mixed-complementarity constraints, and as such
+is a good demo for MathOptComplements' capabilities.
+
+As a reference, we use the solution returned by the [PATH solver](https://pages.cs.wisc.edu/~ferris/path.html):
+```@example tutorial_equilibrium
+using PATHSolver
+model = equilibrium_model()
+JuMP.set_optimizer(model, PATHSolver.Optimizer)
+JuMP.optimize!(model)
+nothing
+
+```
+The solution returned by PATH is:
+```@example tutorial_equilibrium
+JuMP.value(model[:x]) # production in MWh
+```
+
+### Solution with a nonlinear solver
+
+We replace the solver PATH by Ipopt. MathOptComplements takes care of
+reformulating the problem automatically with appropriate nonlinear constraints.
+```@example tutorial_equilibrium
+using MathOptComplements
+using Ipopt
+model = equilibrium_model()
+MathOptComplements.Bridges.add_all_bridges(model)
+set_optimizer(model, Ipopt.Optimizer)
+JuMP.optimize!(model)
+nothing
+```
+The solution returned by Ipopt is:
+```@example tutorial_equilibrium
+JuMP.value(model[:x]) # production in MWh
+```
+
+
+

docs/src/index.md

@@ -0,0 +1,26 @@
+
+# MathOptComplements
+
+MathOptComplements is a [MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl/) extension for complementarity constraints.
+
+## Motivation
+
+MOI implements a set for mixed-complementarity constraints through [`MathOptInterface.Complements`](@extref).
+However, few solvers support mixed-complementarity out of the box, often requiring
+users to apply manual reformulations.
+
+MathOptComplements provides a systematic way to handle the complementarity
+constraints in JuMP and MathOptInterface by introducing a new set [`MathOptComplements.ComplementsWithSetType`](@ref).
+The package provides a rich collection of tools for manipulating complementarity constraints, including:
+- Equivalent **reformulations** in forms better suited to the target solver;
+- Automatic **relaxations** that reformulate complementarity constraints as nonlinear constraints.
+
+Under the hood, MathOptComplements extends [the bridge system](https://jump.dev/MathOptInterface.jl/stable/submodules/Bridges/overview/) implemented in MOI to optimally reformulate
+the complementarity constraints within the model.
+
+
+
+## Funding
+We acknowledge support from the [Fondation Mathématiques Jacques Hadamard](https://www.fondation-hadamard.fr/fr/)
+which has funded the PGMO-IROE project "A new optimization suite for large-scale market equilibrium".
+

docs/src/quickstart.md

@@ -0,0 +1,67 @@
+```@meta
+CurrentModule = MathOptComplements
+```
+## Quickstart
+
+The following code shows how to solve a simple Mathematical Program with Complementarity Constraints (MPCC) with Ipopt. This instance is a JuMP translation
+of `scholtes4.mod` in [MacMPEC](https://www.mcs.anl.gov/~leyffer/macmpec/comments.html).
+
+### Basic usage
+
+We start by writing the model with JuMP:
+
+```@example quickstart
+using JuMP
+z0 = [0, 1]
+model = Model()
+@variable(model, z[i=1:2] >= 0.0, start=z0[i])
+@variable(model, z3, start=0.0)
+@objective(model, Min, z[1] + z[2] - z3)
+@constraint(model, -4 * z[1] + z3 <= 0)
+@constraint(model, -4 * z[2] + z3 <= 0)
+@constraint(model, [z[1], z[2]] ∈ MOI.Complements(2))
+model
+```
+
+Solving this instance with Ipopt simply amounts to:
+```@example quickstart
+using MathOptComplements
+using Ipopt
+MathOptComplements.Bridges.add_all_bridges(model)
+set_optimizer(model, Ipopt.Optimizer)
+JuMP.optimize!(model)
+println("Solution: ", JuMP.value.(model[:z]))
+```
+Under the hood, MathOptComplements takes the complementarity
+constraints and reformulate it as a nonlinear constraint using
+the [`ScholtesRelaxation`](@ref) method.
+
+!!! note
+    We recommend setting the following options in Ipopt for optimal performance:
+    ```julia
+    JuMP.set_optimizer_attribute(model, "mu_strategy", "adaptive")
+    JuMP.set_optimizer_attribute(model, "bound_push", 1e-1)
+    JuMP.set_optimizer_attribute(model, "bound_relax_factor", 0.0)
+    ```
+
+
+### Changing the relaxation
+
+The [`ScholtesRelaxation`](@ref) is used by default, but the user
+has the freedom to use any of the relaxations implemented in MathOptComplements.
+Replacing the [`ScholtesRelaxation`](@ref) by the classical [`FischerBurmeisterRelaxation`](@ref) simply amounts to
+```@example quickstart
+JuMP.set_optimizer(model, () -> MathOptComplements.Optimizer(Ipopt.Optimizer()))
+MOI.set(model, MathOptComplements.DefaultComplementarityReformulation(), MathOptComplements.FischerBurmeisterRelaxation(1e-8))
+JuMP.optimize!(model)
+
+println("Solution: ", JuMP.value.(model[:z]))
+```
+
+Observe that the solution is here closer to the true solution `(0, 0)`
+than the solution returned by the Scholtes relaxation.
+
+!!! warning
+    MPCCs are nonconvex problems and they rarely have a unique solution.
+    In general, changing the relaxation method can yield a different local solution.
+

src/Bridges/Bridges.jl

@@ -20,15 +20,15 @@ include("to_sos1_bridge.jl")
     add_all_bridges(model::MOI.ModelLike, ::Type{T} = Float64)
 
 Add all `MathOptComplements` bridges to `model`. The model is typically a
-[`MOI.Bridges.LazyBridgeOptimizer`](@ref) so that the bridge graph is
+[`MathOptInterface.Bridges.LazyBridgeOptimizer`](@extref) so that the bridge graph is
 extended with the bridges needed to reformulate
-[`MathOptComplements.ComplementsWithSetType`](@ref) and [`MOI.Complements`](@ref)
+[`MathOptComplements.ComplementsWithSetType`](@ref) and [`MathOptInterface.Complements`](@extref)
 constraints.
 
 When used with a `LazyBridgeOptimizer`, the [`NonlinearBridge`](@ref) uses
 the default [`ScholtesRelaxation`](@ref) because the
 [`MathOptComplements.DefaultComplementarityReformulation`](@ref) optimizer
-attribute is only supported by [`MathOptComplements.Optimizer`](@ref).
+attribute is only supported by `MathOptComplements.Optimizer`.
 """
 function add_all_bridges(model::MOI.ModelLike, ::Type{T} = Float64) where {T}
     MOI.Bridges.add_bridge(model, SpecifySetTypeBridge{T})

src/Bridges/complements_vectorize_bridge.jl

@@ -16,17 +16,17 @@ where `T` is the coefficient type.
 
 `ComplementsVectorizeBridge` supports:
 
-  * [`MOI.VectorOfVariables`](@ref) in [`ComplementsWithSetType{S}`](@ref)
-    where `S` is [`MOI.GreaterThan{T}`](@ref), [`MOI.LessThan{T}`](@ref), or
-    [`MOI.EqualTo{T}`](@ref)
+  * [`MathOptInterface.VectorOfVariables`](@extref) in [`ComplementsWithSetType{S}`](@ref)
+    where `S` is [`MathOptInterface.GreaterThan`](@extref), [`MathOptInterface.LessThan`](@extref), or
+    [`MathOptInterface.EqualTo`](@extref)
 
 ## Target nodes
 
 `ComplementsVectorizeBridge` creates:
 
   * `F` in [`ComplementsWithSetType{SV}`](@ref), where `SV` is
-    [`MOI.Nonnegatives`](@ref), [`MOI.Nonpositives`](@ref), or
-    [`MOI.Zeros`](@ref) depending on the input set type
+    [`MathOptInterface.Nonnegatives`](@extref), [`MathOptInterface.Nonpositives`](@extref), or
+    [`MathOptInterface.Zeros`](@extref) depending on the input set type
 
 """
 struct ComplementsVectorizeBridge{T,F,S,SV} <: MOI.Bridges.Constraint.AbstractBridge

src/Bridges/nonlinear.jl

@@ -19,7 +19,7 @@ The relaxation method is determined by the
 
 `NonlinearBridge` supports:
 
-  * [`MOI.VectorOfVariables`](@ref) in [`ComplementsWithSetType{S}`](@ref)
+  * [`MathOptInterface.VectorOfVariables`](@extref) in [`ComplementsWithSetType{S}`](@ref)
 
 ## Target nodes
 

src/Bridges/specify_set_type_bridge.jl

@@ -3,7 +3,7 @@
 
 `SpecifySetTypeBridge` implements the following reformulation:
 
-  * `(x₁, x₂)` in [`MOI.Complements`](@ref) into `(x₁, x₂)` in
+  * `(x₁, x₂)` in [`MathOptInterface.Complements`](@extref) into `(x₁, x₂)` in
     [`ComplementsWithSetType{S}`](@ref)
 
 where `S` is determined by the bounds of `x₂`:
@@ -13,7 +13,7 @@ where `S` is determined by the bounds of `x₂`:
   * `x₂ ≤ 0` gives `S = MOI.Nonpositives`
   * `x₂ ≤ ub` (ub ≠ 0) gives `S = MOI.LessThan{T}`
   * `lb ≤ x₂ ≤ ub` gives `S = MOI.Interval{T}`
-  * `x₂` free gives `S = MOI.Zeros`
+  * `x₂` free gives `S = MOI.Real`
 
 The bridge also adds the appropriate bound on the activity variable `x₁`
 (for example, `x₁ ≥ 0` when `x₂` has a lower bound).
@@ -22,15 +22,15 @@ The bridge also adds the appropriate bound on the activity variable `x₁`
 
 `SpecifySetTypeBridge` supports:
 
-  * [`MOI.VectorOfVariables`](@ref) in [`MOI.Complements`](@ref)
+  * [`MathOptInterface.VectorOfVariables`](@extref) in [`MathOptInterface.Complements`](@extref)
 
 ## Target nodes
 
 `SpecifySetTypeBridge` creates:
 
-  * [`MOI.VectorOfVariables`](@ref) in [`ComplementsWithSetType{S}`](@ref)
-  * [`MOI.VariableIndex`](@ref) in [`MOI.GreaterThan{T}`](@ref) or
-    [`MOI.LessThan{T}`](@ref) (bounds on `x₁`)
+  * [`MathOptInterface.VectorOfVariables`](@extref) in [`ComplementsWithSetType{S}`](@ref)
+  * [`MathOptInterface.VariableIndex`](@extref) in [`MathOptInterface.GreaterThan`](@extref) or
+    [`MathOptInterface.LessThan`](@extref) (bounds on `x₁`)
 
 """
 mutable struct SpecifySetTypeBridge{T} <: MOI.Bridges.Constraint.AbstractBridge

src/Bridges/split_interval_bridge.jl

@@ -18,18 +18,18 @@ must be a variable.
 
 `SplitIntervalBridge` supports:
 
-  * [`MOI.AbstractVectorFunction`](@ref) in
+  * [`MathOptInterface.AbstractVectorFunction`](@extref) in
     [`ComplementsWithSetType{MOI.Interval{T}}`](@ref)
 
 ## Target nodes
 
 `SplitIntervalBridge` creates:
 
-  * [`MOI.VectorOfVariables`](@ref) in
+  * [`MathOptInterface.VectorOfVariables`](@extref) in
     [`ComplementsWithSetType{MOI.GreaterThan{T}}`](@ref)
-  * [`MOI.VectorOfVariables`](@ref) in
+  * [`MathOptInterface.VectorOfVariables`](@extref) in
     [`ComplementsWithSetType{MOI.LessThan{T}}`](@ref)
-  * `G` in [`MOI.EqualTo{T}`](@ref) (the splitting equality)
+  * `G` in [`MathOptInterface.EqualTo`](@extref) (the splitting equality)
 
 where `G` is the scalar function type of the first component.