Complex Dynamics Constraints

[tpresser570]

transfer_example.txt
Hello,
I am trying to solve the optimal control problem for a spacecraft with thrust in a two-body problem. This should closely follow the hovercraft examples but the optimizer fails. My objective function is to minimize the energy (full script below):

@objective(model, Min, ∫(u[1]^2 + u[2]^2 + u[3]^2, t))

My variables are:
@variables(model, begin # state variables x[1:3], Infinite(t) v[1:3], Infinite(t) # control variables u[1:3], Infinite(t) end)

I provide two waypoints, the start and end positions
@constraint(model, [i = 1:3, j = eachindex(tw)], x[i](tw[j]) == xw[i, j])

as well as use the two body dynamics as the constraints on the derivative:

#compute radius ^ 3/2 @expression(model, expr, ((x[1]^2+x[2]^2+x[3]^2)^(3/2))) #ODE equations @constraint(model, [i = 1:3], ∂(x[i], t) == v[i]) @constraint(model, [i = 1:3], ∂(v[i], t) == -μ*x[i]/expr + u[i])

But this results in IPOPT with a 0 Jacobian.

Number of nonzeros in equality constraint Jacobian...: 384
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 210

The Jacobian for the equality constraints contains an invalid number

Has anyone encountered this problem before??

[pulsipher]

Can you please provide a minimum working example (i.e., the full script to reproduce this)?

[tpresser570]

Sure the example transfer is Mars to asteroid Psyche with a constant mass spacecraft with unbounded thrust (.txt file at the bottom):

`using InfiniteOpt
using Ipopt

#----------- Model Setup Mars/Psyche Transfer-----------
#sun gravitational constant
μ = 1.3271645321e11

#mars inital position and velocity
X0_mars = [-5.372189524343949e7;
4.135366207409543e8;
-1.7926248044321448e7;
-18.23019366914193;
-0.014994849925857007;
0.49411154199716095] #[rx,ry,rz,vx,vy,Vz_constr] in km and km/s

#psyche final target position
Xf_psyche= [ -5.372189524343949e7;
4.135366207409543e8;
-1.7926248044321448e7] #[rx,ry,rz] in km

xw = [X0_mars[1] Xf_psyche[1];
X0_mars[2] Xf_psyche[2];
X0_mars[3] Xf_psyche[3]] # positions

tof = 1026 * 86400 #seconds from [NP]
tw = [0, tof] # times

model = InfiniteModel(Ipopt.Optimizer);

@infinite_parameter(model, t in [0, tw[2]], num_supports = 10,
derivative_method = OrthogonalCollocation(2))

X0 = X0_mars

#variables
@variables(model, begin
# state variables
x[1:3], Infinite(t)
v[1:3], Infinite(t)
# control variables
u[1:3], Infinite(t)
end)

#objective
@objective(model, Min, ∫(u[1]^2 + u[2]^2 + u[3]^2, t))

#intial conditions
@constraint(model, [i = 1:3], vi == X0[i])

#compute radius ^ 3/2
@expression(model, expr, ((x[1]^2+x[2]^2+x[3]^2)^(3/2)))
#ODE equations
@constraint(model, [i = 1:3], ∂(x[i], t) == v[i])
@constraint(model, [i = 1:3], ∂(v[i], t) == -μ*x[i]/expr + u[i])

@constraint(model, [i = 1:3, j = eachindex(tw)], xi == xw[i, j])

print(model)
optimize!(model)

x_opt = value.(x);

using Plots
scatter(xw[1,:], xw[2,:], xw[3,:], label = "Waypoints", background_color = :transparent)
plot!(x_opt[1], x_opt[2], x_opt[3], label = "Trajectory")
xlabel!("x_1")
ylabel!("x_2")`

[tpresser570]

This is the script as a txt file

transfer_example.txt

[pulsipher]

First off, welcome to the forum!

With respect to the formulation a few observations:

• The initial condition syntax should be @constraint(model, [i = 1:3], v[i](0) == X0[i]) to create point variables correctly
• The waypoint constraint should be @constraint(model, [i = 1:3, j = eachindex(tw)], x[i](0) == xw[i, j]) to create point variables correctly
• The inequality Jacobian has 0 nonzero elements because there are no inequality constraints, so this is expected (the equality constraint Jacobian does have a number of nonzeros that makes sense)
• The Invalid number in NLP function or derivative detected. error that occurs with the corrected code occurs because x can take on 0 or negative values which is mathematically undefined for (x[1]^2+x[2]^2+x[3]^2)^(3/2)
• This model is poorly scaled with numbers near 0 and others many orders of magnitude larger

To address the above concerns, I would recommend the following:

• Rescaling the constants to be less extreme in their differences if at all possible (I know inter-plant distances are large though)
• Reformulate the equations to avoid invalid values or add bounds to the variables to avoid incurring problematic values
• Adding reasonable bounds to the variables can be quite beneficial in solving the problem with Ipopt
• Nonlinear models like this one typically benefit a lot from good guess values/trajectories for the variables (as discussed in Minimum time manouevering problem + boundaries #219)
• Once, the model is working, you will want to use more than just 10 support points (that is a very course discretization)