Rocket Challenge 2

This competition is now closed, but we still welcome all attempts at this challenge!
Please send your code to greg@bodekerscientific.com and we'd be delighted to run and test your code, and let you know how you fared compared with others.

The Challenge

Consider a spacecraft 100 km above the surface of some planet. It has just been released from the mothership and needs to descend to the surface carrying two astronauts. The spacecraft weighs 2000 kg and has 1000 kg of fuel onboard. The planet has a gravitational acceleration of 12 m.s^-2 (on Earth it is 9.8 m.s^-2). Burning 1 kg of fuel in one second generates 6000 N of upward thrust.

Your job is to write some python code that decides how much fuel to burn, every second, so that the spacecraft lands on the surface under the following constraints:

• The upward acceleration must never be more than 3g otherwise the astronauts will pass out.

• The vertical velocity as it touches down must be greater than 0 m.s^-1 and less than 2 m.s^-1 when the spacecraft lands otherwise you destroy it.

The goal is to land on the surface, under these constraints, with as much fuel in the tank as possible.

The twist:

The planet has an atmosphere. The air is more dense closer to the ground and less dense higher up. This air density will slow your spacecraft because of air friction. As the air becomes more dense so the friction will increase and slow your spacecraft even more. The frictional force depends on the square of the spacecraft velocity.

Ah, but if only life were so simple! This planet has a very turbulent atmosphere with lots of random waves. These waves create a rather bumpy ride down to the surface for your astronauts but also make it that your solution needs to adapt to the environment. Your algorithm needs to figure out what to do given the current altitude, speed and fuel left. So you will need to think.

Because the situation is slightly different each time you run the code, we will run your solution 1000 times and will average the fuel left in the tank across all of the runs to get your 'score'. If your algorithm fails more than 50 times out of the 1000 runs you will be disqualified. The person with the highest score will win.

Most of the job has already been done for you. Here is some sample code. All you need to do is add code where it is indicated below to generate the fuel to burn (ftb):

import matplotlib.pyplot as plt

import numpy as np

thrust_per_kg = 6000

g = 12.0

mass_of_fuel = 1000

mass_of_crew_module = 2000

time = [0]

alt = [100000.0]

actual_speed = [0.0]

acc = [g]

fuel_to_burn = [0.0]

def get_fuel_to_burn(time, alt, speed, fuelmass, density):

ftb = 0.0

# Students should implement their solution between here

ftb = 2.0

# and here

# The line below prevents more fuel being burnt than is available

ftb = ftb if ftb < fuelmass else fuelmass

ftb = max(ftb, 0)

return ftb

def generate_atmos_density():

"""

This function generates the atmospheric density profile for the planet. The density profile decays exponentially

from the surface to the top of the atmosphere of the planet. The density profile will affect the

frictional force that the spacecraft experiences. The profile includes some random waves.

:return: an array of 100,001 elements specifying the density from the surface to 100 km altitude every metre.

"""

altitude = np.arange(0,100001,1)

dens = 100000 * np.exp(-altitude/7000)

# Now generate four 'waves' in the density profile. Each wave will have a randomly generated amplitude, phase and

# wavelength within some bounds. These waves are added up and then applied to the density profile

waves = np.zeros((100001,), dtype = np.float )

for wv in range(1, 5):

amplitude = np.random.normal(wv * 0.1, wv * 0.05, 1)

# The line below randomly selects the wavelength in metres

wavelength = np.random.normal(wv * 4000, wv * 1000, 1)

# The line below randomly generate a phase offset

phase = wavelength * np.random.random()

# Here is the equation that generates the wave

wave = amplitude * np.cos(2 * np.pi * (altitude - phase)/wavelength)

# Add the four waves

waves = waves + wave

# If you are in debugging model it can be interesting to plot 'waves' at this point

# Now apply the waves

dens = dens *(1 + waves)

return dens

#-----------------------------------------------------------------------------------------------------------------------

# Main program starts here

#-----------------------------------------------------------------------------------------------------------------------

dens = generate_atmos_density()

print('Time, Altitude, Speed, Accel., G-force, Fuel left')

print('{:4d}, {:8.1f}, {:6.2f}, {:6.2f}, {:6.2f}, {:6.2f}'.format(time[-1],

alt[-1], actual_speed[-1], acc[-1], (g - acc[-1]) / 9.8, mass_of_fuel))

try:

while alt[-1] > 0:

# Find out how much fuel we will burn in the next second.

fuel_to_burn.append(get_fuel_to_burn(time[-1], alt[-1], actual_speed[-1], mass_of_fuel, dens[int(alt[-1])]))

# Calculate the deceleration from the fuel burnt in the previous time step

acc_from_burn = (thrust_per_kg * fuel_to_burn[-1]) / (mass_of_fuel + mass_of_crew_module)

# Calculate the deceleration from friction which depends on the density profile. The frictional force depends on

# velocity squared and the density of the atmosphere: acc = force/mass

acc_from_friction = actual_speed[-1]**2 * dens[int(alt[-1])] * 5E-7 / (mass_of_fuel + mass_of_crew_module)

# Update the mass of fuel

mass_of_fuel = mass_of_fuel - fuel_to_burn[-1]

# Work out the acceleration at the end of the 1 second period.

acc.append(g - acc_from_burn - acc_from_friction)

# Work out new speed: v = u + a*t

actual_speed.append(actual_speed[-1] + acc[-1])

# Work out new altitude: s = ut + 0.5*a*t^2

alt.append(alt[-1] - (actual_speed[-2] + 0.5 * acc[-1]))

# Update the time

time.append(time[-1] + 1)

# Print the diagnostics

print('{:4d}, {:8.1f}, {:6.2f}, {:6.2f}, {:6.2f}, {:6.2f}'.format(time[-1],

alt[-1], actual_speed[-1], acc[-1], (g-acc[-1])/9.8, mass_of_fuel))

if (g-acc[-1])/9.8 > 3.0:

raise Exception('G-force exceeded 3.0.')

landing_speed = actual_speed[-1]- (actual_speed[-2] - actual_speed[-1])/ (alt[-2] - alt[-1]) * alt[-1]

if (landing_speed < 0.0) or (landing_speed > 2.0):

raise Exception('speed on landing was too high.')

print('Well done! you landed successfully.')

print('Landing speed was {:6.4f} m.s^-1'.format(landing_speed))

print('Final fuel left in tank={:6.3f} kg'.format(mass_of_fuel))

except Exception as err:

print('Sorry your solution failed because {}'.format(err))

# Make plot of altitude and speed

fig, ax = plt.subplots()

ax.plot(time, alt, color='red')

ax.set_xlabel('Time (seconds)', fontsize=12)

ax.set_ylabel('Altitude (m)', color='red', fontsize=12)

ax2=ax.twinx()

ax2.plot(time, actual_speed, color='blue')

ax2.set_ylabel('Speed (m/s)', color='blue', fontsize=12)

plt.show()

plt.close()

# Make plot of acceleration and fuel burnt

fig, ax = plt.subplots()

ax.plot(time, acc, color='red')

ax.set_xlabel('Time (seconds)', fontsize=12)

ax.set_ylabel('Acceleration (m/s^2)', color='red', fontsize=12)

ax2=ax.twinx()

ax2.plot(time, fuel_to_burn, color='blue')

ax2.set_ylabel('Fuel burnt (kg)', color='blue', fontsize=12)

plt.show()

plt.close()

#-----------------------------------------------------------------------------------------------------------------------

# Main program ends here

#-----------------------------------------------------------------------------------------------------------------------

Clearly what is written above is not a successful algorithm - the spacecraft hits the surface at well over 2 m.s^-1 with the result that the astronauts need to be scraped off the floor with a spatula. Your job is to modify the algorithm between the lines indicated in the code above so that the landing occurs within the constraints detailed above. The winner will be the person whose algorithm lands with the most fuel left in the tank, within the constraints.