CentralPath.jl

This package uses JuMP, PolyJuMP and HomotopyContinuation.jl to compute the central paths associated to a given algebraic optimization problem.

The algorithm proceeds as follows:

• it transcribes the central path equations as a semi-algebraic set, parameterized by a barrier parameter mu;
• the semi-algebraic set is passed to HomotopyContinuation.jl, which computes the set of solutions for multiple barrier parameters using homotopy continuation.

The algorithm works only for very small instances, and has mostly a pedagogical purpose. Furthermore, the algorithm works only if the KKT conditions of the problem writes as a system of polynomial equations defining a semi-algebraic set. This is the case for LPs, QPs, and NLPs formulated with polynomial expressions.

Examples

LP

We take the LP from this article (Figure 2), whose central path exhibits a very sharp turn close to the solution.

using JuMP
using PolyJuMP
using HomotopyContinuation
using CentralPath

model = Model(optimizer_with_attributes(
    PolyJuMP.KKT.Optimizer,
    "solver" => HomotopyContinuation.SemialgebraicSetsHCSolver(),
))

# Define LP with JuMP
eps = 3e-1
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@constraint(model, x + 2*y >= eps)
@objective(model, Min, 2*x + 5*y)

# Barrier values
mu_vals = [10^f for f in -9:0.1:5]
# Solve with homotopy continuation
MOI.Utilities.attach_optimizer(model)
solver = JuMP.unsafe_backend(model)
results = CentralPath.compute_central_path(solver, mu_vals)

NLP

We consider the instance HS15 from the Hock & Schittkowski collection. The problem is non-convex but is semi-algebraic. It has two distinct central paths.

using JuMP
using PolyJuMP
using HomotopyContinuation
using CentralPath

model = Model(optimizer_with_attributes(
    PolyJuMP.KKT.Optimizer,
    "solver" => HomotopyContinuation.SemialgebraicSetsHCSolver(),
))
# Define non-convex NLP with JuMP
@variable(model, x1 <= 0.5)
@variable(model, x2)
@objective(model, Min, 100.0 * (x2 - x1^2)^2 + (1.0 - x1)^2)
@constraint(model, x1 * x2 >= 1.0)
@constraint(model, x1 + x2^2 >= 0.0)

# Barrier values
mu_vals = [10^f for f in -9:0.1:9]
# Solve with homotopy continuation
MOI.Utilities.attach_optimizer(model)
solver = JuMP.unsafe_backend(model)
results = CentralPath.compute_central_path(solver, mu_vals)