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 JuMP, and then by specifying the derivatives manually in NLPModels.
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 has one bound constraint ($x_1 \leq 0.5$) and two generic constraints. We display in the following picture the feasible set of HS15. The problem has two local solutions: $(x_1, x_2) = (0.5, 2)$ and $(x_1, x_2) \approx (-0.792, -1.262)$.
Using MadNLP with JuMP
JuMP is the easiest way to implement a nonlinear problem. In JuMP, the user just has only 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 * ((-x1² + x2) ^ 2)) + (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.9.1, running with MUMPS v5.8.2
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 1.11e-16 1.19e-13 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......: 1.1912602012132843e-13 1.1912602012132843e-13
Constraint violation....: 1.1102230246251565e-16 1.1102230246251565e-16
Complementarity.........: 1.5755252336099267e-09 1.5755252336099267e-09
Overall NLP error.......: 1.5755252336099267e-09 1.5755252336099267e-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
Number of KKT factorizations = 24
Number of KKT backsolves = 25
Total wall secs in initialization = 2.393 s
Total wall secs in linear solver = 0.003 s
Total wall secs in NLP function evaluations = 8.973 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 11.410 s
Total wall secs = 22.779 s
EXIT: Optimal Solution Found (tol = 1.0e-08).Under the hood, JuMP builds a sparse AD backend to evaluate the first and second-order derivatives of the objective and the constraints. Then, MadNLP takes as input the callbacks generated by JuMP and wraps them inside a MadNLP.MOIModel.
Using MadNLP with NLPModels
Alternatively, we can compute the derivatives manually in an AbstractNLPModel. This second option, although more complicated, gives us more flexibility and comes without boilerplate.
We define the HS15Model 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 an 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 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
endThe Jacobian is defined as
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:
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, 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)
nothingThis is MadNLP version v0.9.1, running with MUMPS v5.8.2
Number of nonzeros in constraint Jacobian............: 4
Number of nonzeros in Lagrangian Hessian.............: 3
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 1.11e-16 1.19e-13 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......: 1.1912602012132843e-13 1.1912602012132843e-13
Constraint violation....: 1.1102230246251565e-16 1.1102230246251565e-16
Complementarity.........: 1.5755252336101401e-09 1.5755252336101401e-09
Overall NLP error.......: 1.5755252336101401e-09 1.5755252336101401e-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
Number of KKT factorizations = 24
Number of KKT backsolves = 25
Total wall secs in initialization = 2.422 s
Total wall secs in linear solver = 0.002 s
Total wall secs in NLP function evaluations = 0.000 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 7.724 s
Total wall secs = 10.149 s
EXIT: 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-9If we come back at the picture displayed at the beginning of this page, the solution found is the one at the bottom left. At this solution, only one constraint is active, explaining why only the first element in results.multipliers is non-null (up to floating point approximation).
The dual multipliers associated to the variable's lower bounds (resp. upper bounds) are stored in results.multipliers_L (resp. results.multipliers_U):
results.multipliers_L2-element Vector{Float64}:
0.0
0.0By changing the initial point in HS15Model, MadNLP converges to the second local solution of the non-convex problem HS15:
x0 = ones(2)
nlp = HS15Model(x0)
results = madnlp(nlp)
nothingThis is MadNLP version v0.9.1, running with MUMPS v5.8.2
Number of nonzeros in constraint Jacobian............: 4
Number of nonzeros in Lagrangian Hessian.............: 3
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 5.8004900e+01 5.20e-01 8.04e+01 1.49e+00 -1.0 - 0.00e+00 1 0
1 5.5900651e+01 5.15e-01 7.60e+01 1.43e+00 -1.0 - 9.17e-03 1 1h
2 5.3726392e+01 5.10e-01 8.07e+01 2.01e+00 -1.0 - 1.04e-02 1 1h
3 5.3797218e+01 5.09e-01 1.75e+02 1.85e-02 -1.0 - 1.25e-03 1 1h
4 3.0774740e+02 2.03e-03 1.67e+02 8.90e-02 -1.0 - 1.00e+00 1 1H
5 3.0672128e+02 9.77e-04 3.22e+01 2.61e-02 -1.0 - 5.21e-01 1 1h
6 3.0680411e+02 2.38e-08 2.45e-03 1.86e-02 -1.0 - 1.00e+00 1 1h
7 3.0650084e+02 6.07e-07 8.84e-05 5.88e-05 -3.8 - 9.99e-01 1 1h
8 3.0649998e+02 4.73e-12 1.41e-10 3.43e-07 -5.7 - 1.00e+00 1 1h
9 3.0649998e+02 8.88e-16 5.48e-14 1.69e-10 -9.0 - 1.00e+00 2 1h
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 2.0167125901139678e+02 3.0649997549276407e+02
Dual infeasibility......: 5.4833787889654930e-14 8.3336389759955301e-14
Constraint violation....: 8.8817841970012523e-16 8.8817841970012523e-16
Complementarity.........: 1.1127465637341613e-10 1.6911522057533453e-10
Overall NLP error.......: 1.6911522057533453e-10 1.6911522057533453e-10
Number of objective function evaluations = 11
Number of objective gradient evaluations = 10
Number of constraint evaluations = 11
Number of constraint Jacobian evaluations = 10
Number of Lagrangian Hessian evaluations = 9
Number of KKT factorizations = 10
Number of KKT backsolves = 12
Total wall secs in initialization = 0.001 s
Total wall secs in linear solver = 0.001 s
Total wall secs in NLP function evaluations = 0.000 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 0.001 s
Total wall secs = 0.003 s
EXIT: Optimal Solution Found (tol = 1.0e-08).results.solution2-element Vector{Float64}:
0.500000009999211
1.9999999400071047Note that we don't know beforehand towards each solution MadNLP is converging.