Adding NL expressions to an AffExpr type Vector

Specific Domains Optimization (Mathematical)
,
1

Hello again everybody,

I’m trying to generate a Vector where which element is a NL expression, here is the code where I’m trying to implement it.

register(rocket2d_hs_2, :StateFunc, 10, StateFunc; autodiff=true)
x_c = Array{AffExpr}(undef,7,n-1)

for j in 1:n-1
    for i in 1:7
        x_c[i,j] = AffExpr(0.0)
        if i<6
            val = value(0.5*(X[i,j]+X[i,j+1]) + (Δt/8)*(StateFunc(g_0,c,i,X[:,j]...)-StateFunc(g_0,c,i,X[:,j+1]...)))
        else
            val = (X[i,j]+X[i,j+1])/2
        end
        add_to_expression!(x_c[i,j],val)
    end
    for i in 1:5
        @NLconstraints(rocket2d_hs_2,begin
                X[i,j+1]-X[i,j] == (Δt/6)*(StateFunc(g_0,c,i,X[:,j]...) + StateFunc(g_0,c,i,x_c[:,j]...) + StateFunc(g_0,c,i,X[:,j+1]...))
            end)
    end
end

Where X (7xn dim Array) is defined as a @variable and StateFunc is a user.defined function.
I get the following error:

ERROR: MethodError: no method matching /(::NonlinearExpression, ::Int64)
Closest candidates are:
  /(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}, ::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8}) at int.jl:93
  /(::StridedArray{P, N} where N, ::Real) where P<:Dates.Period at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\Dates\src\deprecated.jl:44
  /(::Union{SparseArrays.SparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}}, false} where var"#s814"<:SparseArrays.AbstractSparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s814", Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false} where var"#s814"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}} where {Tv, Ti}, ::Number) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SparseArrays\src\sparsevector.jl:1450

The error occurs when trying to evaluate “Val =” for the first time.

Maybe I’m not perfectly understanding the syntaxis for an AffExpr Array type given in the documentation. In case it’s not posible, there is another way to do something similar?

Thanks for the help and sorry for any inconvenience!

2

Hello again :smile: this question is follow-up of How to create an user-defined function with multiple vector as inputs.

A few things:

  1. AffExpr is a type for affine expressions. You cannot add nonlinear terms to it
  2. You cannot call a user-defined function outside of a macro
  3. You cannot call value before the model has been solved

It’s hard to offer complete advice without a reproducible example that I can copy-and-paste (what is StateFunc, what are X, g_0, c, and Δt?), but this should point you in the right direction:

register(rocket2d_hs_2, :StateFunc, 10, StateFunc; autodiff=true)

for j in 1:n-1
    x_c = map(1:7) do i
        if i < 6
            @NLexpression(
                rocket2d_hs_2,
                0.5*(X[i,j]+X[i,j+1]) + (Δt/8)*(StateFunc(g_0,c,i,X[:,j]...)-StateFunc(g_0,c,i,X[:,j+1]...)),
            )
        else
            @NLexpression(rocket2d_hs_2, (X[i,j]+X[i,j+1])/2)
        end
    end
    @NLconstraint(
        rocket2d_hs_2,
        [i in 1:5],
        X[i,j+1] - X[i,j] == (Δt/6) * (StateFunc(g_0,c,i,X[:,j]...) + StateFunc(g_0,c,i,x_c...) + StateFunc(g_0,c,i,X[:,j+1]...))
    )
end
3

Thank you for your time, what does map(1:7) exactly do?
Also, answering your questions, here is the full code for better context:

using JuMP
import Ipopt
import Plots
import LinearAlgebra
rocket2d_hs_2 = Model(Ipopt.Optimizer)
set_attribute(rocket2d_hs_2, "print_level", 4)
set_attribute(rocket2d_hs_2, "print_timing_statistics","yes")

#----Variables----#

# Initial conditions
m_0 = 30 #Kg
u_0 = 0 #m/s   Initial velocity x-axis
v_0 = 0 #m/s   Initial velocity y-axis
x_0 = 1000 # m  
y_0 = 100 # m

# Other conditions
g_0 = 9.81 #m/s^2
c = 2000
T_max = 1000 #N

#Final conditions
x_f=0
y_f=0
u_f=0
v_f=0

n = 300 #time steps

#---- variables-----#
@variables(rocket2d_hs_2, begin
      1 ≤ t_f ≤ 100
      X[1:7 , 1:n]
      
end)
@constraints(rocket2d_hs_2, begin
      0 .<= X[1,1:n]
      0 .<= X[2,1:n]
      0 .<= X[5,1:n]
      0 .<= X[6,1:n] .<= T_max
     deg2rad(-20) .≤ X[7,1:n] .≤ deg2rad(20)
end)
#----Initial conditions----#
fix(X[1,1], x_0; force=true)
fix(X[2,1], x_0; force=true)
fix(X[3,1], u_0; force=true)
fix(X[4,1], v_0; force=true)
fix(X[5,1], m_0; force=true)

#----Target conditions----#
fix(X[1,n], x_f; force=true)
fix(X[2,n], y_f; force=true)
fix(X[3,n], u_f; force=true)
fix(X[4,n], v_f; force=true)
@NLexpression(rocket2d_hs_2, Δt, t_f/(n-1))

function StateFunc(g_0, c, i, X...)
    if i == 1
        δ = sum(X[3])
    elseif i == 2
        δ = sum(X[4])
    elseif i==3
        δ = sum((X[6]*sin(X[7]))/X[6])
    elseif i==4
        δ = sum(-g_0 + (X[6]*cos(X[7]))/X[6])
    elseif i==5
        δ = sum(-sqrt((X[6]*sin(X[7]))^2 + (X[6]*cos(X[7]))^2 )/c)
    end
    return δ
end
register(rocket2d_hs_2, :StateFunc, 10, StateFunc; autodiff=true)
x_c = Array{AffExpr}(undef,7,n-1)

for j in 1:n-1
    for i in 1:7
        x_c[i,j] = AffExpr(0.0)
        if i<6
            val = value(0.5*(X[i,j]+X[i,j+1]) + (Δt/8)*(StateFunc(g_0,c,i,X[:,j]...)-StateFunc(g_0,c,i,X[:,j+1]...)))
        else
            val = (X[i,j]+X[i,j+1])/2
        end
        add_to_expression!(x_c[i,j],val)
    end
    for i in 1:5
        @NLconstraints(rocket2d_hs_2,begin
                X[i,j+1]-X[i,j] == (Δt/6)*(StateFunc(g_0,c,i,X[:,j]...) + StateFunc(g_0,c,i,x_c[:,j]...) + StateFunc(g_0,c,i,X[:,j+1]...))
            end)
    end
end
@objective(rocket2d_hs_2, Min, t_f)
optimize!(rocket2d_hs_2)
println("Tf optimo es: ",objective_value(rocket2d_hs_2))
4

It’s not map(1:7), it’s map(1:7) do i, see:

https://docs.julialang.org/en/v1/manual/functions/#Do-Block-Syntax-for-Function-Arguments

The do x syntax creates an anonymous function with argument x and passes it as the first argument to map.

2 Likes
5

I would write your model like

using JuMP
import Ipopt
# Number of time steps
n = 300
# Initial contitions
x_0, y_0, u_0, v_0, m_0 = 1000, 100, 0, 0, 30
# Target conditions
x_n, y_n, u_n, v_n = 0, 0, 0, 0
# State lower bounds
lb = [0.0, 0.0, -Inf, -Inf, 0.0, 0.0, deg2rad(-20)]
# State upper bounds
ub = [Inf, Inf, Inf, Inf, Inf, 1000, deg2rad(20)]

g_0 = 9.81  # m/s^2
c = 2000

function state_expression(i::Int, X::Vector)
    if i == 1
        return X[3]
    elseif i == 2
        return X[4]
    elseif i == 3
        return @NLexpression(model, X[6] * sin(X[7]) / X[6])
    elseif i == 4
        return @NLexpression(model, -g_0 + X[6] * cos(X[7]) / X[6])
    else
        @assert i == 5
        return @NLexpression(model, -sqrt((X[6] * sin(X[7]))^2 + (X[6] * cos(X[7]))^2) / c)
    end
end

model = Model(Ipopt.Optimizer)
@variable(model, 1 <= t_f <= 100)
@variable(model, lb[i] <= X[i = 1:7 , 1:n] <= ub[i])
fix.(X[1:5, 1], [x_0, y_0, u_0, v_0, m_0]; force = true)
fix.(X[1:4, n], [x_n, y_n, u_n, v_n]; force = true)

state_func_x = [state_expression(i, X[:, j]) for i in 1:5, j in 1:n]
@expression(model, Δt, t_f / (n - 1))

for j in 1:n-1
    x_c = map(1:7) do i
        if i <= 5
            return @NLexpression(
                model,
                (X[i,j] + X[i, j+1]) / 2 + (Δt / 8) * (state_func_x[i, j] - state_func_x[i, j+1]),
            )
        else
            return (X[i, j] + X[i, j+1]) / 2
        end
    end
    state_func_c = [state_expression(i, x_c) for i in 1:5]
    @NLconstraint(
        model,
        [i in 1:5],
        X[i,j+1] - X[i,j] == Δt / 6 * (state_func_x[i, j] + state_func_c[i] + state_func_x[i, j+1])
    )
end
@objective(model, Min, t_f)
optimize!(model)
objective_value(model)

But I may have made a typo somewhere, because Ipopt converges to a locally infeasible point.

The clunkiness of the current nonlinear interface in JuMP is a known problem. We’re in the middle of rewriting the entire nonlinear API that should make this a lot nicer in the future: https://github.com/jump-dev/JuMP.jl/pull/3106.

1 Like
6

Thank you very much Odow and Nils, everything is clearer now!