# Optimizing using objective with an integral in it

**URL:** <https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [February 27, 2022, 7:09pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162 "2022-02-27T19:09:09Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [February 27, 2022, 7:09pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/1 "2022-02-27T19:09:09Z")

</div>

I’m new to JuMP and I’d like to minimize a function that has an integral in it. Here’s an example f(y) = \int\_0^y y dy - y. Of course, this is just f(y)=(1/2)y^2-y which is minimized at y=1, but this is just a minimal example as a step towards something more complicated.

I don’t know how to pass this objective f into the JuMP stuff. I believe I need to register f as a “user-defined nonlinear function”, but I don’t know how to do that. According to the [docs](https://jump.dev/JuMP.jl/stable/manual/nlp/#User-defined-Functions), “in addition to this list of functions, it is possible to register custom _user-defined_ nonlinear functions.” But I don’t see how to do that and it doesn’t appear to say anything else.

Here’s the code I tried so far. The error message was: Unexpected object #89 (of type var"#89#90" in nonlinear expression.

Thanks!

```julia
using QuadGK
f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y

```

```julia
using JuMP

function example()
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @NLobjective(model, Min, f)
    optimize!(model)

    return
end

example()
# -> Unexpected object #89 (of type var"#89#90" in nonlinear expression.

```

---

<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 27, 2022, 7:24pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/2 "2022-02-27T19:24:13Z")

</div>

You should take a look at the JuMP documentation for some examples:

- [Introduction · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/introduction/)
- [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nlp/)

You need to define a variable `@variable(model, y)`, you need to call the function `f(y)`, and you’ll need to register the user-defined function `register(model, :f, 1, f; autodiff = true)`.

I haven’t run the code so there might be a typo etc, but this should point you in the right direction:

```julia
function example()
    model = Model(Ipopt.Optimizer)
    register(model, :f, 1, f; autodiff = true)
    set_silent(model)
    @variable(model, y)
    @NLobjective(model, Min, f(y))
    optimize!(model)
    return
end

```

Depending on what you are trying to do, you might also want to take a look at [https://github.com/pulsipher/InfiniteOpt.jl](https://github.com/pulsipher/InfiniteOpt.jl).

---

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [February 27, 2022, 9:06pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/3 "2022-02-27T21:06:16Z")

</div>

Thanks. From the docs [here](https://jump.dev/JuMP.jl/stable/reference/nlp/#User-defined-functions), I adjusted your code to replace the `y`’s in your second line with `foo` and `f`. (using `f` instead of `foo` works too). This worked when I use the analytically defined function such as `f(y) = 0.5*y^2 - y` and gives a solution of `y=1`. But when I used the equivalent function (up to numerical approx) `f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y` I get the following error saying something “repeats 79984 times”:

```julia
StackOverflowError:

Stacktrace:
 [1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, n::Int64) (repeats 79984 times)
   @ QuadGK C:\Users\alex\.julia\packages\QuadGK\ENhXl\src\gausskronrod.jl:253

```

Here’s my code. Any thoughts on why it’s not working when I use the numerical integral version of the objective?

```julia
f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y 
g(y) = 0.5*y^2 - y # g = f up to integration approx error
y = 0:0.1:10

function example()
    model = Model(Ipopt.Optimizer)
    @variable(model, y)
    register(model, :foo, 1, f; autodiff = true)
    set_silent(model)
    @NLobjective(model, Min, foo(y))
    optimize!(model)
    @show value(y)
    return
end

example()

```

---

<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 27, 2022, 9:39pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/4 "2022-02-27T21:39:54Z")

</div>

> [@samerb](#):
>
> I adjusted your code to replace the `y` 's in your second line with `foo` and `f` .

Ha. That was an obvious typo. I’ve edited the post.

The issue is that your user-defined function needs to be differentiable:

```julia
julia> using ForwardDiff, QuadGK

julia> f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1]
f (generic function with 1 method)

julia> ForwardDiff.derivative(f, 0.0)
ERROR: StackOverflowError:
Stacktrace:
[1] cachedrule( #unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, n::Int64 ) (repeats 79984 times)
@ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/gausskronrod.jl:253

```

Take a read of the posts here:

> [@Array issue with optimzing a function with quadgk](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837):
>
> I’m trying to minimize a function that involves an integral, the integrand of which is one of the optimizing parameters. I can compute the function and get the gradient with Zygote, but I’m getting an odd conversion error when I try to run an optimization routine. Might anyone know what is going on? using Distributions using Random using ForwardDiff using Zygote using QuadGK using Optim theta\_lower = 1.0 theta\_upper = 5.0 uniform\_dist = Uniform(theta\_lower, theta\_upper) function linearCost(θ)…

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [February 27, 2022, 10:47pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/5 "2022-02-27T22:47:30Z")

</div>

QuadGK will not work like this with JuMP. You might try Quadrature.jl

---

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [March 5, 2022, 11:05pm UTC](https://discourse.julialang.org/t/optimizing-using-objective-with-an-integral-in-it/77162/7 "2022-03-05T23:05:41Z")

</div>

It appears that Quadrature.jl doesn’t yet support differentiation of the bounds of the integrals. [https://github.com/SciML/Quadrature.jl/issues/56](https://github.com/SciML/Quadrature.jl/issues/56)

According to [Chirs R] Hasn’t been added yet. Someone would just need to extend [Quadrature.jl/Quadrature.jl at master · SciML/Quadrature.jl · GitHub](https://github.com/SciML/Quadrature.jl/blob/master/src/Quadrature.jl#L583) which wouldn’t be too hard.)