Quickstart
This page presents a quickstart guide to solve a nonlinear problem with MadNLP.
As a demonstration, we show how to implement the HS15 nonlinear problem from the Hock & Schittkowski collection, first by using a nonlinear modeler and then by specifying the derivatives manually.
The HS15 problem is defined as:
\[\begin{aligned} \min_{x_1, x_2} \; & 100 \times (x_2 - x_1^2)^2 +(1 - x_1)^2 \\ \text{subject to} \quad & x_1 \times x_2 \geq 1 \\ & x_1 + x_2^2 \geq 0 \\ & x_1 \leq 0.5 \end{aligned} \]
Despite its small dimension, its resolution remains challenging as the problem is nonlinear nonconvex. Note that HS15 encompasses one bound constraint ($x_1 \leq 0.5$) and two generic constraints.
Using a nonlinear modeler: JuMP.jl
The easiest way to implement a nonlinear problem is to use a nonlinear modeler as JuMP. In JuMP, the user just has to pass the structure of the problem, the computation of the first- and second-order derivatives being handled automatically.
Using JuMP's syntax, the HS15 problem translates to
using JuMP
model = Model()
@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)
println(model)Min (100.0 * ((-x1² + x2) ^ 2.0)) + (x1² - 2 x1 + 1)
Subject to
x1*x2 ≥ 1
x2² + x1 ≥ 0
x1 ≤ 0.5Then, solving HS15 with MadNLP directly translates to
using MadNLP
JuMP.set_optimizer(model, MadNLP.Optimizer)
JuMP.optimize!(model)This is MadNLP version v0.8.12, running with MUMPS v5.8.1
Number of nonzeros in constraint Jacobian............: 4
Number of nonzeros in Lagrangian Hessian.............: 5
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 1
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 2
inequality constraints with only lower bounds: 2
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
0 1.0000000e+00 1.01e+00 1.00e+00 5.00e-01 -1.0 - 0.00e+00 1 0
1 9.9758855e-01 1.00e+00 4.61e+01 3.80e-01 -1.0 - 9.80e-03 1 1h
2 9.9664309e-01 1.00e+00 5.00e+02 2.35e-03 -1.0 - 9.93e-05 1 1h
3 1.3615174e+00 9.99e-01 4.41e+02 2.44e-03 -1.0 - 4.71e-04 1 1H
4 1.3742697e+00 9.99e-01 3.59e+02 1.40e-01 -1.0 - 2.68e-05 1 1h
5 1.4692139e+00 9.99e-01 4.94e+02 2.31e+00 -1.0 - 1.46e-04 1 1h
6 3.1727722e+00 9.97e-01 3.76e+02 2.36e+00 -1.0 - 9.77e-04 1 11h
7 3.1746542e+00 9.96e-01 1.53e+02 3.13e-03 -1.0 4.0 7.98e-04 1 1h
8 8.2322585e+00 9.85e-01 2.47e+02 5.06e-03 -1.0 - 7.81e-03 1 8h
9 8.2886187e+00 9.84e-01 4.78e+02 1.51e-04 -1.0 - 2.48e-04 1 1h
iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
10 4.0285541e+01 8.72e-01 4.57e+02 6.32e-04 -1.0 - 6.25e-02 1 5h
11 2.8602305e+02 2.66e-01 4.91e+02 5.99e-04 -1.0 - 5.00e-01 1 2h
12 3.9870638e+02 6.41e-03 4.21e+02 3.32e-03 -1.0 - 1.00e+00 1 1h
13 3.4928112e+02 2.93e-02 2.62e+01 1.02e+00 -1.0 - 1.00e+00 1 1h
14 3.5909050e+02 6.98e-03 1.66e+00 6.15e-02 -1.0 - 1.00e+00 1 1h
15 3.6047026e+02 6.05e-05 1.44e-02 6.46e-02 -1.0 - 1.00e+00 1 1h
16 3.6038251e+02 2.50e-07 7.99e-05 1.89e-03 -2.5 - 1.00e+00 1 1h
17 3.6037976e+02 2.63e-10 7.93e-08 1.27e-06 -5.7 - 1.00e+00 1 1h
18 3.6037976e+02 2.22e-16 0.00e+00 1.58e-09 -8.6 - 1.00e+00 2 1h
Number of Iterations....: 18
(scaled) (unscaled)
Objective...............: 3.6037976240508465e+02 3.6037976240508465e+02
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16
Complementarity.........: 1.5755252336095001e-09 1.5755252336095001e-09
Overall NLP error.......: 1.5755252336095001e-09 1.5755252336095001e-09
Number of objective function evaluations = 46
Number of objective gradient evaluations = 19
Number of constraint evaluations = 46
Number of constraint Jacobian evaluations = 19
Number of Lagrangian Hessian evaluations = 18
Total wall secs in initialization = 0.726
Total wall secs in linear solver = 0.002
Total wall secs in NLP function evaluations = 4.558
Total wall secs in solver (w/o init./fun./lin. alg.) = 6.667
Total wall secs = 11.953
EXIT: Optimal Solution Found (tol = 1.0e-08).Under the hood, JuMP builds a nonlinear model with a sparse AD backend to evaluate the first and second-order derivatives of the objective and the constraints. Internally, MadNLP takes as input the callbacks generated by JuMP and wraps them inside a MadNLP.MOIModel.
Specifying the derivatives by hand: NLPModels.jl
Alternatively, we can compute the derivatives manually and define directly a NLPModel associated to our problem. This second option, although more complicated, gives us more flexibility and comes without boilerplate.
We define a new NLPModel structure simply as:
struct HS15Model <: NLPModels.AbstractNLPModel{Float64,Vector{Float64}}
meta::NLPModels.NLPModelMeta{Float64, Vector{Float64}}
counters::NLPModels.Counters
end
function HS15Model(x0)
return HS15Model(
NLPModels.NLPModelMeta(
2, #nvar
ncon = 2,
nnzj = 4,
nnzh = 3,
x0 = x0,
y0 = zeros(2),
lvar = [-Inf, -Inf],
uvar = [0.5, Inf],
lcon = [1.0, 0.0],
ucon = [Inf, Inf],
minimize = true
),
NLPModels.Counters()
)
endMain.HS15ModelThis structure takes as input the initial position x0 and generates a AbstractNLPModel. NLPModelMeta stores the information about the structure of the problem (variables and constraints' lower and upper bounds, number of variables, number of constraints, ...). Counters is just a utility storing the number of time each callbacks is being called.
The objective function takes as input a HS15Model instance and a vector with dimension 2 storing the current values for $x_1$ and $x_2$:
function NLPModels.obj(nlp::HS15Model, x::AbstractVector)
return 100.0 * (x[2] - x[1]^2)^2 + (1.0 - x[1])^2
endThe corresponding gradient writes (note that we update the values of the gradient g inplace):
function NLPModels.grad!(nlp::HS15Model, x::AbstractVector, g::AbstractVector)
z = x[2] - x[1]^2
g[1] = -400.0 * z * x[1] - 2.0 * (1.0 - x[1])
g[2] = 200.0 * z
return g
endSimilarly, we define the constraints
function NLPModels.cons!(nlp::HS15Model, x::AbstractVector, c::AbstractVector)
c[1] = x[1] * x[2]
c[2] = x[1] + x[2]^2
return c
endand the associated Jacobian
function NLPModels.jac_structure!(nlp::HS15Model, I::AbstractVector{T}, J::AbstractVector{T}) where T
copyto!(I, [1, 1, 2, 2])
copyto!(J, [1, 2, 1, 2])
end
function NLPModels.jac_coord!(nlp::HS15Model, x::AbstractVector, J::AbstractVector)
J[1] = x[2] # (1, 1)
J[2] = x[1] # (1, 2)
J[3] = 1.0 # (2, 1)
J[4] = 2*x[2] # (2, 2)
return J
endIt remains to implement the Hessian of the Lagrangian for a HS15Model. The Lagrangian of the problem writes
\[L(x_1, x_2, y_1, y_2) = 100 \times (x_2 - x_1^2)^2 +(1 - x_1)^2 + y_1 \times (x_1 \times x_2) + y_2 \times (x_1 + x_2^2)\]
and we aim at evaluating its second-order derivative $\nabla^2_{xx}L(x_1, x_2, y_1, y_2)$.
We first have to define the sparsity structure of the Hessian, which is assumed to be sparse. The Hessian is a symmetric matrix, and by convention we pass only the lower-triangular part of the matrix to the solver. Hence, we define the sparsity structure as
function NLPModels.hess_structure!(nlp::HS15Model, I::AbstractVector{T}, J::AbstractVector{T}) where T
copyto!(I, [1, 2, 2])
copyto!(J, [1, 1, 2])
endNow that the sparsity structure is defined, the associated Hessian writes down:
function NLPModels.hess_coord!(nlp::HS15Model, x, y, H::AbstractVector; obj_weight=1.0)
# Objective
H[1] = obj_weight * (-400.0 * x[2] + 1200.0 * x[1]^2 + 2.0)
H[2] = obj_weight * (-400.0 * x[1])
H[3] = obj_weight * 200.0
# First constraint
H[2] += y[1] * 1.0
# Second constraint
H[3] += y[2] * 2.0
return H
endOnce the problem specified in NLPModels's syntax, we can create a new MadNLP instance and solve it:
x0 = zeros(2) # initial position
nlp = HS15Model(x0)
solver = MadNLP.MadNLPSolver(nlp; print_level=MadNLP.INFO)
results = MadNLP.solve!(solver)"Execution stats: Optimal Solution Found (tol = 1.0e-08)."MadNLP converges in 19 iterations to a (local) optimal solution. MadNLP returns a MadNLPExecutionStats storing all the results. We can query the primal and the dual solutions respectively by
results.solution2-element Vector{Float64}:
-0.7921232178470456
-1.2624298435831804and
results.multipliers2-element Vector{Float64}:
-477.17046873911966
-3.1262000119573344e-9