# Optimization with integrals in both objective and the constraint

**URL:** <https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838>\
**Category:** Optimization (Mathematical)\
**Tags:** nlsolve, quadgk, julia\
**Created:** [May 4, 2024, 3:53pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838 "2024-05-04T15:53:07Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [May 4, 2024, 3:53pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/1 "2024-05-04T15:53:08Z")

</div>

I’m trying to optimize a function that has an integral in the objective. While I can perform unconstrained optimization, I’m unsure how to use the Julia optimizer. I am trying to do this to add a constraint later, which itself is an integral.

```
using JuMP, Ipopt, Optim, ForwardDiff, QuadGK, NLsolve, Roots, LinearAlgebra, Base.Threads, Distributions

```

##### Set parameters

```
begin 
	r = 0.2 
	α = 0.85
	s = 0.5
	β = 0.9
	A = 1.0
	w = 10.0
	γ = 0.0
	ll = 0.0 # lower limit for optimization 
	ul = 1500.0 # upper limit for optimization 
	Dθ = Uniform(0.0, 5.0)
end

```

##### Define second period functions

```
begin
	θ̃(b, α, w, A, r, s) = min(max(R̄(b,r)/f(b,α),0.0),5.0)
	f(b,α) = b^α
	R̄(b,r) = (1+r)*b
	endY₂(b,θ,α) = θ*f(b,α)
	end#C₂(R₂,b,θ,α,w) = Y₂(b,θ,α) - R₂ + w
end

```

## Objective

```
begin
	Objhᵦ(b, θ, α, w, A, r, s, β) = β * (A*log(w) - s*(R̄(b,r)-Y₂(b,θ,α)))
	Objhᵪ(b, θ, α, w, A, r, s, β) = β * (A*log(Y₂(b,θ,α) - R̄(b,r) + w))
end

```

##### Integration over Dθ

```
begin 
	function FPopthᵦ(D, b, α, w, A, r, s, β)
		θ̲ = minimum(D)
		θ̄ = θ̃(b, α, w, A, r, s)
		return quadgk(y -> Objhᵦ(b, y, α, w, A, r, s, β)*pdf(D, y), θ̲, θ̄)[1]
	end
	
	function FPopthᵪ(D, b, α, w, A, r, s, β)
		θ̲ = θ̃(b, α, w, A, r, s)
		θ̄ = maximum(D)
		return quadgk(y -> Objhᵪ(b, y, α, w, A, r, s, β)*pdf(D, y), θ̲, θ̄)[1]
	end

	### objective function
	FPθh(Dθ, b, α, w, A, r, s, β) = FPopthᵦ(Dθ, b, α, w, A, r, s, β) + FPopthᵪ(Dθ, b, α, w, A, r, s, β)
end

```

##### Outcome

```
optimize(b -> -FPθh(Dθ, b, α, w, A, r, s, β), ll, ul).minimizer

```

## How to use Julia optimizer for the problem stated above?

```
function Opt()
    model = Model(Ipopt.Optimizer)
    @variable(model, 1 <= b <= 2, start = 1)
    register(model, :FP, 7, FPθh; autodiff = true)
    @NLobjective(model, Max, FPθ(Dθ, b, α, w, A, r, s, β))
    optimize!(model)
    return (b = value(b), value = objective_value(model))
end

```

##### This doesn’t work:

```
    Opt()

```

### Finally, how to add the constraint?

```
begin
	function πbh(D, b, α, w, A, r, s)
		θ̲ = minimum(D)
		θ̄ = θ̃(b, α, w, A, r, s)
		return quadgk(y -> Y₂(b,y,α)*pdf(D, y), θ̲, θ̄)[1]
	end
	
	function πch(D, b, α, w, A, r, s)
		θ̲ = θ̃(b, α, w, A, r, s)
		θ̄ = maximum(D)
		return quadgk(y -> R̄(b,r)*pdf(D, y), θ̲, θ̄)[1]
	end

	πθh(Dθ, b, α, w, A, r, s, γ) = πbh(Dθ, b, α, w, A, r, s)+πch(Dθ, b, α, w, A, r, s) - (1+γ)*b
end

```

##### Constraint: πθh(Dθ, b, α, w, A, r, s, γ) ≥ -1

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 5, 2024, 5:19am UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/2 "2024-05-05T05:19:15Z")

</div>

> [@Ayesha](#):
>
> This doesn’t work:

Can you share what the error says?

> [@Ayesha](#):
>
> ` @NLobjective(model, Max, FPθ(Dθ, b, α, w, A, r, s, β))`

You are using the old version of the nonlinear API, starting with JuMP.jl you don’t need `@NL` in the macros anymore. Check out

> **[Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nonlinear/)**
>
> Documentation for JuMP.

---

<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:** [May 5, 2024, 8:28pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/3 "2024-05-05T20:28:06Z")

</div>

If I fix up your syntax slightly, and make a single reproducible example, I get:

```Julia
julia> using JuMP, Ipopt, QuadGK, Distributions

julia> r = 0.2
0.2

julia> α = 0.85
0.85

julia> s = 0.5
0.5

julia> β = 0.9
0.9

julia> A = 1.0
1.0

julia> w = 10.0
10.0

julia> γ = 0.0
0.0

julia> ll = 0.0
0.0

julia> ul = 1500.0
1500.0

julia> Dθ = Uniform(0.0, 5.0)
Uniform{Float64}(a=0.0, b=5.0)

julia> θ̃(b, α, w, A, r, s) = min(max(R̄(b,r)/f(b,α),0.0),5.0)
θ̃ (generic function with 1 method)

julia> f(b, α) = b^α
f (generic function with 1 method)

julia> R̄(b, r) = (1+r)*b
R̄ (generic function with 1 method)

julia> Y₂(b,θ,α) = θ*f(b,α)
Y₂ (generic function with 1 method)

julia> Objhᵦ(b, θ, α, w, A, r, s, β) = β * (A*log(w) - s*(R̄(b,r)-Y₂(b,θ,α)))
Objhᵦ (generic function with 1 method)

julia> Objhᵪ(b, θ, α, w, A, r, s, β) = β * (A*log(Y₂(b,θ,α) - R̄(b,r) + w))
Objhᵪ (generic function with 1 method)

julia> function FPopthᵦ(D, b, α, w, A, r, s, β)
           θ̲ = minimum(D)
           θ̄ = θ̃(b, α, w, A, r, s)
           return quadgk(y -> Objhᵦ(b, y, α, w, A, r, s, β)*pdf(D, y), θ̲, θ̄)[1]
       end
FPopthᵦ (generic function with 1 method)

julia> function FPopthᵪ(D, b, α, w, A, r, s, β)
           θ̲ = θ̃(b, α, w, A, r, s)
           θ̄ = maximum(D)
           return quadgk(y -> Objhᵪ(b, y, α, w, A, r, s, β)*pdf(D, y), θ̲, θ̄)[1]
       end
FPopthᵪ (generic function with 1 method)

julia> FPθh(Dθ, b, α, w, A, r, s, β) = FPopthᵦ(Dθ, b, α, w, A, r, s, β) + FPopthᵪ(Dθ, b, α, w, A, r, s, β)
FPθh (generic function with 1 method)

julia> function Opt()
           model = Model(Ipopt.Optimizer)
           @variable(model, 1 <= b <= 2, start = 1)
           FP(b) = FPθh(Dθ, b, α, w, A, r, s, β)
           @operator(model, op_FP, 1, FP)
           @objective(model, Max, op_FP(b))
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return (b = value(b), value = objective_value(model))
       end
Opt (generic function with 1 method)

julia> Opt()
ERROR: Unable to register the function :op_FP.
... many more lines ...

```

The issue is that JuMP does not support automatically differentiating your `FPθh` function because of the `quadgk`.

> **[Should you use JuMP? · JuMP](https://jump.dev/JuMP.jl/stable/should_i_use/#You-want-to-optimize-a-complicated-Julia-function)**
>
> Documentation for JuMP.

I’m not sure why Optim works? It might be doing something different:

```Julia
using Optim
optimize(b -> -FPθh(Dθ, b, α, w, A, r, s, β), ll, ul).minimizer

```

---

<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:** [May 5, 2024, 8:32pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/4 "2024-05-05T20:32:22Z")

</div>

I guess Optim is doing finite differencing, in which case, you could do:

```julia
function Opt()
    model = Model(Ipopt.Optimizer)
    @variable(model, 1 <= b <= 2, start = 1)
    FP(b) = FPθh(Dθ, b, α, w, A, r, s, β)
    finite_diff(f, h = 1e-8) = x -> (f(x + h) - f(x)) / h
    ∇FP = finite_diff(FP)
    ∇²FP = finite_diff(∇FP)
    @operator(model, op_FP, 1, FP, ∇FP, ∇²FP)
    @objective(model, Max, op_FP(b))
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return (b = value(b), value = objective_value(model))
end

Opt()

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 6, 2024, 6:14am UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/5 "2024-05-06T06:14:56Z")

</div>

You can also use ForwardDiff.jl, which QuadGK.jl natively supports.

---

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [May 9, 2024, 6:26am UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/6 "2024-05-09T06:26:49Z")

</div>

The code using finite difference approximation is working! Thank you so much!

- I have encountered another issue: when I set a different boundary for the choice variable, for example, `1 <= b <= 4`, it shows an AssertionError. The same error occurs when the upper limit is set to 5, 7, 8, or 9, but it works when set to 6, 10, etc. Do you have any idea why this could be happening?

- Finally, is it possible to add a constraint that involves an integral? I have tried using a constraint that’s not an integral, and it works.  
`@constraint(model, πθh(Dθ, b, α, w, A, r, s, γ) >= -1)`

---

<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:** [May 9, 2024, 7:51am UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/7 "2024-05-09T07:51:09Z")

</div>

> it shows an AssertionError

What is the error?

> is it possible to add a constraint that involves an integral

JuMP does not have anything built-in for integrals. You’ll need to implement the necessary constraint.

---

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [May 9, 2024, 4:45pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/8 "2024-05-09T16:45:09Z")

</div>

> [@odow](#):
>
> What is the error?

It throws this:

```julia
AssertionError: is_solved_and_feasible(model)

Opt()@Other: 12
top-level scope@Local: 1[inlined]

```

The log says that the problem is solved to an acceptable level.

```julia
Number of Iterations....: 43

                                   (scaled) (unscaled)
Objective...............: -2.1237353163035482e+00 2.1237353163035482e+00
Dual infeasibility......: 4.4408920985006262e-07 4.4408920985006262e-07
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.6505208684905717e-08 2.6505208684905717e-08
Overall NLP error.......: 4.4408920985006262e-07 4.4408920985006262e-07

Number of objective function evaluations = 174
Number of objective gradient evaluations = 44
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 43
Total seconds in IPOPT = 0.136

EXIT: Solved To Acceptable Level.

```

But it doesn’t show the maximizer. When I plot it as `b(s,r)`, it looks like this:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/b/5b734aed66ec980719db198bdd451c820ca7a5f7.png)

> [@odow](#):
>
> You’ll need to implement the necessary constraint.

Could you please give a bit of direction on how to do this?

---

<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:** [May 9, 2024, 8:04pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/9 "2024-05-09T20:04:59Z")

</div>

See the documentation at: [Solutions · JuMP](https://jump.dev/JuMP.jl/stable/manual/solutions/).

`solution_summary(model)` will provide more information.

You may need:

```Julia
@assert is_solved_and_feasible(model; allow_almost = true)

```

This says that the problem is _not_ solved and feasible, but that it was nearly. There are likely some minor issues with tolerances, etc. As one example, the overall NLP error is 4e-7, but the default tolerance is 1e-8. So things are “nearly” solved.

> is it possible to add a constraint that involves an integral  
> Could you please give a bit of direction on how to do this?

This is a very broad question with no particular answer. It depends on your model and the constraint you want to add.

Presumably, you want to register another operator for your `πθh` function, and then add `@constraint(model, op_πθh(b) >= -1)`.

---

<div class="post-metadata">

**Author:** ![Ayesha](https://avatars.discourse-cdn.com/v4/letter/a/c0e974/32.png) [@Ayesha](https://discourse.julialang.org/u/Ayesha)\
**Post date:** [May 9, 2024, 11:54pm UTC](https://discourse.julialang.org/t/optimization-with-integrals-in-both-objective-and-the-constraint/113838/10 "2024-05-09T23:54:36Z")

</div>

Thanks for your help and for sharing the links.  
This seems to be working for now:

```julia
function Opt(α, w, A, r, s, β)
    model = Model(Ipopt.Optimizer)
    @variable(model, 1 <= b <= 10, start = 1)
    FP(b) = FPθh(Dθ, b, α, w, A, r, s, β)
    finite_diff(f, h = 1e-6) = x -> (f(x + h) - f(x)) / h
    ∇FP = finite_diff(FP)
    ∇²FP = finite_diff(∇FP)
    @operator(model, op_FP, 1, FP, ∇FP, ∇²FP)

    πh(b) = πθh(Dθ, b, α, w, A, r, s, γ)
    ∇πh = finite_diff(πh)
    ∇²πh = finite_diff(∇πh)
    @operator(model, op_πh, 1, πh, ∇πh, ∇²πh)
    @constraint(model, op_πh(b) >= 0.0)

    @objective(model, Max, op_FP(b))
    optimize!(model)
    @assert is_solved_and_feasible(model; allow_almost = true)
    return (b = value(b), value = objective_value(model), sum = solution_summary(model))
end

```

![image](https://global.discourse-cdn.com/julialang/original/3X/a/0/a03aa840ea310de22cbc97a721ebfd57d21c3d38.png)