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

Specific Domains Optimization (Mathematical)
4

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 :upside_down_face:

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

@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

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. 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:

@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.