Nonlinear Modeling
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 @expression
s.
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 function | Julia function |
---|---|
op_ifelse | ifelse |
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.
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^2
square (generic function with 1 method)
julia> f(x, y) = (x - 1)^2 + (y - 2)^2
f (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:
- The model to which the operator is added.
- 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.
- The number of scalar input arguments that the function takes.
- A Julia method which computes the function.
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^2
square (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^2
square (generic function with 1 method)
julia> f(x, y) = (x - 1)^2 + (y - 2)^2
f (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_square
NonlinearOperator(square, :op_square)
julia> op_f = add_nonlinear_operator(model, 2, f; name = :op_f)
NonlinearOperator(f, :op_f)
julia> model[:op_f] = op_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]))
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^2
f (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 optional
NonlinearOperator(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)^2
f (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 optional
NonlinearOperator(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.
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
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