gekkoorbital-mechanics

constrain the initial and final values of a GEKKO ```Var``` to a data-based curve

I am trying to solve a low thrust optimal control problem for Earth orbits, i.e. going from one orbit to another. The formulation of the problem includes six states (r_1, r_2, r_3, v_1, v_2, v_3) and 3 controls (u_1, u_2, u_3) with a simplified point model of gravity. When I specify the full initial state and half of the final state, the solver converges and yields a good solution. When I try the full final state, the problem is over constrained.

My thought on how to remedy this is to allow the trajectory to depart the initial orbit at any point along the orbital curve and join the final orbit an any point along the final orbital curve, giving it more degrees of freedom. Is there a way to constrain the initial and final values of all 6 states to a cspline curve? This is what I have tried so far:

m = GEKKO(remote=True)
m.time = np.linspace(0, tof, t_steps)

theta_i = m.FV(lb = 0, ub = 360)
theta_f = m.FV(lb = 0, ub = 360)

rx_i = m.Param()
ry_i = m.Param()
rz_i = m.Param()

vx_i = m.Param()
vy_i = m.Param()
vz_i = m.Param()

rx_f = m.Param()
ry_f = m.Param()
rz_f = m.Param()

vx_f = m.Param()
vy_f = m.Param()
vz_f = m.Param()

m.cspline(theta_i, rx_i, initial.theta, initial.r[:, 0])
m.cspline(theta_i, ry_i, initial.theta, initial.r[:, 1])
m.cspline(theta_i, rz_i, initial.theta, initial.r[:, 2])

m.cspline(theta_i, vx_i, initial.theta, initial.v[:, 0])
m.cspline(theta_i, vy_i, initial.theta, initial.v[:, 1])
m.cspline(theta_i, vz_i, initial.theta, initial.v[:, 2])

m.cspline(theta_f, rx_f, final.theta, final.r[:, 0])
m.cspline(theta_f, ry_f, final.theta, final.r[:, 1])
m.cspline(theta_f, rz_f, final.theta, final.r[:, 2])

m.cspline(theta_f, vx_f, final.theta, final.v[:, 0])
m.cspline(theta_f, vy_f, final.theta, final.v[:, 1])
m.cspline(theta_f, vz_f, final.theta, final.v[:, 2])

r1 = m.Var(rx_i)
r2 = m.Var(ry_i)
r3 = m.Var(rz_i)

r1dot = m.Var(vx_i)
r2dot = m.Var(vy_i)
r3dot = m.Var(vz_i)

u1 = m.Var(lb = -max_u, ub = max_u)
u2 = m.Var(lb = -max_u, ub = max_u)
u3 = m.Var(lb = -max_u, ub = max_u)

m.Equation(r1.dt() == r1dot)
m.Equation(r2.dt() == r2dot)
m.Equation(r3.dt() == r3dot)

r = m.Intermediate(m.sqrt(r1**2 + r2**2 + r3**2))
v = m.Intermediate(m.sqrt(r1dot**2 + r2dot**2 + r3dot**3))

m.Equation(-mu*r1/r**3 == r1dot.dt() + u1)
m.Equation(-mu*r2/r**3 == r2dot.dt() + u2)
m.Equation(-mu*r3/r**3 == r3dot.dt() + u3)

m.fix_final(r1, rx_f)
m.fix_final(r2, ry_f)
m.fix_final(r3, rz_f)

m.fix_final(r1dot, vx_f)
m.fix_final(r2dot, vy_f)
m.fix_final(r3dot, vz_f)

m.Minimize(m.integral(u1**2 + u2**2 + u3**2))

m.options.IMODE = 6
m.options.solver = 3
#m.options.ATOL = 1e-3
m.options.MAX_ITER = 300
m.solve(disp=True)    # solve

With cspline, I am trying to allow GEKKO to pick a fixed value of the true anomaly (the parameter that determines how far along the orbit you are) from data that I have generated about the states at sampled true anomalies that would have associated position and velocity states. Any help would be much appreciated!

SOLUTION:

I implemented the end constraints as "soft constraints", i.e. quantities to be minimized, instead of hard constraints which would use fix_final. Specific implementation is as follows.

final = np.zeros(len(m.time))
final[-1] = 1
final = m.Param(value=final) 

m.Obj(final*(r1-final_o.r[0, 0])**2)
m.Obj(final*(r2-final_o.r[0, 1])**2)
m.Obj(final*(r3-final_o.r[0, 2])**2)

m.Obj(final*(r1dot-final_o.v[0, 0])**2)
m.Obj(final*(r2dot-final_o.v[0, 1])**2)
m.Obj(final*(r3dot-final_o.v[0, 2])**2)

final makes the solver only consider the last item (the end boundary) when multiplied.

Solution

  • It is generally much harder for an optimizer to exactly reach a fixed endpoint, especially when it depends on a complex sequence of moves. This often leads to infeasible solutions. An alternative is to create a soft constraint (objective minimization) to penalize deviations from the final trajectory. Here is an example that is similar:

    Inverted Pendulum

    import matplotlib.animation as animation
import numpy as np
from gekko import GEKKO

#Defining a model
m = GEKKO()

#################################
#Weight of item
m2 = 1
#################################

#Defining the time, we will go beyond the 6.2s
#to check if the objective was achieved
m.time = np.linspace(0,8,100)
end_loc = int(100.0*6.2/8.0)

#Parameters
m1a = m.Param(value=10)
m2a = m.Param(value=m2)
final = np.zeros(len(m.time))
for i in range(len(m.time)):
    if m.time[i] < 6.2:
        final[i] = 0
    else:
        final[i] = 1
final = m.Param(value=final)

#MV
ua = m.Var(value=0)

#State Variables
theta_a = m.Var(value=0)
qa = m.Var(value=0)
ya = m.Var(value=-1)
va = m.Var(value=0)

#Intermediates
epsilon = m.Intermediate(m2a/(m1a+m2a))

#Defining the State Space Model
m.Equation(ya.dt() == va)
m.Equation(va.dt() == -epsilon*theta_a + ua)
m.Equation(theta_a.dt() == qa)
m.Equation(qa.dt() == theta_a -ua)

#Define the Objectives
#Make all the state variables be zero at time >= 6.2
m.Obj(final*ya**2)
m.Obj(final*va**2)
m.Obj(final*theta_a**2)
m.Obj(final*qa**2)

m.fix(ya,pos=end_loc,val=0.0)
m.fix(va,pos=end_loc,val=0.0)
m.fix(theta_a,pos=end_loc,val=0.0)
m.fix(qa,pos=end_loc,val=0.0)
#Try to minimize change of MV over all horizon
m.Obj(0.001*ua**2)

m.options.IMODE = 6 #MPC
m.solve() #(disp=False)

    This example uses a combination of soft and hard constraints to help it find the optimal solution.