Does anyone see where the error in the following GEKKO-IPOPT nonlinear optimization problem is?

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?

I'm wondering if [(-sin(v6))(v2+cos(v6))0] is the problem. It almost looks like an array element, and is repeated three times in. But I'm not sure how GEKKO is manipulating things to get that. Or what it is manipulating. Any thoughts? — pbhuter, Jun 14 '23 at 15:34
Yes, that is a problem - the variables used in equations can't be lists. I'll create a more complete answer below. — John Hedengren, Jun 30 '23 at 21:12

score 2 · Answer 1 · answered Jun 19 '23 at 13:56

2

The error looks similar to what I am dealing with gekko line break in equation. Looking at your equations in the model.apm there are line breaks interrupting the mathematical relationship of each equation. From just looking at one line the solver then can't find an equality (=) or inequality (>,<). For long equations gekko can split them into smaller parts using placeholder variables gekko shorten equation. I wonder whether there might be an issue in that procedure.

answered Jun 19 '23 at 13:56

The_One_And_Only

164
5

1

I looked into Intermediate variables in GEKKO, and ended up getting that same error, only listing the intermediate variables as having the issue. As I kept playing around, I noted that it might be something to do with my manipulation variables being used in one of the functions. Any thoughts? – pbhuter Jun 24 '23 at 00:26
Good suggestion `The_One_And_Only`. The maximum equation length is 15,000 characters and it typically displays that error, if it runs into that problem. It is good to check though. – John Hedengren Jun 30 '23 at 21:10
The use of the intermediate variables helped me solve the problem of hitting the character limit when applying what I got using the solution below. Thanks! – pbhuter Jul 02 '23 at 15:41

score 0 · Accepted Answer · answered Jun 30 '23 at 21:35

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]))

Does anyone see where the error in the following GEKKO-IPOPT nonlinear optimization problem is?

2 Answers2

Linked