Constraints lead to invalid no. derivative in model

I am using Juniper to solve an orbital insertion problem for a LEO satellite. Treating it as a non-linear and non-convex optimisation problem. However introducing the first-order hold dynamics constraints for orbital mechanics leads to the problem giving the following error:

nl_solver   : MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 3])
mip_solver  : MathOptInterface.OptimizerWithAttributes(HiGHS.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("output_flag") => true])
log_levels  : [:Options, :Table, :Info]

#Variables: 46
#IntBinVar: 0
Obj Sense: Min

Total number of variables............................:       46
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       36
Total number of inequality constraints...............:        5
        inequality constraints with only lower bounds:        5
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

Number of Iterations....: 3000

                                   (scaled)                 (unscaled)
Objective...............:   1.4643330105627386e+10    1.4643330105627386e+10
Dual infeasibility......:   1.7583873920294525e+15    1.7583873920294525e+15
Constraint violation....:   6.4624646118998332e+06    6.4624646118998332e+06
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   5.3678697801909525e+00    5.3678697801909525e+00
Overall NLP error.......:   7.8302837315459702e+10    1.7583873920294525e+15

Number of objective function evaluations             = 4311
Number of objective gradient evaluations             = 22
Number of equality constraint evaluations            = 4333
Number of inequality constraint evaluations          = 4332
Number of equality constraint Jacobian evaluations   = 3021
Number of inequality constraint Jacobian evaluations = 3021
Number of Lagrangian Hessian evaluations             = 3000
Total seconds in IPOPT                               = 1.743

EXIT: Maximum Number of Iterations Exceeded.

Number of Iterations....: 0

Number of objective function evaluations             = 0
Number of objective gradient evaluations             = 0
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 1
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 0.000

EXIT: Invalid number in NLP function or derivative detected.
Status of relaxation: INVALID_MODEL

Pasting the code to reproduce the issue below:

begin
    using JuMP
    using Juniper,Ipopt, HiGHS

    R_EARTH     = 6.378136300e6           # [m] GGM05s Value
        GM_EARTH    = 3.986004415e14          # [m^3/s^2] GGM05s Value
        J2_EARTH    = 0.0010826358191967      # [] GGM05s value
        OMEGA_EARTH = 7.292115146706979e-5    
        ωₑ = OMEGA_EARTH;
        μ = GM_EARTH;
        Rₑ = R_EARTH;
        J₂ = J2_EARTH;

    # Coast phase state transition matrix
    function dT(x...)
        (a,e,i,ω,Ω,θ) = x
        p = a*(1-e^2)
        r = p/(1+e*cos(θ))
        h = a*√(μ/a)
        ωₛ = √(μ/a^3)

        A = replace([
                .0,
                .0,
                .0,
                .0,
                (-3/2)*(Rₑ^2/(a*(1-e^2))^2)*J₂*ωₛ*cos(i),
                h/r^2]
            ,NaN=>0, Inf=>0, -Inf=>0)

        return A
    end

    # Firing phase state transition matrix
    function dU(x...)
        (a,e,i,ω,Ω,θ) = x
        p = a*(1-e^2)
        r = p/(1+e*cos(θ))
        h = a*√(μ/a)
        ωₛ = √(μ/a^3)

        B = replace(
            [(2*a^2/h)*(e*sin(θ)) (2*a^2/h)*(p/r) 0
            (1/h)*(p*sin(θ)) (1/h)*((p+r)*cos(θ) + r*e) 0
            0 0 (r*cos(θ+ω))/h
            -p*cos(θ)/(e*h) (p+r)*sin(θ)/(e*h) -(r*sin(θ+ω)*cos(i))/(h*sin(i))
            0 0 (r*sin(θ+ω))/(h*sin(i))
            (p*cos(θ))/(e*h) -((p+r)*sin(θ))/(e*h) 0]
            , NaN=>0, Inf=>0, -Inf=>0)

        return B
    end

    h = 4
    n = 6
    m = 3

    ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3)
    highs = optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => true)
    model = Model(
        optimizer_with_attributes(
            Juniper.Optimizer,
            "nl_solver" => ipopt,
            "mip_solver" => highs,
        ),
    )

    @variables model begin
        X[1:6, 1:h+1]
        U[1:3, 1:h]
        T[1:h]
    end
    @constraint(model, X[1,:] .>= 6.5e6) # Minimum altitude

    for i=1:2:h # Coast Phase
        @constraint(model, X[:,i+1] .== dT(X[:,i]...)*T[i])
    end
    for i=2:2:h # Firing Phase
        @constraint(model, X[:,i+1] .== X[:,i] + dU(X[:,i]...)*U[:,i])
    end
    @constraint(model, X[:,begin] == [7e6, 1e-4, pi/3, pi/6, 0, 0]) # Initial boundary condition
    @constraint(model, X[:,end] == [7.01e6, 1e-4, pi/3, pi/7, 0, 0]) # Final boundary condition
    @objective(model, Min, sum(U.^2)) # Minimise fuel use
    JuMP.optimize!(model)
end

lanl-ansi / Juniper.jl

Constraints lead to invalid no. derivative in model #267