Trajectory optimisation for a spacecraft using IPOPT GEKKO

Question

So I am trying to solve a nonlinear differential optimisation problem in IPOPT GEKKO (python). I need to minimize time for a spacecrafts trajectory where equations of motion are described by circular restricted three-body dynamics. The trajectory is from earth to L1 between Sun-Earth-system. I have fixed initial and endpoint position and velocity coordinates and I need to integrate differential equations such that the time taken from earth to final location is minimized. In the dynamics of CRTBP I have added nondimensional thrust in x, y, z, direction and are the control variables of this optimisation problem. This leads to a path constraint such that at each thrusting segment the total thrust magnitude must be less then or equal sum of each thrust-components magnitude squared (see in the problem formulation image below).

I am quite new to both optimal control related problems and softwares like IPOPT and GEKKO. I have tried to code this problem in GEKKO using IPOPT as the optimisation software but it never seems to converge. I have tried to follow this trajectory optimisation example especially the 2D-part: https://transport-systems.imperial.ac.uk/tf/60008_21/n7_1_introduction_to_trajectory_optimisation.html

Obviously in my problem, I have more complicated dynamics as well as a 3D-trajectory optimisation. But considering I have fixed starting and end points and also given initial guesses for each variable (close to solution), in the end it is just a matter of finding the best route between A and B given three-body dynamics. Although I keep getting this error,

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

Here is my code:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
#initialize model
m = GEKKO()

# optional solver settings with APOPT
Nsim = 100 #steps when thrust is constant
# tf = 12 days /tau_star nondim

m.time = np.linspace(0, 0.206, Nsim)
#constants
mu = 3.003481e-6
l_star = 1.49597871e11
G_c = 6.67408 * 10**-11
m_star = 1.9885e30
tau_st = (l_star**3/G_c/m_star)**0.5 # circa 58 days
Isp = 3000 #  impulse per second
g0 = 9.80665* 10**-3 # km/s^2 => g = 9.8
T_max = 0.1 # {N, 100 mN}
m0 = 500 #kg
f = (T_max*(tau_st**2)/(l_star * m0)) #maximal nondim thrust
sub_L1 = 0.0158
x_final = sub_L1 - (1-mu)

#manipulating variables thrust components in x,y,z direction
#initial values represent full-thrust in x-direction as well as downwards
ux = m.MV(0.99*f, lb=-f, ub=f)
ux.STATUS = 1
uy = m.MV(0, lb=-f, ub=f)
uy.STATUS = 1
uz = m.MV(m.sqrt(1 - 0.99**2)*f, lb=-f, ub=f)
uz.STATUS = 1

#variables + initial guesses??
x1 = m.Var(value=x_final, lb=0, ub=1)
x2 = m.Var(value= -0.0023, lb=0, ub=1)
x3 = m.Var(value= 0, lb=0, ub=1)
x4 = m.Var(value= 0.3, lb=-1, ub=1)
x5 = m.Var(value= 0.3, lb=-0.5, ub=0.5)
x6 = m.Var(value= 0.3, lb=-0.5, ub=0.5)
#constraint on objetive function
# guess value for tf is 9 days meaning for an optimal path less time is needed to get to x_final
tf = m.FV(value = 0.15, lb=0.03, ub=0.5) # time upper and lower constraints
tf.STATUS = 1


#defining r1 and r2 as equations to be solved implicity together
a = m.Var(value=0.99)
b = m.Var(value=0.004)
m.Equation(((x1 - mu)**2 + x2**2 + x3**2)**(3/2) == a) #r1
m.Equation(((x1 + 1 - mu)**2 + x2**2 + x3**2)**(3/2) == b) #r2
#dynamics of three-body-problem
m.Equation( x1.dt() == x4*tf)
m.Equation( x2.dt() == x5*tf)
m.Equation( x3.dt() == x6*tf)
m.Equation( x4.dt() == tf*(2*x5 + x1 - ((1-mu)*(x1-mu)/a) - (mu*(x1+1-mu))/b + ux))
m.Equation( x5.dt() == tf*(-2*x4 + x2 - ((1-mu)*x2/a) - (mu*x2)/b + uy))
m.Equation( x6.dt() == tf*(-(1-mu) * x3 / a - (mu*x3)/b + uz))

#path constraints
m.Equation(x1 <= -0.1)
eq = m.Param(value=f)
m.Equation(ux**2 + uy**2 + uz**2 <= eq**2)
#starting constraints, starting x-value is from earth's escape
m.fix(x1, pos=0,val=-0.9952)
m.fix(x2, pos=0,val=-0.0023)
m.fix(x3, pos=0,val=-0.0010)
#m.fix(x4, pos=0,val=1.3890)
m.fix(x4, pos=0,val=0.002)
#m.fix(x5, pos=0,val=1.0585)
m.fix(x5, pos=0,val=0.001)
m.fix(x6, pos=0,val=0.0273)
#boundary constraints
m.fix(x1, pos=len(m.time)-1,val=x_final) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x2, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x3, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x4, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x5, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x6, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.Obj(tf) # minimize final time
m.options.IMODE = 6 # non-linar model
#m.options.NODES = 3 # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.options.MAX_ITER = 15000
m.options.RTOL = 1e-3
m.options.OTOL = 1e-3
m.solve() # Solve
print('Optimal time: ' + str(tf.value[0]))
m.solve(disp=True)
m.open_folder(infeasibilities.txt)

Answer 1

One infeasible condition is the equation m.Equation(x1 <= -0.1) but the variable is defined with lower bound of 0 and upper bound of 1 with x1 = m.Var(value=x_final, lb=0, ub=1). Try removing inequality constraints on the state variables until the problem becomes feasible. Successively add the constraints until it becomes infeasible again to identify which constraint(s) is (are) causing the infeasible solution.

Sometimes the problem is feasible but the solver can't find the solution or gets stuck in a non-feasible region. The end-point hard constraints can be challenging for solvers:

m.fix(x1, pos=len(m.time)-1,val=x_final) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x2, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x3, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x4, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x5, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x6, pos=len(m.time)-1,val=0.0) # stationary in sub-L1

One strategy is to combined soft (objective) constraints with the hard constraints to help guide the solver. Check out the strategy taken with the inverted pendulum problem.

It solves better with the combination such as:

# soft constraints
m.Minimize(final*ya**2)
m.Minimize(final*va**2)
m.Minimize(final*theta_a**2)
m.Minimize(final*qa**2)

# hard constraints
m.fix(ya,pos=end_loc,val=0.0)
m.fix(va,pos=end_loc,val=0.0)
m.fix(theta_a,pos=end_loc,val=0.0)
m.fix(qa,pos=end_loc,val=0.0)

Another problem that may be helpful is the HALE aircraft optimization problem.

• Martin, R.A., Gates, N., Ning, A., Hedengren, J.D., Dynamic Optimization of High-Altitude Solar Aircraft Trajectories Under Station-Keeping Constraints, Journal of Guidance, Control, and Dynamics, 2018, doi: 10.2514/1.G003737. Preprint | Article | Source Code