Solving Maximum Radius Orbit Transfer OCP in GEKKO

user_19392856 0
Viewed 45

I've just started looking at GEKKO for solving OCPs, and I've run into an issue. How would one go about solving the maximum radius orbit transfer problem as seen in Bryson and Ho (also here). The main question is adding the terminal constraint on the velocity v(tf), which is a function of the final position r(tf).

Any help is greatly appreciated!

Thank you!

import numpy as np
import matplotlib.pylab as plt
from gekko import GEKKO


if __name__ == '__main__':

    m = GEKKO()

    # Constants
    nt = 101
    m.time = np.linspace(0, 193, nt)
    thr = 0.85
    mu = (3.32/193)**2
    m0 = 10000
    mdot = 12.9

    # Params
    u1 = m.MV(value=0)
    u1.STATUS = 1 
    u1.DCOST = 0

    u2 = m.MV(value=np.pi, fixed_initial=False)
    u2.STATUS = 1 
    u2.DCOST = 0

    # Vars
    t = m.Var(value=0, fixed_initial=True)
    r = m.Var(value=1, fixed_initial=True)
    u = m.Var(value=0, fixed_initial=True)
    v = m.Var(value=np.sqrt(mu), fixed_initial=True)

    # Find Index at Final Time
    p = np.zeros(nt)
    p[-1] = 1.0
    final = m.Param(value=p)

    # Equations
    m.Equation(t.dt() == 1)
    m.Equation(r.dt() == u)
    m.Equation(u.dt() == v**2/r - mu/r**2 + thr*u1/(m0-mdot*t))
    m.Equation(v.dt() == -u*v/r + thr*u2/(m0-mdot*t))

    # Constraints
    m.fix(u, pos=nt-1, val=0.0)
    #m.fix(v, pos=nt-1, val=m.sqrt(mu/r[nt-1])) # <-- How do I correctly add this constraint?

    # Objective Function -- Max Radius at Final Time
    m.Obj(-r*final)

    # Solve!
    m.options.IMODE = 6
    m.options.NODES = 4
    m.options.MV_TYPE = 1
    m.options.SOLVER = 3
    m.solve()
gekko python
1 Answers

There were a couple differences with the original link so the problem is reverted to the original equations. It is easiest for the solver when using soft terminal constraints, but hard terminal constraints (with m.fix()) can also work.

m.Minimize(100*final*(u**2))
m.Minimize(100*final*(v-m.sqrt(mmu/r))**2)

solution

The solution is similar but there are some differences. A small DCOST value of 1e-4 avoids excessive controls movement.

import numpy as np
import matplotlib.pylab as plt
from gekko import GEKKO
m = GEKKO()
nt = 101
m.time = np.linspace(0,1,nt)
r0 = 1; mmu = 11; th = 1.55
m0 = 1; rm0 = -0.25

# Params
u1 = m.MV(value=0, lb=0)
u1.STATUS = 1 
u1.DCOST = 1e-4

# Vars
t = m.Param(value=m.time)
r = m.Var(value=r0)
u = m.Var(value=0)
v = m.Var(value=np.sqrt(mmu/r0))

# Find Index at Final Time
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(r.dt() == u)
m.Equation(r**2 * u.dt() == r*v**2 - mmu + r**2 * th*m.sin(u1)/(m0+rm0*t))
m.Equation(r*v.dt() == -u*v + r*th*m.cos(u1)/(m0+rm0*t))

# Soften Constraints
m.Minimize(100*final*(u**2))
m.Minimize(100*final*(v-m.sqrt(mmu/r))**2)

# Objective Function -- Max Radius at Final Time
m.Maximize(r*final)

# Solve!
m.options.IMODE = 6
m.options.NODES = 3
m.options.MV_TYPE = 1
m.options.SOLVER = 3
m.solve()

# Plot results
plt.figure()
plt.subplot(2,1,1)
plt.title('Spacecraft states')
plt.plot(m.time,r,'b-',label='r')
plt.plot(m.time,u,'g--',label='u')
plt.plot(m.time,v,'r:',label='v')
plt.legend()

plt.subplot(2,1,2)
plt.title('Spacecraft controls')
plt.plot(m.time,u1,'b-',label='u1')
plt.xlabel('Time')
plt.legend()
plt.show()
user_23669410
Related