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