Vertical Shoot Problem

Given a rocket and a fixed amount of fuel, what thrust curve gets the most altitude?

In each of our Monte Carlo runs the rocket is the exact same size and mass. The amount of propellent is the same as well as the Isp. The only thing is that changes is the thrust curve.



In [ ]:

    
# Chart
import csv
import statistics
import matplotlib.pyplot as plt
from glob import glob
%matplotlib inline

FPS2M = 0.3048
LBF2N = 4.44822
LBS2KG = 0.453592

n_times = len(glob('data/sim-*.csv'))

max_alt = 0
max_alt_time = 0
max_alts = []
sim_times = []
sim_vel_ups = []
sim_thrusts = []
sim_alts = []
sim_accs = []
sim_qbars = []
for i in range(n_times):
    # Read data from JSBSim
    sim_time = []
    sim_vel_up = []
    sim_alt = []
    sim_thrust = []
    sim_acc = []
    sim_qbar = []
    this_max_alt = 0
    with open("data/sim-%05d.csv" % i) as datafile:
        reader = csv.reader(datafile, delimiter=',')
        for row in reader:
            # ignore first line
            if row[0][0] == 'T':
                continue
            time      = float(row[0])           # s
            alt       = float(row[1])           # m
            vel_down  = float(row[2]) * FPS2M   # m/s
            thrust    = float(row[3]) * LBF2N   # N
            qbar      = float(row[4])           # psf
            acc       = float(row[5]) * FPS2M   # m/s/s
            
            sim_time.append(time)
            sim_vel_up.append(-vel_down)
            sim_alt.append(alt)
            sim_thrust.append(thrust)
            sim_qbar.append(qbar)
            sim_acc.append(acc)

            # max alt
            if alt > max_alt:
                max_alt = alt
                max_alt_time = time
            if alt > this_max_alt:
                this_max_alt = alt
    max_alts.append((i, this_max_alt))
    sim_times.append(sim_time)
    sim_vel_ups.append(sim_vel_up)
    sim_alts.append(sim_alt)
    sim_thrusts.append(sim_thrust)
    sim_qbars.append(sim_qbar)
    sim_accs.append(sim_acc)

sorted_max = sorted(max_alts, key=lambda tup: tup[1])
mlist = [tup[1] for tup in sorted_max]

print("I ran %d simulations.\n" % n_times)

print("Highest sim went to:         %5.1f km\n" % (mlist[-1] / 1000.0))
print("Lowest sim went to:          %5.1f km\n" % (mlist[0] / 1000.0))

mean = statistics.mean(mlist)
median = statistics.median(mlist)

max_i = sorted_max[-1][0]
min_i = sorted_max[0][0]
median_i = min(sorted_max, key=lambda tup: abs(tup[1]-median))

print("The median sim altitude was: %5.1f km\n" % (median / 1000.0))



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Simulated Rocket Thrust: All Runs")
plt.ylabel(r"Thrust [N]")
plt.xlabel(r"Time [s]")

for i in range(len(max_alts)):
    plt.plot(sim_times[i], sim_thrusts[i], 'b-', lw=1.0, alpha=0.005)

for i in range(10):
    plt.plot(sim_times[sorted_max[-i][0]], sim_thrusts[sorted_max[-i][0]], 'c-', lw=1.2, alpha=0.3)

    
plt.plot([-1,-2], [0,0], 'c-', lw=1.2, alpha=0.3, label="Top 10 Simulations")
plt.plot(sim_times[max_i], sim_thrusts[max_i], 'r-', lw=1.4, alpha=0.8, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_thrusts[min_i], 'g-', lw=1.2, alpha=0.5, label="Minimum Altitude Sim")

#plt.ylim([0, 600])
plt.xlim([-1, 60])
plt.legend(loc=1)
plt.show()



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Simulated Rocket Altitude: All Runs")
plt.ylabel(r"Altitude MSL [meters]")
plt.xlabel(r"Time [s]")
    
for i in range(len(max_alts)):
    plt.plot(sim_times[i], sim_alts[i], 'b-', lw=1.0, alpha=0.005)
    
for i in range(10):
    plt.plot(sim_times[sorted_max[-i][0]], sim_alts[sorted_max[-i][0]], 'c-', lw=1.2, alpha=0.3)

    
plt.plot([0,0], [0,0], 'c-', lw=1.2, alpha=0.3, label="Top 10 Simulations")
plt.plot(sim_times[min_i], sim_alts[min_i], 'g-', lw=1.2, alpha=0.4, label="Minimum Altitude Sim")
plt.plot(sim_times[max_i], sim_alts[max_i], 'r-', lw=1.2, alpha=0.8, label="Maximum Altitude Sim")

#plt.ylim([0, max_alt + 1000])
#plt.xlim([0, 60])
plt.legend(loc=2)
plt.show()



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Simulated Rocket Velocity: All Runs")
plt.ylabel(r"Velocity [m/s]")
plt.xlabel(r"Time [s]")
    
for i in range(len(max_alts)):
    plt.plot(sim_times[i], sim_vel_ups[i], 'b-', lw=1.0, alpha=0.005)

plt.plot(sim_times[max_i], sim_vel_ups[max_i], 'r-', lw=1.2, alpha=0.8, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_vel_ups[min_i], 'g-', lw=1.2, alpha=0.4, label="Minimum Altitude Sim")

#plt.ylim([0, 600])
#plt.xlim([0, max_alt_time])
plt.legend(loc=1)
plt.show()



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Simulated Rocket Acceleration: All Runs")
plt.ylabel(r"Acceleration [m/s/s]")
plt.xlabel(r"Time [s]")
    
for i in range(len(max_alts)):
    plt.plot(sim_times[i], sim_accs[i], 'b-', lw=1.0, alpha=0.005)

plt.plot(sim_times[max_i], sim_accs[max_i], 'r-', lw=1.2, alpha=0.8, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_accs[min_i], 'g-', lw=1.2, alpha=0.4, label="Minimum Altitude Sim")

#plt.ylim([0, 200])
plt.xlim([0, 60])
plt.legend(loc=1)
plt.show()



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Simulated Rocket QBar: All Runs")
plt.ylabel(r"Qbar [psf]")
plt.xlabel(r"Time [s]")
    
for i in range(len(max_alts)):
    plt.plot(sim_times[i], sim_qbars[i], 'b-', lw=1.0, alpha=0.005)

plt.plot(sim_times[max_i], sim_qbars[max_i], 'r-', lw=1.2, alpha=0.8, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_qbars[min_i], 'g-', lw=1.2, alpha=0.4, label="Minimum Altitude Sim")

#plt.ylim([0, 600])
plt.xlim([0, 80])
plt.legend(loc=1)
plt.show()