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!
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).