Constraints lead to invalid no. derivative in model

[GNCGenie]

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

[odow]

The replace doesn't do anything because it acts on the JuMP expressions, not the numeric values.

You need to formulate the problem so that you do not have the undefined values. For example, by adding appropriate bounds and scaling variables appropriately.

As one example, instead of @constraint(model, X[1,:] .>= 6.5e6) # Minimum altitude do

set_lower_bound.(X[1,:], 6.5e6) # Minimum altitude

This question might be better asked on the forum: https://jump.dev/forum. It is not a bug in Juniper.