Monte Carlo

Another attempt at a monte carlo analysis, this time with LV2-like numbers.

This directory has an attempt at a managing JSBSim sessions that can spin up multiple threads to run simulations in parallel. It also kills the sim at apogee, whenever that may be. The upshot is that I can run many more simulations per minute than the naive approach.

This setup:

# starting values, real numbers for an N2501 motor
isp     =   179    # s
thrust  =  2510    # N
impulse = 15280    # N*s
CD      =     0.6

# Randomize
isp = random.gauss(isp, isp*0.03)                       # 3% spread
thrust = random.gauss(thrust, thrust*0.10)              # 10% spread
impulse = random.uniform(impulse - 200, impulse + 200)  # +/- 200 N*s
cd_avg = random.gauss(0.6, 0.1)                         # ~16% spread



In [ ]:

    
print("Ran 4000 indapendent simulations in 12.8 miuntes.")
print("That's about %0.2f seconds per simulation." % (720/4000))



In [ ]:

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

FPS2M = 0.3048
LBF2N = 4.44822
LBS2KG = 0.453592

n_times = 4000

max_alt = 0
max_alt_time = 0
max_alts = []
sim_times = []
sim_vel_ups = []
sim_thrusts = []
sim_alts = []
for i in range(n_times):
    # Read data from JSBSim
    sim_time = []
    sim_vel_up = []
    sim_alt = []
    sim_thrust = []
    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
            
            sim_time.append(time)
            sim_vel_up.append(-vel_down)
            sim_alt.append(alt)
            sim_thrust.append(thrust)

            # 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)

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

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 mean sim altitude was:   %5.1f km\n" % (mean / 1000.0))
print("The median sim altitude was: %5.1f km\n" % (median / 1000.0))



In [ ]:

    
fig, ax1 = plt.subplots(figsize=(18,7))
plt.title(r"Distribution of Maximum Altitudes")
plt.xlabel(r"Maximum Altitude [m]")

n, bins, patches = plt.hist(mlist, 40, histtype='step', normed=1, alpha=0.4, linewidth=1, fill=True)

#plt.ylim([0, 35])
#plt.xlim()
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)


plt.plot(sim_times[max_i], sim_alts[max_i], 'r-', lw=1.2, alpha=0.6, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_alts[min_i], 'g-', lw=1.2, alpha=0.6, label="Minimum Altitude Sim")
plt.plot(sim_times[median_i[0]], sim_alts[median_i[0]], 'k-', lw=1.4, alpha=0.9, label="Median Altitude Sim")

#plt.ylim([0, max_alt + 1000])
#plt.xlim([0, max_alt_time])
plt.legend(loc=2)
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_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.6, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_vel_ups[min_i], 'g-', lw=1.2, alpha=0.6, label="Minimum Altitude Sim")
plt.plot(sim_times[median_i[0]], sim_vel_ups[median_i[0]], 'k-', lw=1.4, alpha=0.9, label="Median 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 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_thrusts[i], 'b-', lw=1.0, alpha=0.005)

plt.plot(sim_times[max_i], sim_thrusts[max_i], 'r-', lw=1.2, alpha=0.6, label="Maximum Altitude Sim")
plt.plot(sim_times[min_i], sim_thrusts[min_i], 'g-', lw=1.2, alpha=0.6, label="Minimum Altitude Sim")
plt.plot(sim_times[median_i[0]], sim_thrusts[median_i[0]], 'k-', lw=1.4, alpha=0.9, label="Median Altitude Sim")

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