Question: I am trying to propagate a spacecraft to optimize the time of flight using IPOPT in GEKKO/Python. Here is the code for my GEKKO model: m = GEKKO() #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, name = 'sma') e = m.Var(value = oe_i[1], lb = 0, ub = 1, name = 'ecc') i = m.Var(value = oe_i[2], lb = 0, ub = math.radians(90), name = 'inc') Om = m.Var(value = oe_i[3], lb = 0, ub = math.radians(360), name = 'raan') om = m.Var(value = oe_i[4], lb = 0, ub = math.radians(360), name = 'ap') nu = m.Var(value = oe_i[5], lb = 0, ub = math.radians(360), name = 'ta') mass = m.Var(value = m0, lb = 0, ub = m0, name = 'mass') #objective function tf = m.FV(1.2 * ((m0 - mass)/dm), lb = 0, ub = t_max) tf.STATUS = 1 #propagate t = 0 while t <= tf: deltas, Tp = propagate(a, e, i, Om, om, nu, mass) 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])) t = t + delta_t #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.time = np.linspace(0,0.2,100) m.Obj(tf) 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.open_folder() m.solve(display=false) # Solve print('Optimal time: ' + str(tf.value[0])) m.solve() m.open_folder(infeasibilities.txt)  I know that the problem I am having is related to the propagate part. I want to propagate from the orbit corresponding to oe_i in 3600 s increments (delta_t) from time 0 to the final time (which is my objective function), achieving the orbit corresponding to oe_f, using the propagate function, which relies on the variations in the manipulation variables. I had originally tried propagating without any sort of loop to go from 0 to the end time, and the model ran fine, but never found a solution. Looking at that code, I realized that it was not actually propagating the orbit over a long period of time, which is why I added the loop. I tried a for loop first, but had similar problems with errors about tf not being an int. I tried looping time to reach the calculated value for tf (1.2 * ((m0 - mass)/dm)), but had the problem with mass not being able to be used in the calculations. If anyone is able to point out where I am going wrong in trying to do my propagation or to an example similar to what I am attempting, I would appreciate it. Thanks!

Answer: Conditional statements can't be used to define the equations of the model such as: while t <= tf: deltas, Tp = propagate(a, e, i, Om, om, nu, mass) 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])) t = t + delta_t  because the value of tf is determined by the optimizer. In Gekko, there is a model building phase where the variables and equations are defined. After the model building phase, there is a model compilation phase that converts the model to byte-code and it gives the model to the solver. There are no callbacks to the Python code. Use m.if3() to create a switching variables such as: p = m.if3(t-tf,1,0) deltas, Tp = propagate(a, e, i, Om, om, nu, mass) m.Equation(p * Tp * a.dt() == p*(deltas[0] * delta_t * deltas[7])) m.Equation(p * Tp * e.dt() == p*(deltas[1] * delta_t * deltas[7])) m.Equation(p * Tp * i.dt() == p*(deltas[2] * delta_t * deltas[7])) m.Equation(p * Tp * Om.dt() == p*(deltas[3] * delta_t * deltas[7])) m.Equation(p * Tp * om.dt() == p*(deltas[4] * delta_t * deltas[7])) m.Equation(p * nu.dt() == p*deltas[5] * delta_t) m.Equation(p * Tp * mass.dt() == p*(deltas[6] * delta_t * deltas[7]))  This turns on the equation with p=1 when t<tf and the equation if off when p=0 and t>tf. This way, the optimizer can choose the value of p and the compiled equations automatically use this input. There is an example in the paper:  Beal, L., Park, J., Petersen, D., Warnick, S., Hedengren, J.D., Combined Model Predictive Control and Scheduling with Dominant Time Constant Compensation, 2017, doi: 10.1016/j.compchemeng.2017.04.024.  that is related to this problem of setting arrival constraints mid-way through the time horizon (see Section 4).