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

[BuRanran]

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.

[odow]

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.

[BuRanran]

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?

[odow]

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)

[BuRanran]

It works. Thanks so much!

[anand.kumar]

Have something similar and was trying to replicate this. However, it seems syntax has updated for SymPy. Got it working, but convert doesn’t seem to work.

update: realised func(z) returns float64. Thus, I may not need convert (maybe?). However, adding this in objective function yields back erf in the expression and I’m back to square -

In my case, objective function involves f(u,w) + g(u)*erf(u/w). Since, erf isn’t defined in Julia natively, I’m hoping, resorting to your approach could work for my case too. I’m using standard solver e.g. Gurobi/Knitro. Any other insights would be helpful.

Thank You.

[odow]

Hi @anand.kumar,

erf is supported as part of SpecialFunctions.jl

See the JuMP documentation section: Nonlinear Modeling · JuMP

If you’re still stuck, please start a new thread and we can discusss there.