Oregon Undergraduate Research Journal: McNair Special Issue 21.3 (2023) ISSN: 2160-617X (online) ourj.uoregon.edu Multidisciplinary Design Optimization of Portland State Aerospace Society (PSAS) Launch Vehicle 4 Aaron Casserly* Abstract Multidisciplinary Design Optimization is a field that enables the solution of challenging engineering problems involving multiple technical specializations and design/performance constraints. In this work, I optimize the design of the PSAS Launch Vehicle 4 (LV4). To that end, I evaluate different optimization approachesβ€”such as RBFOpt Global Optimization, Nelder- Mead minimization, and Simplicial Homology Global Optimization with Nelder-Mead and COBYLA local minimization techniques, calculate structural analysis information for different stages of flight, outline a method of simulating fin β€œstaging”—the dropping of a larger initial fin can at a certain altitude to reduce the required engine thrust and drag in the upper atmosphere and optimize fin parameters. I converged on the ideal design vector. This led to an apogee of 107 km with a 9.8 kN engine (realized with two 5 kN engines). Further debugging is required to resolve the apparent 120 km vehicle drift. 1. Introduction simulation code and to converge on a design vector satisfying engineering and performance The Portland State Aerospace Society (PSAS) constraints. Launch Vehicle 4 is the fourth iteration of their 2. Methodology student-built rocket. PSAS is a feeder organization for those interested in working in Our approach minimizes the gross lift-off weight spaceflight and we care about optimizing it to (GLOW) of LV4 without sacrificing apogee. Given give young professionals experience in relevant the Tsiolkovsky rocket equation 𝛿𝛿𝛿𝛿 = 𝑙𝑙𝑙𝑙 π‘šπ‘š0 , where problem solving to real-world problems, and to 𝛿𝛿𝑒𝑒 π‘šπ‘šπ‘“π‘“ achieve the above 100 km apogee target while 𝑣𝑣𝑒𝑒 represents the effective exhaust velocity, 𝛿𝛿𝑣𝑣 meeting performance constraints. The the total change in velocity, π‘šπ‘š0 the GLOW and π‘šπ‘šπ‘“π‘“ constraints given to the optimizer are real-world the final dry mass of the rocket, we isolate π‘šπ‘š0 and engineering constraints that must be met for this construct an objective function. This approach to be realized. Previous optimization efforts has previously been characterized: included a design vector involving various components of the rocket, and my approach to [M]inimizing GLOW while demanding a include fin parameters and scale them down certain apogee is equivalent to when a near optimal design was converged on simultaneously minimizing structural mass, was an improvement. The goals of this research maximizing engine performance, and project were to improve and extend the existing balancing the conflicting goals of minimizing Multidisciplinary Design Optimization (MDO) losses due to gravity and aerodynamics. Note *Aaron Casserly (casser.aero@gmail.com) is a Jamaican immigrant and lawful permanent resident of the United States. He arrived in Oregon in April of 2019 and enrolled at the University of Oregon later that year. He is now a graduate with a Bachelor of Science in Mathematics and double minor in Physics and Computer Science with final GPA of 3.7/4. He plans to attend Northwestern University for a master’s degree in electrical engineering in the Fall of 2024. Oregon Undergraduate Research Journal: McNair Special Issue Casserly that the sources of this conflict are the so that it may find a suitable neighbourhood and incentive to expel propellant rapidly to avoid then restricts its freedom once it has done so. the cost of carrying propellant in a The All-at-Once (AAO) problem statement is a gravitational field and the incentive to reduce fundamental optimization problem from which velocity in lower atmosphere since drag is all others may be derived. It includes an proportional to air pressure and the square of objective/pseudo-objective function to be velocity. (MDO Jupyter Notebook) minimized with respect to a design vector and subject to certain constraints. For this problem, To represent our problem constraints, we use we are estimating the optimal design vector, ?Μ…?π‘₯ barrier and penalty functions. Both bind the according to π‘™π‘™π‘™π‘™π‘šπ‘š μ𝑛𝑛 = 0, π‘™π‘™π‘™π‘™π‘šπ‘š ρ = ∞, 𝑓𝑓 (?Μ…?π‘₯) =π‘›π‘›β†’βˆž π‘›π‘›β†’βˆž 𝑛𝑛 𝑛𝑛 generated rockets to a feasible region in the π‘šπ‘šπ‘œπ‘œπ‘œπ‘œπ‘œπ‘œ(?Μ…?π‘₯) + μ𝑛𝑛 βˆ‘π‘–π‘– β„Žπ‘–π‘–(?Μ…?π‘₯) + ρ𝑛𝑛 βˆ‘π‘œπ‘œ π‘”π‘”π‘œπ‘œ(?Μ…?π‘₯). design space: barrier functions use an absolute xοΏ½n = π‘šπ‘šπ‘™π‘™π‘™π‘™?Μ…?π‘₯ 𝑓𝑓𝑛𝑛(?Μ…?π‘₯), π‘™π‘™π‘™π‘™π‘šπ‘š xοΏ½n = xοΏ½βˆ—, where ?Μ…?π‘₯ = constraint that may not be violated under any π‘›π‘›β†’βˆž circumstances, whereas penalty functions (π‘šπ‘šπ‘π‘π‘π‘π‘œπ‘œπ‘π‘, ?Μ‡?π‘š, 𝑝𝑝𝑒𝑒) is the design vector containing the slightly disincentivize convergence in sections of total propellant mass, unadjusted propellent the design space far from a set of more lenient mass, mass flow rate, nozzle exit pressure, total constraints. tankage length, airframe diameter, airframe total Adding the objective, barrier, and penalty length, GLOW, ballast mass, conical component functions, we construct a pseudo-objective merit of nosecone length, fin root chord, fin tip chord, function, which takes in an array of values fin sweep angle, fin span, and fin thickness. This sufficient to describe the mathematical model of information is necessary for the evaluation of the the rocket and its performance. Given that each constraint functions, β„Žπ‘œπ‘œπ‘π‘π‘π‘π‘π‘π‘–π‘–π‘’π‘’π‘π‘ = 108401π‘šπ‘š < β„Ž < evaluation of this merit function consists of 151401π‘šπ‘š, and 𝑔𝑔𝑝𝑝𝑒𝑒𝑛𝑛𝑏𝑏𝑝𝑝𝑝𝑝𝑝𝑝 = �𝐹𝐹 ≀ 6π‘˜π‘˜π‘˜π‘˜, 𝐿𝐿𝐿𝐿 β‰₯ simulating the rocket’s trajectory, we are unable 22π‘šπ‘š/𝑠𝑠, π‘π‘π‘šπ‘šπ‘šπ‘šπ‘šπ‘š ≀ 15𝑔𝑔′𝑠𝑠,𝑇𝑇𝑇𝑇𝑇𝑇 β‰₯ 2, 𝐿𝐿/𝐷𝐷 ≀ 21, 𝑝𝑝𝑒𝑒 β‰₯ 𝑔𝑔0 π‘π‘π‘šπ‘š to use an optimization algorithm involving 0.35οΏ½. β„Ž is the strict apogee constraint, and differentiation or a finite difference method. To π‘œπ‘œπ‘π‘π‘π‘π‘π‘π‘–π‘–π‘’π‘’π‘π‘ navigate this limitation, we use a Nelder-Mead the looser penalty constraints are: Thrust (F)β€”the simplex method: a geometry-based optimization Electric Feed System (EFS) that deals with algorithm that does not perform well for higher pressurizing before propellant injection is not dimensions but is satisfactory for our purposes. A feasible for powerful engines, Launch Speed genetic algorithm, which can handle higher- (LS)β€”ensuring stable take-off, Thrust to Weight dimensional spaces, would also serve this Ratio (TWR), Length to Diameter Ratio (L/D), function, but it incurs an additional maximum acceleration, and nozzle over- computational cost. expansion. The barrier and penalty functions are The unaltered optimization code produces an weighted by user-selected parameters before .ork rocket file for further testing in OpenRocket, their addition to the objective function. Given as seen in Figure 1. that an overly low weighting will lead to the optimizer’s neglect of the constraints and an overly high weighting to the optimizer ignoring the objective, we run an iterative sequence of Nelder-Mead optimizations, starting with very low weights and increasing them for every successive optimization. This gives the optimizer more global coverage of the design space early on Figure 1. Example generated rocket. Multidisciplinary Design Optimization 9 Oregon Undergraduate Research Journal: McNair Special Issue Casserly It is then possible to run the simulation with operational environment; robustnessβ€”measure settings emulating the launch site (WGS84 of insensitivity to variations in both the system ellipsoid for Geodetic calculations), pictured in and environment; reliabilityβ€”likelihood of a Figure 2. component/system to perform intended function for a given period of time under the determined operating conditions; deterministic design optimizationβ€”process of obtaining optimal designs with all variables, models, parameters and simulations involved being deterministic; robust design optimization (RDO)β€”optimizing design such that it is insensitive to many variations; and reliability-based design optimization (RBDO)β€”obtaining optimal design while meeting reliability constraints. The combination of these last two, RDO and RBDO, is the basis for reliability-based robust design Figure 2. Settings to emulate the launch site. optimization (RBRDO): Find 𝒙𝒙 minimizing 𝑓𝑓(𝒙𝒙,𝒑𝒑) = 𝐹𝐹 �μ𝑓𝑓(𝒙𝒙,𝒑𝒑),σ𝑓𝑓(𝒙𝒙,𝒑𝒑)οΏ½ subject to (s.t.) These calculations produce a launch data 𝑃𝑃[π’ˆπ’ˆ(𝒙𝒙,𝒑𝒑) ≀ 0] β‰₯ 𝑹𝑹, 𝒙𝒙𝐿𝐿 ≀ 𝒙𝒙 ≀ π’™π’™π‘ˆπ‘ˆ, where 𝒙𝒙 graphic, pictured in Figure 3. represents the design variable vector, 𝒑𝒑 represents the system constant parameter vector, 𝒙𝒙𝐿𝐿 and π’™π’™π‘ˆπ‘ˆ define the boundaries of the design space, μ𝑓𝑓 and σ𝑓𝑓 are the mean and standard deviation of the original optimization objective function, F is the reformulated optimization function with respect to μ𝑓𝑓 and σ𝑓𝑓, g is the unequal constraint vector, P is the probability of the statement in brackets to be true, and R is the Figure 3. OpenRocket launch data graphic. reliability vector related to this. Yao et al. provide illustrations related to RDO and RBDO, seen in Review of uncertainty-based multidisciplinary Figures 4 and 5. design optimization methods for aerospace vehicles (Yao et al.) covers Uncertainty-Based Multidisciplinary Design Optimization (UMDO) theory and the cutting edge methods of that time. Throughout the lifecycle of the aerospace vehicle (design, manufacture, operation, disposal/repurposing), there exist many uncertainties related to the vehicle system itself, along with environmental and operational conditions. Before describing the UMDO procedure, several important definitions are given: uncertaintyβ€”incompleteness in knowledge Figure 4. Graphical illustration of RDO. and inherent variability of the system and Multidisciplinary Design Optimization 10 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Figure 6. The coupling relationship of a three-discipline UMDO problem. Figure 5. Graphical illustration of RBDO. Yao et al. also provide an illustration of the The actual UMDO procedure organizes the conventional double-loop UMDO procedure, elements involved in uncertainty-based design pictured in Figure 7. optimization: system optimization, system analysis, disciplinary analysis, and uncertainty analysis. The key to realizing UMDO for a large, complex system is efficiently arranging these elements into an execution sequence so that it may be implemented on a computer; the coupling relationship related to disciplinary analysis and computationally intensive system analysis make for a very time-consuming procedure. The computational burden of UMDO can be understood by the following modification Figure 7. The conventional double-loop UMDO procedure. to the RBRDO formulation: Find 𝒙𝒙 minimizing 𝑓𝑓(𝒙𝒙,𝒑𝒑,π’šπ’š) = 𝐹𝐹 οΏ½ΞΌ (𝒙𝒙,𝒑𝒑,π’šπ’š),Οƒ (𝒙𝒙,𝒑𝒑,π’šπ’š)οΏ½ s.t. I set out to review the existing open-source 𝑓𝑓 𝑓𝑓 optimization code to analyze the dependencies 𝑃𝑃[𝑔𝑔𝑖𝑖(𝒙𝒙,𝒑𝒑,π’šπ’š) ≀ 0] β‰₯ π‘Ήπ‘Ήπ’Šπ’Š, 𝑙𝑙 = 1,2, … ,𝑙𝑙𝑔𝑔, 𝒙𝒙𝐿𝐿 ≀ 𝒙𝒙 ≀ π’™π’™π‘ˆπ‘ˆ. and components of these procedures such that I π’šπ’š represents the intermediate state variables of would be able to isolate areas for potential the multidisciplinary analysis. We denote the improvements. In making modifications, I set out output vector of disciplinary analysis i as π’šπ’šπ‘–π‘–, the to compare simulated trajectory results and the coupling state vector output from disciplinary optimized design vectors. The following is a list of analysis i and input into disciplinary analysis j as planned improvements: π’šπ’šπ‘–π‘–π‘œπ‘œ, the complete set of output vectors from discipline 𝑙𝑙 coupled with other disciplines π’šπ’šπ‘–π‘–., and 1. The majority of the Aerodynamics model the complete set of coupling state vectors input is based on OpenRocket’s source code, into disciplinary analysis i as π’šπ’š.𝑖𝑖. With these which is an oversimplification of reality. conventions, we have π’šπ’š = [π’šπ’šπ‘–π‘– , 𝑙𝑙 = 1,2, … ,𝑙𝑙𝐷𝐷], π’šπ’š.𝑖𝑖 = The first instance of a possible οΏ½π’šπ’šπ‘œπ‘œπ‘–π‘– , 𝑗𝑗 = 1,2, … ,𝑙𝑙𝐷𝐷, 𝑗𝑗 β‰  𝑙𝑙�, π’šπ’šπ‘–π‘–. = π’šπ’šπ‘–π‘–.(𝒙𝒙𝑖𝑖 ,𝒑𝒑,π’šπ’š.𝑖𝑖) and improvement to the Jupyter Notebook 𝑦𝑦𝑖𝑖. = οΏ½π‘¦π‘¦π‘–π‘–π‘œπ‘œ , 𝑗𝑗 = 1,2, … ,𝑙𝑙𝐷𝐷, 𝑗𝑗 β‰  𝑙𝑙�. 𝒙𝒙𝑖𝑖 is the local detailing the aerodynamics model is design variable vector of discipline i, and 𝒑𝒑 is the related to the fin-body interference system parameter vector. The paper provides a coefficient. This can be improved using figure with information related to the coupling MIL-HDBK-762 (Handbook). relationship for a three-discipline system, 2. The optimization program generates pictured in Figure 6. designs based on a template OpenRocket Multidisciplinary Design Optimization 11 Oregon Undergraduate Research Journal: McNair Special Issue Casserly file. It will be necessary to update this to account for the additional complexity of reality. the latest airframe design, as the design The first improvement I made in relation to this fed into the optimizer is still based on the was to adjust the fin-body interference coefficient airframe of a previous iteration. using MIL-HDBK-762. Using the slender-body 3. Implement fin β€œstaging”; this will be done theory approach (where the slenderness of the by dropping a large fin can. Essentially, modelled body is used to create an approximation larger fins enable a reduction in the to the field surrounding it), the ratio of the fin required launch velocity, and, therefore, normal force gradientβ€”the resulting corrective the engine size may be reduced. force perpendicular to the z-axis of the rocketβ€”in Dropping the can mid-flight will reduce the presence of a cylindrical body compared to the drag due to large fins. I will develop a that of an isolated fin is given by 𝐾𝐾𝐹𝐹(π‘œπ‘œ) = method of simulating this effect. 2 𝑑𝑑4Ο€ 1 1 π‘œπ‘œ 𝑑𝑑 Ο€ 4. Add to MDO capability via reading-in a 2 οΏ½οΏ½1 + 4οΏ½ οΏ½ π‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘™π‘™ οΏ½ οΏ½ βˆ’ οΏ½οΏ½ + οΏ½ βˆ’οΏ½1βˆ’π‘‘π‘‘οΏ½ π‘œπ‘œ 2 2 𝑝𝑝 π‘œπ‘œ 4𝑏𝑏 database of aerodynamic coefficients 𝑑𝑑2 π‘œπ‘œ created by CFD. 2 οΏ½οΏ½ βˆ’ 𝑑𝑑� + 2π‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘Žπ‘™π‘™ �𝑑𝑑��� , where π‘Žπ‘Ž is half the π‘œπ‘œ 𝑑𝑑 π‘œπ‘œ π‘œπ‘œ 5. Improve the β€œUI” such that those fin span (𝑏𝑏) and 𝑑𝑑 is the body diameter, as unfamiliar with the code may set pictured in Figure 8. parameters and understand the process. Work on documentation and comments in relation to better user interaction. 6. Add structural analysis output such as weight, stresses, acceleration, acceleration in propellants, axial load down the rocket, and heatmap to show where the axial load is the highest. 7. If efforts to reduce engine weight have Figure 8. Fin span and body diameter of LV4. failed and we need a large engine (on the order of 10 kN), the use two 5 kN engines The trajectory simulation component of the or four 2.5 kN engines would be optimal. open-source code outputs the angle of attackβ€” the difference between the rocket’s z-axis and This would reduce chamber pressure; relative velocity vectorβ€”of the vehicle as a additionally, lower-thrust engines are more feasible to build. function of time, seen in Figure 9. 8. Compare efficiency/quality (merit function evaluation) of optimization approaches such as global optimization using RBFOpt, iterative Nelder-Mead, and Simplicial Homology Global Optimization (SHGO). 3. Results 3.1. Aerodynamics Figure 9. Angle of attack as a function of time. The majority of the Aerodynamics model is based The optimized design vector before and after on source code from OpenRocket, which does not this modification is represented in Figure 10. Multidisciplinary Design Optimization 12 Oregon Undergraduate Research Journal: McNair Special Issue Casserly expansion is not fully converted into thrust by the nozzle inner wall. Figure 13 demonstrates an under-expanded nozzle. Figure 11. Nozzle and combustion chamber of LV4. Figure 10. Optimized design vector before and after fin-body interference modification. The increased total propellant mass and gross lift-off weight indicate that increased resistance to the fin normal force is required. The higher nozzle exit pressure is related to the optimization program minimizing nozzle over-expansion (Monte). An over/under-expanded nozzle is one in which the exit pressure is greater or lower than the atmospheric pressure. The combustion chamber generates high pressure, high Figure 12. Over-expanded nozzle of LV4. temperature gas, and the ideal nozzle (shape and length optimized) converts this thermal energy into thrust as in Figure 11. An over-expanded nozzle is one in which the atmospheric pressure is greater than the exit pressure, which causes a pinching effect and decreases the efficiency of the nozzle as sections of the nozzle inner wall are not used to produce thrust. Figure 12 demonstrates an over-expanded nozzle. Under- expansion is the opposite: the atmospheric pressure is less than the nozzle exit pressure, which causes the flow to fan out after exiting the nozzle and results in inefficiency, as the Figure 13. Under-expanded nozzle of LV4. Multidisciplinary Design Optimization 13 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Ideally, we would like nozzle exit pressure to construct each variable’s list of values, with equal atmospheric pressure, and minimizing which we populate the dictionary. We are now over-expansion is our best option. The higher exit ready to perform interpolation, as shown in pressure after increasing the accuracy of the fin- Figure 16. body interference coefficient indicates that the atmospheric pressure is higher, likely due to an increased fin normal force gradient. In relation to enabling the MDO to read CFD data, I obtained a demo aerodynamics databaseβ€” shown in Figure 14β€”with the goal of building an interpolation space in the variables. Figure 14. Aerodynamic coefficients database. I created a new notebook to complete this task and used the SciPy interpolate module Figure 16. Cubic interpolation in the qbar variable. (Scipy.Interpolate.Interp1d β€” SciPy v1.8.1 Manual) We use the β€œinterp1d” function to along with the built-in Python CSV module. approximate a continuous function given our Figure 15 displays the first from-scratch code discrete data points. Figure 17 displays another contribution I made. example, this time with linear interpolation. Figure 15. Code for processing database of Aerodynamic coefficients. Figure 17. Linear interpolation in the Cx variable. I initialize a list β€œmatrix” to store the data along with a dictionary β€œDict” to index into a 3.2. Optimization Method variable’s list of values. Using the csv module, we Efficiency/Quality Analysis read in the data and construct a 2-dimensional matrix. We then use a nested structure to The three options for optimization that I considered were global optimization using Multidisciplinary Design Optimization 14 Oregon Undergraduate Research Journal: McNair Special Issue Casserly RBFOpt (black-box optimization), 2. Repeat until termination test is satisfied. scipy.optimize’s Nelder-Mead minimization, and a. Calculate termination test Simplicial Homology Global Optimization information (absolute error). (SHGO). It is important to note that our merit b. If termination test is satisfied, function is a combination of an objective function transform the working simplex. and a penalty function. This means that the 3. Return the best vertex of the current constraints are captured by the penalty simplex 𝐿𝐿 and the merit function component and are not passed into the evaluation. optimization methods directly. Therefore, we have in each case an unconstrained global We can construct the initial simplex by optimization which approximates a constrained generating 𝑙𝑙 + 1 vertices (π‘₯π‘₯0, … , π‘₯π‘₯𝑛𝑛) around some optimization. input point in 𝑹𝑹𝑛𝑛. For practical purposes, we use The Nelder-Mead algorithm is designed to π‘₯π‘₯0 so that the algorithm may be restarted. The minimize a non-linear function 𝑓𝑓:𝑹𝑹𝑛𝑛 β†’ 𝑹𝑹 using remaining n vertices are then generated to obtain function values at a few points in 𝑹𝑹𝑛𝑛. It can be a regular simplex, with all edges having the same viewed as a simplex-based search algorithm. A length. simplex in 𝑹𝑹𝑛𝑛 is defined as the convex hull of 𝑙𝑙 + A key component is the simplex 1 vertices. For example, a simplex in 𝑹𝑹2 is a transformation algorithm, which consists of three triangle, while 𝑹𝑹3 would be a tetrahedron, shown stages: in Figure 18. 1. Ordering to determine the worst (h), second-worst (s), and best (l) vertices in the current working simplex: π‘“π‘“β„Ž = π‘šπ‘šπ‘Žπ‘Žπ‘₯π‘₯π‘œπ‘œπ‘“π‘“π‘œπ‘œ, 𝑓𝑓𝑠𝑠 = π‘šπ‘šπ‘Žπ‘Žπ‘₯π‘₯π‘œπ‘œβ‰ β„Žπ‘“π‘“π‘œπ‘œ, 𝑓𝑓𝑝𝑝 = π‘šπ‘šπ‘™π‘™π‘™π‘™π‘œπ‘œβ‰ β„Žπ‘“π‘“π‘œπ‘œ. 2. Calculation of the centroid of the best side, which is opposite the h-vertex (π‘Žπ‘Ž ≔ 1 βˆ‘ 𝑛𝑛 π‘œπ‘œβ‰ β„Ž π‘₯π‘₯π‘œπ‘œ). 3. Computation of the new working simplex via transforming the current. Figure 18. Simplexes in two- and three-dimensional space. The method begins with a set of points We attempt to replace the worst vertex using π‘₯π‘₯ , … , π‘₯π‘₯ ∈ 𝑹𝑹𝑛𝑛0 𝑛𝑛 , which are the vertices of our reflection, contraction or expansion with respect simplex, and their merit function evaluations. to the best side. The test points lie on the line The algorithm will then perform a series of from the worst point (π‘₯π‘₯β„Ž) to the centroid of the transformations on the working simplex with the best side, as previously calculated. At most, two goal of decreasing the merit function evaluation such points are calculated in each iteration. If at the vertices. This process is terminated when successful, this accepted point becomes the new the absolute errors in the optimal design vector vertex of our working simplex. Otherwise, we and its function evaluation are sufficiently small shrink the simplex towards the best vertex (π‘₯π‘₯𝑝𝑝), (Minimize(Method=’Nelder-Mead’) β€” SciPy v1.8.0 and it is necessary to compute 𝑙𝑙 new vertices. Manual). A simplification of the algorithm is the On testing the Nelder-Mead approach, I following: reduced the termination conditions to an absolute difference between optimal design 1. Construct initial working simplex 𝐿𝐿. vectors of 1 and an absolute difference between Multidisciplinary Design Optimization 15 Oregon Undergraduate Research Journal: McNair Special Issue Casserly merit function evaluations of 0.1. The results indicate that this will produce a locally feasible design; however, this is not ideal for the global optimization that we desire. Combining this with Figure 21. RBFOpt Global Optimization code snippet. SHGO should improve results. The output of this is shown in Figure 19. This design variable range works in conjunction with Bonmin (Bonmin) (Basic Open- source Nonlinear Mixed Integer programming) to find a design vector that minimizes the merit function. With this approach, the majority of our constraints are satisfied, as seen in Figure 22. Figure 19. Iterative Nelder-Mead optimization output. We can see here that the peak altitude is less than desired. Overall, this approach is time- intensive, given a more stringent termination condition, and will at best produce locally feasible designs. The limited design space exploration of the Nelder-Mead algorithm can be understood via Figure 20. Figure 22. RBFOpt Global Optimization Constraint Satisfaction. Simplicial Homology Global Optimization (SHGO), in conjunction with Nelder-Mead minimization, provides good results with a completion time of about 2 hours. The theoretical advantages of SHGO are guaranteed when the objective function is Lipschitz smooth (objective function is continuous, convex, and smooth); however, if this is not the case, the algorithm will converge to the global optimum if the default β€œsimplicial” sampling method is used Figure 20. Iterative Nelder-Mead design space exploration. (Scipy.Optimize.Shgo β€” SciPy v1.8.0 Manual). SHGO is a general-purpose global optimization The RBFOpt (Coin-or/Rbfopt) global algorithm that approximates the homology optimization method provides excellent coverage groups of a complex built on a hypersurface that of the design space and runs to completion in is homeomorphic (similar in form) to a complex under two hours. To perform this optimization, on the objective function. This facilitates we construct a black box using the approximations of locally convex subdomains in RbfoptUserBlackBox class and execute the search space (multidimensional space RbfoptAlgorithm on it. The values comprising the consisting of design vector parameters and their two arrays in the definition of the black box were constraints) and provides an excellent visual tool set based on a feasible range for the design for characterising and solving higher- variables. Figure 21 displays a snippet of our dimensional black box optimization problems. code. The complex is created using sampling points Multidisciplinary Design Optimization 16 Oregon Undergraduate Research Journal: McNair Special Issue Casserly within the feasible search space as vertices. The hyperrectangle (rectangle generalized to algorithm is best suited to finding all local higher dimensions). minima of an objective function with a 2. The set 𝒫𝒫 = [π’™π’™πœ–πœ–πœ’πœ’|π’ˆπ’ˆ(𝒙𝒙) β‰₯ 0] describes a computationally expensive evaluation (such as set of points within the feasible set Ξ©. ours, which involves simulating the trajectory of 3. Given an objective function 𝑓𝑓,β„± the design). represents the set of scalar outputs Using the sampled points of an objective mapped by the objective function 𝑓𝑓: 𝒫𝒫 ⟢ function as vertices, this method constructs a β„± in relation to a sampling set 𝒫𝒫 βŠ† Ξ© βŠ† simplicial complex. The resulting directed ℝ𝑛𝑛. subgraph contains the set of all 1-chains from the 4. If β„‹ is a directed simplicial complex, elements of β„‹1 ∈ β„‹ and enables the finding of then β„‹0 ≔ 𝒫𝒫 is the set of all vertices of minimizer pools (Endres et al.) Sperner’s lemma β„‹. enables us to approximate the domains of 5. Given a set of vertices β„‹0, we construct stationary points for our objective function in the the simplicial complex β„‹ by a feasible search space, denoted by Ξ©. The triangulation connecting every vertex in homology groups produced from the β„‹0. This supplies a set of undirected construction of β„‹ will be invariant given an edges 𝐸𝐸. adequate sampling set. It follows that for the 6. β„‹1 is a set constructed by directing every given sampling set of vertices β„‹0 ∈ β„‹, we are edge in 𝐸𝐸. This is done by selecting a guaranteed to extract the optimal minimiser pool. vertex 𝑣𝑣𝑖𝑖 ∈ β„‹0 and connecting to another The algorithm has four steps: vertex π‘£π‘£π‘œπ‘œ by an edge within 𝐸𝐸. This edge is directed as π‘£π‘£π‘–π‘–π‘£π‘£π‘œπ‘œ from 𝑣𝑣𝑖𝑖 to π‘£π‘£π‘œπ‘œ if and only if 1. Sampling point generation of N vertices the merit function evaluation at the in the search space from which 0-chains former is lesser than the latter. It is of β„‹0 are constructed. directed as π‘£π‘£π‘œπ‘œπ‘£π‘£π‘–π‘– from π‘£π‘£π‘œπ‘œ to 𝑣𝑣𝑖𝑖 if and only if 2. Triangulation of the vertices to construct the merit function evaluation at the the directed simplicial complex β„‹. former is greater than the latter. In these 3. Construction of the minimiser pool using cases, we have πœ•πœ•οΏ½π‘£π‘£ 𝑣𝑣 οΏ½ = 𝑣𝑣 βˆ’ 𝑣𝑣 and Sperner’s lemma. 𝑖𝑖 π‘œπ‘œ π‘œπ‘œ 𝑖𝑖 πœ•πœ•οΏ½π‘£π‘£ 𝑣𝑣 οΏ½ = 𝑣𝑣 βˆ’ 𝑣𝑣 . The case in which 4. Local minimization using the starting π‘œπ‘œ 𝑖𝑖 𝑖𝑖 π‘œπ‘œ points defined in the minimiser pool 𝑓𝑓(𝑣𝑣𝑖𝑖) = π‘“π‘“οΏ½π‘£π‘£π‘œπ‘œοΏ½, with neither 𝑣𝑣𝑖𝑖 nor π‘£π‘£π‘œπ‘œ (Nelder-Mead method). already being a minimizer, we use the rule that β€œthe incidence direction of the Given a set of sampling points 𝒫𝒫, we wish to connecting edge is always directed describe a discrete mapping β„Ž: 𝒫𝒫 βŸΆβ„‹ that will towards the vertex that was generated provide a simplicial approximation for the earliest by the sampling point sequence.” surface of the merit function. To begin, we need If 𝑣𝑣𝑖𝑖 is not connected to another vertex to formally define the set of vertices forming the π‘£π‘£π‘˜π‘˜, then our convention will be to leave 0-chains of the simplicial complex and the edges π‘£π‘£π‘–π‘–π‘£π‘£π‘˜π‘˜ undefined, with πœ•πœ•(π‘£π‘£π‘–π‘–π‘£π‘£π‘˜π‘˜) = 0. The π‘˜π‘˜ forming the 1-chains of β„‹. The following are higher dimensional simplices β„‹ ,π‘˜π‘˜ = useful definitions: 2, 3, …𝑙𝑙 + 1 may be directed in an arbitrary direction to complete the 1. πœ’πœ’ is the set of sampling points created by construction of the complex β„Ž ∢ 𝒫𝒫 βŸΆβ„‹. a sampling sequence in a bounded This will be used to find the minimiser pool for the local minimization starting Multidisciplinary Design Optimization 17 Oregon Undergraduate Research Journal: McNair Special Issue Casserly points required by the algorithm. 7. 𝑣𝑣𝑖𝑖 is a minimiser if and only if all edges connected to 𝑣𝑣𝑖𝑖 are directed away from 𝑣𝑣𝑖𝑖; formally, that is βˆ‚οΏ½π‘£π‘£π‘–π‘–π‘£π‘£π‘œπ‘œοΏ½ = οΏ½π‘£π‘£π‘œπ‘œβ‰ π‘–π‘– βˆ’ 𝑣𝑣𝑖𝑖� ∨ 0 βˆ€v 0jβ‰ i ∈ β„‹ . The set of all minimisers is the minimiser pool β„³. 8. The star of a vertex 𝑣𝑣𝑖𝑖 [π‘ π‘ π‘Žπ‘Ž(𝑣𝑣𝑖𝑖)] is the set of all points 𝒬𝒬 s.t. every simplex containing 𝒬𝒬 contains 𝑣𝑣𝑖𝑖. 9. The π‘˜π‘˜-chain 𝐢𝐢(β„‹π‘˜π‘˜), π‘˜π‘˜ = 𝑙𝑙 + 1 of simplices in st (𝑣𝑣𝑖𝑖) results in a boundary cycle βˆ‚οΏ½πΆπΆ(ℋ𝑛𝑛+1)οΏ½ with βˆ‚ οΏ½βˆ‚οΏ½πΆπΆ(ℋ𝑛𝑛+1)οΏ½οΏ½ = βˆ…. The bounds of the domain defined by Figure 24. Objective function values using sampling points s.t. (𝑣𝑣𝑖𝑖) form the faces of βˆ‚(ℋ𝑛𝑛+1). from Sobol sequence. To place these constructions in a practical From Definition 4 above, we have β„‹0 from 𝒫𝒫. context, we minimize the Ursem01 function in Definition 5 enables us to construct β„‹ using two dimensions, which is defined as: Delaunay triangulation to find a set of connected edges. The edges are then directed according to min 𝑓𝑓(𝒙𝒙) = βˆ’π‘ π‘ π‘™π‘™π‘™π‘™(2π‘₯π‘₯1 βˆ’ 0.5Ο€) βˆ’ 3π‘Žπ‘Žπ‘π‘π‘ π‘ (π‘₯π‘₯2) βˆ’ Definition 6. Definition 7 enables us to find the 0.5π‘₯π‘₯1, π‘₯π‘₯ ∈ Ξ© = [0,9] Γ— [βˆ’2.5, 2.5] minimiser set, which in this case is β„³ = {𝑣𝑣1,𝑣𝑣7, 𝑣𝑣13}. Figure 25 is the resulting structure, A plot of this function with its three local which highlights the domain of s.t. (𝑣𝑣1). minima is shown in Figure 23. Figure 23. The Ursem01 Function. Figure 25. A directed complex β„‹β€”asimplicial approximation for an objective function. The set 𝒫𝒫 contains π‘˜π‘˜ = 15 sampling points from the 2-dimensional Sobol sequence. Figure Increasing the sampling size to π‘˜π‘˜ = 150 and 24 contains a mapping of the objective function repeating the procedure produces the complex in values. Figure 26. Multidisciplinary Design Optimization 18 Oregon Undergraduate Research Journal: McNair Special Issue Casserly design vector and its merit function evaluation between iterations such that the optimization will terminate. Simplicial Homology is my preferred method as it finds approximations to local minima and then uses iterative Nelder-Mead at each of these to find the global minimum. This method is theoretically guaranteed to find the global minimum when using the 'simplicial' sampling method. However, for a merit function as complex as ours (involving trajectory simulation), it is inefficient. The β€œsobol” sampling method will approximate the global minimum Figure 26. A directed complex β„‹β€”asimplicial approximation with an execution time of about 2 hours. I have for an objective function with 150 vertices. also added comments to the code related to setting (black box/design vector) boundaries. This has different minimiser vertices that are better approximations to the local minima, but 3.4. Initial Design Modifications |β„³| is unchanged. This points to the SHGO property: if the number of initial sampling points In coordination with Hayden Reinhold from the is adequate, |β„³| ceases to grow with increasing PSAS airframe team, I have updated the initial π‘˜π‘˜, which provides a heuristic for the number of template.ork OpenRocket file to approximate the sampling points needed to approximately map current design. This involved modifying the local minima of a merit function. component weights and lengths, along with using an approximate thickness to model our isogrid 3.3. User Interface and Documentation plate bulkheads as having uniform density, Improvements demonstrated in Figure 27. In relation to making the workings of the code more understandable, I added the following prior to the code block containing the three main techniques: The optimization approaches are RBFOpt Global Optimization, Iterative Nelder- Mead, and Simplicial Homology Global Optimization using Nelder-Mead at the local minima. RBFOpt produces results in about 2 hours, depending on your machine. The two array arguments passed to the Figure 27. Isogrid plate bulkhead. RbfoptUserBlackBox class define the bounds of The updated initial design fed into the the black box and correspond to minimum and optimizer produced the diagram found in Figure maximum feasible values for the design vector. 28. Iterative Nelder-Mead does take a while; however, in the iterate function in the above code block, you may change the β€œxatol” and β€œfatol” parameters to relax the termination condition. These correspond to the absolute error in the Figure 28. Updated initial design. Multidisciplinary Design Optimization 19 Oregon Undergraduate Research Journal: McNair Special Issue Casserly The main changes were made to the weight/length of Nosecone, Electrical Recovery System (ERS), N2 tank/Reaction Control System (RCS), Avionics/Camera module, Liquid Oxygen (LOX) tank, and fin can. A SHGO simulation with this updated model resulted in global coverage of the design space, as in Figure 29. Figure 31. LV4 Trajectory Information with SHGO Approach. 3.5. Structural Analysis Output Related to the goal of improving the MDO via providing structural analysis output, I added a new notebook that ported relevant code from the structural model notebook. This code contains a structural analysis function that calculates the axial and lateral loads along with bending Figure 29. SHGO coverage of design space with updated design. moments at launch (tip-off), maximum aerodynamic pressure (max Q), and before and We met the majority of our constraints; after engine burnout. The axial forces consist of however, manifesting a 10.6 kN engine poses a friction along the body, parasitic drags related to problem, displayed in Figure 30. each passthru module, and drag coefficient contributions. Lateral load is calculated via summing normal forces, and the bending moment is calculated considering the shear at the top of each component along with the lateral load at the middle. Using the structural plot function in the Display_Information notebook, the main MDO notebook now outputs the graphs seen in Figure 32. Figure 30. SHGO constraint satisfaction with updated design. Our apogee estimate is slightly conservative, so 95.6 km is excellent. Trajectory information indicates a successful launch is possible with the design in Figure 31. Multidisciplinary Design Optimization 20 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Figure 32. Structural analysis information added to MDOoutput. 3.6. Other Modifications In order to make weight reductions easier, Peter McCloud (a scientist affiliated with NASA Aerothermodynamics) suggested that a pie chart with the mass breakdown be added to the MDO output. Figure 33 displays a code snippet and Figure 34 the result. Figure 33. Mass breakdown code added. Figure 34. Mass breakdown pie chart. The only non-trivial component of this is calculating the airframe mass. The volume of the tube was first calculated given the existing airframe length, diameter, and thickness Multidisciplinary Design Optimization 21 Oregon Undergraduate Research Journal: McNair Special Issue Casserly variables. Given the 2700 π‘˜π‘˜π‘”π‘”/π‘šπ‘š3 density of 6061- T6 aluminum, the mass is calculated. The pie chart was constructed with reference to the matplotlib documentation (Basic Pie Chart β€” Figure 37. Launch Parameters with reduced launch rail Matplotlib 3.5.2 Documentation). length. Another issue is related to the engine thrust constraint. It is not feasible for the engine to have The optimal CONS_THRUST parameter in a thrust much greater than 6 kN. Since the thrust this case is somewhere between 5599.6 and constraint is lenient, I experimented with setting 5599.7. Figure 38 displays the results for those the CONS_THRUST parameter to 5.5 kN. Figure two options. 35 contains the results using Simplicial Homology with the simplicial sampling method (guaranteed optimality) and COBYLA local minimization. Figure 35. Constraint satisfaction with 5.5 kN engine thrust constraint. The result is a roughly 9 kN engine with 90 km apogee. Reducing the constraint to 5.4 kN, I Figure 38. Constraint satisfaction to find sweet spot of found that the optimal apogee was a scant 14 km. CONS_THRUST parameter. The sweet spot for the constraint is 5.49698 kN, which produces an apogee of 83 km with an The apogee is dismal for the first and the engine thrust of 9 kN, shown in Figure 36. engine thrust is too high for the second. With a 9.8 m launch rail, along with the CONS_THRUST parameter set to 5599.67, we have an apogee of roughly 83 km and engine power of 9.8 kN. The result is still too powerful of an engine, so we need to experiment with the initial fin parameters to reduce the required thrust. We will investigate the velocity off the launch rail related to the engine thrust. The β€œevent_manager” function shown in Figure 39 contains the desired Figure 36. Constraint satisfaction with optimized engine thrust constraint. parameter and Figure 40 shows the related launch velocity analysis code. Building a launch rail greater than 10 meters is difficult; therefore, I ran the simulation with the LAUNCH_TOWER parameter set to 9.8 m, shown in Figure 37. Figure 39. Velocity off launch rail (rkt.launch_speed parameter). Multidisciplinary Design Optimization 22 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Figure 42. Fin resizing function. Figure 40. Launch velocity analysis code. In the trajectory simulation notebook, I have With the standard fins and unoptimized modified the event_manager, time_step, initial design vector, we have a launch velocity of integration, and trajectory functions to 23.47 m/s and engine thrust of 14 kN. We would accommodate fin staging, in particular the code like to reduce the engine thrust to around 10 kN, seen in Figure 43. which means the optimized design should be < 7 kN. The plan was to drop a rectangular fin can, meaning that only the fin root and tip chord Figure 43. Fin staging condition in event manager function. lengths would be changed. Increasing the fin root and tip chords by 8 m, the velocity off the launch The stage_drop_ECEFβ€”Earth-centred, Earth- rail decreases to 19.3 m/s; however, the engine fixed coordinate system (β€˜Earth-Centred, Earth- thrust remains unchanged at 14 kN. Fixed Coordinate System’)β€”parameter determines the altitude at which we drop the fin 3.7. Fin Staging/Optimization can module. To convert from altitude to ECEF, the function pictured in Figure 44 is used. Fin geometry is defined by the root, tip, sweep angle, semispan, and thickness, as in Figure 41. Figure 44. Altitude to ECEF for entry in trajectory function. The trajectory function now includes parameters related to fin staging, seen in Figure 45. Figure 45. Trajectory function arguments. Figure 46 displays the documentation added Figure 41. Typical Fin Nomenclature (Barrowman Package β€” to the main MDO notebook. Barrowman 0.0.1 Documentation). To simulate dropping a fin can mid-flight, we create a second rocket object with new fin parameters. The idea here is to reduce the required engine thrust by starting out with larger fins, and then the smaller fins will reduce drag in the upper atmosphere. The function pictured in Figure 46. Fin staging simulation instructions. Figure 42 performs fin resizing. Multidisciplinary Design Optimization 23 Oregon Undergraduate Research Journal: McNair Special Issue Casserly To test the code, we drop the ECEF at 50 km and leave the fin parameters unchanged, with zero mass reduction, as seen in Figure 47. Figure 47. Fin staging test arguments. The simulation produces the same output as Figure 51. Displaying fin drop altitude and the Mach number at that point. was produced prior to adding this functionality, indicating that the code is sound. In Figure 52, a plot of the stability margin To find the ideal altitude for dropping the versus time has been added in order to double- larger fin can, we must find the point at which check the fin drop logic and find ideal fin sizing. the smaller fins are passively stable. McCloud This process, shown in Figure 53, involved suggested that this be done when the ratio of the modifications to the Trajectory_Simulation, center of pressure to the center of mass is around Display_Information, and MDO notebooks, as the 2. The stability_margin in the MDO calculates this rocket_plot function required an additional information, shown in Figure 48. parameter β€œstability” for the stability margin list. Figure 48. Center of mass vs. center of pressure metric. Figure 49 contains the code to find the ideal Figure 52. Calibers vs Time plot in MDO Notebook. drop time. Figure 53. Modifications to Display_Information enabling Figure 49. Time after launch to drop the fin can. stability (calibers) vs. time plot. Now, the MDO has a component that The initial fluctuation (especially that of the simulates performance with the smaller fins (post negative section of the graph) does not seem to be dropping) to find the altitude at which they are physical; rather, when the wind is turned on in passively stable, as seen in Figure 50. the trajectory simulation, and before the rocket has left the launch rail, there is an interval during which the wind velocity dominates the rocket’s velocity, meaning that the angle of attack is ~ 45 degrees or more. This interferes with the update [1][5][6] (indexing into multidimensional array) parameter in Fig. 48, which is the center of pressure and the seventh element of the array returned by the aero function in the Figure 50. Finding the ideal drop altitude. Aerodynamics_Model notebook, as in Figure 54. The code shown in Figure 51 prints the altitude at which the fins are dropped, as well as Figure 54. Array returned by the aero function in the Aerodynamics Model. the Mach number at that point. Multidisciplinary Design Optimization 24 Oregon Undergraduate Research Journal: McNair Special Issue Casserly This explanation has been verified by turning by 0.2 m. The result, shown in Figure 57, is a off the wind and viewing the stability plot, seen in minimum just above 2. Figure 55. Figure 57. Calibers vs. time with reduced fin length. On adding the fin root and tip chords to the design vector along with boundaries for these values incentivizing a larger root than tip chord, SHGO produces the design parameters pictured in Figure 58. Figure 55. Stability plot with wind turned off (now strictly positive). A plot of the first five seconds of the stability metric was added, given the fluctuating behaviour, with a line indicating the time at which the rocket has left the launch rail, as in Figure 56. Figure 58. Optimized fin root and tip chords. Figure 56A. Stability metric code. The stability was still quite high for this iteration, so McCloud recommended that the semispan be added to the design vector, shown in Figure 59. Figure 59. Modified design vector variables. Figure 60 displays the boundaries given to SHGO. Figure 60. Boundaries on design vector values given to SHGO. Figure 61A–C contains the results. Figure 56B. Stability metric: first five seconds post-launch. Given that the stability is quite high without fin staging, I reduced the fin root and tip chords Multidisciplinary Design Optimization 25 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Figure 61B. Design vector: optimized fins. Figure 61C. Constraint satisfaction information indicating acceptable apogee with 9 kN engine: optimized fins. Here, we achieve feasible fins and a 95 km apogee with a 9 kN engine. The root and tip chords are relatively small, making for a reasonable altitude (less drag), and the larger span corresponds to greater stability, which contributes to a reduced engine thrust. The minimum stability is still quite high, so we scale down the optimized fin parameters to 80%, reducing the fin area. The result is an apogee of 107 km with a 9.8 kN engine, seen in Figure 62. Figure 61A. Trajectory and stability metric: optimized fins. Multidisciplinary Design Optimization 26 Oregon Undergraduate Research Journal: McNair Special Issue Casserly Figure 62B. Design vector and constraint satisfaction information indicating ideal apogee with 9.8 kN engine: optimized fins scaled down to 80%. This is the ideal design. The 9.8 kN thrust may be realized by using two 5 kN engines. What remains unexplained is the 120 km vehicle drift (calculated by considering the coordinates at launch and apogee, as in Figure 63). Figure 63. Coordinates at launch and apogee. Setting the wind parameter to False in the trajectory function increases drift to 121 km, which means there is either an error in the MDO or the wind does not contribute to the drift. Therefore, I have added a plot of the angle of attack from when the rocket leaves the launch rail to just before engine burnout, pictured in Figure 64. Figure 62A. Trajectory and stability metric: optimized Figure 64. Angle of attack from leaving the launch rail to just fins scaled down to 80%. before engine burnout. Pictured is a non-zero angle of attack; however, it is not significant enough to constitute a 120 km drift. One potential hypothesis is that the drift is due to the combination of the Coriolis acceleration and the wind; however, upon setting these variables to zero and false, respectively, the drift actually increased to 126 km, which means that there must be some computational error in the MDO notebook. Multidisciplinary Design Optimization 27 Oregon Undergraduate Research Journal: McNair Special Issue Casserly 4. Conclusion Design of Aerodynamically Stabilized Free Rockets, Military Handbook. In this work, I set out to improve and extend the β€˜Earth-Centered, Earth-Fixed Coordinate System’. existing open-source Multidisciplinary Design Wikipedia, 6 Mar. 2022. Wikipedia, Optimization (MDO) simulation code. I evaluated https://en.wikipedia.org/w/index.php?title=E different optimization approaches such as arth-centered,_Earth- standard global optimization with RBFOpt, fixed_coordinate_system&oldid=1075493103. Nelder-Mead with different local minimization Endres, Stefan C., et al. β€˜A Simplicial Homology approaches, and Simplicial Homology Global Algorithm for Lipschitz Optimisation’. Journal Optimization. I found the SHGO approach to be of Global Optimization, vol. 72, no. 2, Oct. the most effective and converged on the ideal 2018, pp. 181–217. Springer Link, design vector upon including fin parameters and https://doi.org/10.1007/s10898-018-0645-y. scaling down given stability margin information. Martins, Joaquim R. R. A., and Andrew B. Lambe. The 120 km vehicle drift was caused by some β€˜Multidisciplinary Design Optimization: A computational error in the notebook. This has Survey of Architectures’. AIAA Journal, vol. 51, no. 9, 2013, pp. 2049–75. American Institute been confirmed by running trajectory of Aeronautics and Astronautics, simulations with the design variables returned by https://doi.org/10.2514/1.J051895. the MDO, which indicate that the wind is the Minimize(Method=’Nelder-Mead’) β€” SciPy v1.8.0 cause of the large vehicle drift. Further Manual. https://docs.scipy.org/doc/scipy/ debugging and testing are required. However, the reference/optimize.minimize-neldermead. convergence of that ideal design vector serves as html#optimize-minimize-neldermead. a theoretical guidance to the PSAS engineering Monte, Vaughn. Over-Under Expanded Nozzle - team going forward. Propulsion 1 - Aerospace Notes. https://aero spacenotes.com/propulsion-1/over-under- Acknowledgements expanded-nozzle/. Scipy.Interpolate.Interp1d β€” SciPy v1.8.1 Manual. Thank you to PSAS mentors Andrew Greenberg, https://docs.scipy.org/doc/scipy/reference/ge Cory Gillette, Peter McCloud, and Max Eltzroth nerated/scipy.interpolate.interp1d.html. for their guidance and expertise throughout the duration of the project. Thanks also to Scipy.Optimize.Shgo β€” SciPy v1.8.0 Manual. Christabelle Dragoo and Denise Elder from the https://docs.scipy.org/doc/scipy/reference/ge UO McNair team for their support. nerated/scipy.optimize.shgo.html. Shyy, Wei, et al. β€˜Global Design Optimization for Bibliography Aerodynamics and Rocket Propulsion Components’. Progress in Aerospace Sciences, Barrowman Package β€” Barrowman 0.0.1 vol. 37, no. 1, Jan. 2001, pp. 59–118. Documentation. https://open-aerospace. ScienceDirect, https://doi.org/10.1016/S0376- github.io/barrowman/barrowman.html. 0421(01)00002-1. Basic Pie Chart β€” Matplotlib 3.5.2 Documentation. Yao, Wen, et al. β€˜Review of Uncertainty-Based https://matplotlib.org/stable/gallery/pie_and_ Multidisciplinary Design Optimization polar_charts/pie_features.html. Methods for Aerospace Vehicles’. Progress in Bonmin. 2014. COIN-OR Foundation, 2022. Aerospace Sciences, vol. 47, no. 6, Aug. 2011, GitHub, https://github.com/coin-or/Bonmin. pp. 450–79. ScienceDirect, Coin-or/Rbfopt. 2015. COIN-OR Foundation, 2022. https://doi.org/10.1016/j.paerosci.2011.05.001. GitHub, https://github.com/coin-or/rbfopt. Multidisciplinary Design Optimization 28