Lunar Ascent Trajectory Optimization

Program lascent_ocs
Lunar Ascent Trajectory Optimization
This document is the user’s manual for a Fortran computer program called lascent_ocs that uses the
Sparse Optimization Suite distributed by Applied Mathematical Analysis to solve the single-stage,
finite-burn lunar ascent trajectory optimization problem. The trajectory is modeled with a user-defined
vertical rise phase, pitch-over phase and initial and final boundary conditions for an ascent-to-orbit
phase. This computer program attempts to maximize the final mass placed into lunar orbit or minimize
the time to ascend to the final orbit. The lunar mission orbit can be circular or elliptical.
The important features of this scientific simulation are as follows:
 user-defined vertical rise and pitch-over phases

 angle-of-attack, bank angle and throttle control variables
 3-DOF flight path equations of motion relative to a spherical, rotating Moon
 user-defined vehicle mass and propulsion properties
 user-defined mission orbit semimajor axis, eccentricity and inclination
The Sparse Optimization Suite is a direct transcription method that can be used to solve a variety of
trajectory optimization problems using the following combination of numerical methods:
 collocation and implicit integration

 adaptive mesh refinement
 sparse nonlinear programming
Additional information about the mathematical techniques and numerical methods used in the Sparse
Optimization Suite can be found in the book, “Practical Methods for Optimal Control and Estimation
Using Nonlinear Programming” by John. T. Betts, SIAM, 2010.
The lascent_ocs software consists of Fortran routines that perform the following tasks:
 set algorithm control parameters and call the transcription/optimal control subroutine
 define the problem structure and perform initialization related to scaling, lower and upper
bounds, initial conditions, etc.
 compute the right-hand-side differential equations
 evaluate any point and path constraints
 display the optimal solution results and create an output file
The Sparse Optimization Suite will use this information to automatically transcribe the user’s problem
and perform the optimization using a sparse nonlinear programming method. The software allows the
user to select the type of initial guess, collocation method and other important algorithm control
parameters.
page 1
Program execution
An input file created by the user can be run from the command line or a simple batch file with a
statement similar to the following:
lascent_ocs lascent1.in
If the software is executed without an input file on the command line, the computer program will display
the following information screen and file name prompt:
*************************************
* program lascent_ocs *
* *
* lunar ascent trajectory *
* optimization *
* *
* March 22, 2011 *
*************************************
please input the name of the simulation definition file
The user should respond to this prompt with the name of a compatible input data file including the
filename extension.
Input file format and contents
The lascent_ocs software is “data-driven” by a user-created text file. The following is a typical input
file used by this computer program. In the following discussion the actual input file contents are in
courier font and all explanations are in times italic font.
This data file defines a typical maximum payload simulation with throttling enabled and a circular final
mission orbit with an inclination constraint. For this simulation, the initial azimuth is bounded so that
the software can determine the correct azimuth after the pitch-over maneuver.
Each data item within an input file is preceded by one or more lines of annotation text. Do not delete
any of these annotation lines or increase or decrease the number of lines reserved for each comment.
However, you may change them to reflect your own explanation. The annotation line also includes the
correct units and when appropriate, the valid range of the input. ASCII text input is not case sensitive
but must be spelled correctly.
The first six lines of any input file are reserved for user comments. These lines are ignored by the
software. However the input file must begin with six and only six initial text lines.
**********************************
** lunar launch vehicle simulation
** lv1.in
** November 15, 2005
** maximize payload to orbit
**********************************
The first input is an integer that defines the type of simulation. Option 2 is equivalent to a full throttle or
constant thrust simulation.
page 2
type of trajectory optimization
*******************************
1 = maximize final mass
2 = minimize flight time
------------------------
1
The next three inputs allow the user to specify the duration of the vertical rise and pitch-over maneuvers
along with the total pitch angle.
vertical rise time (seconds)
10
pitch over time duration (seconds)

30
pitch angle change during pitch over (degrees)

1.25
These next two inputs define the propulsion characteristics for the simulation. Here the user can specify
the maximum thrust and specific impulse.
**************************
propulsion characteristics
**************************
maximum thrust (newtons)

5000.0
specific impulse (seconds)

350.0
The initial vehicle mass and a guess for the final mass are defined by the next two data inputs.
****************************
vehicle mass characteristics
****************************
initial vehicle mass (kilograms)

3000.0
initial guess for final vehicle mass (kilograms)

1500.0
The following series of data items are reserved for user-defined initial conditions. To fix one or more
initial conditions, the user should input identical lower and upper bounds. To free or un-constrain one
or more initial states, set the lower and/or upper bounds to 1.0d99. Please note the units and valid data
range for each item.
************************************
initial flight conditions and bounds
************************************
NOTE 1: set upper and lower bounds to
the initial value to constrain or
"fix" the flight azimuth.
NOTE 2: set bound to 1.0d99 to ignore
-------------------------------------
initial flight azimuth (0 <= azimuth <= 360; degrees)

90.0
upper bound for initial flight azimuth (0 <= azimuth <= 360; degrees)
page 3
360.0
lower bound for initial flight azimuth (0 <= azimuth <= 360; degrees)
0.0
initial selenocentric latitude (-90 <= latitude <= + 90; degrees)

20.0
initial selenographic east longitude (0 <= longitude <= 360; degrees)

40.0
The following series of data items allow the user to define the final flight conditions. To fix one or more
conditions, the user should input identical lower and upper bounds. To free or un-constrain one or
more final states set the lower and/or upper bounds to 1.0d99. Please note the units and valid data
range for each item.
****************************************************
initial guess for final flight conditions and bounds
****************************************************
NOTE 1: set upper and lower bounds
to the final value to constrain or
"fix" a flight condition.
-------------------------------------
final time (seconds)

800.0
upper bound for final time (seconds)

1.0d99
lower bound for final time (seconds)

1.0d99
final altitude (meters)

100000.0
upper bound for final altitude (meters)

120000.0
lower bound for final altitude (meters)

10000.0
final velocity (meters/second)

1600.0
upper bound for final velocity (meters/second)

2000.0
lower bound for final velocity (meters/second)

1000.0
final flight path angle (-90 <= fpa <= +90; degrees)
0.0
upper bound for final flight path angle (-90 <= fpa <= +90; degrees)
+90.0
lower bound for final flight path angle (-90 <= fpa <= +90; degrees)
-90.0
page 4
final flight azimuth (0 <= azimuth <= 360; degrees)
90.0
upper bound for final flight azimuth (0 <= azimuth <= 360; degrees)
360.0
lower bound for final flight azimuth (0 <= azimuth <= 360; degrees)
0.0
final declination (-90 <= declination <= +90; degrees)

0.0
upper bound for final declination (-90 <= declination <= +90; degrees)
+90.0
lower bound for final declination (-90 <= declination <= +90; degrees)
-90.0
final east longitude (0 <= longitude <= 360; degrees)

55.0
upper bound for final east longitude (0 <= longitude <= 360; degrees)
360.0
lower bound for final east longitude (0 <= longitude <= 360; degrees)
0.0
The next series of data inputs define lower and upper bounds on the state variables during the ascent-to-
orbit phase. To free or unconstrain one or more states, set the lower and/or upper bounds to the value
1.0d99.
*************************************************************
upper and lower bounds on the flight conditions during ascent
*************************************************************
NOTE: set bound to 1.0d99 to ignore
-----------------------------------
upper bound for altitude (meters)

120000.0
lower bound for altitude (meters)

1.0
upper bound for velocity (meters/second)

2000.0
lower bound for velocity (meters/second)

0.1
upper bound for flight path angle (-90 <= fpa <= +90; degrees)
+90.0
lower bound for flight path angle (-90 <= fpa <= +90; degrees)
-90.0
upper bound for flight azimuth (0 <= azimuth <= 360; degrees)
360.0
page 5
lower bound for flight azimuth (0 <= azimuth <= 360; degrees)
0.0
upper bound for declination (-90 <= declination <= +90; degrees)
+90.0
lower bound for declination (-90 <= declination <= +90; degrees)
-90.0
upper bound for east longitude (0 <= longitude <= 360; degrees)
360.0
lower bound for east longitude (0 <= longitude <= 360; degrees)
0.0
This next section of the input file defines the initial guesses and bounds for the control variables. To
free or unconstrain one or more control variables set the lower and/or upper bounds to 1.0d99.
**********************************
initial flight controls and bounds
**********************************
to the initial value to constrain
or "fix" a flight control.
-------------------------------------
initial angle-of-attack (degrees)

0.0
upper bound for initial angle-of-attack (degrees)

+90.0
lower bound for initial angle-of-attack (degrees)

-90.0
initial bank angle (degrees)

0.0
upper bound for initial bank angle (degrees)

+90.0
lower bound for initial bank angle (degrees)

-90.0
initial throttle setting

1.0
upper bound for initial throttle setting

1.0
lower bound for initial throttle setting

0.5
This section of the input file defines the final guesses and bounds for the control variables. To free one
or more control variables set the lower and/or upper bounds to 1.0d99.
********************************
final flight controls and bounds
********************************
page 6
to the final value to constrain
or "fix" a flight control.
-------------------------------------
final angle-of-attack (degrees)

10.0
upper bound for final angle-of-attack (degrees)

+90.0
lower bound for final angle-of-attack (degrees)

-90.0
final bank angle (degrees)

0.0
upper bound for final bank angle (degrees)

+90.0
lower bound for final bank angle (degrees)

-90.0
final throttle setting

0.75
upper bound for final throttle setting

1.0
lower bound for final throttle setting

0.5
This section of the input file defines the lower and upper bounds for the control variables during the
ascent-to-orbit phase. To free or un-constrain one or more control variables, set the lower and/or upper
bounds to 1.0d99.
**********************************************************
upper and lower bounds on the flight controls during phase
**********************************************************
NOTE: set bound to 1.0d99 to ignore
-----------------------------------
upper bound for angle-of-attack (degrees)

+90.0
lower bound for angle-of-attack (degrees)

-90.0
upper bound for bank angle (degrees)

+90.0
lower bound for bank angle (degrees)

-90.0
upper bound for throttle setting

1.0
lower bound for throttle setting

0.5
page 7
These next three inputs allow the user to specify the final orbit characteristics. Please note that an input
of 1.0d99 for the final orbital inclination will tell the software to ignore the orbital inclination
constraint.
***********************
final orbit constraints
***********************
final semimajor axis (kilometers)

1838.0
final orbital eccentricity (non-dimensional)

0.0d0
final orbital inclination (degrees; enter 1.0d99 to ignore this constraint)

25.0
The software also creates a comma-separated-variable (csv) ascii data file that contains the optimal
control solution and other flight parameters. This next integer input specifies the type of solution data
file to create.
******************************************
type of comma-delimited solution data file
******************************************
1 = OC-defined nodes
2 = user-defined nodes
3 = user-defined step size
---------------------------
1
For options 2 or 3, this next input defines either the number of data points or the time step size of the
data output in the solution file.
number of user-defined nodes or print step size in solution data file
1.0
The name of this output file is specified in the next line of information. Please consult Appendix B for
information about the contents of this file.
name of comma-delimited solution data file
lv1.csv
The next series of program inputs are algorithm control options and parameters for the software. The
first input is an integer that specifies the type of collocation method to use during the solution process.
For most simulations, the trapezoidal method is recommended.
*********************************
discretization/collocation method
*********************************
1 = trapezoidal
2 = separated Hermite-Simpson
3 = compressed Hermite-Simpson
-------------------------------
1
The next integer defines the number of initial grid points to use in the collocation modeling of the ascent
trajectory.
number of initial guess grid points to use
25
page 8
The next series of program inputs are algorithm control options and parameters for the Sparse
Optimization Suite.
****************************
algorithm control parameters
****************************
This input defines the relative error in the objective function.

relative error in the objective function (performance index)
1.0d-5
The next input defines the relative error in the solution of the differential equations.
relative error in the solution of the differential equations
1.0d-7
The next input is an integer that defines the maximum number of mesh refinement iterations.
maximum number of mesh refinement iterations
20
The next input is an integer that defines the maximum number of function evaluations.
maximum number of function evaluations
100000
The next input is an integer that defines the maximum number of algorithm iterations.
maximum number of algorithm iterations
10000
The level of output from the NLP algorithm is controlled with the following integer input.
***************************
sparse NLP iteration output
***************************
1 = none
2 = terse
3 = standard
4 = interpretive
5 = diagnostic
---------------
2
The level of output from the optimal control algorithm is controlled with the following integer input.
Please note that option 4 will create lots of information.
**********************
optimal control output
**********************
1 = none
2 = terse
3 = standard
4 = interpretive
-----------------
1
The level of output from the differential equations algorithm is controlled with the following integer
input. Please note that option 5 will create lots of information.
page 9
****************************
differential equation output
****************************
1 = none
2 = terse
3 = standard
4 = interpretive
5 = diagnostic
---------------
1
The level of output can be further controlled by the user with this final text input. This program option
sets the value of the SOCOUT character variable described in the Sparse Optimization Suite user’s
manual. To ignore this special output control, input the simple character string no.
*******************
user-defined output
-------------------
input no to ignore
*******************
a0b0c0d0e0f0g0h0i0j2k0l0m0n0o0p0q0r0
Optimal control solution
Before attempting to optimize the trajectory, the software will display the conditions at the beginning
and end of the vertical rise as well as the flight conditions at the end of the pitch-over maneuver. The
following is the program output for this example.
liftoff
-------
time 0.000000000000000E+000 seconds

altitude 0.000000000000000E+000 meters
longitude 40.0000000000000 degrees
declination 20.0000000000000 degrees
velocity 0.000000000000000E+000 meters/second
fpa 0.000000000000000E+000 degrees
azimuth 0.000000000000000E+000 degrees
mass 3000.00000000000 kilograms
end of vertical rise

--------------------
time 10.0000000000000 seconds

altitude 2.31424039043486 meters
velocity 0.476409333678226 meters/second
delta-v 16.7072630675081 meters/second
end of pitch-over
-----------------
time 40.0000000000000 seconds

velocity 2.45921128624412 meters/second
fpa 76.9988433933634 degrees
azimuth 89.9999999999999 degrees
aoa 11.7511558278279 degrees
pitch 1.25000000000000 degrees
page 10
After convergence, the software will display the flight path conditions immediately after the pitch-over
and the flight path conditions and classical orbital elements at mission orbit injection. The following is
the program output for this example.
program lascent_ocs
--------------------
input data file ==> lascent1_max_payload.in
maximize final mass
flight path coordinates after pitch over

----------------------------------------
time 40.0000000000000 seconds
relative velocity 2.45921128624421 meters/second
relative fpa 76.9988433933634 degrees
relative azimuth 105.173026738430 degrees
vehicle mass 2941.73050211555 kilograms
final flight path coordinates

-----------------------------
time 1275.75471116451 seconds
relative fpa 1.785475748055718E-014 degrees
final mass and pitch-over to injection delta-v

----------------------------------------------
propellant mass 1551.04869893348 kilograms
mission orbit - orbital elements & state vector

-----------------------------------------------
sma (km) eccentricity inclination (deg) argper (deg)

0.183800000000D+04 0.123339131977D-14 0.250000000000D+02 0.000000000000D+00
raan (deg) true anomaly (deg) arglat (deg) period (min)

0.271286134821D+03 0.154251398995D+03 0.154251398995D+03 0.117848695428D+03
rx (km) ry (km) rz (km) rmag (km)

0.686318804299D+03 0.167132861969D+04 0.337447986892D+03 0.183800000000D+04
vx (kps) vy (kps) vz (kps) vmag (kps)

-.134883302933D+01 0.679412435911D+00 -.621701669512D+00 0.163323751027D+01
page 11
The following are trajectory characteristics plots created from the user-defined summary file. The first
plot illustrates the behavior of the angle of attack and bank angle during the ascent trajectory. Notice the
bank reversal as the angle of attack passes through zero. This behavior is explained on page 1192 of
“Advanced Launch System Trajectory Optimization”, by P. Daniel Burkhart and David G. Hull, AAS
96-168. An examination of the azimuth differential equation given in the Technical Discussion later in
this document will help explain this behavior.
The time evolution of altitude and velocity are illustrated in this next plot.
page 12
This plot illustrates the behavior of the orbital eccentricity and inclination during the ascent to the final
circular mission orbit.
The evolution of the relative flight path and azimuth angles is shown in this next plot. These angles are
computed with respect to a rotating Moon.
This final plot illustrates how the throttle setting and accumulated delta-v change during the ascent from
the end of the pitch-over maneuver to injection into the final mission orbit.
page 13
Verification of the optimal control solution
The optimal solution determined by the Sparse Optimization Suite can be verified by numerically
integrating the spacecraft’s selenocentric equations of motion using the optimal control-computed initial
conditions and the optimal control angles and throttle setting determined by the algorithm. This
computer program uses a Runge-Kutta-Fehlberg 7(8) variable step size method to explicitly integrate the
flight path equations of motion.
The following is the program output summarizing the final flight path coordinates, classical orbital
elements and state vector computed using this explicit numerical integration method. The vehicle mass
and accumulated delta-v are computed by including the instantaneous propellant flow rate and thrust
acceleration in the first-order equations of motion.
********************************************
* verification of optimal control solution *
********************************************

-----------------------------
time 1275.75471116451 seconds
relative fpa -1.066147781771430E-005 degrees

----------------------------------------------
page 14
Creating an initial guess
The software uses a linear interpolation guess based on the flight conditions after the pitch-over
maneuver and the final flight conditions provided by the user. This is accomplished within the software
by setting init = 1. The Sparse Optimization Suite attempts to obtain problem feasibility before
trying to solve the optimal control problem. If the software cannot get feasible, the user will have to
provide a better initial guess.
Problem setup
This section provides additional details about the Sparse Optimization Suite software implementation. It
briefly explains such things as point constraints and the performance index options. The lascent_ocs
simulation consists of a vertical rise followed by a pitch-over maneuver and a final, optimized ascent-to-
orbit propulsive phase.
(1) Performance index
The performance index for the maximize final mass option is the final spacecraft mass given by
J  m f  m t f 
where m f is the spacecraft mass at the final time.
For the minimize flight time program option, the performance index is
J  tf
where t f is the simulation time at mission orbit injection.
The value of the maxmin indicator in the Sparse Optimization Suite algorithm tells the software whether
the user is minimizing or maximizing the performance index.
(2) Path constraints
The path constraints enforced by the software are based on the lower and upper bounds for the dynamic
variables provided by the user. For example,
 L    U
 L    U
and so forth. The subscripts L and U correspond to lower and upper bounds provided by the user.
page 15
(3) Point constraints
The number and type of point functions enforced by the software is a function of the orbital elements of
the final mission orbit. This section discusses the point constraints or final boundary conditions
implemented in the lascent_ocs computer program.
Circular orbit
For a final circular mission orbit, the user provides the semimajor axis of this orbit. The software then
uses the following three point constraints or boundary conditions to “target” a final circular orbit.
 rf v f
rf  au vf  sin  f  0
au rf v f
If the user also specifies an orbital inclination for the final mission orbit, the software enforces the
following point constraint at the final time
sin i f 
h fz

r f  v f z
 sin iu
hf rf  v f
In these equations, the f subscript refers to mission characteristics computed by the software at the final
time and the u subscript refers to user-provide final or “target” orbital conditions.
Elliptical orbit
For a user-defined elliptical mission orbit, the software constrains the final specific orbital energy and
final angular momentum magnitude according to
2 2
 f  v 2f    u  vu2 
rf ru
rf  v f ru  v u
hf   hu 
rf  v f ru  v u
If the user also specifies an orbital inclination for the final mission orbit, the software enforces the
following point constraint
hf
sin i f  z 
rf  v f z  sin i
rf  v f
u
hf
In these equations, the f subscript refers to mission characteristics computed by the software at the final
time and the u subscript refers to user-provided final orbital conditions. The user-defined position
vector ru and velocity vector v u are determined from the semimajor axis, orbital eccentricity and
inclination input by the user. All other angular orbital elements are set to 0.
page 16
Technical Discussion
This section describes the equations of motion implemented in this computer program along with other
important numerical algorithms.
Flight path equations of motion
The first-order flight path equations of motion relative to a rotating, spherical Moon are as follows:
selenocentric radius
dr
r  V sin 
dt
selenographic longitude
d cos  sin
 V
dt r cos 
selenocentric declination
d cos  cos
 V
dt r
speed
dV T cos 
V   g sin    2 r cos   sin  cos   sin  cos  cos 
dt m
flight path angle
d V  T sin   g cos 
  cos     cos    2 sin cos 
dt r  mV  V
r
 2 cos   cos sin  sin   cos  cos  
V
flight azimuth
d V  T sin  
  tan  sin cos     sin 
dt r  mV cos  
r
 2  sin   cos cos  tan     2 sin cos  sin 
V cos 
where
page 17
r  selenocentric radius
V  relative speed
  flight path angle
  selenocentric declination
  selenographic longitude   east 
  flight azimuth   clockwise from north 
  bank angle   for a right turn 
  angle of attack
req  lunar equatorial radius
  lunar inertial rotation rate
  lunar gravitational constant
g  acceleration of gravity   r 2
T  propulsive thrust
m  spacecraft mass
Vertical rise phase
During the vertical rise part of the trajectory, the software uses a subset of the full equations of motion.
This set of equations includes r, v,  ,  with     0 . This is necessary in order to avoid singularities
in the flight path and azimuth differential equations during the trajectory evolution. Furthermore, during
the numerical integration of the vertical rise differential equations,   0 ,   0 ,   90 ,   0 and the
throttle setting is equal to 1.
Pitch-over phase
During the pitch-over part of the trajectory, the software also uses a subset of the full equations of
motion. This set of equations includes r, v,  ,  ,  with   0 . This is necessary because azimuth is
poorly defined until the end of the pitch-over part of the ascent trajectory. Furthermore, during the
numerical integration of the pitch-over differential equations,   0 ,   0 ,   90 ,  0 and the
throttle setting is equal to 1.
The angle of attack during the pitch-over phase is determined from   90      t  t p  where  is
the pitch rate, t is the time since launch and t p is the time since launch at which the pitch-over maneuver
begins. The constant pitch rate is determined from the pitch angle and pitch-over maneuver duration
provided by the user.
For a realistic and contiguous ascent trajectory, the launch vehicle should pitch over to the azimuth
direction determined by the ascent-to-final orbit optimization.
page 18
Crossrange and downrange calculations
The crossrange angle is determined from the following expression which can be derived from spherical
trigonometry relationships:
sin   sin1 sin 2 cos1 cos   cos1 cos1 sin   sin1 cos2 sin 1
The downrange angle is determined from the following three equations:
sin    cos 1 sin 2 cos 1 cos   sin 1 cos 1 sin   cos 1 cos 2 sin 1
cos   cos 2 cos 1 cos   sin 2 sin 1
  tan 1  sin  , cos  
where
1  selenocentric declination of the initial point
 1  flight azimuth at the initial point
2  selenocentric declination of the final point
  2  1
1  east longitude of the initial point

2  east longitude of the final point
Please note that the inverse tangent above is a four quadrant calculation.
Conversion of flight path coordinates to inertial state vector
The software uses the following expressions to calculate the selenocentric, inertial position and velocity
vectors of the spacecraft at the final time:
rx  r cos  cos 
ry  r cos  sin 
rz  r sin 
v x  v cos    cos sin  sin   cos  cos    sin sin  sin  
v y  v sin    cos sin  sin   cos  cos    sin sin  cos  
vz  v  cos cos  sin   cos  cos  
where the correction for lunar rotation is
page 19
v x  v x   ry
v y  v y   rx
The inertial speed can also be computed from the following expression
vi  v 2  2v r cos  sin cos   r 2 2 cos2 
and the inertial flight path angle from
v 2 cos2   2vr cos  cos cos   r 2 2 cos2 

cos  i 
v 2  2vr cos  cos cos   r 2 2 cos2 
where all coordinates on the right-hand-side of these equations are relative to a rotating Moon.
page 20
References and Bibliography
(1) “Direct Trajectory Optimization Using Nonlinear Programming and Collocation”, C. R. Hargraves
and S. W. Paris, AIAA Journal of Guidance, Control and Dynamics, Vol. 10, No. 4, July-August, 1987,
pp. 338-342.
(2) “Optimal Finite-Thrust Spacecraft Trajectories Using Direct Transcription and Nonlinear
Programming”, Paul J. Enright, Ph.D. Thesis, University of Illinois at Urbana-Champaign, 1991.
(3) “Using Sparse Nonlinear Programming to Compute Low Thrust Orbit Transfers”, John T. Betts, The
Journal of the Astronautical Sciences, Vol. 41, No. 3, July-September 1993, pp. 349-371.
(4) “Optimization of Thrust Direction Histories and Vehicle Parameters for Three-dimensional Ascent
Trajectories”, C. Tucker Battle and Robert G. Gottlieb, AIAA 64-663.
(5) “Improved Collocation Methods with Application to Direct Trajectory Optimization”, Albert L.
Herman, Ph.D. Thesis, University of Illinois at Urbana-Champaign, 1995.
(6) “Optimal Low Thrust Interplanetary Trajectories by Direct Method Techniques”, Craig A. Kluever,
The Journal of the Astronautical Sciences, Vol. 45, No. 3, July-September 1997, pp. 247-262.
(7) “Survey of Numerical Methods for Trajectory Optimization”, John T. Betts, AIAA Journal of
Guidance, Control and Dynamics, Vol. 21, No. 2, March-April 1998, pp. 193-207.
(8) “Design and Evaluation of a 3-D Optimal Ascent Guidance Algorithm”, Anthony J. Calise, Nahum
Melamed and Seungjae Lee, AIAA 97-3707.
(9) “Optimization of Launch Vehicle Ascent Trajectories with Path Constraints and Coast Arcs”, C. P.
Gath and A. J. Calise, AIAA Journal of Guidance, Control and Dynamics, Vol. 24, No. 2, March-April,
2001, pp. 296-304.
(10) “Advanced Launch System Trajectory Optimization Using Suboptimal Control”, Douglas A.
Shaver and David G. Hull, AIAA 90-3413.
(11) “Solving the Optimal Control Problem Using a Nonlinear Programming Technique, Part 2: Optimal
Shuttle Ascent Trajectories”, T. Bauer, J. Betts, W. Huffman and K. Zondervan, AIAA 84-2038.
(12) “Semi-Analytical Formulas for Launcher Performance Evaluation”, Emanuele Di Sotto and Paolo
Teofilatto, AIAA Journal of Guidance, Control and Dynamics, Vol. 25, No. 3, May-June, 2002, pp. 538-
545.
page 21
APPENDIX A
Contents of the Simulation Summary CSV File
This appendix is a brief summary of the information contained in the CSV data file produced by the
lascent_ocs software. The comma-separated-variable disk file is created by the odeprt subroutine
and contains the following information:
time (seconds) = simulation time since launch in seconds
altitude (meters) = altitude relative to a spherical Moon in meters
velocity (mps) = Moon-relative velocity in meters per second
flight path angle (d) = Moon-relative flight path angle in degrees
azimuth (deg) = Moon-relative azimuth angle in degrees
declination (deg) = selenocentric declination in degrees
longitude (deg) = selenographic east longitude in degrees
angle of attack (deg) = angle-of-attack in degrees
bank angle (deg) = bank angle in degrees
throttle = throttle setting
thrust (newtons) = thrust in newtons
load factor (g) = load factor in Earth g’s
crossrange (km) = crossrange distance in kilometers
downrange (km) = downrange distance in kilometers
vehicle mass (kg) = vehicle mass in kilograms
fuel flowrate (kg/s) = fuel flowrate in kilograms /second
rmag (meters) = selenocentric distance in meters
vmag (mps) = selenocentric speed in meters/second
semimajor axis (km) = semimajor axis in kilometers
eccentricity = orbital eccentricity (non-dimensional)
inclination (deg) = orbital inclination in degrees
arg of periapsis (deg) = argument of periapsis in degrees
raan (deg) = right ascension of the ascending node in degrees
true anomaly (deg) = true anomaly in degrees
periapsis radius (km) = selenocentric periapsis radius in kiometers
apoapsis radius (km) = selenocentric apoapsis radius in kilometers
eci fpa (deg) = inertial flight path angle in degrees
page 22
APPENDIX B
Fortran Functions and Subroutines
This appendix is a brief summary of the major Fortran functions and subroutines included in the
lascent_ocs computer program.
lascent_ocs.f - main executive program
atan3.for - four quadrant inverse tangent function
crdr.for - subroutine that calculates crossrange and downrange
eci2orb.for - subroutine that converts an inertial position and velocity vector to

classical orbital elements
fpc2car.for - subroutine that converts flight path coordinates to inertial

position and velocity vectors
fpeqms1.for - subroutine that evaluates flight path equations of motion during the
vertical rise phase
fpeqms2.for - subroutine that evaluates flight path equations of motion during the
pitch-over phase
fpeqms3.for - subroutine that evaluates flight path equations of motion for the
verification computations
odeinp.for - simulation input subroutine
odepf.for - point functions subroutine
odeprt.for - print subroutine – creates comma-separated-variable file
oderhs.for - subroutine that evaluates the equations of motion and any algebraic
equations
readfpn.for - read and echo floating point number from an input file subroutine
readint.for - read and echo an integer from an input file subroutine
readtext.for - read and echo text from an input file subroutine
rkf78.for - Runge-Fehlberg-Kutta (RKF78) numerical integration subroutine
rkf78cn.for - evaluate RKF78 integration coefficients subroutine
utility.for - number and text manipulation functions and subroutines
uvector.for - unit vector subroutine
vcross.for - vector cross product subroutine
vdot.for - vector dot product subroutine
vecmag.for - vector scalar magnitude function
xmod.for - modulo 2 pi function
page 23
APPENDIX C
Example Fortran Subroutine
This appendix contains the source code for a single Fortran 77 routine and illustrates typical
programming conventions used in the lascent_ocs software. This subroutine is the differential-
algebraic equations routine required by the Sparse Optimization Suite.
subroutine oderhs(iphase, t, ydyn, ny, p, np, frhs, nf, iferr)
c computes the right hand sides of the

c differential-algebraic (dae) equations
c input
c state variables
c rmag = selenocentric radius (meters) = ydyn(1)

c elon = selenographic longitude (radians) = ydyn(2)
c dec = selenocentric declination (radians) = ydyn(3)
c vrel = relative velocity (meters/second) = ydyn(4)
c fpa = relative flight path angle (radians) = ydyn(5)
c azim = relative azimuth (radians) = ydyn(6)
c mass = vehicle mass (kilograms) = ydyn(7)
c alfa = angle of attack (radians) = ydyn(8)
c beta = bank angle (radians) = ydyn(9)
c throttle = throttle setting (nd) = ydyn(10)
c output
c frhs = right hand sides of differential-algebraic system
c ************************************
implicit double precision (a-h, o-z)
include 'socxcom1.inc'
dimension ydyn(ny), p(np), frhs(nf), xtmp(2)
iferr = 0
c selenocentric radius (meters)
rmag = ydyn(1)
c selenographic longitude (radians)
elon = ydyn(2)
c selenocentric declination (radians)
dec = ydyn(3)
c relative speed (meters/second)
vrel = ydyn(4)
c flight path angle (radians)
fpa = ydyn(5)
c flight azimuth (radians)
page 24
azim = ydyn(6)
c current vehicle mass (kilograms)
xmass = ydyn(7)
c angle of attack (radians)
alfa = ydyn(8)
c bank angle (radians)
bank = ydyn(9)
c throttle setting (non-dimensional)
throttle = ydyn(10)
c compute gravitational acceleration (meters/second**2)
agrav = 1.0d9 * xmu / rmag**2
c *********************
c propulsion properties
c *********************
c propellant flow rate (kilograms/second)
xmdot = throttle * thrmax / vexhaust
c propulsive thrust (newtons)
thrust = throttle * thrmax
c *******************
c equations of motion
c *******************
c compute common trigonometric terms
sfpa = sin(fpa)
cfpa = cos(fpa)
sazim = sin(azim)
cazim = cos(azim)
salfa = sin(alfa)
calfa = cos(alfa)
sbank = sin(bank)
cbank = cos(bank)
sdec = sin(dec)
cdec = cos(dec)
c ****************************
c state differential equations
c ****************************
c no aerodynamic lift or drag
xlift = 0.0d0
drag = 0.0d0
page 25
rdot = vrel * sfpa
elondot = vrel * cfpa * sazim / (rmag * cdec)
decdot = vrel * cfpa * cazim / rmag
tmp1 = omega**2 * rmag * cdec * (sfpa * cdec

& - sdec * cfpa * cazim)
tmp2 = (thrust * calfa - drag) / xmass
tmp3 = -agrav * sfpa
vdot = tmp1 + tmp2 + tmp3
tmp1 = vrel * cfpa / rmag
tmp2 = 2.0d0 * omega * sazim * cdec
tmp3 = omega**2 * (rmag / vrel) * cdec * (cazim

& * sfpa * sdec + cfpa * cdec)
tmp4 = (thrust * salfa + xlift) * cbank / (xmass * vrel)
tmp5 = -agrav * cfpa / vrel
gamdot = tmp1 + tmp2 + tmp3 + tmp4 + tmp5
tmp1 = (vrel / rmag) * tan(dec) * sazim * cfpa
tmp2 = 2.0d0 * omega * (sdec - cazim * cdec * tan(fpa))
tmp3 = (rmag / (vrel * cfpa)) * omega**2

& * sazim * cdec * sdec
tmp4 = ((thrust * salfa + xlift) * sbank) / (xmass * vrel * cfpa)
azidot = tmp1 + tmp2 + tmp3 + tmp4
c *******************************
c first-order equations of motion
c *******************************
frhs(1) = rdot
frhs(2) = elondot
frhs(3) = decdot
frhs(4) = vdot
frhs(5) = gamdot
frhs(6) = azidot
frhs(7) = -xmdot
return
end
page 26
APPENDIX D
Minimum Time, Elliptical Final Orbit Example
This appendix summarizes the lascent_ocs solution and trajectory graphics for a minimum time
ascent to an elliptical final mission orbit.
page 27
Since this is a minimum time solution, the throttle setting during the pitch-over to final mission injection
phase is always 100 percent.
The input file for this example is named lascent1_min_time_elliptical.in.
page 28
Here’s the program output for this example.
program lascent_ocs
--------------------
input data file ==> lascent1_min_time_elliptical.in
minimize total time
flight path coordinates after pitch over

----------------------------------------
time 40.0000000000000 seconds
relative fpa 76.9988433933635 degrees

-----------------------------
time 1012.03207649040 seconds
relative fpa -0.581594381564705 degrees

----------------------------------------------

-----------------------------------------------

0.183800000000D+04 0.500000000000D-01 0.300000000000D+02 0.169425377108D+03

0.259056128591D+03 0.347734565247D+03 0.157159942355D+03 0.117848695428D+03

0.882748192401D+03 0.147009109148D+04 0.339251840931D+03 0.174800000000D+04

-.122640620953D+01 0.898934666349D+00 -.793719576539D+00 0.171526863533D+01
page 29
********************************************
* verification of optimal control solution *
********************************************

-----------------------------
time 1012.03207649040 seconds
relative fpa -0.581726988142458 degrees

----------------------------------------------

-----------------------------------------------

0.183801231387D+04 0.500080360095D-01 0.299999991538D+02 0.169426460294D+03

0.259056130199D+03 0.347733635676D+03 0.157160095970D+03 0.117849879739D+03

0.882743567977D+03 0.147009151140D+04 0.339249188092D+03 0.174799750295D+04

-.122641590194D+01 0.898931343902D+00 -.793724668909D+00 0.171527618057D+01
page 30

Lunar Ascent Trajectory Optimization

Uploaded by

Copyright:

Lunar Ascent Trajectory Optimization

Uploaded by

Document Information

Copyright

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Lunar Ascent Trajectory Optimization

Uploaded by

Copyright:

Program lascent_ocs

Lunar Ascent Trajectory Optimization

The important features of this scientific simulation are as follows:

 user-defined vertical rise and pitch-over phases

 collocation and implicit integration

please input the name of the simulation definition file

Input file format and contents

pitch over time duration (seconds)

pitch angle change during pitch over (degrees)

maximum thrust (newtons)

specific impulse (seconds)

initial vehicle mass (kilograms)

initial guess for final vehicle mass (kilograms)

initial flight azimuth (0 <= azimuth <= 360; degrees)

initial selenocentric latitude (-90 <= latitude <= + 90; degrees)

initial selenographic east longitude (0 <= longitude <= 360; degrees)

final time (seconds)

upper bound for final time (seconds)

lower bound for final time (seconds)

final altitude (meters)

upper bound for final altitude (meters)

lower bound for final altitude (meters)

final velocity (meters/second)

upper bound for final velocity (meters/second)

lower bound for final velocity (meters/second)

final declination (-90 <= declination <= +90; degrees)

final east longitude (0 <= longitude <= 360; degrees)

upper bound for altitude (meters)

lower bound for altitude (meters)

upper bound for velocity (meters/second)

lower bound for velocity (meters/second)

initial angle-of-attack (degrees)

upper bound for initial angle-of-attack (degrees)

lower bound for initial angle-of-attack (degrees)

initial bank angle (degrees)

upper bound for initial bank angle (degrees)

lower bound for initial bank angle (degrees)

initial throttle setting

upper bound for initial throttle setting

lower bound for initial throttle setting

final angle-of-attack (degrees)

upper bound for final angle-of-attack (degrees)

lower bound for final angle-of-attack (degrees)

final bank angle (degrees)

upper bound for final bank angle (degrees)

lower bound for final bank angle (degrees)

final throttle setting

upper bound for final throttle setting

lower bound for final throttle setting

upper bound for angle-of-attack (degrees)

lower bound for angle-of-attack (degrees)

upper bound for bank angle (degrees)

lower bound for bank angle (degrees)

upper bound for throttle setting

lower bound for throttle setting

final semimajor axis (kilometers)

final orbital eccentricity (non-dimensional)

final orbital inclination (degrees; enter 1.0d99 to ignore this constraint)

This input defines the relative error in the objective function.

Optimal control solution

time 0.000000000000000E+000 seconds

end of vertical rise

time 10.0000000000000 seconds

time 40.0000000000000 seconds