39.1 The Process and Analysis

Yes, it is rocket science; although simplified here, it is adequately rich. The objective is to schedule thrust (e.g., change thrust with time) to maximize the altitude a rocket can achieve. You could have full thrust until the fuel burns out, and then let the rocket coast up until gravity pulls it back. But this plan makes the rocket go fast in the lower altitude where the air resistance is high and wastes energy. A better plan to minimize the impact of air resistance is to move slowly in the low altitudes and then full throttle to fast speed in the high altitudes. But moving slow in the low altitudes, in the high gravitational field, means wasting fuel to fight gravity. In a limit of not enough thrust to overcome gravity, the rocket never rises; it burns all its fuel on the launch pad. This optimization is known as the “Goddard problem” in honor of rocket scientist Robert Goddard.

This will consider a simple situation of a single‐stage vertical rocket (no changes in drag coefficient), moving vertically from a stationary flat earth (no Coriolis effects) through an atmosphere with nominal density dependence on height. Even with such idealizations, it will be sufficiently complex!

One generally accepted solution is to accelerate the rocket at full thrust to terminal velocity (the speed at which it would fall if pulled down by gravity when air drag resistance exactly balanced gravity) and then reduce thrust to maintain v_terminal. As the rocket rises both air density and gravitational pull reduce, so thrust would be scheduled with the local v_terminal. Then coast after burnout until gravity and air drag stop upward motion. Although the optimum velocity trajectory seems to be a bit lower than terminal velocity, this rule reveals the diversity of variables that could be used to schedule thrust.

There are many versions for this exercise. Considering options on the DV, thrust (rate of fuel consumption) could be scheduled w.r.t. any number of parameters that change in time: time, elevation, terminal velocity, speed, remaining fuel, or air density. The scheduling could be a continuum relation (e.g., ) or a rule‐based type of model or .

Further, the OF could be any of several: Perhaps maximize elevation at the top of its fuel‐depleted upward drift or perhaps reach a given elevation with minimum fuel consumption.

Do you schedule the throttle with height, or time, or remaining fuel, or terminal velocity? This will be your choice. Let’s consider height, and several options. (i) Should you section height into 10,000 ft intervals and find an optimal throttle position for each interval? This would step the throttle to the new position at each 10k ft interval. (ii) Why should the transition heights be in 10k ft intervals, and why should the intervals be uniform? Could the optimizer find both the transition heights and the thrust values? (iii) Or should you have an equation that defines throttle position with elevation and use the optimizer to determine the coefficients in the equation? This would provide for smooth throttle transitions as elevation increases. What equation complexity would you use? What would you use for an equation form to ensure that throttle is bounded between 0 and 100%? (iv) Should you choose a combination of scheduling techniques, for instance, set the initial thrust as 100%, and then trigger an equation after some height? Let the optimizer determine the transition height and the post‐transition equation coefficients.

This set of questions about the DV schedule will be applicable to any schedule parameter, of either elevation, or time, or fuel content. In any case, what would you do when fuel runs out, but the equation or schedule says go to 83% thrust? And how can you know how to define stages with elevation until you know how high the rocket will rise?

The model of the rocket can be found at http://bocop.saclay.inria.fr/?page_id=465. Thanks to Thomas Hays for pointing me to it while he was a PhD candidate in aerospace engineering. The situation can be illustrated as follows: Figure 39.1 shows the Earth with radius r₀ and the rocket rising vertically from the surface with an elevation (from the Earth center) as r. The rocket path is solid during the time when it had fuel and was thrusting and dashed when the fuel is spent and it is coasting away due to momentum. The elevation hits a peak when the combined gravitational pull and atmosphere drag reduces the upward v to zero and it begins to free fall back to Earth. The figure is not to scale, represents a concept, and illustrates a rocket that will reach more than an Earth diameter from the surface. However, in this case study, the rocket only reaches about 50 miles up, which is only about six‐thousandths of a diameter.

Figure 39.2 illustrates how thrust might be scheduled with elevation. There are many ways to do this; for example, thrust could be scheduled with time, velocity, fuel content, etc. This schedule starts the thrust at 100% and then steps the thrust to certain values at certain heights (waypoints) and holds that thrust value until the next waypoint is reached. After the third waypoint, the u₄ thrust is held until fuel is depleted. When the fuel is spent, the thrust is zero. Here, the optimizer could choose the three waypoint elevations and the four thrust values, a total of 7 DVs. But it could be divided into 10 sections, giving the optimizer 19 DVs. Or the parameter for the schedule could be any other variable. Here the schedule is discontinuous. Alternately, the schedule could be defined by an equation and based on rocket fuel inventory, such as , or a ratio of polynomials such as . In either polynomial case, the optimizer would have four DVs, the schedule coefficients a, b, c, d.

The model can be developed as follows: Start with Newton’s law of motion for the rocket . Note that both F and m will change in time.

The three forces (thrust, gravity, and drag) will change in time. Thrust, , will be a fraction of maximum, full thrust, defined by the control signal u, which changes in time. Gravitational pull will be dependent on the rocket mass and distance from the Earth center, . Air drag will use the conventional (and elementary) turbulent drag on a body immersed in a large fluid field, , using a conventional (elementary) model for how air density depends on elevation , where elevation and ρ₀ is the air density at r₀ sea level. Note that h and v both change in time.

The rocket mass changes with fuel consumption, ideally modeled as a linear response to the control signal , .

The velocity and acceleration are related, and the initial velocity is zero. , . Similarly, position and velocity are related. , , or .

Combining relations in Newton’s law of motion and restating the position and mass models,

(39.1)

(39.2)

(39.3)

The v|v| term could be expressed as v². As v|v| the sign of the drag force works whether the velocity is upward or downward. Since this study only models the rocket on its upward journey, I will use v².

For convenience and generality, convert to scaled variables. Some are logically scaled by their initial values and . One is logically scaled by the maximum value . Velocity and time can be similarly scaled, and , but here the scaling values are not the initial values of zero. Convenient scaling values for v₀ and t₀ can be revealed when the three ODEs are converted to scaled variables. Choosing and , the model becomes

(39.4)

(39.5)

(39.6)

According to the Bocop site, appropriate coefficient values for a rocket on Earth are . Here are some un‐scaling values: miles, ft/s², s, and ft/s. Use the initial fuel mass as 90% of the full rocket. This means that when , the fuel is spent.

There are many ways to solve the three ODEs. In my solution, I use the Euler method, a forward finite difference explicit approximation to the derivative with . A more proper value is . The small integration time step is needed to make the solution independent of it; however, it takes a long time to run, and given the idealizations in the model and precision on the model coefficients, the coarse consequence of the numerical method is probably inconsequential. I have also investigated primitive predictor–corrector and Runge–Kutta techniques to solve the ODEs, but again, I don’t think that the model and coefficient values justify numerical methods that seek analytically exact solutions. So, I defaulted to the simple Euler method. The numerical solution is

(39.7)

(39.8)

(39.9)

where i is the time step counter.

Start with the initial conditions, and do until (solve the ODE set until) , which is the indication that the rocket has reached its peak and started to fall back. At that point use as the OF value.

Function 79 (in Appendix E) and Functions 77 and 78 (in the 2‐D Optimization Examples file on the companion site) explore three variations on the problem. They also show the finite difference‐approach to modeling the solution to the three time‐dependent differential equations (which are coupled and nonlinear). These functions have only 2 DVs. In Function 77 the DVs are initial and second‐stage thrust. The second‐stage transitions at 40k ft and lasts until fuel runs out. In Function 78 the initial thrust is 100% and the decision variables are the elevation to switch to the second thrust level and that level. Function 79 has an alternate objective—reach a desired elevation (about 100k ft), minimizing fuel content. The thrusts are 100% and then 20%, and the DVs are the elevations to transition from first stage to second and then second to thrust off (coast).

The OF has flat spots, cliff discontinuities, ridge discontinuities due to numerical discretization, and hard constraints.

39.2 Pre‐Assignment Note

I explored several approaches to the DVs and would offer this guide, which minimizes the number of variables and finds nearly the max height with the minimum number of function evaluations:

• DV1—the initial thrust, which should be 1.
• DV2—the rocket mass when switching to a new thrust, the trigger to change from the initial thrust to the schedule. My optimizations place this in the m′ = 0.6 region.
• DV3 and DV4—the coefficients in a linear relation for thrust as a function of rocket mass after the initial thrust: u′ = DV3 + DV4 * (DV2 − m׳). Here DV3 would be the jump‐to thrust at the point of change.

I chose to schedule the thrust with respect to rocket mass, because I know that it starts at m׳ = 1 and ends at m׳ = 0.1. Although the thrust could be scheduled w.r.t. either height or time, you don’t know what the range of r׳ or t׳ is. Perhaps the thrust could be scheduled with terminal velocity. As I have explored higher‐order models for thrust as a function of mass, I find that they describe nearly a linear relation.

The choice of DVs progressed with my investigation. My first attempts were to have 3 DVs representing thrust step and holds with height. When it was working, I extended it to 10 and then to 20 thrusts. This made me realize that I could not know the height range for the thrusts when the fuel runs out. So, I changed to schedule thrust with fuel content. The thrust patterns showed equivalently high thrust for the first several stages and then a sharp drop to lower thrusts that seemed all about the same. Looking at this pattern, I decided to use the DVs as aforementioned. I then decided to increase the thrust equation to a cubic, which resulted in 6 DVs. Then it seemed that the optimizer indicated that the initial thrust should be 100%. So, I set it to 100% and eliminated a DV. Very pleased with progress, I explored various models and began to realize that a linear schedule seems about as effective as the cubic, which reduces the number of DVs to 3. Sketch, evaluate, erase, sketch, evaluate, erase, …, ink.

The best strategy is to use full throttle to accelerate to a “high” velocity and then back off power; and as elevation increases, then increase throttle as air resistance diminishes.

My optimum schedule gets the rocket to an altitude of about 50 miles, which is very small relative to the 3960‐miles radius of the Earth, making the Cartesian viewpoint seems reasonable. It peaks at 3.5 min. A characteristic velocity is about 1500 ft/s, which means that the rocket moves about 12 ft in one simulation discretized time step.

Doing this in a simulation, I have not had to make plans for where the rocket will land when it free‐falls back.

39.3 Exercises

Use at least 4 DVs. Perhaps start with just 2 DVs and use the 2‐D program to practice on. Then use the Generic LF optimization software to determine a solution for 4 or more DVs.

1. Explain the situation and basic equations in your own words.
2. Explore the surface in a 2‐D simulation and comment on the several aberrations.
3. Explain your choice for the scheduling of the DVs: Discrete values staged by regular intervals of height, fuel, velocity, mass, or time; or variable stage interval defined by the optimizer; or an equation to schedule the DV continuously w.r.t. some parameter; or some combination.
   1. Choose a 2‐DV simple situation.
   2. Choose an N‐DV situation.
4. Explain your finite difference model and the triggers that you will use and actions that the code will take to recognize when the rocket hits the peak altitude and when the fuel is spent.
5. Implement your N‐DV optimization in my base case version “Generic LF Optimization Goddard Problem.” Clearly reveal your unique code and your unique DV schedule. Discuss the results.
6. Make the simulation go to the next step of rigor by using a better model for air density and/or air drag.
7. Explore the impact of uncertainty on model coefficient values and air density (which changes day to day with barometric pressure, humidity, and temperature) on the optimum height.
8. Investigate the OF response to DV values for any pair of DVs in the vicinity of the optimum. Explore a DV range of of the DV* value. You will see that the OF surface has flat spots. Explain why.