This w ork presen ts a mo del of a so da-b ottle w ater ro c k et dev elop ed with NASA's Op enMD A O Dymos optimal con trol m ultidisciplinary framew ork. This is an acces- sible example that is able to highligh t man y of the b ene ts and c hallenges of m ulti- disciplinary optimization and of collo cation metho ds. Optimization results for igh t range and heigh t at ap ogee with resp ect to empt y mass, initial w ater v olume and launc h angle are presen ted. Nomenclature ()0 v alue of () att= 0 _ time deriv ativ e igh t path angle w w ater densit y A out nozzle area D drag g acceleration of gra vit y h altitude k p olytropic constan t m mass m empt y empt y massmw w ater mass p pressure in the ro c k et pa am bien t pressure q dynamic pressure r range S cross sectional area T thrust t time v igh t sp eed v out w ater exit sp eed at the nozzle Vb in ternal v olume of the ro c k et Vw w ater v olume in the ro c k et 1 Intro duction The so da-b ottle w ater ro c k et, g. 1, is a great problem to in tro duce basic ph ysics prin- ciples to studen ts, suc h as Newton's la ws, ideal gas dynamics, Bernoulli's equation and aero dynamics [1]. The design of these ro c k ets for optimal heigh t or range, ho w ev er, is not a trivial problem. It in v olv es the n umerical in tegration of a non-linear set of ordinary dif- feren tial equations (ODE) and couples the disciplines of igh t dynamics and propulsion, represen ted b y the uid dynamics inside the ro c k et b ottle. In spite of that, there are sev eral w ater ro c k et comp etitions w orldwide, so it is in- teresting to nd a solution for this problem. In this w ork, w e optimize a w ater ro c k et 1 CC BY 2.5 a v ailable at https://commons.wikimedia.org/wiki/File:Launchbottle.jpg Figure 1: W ater Ro c k et for range and heigh t at the ap ogee, using design v ariables that are easily mo di able just b efore launc h: the empt y mass, whic h can b e c hanged with the use of ballast, the initial w ater v olume and the launc h angle. W e use this v ery accessible problem to sho w a simple but real w orld application of state of the art m ultidisciplinary optimization framew ork Op enMD A O [2] and its optimal con trol extension, Dymos [3]. This rep ort is organized as follo ws. First, w e presen t the mo del form ulation with its division in comp onen ts and implemen tation in Op enMD A O. W e discuss the constan ts required b y the mo del and our c hoice of optimization algorithm. Then, w e pro ceed to in v estigate the mo del with the use of parameter sw eeps and con tour plots. Finally , w e presen t and analyse the optimization results for b oth range and heigh t at the ap ogee. 2 Problem F o rmulation A natural ob jectiv e function for a w ater ro c k et is the maxim um heigh t ac hiev ed b y the ro c k et during igh t, or the horizon tal distance it tra v els, i.e. its range. The design of a w ater ro c k et is somewhat constrained b y the so da b ottle used as its engine. This means that the v olume a v ailable for w ater and air is xed, the initial launc h pressure is limited b y the b ottle's strength (since the pressure is directly related to the energy a v ailable for the ro c k et, it is easy to see that it should b e as high as p ossible) and the nozzle throat area is also xed. Giv en these man ufacturing constrain ts, the design v ariables w e are left with are the empt y mass (it can b e easily c hanged through adding ballast), the w ater v olume at the launc h, and the launc h angle. With this considerations in mind, a natural form ulation for the w ater ro c k et problem is 2 maximize range or heigh t w.r.t. empt y mass, initial w ater v olume, launc h angle, tra jectory sub ject to igh t dynamics (eqs. (5) to (8)) uid dynamics inside the ro c k et (table 1) 0< initial w ater v olume < v olume of b ottle 0< launc h angle <90 0< empt y mass This problem is single ob jectiv e (although the ob jectiv e can b e toggled), con tin uous, constrained, non-linear and deterministic. 2.1 Mo del There are t w o main approac hes to solving an optimization problems whic h in v olv es time in tegration of a dieren tial equation: direct in tegration (i.e. time marc hing) and col- lo cation metho ds. In comparativ e study , Hw ang and Munster [4] found that high order collo cation metho ds suc h as Legendre-Gauss-Lobatto (LGL) are more robust and e cien t than the alternativ es. F or this w ork, w e use the LGL collo cation metho d to discretize the ordinary dieren tial equations. In order to do this, eac h phase is discretized in to a 21 no de LGL segmen t, and the v alues of the states at the ev en collo cation p oin ts are added as a design v ariable. The resulting appro ximation is dieren tiated with resp ect to time, and compared with the time deriv ativ e pro vided b y the mo del. The dierence b et w een these v alues, called the c ol lo c ation defe cts are added to the problem as c ol lo c ation c onstr aints , whic h the optimizer m ust driv e to zero. F urthermore, the dierence of the states in the end and b eginning of subsequen t phases are also added to the problem as linkage c onstr aints [3]. This approac h is a form of the Sim ultaneous Analysis and Design (SAND) approac h, meaning that the optimizer is also resp onsible for solving the ph ysics of the problem. This approac h leads to a fairly large problem, in this case with 151 design v ariables and 151 functions of in terest (1 ob jectiv e function, 140 collo cation constrain ts, and 10 link age constrain ts). In order to solv e this problem form ulation, w e need a mo del capable of, giv en a tra- jectory and the remaining design v ariables, calculate the ro c k et's igh t dynamic and uid dynamic collo cation and link age defects and extract the maxim um heigh t and range ac hiev ed during igh t. Suc h a mo del can b e divided in to three basic comp onen ts: a w ater engine, resp onsible for mo delling the uid dynamics inside the ro c k et and returning its thrust; the aero dynamics, resp onsible for calculating the atmospheric drag of the ro c k et; and the equations of motion, resp onsible for propagating the ro c k et's tra jectory in time, using Newton's la ws and the forces pro vided b y the other t w o comp onen ts. In order to in tegrate these three basic comp onen ts, some additional in terfacing com- p onen ts are necessary: an atmospheric mo del to pro vide v alues of am bien t pressure for the w ater engine and air densit y to the calculation of the dynamic pressure for the aero- dynamic mo del, and a comp onen t that calculates the instan taneous mass of the ro c k et 3 T able 1: Constitutiv e equations of the w ater engine mo del Comp onen t Equation water_exhaust_speed v out=p 2(ppa)=w (1) water_flow_rate _Vw=v outA out (2) pressure_rate _p=kp_Vw (VbVw)(3) water_thrust T= (wv out)(v outA out) (4) b y summing the w ater mass with the ro c k et's empt y mass. The high lev el la y out of this mo del is sho wn in g. 2. In the next section, the form ulation of eac h comp onen t will b e discussed. All mo dels are analytically dieren tiated with Op enMD A O's framew ork. 2.1.1 W ater engine The w ater engine is mo delled b y assuming that the air expansion in the ro c k et follo ws an adiabatic pro cess and the w ater o w is incompressible and in viscid, i.e. it follo ws Bernoulli's equation. W e also mak e the follo wing simplifying assumptions: 1. The thrust dev elop ed after the w ater is depleted is negligible 2. The area inside the b ottle is m uc h smaller than the nozzle area 3. The inertial forces do not aect the uid dynamics inside the b ottle This simpli ed mo delling can b e found in Prusa [5]. A more rigorous form ulation, whic h drops all these simplifying assumptions can b e found in Wheeler [6], Gommes [7], Barrio- P erotti et al. [8]. Assumption 1 is kno wn to underestimate the ro c k et p erformance, since the air left in the b ottle after it is out of w ater is kno wn to generate appreciable thrust [9]. This simpli ed mo del, ho w ev er, pro duces ph ysically meaningful results as will b e sho wn in this text. There are t w o states in this dynamical mo del, the w ater v olume in the ro c k et Vw and the gauge pressure inside the ro c k et p . The constitutiv e equations are sho wn in table 1 and a diagram sho wing the mo del organization is in g. 3. 2.1.2 Equation of motion comp onent ( eom ) The mo delling of the ro c k et igh t dynamics consist of a simple 2D mo del, describ ed b y the follo wing equations 4 modelatmos kinetic_energy water_enginewater_exhaust_speed water_ﬂow_rate pressure_rate water_thrust mass_adder dynamic_pressure aero eomFigure 2: N2 diagram for the w ater ro c k et mo del modelwater_exhaust_speedrho_w p p_a v_out water_ﬂow_rateA_out v_out Vdot pressure_ratep k V_b V_w Vdot pdot water_thrustrho_w A_out v_out F Figure 3: N2 diagram for the w ater engine group 5 _v=TD mgsin (5) _ =T mvgcos v(6) _h=vsin (7) _r=vcos (8) The mo del used b y this w ork w as implemen ted b y F alc k and Gra y [3]. It is imp ortan t to notice that this mo del has a singularit y in the igh t path angle deriv ativ e when the igh t sp eed is zero. This happ ens b ecause the rotational dynamics are not mo delled. 2.1.3 Auxilia ry comp onents ( aero, dynamic_presure, mass_adder, atmos ) The atmos group calculates the air densit y for use in the dynamic pressure calculation and the am bien t pressure for the thrust calculation follo wing the US 1976 mo del. It consists of sev eral lo okup tables and w as dev elop ed b y F alc k and Gra y [3], using data from https://www.digitaldutch.com/atmoscalc/index.htm . Since the range of altitudes ac hiev ed b y the w ater ro c k et is v ery small ( <100m ), the air densit y and pressure are practically constan t, th us the use of an atmospheric mo del is not necessary . It w as, ho w ev er, included b ecause it w as readily a v ailable within Dymos' framew ork, and this w as in fact easier than setting the v alues to constan ts. The dynamic pressure comp onen t implemen ts the dynamic pressure equation q=1 2v2(9) the aero dynamic mo del calculates drag via D=qSC D (10) and the mass adder calculates the ro c k et's instan taneous mass b y summing the w ater mass with the ro c k ets empt y mass, i.e. m=m empt y +wVw (11) 2.1.4 Phases The igh t of the w ater ro c k et is split in three distinct phases: prop elled ascen t, ballistic ascen t and ballistic descen t. If the simpli cation of no thrust without w ater w ere lifted, there w ould b e an extra air prop elled ascen t phase b et w een the prop elled ascen t and ballistic ascen t phases. Prop elled ascen t is the igh t phase where the ro c k et still has w ater inside, and hence it is pro ducing thrust. The thrust is giv en b y eq. (4) from the w ater engine mo del, and fed in to the igh t dynamic equations eqs. (5) to (8). It starts at launc h and nishes when the w ater is depleted, i.e. Vw= 0 . 6 T able 2: V alues for parameters in the w ater ro c k et mo del P arameter V alue Unit Reference CD 0.3450 - [10] S  1062=4 mm2[10] k 1.2 - [9, 11, 12] A out 222=4 mm2[13] Vb 2L - w 1000 kg=m3- p0 6.5 bar - v0 0.1 m=s - h0 0m - r0 0m - Ballistic ascen t is the igh t phase where the ro c k et is ascending ( >0 ) but pro duces no thrust. This phase b egins at the end of the prop elled ascen t phase and ends at the ap ogee, de ned b y = 0 . Descen t is the phase where the ro c k et is descending without thrust. It b egins at the end of the ballistic ascen t phase and ends with ground impact, i.e. h= 0 . 2.2 Mo del pa rameters The mo del requires a few constan t parameters. The v alues used for those are in table 2. V alues for the b ottle v olume Vb , its cross-sectional area S and the nozzle area A out are determined b y the so da b ottle that mak es the ro c k et primary structure, and th us are not easily mo di able b y the designer. The p olytropic co e cien t k is a function of the moist air c haracteristics inside the ro c k et. The initial sp eed v0 m ust b e set to a v alue higher than zero, otherwise eq. (6) b ecomes singular. This issue arises from the angular dynamics of the ro c k et not b eing mo delled. The drag co e cien t CD is sensitiv e to the aero dynamic design, but can b e optimized b y a single discipline analysis. The initial pressure p0 should b e maximized in order to obtain the maxim um range or heigh t for the ro c k et. It is limited b y the structural prop erties of the b ottle, whic h are mo di able b y the designer, since the b ottle needs to b e a v ailable commercially . Finally , the starting p oin t of the ro c k et is set to the origin. 2.3 Optimization algo rithm and derivatives computation Due to the size of the resulting problem and the presence of man y equalit y constrain ts, it is imp erativ e to use a gradien t based optimization algorithm. W e c ho ose SciPy's SLSQP algorithm, due to its a v ailabilit y and adequate p erformance for mo deratelly large, equalit y constrained problems. This optimizer is a sequen tial quadratic programming optimizer. It should b e noted that since the problem has as man y design v ariables as functions of in terest, it do es not b ene t from the adjoin t metho d for calculating deriv ativ es. Therefore, the total deriv ativ es w ere calculate from the partial deriv ativ es using the direct metho d. 7 3 Investigation on the coupled mo dels T o qualitativ ely ev aluate the mo del, parameters sw eeps with resp ect to empt y mass (50 300g) and initial w ater v olume(0.5-1.5L) are sho wn in gs. 4 and 5, with con tours of nal range, ap ogee heigh t and optim um launc h angle for maxim um range. The range of the sw eeps w ere c hosen to co v er a represen tativ e range of the searc h space. The empt y mass is limited from b elo w b y zero and unlimited from ab o v e, while the initial w ater v olume is limited from b elo w b y zero and from ab o v e b y the b ottle v olume, in this case 2L. By doing a though t exp erimen t, it is easy to see that with zero or maxim um w ater v olume there will b e no thrust generated and therefore b oth range and heigh t will b e zero. The optim um m ust lie b et w een these extremes. Lik ewise, with a v ery high mass the w eigh t will b e higher than the thrust and the ro c k et will not lea v e the ground. The v ery lo w mass case is sligh tly more complicated, but if w e consider that the prop elled phase is v ery short, ev en if the ro c k et gains a lot of sp eed in this phase it will not co v er m uc h distance. A v ery lo w mass means that the motion will b e quic kly halted b y atmospheric drag once the prop elled phase is o v er. Therefore, w e do not exp ect the zero empt y mass case to b e optimal. These engineering considerations are con rmed b y the con tour plots, whic h suggest that the problem is unimo dal as exp ected. This means that the c hoice of a lo cal optimizer suc h as the SciPy's SLSQP is adequate. These plots also sho w that the con tours are smo oth, whic h suggests lo w solv er noise. This is con rmed b y the go o d con v ergence tolerance ac hiev ed, in the order of 108at mo derate initial igh t path angles and 106up to 89. F or launc h angles higher than this, con v ergence is v ery hard b ecause the igh t sp eed at the ap ogee b ecomes v ery lo w th us the igh t path angle dieren tial equation 6 b ecomes ill conditioned. In ligh t of this, the v alue of w as b ounded to 89for the optimization. W e found that the c hoice of transcription is the factor that most impact the optimizer p erformance. This happ ens b ecause the transcription con trols not only the size of the problem b y in tro ducing man y collo cation design v ariables and constrain ts, but also its stiness and its coupling. The size of the problem is determined b y the n um b er of col- lo cation p oin ts in the transcription; the coupling b y the order of the transcription, since ev ery p oin t in the one transcription segmen t is coupled to the others; and the stiness of the problem w as found to b e related to ho w accurate the c hosen transcription can rep- resen t the true solution: if the transcription is to o course the optimizer will ha v e a hard time reducing the collo cation defects. It w as found that high order segmen ts are ideal for ha ving high accuracy in the represen tation of the solution while using few er p oin ts than a lo w order transcription of similar accuracy , and the reduced problem size comp ensated the higher coupling. Ho w ev er, it w as observ ed empirically that this came at the cost of more sensitivit y to initial guesses, although more w ork is required to quan tify these observ ations. 4 Optimization results The optim um design v ariables, optimizer p erformance parameters and other selected pa- rameters of in terest can b e seen in table 3. These results sho w go o d agreemen t with the parameters sw eeps sho wn in section 3 It is in teresting to notice that b oth the optimal 8 50 100 150 200 250 300 Empty mass (g)0.60.81.01.21.4Water volume (L) 54 m 53 m 50 mFigure 4: Heigh t at ap ogee con tours for the w ater ro c k et problem with launc h angle xed to89 50 100 150 200 250 300 Empty mass (g)0.60.81.01.21.4Water volume (L) 85 m 80 m 75 m 45 40 44 46 50 Figure 5: Range con tours in green-y ello w for the w ater ro c k et problem with optimal launc h angle. Launc h angle con tours are in grey . 9 T able 3: Optimization results Maxim um range Maxim um heigh t Empt y mass 0:1861 kg 0:1429 kg W ater v olume 0:9045 L 0:8704 L Launc h angle 45:029289:0000 Fligh t angle at end of prop elled phase 38:340288:8437 Maxim um range 86:4126 m 3:2027 m Maxim um heigh t 23:5110 m 54:6032 m Maxim um v elo cit y 42:4133 m/s 47:2591 m/s Iterations 105 419 F unction ev aluations 584 3988 Gradien t ev aluations 105 419 T olerance 11085105 empt y mass and w ater v olume are larger for the range optimization than for the heigh t optimization. This suggests that the problem is truly m ultidisciplinary: if the w ater en- gine could b e optimized b y itself, for example b y maximizing its impulse, the w ater v olume w ould b e the same for b oth cases. Instead, w e see that this is not the case b ecause the propulsion and igh t dynamics disciplines are couple through mass. Due to the singularit y at = 90, the con v ergence for the heigh t optimization is m uc h w orse than for the range optimization, b oth in n um b er of iterations and function ev aluations, ev en with a lo w er tolerance set for the former case. It is in teresting to notice that for the maxim um range optimization, the launc h angle is not smaller than 45as w ould b e exp ect for a ballistic igh t with drag. This happ ens b ecause at the small sp eeds of the b eginning of the prop elled ascen t phase, the igh t path angle c hanges v ery fast as sho wn in g. 9. Because of this eect, at the end of the prop elled phase the igh t path angle is smaller than 45, whic h agrees with the in tuition from the ballistic case. This high non-linearit y in the igh t path angle also explains the need for a high order transcription. The tra jectory of the ro c k et optimized for range can b e seen in g. 8, and the tra jectory and selected states for the prop elled phase for the ro c k et optimized for heigh t can b e seen in gs. 6 and 7 resp ectiv ely . In the latter case, notice ho w fast the igh t path angles c hange at the ap ogee. This due to the drastic reduction in sp eed at the ap ogee that is c haracteristic of launc h angles close to v ertical. Once again, these high non-linearities justify the c hoice of a high order transcription 5 Conclusion In this w ork a mo del of a w ater ro c k et implemen t using the state of the art Op enMD A O Dymos framew ork w as presen ted and optimized for b oth range and heigh t at the ap ogee. This is a simple mo del that is able to illustrate w ell v arious concepts and c hallenges of m ultidisciplinary optimization and ph ysical principles. 10 0123r (m) Propelled Ascent Ballistic Ascent Descent 02040h (m) 010203040v (m/s) 0 2 4 6 t (s)90 45 04590 (deg) (a) igh t dynamics states 0.02.5 r (m)01020304050h (m) (b) tra jectory Figure 6: Optimized tra jectory for maxim um heigh t at ap ogee 11 0246p (bar) 0.00.51.01.52.0Vw (L) 0200400T (N) 02040v (m/s) 0.00 0.02 0.04 0.06 0.08 t (s)88.84088.88588.93088.975 (deg) Figure 7: Thrust and selected states during prop elled ascen t phase optimized for heigh t. 12 0255075r (m) Propelled Ascent Ballistic Ascent Descent 01020h (m) 02040v (m/s) 0 1 2 3 4 t (s)60 45 30 15 0153045 (deg) (a) igh t dynamics states 0 20 40 60 80 r (m)01020h (m) (b) tra jectory Figure 8: T ra jectory for maxim um range ro c k et 13 0246p (bar) 0.00.51.01.52.0Vw (L) 0200400T (N) 02040v (m/s) 0.00 0.02 0.04 0.06 0.08 t (s)394245 (deg) Figure 9: Thrust and selected states during prop elled ascen t phase optimized for range. 14 It w as found that the problem is truly m ultidisciplinary in nature. It w as also sho wn ho w the c hoice of optimization function aects the optim um and the con v ergence qualities of the problem. F urthermore, the usefulness of high order transcriptions for dynamic mo dels w as highligh ted. In spite of ha ving a simple mo delling, the w ater ro c k et can generate a fairly large optimization problem if collo cations metho ds are used. In order to ac hiev e acceptable accuracy , the problem m ust ha v e h undreds of design v ariables and functions. Therefore, a gradien t based optimization metho d with high qualit y deriv ativ es is needed. This w ork highligh ts ho w an organized framew ork can ease the dev elopmen t eort of a mo del with suc h requiremen ts. Finally , w e exp ect that this accessible example of the use of Op enMD A O and Dymos will help to increase the reac h of these tec hnologies, an w e exp ect this example to b e incorp orated in to Dymos' source co de in the near future.