{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Design Optimization of the LV4 Sub-orbital Rocket\n",
"## Motivation\n",
"Clean-sheet conceptual design and trajectory optimization of launch vehicles is a classically difficult problem. This problem arises for two reasons. The first is that the 1-D trajectory equation is a 2nd order non-linear ordinary differential equation with non-constant coefficients, and has no closed-form solution. Secondly, detailed design choices in propulsion, structures/weights, aerodynamics, and guidance and control which ultimately all appear as variables in the governing equation are both highly coupled and are non-hierarchical. Schematically the coupling between the variables is presented in the figure below. The traditional approach to dealing with this problem has been to evaluate vehicle performance by using \"one variable-at-a-time\" trade-off analysis. However there are several key limitations to this approach:\n",
"\n",
"* The nature of the coupling between the design variables is poorly understood\n",
"* Early conceptual design is usually carried out with low-fidelity models\n",
"* Optimizing individual design variables does not usually guarantee optimality at the overall system level\n",
"\n",
"In practice this results in a highly iterative design process and concomitant requirement mismatches, developmental dead-ends, and sub-optimal final design. The problem with sub-optimal design is magnified by the ramifications of Konstantin Tsiolkovsy's equation, with its severe implication of exponential growth in design requirements for linear increases in rocket dry mass. Since ultimately all design decisions impact the rocket dry mass in some way it is imperative to understand these trade-offs and compromises as early in the design process as possible to reduce the potential for technical, schedule, and cost risks. This is especially the case for time, technical expertise, and funding constrained student organizations. Thus there is a strong motivation to treat the conceptual design parameters \"all-at-once\" using applied convex optimization techniques. A numerical design optimization approach allows us to start to explore this extremely large trade space and point us towards a good \"initial guess\" for the iterative design process."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import SVG\n",
"SVG(filename='coupling.svg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some commercial and governmental tools exist for launch vehicle design optimization. These include\n",
"codes such as FASTPASS, and SWORD. There are also open-source design and optimization tools that\n",
"can be applied to high-powered rocketry such as Open Rocket, and JSBSim, however these tools either\n",
"cannot be run in a batch mode or lack I/O tools to support direct numerical optimization. In the context\n",
"of this developing interest in clean-sheet small launch vehicle design we can identify a need for a design tool\n",
"high-level enough for simplicity, speed, and ease of use, but which captures enough of the dimensions of the\n",
"optimization problem to still be useful as a guide in the early conceptual design phase."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Methods\n",
"We expect overall project costs to scale with the initial mass of the launch vehicle (including the mass of loaded propellants), which is usually defined by the Gross Lift-Off Weight (GLOW). Therefore the objective of the optimization is to minimize GLOW. \n",
"\n",
"For a top-level early concept development we can distill the key design variables to just thrust, nozzle expansion ratio, airframe diameter and propellant tank lengths. There are also a large number of design constants present in the models, which include specific impulse and combustion chamber pressure, propellant mixture ratio, and others, but which for the sake of brevity will be largely ignored in the following discussion. The model currently includes 6 inequality constraints:\n",
"\n",
"* Apogee $h$: The maximum altitude on the trajectory\n",
"* Thrust-to-Weight Ratio $TWR$: Trade-off between gravity loss and aerodynamic stability\n",
"* Length to Diameter Ratio $L/D$: Trade-off between aerodynamic stability and mechanical (non-rigid body) resonance modes\n",
"* Maximum acceleration $\\frac{a_{max}}{g_0}$: Set by the material limits of various launch vehicle subsystems\n",
"* Diameter $D$: Manufacturability (current process scale-able from \\diameter6\" to up to \\diameter14\") \n",
"* Nozzle over-expansion $\\frac{p_e}{p_a}$: Prevents the assumption made to help linearize thrust model from becoming invalid\n",
"\n",
"Mathematically the problem can be stated as\n",
"\n",
"$$ \\mbox{min} \\hspace{2 mm} m_{initial} = m_{dry}(h) + m_{propellants, initial}(h) $$\n",
"\n",
"\n",
"\\begin{eqnarray*} \\mbox{where} \\hspace{5 mm} \\begin{split}x = \\left\\{ \\begin{array}{ll} & \\dot{m}\\\\\n",
" & D\\\\\n",
" & L\\\\\n",
" & p_e \\end{array} \\right. \n",
" \\end{split} \\hspace{2 mm} \\mbox{subject to} \\hspace{2 mm} \\begin{split}\n",
" g = \\left\\{ \\begin{array}{ll}\n",
" & 5 \\leq L/D \\leq 15\\\\\n",
" & TWR \\geq 2\\\\\n",
" & \\frac{p_e}{p_a} \\geq 0.5\\\\\n",
" & 6 \\leq D \\leq 14\\\\\n",
" & h \\geq 100000\\\\\n",
" & \\frac{a_{max}}{g_0} \\leq 15 \\end{array} \\right. \\end{split}\\end{eqnarray*}\n",
"\n",
"While the objective function itself is given by a very simple function of the design variables $D$ and $L$, the constraints require a much more complex model. The following section discusses the trajectory simulation model required to capture values of these constraints for any particular design vector.\n",
"\n",
"### Trajectory model\n",
"The trajectory of the rocket, in a single degree of freedom, can be best described by altitude $h(t)$, velocity $V(t)$, and total mass $m(t)$ state variables which are functions of time. The initial values of the state variables are given by $h_0$, $V_0$, and $m_0$ respectively. The governing equation of the rocket trajectory is Newton's 2nd law of motion \n",
"\n",
"$$F = \\frac{d(mV)}{dt}.$$\n",
"\n",
"This can be expanded in terms of the sum of forces acting on the rocket during free flight going up to apogee\n",
"\n",
"$$Thrust - Drag - m g_0 = \\frac{d}{dt} \\left(m \\frac{dh}{dt}\\right).$$\n",
"\n",
"We usually express drag as\n",
"$$F_d = \\frac{1}{2} \\rho V^2 C_d A$$\n",
"\n",
"where $\\rho$ is the local air density, $A$ is the frontal area, $C_d$ is the drag coefficient which for a given airframe geometry is a function of angle-of-attack (which for a 1-DoF system is identically zero), and Mach number. The Mach number is a function given by\n",
"\n",
"$$ \\frac{V}{c} = \\frac{V}{\\sqrt{kRT}} $$\n",
"\n",
"where $c$ is the local speed of sound, $k$ is the gas ratio of specific heats, $R$ is the specific gas constant (assumed constant in altitude), and $T$ is the gas temperature. These numbers, as well as $\\rho$, and $p_a$ the ambient pressure, are supplied by the 1976 U.S. Standard Atmosphere model. For simplicity, drag coefficients are interpolated from published aerodynamic data from the 1950's era NACA/USN/NASA Aerobee-150 sounding rocket\\cite{aerobee}, which is expected to be dimensionally similar (though not identical to) to the LV4 sounding rocket.\n",
"\n",
"The rocket thrust force, given a huge number of simplifying assumptions, is a function of the mass flow rate design variable, the chamber pressure constant, and the ambient pressure state variable. It is given by\n",
"$$F_t = \\dot{m}V_e + A_e(p_e - p_a) $$\n",
"\n",
"where $V_e$ is the rocket exhaust velocity, $A_e$ is the rocket nozzle exit area. The exhaust velocity, again given a large number of simplifying assumptions, is given by\n",
"\n",
"$$V_e = \\sqrt{\\frac{2kRT_{ch}}{k - 1}\\left(1 - \\left(\\frac{p_e}{p_{ch}}\\right)^{((k - 1)/k)}\\right)}$$\n",
"\n",
"where $T_{ch}$ is the chamber pressure, which along with the combustion chamber gas ratio of specific heats $k$, and gas constant $R$ are functions of the mixture ratio of the propellants, and the chamber pressure $p_{ch}$ determined using the NASA Chemical Equilibrium Analysis (CEA)\\cite{CEA} tool and are design constants. It should be noted that the $p_e$ is a design variable tied to the expansion ratio of the rocket engine, but that $p_a$ is constantly changing with altitude. The thrust function is maximized when $p_e = p_a$, however this occurs at one, and only one, altitude. Thus the design expansion ratio selection directly benefits from trajectory optimization. Some additional information about PSAS's liquid propulsion project can be found [here](https://github.com/psas/liquid-engine-capstone-2015).\n",
"\n",
"By finite-differencing the derivatives of the governing equation a first-order Runge-Kutta (e.g. Forward Euler) integration of the governing equation is used to determine the rocket trajectory. This approach was chosen only for simplicity; in principle central differencing or trapezoidal integration are more generally accurate discretization approaches. Forward Euler difference equations always have the form\n",
"\n",
"$$y_{i+1} = y_i + \\left(y\\hspace{-1 mm}-\\hspace{-1 mm}\\mbox{rate at } t_i \\right)(t_{i+1} - t_i) $$\n",
"\n",
"where $y(t)$ is the state variable being integrated, and $i$ is the time index superscript. In this way the state variables $h$, $V$, and $m$ can be defined at any discrete time $t_0$, $t_1$, $t_2$ ... $t_i$. Given a vector of the four design variables the trajectory function returns numerical values of the objective function and constraint vector. This trajectory code was benchmarked against known trajectory and design data for the Aerobee-150 (sustainer only), and Armadillo Stig-B sounding rockets, and the results compare favorably. This grants us at least some confidence in the assumptions made in the model."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import csv\n",
"import numpy as np\n",
" \n",
"# Drag coefficient look up\n",
"C_d_t = []\n",
"Ma_t = []\n",
"f = open('CD_sustainer_poweron.csv') # Use aerobee 150 drag data\n",
"aerobee_cd_data = csv.reader(f, delimiter=',')\n",
"for row in aerobee_cd_data:\n",
" C_d_t.append(row[1])\n",
" Ma_t.append(row[0])\n",
"\n",
"# Plot the results\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import rc\n",
"%config InlineBackend.figure_formats=['svg']\n",
"%matplotlib inline\n",
"rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})\n",
"rc('text', usetex=True)\n",
"\n",
"Ma_t = np.array(Ma_t, dtype=float)\n",
"C_d_t = np.array(C_d_t, dtype=float)\n",
"plt.plot(Ma_t, C_d_t)\n",
"plt.xlabel('Mach number')\n",
"plt.ylabel('Cd')\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Airframe Section Parameters\n",
"\n",
"These are parameters from the current 18\" carbon composite airframe modules used in the dry mass model. Further information about this project can be found [here](https://github.com/psas/lv3.0-airframe)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Module planar density = 1.92 kg/m^2\n"
]
}
],
"source": [
"# Airframe Section Parameters\n",
"from math import pi\n",
"m = 0.463 # Single LV3 airframe module mass [kg], 18 in.\n",
"Do = 0.168 # Outer diameter [m]\n",
"L = 0.457 # Length [m]\n",
"A_fiber = pi * Do *L\n",
"rho_fiber = m/A_fiber # Planar density [kg/m^2]\n",
"print(\"Module planar density = {0:.2f} kg/m^2\".format(rho_fiber))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Trjaectory Simulation"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true,
"scrolled": true
},
"outputs": [],
"source": [
"from math import sqrt, pi, exp, log, cos\n",
"import numpy as np\n",
"\n",
"\n",
"# A simple forward Euler integration for rocket trajectories\n",
"def dry_mass(L, dia):\n",
" m_avionics = 3.3 # Avionics mass [kg]\n",
" m_recovery = 4 # Recovery system mass [kg]\n",
" m_payload = 1 + 3 # Payload mass [kg]\n",
" m_tank_fudgefactor = 1.5\n",
" m_tankage = rho_fiber * m_tank_fudgefactor * L * dia * pi + 3 # Tank mass Estimation [kg]\n",
" m_engine = 3 # Engine mass [kg]\n",
" m_feedsys = 10 # Feed system mass [kg]\n",
" len_airframe = 1 # Length of Airframe (sans tankage) [m]\n",
" m_airframe = rho_fiber * m_tank_fudgefactor * len_airframe * dia * pi + 1 # Airframe mass [kg]\n",
" return (m_avionics + m_recovery + m_payload + m_tankage \n",
" + m_engine + m_feedsys + m_airframe) # Total dry mass [kg]\n",
" \n",
"def propellant_mass(A, L, OF=1.3):\n",
" rho_alc = 852.3 # Density, ethanol fuel [kg/m^3]\n",
" rho_lox = 1141.0 # Density, lox [kg/m^3]\n",
" L_lox = L/(rho_lox/(rho_alc*OF) + 1)\n",
" m_lox = rho_lox*L_lox*A # Oxidizer mass [kg]\n",
" m_alc = rho_alc*(L-L_lox)*A # Fuel mass [kg]\n",
" return m_alc + m_lox # Propellant Mass [kg]\n",
"\n",
"def std_at(h): # U.S. 1976 Standard Atmosphere\n",
" if h < 11000:\n",
" T = 15.04 - 0.00649*h\n",
" p = 101.29*((T + 273.1)/288.08)**5.256\n",
"\n",
" elif 11000 <= h and h <25000:\n",
" T = -56.46\n",
" p = 22.65*exp(1.73 - 0.000157*h)\n",
"\n",
" else:\n",
" T = -131.21 + 0.00299*h\n",
" p = 2.488 * ((T + 273.1)/216.6)**(-11.388)\n",
"\n",
" rho = p/(0.2869*(T + 273.1)) # Ambient air density [kg/m^3]\n",
" p_a = p*1000 # Ambient air pressure [Pa]\n",
" T_a = T + 273.1 # Ambient air temperature [K]\n",
" return p_a, rho, T_a\n",
"\n",
"def thrust(x, p_ch, T_ch, p_e, ke, Re, mdot):\n",
" p_a = std_at(x)[0] # Ambient air pressure [Pa]\n",
" p_t = p_ch * (1 + (ke - 1)/2)**(-ke/(ke - 1)) # Throat pressure [Pa]\n",
" T_t = T_ch * (1/(1 + (ke - 1)/2)) # Throat temperature [K]\n",
" A_t = (mdot / p_t) * sqrt(Re * T_t/ke) # Throat area [m^2]\n",
" A_e = A_t * (2/(ke + 1))**(1/(ke - 1)) * (p_ch/p_e)**(1/ke) * \\\n",
" 1/sqrt((ke + 1)/(ke - 1) * (1 - (p_e/p_ch)**((ke - 1)/ke))) # Exit area [m^2]\n",
" \n",
" ex = A_e/A_t # Expansion ratio\n",
" #alpha_t = [14, 11, 10, 9] # Lookup table of divergence angles, assuming 80% bell length\n",
" #ex_t = [5, 10, 15, 20] # Lookup table of expansion ratios from alpha_t\n",
" #alpha= np.interp(ex, ex_t, alpha_t)\n",
" #lam = 0.5 * (1 + cos(alpha * pi/180)) # Thrust cosine loss correction (not an O(1) effect, even in extreme cases) \n",
" lam = 1\n",
" Ve = lam * sqrt(2 * ke/(ke - 1) * Re * T_ch * \\\n",
" (1 - (p_e/p_ch)**((ke - 1)/ke))) * 0.9 # Exhaust velocity [m/s]\n",
" \n",
" F = mdot * Ve + (p_e - p_a) * A_e # Thrust force, ignoring that isp increases w/ p_ch [N]\n",
" return F, A_t, A_e, Ve\n",
"\n",
"def drag(x, v, A, Ma, C_d_t, Ma_t):\n",
" (p_a, rho, T_a) = std_at(x)\n",
" #C_d_t = [0.15, 0.15, 0.3, 0.45, 0.25, 0.2, 0.175, .15, .15] # V2 rocket drag coefficient lookup table\n",
" #Ma_t = [0, 0.6, 1.0, 1.1, 2, 3, 4, 5, 5.6] # V2 rocket Mach number lookup table\n",
" C_d = np.interp(Ma, Ma_t, C_d_t) # Drag coefficient function\n",
" q = 0.5 * rho * v**2 # Dyanmic pressure [Pa]\n",
" D = q * C_d * A # Drag force [N]\n",
" return D, q\n",
"\n",
"def trajectory(L, mdot, dia, p_e, p_ch=350, T_ch=3500, ke=1.3, Re=349, x_init=0):\n",
" # Note combustion gas properties ke, Re, T_ch, etc, determined from CEA\n",
" # Physical constants\n",
" g_0 = 9.81 # Gravitational acceleration [m/s^2]\n",
" dt = 1 # Time step [s]\n",
" ka = 1.4 # Ratio of specific heats, air \n",
" Ra = 287.1 # Avg. specific gas constant (dry air)\n",
" \n",
" # LV4 design variables\n",
" dia = dia*0.0254 # Convert in. to m\n",
" A = pi*(dia/2)**2 # Airframe frontal area projected onto a circle of diameter variable dia\n",
" m_dry = dry_mass(L, A) # Dry mass, call from function dry_mass()\n",
" mdot = mdot # Mass flow rate [kg/s]\n",
" p_ch = p_ch*6894.76 # Chamber pressure, convert psi to Pa\n",
" p_e = p_e*1000 # Exit pressure, convert kPa to Pa\n",
"\n",
" # Initial conditions\n",
" x = [x_init]\n",
" v = [0.]\n",
" a = [0.]\n",
" t = [0.]\n",
" rho = [std_at(x[-1])[1]]\n",
" p_a = [std_at(x[-1])[0]]\n",
" T_a = [std_at(x[-1])[2]]\n",
" m_prop = [propellant_mass(A, L)]\n",
" m = [m_dry + m_prop[-1]]\n",
" (F, A_t, A_e, Ve) = thrust(x[-1], p_ch, T_ch, p_e, ke, Re, mdot)\n",
" F = [F]\n",
" D = [0.]\n",
" Ma = [0.]\n",
" q = [0.]\n",
" r = (m_prop[0] + m_dry)/m_dry # Mass ratio\n",
" dV1 = Ve*log(r)/1000 # Tsiolkovsky's bane (delta-V)\n",
"\n",
" while True:\n",
" p_a.append(std_at(x[-1])[0])\n",
" rho.append(std_at(x[-1])[1])\n",
" T_a.append(std_at(x[-1])[2])\n",
" # Check of the propellant tanks are empty\n",
" if m_prop[-1] > 0:\n",
" (Fr, A_t, A_e, Ve) = thrust(x[-1], p_ch, T_ch, p_e, ke, Re, mdot)\n",
" F.append(Fr)\n",
" m_prop.append(m_prop[-1] - mdot*dt)\n",
" mdot_old = mdot\n",
" else:\n",
" Ve = thrust(x[-1], p_ch, T_ch, p_e, ke, Re, mdot_old)[3]\n",
" F.append(0)\n",
" mdot = 0\n",
" m_prop[-1] = 0\n",
" q.append(drag(x[-1], v[-1], A, Ma[-1], C_d_t, Ma_t)[1])\n",
" D.append(drag(x[-1], v[-1], A, Ma[-1], C_d_t, Ma_t)[0])\n",
" a.append((F[-1] - D[-1])/m[-1] - g_0)\n",
" v.append(a[-1]*dt + v[-1])\n",
" x.append(v[-1]*dt + x[-1]) \n",
" Ma.append(v[-1]/sqrt(ka*Ra*T_a[-1]))\n",
" t.append(t[-1] + dt)\n",
" m.append(m_dry + m_prop[-1])\n",
" TWR = a[1]/g_0 # Thrust-to-weight ratio constraint\n",
" ex = A_e/A_t\n",
" S_crit = p_e/p_a[0] # Sommerfield criterion constraint\n",
" if v[-1] <= 0:\n",
" x = np.array(x)\n",
" a = np.array(a)\n",
" F = np.array(F)\n",
" D = np.array(D)\n",
" q = np.array(q)\n",
" return x, v, a, t, F, D, Ma, rho, p_a, T_a, TWR, ex, Ve, A_t, dV1, m, S_crit, q, m_prop\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Optimization Scheme\n",
"There are several general approaches to design optimization for launch vehicle design found in literature:\n",
"\n",
"* Design of Experiments methods (Taguchi Methods`[1]`, Response Surface Methods`[2]`)\n",
"* Gradient methods (steepest descent)\n",
"* Stochastic methods (genetic algorithm`[3]`, simulated annealing)\n",
"\n",
"While there are certain advantages to each of these approaches, I only have experience with gradient based methods for the most part. Furthermore the difficulty of differentiating the governing equations, the expected multi-modal nature of the response surface and the discrete nature of some of the variables ($m$, $C_d$, etc.) argue against gradient based methods. We therefore use a fast, gradient free, non-stochastic approach: the Nelder-Mead method (e.g. Simplex Search). The limitations of this choice may become evident for high-dimension problems, but for $x \\in \\textbf{R}^4$, we expect the algorithm to have no major problems with convergence. Inequality constraints were handled by using an exterior penalty function. To desensitize the results to arbitrary values of numerically large constraints the pseudo-objective function has been non-dimensionalized."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true,
"scrolled": true
},
"outputs": [],
"source": [
"# Pseudo-objective function\n",
"global cons_TWR, cons_S_crit, cons_accel, cons_LD, cons_alt\n",
"\n",
"#CHANGE CONSTRAINTS HERE\n",
"cons_TWR = 8 # TWR constraint\n",
"cons_S_crit = 0.35 # Critical pressure ratio constraint\n",
"cons_accel = 15 # Max acceleration constraint\n",
"cons_LD = 15 # L/D ratio constraint\n",
"cons_alt = 100000 # Min altitude constraint\n",
"\n",
"def f(x, p_ch=350, rp=50): ##CHANGE CHAMBER PRESSURE HERE\n",
" L = x[0] # Rocket length (m)\n",
" mdot = x[1] # Propellant mass flow rate (kg/s)\n",
" dia = x[2] # Rocket diameter (in)\n",
" p_e = x[3] # Pressure (kPa)\n",
" (alt, v, a, t, F, D, Ma, rho, p_a, T_a, TWR, ex, Ve, A_t, dV1, m, S_crit, q, m_prop) \\\n",
" = trajectory(L, mdot, dia, p_e, p_ch)\n",
"\n",
"\n",
" obj_func = m[0]/m0 + rp*(max(0, (L+2)/(dia*0.0254*cons_LD) - 1)**2 + \\\n",
" max(0, -TWR/cons_TWR + 1)**2 + max(0, -S_crit/cons_S_crit + 1)**2 + \\\n",
" max(0, -alt[-1]/cons_alt + 1)**2 + max(0, max(abs(a))/(cons_accel*9.81) - 1)**2)\n",
" return obj_func\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"##CHANGE INITIAL DESIGN GUESS HERE\n",
"L = 1.5 # Rocket length (m)\n",
"mdot = 0.453592 * 0.9 * 5 # Propellant mass flow rate (kg/s)\n",
"dia = 9 # Rocket diameter (in)\n",
"p_e = 40 # Pressure (kPa)\n",
"X0 = np.array([L, mdot, dia, p_e])\n",
"m0 = trajectory(X0[0], X0[1], X0[2], X0[3])[-4][0] # Initial guess GLOW\n",
"\n",
"from scipy.optimize import minimize\n",
"res = minimize(f, X0, method='nelder-mead') "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"# Plot the results\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import rc\n",
"import pylab\n",
"%config InlineBackend.figure_formats=['svg']\n",
"%matplotlib inline\n",
"rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})\n",
"rc('text', usetex=True)\n",
"\n",
"pylab.rcParams['figure.figsize'] = (10.0, 10.0)\n",
"\n",
"def plot_trajectory(res):\n",
" p_ch = 350 # Chamber pressure [psi] **DONT FORGET TO CHANGE THE VALUE IN THE OBJECTIVE FUNCTION IN def f()**\n",
" (alt, v, a, t, F, D, Ma, rho, p_a, T_a, TWR, ex, Ve, A_t, dV1, m, S_crit, q, m_prop) = \\\n",
"trajectory(res.x[0], res.x[1], res.x[2], res.x[3], p_ch)\n",
" # Rename the optimized output\n",
" L = res.x[0] \n",
" mdot = res.x[1]\n",
" dia = res.x[2]\n",
" p_e = res.x[3]\n",
"\n",
" f, (ax1, ax2, ax3, ax4, ax6, ax7, ax8) = plt.subplots(7, sharex=True)\n",
" for n in (ax1, ax2, ax3, ax4, ax6, ax7, ax8):\n",
" n.spines['top'].set_visible(False)\n",
" n.spines['right'].set_visible(False)\n",
" n.yaxis.set_ticks_position('left')\n",
" n.xaxis.set_ticks_position('bottom')\n",
" n.yaxis.labelpad = 20\n",
" \n",
" ax1.plot(t, alt/1000, 'k')\n",
" ax1.set_ylabel(\"Altitude (km)\")\n",
" ax1.yaxis.major.locator.set_params(nbins=6)\n",
" #ax1.set_title('LV4 Trajectory')\n",
" ax2.plot(t, v, 'k')\n",
" ax2.yaxis.major.locator.set_params(nbins=6)\n",
" ax2.set_ylabel(\"Velocity (m/s)\")\n",
" ax3.plot(t, a/9.81, 'k')\n",
" ax3.yaxis.major.locator.set_params(nbins=10)\n",
" ax3.set_ylabel(\"Acceleration/g0\")\n",
" ax4.plot(t, F/1000, 'k')\n",
" ax4.yaxis.major.locator.set_params(nbins=6)\n",
" ax4.set_ylabel(\"Thrust (kN)\")\n",
" ax6.plot(t, q/1000, 'k')\n",
" ax6.yaxis.major.locator.set_params(nbins=6)\n",
" ax6.set_ylabel(\"Dynamic Pressure (kPa)\")\n",
" ax7.plot(t, Ma, 'k')\n",
" ax7.yaxis.major.locator.set_params(nbins=6) \n",
" ax7.set_ylabel(\"Mach number\")\n",
" ax7.set_xlabel(\"t (s)\")\n",
" ax8.plot(t, np.array(m)*0.666*np.array(a), 'k')\n",
" #ax8.yaxis.major.locator.set_params(nbins=6) \n",
" ax8.set_ylabel(\"LOX Tank Axial Load\")\n",
" ax8.set_xlabel(\"t (s)\")\n",
" plt.savefig(traj.svg)\n",
" plt.show()\n",
" \n",
" \n",
" np.set_printoptions(precision=3)\n",
" print('\\n')\n",
" print('OPTIMIZED DESIGN VECTOR')\n",
" print('-----------------------------')\n",
" print('x_optimized = ', res.x)\n",
" print('x_initial_guess = ', X0)\n",
" print('design tankage length = {0:.2f} m'.format(res.x[0]))\n",
" print('design mass flow rate = {0:.2f} kg/s'.format(res.x[1]))\n",
" print('design airframe diameter = {0:.2f} in.'.format(res.x[2]))\n",
" print('design nozzle exit pressure = {0:.2f} kPa'.format(res.x[3]))\n",
" print('iterations =', res.nit)\n",
" print('design GLOW = {0:.1f} kg'.format(m[0]))\n",
" print('x0 GLOW = {0:.1f} kg'.format( \\\n",
" trajectory(X0[0], X0[1], X0[2], X0[3], p_ch)[-4][0]))\n",
"\n",
" print('\\n')\n",
" print('CONSTRAINTS')\n",
" print('-----------------------------')\n",
" print('L/D ratio (c.f. < {}) = {:.2f}'.format(cons_LD, (L+2)/(dia*0.0254)))\n",
" print('Sommerfield criterion (c.f. pe/pa >= {}) = {:.1f}'.format(cons_S_crit, S_crit))\n",
" print(\"Max acceleration (c.f. < {}) = {:.2f} g's\".format(cons_accel, max(abs(a))/9.81))\n",
" print('TWR at lift off (c.f. > {}) = {:.2f}'.format(cons_TWR, TWR))\n",
" print('altitude at apogee = {0:.1f} km'.format(alt[-1]/1000))\n",
"\n",
" print('\\n')\n",
" print('ADDITIONAL INFORMATION')\n",
" print('-----------------------------')\n",
" print('mission time at apogee = {0:.1f} s'.format(t[-1]))\n",
" print('design total propellant mass = {0:.3f} kg'.format(m_prop[0]))\n",
" print('design thrust (sea level) = {0:.1f} kN'.format(F[0]/1000))\n",
" j = 0\n",
" for thing in F:\n",
" if thing == 0:\n",
" fdex = j\n",
" break\n",
" j += 1\n",
" print('design thrust (vacuum) = {0:.1f} kN'.format(F[fdex - 1]/1000))\n",
" print('design burn time = {} s'.format(fdex))\n",
" print('design expansion ratio = {0:.1f}'.format(ex))\n",
" print('design throat area = {0:.1f} in.^2'.format(A_t/0.0254**2))\n",
" print('design isp = {0:.1f} s'.format(Ve/9.81))\n",
" print('design chamber pressure = {0:.1f} psi'.format(p_ch))\n",
" print('design total impulse = {0:.1f} kN*s'.format(j*(F[fdex - 1]/1000 + F[0]/1000)/2))\n",
" print('design dV = {0:.1f} km/s'.format(dV1))\n",
" print('estimated minimum required dV = {0:.1f} km/s'.format( sqrt(2*9.81*alt[-1])/1000))"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"OPTIMIZED DESIGN VECTOR\n",
"-----------------------------\n",
"x_optimized = [ 1.337 3.444 11.306 107.933]\n",
"x_initial_guess = [ 1.5 2.041 9. 40. ]\n",
"design tankage length = 1.34 m\n",
"design mass flow rate = 3.44 kg/s\n",
"design airframe diameter = 11.31 in.\n",
"design nozzle exit pressure = 107.93 kPa\n",
"iterations = 464\n",
"design GLOW = 115.8 kg\n",
"x0 GLOW = 90.5 kg\n",
"\n",
"\n",
"CONSTRAINTS\n",
"-----------------------------\n",
"L/D ratio (c.f. < 15) = 11.62\n",
"Sommerfield criterion (c.f. pe/pa >= 0.35) = 1.1\n",
"Max acceleration (c.f. < 15) = 15.90 g's\n",
"TWR at lift off (c.f. > 8) = 5.40\n",
"altitude at apogee = 97.8 km\n",
"\n",
"\n",
"ADDITIONAL INFORMATION\n",
"-----------------------------\n",
"mission time at apogee = 153.0 s\n",
"design total propellant mass = 86.106 kg\n",
"design thrust (sea level) = 7.3 kN\n",
"design thrust (vacuum) = 8.0 kN\n",
"design burn time = 26 s\n",
"design expansion ratio = 3.5\n",
"design throat area = 3.7 in.^2\n",
"design isp = 213.5 s\n",
"design chamber pressure = 350.0 psi\n",
"design total impulse = 198.6 kN*s\n",
"design dV = 2.9 km/s\n",
"estimated minimum required dV = 1.4 km/s\n"
]
}
],
"source": [
"from IPython.core.interactiveshell import InteractiveShell\n",
"InteractiveShell.ast_node_interactivity = \"all\"\n",
"\n",
"plot_trajectory(res)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Fine Tuning\n",
"#res.x = [ 1.585, 2.58, 8.987, 60]\n",
"\n",
"#res.x = [ 1.259, 3.11, 10.847, 150.025] #Is an especially good design vector.\n",
"#plot_trajectory(res)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Discussion of Preliminary Results\n",
"As we would expect using Nelder-Mead, the optimized design vector is quite sensitive to the initial guesses. The results appear to always be \"feasible\" designs, but often several of the constrants appear to be simultaneously active. We attempted to find a more global optima by brute forcing a large array of initial guesses using Montecarlo techniques, and examining the frequency spectra to get a sense of the modality of the response surface. We noticed that there didn't appear to be any regimes of feasible designs with significant improvements in GLOW. There is a multi-modality of the respose surface, but given the parameters in the model the best LV configurations all have GLOW on the order of 100 kg or higher.\n",
"\n",
"By changing values in the `dry_mass` function we can see that there is an enormous sensitivty of GLOW to tank mass and inert dry mass, as we would have expected. It is also worthy of note that our optimized designs usually yeild a GLOW of less than 200 kg and total impulse of less than 300 kN s, which would surely make LV4 the smallest ever suborbital rocket to use liquid bi-propellant propulsion. This is a direct result of using our extremely lightweight composite modules as propellant tanks (with tank factors ~10,000). While increasing the chamber pressure to 500 psi only increases performance by about 10% or so (by various metrics), the real advantage of the electric feed system project is that we get more pressure margin for regenerative cooling while simultaneously buying down risk on the composite propellant tanks (which now might only require 2-3 bar ullage pressurization vs. 20 or so), though this obviously depends on how much dry mass is added by the motors, centrifugal impellers, volutes, batteries and power conditioning system (likely a 100+ V DC bus). Another obvious takeaway is in the TWR constraint. At TWR of 2, LV4 will have $<0$ caliber stability margin off the launch rail. This has the advantage of drastically reducing our propulsion requirements. The significant downside is that full 6-DoF attitude control is now a critical path tech requirement (this would need to be achieved using thrust vanes or by engine gimballing).\n",
"\n",
"## A Final Warning\n",
"At present there are some substantial implied limitations to the model as a useful design tool for LV4. The dry mass is based on some very rough estimations, and the optimization results are very sensitive to it. I've tried to deal with this by being conservative in my estimates. The tank mass model is probably optimistic in the sense that it uses a scalar multiple (1.5x) of the current LV3 carbon fiber planar density, plus some allowances for the mass of the machined endcaps as a guide to LV4 tank weights. It also does not take into account ullage volume and other neccessary wasted space that nonetheless have the effect of increasing the overall tank mass. We presently do not know if the LV3 CF technology will work as supposed as a pressure vessel in the LV4 diamaters, post-hoc FEA (and probably some empirical work) is needed to verify this assumption. Another issues is in the drag coeficients used in the drag model. These presently use historical data from the Aerobee-150 sustainer stage. While this may work as a first approximation, it is likely a significant source of error. This can be remedied by using a substantially more complex aerodynamics model in this simulation (perhaps using Barrowmans method) or otherwise post-hoc analysis in Star-CCM+ or another CFD program. The drag coefficients are expected to be largely functions of the fin geometry, which itself depends on static stability analysis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Alternative Optimization Code\n",
"Our own homebrewed Nelder-Mead algorithem. It returns similar design vectors as the Nelder-Mead scheme built into `scipy.optimize` but is much less 'optimized' for computation time (sorry for the pun!)."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Class simplex: \n",
"# Nelder-Mead simplex search\n",
"import numpy as np\n",
"from math import sqrt, pi, exp, log, cos\n",
"import math as m\n",
"\n",
"def search(f, x_start, max_iter = 100, gamma = 5, beta = 0.5, rp=100, a=10, epsilon = 1E-6):\n",
" \"\"\"\n",
" parameters of the function:\n",
" f is the function to be optimized\n",
" x_start (numpy array) is the initial simplex vertices\n",
" epsilon is the termination criteria\n",
" gamma is the contraction coefficient\n",
" beta is the expansion coefficient\n",
" \"\"\"\n",
" # Init Arrays\n",
" N = len(x_start) # Amount of design variables\n",
" fb = [] # Empty function matrix\n",
" xnew = [] # Empty re-write for design variables\n",
" x = [] # Empty x matrix\n",
" C = [[0]*N]*(N+1) # Empty center point matrix #####CHANGED\n",
" \n",
" # Generate vertices of initial simplex\n",
" x0 = (x_start) # x0 Value for x Matrix\n",
" x1 = [x0 + [((N + 1)**0.5 + N - 1.)/(N + 1.)*a, 0., 0., 0.]]\n",
" x2 = [x0 + [0., ((N + 1)**0.5 - 1.)/(N + 1.)*a, 0., 0.]]\n",
" x3 = [x0 + [0., 0., ((N + 1)**0.5 - 1.)/(N + 1.)*a, 0.]]\n",
" x4 = [x0 + [0., 0., 0., ((N + 1)**0.5 - 1.)/(N + 1.)*a]]\n",
" x = np.vstack((x0, x1, x2, x3, x4))\n",
"\n",
" # Simplex iteration\n",
" while True:\n",
" # Find best, worst, 2nd worst, and new center point\n",
" f_run = np.array([f(x[0], rp), f(x[1], rp), f(x[2], rp), f(x[3], rp), f(x[4], rp)]).tolist() \\\n",
" # Func. values at vertices\n",
" xw = x[f_run.index(sorted(f_run)[-1])] # Worst point\n",
" xb = x[f_run.index(sorted(f_run)[0])] # Best point\n",
" xs = x[f_run.index(sorted(f_run)[-2])] # 2nd worst point \n",
" for i in range(0, N+1):\n",
" if i == f_run.index(sorted(f_run)[-1]):\n",
" C[i] = [0,0,0,0]\n",
" else:\n",
" C[i] = x[i].tolist()\n",
" xc = sum(np.array(C))/(N) # Center point\n",
" xr = 2*xc - xw # Reflection point\n",
" fxr = f(xr, rp)\n",
" fxc = f(xc, rp)\n",
" \n",
" # Check cases\n",
" # f(xr, rp) < f(xb, rp): # Expansion\n",
" if fxr < f_run[f_run.index(sorted(f_run)[0])]: \n",
" xnew = (1 + gamma)*xc - gamma*xr\n",
" # f(xr, rp) > f(xw, rp): # Contraction 1\n",
" elif fxr > f_run[f_run.index(sorted(f_run)[-1])]:\n",
" xnew = (1 - beta)*xc + beta*xw\n",
" # f(xs, rp) < f(xr, rp) and f(xr, rp) < f(xw, rp): # Contraction 2\n",
" elif f_run[f_run.index(sorted(f_run)[-2])] < fxr and fxr < f_run[f_run.index(sorted(f_run)[-1])]: \n",
" xnew = (1 + beta)*xc - beta*xw\n",
" else:\n",
" xnew = xr\n",
" \n",
" # Replace Vertices\n",
" x[f_run.index(sorted(f_run)[-1])] = xnew\n",
" #x[f_run.index(sorted(f_run)[1])] = xb # Replace best\n",
" #x[f_run.index(sorted(f_run)[2])] = xs # Replace second best\n",
" fb.append(f(xb, rp))\n",
" print('Current optimum = ', fb[-1])\n",
" \n",
" # Break if any termination critera is satisfied\n",
" if len(fb) == max_iter: #or term_check(x, xc, xw, N, rp, f_run) <= epsilon:\n",
" (alt, v, a, t, F, D, Ma, rho, p_a, T_a, TWR, ex, Ve, A_t, dV1, m, S_crit, q, m_prop, p_ch) \\\n",
" = trajectory(xb[0], xb[1], xb[2], xb[3])\n",
" return f(x[f_run.index(sorted(f_run)[0])], rp), x[f_run.index(sorted(f_run)[0])], len(fb)\n",
" \n",
"def term_check(N, rp, f_run, fxc): # Termination critera\n",
" M = [0]*(N + 1)\n",
" for i in range(0, N + 1):\n",
" if i == f_run.index(sorted(f_run)[-1]): # Avoid worst point\n",
" M[i] = 0\n",
" else:\n",
" M[i] = (f_run[i] - fxc)**2\n",
" return m.sqrt(sum(M)/N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"1. Stanley, D. O., Unal, R., and Joyner, C. R., *``Application of Taguchi Methods to Dual Mixture Ratio Propulsion System Optimization for SSTO Vehicles\"*, Journal of Spacecraft and Rockets, Vol. 29, No. 4, 1992, pp. 453-459.\\\\\n",
"\n",
"2. Stanley, D. O., Engelund. W. C., Lepsch. R. A., McMillin, M. L.Wt K. E.. Powell. R. W., Guinta. A. A., and Unal, R., *``Rocket-Powered Single Stage Vehicle Configuration Selection and Design\"*, Journal of Spacecraft and Rockets,Vol. 31, No. 5, 1994. pp. 792-798; also AIAA Paper93-Feb. 1993.\\\\\n",
"\n",
"3. Anderson, M., Burkhalter J., and Jenkins R., *``Multidisciplinary Intelligence Systems Approach To Solid Rocket Motor Design, Part I: Single And Dual Goal Optimization\"*. AIAA 2001-3599, July, 2001.\\\\"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}