# JuMP(+Juniper) Error registering user-defined function

**URL:** <https://discourse.julialang.org/t/jump-juniper-error-registering-user-defined-function/94909>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [February 20, 2023, 10:32am UTC](https://discourse.julialang.org/t/jump-juniper-error-registering-user-defined-function/94909 "2023-02-20T10:32:38Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![PrimulaBunce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/primulabunce/32/46933_2.png) [@PrimulaBunce](https://discourse.julialang.org/u/PrimulaBunce)\
**Post date:** [February 20, 2023, 10:32am UTC](https://discourse.julialang.org/t/jump-juniper-error-registering-user-defined-function/94909/1 "2023-02-20T10:32:38Z")

</div>

Hi,

I’m trying to register a user-defined function in JuMP and get the following error:

```julia
ERROR: Unable to register the function :f_objective.

Common reasons for this include:

 * The function takes `f(x::Vector)` as input, instead of the splatted
   `f(x...)`.
 * The function assumes `Float64` will be passed as input, it must work for any
   generic `Real` type.
 * The function allocates temporary storage using `zeros(3)` or similar. This
   defaults to `Float64`, so use `zeros(T, 3)` instead.

As examples, instead of:
`my_function(x::Vector) = sum(x.^2)`
use:
`my_function(x::T...) where {T<:Real} = sum(x[i]^2 for i in 1:length(x))`

Instead of:
`function my_function(x::Float64...)
    y = zeros(length(x))
    for i in 1:length(x)
        y[i] = x[i]^2
    end
    return sum(y)
end`
use:
`function my_function(x::T...) where {T<:Real}
    y = zeros(T, length(x))
    for i in 1:length(x)
        y[i] = x[i]^2
    end
    return sum(y)
end`

Review the stacktrace below for more information, but it can often be hard to
understand why and where your function is failing. You can also debug this
outside of JuMP as follows:
`import ForwardDiff

# If the input dimension is 1
x = 1.0
my_function(a) = a^2
ForwardDiff.derivative(my_function, x)

# If the input dimension is more than 1
x = [1.0, 2.0]
my_function(a, b) = a^2 + b^2
ForwardDiff.gradient(x -> my_function(x...), x)`

Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] _validate_register_assumptions(f::typeof(f_objective), name::Symbol, dimension::Int64)
   @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:267
 [3] (MathOptInterface.Nonlinear._MultivariateOperator{20, F′, F′′} where {F′, F′′})(op::Symbol, f::Function)
   @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:297
 [4] register_operator(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, nargs::Int64, f::Function)
   @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:350
 [5] register_operator(model::MathOptInterface.Nonlinear.Model, op::Symbol, nargs::Int64, f::Function)      
   @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\model.jl:217
 [6] register(model::Model, op::Symbol, dimension::Int64, f::Function; autodiff::Bool)
   @ JuMP C:\Users\kjaya\.julia\packages\JuMP\yYfHy\src\nlp.jl:677
 [7] top-level scope
   @ g:\Other computers\PC\PhD\Courses\D1000 Linear Control Theory\PhD paper 1\Git\latestOrganizedCode.jl:87
caused by: MethodError: no method matching schur!(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}})
Closest candidates are:
  schur!(::StridedMatrix{var"#s814"} where var"#s814"<:Union{Float32, Float64, ComplexF32, ComplexF64}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\schur.jl:97
  schur!(::StridedMatrix{T}, ::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\schur.jl:288
Stacktrace:
  [1] schur(A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\schur.jl:143
  [2] lyapc(A::Adjoint{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}}}, C::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}})
    @ MatrixEquations C:\Users\kjaya\.julia\packages\MatrixEquations\pS06e\src\lyapunov.jl:68
  [3] f_objective(::ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}, ::Vararg{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}, N} where N)
    @ Main g:\Other computers\PC\PhD\Courses\D1000 Linear Control Theory\PhD paper 1\Git\latestOrganizedCode.jl:82
  [4] (::MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)})(x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}})
    @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:263
  [5] chunk_mode_gradient(f::MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}}})
    @ ForwardDiff C:\Users\kjaya\.julia\packages\ForwardDiff\QdStj\src\gradient.jl:150
  [6] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}}}, ::Val{true})     
    @ ForwardDiff C:\Users\kjaya\.julia\packages\ForwardDiff\QdStj\src\gradient.jl:21
  [7] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{MathOptInterface.Nonlinear.var"#15#16"{typeof(f_objective)}, Float64}, Float64, 10}}}) (repeats 2 times) @ ForwardDiff C:\Users\kjaya\.julia\packages\ForwardDiff\QdStj\src\gradient.jl:17
  [8] _validate_register_assumptions(f::typeof(f_objective), name::Symbol, dimension::Int64)
    @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:263
  [9] (MathOptInterface.Nonlinear._MultivariateOperator{20, F′, F′′} where {F′, F′′})(op::Symbol, f::Function)
    @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:297
 [10] register_operator(registry::MathOptInterface.Nonlinear.OperatorRegistry, op::Symbol, nargs::Int64, f::Function)
    @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\operators.jl:350
 [11] register_operator(model::MathOptInterface.Nonlinear.Model, op::Symbol, nargs::Int64, f::Function)     
    @ MathOptInterface.Nonlinear C:\Users\kjaya\.julia\packages\MathOptInterface\NCblk\src\Nonlinear\model.jl:217
 [12] register(model::Model, op::Symbol, dimension::Int64, f::Function; autodiff::Bool)
    @ JuMP C:\Users\kjaya\.julia\packages\JuMP\yYfHy\src\nlp.jl:677
 [13] top-level scope
    @ g:\Other computers\PC\PhD\Courses\D1000 Linear Control Theory\PhD paper 1\Git\latestOrganizedCode.jl:87

```

Here’s a minimal working example:

```julia
using JuMP, Ipopt, HiGHS, Juniper, LinearAlgebra, MatrixEquations

# Setting up the solver
nlp_solver = optimizer_with_attributes(Ipopt.Optimizer, MOI.Silent() => true, "sb" => "yes", "max_iter" => Int(1E4), );
mip_solver = optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "on", "log_to_console" => false, );
minlp_solver = optimizer_with_attributes(Juniper.Optimizer, MOI.Silent() => true, "mip_solver" => mip_solver, "nl_solver" => nlp_solver, );

m = Model(minlp_solver);

@variable(m, x[1:20], Bin);

function f_objective(x::T...) where {T<:Real}
    A = diagm( 0 => (-1-1)*ones(T,20), 
                1 => (2-1)*ones(T,20-1), 
                -1 => (3-1)*ones(T,20-1),
                2 => (-5-1)*ones(T,20-2),
                -2 => (-4-1)*ones(T,20-2)) + ones(T,20,20);

    x = [x...]

    C = diagm(x)

    Wo = MatrixEquations.lyapc(A', A'*C)

    return(tr(Wo))
end
register(m, :f_objective, 20, f_objective; autodiff=true)

```

The problem seems to be in `Wo = MatrixEquations.lyapc(A', A'*C)`…

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [February 20, 2023, 4:17pm UTC](https://discourse.julialang.org/t/jump-juniper-error-registering-user-defined-function/94909/2 "2023-02-20T16:17:45Z")

</div>

The problem seems to be in LinearAlgebra and how it uses the shur!() function. This runs without issues:

```julia
using JuMP, Ipopt, HiGHS, Juniper, MatrixEquations

# Setting up the solver
nlp_solver = optimizer_with_attributes(Ipopt.Optimizer, MOI.Silent() => true, "sb" => "yes", "max_iter" => Int(1E4), );
mip_solver = optimizer_with_attributes(HiGHS.Optimizer, "presolve" => "on", "log_to_console" => false, );
minlp_solver = optimizer_with_attributes(Juniper.Optimizer, MOI.Silent() => true, "mip_solver" => mip_solver, "nl_solver" => nlp_solver, );

m = Model(minlp_solver);
# Setting up the solver

m = Model();

@variable(m, x[1:20], Bin);

function f_objective(x::T...) where {T<:Real}
    A = diagm( 0 => (-1-1)*ones(T,20), 
                1 => (2-1)*ones(T,20-1), 
                -1 => (3-1)*ones(T,20-1),
                2 => (-5-1)*ones(T,20-2),
                -2 => (-4-1)*ones(T,20-2)) + ones(T,20,20);

    x = [x...]

    C = diagm(x)

    Wo = MatrixEquations.lyapc(A', A'*C)
    return(tr(Wo))
end
register(m, :f_objective, 20, f_objective; autodiff=true)

```

Looking at the stacktrace, it can be the second of the ‘common reasons’, i.e. at some point LinearAlgebra requires a float and cannot handle a real number. Not sure how to solve except reaching out to the creators of LinearAlgebra on github: [here I guess](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/LinearAlgebra.jl)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 20, 2023, 7:46pm UTC](https://discourse.julialang.org/t/jump-juniper-error-registering-user-defined-function/94909/3 "2023-02-20T19:46:36Z")

</div>

As @blob suggests, you need `f_objective` to be differentiable. Unfortunately it isn’t because of the use of `MatrixEquations.lyapc`.

Instead of using a user-defined function, you should write out the equation using linear constraints:

```julia
using JuMP
import HiGHS
import LinearAlgebra
model = Model(HiGHS.Optimizer)
@variable(model, x[1:20], Bin)
@variable(model, X[1:20, 1:20], Symmetric)
A = LinearAlgebra.diagm(
    0 => (-1 - 1) * ones(20), 
    1 => (2 - 1) * ones(20 - 1), 
    -1 => (3 - 1) * ones(20 - 1),
    2 => (-5 - 1) * ones(20 - 2),
    -2 => (-4 - 1) * ones(20 - 2)
) + ones(20, 20)
C = LinearAlgebra.diagm(x)
@constraint(model, A' * X + X * A + A' * C .== 0)
@objective(model, Min, LinearAlgebra.tr(X))

```

When I tested I only got the solution `x = 0`, so I might have made a mistake, but it should point you in the right direction.