Question: In my code, I get the following error when running: Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) ((((((((((-cos(v4)))*(sin(v5)))-((((sin(v4))*(cos(v5))))*(cos(v3)))))*(((sqrt(( 398574405096000.0/((v1)*((1-((v2)^(2))))))))*([(-sin(v6))(v2+cos(v6))0])))))^(2 ))+((((((((-sin(v4)))*(sin(v5)))+((((cos(v4))*(cos(v5))))*(cos(v3)))))*(((sqrt( (398574405096000.0/((v1)*((1-((v2)^(2))))))))*([(-sin(v6))(v2+cos(v6))0])))))^( 2)))+((((((cos(v5))*(sin(v3))))*(((sqrt((398574405096000.0/((v1)*((1-((v2)^(2)) ))))))*([(-sin(v6))(v2+cos(v6))0])))))^(2))) STOPPING...  I tried searching my code for variations of the math functions (sqrt, cos, etc.) to see if I could find something that looked like the above equation, but I cannot find it in any way. I assume GEKKO manipulated some things around to get it, likely as part of the solver. My thinking is that the 'v' values are equivalent to my orbital elements, and I see that the value of mu is expressed. I'm hoping someone can put another set of eyes on my code and maybe help me out. Here is my code: from gecko import GEKKO import numpy as np import matplotlib.pyplot as plt import math def oe2rv(oe): a = oe[0] e = oe[1] i = oe[2] Om = oe[3] om = oe[4] nu = oe[5] p = a * (1 - e**2) r = p/(1 + e * m.cos(nu)) rv = np.array([r * m.cos(nu), r * m.sin(nu), 0]) vv = m.sqrt(mu/p) * np.array([-m.sin(nu), e + m.cos(nu), 0]) cO = m.cos(Om) sO = m.sin(Om) co = m.cos(om) so = m.sin(om) ci = m.cos(i) si = m.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 = R * rv vi = R * vv return ri, vi def TwoBody(ri, vi): ri_dot[0] = vi[0] ri_dot[1] = vi[1] ri_dot[2] = vi[2] r_mag = m.sqrt(ri[0]**2 + ri[1]**2 + ri[2]**2) r3 = r_mag**3 c = -mu/r3 vi_dot[0] = c * ri[0] vi_dot[1] = c * ri[1] vi_dot[2] = c * ri[2] return ri_dot, vi_dot def Euler(ri, vi): ri_dot, vi_dot = TwoBody(ri, vi) for i in range(0, 3): ri_new[i] = ri[i] + ri_dot[i] * dt vi_new[i] = vi[i] + vi_dot[i] * dt return ri_new, vi_new def trap(E_en, E_ex, e): dE = (E_en - E_ex)/20 E = np.linspace(E_en, E_ex, 20) nu_new = m.acos((m.cos(E[0]) - e)/(1 - e * m.cos(E[0]))) for i in range(1, 19): nu_new = nu_new + 2 * m.acos((m.cos(E[i]) - e)/(1 - e * m.cos(E[i]))) nu_new = m.acos((m.cos(E[19]) - e)/(1 - e * m.cos(E[19]))) nu_new = (dE/2) * nu_new return nu_new def propagate(a, e, i, Om, om, nu, mass): oe = np.array([a, e, i, Om, om, nu]) ri, vi = oe2rv(oe) r = m.sqrt(ri[0]**2 + ri[1]**2 + ri[2]**2) v = m.sqrt(vi[0]**2 + vi[1]**2 + vi[2]**2) h = np.cross(ri, vi) d1 = m.sqrt(4 * ((al_a * a * v**2)/mu + l_e * (e + m.cos(nu)))**2 + l_e**2 * (r**2/a**2) * m.sin(nu)**2) s_a = (-l_e * (r/a) * m.sin(nu))/d1 c_a = (-2 * (((al_a * a * v**2)/mu) + l_e * (e + m.cos(nu))))/d1 d2 = m.sqrt(l_i**2 * ((r**2 * v**2)/h**2) * (m.cos(om + nu))**2 + (4 * al_a**2 * a**2 * v**4 * c_a**2)/mu**2 + l_e**2 * ((2 * (e + m.cos(nu)) * c_a + r/a * m.sin(nu) * s_a)**2)) s_b = (-l_i * ((r * v)/h) * m.cos(om + nu))/d2 c_b = (((-al_a * 2 * a * v**2)/mu) * c_a - l_e * (2 * (e + m.cos(nu)) * c_a + (r/a) * m.sin(nu) * s_a))/d2 a_n = aT * s_a * c_b a_t = aT * c_a * c_b a_h = aT * s_b n = m.sqrt(mu/a**3) Om_J2 = ((-3 * n * R_E**2 * J2)/(2 * a**2 * (1 - e**2)**2)) * m.cos(i) om_J2 = ((3 * n * R_E**2 * J2)/(4 * a**2 * (1 - e**2)**2)) * (4 - 5 * (m.sin(i))**2) nu_new = trap(E_en, E_ex, e) da_dt = a_t * (2 * a**2 * v)/mu de_dt = (1/v) * (2 * (e + m.cos(nu)) * a_t + (r/a) * a_n * m.sin(nu)) di_dt = (r/h) * a_h * m.cos(om + nu) dOm_dt = (r/(h * m.sin(i))) * a_h * m.sin(om + nu) + Om_J2 dom_dt = (1/(e * v)) * (2 * a_t * m.sin(nu) - (2 * e + (r/a) * m.cos(nu)) * a_n) - (r/(h * m.sin(i))) * a_h * m.sin(om + nu) * m.cos(i) + om_J2 dnu_dt = nu_new - nu dm_dt = (-2 * eta * P)/pow((g * Isp), 2) dt_dE = r/(n * a) Tp = (2 * math.pi/m.sqrt(mu)) * a**(3/2) deltas = np.array([da_dt, de_dt, di_dt, dOm_dt, dom_dt, dnu_dt, dm_dt, dt_dE]) return deltas, Tp #initialize model m = GEKKO() #optional solver settings with APOPT Nsim = 100 #number of steps with constant thrust m.time = np.linspace(0, 0.2, Nsim) #constants mu = 3.98574405096E14 g = 9.81 R_E = 6.2781E6 J2 = 1.08262668E-3 P = 10E3 eta = 0.65 Isp = 3300 m0 = 1200 aT = (2 * eta * P)/(m0 * g * Isp) delta_t = 3600 t_max = 86400 * 200 E_en = math.pi E_ex = -math.pi oe_i = np.array([6927000, 0, math.radians(28.5), 0, 0, 0]) oe_f = np.array([42164000, 0, 0, 0, 0, 0]) v_i = m.sqrt(mu/oe_i[0]) v_f = m.sqrt(mu/oe_f[0]) dv = abs(v_i - v_f) dm = (2 * eta * P)/pow((g * Isp), 2) m_f = m0 * m.exp(-dv/(g * Isp)) #manipulating variables and initial guesses al_a = m.MV(value = -1, lb = -2, ub = 2) al_a.STATUS = 1 l_e = m.MV(value = 0.001, lb = 0, ub = 10**6) l_e.STATUS = 1 l_i = m.MV(value = 1, lb = 0, ub = 10**6) l_i.STATUS = 1 #variables and initial guesses a = m.Var(value = oe_i[0], lb = oe_i[0] - 6378000, ub = oe_f[0] + 6378000) e = m.Var(value = oe_i[1], lb = 0, ub = 1) i = m.Var(value = oe_i[2], lb = 0, ub = math.radians(90)) Om = m.Var(value = oe_i[3], lb = 0, ub = math.radians(360)) om = m.Var(value = oe_i[4], lb = 0, ub = math.radians(360)) nu = m.Var(value = oe_i[5], lb = 0, ub = math.radians(360)) mass = m.Var(value = m0, lb = 0, ub = m0) #objective function tf = m.FV(value = 1.2 * ((m0 - m_f)/dm), lb = 0, ub = t_max) tf.STATUS = 1 #propagation deltas, Tp = propagate(a, e, i, Om, om, nu, mass) m.Equation(a.dt() == (deltas[0] * delta_t * deltas[7])/Tp) m.Equation(e.dt() == (deltas[1] * delta_t * deltas[7])/Tp) m.Equation(i.dt() == (deltas[2] * delta_t * deltas[7])/Tp) m.Equation(Om.dt() == (deltas[3] * delta_t * deltas[7])/Tp) m.Equation(om.dt() == (deltas[4] * delta_t * deltas[7])/Tp) m.Equation(nu.dt() == deltas[5] * delta_t) m.Equation(mass.dt() == (deltas[6] * delta_t * deltas[7])/Tp) #starting constraints m.fix(a, pos = 0, val = oe_i[0]) m.fix(e, pos = 0, val = oe_i[1]) m.fix(i, pos = 0, val = oe_i[2]) m.fix(Om, pos = 0, val = oe_i[3]) m.fix(om, pos = 0, val = oe_i[4]) m.fix(nu, pos = 0, val = oe_i[5]) m.fix(mass, pos = 0, val = m0) #boundary constraints m.fix(a, pos = len(m.time) - 1, val = oe_f[0]) m.fix(e, pos = len(m.time) - 1, val = oe_f[1]) m.fix(i, pos = len(m.time) - 1, val = oe_f[2]) m.fix(Om, pos = len(m.time) - 1, val = oe_f[3]) m.fix(om, pos = len(m.time) - 1, val = oe_f[4]) m.fix(nu, pos = len(m.time) - 1, val = oe_f[5]) m.fix(mass, pos = len(m.time) - 1, val = 0) m.Obj(tf) #minimize final time m.options.IMODE = 6 # non-linear model m.options.SOLVER = 3 # solver (IPOPT) m.options.MAX_ITER = 15000 m.options.RTOL = 1e-7 m.options.OTOL = 1e-7 m.solve(disp=True, debug=True) # Solve print('Optimal time: ' + str(tf.value[0])) m.solve(disp=True) m.open_folder(infeasibilities.txt)  After doing some playing around, I believe the issue is that I am using the manipulating variables ('al_a', 'l_e' and 'l_i') in the 'propagate' function. Does that make sense as a possible problem? If that is the problem, is it possible to use the values of those variables in that function - and, if so, how?

Answer: The square brackets indicate that a list or numpy array was used instead of a scalar value in one of the expressions. Adding names (e.g. name='a') to the variables helps with a more readable model apm file that is in the local run directory m.path. Open the directory with m.open_folder(). #variables and initial guesses a = m.Var(value = oe_i[0], lb = oe_i[0] - 6378000, ub = oe_f[0] + 6378000,name='a') e = m.Var(value = oe_i[1], lb = 0, ub = 1,name='e') i = m.Var(value = oe_i[2], lb = 0, ub = math.radians(90),name='i') Om = m.Var(value = oe_i[3], lb = 0, ub = math.radians(360),name='om1') om = m.Var(value = oe_i[4], lb = 0, ub = math.radians(360),name='om2') nu = m.Var(value = oe_i[5], lb = 0, ub = math.radians(360),name='nu') mass = m.Var(value = m0, lb = 0, ub = m0,name='mass')  This gives an updated version of the error with: @error: Equation Definition Equation without an equality (=) or inequality (>,<) ((((((((((-cos(om1)))*(sin(om2)))-((((sin(om1))*(cos(om2))))*(cos(i)))))*(((sqr t((398574405096000.0/((a)*((1-((e)^(2))))))))*([(-sin(nu))(e+cos(nu))0])))))^(2 ))+((((((((-sin(om1)))*(sin(om2)))+((((cos(om1))*(cos(om2))))*(cos(i)))))*(((sq rt((398574405096000.0/((a)*((1-((e)^(2))))))))*([(-sin(nu))(e+cos(nu))0])))))^( 2)))+((((((cos(om2))*(sin(i))))*(((sqrt((398574405096000.0/((a)*((1-((e)^(2)))) ))))*([(-sin(nu))(e+cos(nu))0])))))^(2))) STOPPING...  One of the errors is in the propagate function where h is a (3,3) numpy array but is used in equations. h = np.cross(ri, vi) s_b = (-l_i * ((r * v)/h) * m.cos(om + nu))/d2  After the equations are corrected, one additional suggestion is to move Tp to the right-hand side of the equation to avoid potential divide-by-zero. m.Equation(Tp*a.dt() == (deltas[0] * delta_t * deltas[7])) m.Equation(Tp*e.dt() == (deltas[1] * delta_t * deltas[7])) m.Equation(Tp*i.dt() == (deltas[2] * delta_t * deltas[7])) m.Equation(Tp*Om.dt() == (deltas[3] * delta_t * deltas[7])) m.Equation(Tp*om.dt() == (deltas[4] * delta_t * deltas[7])) m.Equation(nu.dt() == deltas[5] * delta_t) m.Equation(Tp*mass.dt() == (deltas[6] * delta_t * deltas[7]))