Nonlinear Modeling

Warning

This page describes a new nonlinear interface to JuMP. It replaces the legacy @NL interface, which is documented at Nonlinear Modeling (Legacy). The API described below is stable, and it will not break with future 1.X releases of JuMP. However, solver support may be limited, and there may be gaps in functionality compared with the legacy interface. To report a bug, or request a missing feature, please open an issue.

JuMP has support for nonlinear (convex and nonconvex) optimization problems. JuMP is able to automatically provide exact, sparse second-order derivatives to solvers. This information can improve solver accuracy and performance.

Set a nonlinear objective

Use @objective to set a nonlinear objective.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @objective(model, Min, exp(x[1]) - sqrt(x[2]))
exp(x[1]) - sqrt(x[2])

To modify a nonlinear objective, call @objective again.

Add a nonlinear constraint

Use @constraint to add a nonlinear constraint.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @constraint(model, exp(x[1]) <= 1)
exp(x[1]) - 1.0 ≤ 0

julia> @constraint(model, con[i = 1:2], 2^x[i] >= i)
2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarNonlinearFunction, MathOptInterface.GreaterThan{Float64}}, ScalarShape}}:
 con[1] : (2.0 ^ x[1]) - 1.0 ≥ 0
 con[2] : (2.0 ^ x[2]) - 2.0 ≥ 0

Delete a nonlinear constraint using delete:

julia> delete(model, con[1])

Create a nonlinear expression

Use @expression to create nonlinear expression objects:

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> expr = @expression(model, exp(x[1]) + sqrt(x[2]))
exp(x[1]) + sqrt(x[2])

julia> my_anon_expr = @expression(model, [i = 1:2], sin(x[i]))
2-element Vector{NonlinearExpr}:
 sin(x[1])
 sin(x[2])

julia> @expression(model, my_expr[i = 1:2], sin(x[i]))
2-element Vector{NonlinearExpr}:
 sin(x[1])
 sin(x[2])

A NonlinearExpr can be used in @objective, @constraint, and even nested in other @expressions.

julia> @objective(model, Min, expr^2 + 1)
((exp(x[1]) + sqrt(x[2])) ^ 2.0) + 1.0

julia> @constraint(model, [i = 1:2], my_expr[i] <= i)
2-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarNonlinearFunction, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
 sin(x[1]) - 1.0 ≤ 0
 sin(x[2]) - 2.0 ≤ 0

julia> @expression(model, nested[i = 1:2], sin(my_expr[i]))
2-element Vector{NonlinearExpr}:
 sin(sin(x[1]))
 sin(sin(x[2]))

Use value to query the value of a nonlinear expression:

julia> set_start_value(x[1], 1.0)

julia> value(start_value, nested[1])
0.7456241416655579

julia> sin(sin(1.0))
0.7456241416655579

Nonlinear expressions in detail

Nonlinear expressions in JuMP are represented by a NonlinearExpr object.

Constructors

Nonlinear expressions can be created using the NonlinearExpr constructors:

julia> model = Model();

julia> @variable(model, x);

julia> expr = NonlinearExpr(:sin, Any[x])
sin(x)

or via operator overloading:

julia> model = Model();

julia> @variable(model, x);

julia> expr = sin(x)
sin(x)

Supported arguments

Nonlinear expressions can contain a mix of numbers, AffExpr, QuadExpr, and other NonlinearExpr:

julia> model = Model();

julia> @variable(model, x);

julia> aff = x + 1;

julia> quad = x^2 + x;

julia> expr = cos(x) * sin(quad) + aff
(cos(x) * sin(x² + x)) + (x + 1)

Supported operators

The list of supported operators may vary between solvers. Given an optimizer, query the list of supported operators using MOI.ListOfSupportedNonlinearOperators:

julia> import Ipopt

julia> import MathOptInterface as MOI

julia> MOI.get(Ipopt.Optimizer(), MOI.ListOfSupportedNonlinearOperators())
85-element Vector{Symbol}:
 :+
 :-
 :abs
 :sqrt
 :cbrt
 :abs2
 :inv
 :log
 :log10
 :log2
 ⋮
 :min
 :max
 :&&
 :||
 :<=
 :(==)
 :>=
 :<
 :>

In some univariate cases, the operator is defined in SpecialFunctions.jl. To use these functions, you must explicitly import SpecialFunctions.jl

julia> import Ipopt

julia> op = MOI.get(Ipopt.Optimizer(), MOI.ListOfSupportedNonlinearOperators());

julia> :erfcx in op
true

julia> :dawson in op
true

julia> import SpecialFunctions

julia> model = Model();

julia> @variable(model, x)
x

julia> @expression(model, SpecialFunctions.erfcx(x))
erfcx(x)

julia> @expression(model, SpecialFunctions.dawson(x))
dawson(x)

Limitations

Some nonlinear expressions cannot be created via operator overloading. For example, to minimize the likelihood of bugs in user-code, we have not overloaded comparisons such as < and >= between JuMP objects:

julia> model = Model();

julia> @variable(model, x);

julia> x < 1
ERROR: Cannot evaluate `<` between a variable and a number.
[...]

Instead, wrap the expression in the @expression macro:

julia> model = Model();

julia> @variable(model, x);

julia> expr = @expression(model, x < 1)
x < 1

For technical reasons, other operators that are not overloaded include ||, &&, and ifelse.

julia> model = Model();

julia> @variable(model, x);

julia> expr = @expression(model, ifelse(x < -1 || x >= 1, x^2, 0.0))
ifelse((x < -1) || (x >= 1), x², 0.0)

As an alternative, use the JuMP.op_ functions, which fallback to the various comparison and logical operators:

julia> model = Model();

julia> @variable(model, x);

julia> expr = op_ifelse(
           op_or(op_strictly_less_than(x, -1), op_greater_than_or_equal_to(x, 1)),
           x^2,
           0.0,
       )
ifelse((x < -1) || (x >= 1), x², 0.0)

The available functions are:

JuMP functionJulia function
op_ifelseifelse
op_and&&
op_or||
op_greater_than_or_equal_to>=
op_less_than_or_equal_to<=
op_equal_to==
op_strictly_greater_than>
op_strictly_less_than<

Fields

Each NonlinearExpr has two fields.

The .head field is a Symbol that represents the operator being called:

julia> expr.head
:sin

The .args field is a Vector{Any} containing the arguments to the operator:

julia> expr.args
1-element Vector{Any}:
 x

User-defined operators

In addition to a standard list of univariate and multivariate operators recognized by the MOI.Nonlinear submodule, JuMP supports user-defined Julia operators.

Warning

User-defined operators must return a scalar output. For a work-around, see User-defined operators with vector outputs.

Add an operator

Add a user-defined operator using the @operator macro:

julia> using JuMP
julia> square(x) = x^2square (generic function with 1 method)
julia> f(x, y) = (x - 1)^2 + (y - 2)^2f (generic function with 1 method)
julia> model = Model();
julia> @operator(model, op_square, 1, square)NonlinearOperator(square, :op_square)
julia> @operator(model, op_f, 2, f)NonlinearOperator(f, :op_f)
julia> @variable(model, x[1:2]);
julia> @objective(model, Min, op_f(x[1], op_square(x[2])))op_f(x[1], op_square(x[2]))

The arguments to @operator are:

  1. The model to which the operator is added.
  2. A Julia symbol object which serves as the name of the user-defined operator in JuMP expressions. This name must not be the same as that of the function.
  3. The number of scalar input arguments that the function takes.
  4. A Julia method which computes the function.
Warning

User-defined operators cannot be deleted.

You can obtain a reference to the operator using the model[:key] syntax:

julia> using JuMP
julia> square(x) = x^2square (generic function with 1 method)
julia> model = Model();
julia> @operator(model, op_square, 1, square)NonlinearOperator(square, :op_square)
julia> op_square_2 = model[:op_square]NonlinearOperator(square, :op_square)

Add an operator without macros

The @operator macro is syntactic sugar for add_nonlinear_operator. Thus, the non-macro version of the preceding example is:

julia> using JuMP
julia> square(x) = x^2square (generic function with 1 method)
julia> f(x, y) = (x - 1)^2 + (y - 2)^2f (generic function with 1 method)
julia> model = Model();
julia> op_square = add_nonlinear_operator(model, 1, square; name = :op_square)NonlinearOperator(square, :op_square)
julia> model[:op_square] = op_squareNonlinearOperator(square, :op_square)
julia> op_f = add_nonlinear_operator(model, 2, f; name = :op_f)NonlinearOperator(f, :op_f)
julia> model[:op_f] = op_fNonlinearOperator(f, :op_f)
julia> @variable(model, x[1:2]);
julia> @objective(model, Min, op_f(x[1], op_square(x[2])))op_f(x[1], op_square(x[2]))

Operators with the same name as an existing function

A common error encountered is the following:

julia> using JuMP

julia> model = Model();

julia> f(x) = x^2
f (generic function with 1 method)

julia> @operator(model, f, 1, f)
ERROR: Unable to add the nonlinear operator `:f` with the same name as
an existing function.
[...]

This error occurs because @operator(model, f, 1, f) is equivalent to:

julia> f = add_nonlinear_operator(model, 1, f; name = :f)

but f already exists as a Julia function.

If you evaluate the function without adding it as an operator, JuMP will trace the function using operator overloading:

julia> @variable(model, x);

julia> f(x)
x²

To force JuMP to treat f as a user-defined operator and not trace it, add the operator using add_nonlinear_operator and define a new method which manually creates a NonlinearExpr:

julia> _ = add_nonlinear_operator(model, 1, f; name = :f)
NonlinearOperator(f, :f)

julia> f(x::AbstractJuMPScalar) = NonlinearExpr(:f, Any[x])
f (generic function with 2 methods)

julia> @expression(model, log(f(x)))
log(f(x))

Gradients and Hessians

By default, JuMP will use automatic differentiation to compute the gradient and Hessian of user-defined operators. If your function is not amenable to automatic differentiation, or you can compute analytic derivatives, you may pass additional arguments to @operator to compute the first- and second-derivatives.

Univariate functions

For univariate functions, a gradient function ∇f returns a number that represents the first-order derivative. You may, in addition, pass a third function which returns a number representing the second-order derivative:

julia> using JuMP
julia> f(x) = x^2f (generic function with 1 method)
julia> ∇f(x) = 2x∇f (generic function with 1 method)
julia> ∇²f(x) = 2∇²f (generic function with 1 method)
julia> model = Model();
julia> @operator(model, op_f, 1, f, ∇f, ∇²f) # Providing ∇²f is optionalNonlinearOperator(f, :op_f)
julia> @variable(model, x)x
julia> @objective(model, Min, op_f(x))op_f(x)

Multivariate functions

For multivariate functions, the gradient function ∇f must take an AbstractVector as the first argument that is filled in-place. The Hessian function, ∇²f, must take an AbstractMatrix as the first argument, the lower-triangular of which is filled in-place:

julia> using JuMP
julia> f(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2f (generic function with 1 method)
julia> function ∇f(g::AbstractVector{T}, x::T...) where {T} g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2 g[2] = 200 * (x[2] - x[1]^2) return end∇f (generic function with 1 method)
julia> function ∇²f(H::AbstractMatrix{T}, x::T...) where {T} H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2 # H[1, 2] = -400 * x[1] <-- Not needed. Fill the lower-triangular only. H[2, 1] = -400 * x[1] H[2, 2] = 200.0 return end∇²f (generic function with 1 method)
julia> model = Model();
julia> @operator(model, rosenbrock, 2, f, ∇f, ∇²f) # Providing ∇²f is optionalNonlinearOperator(f, :rosenbrock)
julia> @variable(model, x[1:2])2-element Vector{VariableRef}: x[1] x[2]
julia> @objective(model, Min, rosenbrock(x[1], x[2]))rosenbrock(x[1], x[2])

You may assume the Hessian matrix H is initialized with zeros, and because H is symmetric, you need only to fill in the non-zero lower-triangular terms. The matrix type passed in as H depends on the automatic differentiation system, so make sure the first argument to the Hessian function supports an AbstractMatrix (it may be something other than Matrix{Float64}). Moreover, you may assume only that H supports size(H) and setindex!. Finally, the matrix is treated as dense, so the performance will be poor on functions with high-dimensional input.

User-defined operators with vector inputs

User-defined operators which take vectors as input arguments (for example, f(x::Vector)) are not supported. Instead, use Julia's splatting syntax to create a function with scalar arguments. For example, instead of:

f(x::Vector) = sum(x[i]^i for i in 1:length(x))

define:

f(x...) = sum(x[i]^i for i in 1:length(x))

Another approach is to define the splatted function as an anonymous function:

julia> using JuMP
julia> model = Model();
julia> @variable(model, x[1:5])5-element Vector{VariableRef}: x[1] x[2] x[3] x[4] x[5]
julia> f(x::Vector) = sum(x[i]^i for i in 1:length(x))f (generic function with 1 method)
julia> @operator(model, op_f, 5, (x...) -> f(collect(x)))NonlinearOperator(#6, :op_f)
julia> @objective(model, Min, op_f(x...))op_f(x[1], x[2], x[3], x[4], x[5])

Automatic differentiation

JuMP does not support black-box optimization, so all user-defined operators must provide derivatives in some form. Fortunately, JuMP supports automatic differentiation of user-defined operators.

Info

Automatic differentiation is not finite differencing. JuMP's automatically computed derivatives are not subject to approximation error.

JuMP uses ForwardDiff.jl to perform automatic differentiation of user-defined operators; see the ForwardDiff.jl documentation for a description of how to write a function suitable for automatic differentiation.

Common mistakes when writing a user-defined operator

Warning

Get an error like No method matching Float64(::ForwardDiff.Dual)? Read this section, and see the guidelines at ForwardDiff.jl.

The most common error is that your user-defined operator is not generic with respect to the number type, that is, don't assume that the input to the function is Float64.

f(x::Float64) = 2 * x  # This will not work.
f(x::Real)    = 2 * x  # This is good.
f(x)          = 2 * x  # This is also good.

Another reason you may encounter this error is if you create arrays inside your function which are Float64.

function bad_f(x...)
    y = zeros(length(x))  # This constructs an array of `Float64`!
    for i in 1:length(x)
        y[i] = x[i]^i
    end
    return sum(y)
end

function good_f(x::T...) where {T<:Real}
    y = zeros(T, length(x))  # Construct an array of type `T` instead!
    for i in 1:length(x)
        y[i] = x[i]^i
    end
    return sum(y)
end