# New to Julia -- EAGO Nonlinear model works when finding Min but not Max?

**URL:** https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757
**Category:** Optimization (Mathematical)
**Tags:** question, package
**Created:** [September 24, 2022, 6:53pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757 "2022-09-24T18:53:10Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![lying\_syrup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lying_syrup/32/42933_2.png) [@lying\_syrup](https://discourse.julialang.org/u/lying_syrup)
#### Post date: [September 24, 2022, 6:53pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/1 "2022-09-24T18:53:10Z")

</div>

I’m not sure if a question like this fits more here or on the EAGO GitHub. I’m very new to Julia (started using it three days ago), and I’m also not sure if this is EAGO specific.  
In short, I’m trying to learn to use the JuMP optimizer syntax. When I put in a simple test optimization, the minimum of sin(x) or cos(x) for x in [-2, 2], the optimizer works for both (and their negatives), but not for the maximum. For the maximum, I get a KeyError with a _lot_ of in between stuff. The program is simple:

```julia
using JuMP
using EAGO

model = Model(EAGO.Optimizer)
@variable(model, x)
@NLconstraint(model, -2 <= x <= 2)
@NLobjective(model, Max, cos(x))
println("model = ", model)
optimize!(model)

```

I’m using Julia version 1.8.0, JuMP version 0.21.10, and EAGO version 0.6.1. Has anyone else encountered something like this for EAGO or other nonlinear solvers? What’s going on? Should I just negative anything I want the maxima of?  
The error text is below. Note that for different functions, a different key number is quoted in the error.

```julia
ERROR: LoadError: KeyError: key 2 not found
Stacktrace:
  [1] getindex
    @ ./dict.jl:498 [inlined]
  [2] forward_univariate_tiepnt_2!(k::Int64, op::Int64, child_idx::Int64, setstorage::Vector{MC{1, NS}}, x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, subgrad_tol::Float64, sparsity::Vector{Int64}, tpdict::Dict{Int64, NTuple{4, Int64}}, tp1storage::Vector{Float64}, tp2storage::Vector{Float64}, tp3storage::Vector{Float64}, tp4storage::Vector{Float64}, is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool, ctx::Cassette.Context{nametype(GuardCtx), EAGO.GuardTracker, Nothing, Cassette.var"##PassType#322", Nothing, Nothing})
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/functions/nonlinear/forward_pass.jl:758
  [3] forward_pass_kernel!(nd::Vector{JuMP._Derivatives.NodeData}, adj::SparseArrays.SparseMatrixCSC{Bool, Int64}, x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int64}, setstorage::Vector{MC{1, NS}}, numberstorage::Vector{Float64}, numvalued::Vector{Bool}, tpdict::Dict{Int64, NTuple{4, Int64}}, tp1storage::Vector{Float64}, tp2storage::Vector{Float64}, tp3storage::Vector{Float64}, tp4storage::Vector{Float64}, user_operators::JuMP._Derivatives.UserOperatorRegistry, subexpressions::Vector{EAGO.NonlinearExpression}, func_sparsity::Vector{Int64}, reverse_sparsity::Vector{Int64}, num_mv_buffer::Vector{Float64}, ctx::Cassette.Context{nametype(GuardCtx), EAGO.GuardTracker, Nothing, Cassette.var"##PassType#322", Nothing, Nothing}, is_post::Bool, is_intersect::Bool, is_first_eval::Bool, interval_intersect::Bool, cv_buffer::Vector{Float64}, cc_buffer::Vector{Float64}, treat_x_as_number::Vector{Bool}, subgrad_tol::Float64)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/functions/nonlinear/forward_pass.jl:1062
  [4] forward_pass!(evaluator::EAGO.Evaluator, d::EAGO.NonlinearExpression{MC{1, NS}})
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/functions/nonlinear/forward_pass.jl:1098
  [5] forward_pass!(evaluator::EAGO.Evaluator, d::EAGO.BufferedNonlinearFunction{MC{1, NS}})
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/functions/nonlinear/forward_pass.jl:1115
  [6] set_constraint_propagation_fbbt!(m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/domain_reduction.jl:677
  [7] preprocess!(t::EAGO.DefaultExt, m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/optimize/optimize_nonconvex.jl:595
  [8] preprocess!(m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/optimize/optimize_nonconvex.jl:1050
  [9] global_solve!(m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/optimize/optimize_nonconvex.jl:1085
 [10] optimize!(#unused#::Val{EAGO.MINCVX}, m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/optimize/optimize_nonconvex.jl:1162
 [11] optimize!(m::Optimizer)
    @ EAGO ~/.julia/packages/EAGO/5sz05/src/eago_optimizer/optimize/optimize.jl:37
 [12] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/YDdD3/src/Utilities/cachingoptimizer.jl:252
 [13] optimize!(b::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}}})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/YDdD3/src/Bridges/bridge_optimizer.jl:319
 [14] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.GenericModel{Float64, MathOptInterface.Utilities.ModelFunctionConstraints{Float64}}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/YDdD3/src/Utilities/cachingoptimizer.jl:252
 [15] optimize!(model::Model, optimizer_factory::Nothing; bridge_constraints::Bool, ignore_optimize_hook::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JuMP ~/.julia/packages/JuMP/klrjG/src/optimizer_interface.jl:185
 [16] optimize! (repeats 2 times)
    @ ~/.julia/packages/JuMP/klrjG/src/optimizer_interface.jl:146 [inlined]
 [17] top-level scope
    @ ~/Documents/Julia/test.jl:12

```

---

<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: [September 25, 2022, 12:53am UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/2 "2022-09-25T00:53:22Z")

</div>

This looks like a bug in EAGO.

You might like to try Juniper and the latest version of JuMP: [GitHub - lanl-ansi/Juniper.jl: A JuMP-based Nonlinear Integer Program Solver](https://github.com/lanl-ansi/Juniper.jl)

But if you dont have integer variables, you should first try Ipopt: [GitHub - jump-dev/Ipopt.jl: Julia interface to the Ipopt nonlinear solver](https://github.com/jump-dev/Ipopt.jl)

The jump documentation has a bunch of nonlinear tutorials to get you started: [Introduction · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/introduction/)

---

<div class="post-metadata">

### Author: ![lying\_syrup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lying_syrup/32/42933_2.png) [@lying\_syrup](https://discourse.julialang.org/u/lying_syrup)
#### Post date: [September 26, 2022, 1:41am UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/3 "2022-09-26T01:41:24Z")

</div>

Thanks!

---

<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: [September 26, 2022, 3:54am UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/4 "2022-09-26T03:54:51Z")

</div>

No problem 😄

For future readers, there are a bunch of open issues to update EAGO to the latest version of JuMP: [Add support for MOI 0.10 and JuMP 0.22 · Issue #100 · PSORLab/EAGO.jl · GitHub](https://github.com/PSORLab/EAGO.jl/issues/100), [Upcoming refactoring of JuMP's nonlinear API · Issue #104 · PSORLab/EAGO.jl · GitHub](https://github.com/PSORLab/EAGO.jl/issues/104)

---

<div class="post-metadata">

### Author: ![lying\_syrup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lying_syrup/32/42933_2.png) [@lying\_syrup](https://discourse.julialang.org/u/lying_syrup)
#### Post date: [September 28, 2022, 5:50pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/5 "2022-09-28T17:50:15Z")

</div>

Hello again, I was wondering if you could help me again here: I have been able to get Ipopt doing somewhat what I need it to do, but it’s a local solver. [Graphing Calculator](https://www.desmos.com/calculator/kma0aau4xj)  
To make a long story short, I’m trying to minimize the maximum value of a function’s range over an interval ([0, 1] to be particular). In doing so, I use Ipopt to first find the maximum value given a specific set of parameters like in the desmos link above, and then again use Ipopt to minimize _that_. Unfortunately, Ipopt seems to focus on the smaller corners more often than it finds the actual global result. Are there better methods of finding the global max other than seeding at many different values? My goal is to come away with a generalized function of any amount of input values in [0, 1], so just seeding will get far too slow for more inputs.  
I am fine with Ipopt finding local minima on the second run, but I need to have the global extrema for the first step to work.

The code is below.

```julia

using JuMP
using Ipopt

function find_diameter(n)
    if n%2 == 0
        return 2.0
    end
    x = (n+n%2)/(2*n)
    x_b = ceil(x*n) # x between: x is this and the next vertex
    x_a = (x*n)%1 # x along: how far is x between them (%)
    return sqrt((1-cos((2*π*x_b-2*π)/n)*(1-x_a)-cos((2*π*x_b)/n)*x_a)^2 + (sin((2*π*x_b-2*π)/n)*(1-x_a)+sin((2*π*x_b)/n)*x_a)^2)
end

function ramp_approx(x, stp, enp, slope, eps=10^-14)
    sta = x-stp
    fin = x-enp
    return slope*((sta/2)*(1+(sta/sqrt(sta^2+eps^2)))-((fin/2)*(1+(fin/sqrt(fin^2+eps^2)))))
end

function param(t, n)
    toReturn_x = 0
    toReturn_y = 0
    for k = 1:1:n
        toReturn_x += ramp_approx(t, (k-1)/n, k/n, n*(cos(π*(n-2*k)/n)-cos(π*(n-2*k+2)/n)))
        toReturn_y += ramp_approx(t, (k-1)/n, k/n, n*(sin(π*(n-2*k)/n)-sin(π*(n-2*k+2)/n)))
    end
    return (toReturn_x, toReturn_y)
end

function find_max(w, x, y)
    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "tol", 10^-14)
    #set_optimizer_attribute(model, "print_level", 8)
    @variable(model, 0 <= z <= 1)
    register(model, :dist, 3, dist; autodiff = true)
    @NLobjective(model, Max, (dist(z, w, 3)+dist(z, x, 3)+dist(z, y, 3))/3)
    optimize!(model)
    return objective_value(model)
end

function ∇find_max(g::AbstractVector{T}, w, x, y) where {T<:Real}
    h = 2^-20
    g[1] = (find_max(w+h,x,y)-find_max(w-h,x,y))/(2*h)
    g[2] = (find_max(w,x+h,y)-find_max(w,x-h,y))/(2*h)
    g[3] = (find_max(w,x,y+h)-find_max(w,x,y-h))/(2*h)
    return g
end

function magic_number(n)
    toReturn = 0.0
    if n%2 == 0
        for k = 0:n-1
            toReturn += sqrt(1.5+0.5*cos(2*π/n)-cos(2*k*π/n)-cos(2*(k-1)*π/n))/(2*n)
        end
    elseif n%2 == 1
        for k = 0:n-1
            toReturn += sqrt((1.5+0.5*cos(2*π/n)-cos(2*k*π/n)-cos(2*(k-1)*π/n))/(2-2*cos((n-1)*π/n)))/n
        end
    else
        return 0.0
    end
    return toReturn
end

model = Model(Ipopt.Optimizer)
register(model, :find_max, 3, find_max, ∇find_max)
set_optimizer_attribute(model, "tol", 10^-8)
set_optimizer_attribute(model, "max_wall_time", 60.0)
@variable(model, 0 <= w <= 1)
@variable(model, 0 <= x)
@variable(model, 0 <= y)
@NLconstraint(model, x <= w)
@NLconstraint(model, y <= x)
@NLconstraint(model, y <= w)
@NLobjective(model, Min, find_max(w, x, y))
optimize!(model)
println(" ")
println("model = ", model)
@show termination_status(model)
final_w = JuMP.value.(w)
final_x = JuMP.value.(x)
final_y = JuMP.value.(y)
final_output = find_max(final_w, final_x, final_y)

println("final w = ", final_w)
println("final x = ", final_x)
println("final y = ", final_y)
println("final output = ", final_output, ", divided by diameter: ", final_output/find_diameter(3))
println("target output = ", magic_number(3)*find_diameter(3))

```

---

<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: [September 28, 2022, 6:21pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/6 "2022-09-28T18:21:27Z")

</div>

I can’t run your code because I don’t know what `dist` is. Can you provide a reproducible example?

Your gradient computation is likely to run into issues. Because of tolerances, you could end up getting slightly different gradients each time you call it. You really need some way of computing the analytical gradient from the solution of `find_max`, but how exactly would depend one what `dist` is.

---

<div class="post-metadata">

### Author: ![lying\_syrup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lying_syrup/32/42933_2.png) [@lying\_syrup](https://discourse.julialang.org/u/lying_syrup)
#### Post date: [September 28, 2022, 6:26pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/7 "2022-09-28T18:26:30Z")

</div>

Sorry about that! I thought I’d included everything (I’m storing the functions in a separate file since they don’t need to change or take up space as I run things).

```julia
function dist(a,b,n)
    return ((param(a,n)[1]-param(b,n)[1])^2+(param(a,n)[2]-param(b,n)[2])^2)^0.5
end

```

I had to supply the gradient in order to have Ipopt optimize the results of another optimisation as seen in find\_max, automatic differentiation did not work with that at all. For my purposes, minimizable error is okay, this is for a mathematics project and these computations are only to help me find potential values which I will later prove on paper.  
In fact, the comparison with the output of magic\_number at the bottom is only there because I’m trying to _re-confirm_ already proven results at this stage.

Thanks for your help!

---

<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: [October 1, 2022, 6:58am UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/8 "2022-10-01T06:58:43Z")

</div>

Oops. I missed the reply. Sorry for the delay.

There’s a couple of style things I would change, like your `dist` function can be

```julia
function dist(a,b,n)
    ap, bp = param(a,n), param(b,n)
    return sqrt((ap[1]-pb[1])^2+(ap[2]-bp[2])^2)
end

```

But the bigger problem is that 1) You’re asking for a very tight tolerance inn the inner `find_max`. Ipopt likely can’t find an optimal solution to that tolerance, because when I run it it hits the maximum number of iterations. That’s related to the fact your gradients are finite differences.

Where does the ramp approximation come from? It looks like you’re approximating some function with break points? Is that the true function you’re optimizing over? Stepping back a bit, what is the underlying problem are you trying to solve?

---

<div class="post-metadata">

### Author: ![lying\_syrup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lying_syrup/32/42933_2.png) [@lying\_syrup](https://discourse.julialang.org/u/lying_syrup)
#### Post date: [October 2, 2022, 7:31pm UTC](https://discourse.julialang.org/t/new-to-julia-eago-nonlinear-model-works-when-finding-min-but-not-max/87757/9 "2022-10-02T19:31:03Z")

</div>

I’m trying to find rendezvous numbers/points for 2-d metric spaces. I’m not sure how accessible the small amount of literature out there is, so to sum it up:  
For any compact, connected metric space (X, d), there uniquely exists a number a(X, d) such that for any finite collection of vertices P, there exists another vertex y such that the average of the distance from y to all vertices in P is equal to a(X, d).  
As a result, if we examine the average distance as a function F from X → Positive Reals, we can find at least one collection P that maximizes the minimum of the range of F, and another collection Q that minimizes the maximum of the range of F. This is what I’m trying to do, by optimizing this function F, thereby finding the collections P and Q.  
As for the dist function, I wrote it that way in order for JuMP to accept it. That’s also the reason for the approximation of the step function.  
As of right now, I’m having success with Optim.jl, now I need to tweak parameters and choose the best solvers. I will likely rewrite the parameterization to be more exact again instead of a ramp function approximation since that seems unnecessary with Optim!

Thanks for your help!
