Nonlinear Modeling

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:

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.

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.

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.

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

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