Space Shuttle Reentry

As a next example, we will solve another classic problem outlined by Betts in [1].

Note

We should note that this example is derived entirely from publicly available sources (see [1-4]). Furthermore, this problem is routinely used [1-4] as a benchmark for optimal control packages such as ASSET.

This problem involves maximizing the cross range of the Space Shuttle during reentry. The dynamics are written in spherical coordinates and incorporate gravity and an empirical model for the lift and drag characteristics of the shuttle. The state variables are the altitude, \(h\), latitude, \(\theta\), velocity, \(v\), flight path angle, \(\gamma\), and azimuth, \(\psi\). The controls are the angle of attack, \(\alpha\), and bank angle, \(\beta\), of the vehicle.

\[ \begin{align}\begin{aligned}\dot{h} &= v \sin(\gamma)\\\dot{\theta} &= \frac{v}{r} \cos(\gamma) \cos(\psi)\\\dot{v} &= -\frac{D}{m} - g \sin(\gamma)\\\dot{\gamma} &= \frac{L}{mv}\cos(\beta) + \cos(\gamma)\left( \frac{v}{r} - \frac{g}{v} \right)\\\dot{\psi} &= \frac{L}{mv \cos(\gamma)}\sin(\beta) +\frac{v}{r \cos(\theta)}\cos(\gamma)\sin(\psi)\sin(\theta)\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}D,L &= \frac{1}{2} C_{L,D} S \rho v^2\\g &= \frac{\mu}{r^2}\\r &= R_e + h\\\rho &= \rho_0 e^{-h/h_0}\\C_L &= a_0 + a_1 \hat{\alpha}\\C_D &= b_0 + b_1 \hat{\alpha} + b_2 \hat{\alpha}^2\end{aligned}\end{align} \]

The initial and terminal values of the state variables for this problem are given below.

\[ \begin{align}\begin{aligned}h(0) &=260000\; ft \quad \quad v(0) &= 25600 \; ft/sec \quad \quad \gamma(0)&=-1 ^\circ \quad \quad \theta(0)=0 ^\circ \quad \quad \phi(0)=0 ^\circ\\h(t_f)&=80000 \; ft \quad \quad v(t_f)&= 2500 \; ft/sec \quad \quad \gamma(t_f)&=-5 ^\circ\end{aligned}\end{align} \]

The objective is to maximize the cross range, which for this model and initial conditions is equivalent to maximizing \(\theta(t_f)\).

Similar to Betts, we examine the solution to this problem both with and without a path constraint on the wing leading edge heating rate, \(q\).

\[ \begin{align}\begin{aligned}q &= q_r q_{\alpha} < 70 \frac{BTU}{ft^2 s}\\q_r &= 17700 \sqrt{\rho}(.0001 v)^{3.07}\\q_{\alpha} &= c_0 + c_1 \hat{\alpha} + c_2 \hat{\alpha}^2 + c_3 \hat{\alpha}^3\end{aligned}\end{align} \]

As we have mentioned previously, solving problems in standard units (Miles, Km, fps etc.) is typically very ill conditioned and degrades the performance of an otherwise well posed problem. Therefore, as the first step to solving this problem, we will non-dimensionalize all variables and equations to be of order unity. This is done by defining characteristic length, mass, and time units from which we can define other derived units through dimensional analysis. In this example we select Lstar to be 100,000 feet and tstar to be 60 seconds. Mstar is then set to be the mass of the shuttle. After constructing our derived units, we can simply divide our physical constants by the appropriate unit to get their non-dimensional value.

################### Non Dimensionalize ##################################
g0 = 32.2
W  = 203000

Lstar = 100000.0     ## feet
Tstar = 60.0         ## sec
Mstar = W/g0         ## slugs

Vstar   = Lstar/Tstar
Fstar   = Mstar*Lstar/(Tstar**2)
Astar   = Lstar/(Tstar**2)
Rhostar = Mstar/(Lstar**3)
BTUstar = 778.0*Lstar*Fstar
Mustar  = (Lstar**3)/(Tstar**2)

tmax = 2500            /Tstar
Re = 20902900          /Lstar
S  = 2690.0            /(Lstar**2)
m  = (W/g0)            /Mstar
mu = (0.140765e17)     /Mustar
rho0 =.002378          /Rhostar
h_ref = 23800          /Lstar

a0 = -.20704
a1 = .029244

b0 = .07854
b1 = -.61592e-2
b2 = .621408e-3

c0 =  1.0672181
c1 = -.19213774e-1
c2 = .21286289e-3
c3 = -.10117e-5

Qlimit = 70.0

##############################################################################

Having non-dimensionalized our constants, we can now write the EOM’s as an oc.ODEBase object as we have done in previous examples. For this model, there are five state variables \((h,\theta,v,\gamma,\psi)\) and two control variables \((\alpha,\beta)\).

class ShuttleReentry(oc.ODEBase):
    def __init__(self):

        Xvars = 5
        Uvars = 2

        ############################################################
        XtU  = oc.ODEArguments(Xvars,Uvars)


        h,theta,v,gamma,psi = XtU.XVec().tolist()

        alpha,beta = XtU.UVec().tolist()


        alphadeg = (180.0/np.pi)*alpha

        CL  = a0 + a1*alphadeg
        CD  = b0 + b1*alphadeg + b2*(alphadeg**2)
        rho = rho0*vf.exp(-h/h_ref)
        r   = h + Re

        L   = 0.5*CL*S*rho*(v**2)
        D   = 0.5*CD*S*rho*(v**2)
        g   = mu/(r**2)

        sgam = vf.sin(gamma)
        cgam = vf.cos(gamma)

        sbet = vf.sin(beta)
        cbet = vf.cos(beta)

        spsi = vf.sin(psi)
        cpsi = vf.cos(psi)
        tantheta = vf.tan(theta)

        hdot     = v*sgam
        thetadot = (v/r)*cgam*cpsi
        vdot     = -D/m - g*sgam
        gammadot = (L/(m*v))*cbet +cgam*(v/r - g/v)
        psidot   = L*sbet/(m*v*cgam) + (v/(r))*cgam*spsi*tantheta


        ode = vf.stack([hdot,thetadot,vdot,gammadot,psidot])
        ##############################################################
        super().__init__(ode,Xvars,Uvars)

Additionally, we can express our heating rate constraint as an ASSET VectorFunction.

def QFunc():
    h,v,alpha = Args(3).tolist()
    alphadeg = (180.0/np.pi)*alpha
    rhodim = rho0*vf.exp(-h/h_ref)*Rhostar
    vdim = v*Vstar

    qr = 17700*vf.sqrt(rhodim)*((.0001*vdim)**3.07)
    qa = c0 + c1*alphadeg + c2*(alphadeg**2)+ c3*(alphadeg**3)

    return qa*qr

Next we must define a suitable initial guess for the optimization. Bett’s problem definition places an upper limit of 2500 sec on this problem, but we assume an initial guess of \(t_f\) = 1000 sec as is done in [2,3]. We are given initial and terminal values of the altitude, velocity, and \(\gamma\), so it is natural to construct the initial guess for these state variables as linear functions over the interval (\(0-t_f\)). \(\psi\) and \(\theta\) are only given initial values and we have no good physical intuition for how they will evolve so like [2,3] our initial guess assumes that they are constant. For both control angles, we just assume that they are 0.

tf  = 1000/Tstar

ht0  = 260000/Lstar
htf  = 80000 /Lstar
vt0  = 25600/Vstar
vtf  = 2500 /Vstar


gammat0 = np.deg2rad(-1.0)
gammatf = np.deg2rad(-5.0)
psit0   = np.deg2rad(90.0)


ts = np.linspace(0,tf,200)

TrajIG = []
for t in ts:
    X = np.zeros((8))
    X[0] = ht0*(1-t/tf) + htf*t/tf
    X[1] = 0
    X[2] = vt0*(1-t/tf) + vtf*t/tf
    X[3] = gammat0*(1-t/tf) + gammatf*t/tf
    X[4] = psit0
    X[5] = t
    X[6] =.00
    X[7] =.00
    TrajIG.append(np.copy(X))

With preliminaries completed we can now solve the problem. We first construct our ode and phase object, and use 40 LGL3 segments to discretize the problem. We then enforce our known initial conditions as a boundary value constraint. Next, we apply the given bounds on our states and controls as path constraints and also place the specified upper bound on the final time. Last, we enforce the terminal conditions on altitude, velocity, and flight path angle, and then specify that the objective is to minimize \(-\Delta \theta\). This is equivalent to maximizing \(\Delta \theta\). Given our rather poor initial guess for this problem, PSIOPT is invoked in solve_optimize mode, so that it first finds a feasible solution satisfying all constraints before minimizing the objective. Furthermore, we enable the line search as an extra safe-guard.

ode = ShuttleReentry()

phase = ode.phase("LGL3",TrajIG,40)

phase.addBoundaryValue("Front",range(0,6),TrajIG[0][0:6])
phase.addLUVarBounds("Path",[1,3],np.deg2rad(-89.0),np.deg2rad(89.0),1.0)
phase.addLUVarBound("Path",6,np.deg2rad(-90.0),np.deg2rad(90.0),1.0)
phase.addLUVarBound("Path",7,np.deg2rad(-90.0),np.deg2rad(1.0) ,1.0)
phase.addUpperDeltaTimeBound(tmax,1.0)
phase.addBoundaryValue("Back" ,[0,2,3],[htf,vtf,gammatf])
phase.addDeltaVarObjective(1,-1.0)
phase.setThreads(8,8)

## Our IG is bad, so i turn on line search
phase.optimizer.set_SoeLSMode("L1")
phase.optimizer.set_OptLSMode("L1")
phase.optimizer.set_PrintLevel(1)

## IG is bad, solve first before optimize
phase.solve_optimize()

#Refine to more segments and Reoptimize
phase.refineTrajManual(300)
phase.optimize()

Traj1 = phase.returnTraj()

## Add in Heating Rate Constraint, scale so rhs is order 1
phase.addUpperFuncBound("Path",QFunc(),[0,2,6],Qlimit,1/Qlimit)
phase.optimize()

Traj2 = phase.returnTraj()

print("Final Time:",Traj1[-1][5]*Tstar,"(s) , Final Cross Range:",Traj1[-1][1]*180/np.pi, " deg")
print("Final Time:",Traj2[-1][5]*Tstar,"(s) , Final Cross Range:",Traj2[-1][1]*180/np.pi, " deg")


Plot(Traj1,Traj2)

For this problem, PSIOPT is able to find a feasible solution in 29 iterations of the solve algorithm, and then an optimum solution after another 98 iterations in the optimize algorithm. We then refine the trajectory to a higher number of segments and re-optimize the solution, which converges in only 5 iterations. The total run-time (i9-12900k) is 90 milliseconds. The final objective value for \(\Delta \theta\) is 34.141 degrees, which is exactly that given by Betts in [1]. Next we add the path constraint on leading edge heating rate to the phase and optimize the new problem using the previous solution as the initial guess. Owing to the excellent initial guess, the heat rate limited problem converges in another 24 iterations taking only 60 milliseconds. The additional heating rate constraint reduces the maximum cross range of the shuttle to 30.63 degrees (also identical to Betts). A plot of the converged state and control histories for both problem formulations can be seen below.

The complete code for this example is listed at the bottom of this page.

References

1. Betts, J.T. “Practical methods for Optimal Control and Estimation Using Nonlinear Programming”, Cambridge University Press, 2009

2. Agamawi, Y. M., & Rao, A. V. (2020). Cgpops: A c++ software for solving multiple-phase optimal control problems using adaptive gaussian quadrature collocation and sparse nonlinear programming. ACM Transactions on Mathematical Software (TOMS), 46(3), 1-38.

3. Patterson, M. A., & Rao, A. V. (2014). GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming. ACM Transactions on Mathematical Software (TOMS), 41(1), 1-37.

4. Falck, R. (2022). https://openmdao.github.io/dymos/examples/reentry/reentry.html.

Full Code

import numpy as np
import asset_asrl as ast
import matplotlib.pyplot as plt

vf        = ast.VectorFunctions
oc        = ast.OptimalControl
Args      = vf.Arguments

'''
Space Shuttle Reentry
Betts, J.T. Practical methods for Optimal Control and Estimation Using Nonlinear Programming, Cambridge University Press, 2009
'''

################### Non Dimensionalize ##################################
g0 = 32.2
W  = 203000

Lstar = 100000.0     ## feet
Tstar = 60.0         ## sec
Mstar = W/g0         ## slugs

Vstar   = Lstar/Tstar
Fstar   = Mstar*Lstar/(Tstar**2)
Astar   = Lstar/(Tstar**2)
Rhostar = Mstar/(Lstar**3)
BTUstar = 778.0*Lstar*Fstar
Mustar  = (Lstar**3)/(Tstar**2)

tmax = 2500            /Tstar
Re = 20902900          /Lstar
S  = 2690.0            /(Lstar**2)
m  = (W/g0)            /Mstar
mu = (0.140765e17)     /Mustar
rho0 =.002378          /Rhostar
h_ref = 23800          /Lstar

a0 = -.20704
a1 = .029244

b0 = .07854
b1 = -.61592e-2
b2 = .621408e-3

c0 =  1.0672181
c1 = -.19213774e-1
c2 = .21286289e-3
c3 = -.10117e-5

Qlimit = 70.0

##############################################################################
class ShuttleReentry(oc.ODEBase):
    def __init__(self):

        Xvars = 5
        Uvars = 2

        ############################################################
        XtU  = oc.ODEArguments(Xvars,Uvars)


        h,theta,v,gamma,psi = XtU.XVec().tolist()

        alpha,beta = XtU.UVec().tolist()


        alphadeg = (180.0/np.pi)*alpha

        CL  = a0 + a1*alphadeg
        CD  = b0 + b1*alphadeg + b2*(alphadeg**2)
        rho = rho0*vf.exp(-h/h_ref)
        r   = h + Re

        L   = 0.5*CL*S*rho*(v**2)
        D   = 0.5*CD*S*rho*(v**2)
        g   = mu/(r**2)

        sgam = vf.sin(gamma)
        cgam = vf.cos(gamma)

        sbet = vf.sin(beta)
        cbet = vf.cos(beta)

        spsi = vf.sin(psi)
        cpsi = vf.cos(psi)
        tantheta = vf.tan(theta)

        hdot     = v*sgam
        thetadot = (v/r)*cgam*cpsi
        vdot     = -D/m - g*sgam
        gammadot = (L/(m*v))*cbet +cgam*(v/r - g/v)
        psidot   = L*sbet/(m*v*cgam) + (v/(r))*cgam*spsi*tantheta


        ode = vf.stack([hdot,thetadot,vdot,gammadot,psidot])
        ##############################################################
        super().__init__(ode,Xvars,Uvars)

def QFunc():
    h,v,alpha = Args(3).tolist()
    alphadeg = (180.0/np.pi)*alpha
    rhodim = rho0*vf.exp(-h/h_ref)*Rhostar
    vdim = v*Vstar

    qr = 17700*vf.sqrt(rhodim)*((.0001*vdim)**3.07)
    qa = c0 + c1*alphadeg + c2*(alphadeg**2)+ c3*(alphadeg**3)

    return qa*qr

#############################################################################

def Plot(Traj1,Traj2):
    TT1 = np.array(Traj1).T
    TT2 = np.array(Traj2).T

    fig, axs = plt.subplots(4,1)

    axs[0].plot(TT1[5]*Tstar/60.0,TT1[0]*Lstar/5280,label='No Q limit',color='r')
    axs[0].plot(TT2[5]*Tstar/60.0,TT2[0]*Lstar/5280,label='Q limited',color='b')
    axs[0].set_ylabel("Altitude (Miles)")


    axs[1].plot(TT1[5]*Tstar/60.0,TT1[2]*Vstar,label='No Q limit',color='r')
    axs[1].plot(TT2[5]*Tstar/60.0,TT2[2]*Vstar,label='Q limited',color='b')

    axs[1].set_ylabel(r"Velocity $\frac{ft}{s}$")


    axs[2].plot(TT1[5]*Tstar/60.0,np.rad2deg(TT1[6]),label='Angle of Attack No Q limit',color='r')
    axs[2].plot(TT1[5]*Tstar/60.0,np.rad2deg(TT1[7]),label='Bank Angle No Q limit',color='r',linestyle='dotted')
    axs[2].plot(TT2[5]*Tstar/60.0,np.rad2deg(TT2[6]),label='Angle of Attack  Q limited',color='b')
    axs[2].plot(TT2[5]*Tstar/60.0,np.rad2deg(TT2[7]),label='Bank Angle  Q limited',color='b',linestyle='dotted')

    axs[2].set_ylabel("Angle(deg)")

    qfunc = QFunc().eval(8,[0,2,6])

    qs1 = [qfunc(T)[0] for T in Traj1]
    qs2 = [qfunc(T)[0] for T in Traj2]

    axs[3].plot(TT1[5]*Tstar/60.0,qs1,label='No Q limit',color='r')
    axs[3].plot(TT2[5]*Tstar/60.0,qs2,label='Q limited',color='b')
    axs[3].set_ylabel(r"Q $(\frac{BTU}{ft^2 *s})$")

    axs[3].plot(TT2[5]*Tstar/60.0,np.ones_like(TT2[5])*70,label=r'Q = 70 $\frac{BTU}{ft^2 *s}$',color='k',linestyle='dashed')


    for i in range(0,4):
        axs[i].grid(True)
        axs[i].set_xlabel("Time (min)")
        axs[i].legend()

    fig.set_size_inches(8.0, 11.0, forward=True)
    fig.tight_layout()

    plt.show()




if __name__ == "__main__":
    ##########################################################################
    tf  = 1000/Tstar

    ht0  = 260000/Lstar
    htf  = 80000 /Lstar
    vt0  = 25600/Vstar
    vtf  = 2500 /Vstar


    gammat0 = np.deg2rad(-1.0)
    gammatf = np.deg2rad(-5.0)
    psit0   = np.deg2rad(90.0)


    ts = np.linspace(0,tf,200)

    TrajIG = []
    for t in ts:
        X = np.zeros((8))
        X[0] = ht0*(1-t/tf) + htf*t/tf
        X[1] = 0
        X[2] = vt0*(1-t/tf) + vtf*t/tf
        X[3] = gammat0*(1-t/tf) + gammatf*t/tf
        X[4] = psit0
        X[5] = t
        X[6] =.00
        X[7] =.00
        TrajIG.append(np.copy(X))


    ################################################################

    ode = ShuttleReentry()

    phase = ode.phase("LGL3",TrajIG,40)

    phase.addBoundaryValue("Front",range(0,6),TrajIG[0][0:6])
    phase.addLUVarBounds("Path",[1,3],np.deg2rad(-89.0),np.deg2rad(89.0),1.0)
    phase.addLUVarBound("Path",6,np.deg2rad(-90.0),np.deg2rad(90.0),1.0)
    phase.addLUVarBound("Path",7,np.deg2rad(-90.0),np.deg2rad(1.0) ,1.0)
    phase.addUpperDeltaTimeBound(tmax,1.0)
    phase.addBoundaryValue("Back" ,[0,2,3],[htf,vtf,gammatf])
    phase.addDeltaVarObjective(1,-1.0)
    phase.setThreads(8,8)

    ## Our IG is bad, so i turn on line search
    phase.optimizer.set_SoeLSMode("L1")
    phase.optimizer.set_OptLSMode("L1")
    phase.optimizer.set_PrintLevel(1)

    ## IG is bad, solve first before optimize
    phase.solve_optimize()

    #Refine to more segments and Reoptimize
    phase.refineTrajManual(300)
    phase.optimize()

    Traj1 = phase.returnTraj()

    ## Add in Heating Rate Constraint, scale so rhs is order 1
    phase.addUpperFuncBound("Path",QFunc(),[0,2,6],Qlimit,1/Qlimit)
    phase.optimize()

    Traj2 = phase.returnTraj()

    print("Final Time:",Traj1[-1][5]*Tstar,"(s) , Final Cross Range:",Traj1[-1][1]*180/np.pi, " deg")
    print("Final Time:",Traj2[-1][5]*Tstar,"(s) , Final Cross Range:",Traj2[-1][1]*180/np.pi, " deg")


    Plot(Traj1,Traj2)