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

[blob]

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

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

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

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

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

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

[odow]

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:

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

JuliaCon 2023 Room 32-155 Day 3

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.

[blob]

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

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

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

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

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

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

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

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)

[blob]

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:

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

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.

[odow]

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.

[blob]

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

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)

[ccoffrin]

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!

[odow]

I think you’re after this formulation:

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.

[blob]

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!

[blob]

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
2. The way solvers treat bounded variables and constraints. For ipopt: julia - Why does IPOPT evaluate objective function despite breaching constraints? - Stack Overflow. 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 )

[odow]

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