JuMP 1.15.0 is released
releases ·We are happy to announce the release of JuMP 1.15.0.
This is a very large minor release because it adds an entirely new data structure and API path for working with nonlinear programs.
The previous nonlinear interface remains unchanged and is documented at Nonlinear Modeling (Legacy). The new interface is a treated as a non-breaking feature addition and is documented at Nonlinear Modeling.
Background
The “legacy” nonlinear interface was some of the oldest code in JuMP. When we created MathOptInterface (and rewrote nearly all of JuMP), we explicitly decided to leave the nonlinear interface as it was, and “come back to it later.”
Over the years, the limitations of JuMP’s nonlinear interface have been a
constant source of frustration for users. Why could you only create nonlinear
expressions in the @NL
macros? Why were the @NLconstraint
and @NLobjective
macros needed instead of the normal @constraint
and @objective
? Why could
you not use sum(::Vector)
or compute a dot product inside a nonlinear macro?
The discourse forum is filled
with questions, hacks, and work-arounds for all things related to nonlinear
programming.
Two years ago, and motivated by this public demand, we received funding from Los Alamos National Laboratory to begin a rewrite of JuMP’s nonlinear interface with the goal of bringing first-class nonlinear support to JuMP. This release is the culmination of that work.
Highlights
Here are some highlights of the new features.
Operator overloading
Nonlinear expressions can be defined outside of macros like their quadratic and affine counterparts:
julia> using JuMP
julia> model = Model();
julia> @variable(model, x);
julia> expr = x^2.3 + sin(x)
(x ^ 2.3) + sin(x)
julia> typeof(expr)
NonlinearExpr (alias for GenericNonlinearExpr{GenericVariableRef{Float64}})
julia> expr.head
:+
julia> expr.args
2-element Vector{Any}:
x ^ 2.3
sin(x)
Macro definitions
We can define nonlinear expressions inside of @expression
, @objective
, and
@constraint
. The special @NL
macros are no longer required.
julia> @expression(model, x + x^3 * sin(x)^x)
x + ((x ^ 3) * (sin(x) ^ x))
LinearAlgebra support
We can also do linear algebra operations:
julia> @variable(model, X[1:3, 1:2]);
julia> @variable(model, v[1:3]);
julia> @variable(model, w[1:2]);
julia> @expression(model, v' * X * w)
0.0 + ((X[1,1]*v[1] + X[2,1]*v[2] + X[3,1]*v[3]) * w[1]) + ((X[1,2]*v[1] + X[2,2]*v[2] + X[3,2]*v[3]) * w[2])
including broadcasting:
julia> @variable(model, Z[1:2, 1:2], Symmetric);
julia> @objective(model, Min, sum(Z^4 .- Z.^3))
0.0 + (((+(0.0) + (((+(0.0) + ((Z[1,2]) * (Z[1,1]))) + ...
# ... long expression not shown ...
Function tracing
Because we support operator overloading, calling a regular Julia function with JuMP variables produces a nonnlinear expression:
julia> my_func(y) = 2^y + exp(y^-2.3);
julia> my_func(2.0)
5.225149771643138
julia> my_func(x)
(2.0 ^ x) + exp(x ^ -2.3)
Function tracing should greatly simplify a number of nonlinear models. For example, the model in this Discourse thread becomes:
using JuMP, Ipopt
function example()
Q = -0.8:0.4:0.8
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, -2 <= p[1:5] <= 2)
@variable(model, -1 <= w <= 3)
@variable(model, -1 <= q <= 3)
@objective(model, Min, w)
f(p, q) = (1 / sqrt(2π)) * exp(-((p - q)^2) / 2)
total(p, q) = sum(_p * f(i, q) for (i, _p) in enumerate(p))
l1(p, q) = 1 - total(p, q) + 0.5 * total(p, 0.5)
l2(p, q) = total(p, q) - 1
lhs(p, q, _q) = l1(p, q) - l1(p, _q)
@constraint(model, [_q in Q], w * lhs(p, q, _q) + (1 - w) * l2(p, q) <= 0)
optimize!(model)
return
end
User-defined operators
To prevent tracing, you can create an operator instead:
julia> my_func2(a) = sin(a) * eta(a);
julia> @operator(model, op_my_func, 1, my_func2)
NonlinearOperator(my_func2, :op_my_func)
julia> @expression(model, op_my_func(x)^2)
op_my_func(x) ^ 2.0
Vector-valued nonlinear functions
We now support multi-objective programs with nonlinear objectives:
julia> model = Model();
julia> @variable(model, x[1:4]);
julia> @objective(model, Min, [sum(abs.(x)), sqrt(sum(x.^2))])
2-element Vector{NonlinearExpr}:
(((+(0.0) + abs(x[1])) + abs(x[2])) + abs(x[3])) + abs(x[4])
sqrt(x[1]² + x[2]² + x[3]² + x[4]²)
and we support nonlinear mixed complementarity problems
julia> model = Model();
julia> @variable(model, x >= 0);
julia> @constraint(model, log(x + 0.1) ⟂ x)
[log(x + 0.1), +(x)] ∈ MathOptInterface.Complements(2)
Potential breaking changes
Although the new nonlinear interface is a feature addition, there are two changes which might be breaking for a very small number of users.
- The syntax inside JuMP macros is parsed using a different code path, even for
linear and quadratic expressions. We made this change to unify how we parse
linear, quadratic, and nonlinear expressions. In all cases, the new code
returns equivalent expressions, but because of the different order of
operations, there are three changes to be aware of when updating:
- The printed form of the expression may change, for example from
x * y
toy * x
. This can cause tests which test theString
representation of a model to fail. - Some coefficients may change slightly due to floating point round-off error.
- Particularly when working with a JuMP extension, you may encounter a
MethodError
due to a missing or ambiguous method. These errors are due to previously existing bugs that were not triggered by the previous parsing code. If you encounter such an error, please open a GitHub issue.
- The printed form of the expression may change, for example from
- The methods for
Base.:^(x::VariableRef, n::Integer)
andBase.:^(x::AffExpr, n::Integer)
have changed. Previously, these methods supported onlyn = 0, 1, 2
and they always returned aQuadExpr
, even for the case whenn = 0
orn = 1
. Now:-
x^0
returnsone(T)
, whereT
is thevalue_type
of the model (defaults toFloat64
) -
x^1
returnsx
-
x^2
returns aQuadExpr
-
x^n
where!(0 <= n <= 2)
returns aNonlinearExpr
. We made this change to support nonlinear expressions and to align the mathematical definition of the operation with their return type. (Previously, users were surprised thatx^1
returned aQuadExpr
.) As a consequence of this change, the methods are now not type-stable. This means that the compiler cannot prove thatx^2
returns aQuadExpr
. If benchmarking shows that this is a performance problem, you can use the type-stablex * x
instead ofx^2
.
-
Next steps
Try it out! If you encounter any problems, or have any feedback on the new interface, please start a discourse thread or open a GitHub issue.