# Nonlinear JuMP - vector of expressions+differences between @NL and no @NL

**URL:** <https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [August 2, 2023, 4:20pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405 "2023-08-02T16:20:34Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 2, 2023, 4:20pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/1 "2023-08-02T16:20:34Z")

</div>

Following JuliaCon, I wanted to try out the nonlinear version of JuMP, especially the option of having vectors instead of a list of expressions. In the process I managed to

1. Make the nonlinear version work without the `@NL` macros, but still with a list of expressions,
2. Get three different errors: a classic `no method matching +(::VariableRef, ::Vector{AbstractJuMPScalar})`, an Ipopt error `EXIT: Problem has too few degrees of freedom.`, and `no method matching zero(::Type{Vector{AbstractJuMPScalar}})` if I tried to make a vector
3. Get a different number of iterations for the same problem using @NL macros and the new version, and imposing constraints directly in `@variable` or by `@constraint`.

Here we go.

Technical: I used `od/nlp-expr` `(not sure if correct because I cannot remember which one was mentioned during JuliaCon 😅)

Optimization problem (artificial):

\min z

s.t.  
(x\_1(0)-x\_{10})^2+(x\_2(0)-x\_{20})^2 \leq z, \forall \mu  
\dot{x}\_1(t)= x\_2(t)  
\dot{x}\_2(t)=\mu(1-x^2\_1(t))x\_2(t)-x\_1(t)  
where \mu=1,2, t=0,0.01,\ldots,20, x\_{10}=2, x\_{20}=0. Since we have two values of \mu, let’s say \mu\_1, \mu\_2, we end up with an extended formulation:  
\min z

s.t.  
(x\_{1,\mu\_1}(0)-x\_{10})^2+(x\_{2,\mu\_1}(0)-x\_{20})^2 \leq z  
(x\_{1,\mu\_2}(0)-x\_{10})^2+(x\_{2,\mu\_2}(0)-x\_{20})^2 \leq z  
\dot{x}\_{1,\mu\_1}(t)= x\_{2,\mu\_1}(t)  
\dot{x}\_{2,\mu\_1}(t)=\mu\_1(1-x\_{1,\mu\_1}^2(t))x\_{2,\mu\_1}(t)-x\_{1,\mu\_1}(t)  
\dot{x}\_{1,\mu\_2}(t)= x\_{2,\mu\_2}(t)  
\dot{x}\_{2,\mu\_2}(t)=\mu\_2(1-x^2\_{1,\mu\_2}(t))x\_{2,\mu\_2}(t)-x\_{1,\mu\_2}(t)

The main idea is put the dynamics into a three-dimensional array indexed by time (index `j=1:2001`), version of \mu (index `k=1:2`), and number of equation (index `m=1:2`). It is indexing by the number of equation that fails.

First a version that works, without `@NL` macros:

> **Summary**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 #Final time, exact value doesn't matter
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE - ideally I would like to avoid having every equation separately, but this works
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP, zVar)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> plot(value.(x)[:,1,2],value.(x)[:,2,2])
> plot!(value.(x)[:,1,2],value.(x)[:,2,2])
> 
> ```

The three different errors (I marked in each formulation what is different):

> **Summary**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 #Final time, exact value doesn't matter
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> # ######################################Doesn't work: no method matching +(::VariableRef, ::Vector{AbstractJuMPScalar})
> # @expressions(
> # vdP,
> # begin
> # #initial conditions
> # init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> # #RHS of ODE
> ##########CHANGE BEGINS HERE - one big vector of RHS
> # x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> # x[j,2,k]
> # scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> # ]
> ##########CHANGE ENDS HERE
> # end
> # )
> 
> ##########CHANGE BEGINS HERE - everything is indexed together by j, m, k
> # @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == x[j, m, k] + Δt * x_ODE[j,m, k])
> ##########CHANGE ENDS HERE
> 
> # ######################################
> 
> ######################################Doesn't work: Ipopt error:
> # Exception of type: TOO_FEW_DOF in file "Interfaces/IpIpoptApplication.cpp" at line 655:
> # Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!
> # EXIT: Problem has too few degrees of freedom.
> # @expressions(
> # vdP,
> # begin
> # #initial conditions
> # init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> # #RHS of ODE
> # x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> # x[j,2,k]
> # scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> # ]
> # end
> # )
> ##########CHANGE BEGINS HERE - broadcasting because the error said + couldn't work
> # @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] .== x[j, m, k] .+ Δt * x_ODE[j,m, k])
> ##########CHANGE ENDS HERE
> 
> ######################################
> 
> # ######################################Doesn't work: no method matching zero(::Type{Vector{AbstractJuMPScalar}})
> # @expressions(
> # vdP,
> # begin
> # #initial conditions
> # init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> # #RHS of ODE
> # x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> # x[j,2,k]
> # scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> # ]
> # end
> # )
> ##########CHANGE BEGINS HERE - no broadcasting, but I wanted to reproduce the error from a more complex example so I added constant terms to RHS: 1-1
> # @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == 1.0-1.0+x[j, m, k] + Δt * x_ODE[j,m, k])
> ##########CHANGE ENDS HERE
> 
> # ######################################
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP, zVar >= 0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> ```

And finally, the same problem implemented with @NL macros:

> **Summary**
>
> ```julia
> ######################################WORKS
> 
> tf = 20 #Final time, exact value doesn't matter
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @NLexpressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @NLconstraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @NLconstraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @NLexpression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP, zVar)
> @NLconstraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @NLobjective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> plot!(value.(x)[:,1,2],value.(x)[:,2,2],ls=:dash)
> plot!(value.(x)[:,1,2],value.(x)[:,2,2],ls=:dash)
> 
> ```

If I do not have any constraints on `z`, then the new version needs 146 iterations whereas `@NL` needs 252 iterations. The difference in iterations changes in the other direction if I impose an additional constraint on `z`: `@constraint(vdP,z>=0)`. Then the new version needs 603 iterations whereas `@NL` version needs 461 iterations. If I impose the constraint on `z` directly when creating the variable `@variable(vdP,zVar>=0)`, then the new version needs 1266 iterations whereas `@NL` version needs 442 iterations. I get that the number of iterations may change with the number of constraints, so the differences between no constraints and constraints make sense. I am not sure why I get a difference for two different formulations of the same constraint.

As a last thing, I noticed that this works with/without broadcast . after 1.0:

```julia
@expressions(
    vdP,
    begin
        #initial conditions
        init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
        #RHS of ODE
        x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
        x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0.-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k] ####<- 
    end
)

```

but this works only without broadcast . after 1.0:

```julia
@NLexpressions(
    vdP,
    begin
        #initial conditions
        init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
        #RHS of ODE
        x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
        x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k] ####<- 

    end
)

```

With `.`, JuMP throws an error `Unrecognized function ".-" used in nonlinear expression.` I guess this is related to the rule of “no vectors in nonlinear expressions”, right?

Questions:

1. Is it possible to avoid writing out the equations one by one?
2. What’s different in the `@NL` formulation that makes the number of iterations different?
3. Why having a number in the constraint changes the error if there is already a number in the 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:** [August 2, 2023, 9:37pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/2 "2023-08-02T21:37:41Z")

</div>

I got a bit confused with your edits. Can you make explicit reproducible examples for the different errors?

I might have made a mistake, but here’s how I’d write your model:

```julia
using JuMP, Ipopt
tf, Δt = 20, 0.01
n = round(Int, tf / Δt + 1)
μ = [1.0, 2.0]
K = length(μ)
x_0 = [2.0, 0.0]
model = Model(Ipopt.Optimizer)
@variable(model, x[1:n, 1:2, 1:K])
@variable(model, obj)
function ODE(x, j, m, k)
    if m == 1
        return x[j, 2, k]
    else
        return μ[k] * (1.0 - x[j,1,k] * x[j,1,k]) * x[j,2,k] - x[j,1,k]
    end
end
@constraints(model, begin
    [j = 1:n-1, m = 1:2, k = 1:K], x[j+1, m, k] == x[j, m, k] + Δt * ODE(x, j, m, k)
    [k = 1:K], sum((x[1, m, k] - x_0[m])^2 for m in 1:2) <= obj
end)
@objective(model, Min, obj)
optimize!(model)

```

## No method matching

The no method matching seems like a typical broadcasting error that isn’t specific to JuMP.

```julia
julia> x = [1, 2];

julia> 2 + x
ERROR: MethodError: no method matching +(::Int64, ::Vector{Int64})

```

## Too few degrees of freedom is the

Watch this talk from JuMP-dev:

[![](https://global.discourse-cdn.com/julialang/original/3X/a/0/a0836e3e5ff421f4745d51a4b3ffeff4c0044e56.jpeg "JuliaCon 2023 Room 32-155 Day 3") ](https://www.youtube.com/watch?v=DsC_gb9MFMs&t=18545)

## zero issue

I’d need to dig into this one a bit.

## Broadcast thing

`@NL` expressions support only a limited syntax. They’ll be going away. You cannot broadcast in `@NL` because they don’t support vectors.

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 3, 2023, 7:14am UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/3 "2023-08-03T07:14:33Z")

</div>

Thank you for the detailed response, including the suggestion for the code. I have rarely used functions inside constraints to avoid the whole registration process. But I will have a look now.

**Errors**  
As for the three errors, sorry for getting it all messed up, here are the three examples separately:

Error 1: `MethodError: no method matching zero(::Type{Vector{AbstractJuMPScalar}})`. This error appears if I have `1.0-1.0` in the constraint

> **Summary Error 1**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] 
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> x[j,2,k]
> scenariosMu[k]*(1.0.-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> ]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == 1.0-1.0+x[j, m, k] + Δt * x_ODE[j,m, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> ```

Error 2: no method matching +(::VariableRef, ::Vector{AbstractJuMPScalar}). Yes, this error is not specific to JuMP, but why does it appear here? My `x[j,m,k]` is a scalar and `x_ODE[j,m, k]` is (supposed to be) a scalar as well.

> **Summary Error 2**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> x[j,2,k]
> scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> ]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == x[j, m, k] + Δt * x_ODE[j,m, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar)
> # @constraint(vdP,zVar>=0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> ```

Error 3: Ipopt fails with `EXIT: Problem has too few degrees of freedom.`. Thank you for the video, I couldn’t find the recording from that room 🙂 In my case, I get a different description from ipopt (results from setting the print level to 12):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/0/2049dc73521cde1aca043030f75124ec753bb05e.png)  
I also know that my problem is feasible and has the right number of degrees of freedom if I write the equations separately. So here the error suggests that my attempted “vectorization” messes up with the number of constraints/variables.

> **Summary Error 3**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0]
> noScenarios = length(scenariosMu)
> vdP = Model(Ipopt.Optimizer)
> set_optimizer_attributes(
> vdP,
> "print_level" => 12,
> )
> 
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
> x[j,2,k]
> scenariosMu[k]*(1.0.-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> ]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] .== x[j, m, k] .+ Δt * x_ODE[j,m, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar)
> # @constraint(vdP,zVar>=0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> ```

**Different number of iterations for different implementations:**  
Case 1. No constraints on `z`:

- W/o `@NL`: 146 iterations
- W/ `@NL`: 252 iterations

> **Summary Case 1**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] 
> noScenarios = length(scenariosMu)
> 
> ####No @NL
> vdP = Model(Ipopt.Optimizer)
> 
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar)
> # @constraint(vdP,zVar>=0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> #########@NL
> vdPnl = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdPnl, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @NLexpressions(
> vdPnl,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @NLexpression(vdPnl, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdPnl, zVar)
> # @constraint(vdPnl,zVar>=0)
> @NLconstraint(vdPnl,[k=1:noScenarios],objMin[k]<=zVar)
> @NLobjective(vdPnl, Min, zVar)
> 
> optimize!(vdPnl)
> 
> ```

Case 2. Constraint imposed when creating `z`:

- W/o `@NL`: 1266 iterations
- W/ `@NL`: 442 iterations

> **Summary Case 2**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> #########No @NL
> vdP = Model(Ipopt.Optimizer)
> 
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar>=0)
> # @constraint(vdP,zVar>=0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> #########@NL
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] 
> noScenarios = length(scenariosMu)
> vdPnl = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdPnl, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @NLexpressions(
> vdPnl,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @NLexpression(vdPnl, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdPnl, zVar>=0)
> # @constraint(vdPnl,zVar>=0)
> @NLconstraint(vdPnl,[k=1:noScenarios],objMin[k]<=zVar)
> @NLobjective(vdPnl, Min, zVar)
> 
> optimize!(vdPnl)
> 
> ```

Case 3. Constraint on `z` imposed separately:

- W/o `@NL`: 603 iterations
- W/ `@NL`: 461 iterations

> **Summary Case 3**
>
> ```julia
> using JuMP, Ipopt, Plots
> 
> tf = 20 #Final time, exact value doesn't matter
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] #####Values don't matter
> noScenarios = length(scenariosMu)
> #########No @NL
> vdP = Model(Ipopt.Optimizer)
> 
> x_0 = [2.0, 0.0]
> @variables(vdP, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @expressions(
> vdP,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @expression(vdP, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdP,zVar)
> @constraint(vdP,zVar>=0)
> @constraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)
> @objective(vdP, Min, zVar)
> 
> optimize!(vdP)
> 
> #########@NL
> tf = 20 
> 
> Δt = 0.01
> n = Int64(round(tf / Δt + 1))
> scenariosMu = [1.0 2.0] 
> noScenarios = length(scenariosMu)
> vdPnl = Model(Ipopt.Optimizer)
> x_0 = [2.0, 0.0]
> @variables(vdPnl, 
> begin
> x[j=1:n, m=1:2, k=1:noScenarios] 
> end
> )
> 
> @NLexpressions(
> vdPnl,
> begin
> #initial conditions
> init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
> #RHS of ODE
> x1_ODE[j=1:n-1, k=1:noScenarios], x[j,2,k]
> x2_ODE[j=1:n-1, k=1:noScenarios], scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
> end
> )
> 
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 1, k] == x[j, 1, k] + Δt * x1_ODE[j, k])
> @NLconstraint(vdPnl,[j=1:n-1,k=1:noScenarios], x[j+1, 2, k] == x[j, 2, k] + Δt * x2_ODE[j, k])
> 
> @NLexpression(vdPnl, objMin[k=1:noScenarios], sum(init_x[m,k]*init_x[m,k] for m=1:2))
> @variable(vdPnl, zVar)
> @constraint(vdPnl,zVar>=0)
> @NLconstraint(vdPnl,[k=1:noScenarios],objMin[k]<=zVar)
> @NLobjective(vdPnl, Min, zVar)
> 
> optimize!(vdPnl)
> 
> ```

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 3, 2023, 7:20am UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/4 "2023-08-03T07:20:18Z")

</div>

Having written all this stuff, I realised what the error was, at least for Error 2 and Error 3, and have a guess for Error 1. As usual, the error was between the screen and the chair 🙃

For future reference:  
When I create an expression like this:

```julia
@expression(vdP,
        #RHS of ODE
        x_ODE[j=1:n-1,m=1:2, k=1:noScenarios], [
            x[j,2,k]
            scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
        ]
)

```

I am effectively saying that for every index `j=1:n-1, m=1:2, k=1:noScenarios`, my `x_ODE[j,m,k]` is a vector with two elements. For example:

![image](https://global.discourse-cdn.com/julialang/original/3X/f/2/f2d406bb1b5c0797d5e05b6ef8d2f6cfdbfc1c63.png)

This is why `@constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == x[j, m, k] + Δt * x_ODE[j,m, k])` gives me the error `no method matching +(::VariableRef, ::Vector{AbstractJuMPScalar})`. My `x[j,m,k]` is a scalar (correct), but my `x_ODE[j,m, k]` is a vector (incorrect). Error 2 sorted.

If I do broadcasting `@constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] .== x[j, m, k] .+ Δt * x_ODE[j,m, k])` I assign two different equality constraints to the same `x[j,m,k]` which corresponds to overconstraining my problem. That’s why too few degrees of freedom. Error 3 sorted.

That leaves Error 1. I guess this thread explains this [issue](https://discourse.julialang.org/t/using-literal-1-0-in-jump-constraint-throws-an-error/98496). As my RHS `@constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == 1.0-1.0+ x[j, m, k] + Δt * x_ODE[j,m, k])` is a vector (incorrectly) of abstract elements, the value `1.0-1.0` should also be a vector, but that doesn’t work, as indicated in the thread linked.

An ugly hack to sort out my incorrect indexing is to use:  
`@constraint(vdP,[j=1:n-1,m=1:2,k=1:noScenarios], x[j+1, m, k] == x[j, m, k] + Δt * x_ODE[j,1, k][m])`

Or actually write the expression with correct indexing:

```julia
@expressions(
    vdP,
    begin
        #initial conditions
        init_x[m=1:2, k=1:noScenarios], x[1, m, k] - x_0[m]
        #RHS of ODE - no need to have another index
        x_ODE[j=1:n-1,k=1:noScenarios], [
            x[j,2,k]
            scenariosMu[k]*(1.0-x[j,1,k]*x[j,1,k])*x[j,2,k]-x[j,1,k]
        ]
    end
)

@constraint(vdP,[j=1:n-1,k=1:noScenarios], x[j+1, 1:2, k] .== x[j, 1:2, k] .+ Δt * x_ODE[j,k])

```

That leaves only the different number of iterations in different formulations.

---

<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:** [August 3, 2023, 8:13pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/5 "2023-08-03T20:13:56Z")

</div>

Did you try my formulation above? You can use the function tracing to simplify things.

For the different iteration counts, it is expected that the problems are not identical. We’ve changed the way that some things are represented internally, so it may lead to a different ordering of constraints, or a slightly different expression graph of the constraints. For a convex problem, they should find the same optimal solution, but they might take different iterations to get there. It seems like your problem is non-convex, so anything could happen, including finding a different local minima.

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 4, 2023, 6:37am UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/6 "2023-08-04T06:37:58Z")

</div>

Thank you. These equations have very little to do with my actual problem, and I am happy with local solutions. Is there a way I can read more about the different formulations? I am intrigued by why imposing a constraint when creating a variable results in twice as many iterations.

The code you gave me runs fine! And after a minor adjustment I can make it work for my case (I wanted to avoid listing the equations indexed by m):

```julia
using JuMP, Ipopt
tf, Δt = 20, 0.01
n = round(Int, tf / Δt + 1)
μ = [1.0, 2.0]
K = length(μ)
x_0 = [2.0, 0.0]

# function ODE(x, j, m, k)
# if m == 1 ###I need to avoid writing rhs equations one by one
# return x[j, 2, k]
# else
# return μ[k] * (1.0 - x[j,1,k] * x[j,1,k]) * x[j,2,k] - x[j,1,k]
# end
# end

function ODE(x,μ,j,k)
    return [x[j, 2, k],μ[k] * (1.0 - x[j,1,k] * x[j,1,k]) * x[j,2,k] - x[j,1,k]] 
end

model = Model(Ipopt.Optimizer)
@variable(model, x[1:n, 1:2, 1:K])
@variable(model, obj)
@constraints(model, begin
    [j = 1:n-1, m = 1:2, k = 1:K], x[j+1, m, k] == x[j, m, k] + Δt * ODE(x,μ,j,k)[m]
    [k = 1:K], sum((x[1, m, k] - x_0[m])^2 for m in 1:2) <= obj
end)
@objective(model, Min, obj)
optimize!(model)

```

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [August 4, 2023, 3:52pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/7 "2023-08-04T15:52:11Z")

</div>

Hi @blob, I wanted to thank you for giving the non-linear system a try and reporting your experience in this forum. It is extremely helpful to us in understanding the unexpected challenges that new users may face. Thanks for your input!

---

<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:** [August 5, 2023, 11:19pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/8 "2023-08-05T23:19:09Z")

</div>

I think you’re after this formulation:

```julia
using JuMP, Ipopt
function ODE(x, μ, j, k)
    return [
        x[j, 2, k],
        μ[k] * (1 - x[j, 1, k] * x[j, 1, k]) * x[j, 2, k] - x[j, 1, k]
    ]
end
tf, Δt, μ, x_0 = 20, 0.01, [1.0, 2.0], [2.0, 0.0]
n = round(Int, tf / Δt + 1)
K = length(μ)
model = Model(Ipopt.Optimizer)
@variable(model, x[1:n, 1:2, 1:K])
@variable(model, obj)
@constraints(model, begin
    [j = 1:n-1, k = 1:K], x[j+1, :, k] .== x[j, :, k] .+ Δt * ODE(x, μ, j, k)
    [k = 1:K], sum((x[1, m, k] - x_0[m])^2 for m in 1:2) <= obj
end)
@objective(model, Min, obj)
optimize!(model)

```

Looking at the different number of iterations isn’t very useful. Because the problem is non-convex, small changes in the formulation can lead to dramatically different sequences of iterations.

There are small changes, like going from the nonlinear constraint `@NLconstraint(vdP,[k=1:noScenarios],objMin[k]<=zVar)` to the quadratic constraint using `@constraint`. In the later case, Ipopt exploits the fact that the Hessian is constant.

There are also minor parsing changes in how we construct the expressions. In the `@NL` case, you have something like `x[j+1, 1, k] - (x[j, 1, k] + Δt * x1_ODE[j, k]) == 0`. In the latter case, you end up with something like `(x[j+1, 1, k] - x[j, 1, k]) - Δt * x1_ODE[j, k] == 0`. A very minor change, but it could still lead to slightly different numerical solutions.

To improve performance, you should give all `x` an initial starting solution, and probably give good bounds on values of `x`.

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 7, 2023, 5:46am UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/9 "2023-08-07T05:46:07Z")

</div>

> Blockquote There are also minor parsing changes in how we construct the expressions. In the `@NL` case, you have something like `x[j+1, 1, k] - (x[j, 1, k] + Δt * x1_ODE[j, k]) == 0`. In the latter case, you end up with something like `(x[j+1, 1, k] - x[j, 1, k]) - Δt * x1_ODE[j, k] == 0`. A very minor change, but it could still lead to slightly different numerical solutions.

Wow, didn’t expect that such small changes may have an impact. Thank you!

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [August 7, 2023, 12:33pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/10 "2023-08-07T12:33:40Z")

</div>

Another follow-up, mostly regarding the recommendation

> probably give good bounds on values of `x`.

There are two (that I identified, may be more) sources of different number of iterations:

1. using `@NL` and using the new version → differences due most likely to different formulations as explained by @odow in this thread
2. setting the constraints in two different ways (in either version of JuMP):  
a. directly in `@variable`, e.g. `@variable(model,obj>=0)`  
b. separately as `@constraint` , e.g. `@variable(model,obj)`, `@constraint(model,obj>=0)`

Case 2 seems unrelated to the version of JuMP and the differences here are most likely due to:

1. The way MathOptInterface works with constraints because it uses different formulations (thus see Case 1): [Reference · MathOptInterface](https://jump.dev/MathOptInterface.jl/v0.9.3/apireference/#Variables-1)
2. The way solvers treat bounded variables and constraints. For ipopt: [julia - Why does IPOPT evaluate objective function despite breaching constraints? - Stack Overflow](https://stackoverflow.com/questions/36907064/why-does-ipopt-evaluate-objective-function-despite-breaching-constraints). The bounds cannot be relaxed (the constraints can), and thus the iterations are different. For MadNLP the differences coming from different formulations are even more pronounced: no constraints - works fine, bounded variable - restoration failed, free variable+constraint - runs out of 10000 iterations (but MadNLP didn’t work with the new nonlinear version of JuMP, so the comparison is messy).

I guess I just need to play with the options (or try and formulate my problem as unconstrained optimization 🙃)

---

<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:** [August 7, 2023, 8:19pm UTC](https://discourse.julialang.org/t/nonlinear-jump-vector-of-expressions-differences-between-nl-and-no-nl/102405/11 "2023-08-07T20:19:00Z")

</div>

> [@blob](#):
>
> Wow, didn’t expect that such small changes may have an impact

It might not have been that, there were other changes as well. But in general, the problem we are building is mathematically equivalent, but slightly different in code.

> Case 2 seems unrelated to the version of JuMP and the differences here are most likely due to

Correct. `@variable` adds a variable bound. `@constraint` adds a new row to the matrix.

You shouldn’t expect to get useful information out of analyzing the differences. Ipopt and MadNLP assume that the problem is convex. Your problem isn’t, so anything can happen.

Try adding bounds on `x`, not just on `obj`.