With the following code in Python GEKKO, IPOPT throws an error:
Restoration phase is called at point that is almost feasible, with constraint violation 0.000000e+00. Abort.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math
#constants
mu = 3.98574405096E14
g = 9.81
R_E = 6.3781E6
J2 = 1.08262668E-3
P = 10E3
eta = 0.65
Isp = 3300.0
m0 = 1200.0
t_max = 86400.0 * 365 * 2
E_en = math.pi
E_ex = -math.pi
oe_i = np.array([R_E + 200000.0, 0, math.radians(28.5), math.radians(0.0), math.radians(0.0), math.radians(0.0)])
oe_f = np.array([R_E + 210000.0, 0, math.radians(28.5), math.radians(0.0), math.radians(0.0), math.radians(0.0)])
dm = (2 * eta * P)/((g * Isp)**2)
delta_t = 3600.0
#initialize model
traj = GEKKO()
nt = 200
traj.time = np.linspace(0, 1, nt)
#state variables
a = traj.SV(value = oe_i[0], lb = R_E, ub = oe_f[0] + 5000.0, name = 'sma')
e = traj.SV(value = oe_i[1], lb = 0, ub = 1, name = 'ecc')
i = traj.SV(value = oe_i[2], lb = math.radians(-90), ub = math.radians(90), name = 'inc')
Om = traj.SV(value = oe_i[3], lb = math.radians(-180), ub = math.radians(180), name = 'raan')
om = traj.SV(value = oe_i[4], lb = math.radians(-180), ub = math.radians(180), name = 'ap')
nu = traj.SV(value = oe_i[5], lb = math.radians(-180), ub = math.radians(180), name = 'ta')
m = traj.SV(value = m0, lb = 1000.0, ub = m0, name = 'mass')
t = traj.SV(value = 0.0, lb = 0.0, ub = t_max, name = 'time')
traj.periodic(i)
traj.periodic(Om)
traj.periodic(om)
traj.periodic(nu)
q = np.zeros(nt)
q[-1] = 1.0
final = traj.Param(value = q)
#objective function
tf = traj.FV(1.2 * ((m0 - m)/dm), lb = 0.0, ub = t_max)
#tf = traj.FV(1.0, lb = 0.0, ub = t_max)
tf.STATUS = 1
#manipulating variables and initial guesses
al_a = traj.MV(value = -1.0, lb = -2.0, ub = 2.0, name = 'al_a')
al_a.STATUS = 1
l_e = traj.MV(value = 0.001, lb = 0.0, ub = 1.0E6, name = 'l_e')
l_e.STATUS = 1
l_i = traj.MV(value = 1.0, lb = 0.0, ub = 1.0E6, name = 'l_i')
l_i.STATUS = 1
#equations
p = a * (1 - e**2)
r = p/(1 + e * traj.cos(nu))
rv = (np.array([r * traj.cos(nu), r * traj.sin(nu), 0]))
vv = (np.array([-traj.sin(nu) * traj.sqrt(mu/p), (e + traj.cos(nu)) * traj.sqrt(mu/p), 0]))
cO = traj.cos(Om)
sO = traj.sin(Om)
co = traj.cos(om)
so = traj.sin(om)
ci = traj.cos(i)
si = traj.sin(i)
R = np.array([[(cO * co - sO * so * ci), (-cO * so - sO * co * ci), (sO * si)], [(sO * co + cO * so * ci), (-sO * so + cO * co * ci), (-cO * si)], [(so * si), (co * si), ci]])
ri = np.array([R[0][0] * rv[0] + R[0][1] * rv[1] + R[0][2] * rv[2], R[1][0] * rv[0] + R[1][1] * rv[1] + R[1][2] * rv[2], R[2][0] * rv[0] + R[2][1] * rv[1] + R[2][2] * rv[2]])
vi = np.array([R[0][0] * vv[0] + R[0][1] * vv[1] + R[0][2] * vv[2], R[1][0] * vv[0] + R[1][1] * vv[1] + R[1][2] * vv[2], R[2][0] * vv[0] + R[2][1] * vv[1] + R[2][2] * vv[2]])
r = traj.sqrt(ri[0]**2 + ri[1]**2 + ri[2]**2)
v = traj.sqrt(vi[0]**2 + vi[1]**2 + vi[2]**2)
hi = np.cross(ri, vi)
h = traj.sqrt(hi[0]**2 + hi[1]**2 + hi[2]**2)
s_a = traj.Intermediate((-l_e * (r/a) * traj.sin(nu))/(traj.sqrt(4 * ((al_a * (a * v**2)/mu) + l_e * (e + traj.cos(nu)))**2 + l_e**2 * (r**2/a**2) * (traj.sin(nu))**2)))
c_a = traj.Intermediate((-2 * (((al_a * a * v**2)/mu) + l_e * (e + traj.cos(nu))))/(traj.sqrt(4 * ((al_a * a * v**2)/mu + l_e * (e + traj.cos(nu)))**2 + l_e**2 * (r**2/a**2) * (traj.sin(nu))**2)))
s_b = traj.Intermediate((-l_i * ((r * v)/h) * traj.cos(om + nu))/(traj.sqrt(l_i**2 * ((r**2 * v**2)/h**2) * (traj.cos(om + nu))**2 + ((4 * al_a**2 * a**2 * v**4)/mu**2) * c_a**2 + l_e**2 * ((2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a)**2))))
c_b = traj.Intermediate((((-al_a * 2 * a * v**2)/mu) * c_a - l_e * (2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a))/(traj.sqrt(l_i**2 * ((r**2 * v**2)/h**2) * (traj.cos(om + nu))**2 + ((4 * al_a**2 * a**2 * v**4)/mu**2) * c_a**2 + l_e**2 * ((2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a)**2))))
a_T = (2 * eta * P)/(m * g * Isp)
a_n = a_T * s_a * c_b
a_t = a_T * c_a * c_b
a_h = a_T * s_b
n = traj.sqrt(mu/a**3)
Om_J2 = ((-3 * n * R_E**2 * J2)/(2 * a**2 * (1 - e**2)**2)) * traj.cos(i)
om_J2 = ((3 * n * R_E**2 * J2)/(4 * a**2 * (1 - e**2)**2)) * (4 - 5 * (traj.sin(i))**2)
dt_dE = r/(n * a)
Tp = (2 * math.pi/traj.sqrt(mu)) * a**(3/2)
#deltas
tp = traj.if3(t - tf, 1, 0)
traj.Equation(tp * Tp * a.dt() == (a_t * (2 * v * a**2)/mu) * delta_t * dt_dE)
traj.Equation(tp * Tp * e.dt() == ((1/v) * (2 * (e + traj.cos(nu)) * a_t + (r/a) * a_n * traj.sin(nu))) * delta_t * dt_dE)
traj.Equation(tp * Tp * i.dt() == ((r/h) * a_h * traj.cos(om + nu)) * delta_t * dt_dE)
traj.Equation(tp * Tp * Om.dt() == ((r/(h * traj.sin(i))) * a_h * traj.sin(om + nu) + Om_J2) * delta_t * dt_dE)
traj.Equation(tp * Tp * om.dt() == ((1/(e * v)) * (2 * a_t * traj.sin(nu) - (2 * e + (r/a) * traj.cos(nu)) * a_n) - (r/(h * traj.sin(i))) * a_h * traj.sin(om + nu) * traj.cos(i) + om_J2) * delta_t * dt_dE)
traj.Equation(tp * nu.dt() == (traj.acos((traj.cos((1/dt_dE) * delta_t) - e)/(1 - e * traj.cos((1/dt_dE) * delta_t))) - nu))
traj.Equation(tp * Tp * m.dt() == ((-2 * eta * P)/((g * Isp)**2)) * delta_t * dt_dE)
traj.Equation(t.dt() == delta_t)
traj.Equation(a * final == oe_f[0])
traj.Equation(e * final == oe_f[1])
traj.Equation(i * final == oe_f[2])
traj.Equation(Om * final == oe_f[3])
traj.Equation(om * final == oe_f[4])
traj.Equation(nu * final == oe_f[5])
#solve
traj.Obj(tf)
traj.options.IMODE = 6
traj.options.SOLVER = 3
traj.options.MAX_ITER = 15000
traj.options.RTOL = 1e-6
traj.options.OTOL = 1e-6
#traj.open_folder()
traj.solve()
print('Optimal time: ' + str(tf.value[0]))
traj.solve()
#traj.open_folder(infeasibilities.txt)
The code has been cobbled together from various GEKKO examples and other questions on stackoverflow. I have tried making changes to variable types (Var vs SV), initial MV guesses, changing the objective function, and changing minimizing vs maximizing the objective, but I keep getting the same error.
I am hoping someone else might be able to put another pair of eyes on my code and maybe see a way to get it to converge.
The periodic constraints can be difficult to converge. Selective variable bounds for the trigonometric function arguments (e.g. -π to π), softened (objective-based) terminal constraints, and MPCC version of the conditional statement with m.if2() can help. Also, make sure to use the APOPT solver if you are using m.if3() otherwise the binary switching condition won’t necessarily be respected in the solution. IPOPT can help to find an initial solution, but then use APOPT to converge the integer constraints.
I recommend starting with no periodic constraints and simplify the model to get a solution. Then, add the constraints and better initial guess values based on the constrained solutions. Here are additional initialization strategies.