Adding NL expressions to an AffExpr type Vector

Specific Domains Optimization (Mathematical)
,
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