Ipopt v1.8.0 adds the `Ipopt._VectorNonlinearOracle` set. The purpose of this is…sue is to explain our motivation for adding the set, and what we plan to do about it in the future.
## Current status and future plans
`Ipopt._VectorNonlinearOracle` is experimental. You are encouraged to try it out and report back experiences (see usage instructions below). It is marked as experimental (it begins with an underscore) because we plan to make a change to the API no later than September 30, 2025 (six months after adding it).
We are considering three options:
1. Rename the set to `Ipopt.VectorNonlinearFunction` and publicly support it as part of the v1.X API of Ipopt.jl
2. Move the set to MathOptInterface as part of the public API for MathOptInterface.jl, and add support for it in other packages like KNITRO.jl and NLopt.jl
3. Remove the set in favor of a different approach to vector-valued user defined functions.
If the set is useful for a small number of users who use only Ipopt.jl, and no one shows wider interest, we will choose (1).
If the set is useful for a number of different groups, including by people who want to use it with different NLP solvers, and the API proves to be a good design, we will choose (2).
If issues arise, we reserve the right to delete the set entirely and pursue a different approach in (3).
**Thus, if you use this set please comment below with feedback of how you are using it and whether we could improve it.**
## Compatibility
The `_VectorNonlinearOracle ` feature is experimental. It relies on a private API feature of Ipopt.jl that will change in a future release. If you use this feature, you must pin the version of Ipopt.jl in your `Project.toml` to ensure that future updates to Ipopt.jl do not break your existing code.
A known good version of Ipopt.jl is v1.8.0. Pin the version using:
```
[compat]
Ipopt = "=1.8.0"
```
## Definition
`Ipopt._VectorNonlinearOracle` is defined as:
```julia
_VectorNonlinearOracle(;
dimension::Int,
l::Vector{Float64},
u::Vector{Float64},
eval_f::Function,
jacobian_structure::Vector{Tuple{Int,Int}},
eval_jacobian::Function,
hessian_lagrangian_structure::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[],
eval_hessian_lagrangian::Union{Nothing,Function} = nothing,
) <: MOI.AbstractVectorSet
```
which represents the set:
```math
S = \{x \in \mathbb{R}^{dimension}: l \le f(x) \le u\}
```
where ``f`` is defined by the vectors `l` and `u`, and the callback oracles `eval_f`, `eval_jacobian`, and `eval_hessian_lagrangian`.
### Function
The `eval_f` function must have the signature
```julia
eval_f(ret::AbstractVector, x::AbstractVector)::Nothing
```
which fills $f(x)$ into the dense vector `ret`.
### Jacobian
The `eval_jacobian` function must have the signature
```julia
eval_jacobian(ret::AbstractVector, x::AbstractVector)::Nothing
```
which fills the sparse Jacobian $\nabla f(x)$ into `ret`.
The one-indexed sparsity structure must be provided in the `jacobian_structure` argument.
### Hessian
The `eval_hessian_lagrangian` function is optional.
If `eval_hessian_lagrangian === nothing`, Ipopt will use a Hessian approximation instead of the exact Hessian.
If `eval_hessian_lagrangian` is a function, it must have the signature
```julia
eval_hessian_lagrangian(
ret::AbstractVector,
x::AbstractVector,
μ::AbstractVector,
)::Nothing
```
which fills the sparse Hessian of the Lagrangian $\sum \mu_i \nabla^2 f_i(x)$ into `ret`.
The one-indexed sparsity structure must be provided in the `hessian_lagrangian_structure` argument.
## Use with JuMP
This section provides some usage examples for JuMP.
### Example: f(x) <= 1
A simplest example is to represent a constraint `x^2 + y^2 <= 1` in the problem:
```math
\begin{aligned}
Max & x + y \\
& x^2 + y^2 \le 1
\end{aligned}
```
This can be achieved as follows
```julia
julia> using JuMP, Ipopt
julia> set = Ipopt._VectorNonlinearOracle(;
dimension = 2,
l = [-Inf],
u = [1.0],
eval_f = (ret, x) -> (ret[1] = x[1]^2 + x[2]^2),
jacobian_structure = [(1, 1), (1, 2)],
eval_jacobian = (ret, x) -> ret .= 2.0 .* x,
hessian_lagrangian_structure = [(1, 1), (2, 2)],
eval_hessian_lagrangian = (ret, x, u) -> ret .= 2.0 .* u[1],
);
julia> begin
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@objective(model, Max, x + y)
@constraint(model, [x, y] in set)
optimize!(model)
assert_is_solved_and_feasible(model)
value(x), value(y), 1 / sqrt(2)
end
(0.707106783471979, 0.707106783471979, 0.7071067811865475)
```
### Example: y = F(x)
To model $y = F(x)$, we need to implement the set $[x, y] \\in S$ such that $0 \\le F(x) - y \\le 0$
We can modify our previous example to:
```math
\begin{aligned}
Max & x + y \\
& x^2 + y^2 - z = 0 \\
& z \le 1
\end{aligned}
```
which is implemented as follows:
```Julia
julia> using JuMP, Ipopt
julia> set = Ipopt._VectorNonlinearOracle(;
dimension = 3,
l = [0.0],
u = [0.0],
eval_f = (ret, x) -> (ret[1] = x[1]^2 + x[2]^2 - x[3]),
jacobian_structure = [(1, 1), (1, 2), (1, 3)],
eval_jacobian = (ret, x) -> begin
ret[1:2] .= 2.0 .* x[1:2]
ret[3] = -1.0
end,
hessian_lagrangian_structure = [(1, 1), (2, 2)],
eval_hessian_lagrangian = (ret, x, u) -> ret .= 2.0 .* u[1],
);
julia> begin
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@variable(model, z <= 1)
@objective(model, Max, x + y)
# x^2 + y^2 - z == 0
@constraint(model, [x, y, z] in set)
optimize!(model)
assert_is_solved_and_feasible(model)
value(x), value(y), value(z), 1 / sqrt(2)
end
(0.707106783471979, 0.707106783471979, 1.0000000064625545, 0.7071067811865475)
```
### Example: multiple rows
A more involved example is:
```math
\begin{aligned}
& x^2 - a & = 0 \\
& y^2 + z^3 - b & = 0
\end{aligned}
```
which is implemented as follows:
```Julia
julia> set = Ipopt._VectorNonlinearOracle(;
dimension = 5,
l = [0.0, 0.0],
u = [0.0, 0.0],
eval_f = (ret, x) -> begin
ret[1] = x[1]^2 - x[4]
ret[2] = x[2]^2 + x[3]^3 - x[5]
return
end,
jacobian_structure = [(1, 1), (2, 2), (2, 3), (1, 4), (2, 5)],
eval_jacobian = (ret, x) -> begin
ret[1] = 2 * x[1]
ret[2] = 2 * x[2]
ret[3] = 3 * x[3]^2
ret[4] = -1.0
ret[5] = -1.0
return
end,
hessian_lagrangian_structure = [(1, 1), (2, 2), (3, 3)],
eval_hessian_lagrangian = (ret, x, u) -> begin
ret[1] = 2 * u[1]
ret[2] = 2 * u[2]
ret[3] = 6 * x[3] * u[2]
return
end,
);
julia> begin
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, x[i in 1:3] == i)
@variable(model, y[1:2])
@constraint(model, [x; y] in set)
optimize!(model)
assert_is_solved_and_feasible(model)
value.(x), value.(y)
end
([1.0, 2.0, 3.0], [1.0, 31.0])
```
### Example: no Hessian
The Hessian evaluator is optional:
```julia
julia> using JuMP, Ipopt
julia> set = Ipopt._VectorNonlinearOracle(;
dimension = 2,
l = [-Inf],
u = [1.0],
eval_f = (ret, x) -> (ret[1] = x[1]^2 + x[2]^2),
jacobian_structure = [(1, 1), (1, 2)],
eval_jacobian = (ret, x) -> ret .= 2.0 .* x,
);
julia> begin
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, 0 <= x <= 1)
@variable(model, 0 <= y <= 1)
@objective(model, Max, x + y)
@constraint(model, [x, y] in set)
optimize!(model)
assert_is_solved_and_feasible(model)
value(x), value(y), 1 / sqrt(2)
end
(0.7071067847123265, 0.7071067847123265, 0.7071067811865475)
```