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)
