Multi-Disciplinary Design Optimization

Abstract

This is a program for generating an optimized rocket design based on a feasible initial design, some fixed parameters, physical constants, and a theoretical rocket model. Mass is minimized while satisfying a number of mission constraints.

A feasible design is successfully produced and analyzed, as well as a model of it that is compatible with OpenRocket simulation. This provides a centralized source of authoritative information on the design of Portland State Aerospace Society's 4th launch vehicle design, LV4.

This program and its dependencies were adapted by Cory from the work done by Ondrej, Erin, Ian on the original MDO.

This depends on three other programs, located in the same directory. They should be consulted, as they contain many details which are omitted here.

  • inputs
    • contains all library imports,
    • defines all constants and parameters,
    • sets initial conditions and constraints, and
    • provides a single location to specify all the parts of a rocket design which are not optimized.
  • trajectory
    • models a rocket based on inputs, theoretical rocket equations, and the design variables we optimizing.
    • simulates a trajectory by numerical integration of equations of motion
  • openrocket_interface
    • models an engine system for trajectory based on the design parameters and variables,
    • generates an OpenRocket model and its corresponding engine file for further testing, after conclusion of opitimazion.

End-User Requirements

  • A local copy of this repository
  • Python3
  • Jupyter Notebook
  • SciPy
  • NumPy
  • Matplotlib
  • OpenRocket
    • This depends on Java 8 (not later versions).
    • Technically unnecessary, but a file is created in the folder OpenRocket uses to store engine thrust curves.

Standard Operating Procedure

  1. RTFM.
  2. Review assumptions of supporting notebooks.
  3. Make appropriate modifications (if any) to parameters in inputs and set initial design variables there.
    • Note, garbage in, garbage out applies. If the initial design is infeasible, the optimization will abort.
  4. Run all the code blocks in this notebook in sequential order.
  5. Wait patiently for the algorithm to converge, perhaps 20 minutes depending on your computer.

Motivation

The equations of motion for a rocket, a variable-mass system, are highly non-linear and there are no closed-form analytic solutions. Additionally, launch vehicles are complex systems with multiple interactions and mutual dependencies between their subsystems, each of which in and of itself constitute a difficult engineering project.

This presents a difficulty for clean-sheet design of a launch vehicle to fulfill specific mission requirements. A sufficiently detailed model and numerical integration techniques can be used to circumvent mathematical challenges to evaluating a given rocket design, at the cost of increased computational resources. A well-posed optimization problem and an algorithm implentating it can be used to capture constraints and requirements, to allow a feasible and locally optimal solution to emerge naturally from system dynamics.

There are other approaches to rocket design, such as trade space analysis or design of experiments. We have chosen our approach because it allows us to optimize the whole system all at once, hopefully cutting down on the time and effort that would be spent keeping the system consistent while developing only chunks of it at a time.

Method

We may construct an objective function that represents the quality of our rocket as we see it. We selected Gross Liftoff Weight (GLOW) as our primary figure of merit, since it incorporates many aspects of a mission profile.

The Tsiolkovsky rocket equation $$\frac{\Delta v}{v_e} = \ln\left(\frac{m_0}{m_f}\right)$$ can be rearranged to isolate GLOW, $m_0$. The structure of the rocket and its subsystems determine its final dry mass, $m_f$; the effective exhaust velocity, $v_e$, captures engine performance; and the total change in velocity, $\delta v$, is determined by mission objectives, and can potentially include estimates of gravitational and atmospheric losses. So, minimizing GLOW while demanding a certain apogee is equivalent to simultaneously minimizing structural mass, maximizing engine performance, and balancing the conflicting goals of minimizing losses due to gravity and aerodynamics. Note that the sources of this conflict are the incentive to expel propellant rapidly to avoid the cost of carrying propellant in a gravitational field and the incentive to reduce velocity in lower atmosphere since drag is proportional to air pressure and the square of velocity.

In additional to this high-level objective, there are also several requirements and constraints that a given design must fulfill for a successful mission. These will be detailed in a following section. To represent them, we may use a series of barrier and penalty functions which bound potential rockets to a feasible region in the design space. A barrier function provides an absolute constraint that can not be violated under any circumstance, while a penalty function only disincentives designs if they stray far from a more permeable set of constraints.

By adding our three types of functions together, we obtain a merit function to minimize which encapsulates all of our relevant considerations as a positive real number. The input to our merit function is simply an array of values that are sufficient to determine our mathematical model of a rocket and its performance. However, each evaluation of the function entails that a simulation of the rocket's trajectory be computed, so the merit function itself does not admit to differentiation or finite differences.

There are several ways to approach this issue, two of which are covered in this program. If the design space is between 2 and 4 dimensions (inclusive), we may use the Nelder-Mead simplex method, a directed search that doesn't use derivatives. Note that a simplex is just a triangle generalized to n-dimensions. Performance of this geometry-based algorithm degrades as dimensionality increases, facilitating it becoming stuck in a local minimum. Another option is to use a genetic algorithm to seek a global optimum stochastically, which is much more computationally expensive but not restricted to low-dimensional spaces.

Theoretically, there are other methods possible, such as constructing and optimizing on a response surface model that is itself only a calculus-and-statistics-friendly model of our simulation model. We could also do Monte-Carlo simulation or some other stochastic optimization method, or even some sort of Latin Square trade-analysis. However, our results using both of the methods we selected closely agree and have been satisfactory.

Where we deviate significantly from the previous MDO is in our decision to do an iterative sequence of Nelder-Mead optimizations instead of just one. Before the penalty, and barrier functions are added to the objective function, they are weighed by user-selected parameters. Balancing these weights is an art; if they are too low, the optimizer can ignore or even abuse them, and if they are too high, the optimizer can neglect the objective and rapidly spiral into a somewhat artificial local minimum. Rather than picking arbitrary values that seem to give us the results we want, now we begin with very low weights, and for each successive optimization, we increase them.

This ought to enable more global coverage initially, while restricting the optimizer's freedom of movement in the design space after it has had a chance to find a satisfactory neighborhood. Since the weights become astronomically high, this also gives us a trivial method for determining whether we have articulated a well-posed problem and produced a strictly feasible design, or asked an unsatisfiable problem and received only the closest to being feasible design in a region of even worse alternatives.

Note that if a genetic algorithm is performed instead, an arbitrary weight must be selected intelligently, since a sequence of such optimization problems is exceedingly impractical.

We technically may not speak of "convergence", since it is possible for the Nelder-Mead method to stagnate on non-stationary points if the geometry of the simplex becomes ill-conditioned. However, since we are doing a sequence of optimizations, the simplex geometry is reset each time, so this is likely to not be an issue. The results have been in the same neighborhood as solutions given by a genetic algorithm, so our global coverage of the feasible region of the design space seems satisfactory.

Assumptions

  • If an optimization iteration returns a vector within $10^{-6}$ of the last iteration, there is no point in going further.
  • A simulation time-step of 0.25 s is good enough but not too good.
  • 25 iterations is a reasonable maximum number to perform.
  • Our initial design vector is already feasible (garbage in, garbage out).
  • A 12" diameter is a reasonable stipulation.
    • If this is much wider, then drag becomes horrible, if this is much thinner, then the L/D ratio becomes unreasonable. The standing assumption was already that LV4 would be 12", so we made this a fixed parameter instead of a design variable, in part for conceptual simplicity for members of our organization.
    • Our optimization algorithm also performs poorly in high-dimensional spaces, so this should protect us from the curse of dimensionality. It is trivial to make it a design variable again, but we are moving forward as is.
  • Our constraints are reasonable and satisfiable.
  • SciPy's implementations of the Nelder-Mead and genetic algorithm methods are sufficient for our purposes.

Optimization Conditions

The objective of the optimization is to minimize Gross Lift-Off Weight (GLOW). Overall project costs are likely to scale with GLOW.

We distilled the key design variables to just propellant mass, mass flow rate, and exit pressure. All other relevant variables are determined from those and our design constants and parameters. For the sake of brevity, we won't go into those details in the following discussion, but the interested reader can learn more from the supporting notebooks.

Our optimization routine currently includes 7 inequality constraints:

  • Apogee $h$: There are minimum and maximum acceptable altitudes of the trajectory.
  • Thrust $F$: The electric feed system that pressurizes propellant prior to injection is not feasible for engines that are too powerful.
  • Launch Speed $LS$: When the rocket leaves the launch rail it must be aerodynamically stable.
  • Thrust-to-Weight Ratio $TWR$: Trade-off between gravity loss and aerodynamic stability.
  • Length to Diameter Ratio $L/D$: Trade-off between aerodynamic stability and mechanical (non-rigid body) resonance modes.
  • Maximum acceleration $\frac{a_{max}}{g_0}$: Set by the material limit loads of various launch vehicle subsystems.
  • Nozzle over-expansion $\frac{p_e}{p_a}$: To ensure flow seperation.

Mathematically the problem can be stated as estimating $\bar{x}^*$ to a specified degree of precision according to these conditions,

$$ \lim_{n \to \infty} \mu_{n} = 0$$$$ \lim_{n \to \infty} \rho_{n} = \infty$$$$f_{n}(\bar{x}) = m_{obj}(\bar{x}) + \mu_{n}\sum_i h_i(\bar{x}) + \rho_{n}\sum_j g_j(\bar{x})$$$$ \bar{x}_n = \min_{\bar{x}} f_{n}(\bar{x})$$$$ \lim_{n \to \infty} \bar{x}_n = \bar{x}^{*},$$\begin{eqnarray*} \mbox{where} \hspace{5 mm} \begin{split}\bar{x} = \left\{ \begin{array}{ll} & m_{prop}\\ & \dot{m}\\ & p_e \end{array} \right. \end{split} \hspace{2 mm} \mbox{, subject to} \hspace{2 mm} \begin{split} h_{barrier} = \left\{ \begin{array} & 108401 m < h < 151401 \hspace{1 mm} m \end{array} \right. \end{split} \hspace{2 mm} \mbox{, and subject to} \hspace{2 mm} \begin{split} g_{penalty} = \left\{ \begin{array}{ll} & F \leq 6 \hspace{1 mm} kN\\ & LS \geq 22 \hspace{1 mm} m/s\\ & \frac{a_{max}}{g_0} \leq 15 \hspace{1 mm} g's\\ & TWR \geq 2\\ & L/D \leq 21\\ & \frac{p_e}{p_a} \geq 0.35\\ \end{array} \right. \end{split} \end{eqnarray*}

It is possible that for any given $n$ or in the limit, there is not a unique global minimum. However, a local minimum is attained which is very close to the results of a genetic algorithm.

By sending the limit of $\mu$ to 0, the barrier function is relaxed, which allows designs to approach the absolute window for acceptable apogees. By sending the limit of $\rho$ to $\infty$, the penalty function is strengthened, which prevents designs from leaving the feasible region.


In [1]:
# Nelder-Mead simplex search
%matplotlib inline
%run openrocket_interface.ipynb
%run System_Definition.ipynb
%run Simulation.ipynb

# i'm sorry for global vars...
global allvectors, dbz, allobjval
dbz = 0 # arithmetic error tracker (from crossing boundary constraints)
allvectors = []               # array for all design vecs, global variable
allobjfun = []                # array for tracking objective function evaluations

Functions of Merit

We chose to abstract all of the functions used within the merit function for increased flexibility, ease of reading, and later utility.

  • objective is arbitrarily constructed. It is
    • normalized (by a somewhat arbitrary constant) to bring it into the same range as constraints,
    • squared to reward (or punish) relative to the distance from nominal value, and
    • divided by two so that our nominal value is 0.5 instead of 1.0.
  • exact is less arbitrarily constructed. It is
    • squared to be horizontally symmetric (which could also be obtained by absolute value),
    • determined by the distance from a constant, and
    • divided by two is for aesthetic value.
  • exterior is not particularly arbitrary. It is
    • boolean so that we can specify whether it is minimizing or maximizing its variable,
    • 0 when the inequality is satisfied, otherwise it is just as punishing as exact.
  • barrier comes in two flavors, one of which is not used here. It is
    • boolean so that we can specify whether it is a lower or an upper bound,
    • completely inviolable, unlike exact and exterior penalties.

Technically logarithmic barrier functions allow negative penalties (i.e. rewards), but since we use upper and lower altitude barriers, it is impossible that their sum be less than 0. If the optimizer steps outside of the apogee window, the barrier functions can attempt undefined operations (specifically, taking the logarithm of a negative number), so some error handling is required to return an infinite value in those cases. Provided that the initial design is within the feasible region, the optimizer will not become disoriented by infinite values.


In [2]:
# all of our comparisons are ratios instead of subtractions because
# it's normalized, instead of dependent on magnitudes of variables and constraints

# minimize this, **2 makes it well behaved w.r.t. when var=cons
def objective(var, cons):
    return (var/cons)**2 / 2

# **2 because i like it more than abs(), but that also works
def exact(var, cons):
    return (var/cons - 1)**2 / 2

# this is your basic exterior penalty, either punishes for unfeasibility or is inactive
def exterior(var, cons, good_if_less_than=False):
    if good_if_less_than:
        return max(0, var/cons - 1)**2 / 2
    else:
        return max(0, -(var/cons - 1))**2 / 2

# this barrier function restricts our objective function to the strictly feasible region
# make rockets great again, build that wall, etc, watch out for undefined operations
def barrier(var, cons, int_point=False, good_if_less_than=True):
    global dbz
    try: # just in case we accidentally leave feasible region
        if not int_point:
            if good_if_less_than:
                return -log(-(var/cons - 1))
            else:
                return -log(var/cons - 1)
        elif int_point:
            def interior(g): return 1/g # in case we don't like logarithms, which is a mistake
            if good_if_less_than:
                return -interior(var/cons - 1)
            else:
                return -interior(-(var/cons - 1))
    except:
        dbz += 1 # keep track of arithmetic errors, side effect
        return float('inf') # ordinarily, this is bad practice since it could confuse the optimizer
                            # however, since this is a barrier function not an ordinary penalty, i think it's fine

Optimization Problem

Given a design vector $x$ and the iteration number $n$ our merit function f runs a trajectory simulation and evaluates the quality of that rocket. We keep track of each design and its merit value for later visualization, hence why global variables are used.

We run an iterative sequence of optimization routines for the win. We use the Euclidean distance in the design space between successive optimal designs to decide when it is no longer worth continuing, around the 6th decimal place.


In [3]:
# this manages all our constraints
# penalty parameters: mu -> 0 and rho -> infinity 
def penalty(ls, F, LD, TWR, S_crit, alt, max_g, stbl_margin, ballast, mu, rho):
    b = [barrier(alt, cons_alt, int_point=False, good_if_less_than=False),
         barrier(alt, cons_ceiling, int_point=False, good_if_less_than=True),
         #barrier(ballast, 0.0, int_point=False, good_if_less_than=False)
        ]
    eq = []
    ext = [exterior(F, cons_thrust, good_if_less_than=True),
           exterior(ls, cons_ls, good_if_less_than=False),
           exterior(LD, cons_LD, good_if_less_than=True),
           exterior(TWR, cons_TWR, good_if_less_than=False),
           exterior(S_crit, cons_S_crit, good_if_less_than=False),
           exterior(max_g, cons_accel, good_if_less_than=True),
           #exterior(stbl_margin, cons_stblty, good_if_less_than=False)
          ]
    return mu*sum(b) + rho*(sum(eq) + sum(ext))

# Pseudo-objective merit function
# x is array of design parameters, n is index of penalty and barrier functions
# print blocks are sanity checks so i'm not staring at a blank screen and can see what various tweaks actually do
def f(x, n=8):
    global allvectors, allobjfun
    # get trajectory data
    sim = trajectory(x[0], x[1], x[2],
               throttle_window, min_throttle,
               rcs_mdot, rcs_p_e, rcs_p_ch,
               x0[3], x0[4], x0[5], x0[6], x0[7],
               #ballast, root, tip, sweep, span,
               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,
               dt, True)
    
    obj_func = objective(sim.LV4.GLOW, cons_mass) # minimize GLOW
    # then, calculate penalization from trajectory performance
    pen_func = penalty(sim.launch_speed,
                  sim.thrust[0]/1000,
                  sim.LV4.length/sim.LV4.diameter,
                  sim.TWR,
                  sim.S_crit,
                  sim.alt[-1],
                  sim.max_g_force,
                  sim.min_stability,
                  sim.LV4.ballast,
                  mu_0 / (2**n),   # initial mu and
                  rho_0 * (2**n)) #                 rho are selected for nice behavior
    # add objective and penalty functions
    merit_func = obj_func + pen_func
    allvectors.append(x) # maintains a list of every design, side effect
    allobjfun.append(log(merit_func)) # log plot is cleaner to look at, for later
    return merit_func

# we want to iterate our optimizer for theoretical "convergence" reasons (given some assumptions)
# n = number of sequential iterations
def iterate(f, x_0, n):
    x = x_0[:3] # initial design vector, stripping off throttling variables for now
    global dbz
    designs = []
    for i in range(n):
        print("Iteration " + str(i+1) + ":")
        # this minimizer uses simplex method
        res = minimize(f, x, args=(i), method='nelder-mead', options={'disp': True, 'adaptive':True, 'xatol': 0.001, 'fatol': 0.001})
        x = res.x # feed optimal design vec into next iteration
        
        designs.append(res.x)   # we want to compare sequential objectives 
                                # so we can stop when convergence criteria met
        alt = trajectory(x[0], x[1], x[2],
               throttle_window, min_throttle,
               rcs_mdot, rcs_p_e, rcs_p_ch,
               x0[3], x0[4], x0[5], x0[6], x0[7],
               #ballast, root, tip, sweep, span,
               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,
               dt).alt
        
        print("         Arithmetic errors (from violations of acceptable altitude window): "+str(dbz))
        print("Propellant mass (kg): "+str(x[0]))
        #print("Airframe diameter (in): "+str(dia)) # in case we want diameter to be a design variable
        print("Mass flow rate (kg/s): "+str(x[1]))
        print("Exit pressure (Pa): "+str(x[2]))
        print("Peak Altitude (km): "+str(alt[-1]/1000))
        #print("Lower and Upper Drag Bounds (N): "+str(x[3])+', '+str(x[4]))
        print('')
        dbz=0 # I only care about divisions by zero in each individual iteration, side effect
        if (i > 0) and (np.linalg.norm(designs[-1] - designs[-2]) < delta):
            print("Early termination! The Euclidean distance between the last two designs was < " + str(delta))
            break
    return x

def breed_rockets(f):
    res = differential_evolution(f, [(80, 200), (2.0, 4.5), (45000, 150000)],
                                 strategy='best1bin', popsize=80, mutation=(.1, .8), recombination=.05,
                                 updating='immediate', disp=True, atol=0.1, tol=0.1)
                                 #polish=True,workers=-1, disp=True, atol=0.1, tol=0.1)
    return res.x

Optimization Information and Graphing

It is to the benefit of our intuition if we can visualize the design space and our final trajectory. This block of code simply provides all of the relevant information from a trajectory as text, and displays useful graphs.


In [4]:
# this creates some plots of the phase spaces of all our designs, doesn't save them
def phase_plot(m_prop, mdot_0, p_e):
    fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True)
    
    ax1.plot(m_prop, p_e)
    ax1.set_title('Design Vectors')
    ax1.yaxis.major.locator.set_params(nbins=6)
    ax1.set_xlabel("Propellant (kg)")
    ax1.set_ylabel("Exit pressure (kPa)")
    
    ax2.plot(mdot_0, p_e)
    ax2.yaxis.major.locator.set_params(nbins=6)
    ax2.set_xlabel("Mass flow rate (kg/s)")
    ax2.set_ylabel("Exit pressure (kPa)")
    
    ax3.plot(m_prop, mdot_0)
    ax3.yaxis.major.locator.set_params(nbins=6)
    ax3.set_ylabel("Mass flow rate (kg/s)")
    ax3.set_xlabel("Propellant (kg)")
    
    # we display the first diagram of projected 2d phase portraits
    plt.show()
    
    fig2 = plt.figure()
    ax = fig2.add_subplot(111, projection='3d')
    ax.set_title('Design Space Trajectory')
    ax.plot(m_prop, mdot_0, p_e)
    ax.set_xlabel("Propellant (kg)")
    ax.set_ylabel("Mass flow rate (kg/s)")
    ax.set_zlabel("Exit pressure (kPa)")
    
    # we display the interactive 3d phase portrait
    plt.show()
    # note, we're choosing to not automatically save these, but they can be saved from the interface
    fig3 = plt.figure()
    ax = fig3.add_subplot(111)
    ax.set_title('Design Values as % of optimum design value')
    ax.plot(m_prop/(m_prop[-1]), "r", label="m_prop")
    ax.plot(mdot_0/(mdot_0[-1]), "g", label="mdot_0")
    ax.plot(p_e/(p_e[-1]), "b", label="p_e")
    ax.legend()
    plt.show()

def design_grapher(): # no arguments because I'm using dirty dirty global variables
    # get us some nice plots of the phase space of design vectors
    designplot = [[],[],[]]
    for i in range(0, len(allvectors)):
        for j in range(0, len(designplot)):
            designplot[j].append(allvectors[i][j])
    phase_plot(designplot[0], designplot[1], designplot[2])

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_title('Objective + Constraint convergence')
    ax.plot(allobjfun, "r", label="Objective Function Evaluations (Logarithmic Scale)")
    ax.legend()
    plt.show()

Top-Level of Optimization Routine

Here's where the magic happens. This code block runs the iterative optimization, provides details from our optimized trajectory, uses the OpenRocket interface to make a model rocket and engine for higher-fidelity analysis, and then displays visuals.


In [5]:
# Results, this is the big boi function
if __name__ == '__main__':
    test = trajectory(x0[0], x0[1], x0[2],
               throttle_window, min_throttle,
               rcs_mdot, rcs_p_e, rcs_p_ch,
               x0[3], x0[4], x0[5], x0[6], x0[7],   
               #ballast, root, tip, sweep, span,
               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,
               dt)
    
    if (test.alt[-1] < cons_alt) or (test.alt[-1] > cons_ceiling): # rudimentary error handling to save heartache
        raise Exception('Rocket apogee out of bounds! Apogee {:.3f} km'.format(test.alt[-1]/1000))
        
    # feed initial design into iterative optimizer, get most (locally) feasible design
    res = iterate(f, x0, iterations)
    # probe design space, darwin style. if design space has more than 3 dimensions, you need this. takes forever.
    #res = breed_rockets(f)
    
    print("Optimization done!")
    
    # Rename the optimized output for convenience
    m_prop    = res[0]
    mdot = res[1]
    p_e  = res[2]
    
    # get trajectory info from optimal design
    sim = trajectory(m_prop, mdot, p_e,
               throttle_window, min_throttle,
               rcs_mdot, rcs_p_e, rcs_p_ch,
               x0[3], x0[4], x0[5], x0[6], x0[7],
               #ballast, root, tip, sweep, span,
               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,
               0.05)
    
    textlist = print_results(sim, True)
    # draw pretty pictures of optimized trajectory
    rocket_plot(sim.t, sim.alt, sim.v, sim.a, sim.thrust,
                sim.dyn_press, sim.Ma, sim.m, sim.p_a, sim.drag, sim.throttle, True)
    
    # get/print info about our trajectory and rocket
    for line in textlist:
        print(line)
    
    print('\nMaking an OpenRocket rocket and corresponding engine!')
    # create an openrocket file with matching engine for our design (and print/save trajectory data)
    make_engine(mdot, sim.m_prop, sim.thrust[0:sim.F_index + 1],
                sim.LV4.inr_r, airframe_thickness, sim.LV4.l_o, sim.LV4.l_f, sim.LV4.m_tank_o, sim.LV4.m_tank_f,
                sim.t[sim.F_index], sim.LV4.engine.Ve/g_n,
                sim.LV4.eng_sys_dry_mass, sim.LV4.eng_sys_len, sim.openrocket_CoM,
                sim.LV4.ballast, sim.LV4.fin)
    
    # draw more pretty pictures, but of the optimizer guts
    design_grapher()


Iteration 1:
Optimization terminated successfully.
         Current function value: 0.587137
         Iterations: 209
         Function evaluations: 449
         Arithmetic errors (from violations of acceptable altitude window): 13
Propellant mass (kg): 144.94101497075505
Mass flow rate (kg/s): 2.8788268521459215
Exit pressure (Pa): 73774.67204148213
Peak Altitude (km): 111.78216340369084

Iteration 2:
Optimization terminated successfully.
         Current function value: 0.685156
         Iterations: 198
         Function evaluations: 445
         Arithmetic errors (from violations of acceptable altitude window): 177
Propellant mass (kg): 144.73267806591193
Mass flow rate (kg/s): 2.8711639874533255
Exit pressure (Pa): 74610.00364614137
Peak Altitude (km): 111.30373231073928

Iteration 3:
Optimization terminated successfully.
         Current function value: 0.894289
         Iterations: 61
         Function evaluations: 208
         Arithmetic errors (from violations of acceptable altitude window): 140
Propellant mass (kg): 144.73267806591193
Mass flow rate (kg/s): 2.8711639874533255
Exit pressure (Pa): 74610.00364614137
Peak Altitude (km): 111.30373231073928

Early termination! The Euclidean distance between the last two designs was < 0.001
Optimization done!
DESIGN VECTOR

-----------------------------

design total propellant mass               = 144.589 kg

design mass flow rate                      = 2.871 kg/s

design nozzle exit pressure                = 74610.004 kPa

total tankage length (after adjustment)    = 2.198 m

design airframe diameter                   = 0.321 in.

design airframe total length               = 6.233 m.

design GLOW                                = 259.340 kg

design ballast mass                        = 0.000 kg

design fin root chord                      = 0.700 m

design fin tip chord                       = 0.450 m

design fin sweep angle                     = 60.000 deg

design fin span                            = 0.420 m



CONSTRAINTS

-----------------------------

L/D ratio (c.f. < 20.0)                       = 19.399

Sommerfield criterion (c.f. pe/pa >= 0.35)    = 0.872

max acceleration (c.f. < 15.0)                = 5.639 gs

TWR at lift off (c.f. > 2.0)                 = 1.727

Lowest stability margin caliber (c.f. > 2.0) = 1.588

speed when leaving launch rail (c.f. > 22.0)  = 25.711 m/s

altitude at apogee (c.f. > 111.401)              = 111.564 km

design thrust (ground level) (c.f. < 6.0)    = 6.932 kN



ADDITIONAL INFORMATION

-----------------------------

design thrust (vacuum)                     = 7.89 kN

design total dry mass                      = 114.776 kg

mission time at apogee                     = 179.250 s

mission time at burnout                    = 50.300 s

max dynamic pressure                       = 70.033 kPa

design dV                                  = 2.004 km/s

estimated minimum required dV              = 1.479 km/s



ENGINE SYSTEM DETAILS

-----------------------------

Ox flow:                                   = 1.623 kg/s

Fuel flow:                                 = 1.248 kg/s

Ox mass:                                   = 81.805 kg

Fuel mass:                                 = 62.927 kg

Ox tank length + ullage:                   = 1.081 m

Fuel tank length + ullage:                 = 1.117 m

design chamber pressure                    = 2413.166 kPa

design expansion ratio                     = 5.798

design Exit area                           = 17.783 in.^2

design throat area                         = 3.067 in.^2

design Throat pressure                     = 1398.335 kPa

design Throat temperature                  = 2915.458 K

design Chamber temperature                 = 3097.820 K

design exit velocity                       = 2458.196 m/s

design isp                                 = 250.666 s

design average impulse                     = 372.876 kN*s



POST-FLIGHT MASS BUDGET

-----------------------------

LV4	GLOW:259.3395524909266	Current Mass:114.77644572264897

0.	Fins				Mass: 8.98

  0.		Front					Mass: 2.25

  1.		Back					Mass: 2.25

  2.		Left					Mass: 2.25

  3.		Right					Mass: 2.25

1.	Engine Subsystem		Mass: 56.15

  0.		Engine Module				Mass: 22.72

    0.			Engine Module					Mass: 0.81

    1.			Engine Module					Mass: 0.26

    2.			Engine Module					Mass: 1.7

    3.			Ring&Clamp					Mass: 1.49

    4.			Bulkhead					Mass: 0.63

    5.			Bulkhead					Mass: 0.63

    6.			Engine						Mass: 6.12

    7.			EMS						Mass: 1.0

    8.			EFS						Mass: 10.07

  1.		LOX Tank				Mass: 17.41

    0.			Tank Walls					Mass: 11.74

    1.			Insulation					Mass: 1.75

    2.			Bottom Cap					Mass: 0.52

    3.			Top Cap						Mass: 0.52

    4.			Ring&Clamp					Mass: 1.52

    5.			Bulkhead					Mass: 0.63

    6.			Bulkhead					Mass: 0.63

    7.			Propellant					Mass: 0.1

  2.		IPA/H20 Tank				Mass: 16.02

    0.			Tank Walls					Mass: 12.14

    1.			Bottom Cap					Mass: 0.52

    2.			Top Cap						Mass: 0.52

    3.			Ring&Clamp					Mass: 1.52

    4.			Bulkhead					Mass: 0.63

    5.			Bulkhead					Mass: 0.63

    6.			Propellant					Mass: 0.07

2.	Passthru Module			Mass: 6.22

  0.		Passthru Module				Mass: 0.26

  1.		Passthru Module				Mass: 2.66

  2.		Passthru Module				Mass: 0.54

  3.		Ring&Clamp				Mass: 1.49

  4.		Bulkhead				Mass: 0.63

  5.		Bulkhead				Mass: 0.63

3.	AV Module			Mass: 8.36

  0.		AV Module				Mass: 0.69

  1.		AV Module				Mass: 0.18

  2.		AV Module				Mass: 1.44

  3.		Ring&Clamp				Mass: 1.49

  4.		Bulkhead				Mass: 0.63

  5.		Bulkhead				Mass: 0.63

  6.		AV/360					Mass: 3.3

4.	N2 Module			Mass: 12.74

  0.		N2 Module				Mass: 0.58

  1.		N2 Module				Mass: 0.18

  2.		N2 Module				Mass: 1.21

  3.		Ring&Clamp				Mass: 1.49

  4.		Bulkhead				Mass: 0.63

  5.		Bulkhead				Mass: 0.63

  6.		RCS					Mass: 8.02

    0.			N2 Tank						Mass: 4.0

      0.				Tank Walls						Mass: 1.79

      1.				Bottom Cap						Mass: 0.15

      2.				Top Cap							Mass: 0.15

      3.				Propellant						Mass: 1.9

    1.			Gas Jet						Mass: 4.02

5.	RCS Module			Mass: 6.22

  0.		RCS Module				Mass: 0.26

  1.		RCS Module				Mass: 2.66

  2.		RCS Module				Mass: 0.54

  3.		Ring&Clamp				Mass: 1.49

  4.		Bulkhead				Mass: 0.63

  5.		Bulkhead				Mass: 0.63

6.	ERS Module			Mass: 9.74

  0.		ERS Module				Mass: 0.26

  1.		ERS Module				Mass: 2.66

  2.		ERS Module				Mass: 0.54

  3.		Ring&Clamp				Mass: 1.49

  4.		Bulkhead				Mass: 0.63

  5.		Bulkhead				Mass: 0.63

  6.		ERS					Mass: 3.52

7.	Nosecone			Mass: 6.38

  0.		Cylinder				Mass: 4.33

    0.			Cylinder					Mass: 0.65

    1.			Cylinder					Mass: 0.21

    2.			Cylinder					Mass: 1.36

    3.			Ring&Clamp					Mass: 1.49

    4.			Bulkhead					Mass: 0.63

  1.		Cone					Mass: 0.95

    0.			Cone						Mass: 0.28

    1.			Cone						Mass: 0.08

    2.			Cone						Mass: 0.58

  2.		Parachutes				Mass: 1.1

  3.		Ballast					Mass: 0.0


Making an OpenRocket rocket and corresponding engine!
New rocket generated from template!
Reopen OpenRocket to run simulation.

In [6]:
'''m_prop    = 141.88397561701234
mdot = 2.8728568968933215
p_e  = 71035.22181877575

    
    # get trajectory info from optimal design
sim = trajectory(m_prop, mdot, p_e,
               throttle_window, min_throttle,
               rcs_mdot, rcs_p_e, rcs_p_ch,
               root, tip, sweep, span,
               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,
               0.05)
sim.LV4.read_out()
    
    # get/print info about our trajectory and rocket
    #res_text = print_results(res)
    #for line in res_text:
    #    print(line)
res_text = ""
frame_thickness = 0.008255
inr_radius = sim.LV4.inr_r
print('Engine system details in trajectory log!')
print('\nMaking an OpenRocket rocket and corresponding engine!')
    
    # create an openrocket file with matching engine for our design (and print/save trajectory data)
make_engine(mdot, sim.m_prop, sim.F[0:sim.F_index + 1],
                inr_radius, frame_thickness, sim.LV4.l_o, sim.LV4.l_f, sim.LV4.m_tank_o, sim.LV4.m_tank_f,
                sim.t[sim.F_index], sim.LV4.engine.Ve/g_n, res_text,
                sim.LV4.eng_sys_dry_mass, sim.LV4.eng_sys_len, sim.openrocket_CoM)
    
    # draw pretty pictures of optimized trajectory
#rocket_plot(sim.t, sim.alt, sim.v, sim.a, sim.F, sim.q, sim.Ma, sim.m, sim.p_a, sim.D, sim.throttle)'''


Out[6]:
'm_prop    = 141.88397561701234\nmdot = 2.8728568968933215\np_e  = 71035.22181877575\n\n    \n    # get trajectory info from optimal design\nsim = trajectory(m_prop, mdot, p_e,\n               throttle_window, min_throttle,\n               rcs_mdot, rcs_p_e, rcs_p_ch,\n               root, tip, sweep, span,\n               airfrm_in_rad, OF, p_ch, T_ch, ke, MM,\n               0.05)\nsim.LV4.read_out()\n    \n    # get/print info about our trajectory and rocket\n    #res_text = print_results(res)\n    #for line in res_text:\n    #    print(line)\nres_text = ""\nframe_thickness = 0.008255\ninr_radius = sim.LV4.inr_r\nprint(\'Engine system details in trajectory log!\')\nprint(\'\nMaking an OpenRocket rocket and corresponding engine!\')\n    \n    # create an openrocket file with matching engine for our design (and print/save trajectory data)\nmake_engine(mdot, sim.m_prop, sim.F[0:sim.F_index + 1],\n                inr_radius, frame_thickness, sim.LV4.l_o, sim.LV4.l_f, sim.LV4.m_tank_o, sim.LV4.m_tank_f,\n                sim.t[sim.F_index], sim.LV4.engine.Ve/g_n, res_text,\n                sim.LV4.eng_sys_dry_mass, sim.LV4.eng_sys_len, sim.openrocket_CoM)\n    \n    # draw pretty pictures of optimized trajectory\n#rocket_plot(sim.t, sim.alt, sim.v, sim.a, sim.F, sim.q, sim.Ma, sim.m, sim.p_a, sim.D, sim.throttle)'

Above, we can see the details of each iteration of the optimization sequence.

Then, we have detailed information from the trajectory of our optimal rocket design. Note that there are some low-level details suppressed here, but contained in the saved text file.

There is a plot of that trajectory and some of the state variables we care about, which can be enlightening.

Then there are three projections of the optimization path in the design space, as well as its 3D representation. This informs a sense of how an initial design is improved by the optimizer, and how design variables relate to each other.

The next visualization is of our design variables normalized as a percentage of their optimal values. This reduces the entire optimization into one easily comprehendable image. Note that spikes occur where iterations begin, since the Nelder-Mead method is reinitialized and must relearn how far it can move in any direction.

Finally, there is a semi-log plot of the merit values of our sequence of designs. The spikes are again from the beginning of iterations. If our problem is well-posed, punishment for deviation from feasibility will exponentially increase with each iteration and the minimum of each iteration will be decreasing. An ill-conditioned problem will have an increasing minimum with each iteration, since there would be no feasible design and the punishments would be growing. In either case, gaps appear in the plot where a design failed to satisfy the apogee window.

Conclusion and Next Steps

Prima facie, this program produces viable rocket designs. The next step is to experiment with in OpenRocket with the model produced. It is critical that we make sure the stability requirement set by Base-11 is met (calibers > 2 throughout flight), and it is also an opportunity to experiment with different masses and fins, as well as environmental conditions such as wind or launching at 85 degrees.

Our sequence of merit evaluations is bounded below by 0 but is not bounded above, since some designs have infinite values. Provided that the problem is well-conditioned, the subsequence of minimums from each iteration is monotonically decreasing. So, it seems reasonable to infer that an ideal optimal rocket design exists for a well-posed problem. This program approximates that ideal design, but whether it is a good approximation might be a qualitative judgement.

There are other design and optimization methods, and different choices that could have been made. This code is fairly flexible so it isn't too difficult to change assumptions on the fly. The goal is not to provide a definitive final design of a rocket; rather it is to provide a design tool which will be useful throughout the design, manufacturing, and testing processes of LV4, up to and including the day it launches.

It is important to remember that the difference between design parameters and variables is largely stipulative. With slight effort, some parameters could be freed for optimization and some variables could be frozen. The chamber pressure has been stipulated and determines the thermodynamic properties that made us select the oxygen/fuel ratio we stipulate. There is another script in this repository that can be used to perform a trade-space analysis of chamber pressure and propellant configurations, if we decide to revisit those preliminary assumptions.

We have been conservative in our mass estimates in order to perform a more robust optimization. We could also take the maximum merit evaluation in the neighborborhood of a given design, as bounded by engineering tolerances, instead of simply evaluating a single design. That would obviously be computationally costly, but each neighborhood could be evaluated by parallel distributed processes, and the benefit of implementing a minmax problem instead of a minimization problem could be worth it, considering the inevitable inefficiencies that may arise in a physical system.

It would be a good idea to do some statistical analysis of the design space to produce a heat map as well, so we could characterize the topology of our design space. It's possible that there may be multiple optimal neighborhoods, and this would answer that question. Note that a sensitivity analysis of the current optimization program could also answer that question, by probing the design space with different initial designs.