ERROR: The objective function (erf) is not supported by JuMP

I have a non-linear optimization problem to be solved with JuMP and IPOPT solver. The objective function is a probability distribution function of a normal distribution. It can be optimized by IPOPT directly, but not supported by JuMP. Is there any method to solve this problem? As the constrains are complicated, I hope to write constrains with JuMP. Thanks a lot !
The codes with error are as follows.

func1(x)=1/(2*3.1415926)^0.5/20*exp(-(x-100)^2/2/20^2)
gp = Model(with_optimizer(Ipopt.Optimizer))
@variable(gp, Q)
@objective(gp, Max, integrate(func1)(Q))
optimize!(gp)

ERROR: LoadError: The objective function 0.199471141902022*sqrt(2)*sqrt(pi)*erf(sqrt(2)*(x - 100)/40) is not supported by JuMP.

Here are the docs: https://jump.dev/JuMP.jl/stable/nlp/
You need to use @NLobjective and @NLconstraint.

You will also need to register the integrate function as a user-defined function, so something like:

using JuMP
using Ipopt

gp = Model(Ipopt.Optimizer)
@variable(gp, Q)
func1 = integrate(1/(2*3.1415926)^0.5/20*exp(-(x-100)^2/2/20^2))
register(gp, :func1, 1, func1; autodiff=true)
@NLobjective(gp, Max, func1(Q))
optimize!(gp)

Note that I haven’t tested this because I don’t know where integrate comes from.

Thanks for your reply. I correct the @objective to @NLobjective, it really helps. But there are still errors.
Ps. The integrate() is a integral function in package “SymPy”

using JuMP
using Ipopt
using SymPy

gp = Model(Ipopt.Optimizer)
@variable(gp, Q)
func1(x)=1/(2*3.1415926)^0.5/20*exp(-(x-100)^2/2/20^2)
func2(y)=integrate(func1, -Inf, y)
register(gp, :func2, 1, func2; autodiff=true)
@NLobjective(gp, Max, func2(Q))
optimize!(gp)

ERROR: LoadError: Expected return type of Float64 from a user-defined function, but got Sym.

That the type of func2(constant) is Sym, for example:

julia> func2(100)
0.199471141902022⋅√2⋅√π

julia> typeof(func2(100))
Sym

I want to transfer Sym to float by using function float(func2(Q)), but function float() can not be used in the objective. How can I fix this problem?

You want something like

using JuMP
using Ipopt
using SymPy

func1(x)=1/(2*3.1415926)^0.5/20*exp(-(x-100)^2/2/20^2)
func2 = lambdify(integrate(func1))
func3(y) = convert(Real, func2(y))

gp = Model(Ipopt.Optimizer)
@variable(gp, Q)
register(gp, :func3, 1, func3; autodiff=true)
@NLobjective(gp, Max, func3(Q))
optimize!(gp)

It works. Thanks so much!

1 Like