This piece of work will count for 50% of the overall mark for MATH6120. It is intended that you carry out the analysis using the AMPL modeling software system.

 Please refer to the university calendar for the penalties applied to late assignments.

Your report should be word processed. Ensure that you take frequent and multiple backups of your work, since excuses concerning lost or corrupted files will not be treated sympathetically.

As a very rough guide, I would like about six to ten A4 sides of description, with any extra material (e.g., diagrams, tables or computer output) attached at the end. Your report should include:

• an introduction;

• the formulation of your model with a statement of any assumptions made;

• a brief description of how you solved it;

• a presentation of the results, including any illustrations or graphs that you found useful;

• a summary of recommended actions that should be followed (and possible sugges- tions for further investigations)

• printouts of the AMPL model(s) and corresponding output in one or more appen- dices.

It should be written from the point of view of you, the OR analyst, describing some work that you have done to a space engineer (who will know the mathematical program- ming terminology).

Marks will be awarded on the basis of

• quality of the analysis;

• clarity of the presentation;

• originality.

Your work must be written up independently. The model and report must be your own work. You are reminded of the University policy on plagiarism. Your report must acknowledge clearly all people with whom you have dis- cussed any part of the coursework, as well as any references you might have used.

Shooting Rockets

Elon Musk has a problem. His Falcon Heavy launch vehicle, designed to carry humans into space, does not reach the height necessary for successful operations. Preliminary analysis has indicated that the issue stems from an underperforming first stage, more precisely from a suboptimal flow of fuel into the combustion chamber of the rocket engine. You are tasked with optimising the fuel flow in order to maximise the height of the rocket at the end of the first stage burn.

Figure 1: The Falcon Heavy at launch. Malfunctioning rockets are known to spontaneously disassemble. Especially if the Maths isn’t right.

The underlying dynamics of the flight can be described as follows. Launch time is set to t0  := 0. Let tfinal  be the time when the rocket finishes firing (e. g. because all fuel has been used.). This will be one of many decision variables. We consider the following functions

h : [t0 ,tfinal] −→ ❘ ,

v : [t0 ,tfinal] −→ ❘ ,

m : [t0 ,tfinal] −→ ❘ ,

θ : [t0 ,tfinal] −→ ❘ .

At any time t ∈ [t0 ,tfinal], we have that h(t) is the height of the rocket above ground, v(t) is the speed of the rocket, m(t) is the mass of the rocket, and θ(t) is the thrust of the rocket. The thrust is proportional to how much fuel is pumped into the reaction chamber, so we will use θ as one of the functions we need to determine.

Now, the following equations describe the behaviour of the rocket for all times t ∈ [t0 ,tfinal]:

h′ (t)   =   v(t),                                                                            (1)

v′ (t)   =    −   2 ,                                (2)

m′ (t)   =   −2θ(t)                                                                        (3)

with the auxiliary function

D(h(t),v(t)) := D0  × (v(t))2  × exp  −h0   ,                      (4)

(You do not need to understand all the details of these equations. However, it helps debugging your code if you have a rough understanding of what is going on here. For example, suppose that the time t is measured in seconds and the height h(t) in meters. Then clearly the first derivative of h(t) has dimensions m/s, i. e. h′ (t) measures a velocity. Equation  (1) simply tells us that the rocket is moving away from the surface of the earth, as fast as v(t) prescribes. Equation (3) tells us that the change in mass, m′ (t) is proportional to what what fuel is pumped into the rocket engine, i. e. how much exhaust leaves the rocket. Equation (2) tells us about the rate of change in velocity, v′ (t), i. e. the acceleration of the rocket. The acceleration is negatively effected by the gravitational effect of the earth at height h(t), measured in −(h(0)/h(t))2 . Apart from gravitational effects, the acceleration is inverse proportional to the mass m(t) and proportional to the term θ(t)−D(h(t),v(t)). The first part of this is the fuel flow again, the second part given by the function D describes friction with the atmosphere, which decays exponentially with the height h(t) and is otherwise proportional to the velocity v(t). The constants D0 and h0 are defined below.)

We also have “boundary conditions”

h(0)   =   1,                                                       (6)

m(0)   =   1,                                                       (7)

m(tfinal )   =   0.6,                                                    (8)

that is, the rocket starts at rest (v(0) = 0), we normalised height and mass to one unit, and overall 40% of the mass of the rocket is fuel. Furthermore, there are the simple

bounds

h(t) ≥ h(0),                                                (10)

m(tfinal ) ≤ m(t) ≤ m(0),                                                (11)

0 ≤ θ(t) ≤ θmax                                                                         (12)

that need to hold for all t ∈ [t0 ,tfinal]. Finally, constant parameters are defined as

θmax     :=   3.5,

h0     :=   500,

v0     :=   620,

D0     :=   v0 /2.

With this, the overall problem becomes

max    tfinal,h,v,m,θ

subject to

h(tfinal )

(1)– (12).

However, this is not a problem we can readily solve, as some of the decision variables are functions  h,v,m,θ . We would need to determine the function values at all times t ∈ [t0 ,tfinal]. In a sense, the above is an optimisation problem with an infinite number of variables.

To get around that difficulty, we will discretize the time interval [t0 ,tfinal] into a pre- scribed number of subintervals of same size. Let N be the number of these subintervals, i. e. we have as their boundary points

0                   1                   2                           k                            N

and from now on we are only interested in the function values h(tk ),v(tk ),m(tk ),θ(tk ) at all of these grid points. Collecting these decision variables into vectors

 := (v(t0 ),v(t1 ), . . . ,v(tfinal ))⊤ ∈ ❘N+1                                and similarly for h,m,θ, gives us the unknowns tfinal  ∈ ❘ ,  , ,  , θ¯ ∈ ❘N+1 .

The only problem now is the system (1)– (3). For this, we use the same trick as in the section on numerical differentiation: replace the slope h′ (t) of the tangent of the graph of h at t = tk  by the slope of a secant and replace (1) by

h(t1 ) − h(t0 ) tfinal /N

h(tk+1) − h(tk−1) 2tfinal /N

h(tN ) − h(tN−1) tfinal /N

=   v(t0 ),

=   v(tk )

=   v(tN ),

(k = 1, . . . ,N − 1),

(13)

(14)

(15)

i. e. in terms of our new variables, and slightly simplified

N(1 − 0 )   =   tfinal  × 0 ,

N(k+1 − k−1)   =   tfinal  × 2k            (k = 1, . . . ,N − 1),

N(N  − N−1)   =   tfinal  × N .

The equations (2) and (3) can be replaced in the same fashion.

We end up with a nonlinear optimisation problem of the form

max

tfinal∈❘, ,, ,θ¯∈❘N+1 s.t.

N

N(1 − 0 ) = tfinal  × 0 ,        N(k+1 − k−1) = 2tfinal  × k N(N  − N−1) = tfinal  × N ,

. . . ,

....

(k = 1, . . . ,N − 1),