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)

@NLobjective(model, Min, 100.0 * (x2 - x1^2)^2 + (1.0 - x1)^2)
@NLconstraint(model, x1 * x2 >= 1.0)
@NLconstraint(model, x1 + x2^2 >= 0.0)

Then, solving HS15 with MadNLP directly translates to

JuMP.set_optimizer(model, MadNLP.Optimizer)
JuMP.optimize!(model)

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()
    )
end

This 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
end

The 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
end

We define similarly the constraints

function NLPModels.cons!(nlp::HS15Model, x::AbstractVector, c::AbstractVector)
    c[1] = x[1] * x[2]
    c[2] = x[1] + x[2]^2
end

and 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)
end

It 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])
end

Now 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
end

Once 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)
ips = MadNLP.InteriorPointSolver(nlp)
MadNLP.optimize!(ips)

MadNLP converges in 19 iterations to a (local) optimal solution. We can query the primal and the dual solutions respectively by

ips.x

and

ips.l