Concurrent Trajectory Optimization and Aircraft Design for the Air Cargo Challenge Competition

: A coupled aerostructural aircraft design and trajectory optimization framework is developed for the Air Cargo Challenge competition to maximize the expected score based on cargo carried, altitude achieved and distance traveled. Its modular architecture makes it easily adaptable to any problem where the performance depends not only on the design of the aircraft but also on its ﬂight trajectory. It is based on OpenAeroStruct, an aerostructural solver that uses analytic derivatives for efﬁcient gradient-based optimization. A trajectory optimization module using a collocation method is coupled with the option of using b-splines to increase computational efﬁciency together with an experimentally-based power decay model that accurately determines the aircraft propulsive response to control input depending on the battery discharge level. The optimization problem totaled 206 variables and 283 constraints and was solved in less than 7 h on a standard computer with 12% reduction when using b-splines for trajectory control variables. The results revealed the need to consider the multi-objective total score to account for the different score components and highlighted the importance of the payload level and chosen trajectory. The wing area should be increased within allowable limits to maximize payload capacity, climb to maximum target height should be the focus of the ﬁrst 60 s of ﬂight and full throttle should be avoided in cruise to reduce losses and extend ﬂight distance. The framework proved to be a valuable tool for students to easily obtain guidelines for both the model aircraft design and control to maximize the competition score.


Introduction
The Air Cargo Challenge (ACC) is an international competition targeted at aeronautical engineering students of higher education institutes created by Instituto Superior Técnico in 2003 and since run every two years with the sponsorship of the European Association of Aerospace Students (EUROAVIA).
The participating teams face the complex engineering task of designing, building and flying a radio-controlled aircraft model to carry the highest possible cargo while satisfying the competition rules. The winning team is awarded the right to organize the next event and adjust its regulations.
The ACC 2022 edition is organized by AkaModell Munich, having the rules [1] changed significantly from previous editions. While before the main focus was on maximum cargo, it has evolved to a real-life scenario of transporting medical emergency supplies, in particular blood bags, from one point to another. Consequently, not only cargo but also trajectory are important for the performance score.
The expected flight profile for the ACC 2022 competition is represented in Figure 1, with six segments identified. The goals of the aircraft are to (i) transport as many bags as possible; (ii) as far as possible (free flight loops in segment [3][4] in two minutes, while gaining altitude during 60 s to avoid obstacles (segment 2-3) after take-off (segment 1-2); (iii) as quickly as possible, before landing (segment 4-1). The score of the competition reflects how these three (i)-(iii) goals are met compared to the other teams. In addition, the aircraft design should privilege little space for transport, easy assembly and short take-off.
The ACC 2022 rules impose a fixed-wing, radio-controlled aircraft configuration, where the propulsion system must use a prescribed electric motor and propeller, thus limiting maximum thrust, to create a fair competition. The ACC prototype aircraft designed and built for proof-of-concept in 2021 by the Olissipo Air Team of Instituto Superior Técnico is shown in Figure 2. Notice the adopted tractor and V-tail configuration for low structural weight. Since ACC is a demanding multidisciplinary design problem, no single work has been published covering all relevant coupled disciplines. Among the few published works, excluding academic bachelor and master theses, the most significant scientific documents reveal emphasis on singular disciplines, such as aerodynamic design [2,3], mass prediction models [4], parametric geometric definition [5] or flight dynamics and control [6].
Given the mission goals of the ACC 2022 competition, not only the aircraft design but also the trajectory flown are of utmost importance.
On one hand, the aircraft design should take into account both aerodynamics and structures, aiming to maximize lift for take-off and speed for cruise while actively reducing the overall weight and drag. The aerodynamic and structural disciplines need to be coupled in the form of aerostructural design, as proved necessary in most aircraft design problems from simple configurations to vastly complex, such as multi-operating point conditions [7], morphing wings [8], transonic wings [9] and box-wing configurations [10]. The design process focuses on the conceptual phase, to provide an overview of the aircraft shape, size, weight and performance, given the specified operating conditions [11].
On the other hand, the trajectory followed in both the climb and cruise flight segments are also relevant for the score in the competition, where they are sought to be the fastest and the furthest, respectively. Consequently, the approach needs to account for trajectory optimization. Typically, aircraft trajectory optimizations minimize cruise fuel consumption but they can also be used to improve data collection [12], search efficiency [13], energy efficiency [14], fuel consumption [15] or environmental impact [16]. This work focuses on a single multi-segment trajectory optimization considering take-off, climb and cruise portions while respecting the flight path continuity.
In order to search for the best scoring solution for the ACC 2022 aircraft to perform the specific mission, a concurrent design approach is created, where the aircraft design and trajectory optimization are merged using Multidisciplinary Analysis and Design Optimization (MDAO) techniques [17]. These techniques have been extensively proved successful in aerospace, typically to solve highly complex problems, such as morphing wings [18], tow-steered composite wing structures [19] and tilt-wing take-off trajectory and control scheduling [20]. Recent applications include the design of supersonic aircraft, considering the disciplines of aerodynamic, heat transfer and propulsion, together with a pathdependent formulation [21], in a similar thought process of the one adopted here, where both aircraft design and trajectory are simultaneously considered for the optimal solution.
The purpose of this work is then to develop an MDAO design framework that couples the aerodynamics, structures, propulsion and control disciplines, capable of searching for the optimal conceptual design and trajectory for an aircraft competing at the ACC 2022. This work builds over and offers a considerable extension to OpenAeroStruct [22], an open-source low-fidelity aerostructural analysis and optimization tool. The resulting framework is meant to be modular to easily allow new performance metrics, constraints and design variables to be added, thus not only providing a significant design efficiency improvement to the Olissipo Air Team for future ACC competitions but also becoming an open-source software tailored to aerospace engineering students.

Disciplinary Models
Being the design focus at the conceptual level, the developed optimization framework includes low-fidelity numerical disciplinary models, which provide acceptable physical solutions and satisfactory design trends at an affordable computational cost. Four disciplines are included-aerodynamics, structures, propulsion and trajectory optimization.

Aerodynamics
The Vortex Latice Method (VLM) [23] is the numerical method used to study the aircraft aerodynamics. It is particularly suited for lifting surfaces made of thin airfoils and operating at low angles of attack, which is within the scope of the work. It can handle a variety of wing planforms, even low aspect ratio, swept or delta wings, by modeling the wing as a combination spanwise and chordwise panels. Each panel k has a horseshoe vortex, as illustrated in Figure 3, consisting of a finite bound vortex in the spanwise direction and two semi-infinite trailing vortices that extend into the freestream direction. Each vortex filament segment dl induces a velocity field v at an arbitrary point P distancing r, according to the Biot-Savart law, where Γ k is the vortex circulation strength.
Integrating Equation (1) over the finite and semi-infinite straight vortex filaments, adding the contribution of all panels and imposing flow tangency conditions at the control point in each panel, leads to a linear system of equations, where A ik is the aerodynamic influence coefficient matrix, v ∞ is the freestream velocity vector and n i is the normal vector of panel i. Solving Equation (2) for the circulation, the aerodynamic force acting in each panel can be calculated as where v k is the induced velocity vector at the center of the bound vortex, and l k is the bound vortex vector. Lift and induced drag result from adding these panel forces decomposed into the normal and streamwise directions, respectively. To improve the drag estimation model, a skin friction drag is calculated using flat-plate-based estimates [11]. It should be noted that only the aircraft lifting surfaces (main wing and tail) are modeled using VLM. The influence of all other parts are represented by prescribed additional contributions to total lift and drag.

Structures
The Finite Element Model (FEM) is the numerical method used to study the aircraft structure. Only the main load bearing components are modeled, namely the wing and tail spars. A failure analysis is used to size these spars, ultimately reducing their mass while preventing failing under operating loads.
A beam element consisting of a combination of truss, beam and torsional elements, which simultaneously carry axial, bending and torsional loads, is used. It comprises two nodes that can translate and rotate in every direction, totaling 12 degrees-of-freedom (DOF), as presented in Figure 4. The resulting local stiffness matrix of size 12 × 12 for each element depends on the material properties (Young's modulus E and shear modulus G), spar geometry (element cross-sectional area A and length L) and inertia (polar moment of inertia J and second moments of area I). The global stiffness matrix K ij is assembled by applying transformation matrices to the local stiffness matrices. The (aerodynamic) applied loads are transferred to the nodes i in the form of a force vector f i . The static equilibrium condition then leads to a linear system of equations, The vector of displacements u j , determined by solving Equation (4) is then used to calculate the stresses in each element according to the constitutive laws.
A wingbox spar configuration have been adopted in this study due to their higher load-bearing capability for the same weight [24]. The FEM model requires the skin and spar thicknesses, t skin and t spar , to evaluate the sectional properties, as visualized in Figure 5. The details about the wingbox properties calculation and the stress analysis can be found in [25].
The aircraft maximum take-off weight is given by where W bat is the battery weight and W payload is the total weight of the blood bags. The empty weight W empty is decomposed in wingless and tailless structural weight W 0 , wing weight W wing and tail weight W tail . These last two result from the aerostructural sizing of the corresponding spars. As such, only W 0 is prescribed, being all other contributions to MTOW variable in the design process.

Propulsion
The ACC 2022 regulations [1] impose an electric propulsive system, with a prescribed motor and propeller. The motor is the Axi 2826/10 Gold line v2, a high efficiency brushless motor with reduced maintenance and high power-to-weight ratio, that might be modeled as a transfer function from input electric power to output mechanical power with a loss factor. The propeller is the Aeronaut CAMcarbon Light 10 × 6, a 10-inch diameter by 6-inch pitch propeller, which transforms mechanical power to thrust. Propellers are often modeled with the blade element momentum theory, that can be enhanced by adding a loss factor k to account for the propeller profile drag, blade tip and hub losses, and number of blades [20]. The power is controlled by an electronic speed controller (ESC), freely chosen, using a radio control signal sent by the operator. The ESC can be modeled as a linear relation with a loss factor between operator input δ t and the output power. A lithium-polymer (LiPo) battery must be used, for which the discharge curve during operation must be estimated to determine the maximum available power as function of the remaining energy [26] and rate of discharge.
Rather than using simple mathematical models, the propulsive model is built based on wind tunnel experimental data using the real parts because it is believed to better represent the system. The experiment comprised a force scale apparatus to measure electric voltage and current and the propeller thrust for different wind tunnel speeds. Several experiments produced data for curve fitting of power and thrust (outputs) with airspeed (operating state) and throttle (control input) setting. Figure 6 presents the measured thrust and electric power results using a constant power supply with the speed controller set at full throttle for varying airspeed. It includes the quadratic fitting curve for the experimental data and the theoretical curve with a loss factor k of 1.5. As expected, thrust is maximum at v ∞ = 0 (static) and decreases non-linearly with airspeed due to the reduction in the propeller induced momentum on the flow. The results demonstrate that the momentum theory does not accurately describe the propulsive system, particularly for airspeeds above 10 m s −1 , probably because the introduced linear loss factor k that corrects the initial thrust but fails to capture the existing non-linear loss phenomena. As such, the experimental quadratic fit curve presented in Figure 6a with very good fit R 2 = 0.9973 will be used to describe the available maximum thrust. The quadratic fitting for the electric power consumption in Figure 6b will be used to identify the operating point in the thrust-airspeed curve.
A second experiment studied the dependence of thrust and power with the throttle input from the controller (ESC) still using a constant power supply. Due to the ESC setup, it does not power the motor below 20% throttle setting and it saturates above 70%, so only this active range is presented in Figure 7. The dependency of power on throttle remains approximately the same regardless of airspeed, with differences resulting from the maximum power reduction observed earlier on the full throttle results with increasing airspeed. Figure 8 exhibits the good fitting achieved for the experimental thrust-speed data obtained using the power consumption at each throttle setting and a linear approximation to estimate the thrust with the maximum power for a certain airspeed as where δ t and T(v) are the throttle setting and the thrust experimental data, respectively. To complete the characterization of the propulsive system, a model for the LiPo battery discharge curve was derived based on typical LiPo battery properties and characteristic three-staged behavior. There are three regions in a discharge curve: (1) a rapid exponential decay in voltage, typically from 4.2 V to 3.7 V, during the first 20% capacity; (2) linear slowly decaying voltage to 3.2 V when 80% capacity is used; and (3) when the voltage starts to rapidly fall and should be not used. This led to the generic single-cell LiPo battery voltage model as function of the consumed capacity represented in Figure 9. The voltage at the initial exponential region is defined as with V i , V f , C % , C % i , C % f and A = 10 are the initial and final voltages, the depleted, initial and final capacities and the exponential factor, respectively. The battery discharge model allows for the implemented propulsion model to account for the available maximum power decrease as a function of its usage in flight. The propulsive model uses as input the throttle setting δ t , which determines the current I drawn by the propulsion system using the current-throttle curve ( Figure 10) obtained from the power-throttle curve (Figure 7) divided by the average voltage of the experiments (11.5 V). Because of the ESC characteristics, it was transformed into a linear function to simplify the problem. By keeping track of the elapsed flight time the battery capacity consumed is calculated, which is then used to determine its voltage V from using the battery voltage-capacity curve ( Figure 9). The electric power is then computed as P = V I. Afterwards, the maximum power P max for the given speed is determined using the power-speed curve (Figure 6a), and the full throttle thrust is obtained using the thrust-speed curve (Figure 6b). The output thrust given the actual throttle setting (and battery state) is finally obtained using Equation (6).
The discretization of time, and consequent computation of current, capacity consumed, voltage, power, thrust and energy, corresponds to the number of collocation points used in the trajectory.

Trajectory Optimization Method
To find the best path to be followed by the aircraft to maximize the competition score, a trajectory optimization method must be employed. Following the appropriate equations of motion and the boundary and initial conditions of the problem, the input control variables determine the change in the state variables that define the position and velocity of the aircraft as function of time [27].
The methods for solving trajectory optimization problems can be classified as direct or indirect [28]. In the context of concurrent design and trajectory optimization, the usage of an indirect method would imply an optimization architecture with multi-level optimizers in which the trajectory optimization would be nested inside the design optimizer, thus implying considerable computational cost [29]. As such, this work uses a direct collocation method (also known as a direct transcription method), which is an implicit direct trajectory optimization based methodology that uses the control and state variables as optimization variables and the numerically integrated differential equations of motion as constraints. This option makes it possible to pose the complete problem using a single-level optimizer.
Since the trajectory planned for the ACC aircraft is mostly forward flight with high radius turns, only the longitudinal equations of motion are considered, resulting in the state vector {x, v x , z, v z }, where x and z are the horizontal and vertical position, respectively, and v x and v z the corresponding velocities. The collocation method divides all variables defined for the trajectory, which are functions of time t, in finite intervals of time [t k , t k+1 ] and describes them in a polynomial manner of a specific degree. In each interval, the dynamics described for the problem must be satisfied, so a numerical integration method is used to calculate the imposed constraint as Newton's 2nd law, F = mv, governs the aircraft motion, leading to the set of defect functions to be constrained to zero by the optimizer where F x k and F z k are the force control variables and m is the aircraft mass. Ultimately, the optimizer will collocate the trajectory's control (F x k , F z k ) and state variables ( into the values that optimize the objective function. Typical approaches to aircraft trajectory optimization are multi-point based, meaning the trajectory is divided into several different flight segments, with each segment being individually optimized to achieve a more robust design [31]. However, they do not consider the different aircraft behavior in each segment, which is the ultimate purpose of this work. Using more collocation points in the direct collocation method results in more accurate solutions but also increases the computational cost. To mitigate that, 3rd-order b-spline function approximations are adopted to increase computational efficiency while maintaining good solution convergence and accuracy. The b-splines create fits to the system dynamics criteria with potentially fewer control points (used as optimization variables) than collocation points (used as trajectory points).
The benefits of the aforementioned approach were first assessed using the 1-D Bang-Bang problem [32], whose goal is to minimize the time t f needed for a unit mass body to travel a prescribed distance x f , starting from rest, by applying a force F. The dynamics of the problem were satisfied by enforcing the defects, Equations (10) and (12), to zero. The force was assumed bounded in the interval [−2, 1] and x f = 300 m, which led to the exact analytical solution t f = 30 s.
First, the optimal solution was achieved using the original direct collocation method with trapezoidal implicit integration for an increasing number of collocation points. That increase in collocation points led to significant gain in accuracy until a certain threshold (about 40 in this case, resulting in t f = 29.9 s after a total of 126 function evaluations by the optimizer). Next, b-splines were used to describe the variables of the problem (x(t), v(t) and F(t)) given a set of control points, thus allowing a different number of control points and collocation points. The use of b-splines was tested in an incremental fashion: first applying it only to the control force using 8, 16 and 30 control points (Figure 11a), then using that spline with 16 control points and adding the b-splines for states position and velocity ( Figure 11b). As expected, the b-spline modified approaches also experienced significant gain in accuracy by adding more control points, but always fewer than collocation for net computational cost gain.
Comparing the computational cost reduction in Table 1, the force b-spline version led to savings of about 84% in the number of function evaluations in the optimization process, while the difference of the objective function (time of the trajectory) is less than 1% compared to the original method. When both the control and state variables are controlled by splines, the gains are still considerable, varying between 22 and 48%. Despite the computational efficiency gains being problem specific, the b-spline approach was implemented for both the control and state variables in the trajectory optimization module of the concurrent design framework and demonstrated in the final optimization problem (Section 3.4).

Competition Related Models
The specificity of the ACC competition implies the development of four particular models to be implemented into the design framework: the cargo bay parasite drag prediction model, the geometric rhombus box model, the trim and stability model, and the take-off distance.

Cargo Bay Drag Model
Since OpenAeroStruct [22] only models the lifting surfaces, additional drag contributions must be estimated for all other parts, particularly for the bulky fuselage where the cargo bay is located. The contributions of other parts to the aerodynamic characteristics were lumped as an additional parasite drag coefficient based on data collected by the ACC design team from previous prototypes.
The cargo bay of the aircraft is designed to carry simulated blood bags as payload [1], whose arrangement was assumed to be n bags stacked vertically by m bags lined-up in tandem, as schematically described in Figure 12. The model estimates not only the weight and size of the cargo but also the added drag based on the number of blood bags, based on 300 g blood bags of size product 150 mm × 100 mm × 30 mm [1]. The fuselage parasite drag is calculated as [33] where S re f , C f , F f orm , F inter f and S wet are the reference area, surface skin friction coefficient, form factor, interference factor and surface wetted area, respectively. While a unitary interference factor is considered, the form factor can be estimated for various cross section fuselages [34]. For simplicity, the fuselage was approximated to an ellipsoid body of revolution, being to the form factor given by [35] where the fineness ratio is f = 1/ 4A max π . For fineness ratios in the range 6-12, it correlates well with other possible formulas for the form factor. Since the fuselage is directly behind the propeller wake, a fully turbulent flow is assumed, being the skin friction coefficient given by [36] where Re is the local Reynolds number. The fuselage wetted area depends on the m × n scheme chosen and it is increased lengthwise by the electronics (200 mm) and tail joint (150 mm). The estimated fuselage parasite drag for different arrangements of bags is shown in Figure 13. Fewer bags tend to benefit smaller m arrangements, while more bags tend to benefit higher m arrangements, leading to desired higher length to cross section ratios. The model selects the lower parasite drag configuration for the number of bags being evaluated.

Rhombus Box Model
The ACC 2022 edition [1] limits the size of the assembled aircraft model to fit inside a rhombus box of fixed 1500 mm sides and free inner angle β, as illustrated in Figure 14. A trade-off between high wing span for aerodynamic efficiency and long tail arm for stability and control must be found. The wing tip leading edge was assumed to coincide with the box side, with the wing positioned in the box's diagonal. Moreover, the wing and tail chords must fit inside the box and the propulsion system is ahead of the wing root leading edge. Two geometric gaps are then defined for the wing, where l r , b w and c w tip are the rhombus box lateral side length, wing span and wing tip chord length, respectively, and two geometric gaps for the tail, where l t , c w , c t , c t root , Λ w , c t tip and Λ t represent the tail length, wing mean chord, tail mean chord, tail root chord, wing sweep angle, tail tip chord and tail sweep angle, respectively. These four expressions are used as non-negative inequality constraints.

Trim and Stability Model
The longitudinal static stability of the aircraft can be determined from the static margin. However, numerical studies using XFLR5 [37] showed that similar aircraft exhibit negative C m α slopes, indicating a positive static margin of 10-15% so this stability condition was not explicitly used but only the longitudinal aircraft trim condition at each trajectory collocation point, representing the equilibrium of horizontal and vertical forces and pitching moment.
Since the aircraft has a V-tail design, the vertical fin's area is determined by considering the projected planform of the tail area on the vertical plane S VT , which is estimated from the vertical tail volume coefficient V VT , where l VT is the distance between the wing and tail quarter chords (considering mean geometric chord).

Take-Off Distance
The take-off distance is calculated assuming a constant acceleration [33], with average speed v avg = v LOF / √ 2, being the lift-off speed given by where W, ρ, S and C L max are the aircraft weight, air density, reference area (wing area in this case) and take-off lift coefficient, respectively. Lift L and drag D are obtained from the aerodynamic coefficients estimated for the Olissipo Air Team 2021 prototype at take-off using XFLR5 and the cargo bay drag model, leading to C L = 0.3 and C D = 0.02. Both weight W and thrust T are variables in the design. The friction force is expressed as F wheel = µ track W and, since the ACC competition runway is a grass track (small dry turf), the friction coefficient µ track = 0.07 is used [33].

Flight Segments Integration and Competition Scoring
The three flight segments-take-off, climb and cruise-are handled individually to account for the specific design variables, constraints and overall performance metrics [1].

Take-Off Segment
The take-off segment starts when the electric motor is activated and ends when the aircraft takes-off.
Outside the trajectory model, the take-off distance and speed are evaluated and the payload score computed as where PC team and PC max represent the aircraft payload capacity of the team and the maximum payload achieved overall in the competition, respectively. Take-off must be achieved in less than 60 m for a valid flight.

Climb Segment
The climb segment employs the trajectory collocation method described in Section 2.1.4. The defect constraints that represent the equations of motion in the vertical and horizontal directions must be respected, which imply computing the forces at each collocation point. Figure 15 shows the applied forces in this segment. The weight W is given by the payload and the structural weight, the lift L and the drag D are calculated by the aerodynamic model and the thrust T is given by the propulsion model as a function of time and speed. The climb angle is set as a design variable to allow the optimizer to choose the best flight trajectory.
The climb starts the trajectory and ends when one collocation point reaches at least 60 s of flight time, with the corresponding altitude used to calculate the climb score, where PS altitude,team and PS altitude,max represent the altitude score of the team and the maximum score altitude achieved overall in the competition, respectively. The altitude score is given by being h 60s the altitude achieved after 60 s of climb, which is a concave function with a maximum at 100 m.

Cruise Segment
The cruise segment is handled similarly to the climb segment, but the flight path design angle is set to zero (or negative for descent).
It starts in the last climb collocation point and ends after 120 s, when the total distance traveled is used to calculate the distance score, where D team and D max are the distance traveled by the team and the maximum distance achieved overall in the competition in 120 s of flight time.

Final Score
The ACC 2022 overall performance metric is evaluated as where P to_bonus is a take-off bonus (0.1 if take-off distance is less than 40 m, zero otherwise).
The winning team will be the one that balances each individual metric, implying a trade-off in the design and trajectory. The best aircraft solution is expected to be neither the lightest nor the fastest.

Problem Constraints
Besides the box size restriction in the ACC regulations, the design space of the problem is limited by many other constraints that emerge from aerodynamic, structural, energy, payload, trim and stability, and trajectory.

Box Size Constraints
The competition severely constrains the aircraft size through the rhombus box fitting regulations (Section 2.2.2). These are expressed by the set of inequality constraints The tail length l t is a design variable that not only is used in the last two inequalities but also sets the relative position of the wing and tail discretization meshes at their root as an equality constraint, Other dimensions limited by the transport box are constrained by the design variable bounds.

Aerodynamic Constraints
To avoid lift distributions that might lead to stall, a limitation is imposed on the maximum section lift (2D) coefficients.
Using the wing airfoil designed by the Olissipo team, a XFLR5 simulation at the expected operating conditions (Re ≈ 0.35 × 10 6 ) produced the results shown in Figure 16. Based on these results for the wing and adopting an NACA0009 airfoil for the tail surface, led to the aerodynamic constraints

Structural Constraints
For the stress failure calculations, the equivalent Von Mises stress σ V M is computed for each beam element of the wing and tail spars and compared to an admissible stress (material yield stress σ Y over a safety factor SF) in the form of an inequality constraint.
For each spar, wing and tail, the Kreisselmeier-Steinhauser (KS) function [38] is used to aggregates the several stress constraints into a single one, thus significantly reducing the optimization cost, KS f ailure w < 0 , KS f ailure t < 0 .
A carbon fiber laminate material was chosen for both the wing box and the tubular tail spar, with mechanical properties: Young's modulus E = 70 GPa, shear modulus G = 5 GPa, yield stress σ y = 600 MPa and density ρ = 1400 kg m −3 [39].

Energy Constraints
Using the propulsive model described in Section 2.1.3, the optimizer must ensure the aircraft produces enough power and has enough energy for the propulsion needs.
Firstly, the consumed energy E consumed can not exceed the maximum available energy E max , which depends on the battery mass W bat , specific energy e bat and usable energy factor f usable , E max = W bat e bat f usable > E consumed .
Secondly, the consumed battery capacity C consumed can not exceed the maximum capacity C max , where C bat is the nominal battery capacity. Finally, the nominal battery capacity and its mass are related by 43.2 C bat − W bat e bat = 0 , with C bat expressed in mA h and assuming a battery with 12 V nominal voltage.

Payload Constraints
To comply with the ACC regulations, the payload must not affect the stability by not affecting the aircraft center of gravity (CG). The framework enforces the longitudinal position of the payload CG payload as within a tolerance of 10 −2 , where CG empty is the prescribed empty aircraft CG.

Trim and Stability Constraints
Besides the force balance enforced by the defect constraints in the trajectory model, the pitch trim condition is aldo enforced at each trajectory collocation point k, which is achieved by the optimizer mainly by varying the stabilator angle δ stab . The V-tail vertical stabilizer volume (Equation (21)) is bounded using historical data [33], 0.04 < V VT < 0.09 .

Trajectory Constraints
To ensure flight continuity from take-off to climb segments, that are handled in separate modules, the longitudinal aircraft speed must match, with a chosen tolerance of 10 −2 , whereẋ to andẋ i are the final take-off speed and the initial flight speed, respectively. Notice that continuity between climb and cruise segments is automatically guaranteed by the trajectory module that handles both. A climb helping mechanism is implemented where the maximum possible achievable altitude is registered 60 s after take-off, to allow for the aircraft to cruise with negative path angle, which results in better aircraft's weight contribution to speed. Additionally, the trajectory must comply with the drones flight ceiling defined by European regulations, z ≤ 120 m .

Aircraft and Trajectory Parameters
Following the description of the analysis models implemented in the framework, the complete list of parameters, including both aircraft and trajectory related, are summarized in Tables 2 and 3, respectively. These include all design variables, constraints and constants, and corresponding bounds or fixed values.
The fixed design parameters are defined by the user during problem setup and kept constant during the optimization.
The number of control points of the aircraft variables in Table 2 (when applicable) was chosen to allow enough design freedom for the optimizer to enhance aircraft performance but also kept realistic in terms of model aircraft building difficulty. This led to either constant (n = 1), linear (n = 2) or quadratic (n = 3) functions.
The number of control points for the trajectory variables in Table 3 was defined based on preliminary convergence studies.

MDO Framework Architecture
Having presented in Sections 2.3.2 and 2.3.3 all constraints and design variables, the non-linear programming problem may be formulated in standard form as where x is the vector of all design variables, C is the vector of constraints. The objective function F(x) is defined as the inverse of the competition score function (Equation (28)) so that the desired maximization of the latter is obtained. Each collocation point has 6 design variables (4 state and 2 control variables) and 4 constraints for the trajectory and several other to limit some metrics (stability, C l ) in each flight stage. This translates to a problem with 6 × n + 26 design variables and 9 × n + 13 constraints, where n is the number of control points chosen for the trajectory control.
Every equality constraint is handled as two bounded inequalities with a certain tolerance, as indicated in Tables 2 and 3, to improve the optimizer convergence performance.
The concurrent design approach adopted is a clear example of Multidisciplinary Design Optimization (MDO). While common MDO problems in aircraft design tackle aerostructural phenomena [40,41], this work extends that to also include trajectory optimization.
The architecture of the MDO framework is responsible for the management of the interdependence of the disciplines described in Section 2.1 that need to be considered simultaneously in the design process [42]. A Multidisciplinary Feasible (MDF) monolithic architecture is used, which leads to relatively smaller optimization problems in terms of number of variables and constraints, and relies on a single-level optimizer. As a drawback, it requires a coupled Multidisciplinary Analysis (MDA) at every optimization step, whose convergence must be assured by an iterative method such as a Gauss-Seidel fixed point iteration or a Newton-based method. This process ensures that, after each optimization iteration, the solution is multidisciplinary feasible, even if the optimizer ends prematurely, though it might not satisfy the constraints. The architecture of the developed concurrent design framework is presented in Figure 17 using an XDSM diagram [43]. The optimization framework is based on the OpenMDAO software [44], an opensource framework for efficient multidisciplinary optimization. It is highly scalable and facilitates the integration of different disciplinary analysis modules. By defining the new computation routines together with their corresponding analytic derivatives of output with respect to input arguments, this software enables fast and stable optimization using gradient-based algorithms. As starting ground for development, the OpenAeroStruct, a wing and tail aero-structural design tool [22], was used. Figure 18 comprises the N2 diagram of the design framework, comprising three major groups: Problem Variables; Mission Points Calculations; and Mission Performance and Constraints.

Problem Variables
Includes all state, control, geometric and auxiliary variables that describe the aircraft trajectory, aerodynamic and structural geometry and others related to payload and energy.
Initial values are user-defined, some remain fixed but most are design variables, as listed in Tables 2 and 3. An optional b-spline representation is added for every control and state variable.

Mission Points
Involves aircraft design, aero-structural coupling and propulsion calculations. This group starts with some pre-preprocessing and computation of operating conditions (atmospheric properties and Reynolds number). The take-off distance is calculated as described in Section 2.2.4, and the aero-structural geometric properties and computational meshes are updated.
The framework then handles each collocation point separately, exemplified for a 5-point problem in Figure 18. While the wing shape is controlled outside the mission points, since it is fixed throughout the mission, the tail is handled inside each mission point to allow for a variable stabilator angle to satisfy the pitching moment trim condition at every trajectory point.
In each mission point, the disciplinary modules are evaluated and coupled, namely aerodynamics, structures, propulsion and dynamics (trim conditions). The aero-structural MDA uses the models described in Sections 2.1.1 and 2.1.2, where consistent aerodynamic loads and structural displacements are found. The power and thrust are calculated using the energy decay battery model (Section 2.1.3). The energy decay starts with the take-off model, where the initial consumed battery energy is calculated. Then, the consumed energy is assessed between each mission point and the total remaining battery capacity updated for the next mission point. Finally, some important metrics are calculated to assess the point performance, in particular the horizontal and vertical forces, needed for the collocation method defects calculation to assure the aircraft dynamics (Section 2.1.4), and the pitching moment, needed for the trim calculation.

Mission Performance and Constraints
Evaluates the aircraft and overall performance. This group implements the trajectory collocation method (Section 2.1.4) and evaluates the aircraft mission performance metrics, in particular the mission score (Section 2.3.1). It is also responsible for most of the problem constraints related to trajectory (defects in equations of motion), energy (usable capacity and energy) and geometry (available space inside rhombus box and vertical tail volume).

Results
After stating the selected framework parameters, five different analyses are presented: single-objective solutions for each of the three score contributions-payload, climb and distance; an in-depth analysis of the coupled design and trajectory optimization; and the assessment of the b-spline interpolation numerical efficiency.

Framework Parameters
The MDA solver uses an NLBGS with Aitken's relaxation method [45] to solve the aero-structural coupling, with standard convergence criteria (10 −7 absolute error tolerance and 10 −30 relative error tolerance).
As for the optimizer, the Sequential Least Squares Programming (SLSQP) method was selected, which proved to be a robust gradient-based algorithm. The tolerance was set to 10 −2 as a compromise of engineering converged solutions without excessive computational cost.
Prior to optimization, a mesh convergence study was conducted, running the aerostructural MDA of the baseline aircraft, to select meshes that led to accurate solutions with reduced computational cost. Meshes of size 23 × 5 and 15 × 3 (spanwise × chordwise panels) were selected for the wing and tail discretization, respectively, that led to differences in drag coefficient below 3% when compared to the finest meshes tested.
Similarly for the trajectory module, the effect of the number of collocation points was studied by analyzing the final score (objective function) of the baseline generic feasible aircraft, as shown in Figure 19. A trajectory defined with 30 collocation points was found adequate to accurately define the trajectory while maintaining the computational cost low. This led to optimization problems with 206 design variables and 283 constraints.
Since the final score (Equation (28)) depends on the performance of the best team in the competition, reasonable assumptions had to be made. For the maximum payload P max in Equation (24), it was assumed a cargo of 13 bags of 300 g. The maximum altitude score PS altitude,max in Equation (25) is obtained for a height of 100 m. As for the maximum distance D max in Equation (27), it was assumed an aircraft flying at 24 m s −1 during 120 s.
In all cases, multiple runs were made with different initial aircraft design points to mitigate the fact that gradient-based optimizers converge to local minima. The initial values of all other parameters can found in Tables 2 and 3.
The concurrent optimization framework was run on a desktop computer with an AMD Ryzen™ 5 2600X processor, running at 3.8 GHz, and 16 GB of RAM.

Specialized Design Solutions
As a first design exercise, the framework was used to determine the best solution for each score part (Payload, Climb and Distance) separately, using the inverse of Equations (24), (25) and (27), respectively, as the single-objective function. The resulting optimal aircraft planforms specialized for each score component maximization are depicted in Figure 20.

Payload Score Maximization
To carry the most cargo, the aircraft should employ the largest possible lifting surfaces. The optimal aircraft planform is presented in Figure 20a where, as expected, a very large wing (0.72 m 2 ) and tail (0.15 m 2 ) are designed, extending the most inside the allowable area inside the rhombus box (represented in dashed line), while reaching the allowable upper bound values for both the wing chord (0.40 m) and tail chord (0.25 m). Given the large wing chord, the wingbox is sized using front and rear spar heights that fit inside the smallest thickness-to-chord airfoil ratios (8%), which helps reducing the parasite drag. The complete optimal design variables are listed in Table 4.
The cargo capacity is almost 3 kg, which means the aircraft can carry nine bags. Since both the climb capability and distance traveled are irrelevant for this case, and the climb rate is kept relatively low, reaching a height of just 38 m after 60 s and the speed is reduced, leading to a covered distance of just 2829 m.

Climb Score Maximization
The climb segment lasts the first 60 s of flight when the altitude of the aircraft shall reach 100 m for maximum points. In contrast to the previous case, the wing area should be decreased to obtain a lighter aircraft, as the payload capacity is now irrelevant. This is confirmed by looking at the optimal aircraft planform in Figure 20b where, as expected, a much smaller wing and tail are designed, being the smallest possible while still satisfying the stall constraint.
In this case, the cargo capacity is reduced to just one bag and the total weight is shy of 2 kg. In contrast, the climb capability is enhanced, with focus on the first flight segment so that an impressive 114 m height is achieved at the target 60 s mark. This shift of objective is clearly demonstrated when looking at the optimized vertical trajectory shown in Figure 21a.

Distance Score Maximization
Following the same methodology as before, the optimal planform is presented in Figure 20c for the distance score maximization. This configuration is similar to climb optimization with a small light weight aircraft, but uses a different trajectory that increases the flight distance, as seen in Figure 21b. Since the acceleration dependents on the aircraft mass, the wing area is only enough to produce enough lift to trim the aircraft.
Interestingly, the control setting is far from being capped at full power, as shown in Figure 22a, since this would lead to worse L/D and, consequently, reduced range. Instead, a constant low power setting of about 75 W should be sought to extend the flight distance. The complete results from the specialized single-objective optimizations are summarized in Table 4. The highest total score is the payload configuration (2160.5), despite its relatively poor distance and climb performance because a heavier and slower aircraft still scores significant points in both the climb and distance components. The climb configuration is a strong second (2144.9) with a large penalty in the payload performance, as expected, but with a reasonably good distance score. The distance configuration not only greatly penalizes the payload but also the climb performance, since the aircraft is optimized mostly for aerodynamic efficiency, having the highest wing aspect ratio of all. Taking these three specialized single-objective optimal configurations, the multi-objective utopia score would be 692.3 + 943.6 + 930.1 = 2566.0.

Optimal Design and Trajectory for Maximum Total Score
While the single-objective optimal configurations were as expected, the overall best solution combining the three score components can only be found by considering the multi-objective total score (Equation (28)). Figure 23 presents the initial and final aircraft planform and the score optimization history for the maximization of the total score. The aircraft score is maximized by manipulating the wing geometry to generate higher lift, and by adapting the trajectory to favor the altitude gain in the first 60 s to comply with the competition maximum altitude target. Furthermore, the aircraft payload is significantly increased to eleven bags.
The details of this multi-objective optimal design are included in the penultimate column of Table 4. Compared to the single-objective solutions, the multi-objective solution achieves a total score significantly higher (2714.4), despite being inferior at each of the individual score components, which demonstrates the importance of designing and operating the aircraft considering the total score metric.
The aircraft geometric variables indicate three major aspects of the wing geometric preferences for the optimized solution: (i) the wing area is relatively high; (ii) the sweep is preferred close to zero and; (iii) the wing chord is kept at upper bound (Figure 24), representing the wing surface area maximization for a fixed optimized span. Additionally, the thickness of the wing is kept low, which benefits a faster aircraft and increases the climb and distance score. As for the tail, sweep is kept close to null or negative and the area is kept between 15 and 20% of the wing area.
The structural robustness is demonstrated by the non-positive KS failure function values in Figure 25 for the wing and tail spars. The optimal solution of both wingbox spar and skin thicknesses, t spar and t skin , are at the lower bound 0.15 mm. The margin in the structural integrity of the spars is likely due to a combination of two factors: the mechanical properties of the composite material and the high wingbox stiffness. Following Olissipo Air Team's recent studies, the wingbox has the two spars located at 20% and 80% of the wing chord, which produces very high bending and torsional stiffness. Still, the framework evaluates correctly the empty weight of the aircraft, with the optimization solutions in the range of 2-3 kg, given the Olissipo Air Team's prototype ( Figure 2) has an empty mass of 2.5 kg.
Regarding the optimized trajectory depicted in Figure 21, the optimizer focuses on climbing to maximize the climb score and, afterwards, focuses on increasing the speed to increase the distance score. The vertical speed is highest in the first 60 s and rapidly decreases to zero for the remainder of the flight. Consequently, the horizontal speed increases when the aircraft finishes that climb segment.
As for the propulsion, Figure 26a presents the optimal throttle control and the corresponding thrust output. The propulsion system engages with highest effort in the climb stage, as expected, demonstrated by the near maximum throttle setting during that flight segment. The impact of the realistic voltage declining battery model implemented is clearly seen in this stage, where a continuously increasing throttle setting over time is necessary to output a constant mechanical thrust. The low throttle values for the cruise segment indicate that the trimmed flight occurred for 20-40% of the available thrust, depending on the aircraft used. This is due to the flight region constraint: if no altitude limit was present, the throttle setting could have been higher and the aircraft speed would increase, but would change the aircraft balance, which tends to increase the aircraft altitude (different trim condition).
The trim conditions are continuously satisfied by using the stabilator control scheduling shown in Figure 26b, together with the previously mentioned throttle control scheduling in Figure 26a. However, the obtained throttle data also shows that the cruise flight could be further optimized to achieve a trimmed configuration with higher throttle and speed while maintaining similar climb behavior.
Similarly to the specialized designs, this design also uses a minimum sized battery (lower bound of 200 g, holding a capacity of around 2800 mA h), saving weight that is shifted towards payload. This battery provides more than enough energy for the completion of the short mission, as shown in Figure 22b, where the available capacity never gets below 75%. The battery oversizing can be seen as a built-in safety factor to account for not only modeling errors but also to allow for missed go-around landings.

B-Spline Control Functions
The computational efficiency gain was assessed using b-spline interpolations for the control functions considering the same total score maximization problem. Both throttle δ t and stabilator δ stab variables were represented with 20 point b-splines, keeping the same 30 collocation points. The detailed results are summarized in the last column of Table 4.
A positive impact on the optimization efficiency was experienced by using b-spline functions for the control variables, with a reduction of 12% in computational cost.
The throttle b-spline function produced a smoother thrust scheduling, as shown in Figure 27. As a consequence, the obtained trajectory state function was also smoother, thus closer to the real motion, as exemplified with the aircraft vertical displacement in Figure 28. It is worth mentioning that the optimal design obtained using the b-spline control is significantly different from that in Section 3.3, thus revealing that a new (local) maximum was found by the optimizer. All score components increased, totaling 2907.8, representing an improvement of 7.1%. This was obtained with not only a different trajectory ( Figure 28) but also a different aircraft design ( Figure 29) with a slightly higher aspect ratio wing and payload. Although the number of control points for the b-splines is problem dependent, the use of b-splines to describe the control variables should always be adopted as it can significantly reduce the overall optimization cost. However, it was found that representing the state variables (aircraft position and velocity components) also with b-splines brings additional stiffness to the problem since the defect constraints (Equations (10)- (13)) that enforce the aircraft dynamics can become unfeasible due to the reduction in degrees-of-freedom.

Discussion
A framework to study the Air Cargo Challenge competition optimal design was developed, which included problem specific modules to fully comply with the 2022 edition rules, but it can be easily adapted to any similar problem given its modular architecture. This can be achieved by replacing the different aircraft modules (propulsion, geometry) and providing the required constraint functions.
The final complex multidisciplinary and multiobjective problem with hundreds of design variables and constraints was solved in a few hours on a standard computer, with significant savings if using b-splines for control variables, revealing the capability of the framework.
The optimal design characteristics were determined based on different case studies, leading to the following findings: (i) climb maximum height can be achieved and should be the focus of the first 60 s of flight; (ii) the wing area should be maximized within the rhombus box limits to maximize payload capacity; (iii) wing and tail geometries follow similar trends in all cases, which helps defining a optimal conceptual region; (iv) to maximize traveled distance, full throttle should be avoided in the cruise flight segment to reduce aerodynamic losses.
Despite no prototypes have been built and flown to validate the accuracy of the results presented, it was found that the optimal design cases presented in Sections 3.2 and 3.3 were able to capture the correct problem physics.

Conclusions
The concurrent trajectory optimization and aircraft design framework proved to be a valuable tool for students to easily obtain important guidelines not only for the conceptual and preliminary design of the model aircraft but also for operating it during the Air Cargo Challenge competition in order to maximize the overall team score.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: