# Nonlinear program with JuMP not following all constraints

**URL:** <https://discourse.julialang.org/t/nonlinear-program-with-jump-not-following-all-constraints/123691>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [December 11, 2024, 7:28am UTC](https://discourse.julialang.org/t/nonlinear-program-with-jump-not-following-all-constraints/123691 "2024-12-11T07:28:33Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yuricst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuricst/32/207613_2.png) [@Yuricst](https://discourse.julialang.org/u/Yuricst)\
**Post date:** [December 11, 2024, 7:28am UTC](https://discourse.julialang.org/t/nonlinear-program-with-jump-not-following-all-constraints/123691/1 "2024-12-11T07:28:34Z")

</div>

Hello!  
I am trying to solve a nonlinear nonconvex NLP using JuMP. More specifically, I would like to write a wrapper that accepts a somewhat generic objective/constraint evaluation function `f_fitness` and generates a JuMP model with memoization.

My current issue is that when I have multiple constraints, only the last constraint is followed (even though I see them registered within the `model`).  
For example, I am working with a toy problem with one equality constraint h and two inequality constraints g\_1 and g\_2:

\begin{aligned} \min\_{x,y} \quad& x^2 - y \\ \text{such that}\quad& h: x^3 + y - 2.4 = 0 \\& g\_1 : -0.3x + y - 2 \leq 0 \\& g\_2 : x + y - 5 \leq 0 \end{aligned}

Solving a JuMP model with memoization (please see example code at the end of this post), depending on the order in which h, g\_1, and g\_2 are registered to the `model`, the converged solution changes (in all cases, Ipopt terminates with `LOCALLY_SOLVED`):

 ![Screenshot 2024-12-11 at 16.22.24](https://global.discourse-cdn.com/julialang/original/3X/2/9/29e0b0fef1661aba90a2399ed55961e6d6dcaff4.png)  
 ![Screenshot 2024-12-11 at 16.22.46](https://global.discourse-cdn.com/julialang/original/3X/7/7/77689a28fbcc2b4a602ec57626108592161f1783.png)  
 ![Screenshot 2024-12-11 at 16.22.37](https://global.discourse-cdn.com/julialang/original/3X/a/b/ab5512cff2a641de6a44035759630c93612c30e4.png)

I imagine that there is some silly mistake I am missing/perhaps misunderstanding how the `add_nonlinear_operator` function works, but I would appreciate any help/guidance.

Here’s a minimum working example (below, the constraints are registered in order: h, then g\_1, then g\_2, so only g\_2 is satisfied):

```julia
using GLMakie
using Ipopt
using JuMP

function f_fitness(x::T...) where {T<:Real}
	# objective
    f = x[1]^2 - x[2]
    
    # equality constraints
    h = zeros(T, 1)
    h = x[1]^3 + x[2] - 2.4

    # inequality constraints
    g = zeros(T, 2)
    g[1] = -0.3x[1] + x[2] - 2 # y <= 0.3x + 2
    g[2] = x[1] + x[2] - 5 # y <= -x + 5
    
    fitness = [f; h; g]
    return fitness
end

function memoize(foo::Function, n_outputs::Int)
    last_x, last_f = nothing, nothing
    last_dx, last_dfdx = nothing, nothing
    function foo_i(i, x::T...) where {T<:Real}
        if T == Float64
            if x !== last_x
                last_x, last_f = x, foo(x...)
            end
            return last_f[i]::T
        else
            if x !== last_dx
                last_dx, last_dfdx = x, foo(x...)
            end
            return last_dfdx[i]::T
        end
    end
    return [(x...) -> foo_i(i, x...) for i in 1:n_outputs]
end

# problem dimensions
nx = 2 # number of decision vectors
nh = 1 # number of equality constraints
ng = 2 # number of inequality constraints
lx = -10*ones(nx,)
ux = 10*ones(nx,)
x0 = [-1.2, 10]

# get model
order = 2
diff_f = "forward"
model = Model(Ipopt.Optimizer)

# generate memoized fitness function 
nfitness = 1 + nh + ng
memoized_fitness = memoize(f_fitness, nfitness)

# set variables
@variable(model, lx[i] <= x[i in 1:nx] <= ux[i], start = x0[i])

# set objective
@operator(model, fobj, nx, memoized_fitness[1])
@objective(model, Min, fobj(x...))

# append equality constraints
h_ops = add_nonlinear_operator.(model, nx, memoized_fitness[2:1+nh])
@constraint(model, hs[ih in 1:nh], h_ops[ih](x...) == 0)

# append inequality constraints
g_ops = add_nonlinear_operator.(model, nx, memoized_fitness[2+nh:end])
@constraint(model, gs[ig in 1:ng], g_ops[ig](x...) <= 0)

# run optimizer
optimize!(model)

# checks
@assert is_solved_and_feasible(model)
println(termination_status(model))
println("Decision vector: $(value.(model[:x]))")
println("Objective: $(objective_value(model))")

# plot contour
fig = Figure(size=(500,500))
ax = Axis(fig[1,1], xlabel="x1", ylabel="x2")
xs_grid = LinRange(-11, 11, 100)
ys_grid = LinRange(-11, 11, 100)
contourf!(ax, xs_grid, ys_grid, (x, y) -> f_fitness(x, y)[1], levels=20)
lines!(ax, [x for x in xs_grid], [2.4 - x^3 for x in xs_grid], color=:blue)
fill_between!(ax, xs_grid, maximum(xs_grid) * ones(length(xs_grid)), 
    [0.3x + 2 for x in xs_grid], color=:black, alpha=0.35)
fill_between!(ax, xs_grid, maximum(xs_grid) * ones(length(xs_grid)), 
    [5 - x for x in xs_grid], color=:black, alpha=0.35)
scatter!(ax, [value(model[:x][1])], [value(model[:x][2])], markersize=5, color=:red)
xlims!(ax, minimum(xs_grid), maximum(xs_grid))
ylims!(ax, minimum(ys_grid), maximum(ys_grid))
display(fig)

```

Thank you!

---

<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:** [December 11, 2024, 8:13am UTC](https://discourse.julialang.org/t/nonlinear-program-with-jump-not-following-all-constraints/123691/2 "2024-12-11T08:13:44Z")

</div>

This is a very subtle issue, and it took me a while to figure out! I was getting very confused.

You need to give each operator a unique name:

```julia
for i in 1:nh
    op = add_nonlinear_operator(
        model,
        nx,
        memoized_fitness[1+i];
        name = Symbol("op_h_$i"),
    )
    @constraint(model, op(x...) == 0)
end
for i in 1:ng
    op = add_nonlinear_operator(
        model,
        nx,
        memoized_fitness[1+nh+i];
        name = Symbol("op_g_$i"),
    )
    @constraint(model, op(x...) <= 0)
end

```

Now, this should happen by default with `add_nonlinear_operator`, because the `name` argument is by default inferred from `Symbol(f)` where `f` is the function. Typically, this is sufficient:

```julia
help?> add_nonlinear_operator
search: add_nonlinear_operator add_nonlinear_constraint add_nonlinear_expression

  add_nonlinear_operator(
      model::Model,
      dim::Int,
      f::Function,
      [∇f::Function,]
      [∇²f::Function];
      [name::Symbol = Symbol(f),]
  )

```

However, the anonymous functions returned by `memoize` all have the same name for some reason!

```Julia
julia> memoized_fitness
4-element Vector{var"#333#336"{Int64, var"#foo_i#334"{typeof(f_fitness)}}}:
 #333 (generic function with 1 method)
 #333 (generic function with 1 method)
 #333 (generic function with 1 method)
 #333 (generic function with 1 method)

```

I didn’t expect this. I would have assumed that they should have had unique names.

I wonder if we should throw a warning or error when defining a name that already exists.

---

<div class="post-metadata">

**Author:** ![Yuricst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuricst/32/207613_2.png) [@Yuricst](https://discourse.julialang.org/u/Yuricst)\
**Post date:** [December 11, 2024, 3:09pm UTC](https://discourse.julialang.org/t/nonlinear-program-with-jump-not-following-all-constraints/123691/3 "2024-12-11T15:09:20Z")

</div>

Thank you so much - I was stuck with this for quite a while, so this helps me tremendously!