Multidisciplinary Analysis and Optimization Method for Conceptually Designing of Electric Flying-Wing Unmanned Aerial Vehicles

: Current unmanned aerial vehicles have been designed by applying the traditional approach to aircraft conceptual design which has drawbacks in terms of the individual analysis of each discipline involved in the conception of new aircraft, the reliance on the designer’s experience and intuition, and the inability of evaluating all possible design solutions. Multidisciplinary analysis and optimization focus on solving these problems, by synthesizing all the disciplines involved and accounting for their mutual interaction. This study presents a multidisciplinary analysis and optimization method for conceptually designing electrical ﬂying-wing micro-unmanned aerial vehicles. The conceptual design task was formulated as a non-linear mathematical programming problem. The method considers the trimming of the UAV during each mission proﬁle phase, consisting of the climb, cruise, and descent. We used two algorithms, one for design space exploration and another for optimization. Typical examples of solving conceptual design problems were considered in the work: the modernization of an existing UAV; the effect of the change of the payload and endurance change on the takeoff weight; and the inﬂuence of different static margins on aerodynamic characteristics. The advantages of using this design method are the remotion of additional internal cycles to solve the sizing equation at each optimization step, and the possibility of not only obtaining a unique optimal solution but also a vector of optimal solutions. formal analysis, O.U.E.B.; investigation, O.U.E.B.; resources, E.K.; data curation, O.U.E.B. and E.K.; writing—original draft preparation, O.U.E.B. and O.L.; writing—review and editing, E.K.; visualization, O.U.E.B.; supervi-sion, O.L.; administration, O.U.E.B.;


Introduction
Unmanned aerial vehicles (UAVs) have been an object of interest since the early years of aviation [1][2][3], but they have only recently become widely used. The interest in this aircraft resides in their implementation in either low human interventions or complete human exclusion tasks. For instance, in traffic monitoring [4][5][6], infrastructure monitoring [7,8], environmental monitoring [9][10][11], and surveillance [12][13][14], UAVs are employed because they require long periods of activity and direct human intervention is limited. Agriculture [15][16][17], terrain mapping [18], and data collection (gathering) are examples of tasks in which the use of UAV cuts operational costs, along with military applications [14], emergency rescue, and natural disaster tasks [19,20], in which the deployment of UAVs is in preference due to the dangerous nature of the performed tasks.
Regardless of the task, the search for a design configuration and parameters that fulfill requirements is the priority and is accomplished during the conceptual design phase. Even though during this phase the time consumption and cost are less than the whole design process, many tasks must be completed, and many important decisions must be taken. Several authors have reported the importance of this design phase [21][22][23][24][25][26]; for example, it is estimated that up to 80% of the costs of an aircraft's life cycle are determined during the conceptual and preliminary design phase [27], reminding us of the importance of this phase. Conceptual aircraft design has seen different approaches. First, they were designed by trial and error. Then, after gathering enough information, a statistical methodology was developed, which is based on empirical data from existing aircraft. This has been the traditional design paradigm until today and has been adequate for past demands. Examples of this design paradigm are found in the literature [21][22][23][24][25][28][29][30]. The methods and approaches founded on these studies are primarily based on empirical formulas, decoupling of the involved disciplines for separated analysis, a sequential approach which converges by iteration, and, more importantly, on the designer's experience and intuition. To illustrate this, the studies [31][32][33][34] are related to the design of UAVs by applying the traditional design paradigm. Other works present a slightly different method for conceptual design. For instance, in work [35], a method for the conceptual design of small aircraft with hybrid-electric propulsion systems was proposed by adapting the existing methods. Similarly, [36] proposed the use of neural networks as a surrogate model, due to the high computational requirements required by the genetic algorithm used for the optimization of the wing geometry. Another example was presented in [37], in which not only genetic algorithms were used, but also differential evolution algorithms for performing multiobjective optimization.
However, current needs demand more efficient aircraft [38]; hence, new methods for aircraft design are required. Concurrent design is a design paradigm aimed at significantly reducing product development time and optimizing product performance, which has been gradually replacing the traditional design paradigm [39,40]. One of the main methodologies of concurrent design is multidisciplinary analysis and optimization (MDAO), a new physicsbased design methodology which considers the relation between disciplines involved in aircraft design, has been widely applied to aircraft design [41][42][43][44][45][46][47][48], and lately, to the design of UAVs. For example, in [49], an MDAO-based methodology for designing an electric UAV, which optimizes aerodynamic configurations and propulsion systems, was presented; or in [50], which performed a multi-objective MDAO for maximizing both the payload weight and cruise duration of a tail-sitter UAV.
In this paper, we present an MDAO design method for electric propulsion flying wing UAVs. Aerodynamics, stability and control, propulsion, and weight and balance were the disciplines considered in the conceptual phase. The objective of the current work was to further increase the model fidelity at the conceptual design level, without compromising the design time. This was accomplished by introducing a trimming procedure during the analysis of the mission profile stage, which considers the aircraft pitch angle by performing penalty-based optimization for minimizing the maximum take-off weight (MTOW), varying geometrical and flight conditions parameters subject to stability, balance, aerodynamic constraints, by solving the sizing equation, and reducing the population size within the optimization for reducing the computational time.
The advantages of using this method in design are the remotion of additional internal cycles to solve the sizing equation at each optimization step; the possibility of not only obtaining the single optimal solution, but also the optimal solution vector, the ability to perform a rapid study of the impact of various parameters and flight conditions on flight performance and technical indicators, while ensuring a sufficiently detailed view of the aircraft designed considering the aircraft trimming, given the static margin, as well as the requirements of pitch control.

Materials and Methods
The problem of selecting the optimal parameters of an electric flying wing UAV was formulated in terms of non-linear mathematical programming: . . , n x subject to g j (x) ≤ 0 j = 1, . . . , n g h l (x) = 0 l = 1, . . . , n h The MTOW W 0 was chosen as the objective function. The design variables were the aspect ratio AR, leading-edge sweep angle Λ LE , taper ratio λ, geometric twist τ, wing loading W / S , flight speed v, the first coefficient of the polynomial of the camber line C, the angle of attack α, and the static margin (x CG − x AC ). The proposed method allowed for the consideration of other types and a number of design variables. However, emphasis was placed on these parameters, as they have the greatest influence on the flight and energetic characteristics of small electric UAVs. The constraints were as follows:

•
Zero longitudinal momentum to trim the aircraft, considering the given static margin; • The constraint on the maximum lift coefficient value; • An inequality constraint on the minimum relative elevon moment arm value to ensure enough UAV pitch control; • The equality constraint on the value of the lift coefficient required to ensure horizontal flight.
The objective function, as well as the constraints, are presented as follows: where: x = [AR, Λ, λ, θ, W / S , V] and x' = [C, α, (x CG − x AC )] are the design variables vector; C M is the pitching moment coefficient; δ CM is the pitch coefficient threshold; C L is the lift coefficient; C L max is the maximum lift coefficient; f pc is the fudge factor for pitch control; l el is the relative elevon moment arm; C L eq is the trimming lift coefficient. The method relies on 2 main steps: 1.
Calculation of the objective function and constraints values.

2.
Solution of the optimization problem.
The objective function and the constraints calculation process is presented in Figure 1.

Calculation of the Objective Function
The sizing equation used for calculating the MTOW was the following: where the absolute weights in kg been the payload Wpayload, equipment Wequip, structure Wstruct, and reserve weight Wreserve are represented in the numerator, while the weight ratios been the power unit Wpu /W0 and batteries weight ratio Wbatt /W0 are in the numerator; i is the flight phase. Mathematical models for determining the energy carrier weight of the batteries are different from those for determining the fuel weight of aircraft with a heat engine. For designing aircraft with electric propulsion systems, the battery weight ratio was calculated as the ratio of specific energy consumption SEC to the specific energy of battery E. The power unit weight ratio was calculated as the product of the engine weight-to-thrust ( W /P)pu by the power-to-weight ratio P /W, thus where SEC = P W t in kWh/kg; ηpu is the power unit efficiency; SoC is the battery state of charge; g is the Earth's gravitational acceleration in m/c 2 . The battery and power unit weight ratio depend on the power-to-weight ratio, which was calculated by solving the flight equilibrium equations with respect to and normal to the flight path: T cos α − D − W sin γ = 0, L − W cos γ + T sin α = 0. (4) where T is the thrust in N; D is the drag force; W is the UAV's weight; L is the lift force; α is the angle of attack; γ is the flight path angle. Obtained from the analysis of the free body diagram shown in Figure 2.

Calculation of the Objective Function
The sizing equation used for calculating the MTOW was the following: where the absolute weights in kg been the payload W payload , equipment W equip , structure W struct , and reserve weight W reserve are represented in the numerator, while the weight ratios been the power unit Wpu / W 0 and batteries weight ratio Wbatt / W 0 are in the numerator; i is the flight phase. Mathematical models for determining the energy carrier weight of the batteries are different from those for determining the fuel weight of aircraft with a heat engine. For designing aircraft with electric propulsion systems, the battery weight ratio was calculated as the ratio of specific energy consumption SEC to the specific energy of battery E. The power unit weight ratio was calculated as the product of the engine weight-to-thrust ( W / P ) pu by the power-to-weight ratio P / W , thus where SEC = P W t in kWh/kg; η pu is the power unit efficiency; SoC is the battery state of charge; g is the Earth's gravitational acceleration in m/c 2 .
The battery and power unit weight ratio depend on the power-to-weight ratio, which was calculated by solving the flight equilibrium equations with respect to and normal to the flight path: where T is the thrust in N; D is the drag force; W is the UAV's weight; L is the lift force; α is the angle of attack; γ is the flight path angle.
Obtained from the analysis of the free body diagram shown in Figure 2. The power-to-weight ratio can be written as ( P /W) = T /W V, then the equation for the required power-to-weight ratio considering the propeller efficiency became: The power-to-weight ratio depends on the aerodynamic efficiency, and for its determination AVL [52], a vortex lattice method was used to calculate the lift and inductive wing drag, and empirical formulas were used for calculating the parasite drag [53].
For small/mini-UAVs, the weight of the flying wing structure primarily depends on the manufacturing technology, more than on its strength. Hence, the following equation, depending primarily on the area and material of the wing, was used for its calculation: where: k is a coefficient accounting for the manufacture; ρ is the density of the material; δ is the thickness of the layer; S is the wing area. The payload, equipment, and reserve weight were considered fixed weights. Therefore, their value depends exclusively on the value assigned by the designer.

Calculation of the Constraints
The calculation of the terms of the sizing equation is performed when the UAV is trimmed. The trimming of the flying wing UAV can be ensured by two parameters: the angle of attack and the airfoil shape; both parameters influence the lift and pitching moment coefficients.
The baseline airfoil was the NACA M6, modeled by a similar method of the combination of camber lines and thickness distributions. The third-order polynomial yc = ax 3 − bx 2 + cx, describes the airfoil. The coefficient a of the polynomial of the camber line is responsible for determining the degree of bending in its tail.
Trimming the aircraft, with respect to the UAV center of gravity x̅ CG and the given static margin, was performed by changing the coefficient C of the camber line and the angle of attack (see Figure 3). The trimming was formulated as the following optimization problem within the main optimization (see Figure 1 The power-to-weight ratio can be written as ( P / W ) = T / W V, then the equation for the required power-to-weight ratio considering the propeller efficiency became: The power-to-weight ratio depends on the aerodynamic efficiency, and for its determination AVL [52], a vortex lattice method was used to calculate the lift and inductive wing drag, and empirical formulas were used for calculating the parasite drag [53].
For small/mini-UAVs, the weight of the flying wing structure primarily depends on the manufacturing technology, more than on its strength. Hence, the following equation, depending primarily on the area and material of the wing, was used for its calculation: where: k is a coefficient accounting for the manufacture; ρ is the density of the material; δ is the thickness of the layer; S is the wing area. The payload, equipment, and reserve weight were considered fixed weights. Therefore, their value depends exclusively on the value assigned by the designer.

Calculation of the Constraints
The calculation of the terms of the sizing equation is performed when the UAV is trimmed. The trimming of the flying wing UAV can be ensured by two parameters: the angle of attack and the airfoil shape; both parameters influence the lift and pitching moment coefficients.
The baseline airfoil was the NACA M6, modeled by a similar method of the combination of camber lines and thickness distributions. The third-order polynomial y c = ax 3 − bx 2 + cx, describes the airfoil. The coefficient a of the polynomial of the camber line is responsible for determining the degree of bending in its tail.
Trimming the aircraft, with respect to the UAV center of gravity x CG and the given static margin, was performed by changing the coefficient C of the camber line and the angle of attack (see Figure 3). The trimming was formulated as the following optimization problem within the main optimization (see Figure 1): The optimization problem was solved using the COBYLA method, which was applied using the Python package OpenMDAO [54].
The maximum lift coefficient of the wing was limited to 0.35, recommended for flying wing configurations.
The constraint of the relative elevon arm was formulated by statistical analysis and limited to 20% of the aspect ratio. It was formulated as follows and is geometrically presented in Figure 4.  The optimization problem was solved using the COBYLA method, which was applied using the Python package OpenMDAO [54].
The maximum lift coefficient of the wing was limited to 0.35, recommended for flying wing configurations.
The constraint of the relative elevon arm was formulated by statistical analysis and limited to 20% of the aspect ratio. It was formulated as follows and is geometrically presented in Figure 4.
where xAC is the position of the elevon aerodynamic center in m. The elevon occupies the relative area bounded from 0.8 to 1 along the chord. The elevon aerodynamic center was located at 85% along the semi-wingspan.

Metrics
The aircraft mass growth factor kMGW [56,57] is defined as and is used for assessing the design solutions, i.e., the weight efficiency of the evaluated design. The elevon occupies the relative area bounded from 0.8 to 1 along the chord. The elevon aerodynamic center was located at 85% along the semi-wingspan.

Metrics
The aircraft mass growth factor k MGW [55,56] is defined as and is used for assessing the design solutions, i.e., the weight efficiency of the evaluated design.

Optimization Methods
In this paper, two optimization methods, the stationary Gauss-Seidel method and a differential evolutionary algorithm based on penalty functions, were used for design space exploration and minimization of the objective function, respectively. The Gauss-Seidel method is written as: where r(u (k) ) = b + Au (k) is the residual vector; E −1 is the lower triangular matrix; is the parameters vector. The Gauss-Seidel method was used along with the response surface method for design space exploration. First, the objective function was calculated for a pair of parameter vectors, while the rest of the parameters remain unchanged. Second, the response surface was constructed. Third, the minimum objective function was calculated from the obtained response surface, and the pair of parameters corresponding to this minimum was now fixed while the next iteration was performed for a new pair of parameters. The algorithm for minimizing the objective function using the Gauss-Seidel method, along with the response surface method, is presented in Figure 5. In this paper, two optimization methods, the stationary Gauss-Seidel method and a differential evolutionary algorithm based on penalty functions, were used for design space exploration and minimization of the objective function, respectively.
The Gauss-Seidel method is written as: where r(u (k) ) = b + Au (k) is the residual vector; E −1 is the lower triangular matrix; u (k) = [x (k) , m0 0 ] is the parameters vector. The Gauss-Seidel method was used along with the response surface method for design space exploration. First, the objective function was calculated for a pair of parameter vectors, while the rest of the parameters remain unchanged. Second, the response surface was constructed. Third, the minimum objective function was calculated from the obtained response surface, and the pair of parameters corresponding to this minimum was now fixed while the next iteration was performed for a new pair of parameters. The algorithm for minimizing the objective function using the Gauss-Seidel method, along with the response surface method, is presented in Figure 5. Alternatively, a differential evolutionary algorithm based on penalty functions was proposed to increase the objective function probability to reach the global minimum. The penalty-based function is formulated as follows: minimize L(x) by varying x The adaptative penalty function was taken from [55]: Alternatively, a differential evolutionary algorithm based on penalty functions was proposed to increase the objective function probability to reach the global minimum. The penalty-based function is formulated as follows: The adaptative penalty function was taken from [57]: where: ψ(x) = ∑ m j=1 max 0, g j (x) ; f(x) is objective function; g j (x) is inequality constraint; U * is an upper bound on the constrained global minimum value.
The value of U * is provided by the user. In our method, U * corresponds to 15 times the payload weight, which corresponds to the maximum aircraft mass growth of existing aircraft reported in [55,56].
The SHADE E-PSR evolutionary algorithm, as presented in [58], was used as the optimization algorithm, along with the exponential population size reduction method [59] for decreasing the optimization execution time.
Particular attention was given to the process of solving the sizing equation, which was directly solved in the optimization process ( Figure 6) and did not require a separate cycle for its convergence. The take-off weight, as a design variable, was solved by the sizing equation with the optimization process (without the need for an extra iteration process for the equation convergence) by narrowing the upper and lower take-off parameter values after each iteration. The lower and upper take-off parameter bounds are updated to the minimum and maximum L(x) that fulfill ψ(x) = 0, or else the take-off bounds keep their previous value. ; f(x) is objective function; gj(x) is inequality constraint; U* is an upper bound on the constrained global minimum value.
The value of U* is provided by the user. In our method, U* corresponds to 15 times the payload weight, which corresponds to the maximum aircraft mass growth of existing aircraft reported in [56,57].
The SHADE E-PSR evolutionary algorithm, as presented in [58], was used as the optimization algorithm, along with the exponential population size reduction method [59] for decreasing the optimization execution time.
Particular attention was given to the process of solving the sizing equation, which was directly solved in the optimization process ( Figure 6) and did not require a separate cycle for its convergence. The take-off weight, as a design variable, was solved by the sizing equation with the optimization process (without the need for an extra iteration process for the equation convergence) by narrowing the upper and lower take-off parameter values after each iteration. The lower and upper take-off parameter bounds are updated to the minimum and maximum L(x) that fulfill ψ(x) = 0, or else the take-off bounds keep their previous value.

Verification of the MDAO Method
A mesh independence study of the vortex lattice method was performed. The variable parameters are the number of panels by semi-wingspan (from 10 to 100 in steps of five) and by chord (from 10 to 30 in steps of five). Thus, 95 types of meshes were created, ranging from coarse to dense, to ensure that the simulation results were sufficiently independent of the mesh. As a control value, the variant with the densest mesh was taken as the reference value. The result of the analysis shows that the aerodynamic results are independent of the mesh density, with a mesh consisting of thirty elements along the chord and ten along the semi-wingspan obtaining the following errors: 0.97%, 0%, and 0.88% for the lift, inductive drag, and pitching moment coefficients, respectively.

Verification of the MDAO Method
A mesh independence study of the vortex lattice method was performed. The variable parameters are the number of panels by semi-wingspan (from 10 to 100 in steps of five) and by chord (from 10 to 30 in steps of five). Thus, 95 types of meshes were created, ranging from coarse to dense, to ensure that the simulation results were sufficiently independent of the mesh. As a control value, the variant with the densest mesh was taken as the reference value. The result of the analysis shows that the aerodynamic results are independent of the mesh density, with a mesh consisting of thirty elements along the chord and ten along the semi-wingspan obtaining the following errors: 0.97%, 0%, and 0.88% for the lift, inductive drag, and pitching moment coefficients, respectively.
The results of the aircraft trim were verified by comparing the obtained values with the geometric solution to this problem. Table 1 presents the input data, and Table 2 presents the design variables of the algorithm for trimming the flying wing UAV. The combination that trims the aircraft has parameters α = 1.13 • and C = 0.2080. The found point coincides with the point by using the COBYLA algorithm (see Figure 7). The error of the required angle of attack and the polynomial coefficient of the airfoil camber line was 0% and 0.1%, respectively. The results of the aircraft trim were verified by comparing the obtained values with the geometric solution to this problem. Table 1 presents the input data, and Table 2 presents the design variables of the algorithm for trimming the flying wing UAV. The combination that trims the aircraft has parameters α = 1.13° and C = 0.2080. The found point coincides with the point by using the COBYLA algorithm (see Figure 7). The error of the required angle of attack and the polynomial coefficient of the airfoil camber line was 0% and 0.1%, respectively. Verification of the objective function calculation algorithm. Table 3 presents the input data for the process shown in Figure 5, calculating the objective function and constraint values.  Verification of the objective function calculation algorithm. Table 3 presents the input data for the process shown in Figure 5, calculating the objective function and constraint values.   Table 4 presents the design variable for the process shown in Figure 6. To verify the performance of the objective function calculation algorithm, the response surfaces of MTOW, the lift-to-drag ratio from aspect ratio and leading-edge sweep angle (see Figure 8a), and the specific wing loading and cruise speed (see Figure 8b) at other "frozen" parameters were plotted. The intervals at which the target function has been calculated are for aspect ratio, from 4 to 10; leading-edge sweep angle, from 0 to 50; wing loading, from 3 to 10 daN/m 2 ; and cruise speed, from 10 to 30 m/s. Figure 8b shows the response surfaces at a specific wing loading of 6 daN/kg and a cruise speed of 20 m/s.  Table 4 presents the design variable for the process shown in Figure 6. To verify the performance of the objective function calculation algorithm, the re sponse surfaces of MTOW, the lift-to-drag ratio from aspect ratio and leading-edge swee angle (see Figure 8a), and the specific wing loading and cruise speed (see Figure 8b) a other "frozen" parameters were plotted. The intervals at which the target function ha been calculated are for aspect ratio, from 4 to 10; leading-edge sweep angle, from 0 to 5 wing loading, from 3 to 10 daN/m 2 ; and cruise speed, from 10 to 30 m/s. Figure   The area of the feasible solution is represented by colors. The part of the surface marked in gray does not fall within the area of feasible values according to the presented constraints.
For the given initial data, the most suitable combination of the two parameters is AR = 4 and Λ = 17. The MTOW is 3.202 kg, and the lift-to-drag ratio is 18.75. It should be said that the solution satisfying the minimum MTOW and maximum L / D exists for the same combination of these parameters. This circumstance is characteristic of small UAVs with an electric propulsion system. The MTOW of such vehicles is largely determined by the batteries' weight, which in turn is directly determined by the lift-to-drag ratio, strongly dependent on the aspect ratio and leading-edge sweep angle, while the considered parameters have no significant impact on the design weight of small UAVs, unlike the manufacturing technology and the absolute dimensions of the wing structure.
In the second iteration, a rational combination of wing loading and cruise speed at "freezing" parameters of AR = 4 and Λ = 17 are selected. Figure 9 shows the obtained response surfaces.
Drones 2022, 6, x FOR PEER REVIEW 11 o The area of the feasible solution is represented by colors. The part of the surf marked in gray does not fall within the area of feasible values according to the presen constraints.
For the given initial data, the most suitable combination of the two parameters is = 4 and Λ = 17. The MTOW is 3.202 kg, and the lift-to-drag ratio is 18.75. It should be s that the solution satisfying the minimum MTOW and maximum L /D exists for the sa combination of these parameters. This circumstance is characteristic of small UAVs w an electric propulsion system. The MTOW of such vehicles is largely determined by batteries' weight, which in turn is directly determined by the lift-to-drag ratio, stron dependent on the aspect ratio and leading-edge sweep angle, while the considered rameters have no significant impact on the design weight of small UAVs, unlike the m ufacturing technology and the absolute dimensions of the wing structure.
In the second iteration, a rational combination of wing loading and cruise speed "freezing" parameters of AR = 4 and Λ = 17 are selected. Figure 9 shows the obtain response surfaces. can be explained the fact that the wing loading determines the wing area and significantly affects the str ture weight; despite the need to increase the power unit and batteries weight when creasing the wing loading, reducing the structure weight leads to a decrease in take weight. At the same time, the maximum lift-to-drag ratio is realized at smaller angles attack and requires a lower flight speed and lower wing loading.
The validation of the SHADE E-PSR algorithm was assessed based on the repeata ity of the results obtained. The optimization parameters of the SHADE E-PSL are p sented in Table 5.  The response surfaces of the objective function are shown in Figure 9a,b, and the region of the feasible solution is represented in color. The line along which the surface breaks define the domain of the sizing equation. For the given input data, the best combination for minimization of W 0 = 1.568 kg was (W/S) = 8.7 daN/kg and v = 20 m/s, and for maximization of L / D = 19.97 was W / S = 4.6 daN/kg and v = 20 m/s. The difference of optimal parameter values in the optimization of different objective functions can be explained by the fact that the wing loading determines the wing area and significantly affects the structure weight; despite the need to increase the power unit and batteries weight when increasing the wing loading, reducing the structure weight leads to a decrease in takeoff weight. At the same time, the maximum lift-to-drag ratio is realized at smaller angles of attack and requires a lower flight speed and lower wing loading.
The validation of the SHADE E-PSR algorithm was assessed based on the repeatability of the results obtained. The optimization parameters of the SHADE E-PSL are presented in Table 5.
The algorithm was executed ten times with different initial populations. To ensure the randomness of the design variable vectors, each initial sample was generated with the Latin hypercube sampling method [60] from the SMT (surrogate modelling toolbox) package [61]. Table 6 presents the results of the validation of the algorithm; the taken sample consists of ten vectors, each vector corresponding to the best of the optimal remaining vector after each algorithm run. Analysis of the results in Table 6 confirms the good performance of the algorithm due to the variation coefficient value of 0.31%.
In addition, a statistical analysis of the design variables was performed to investigate the tendency of the differential evolution algorithm to find the global minimum. The taken sample consists of 60 individuals fulfilling at least one of the algorithms' stopping criteria. Table 7 presents the statistical indicators of the optimization results.

Cases of Study
The presented method for conceptually designing electric flying wing aircraft was tested by solving three case studies. First, by modernizing an existing flying wing UAV. Second, by analyzing the effect of different payloads and endurance on the MTOW. Third, by analyzing different static margins on the aircraft performance characteristics.

Modernization of Existing UAV
The UAV to be modernized is a UAV "Boomerang", designed according to the flying wing configuration with a twin-wing vertical tail, and a single-engine propulsion system consisting of a brushless motor. A layout is shown in Figure 10.
The geometric, flight performance, weight, mass, and aerodynamic characteristics of the UAV "Boomerang" are presented in Table 8. Table 3 presents the input data of the optimization algorithm for calculating the objective function. A couple of changes were introduced to match the characteristics of the "Boomerang"; the mass of the equipment is duplicated, and the flight endurance was changed from 1 h to 40 min. Table 4 presents the intervals of the design variables for UAV parameter optimization. The value of the taper ratio was fixed to 1 due to the almost null correlation with either the objective function or the constraints. Table 5 presents the parameters of the SHADE E-PSR evolutionary algorithm for UAV parameter optimization. Drones 2022, 6, x FOR PEER REVIEW 13 of 23 Figure 10. Boomerang UAV appearance.
The geometric, flight performance, weight, mass, and aerodynamic characteristics of the UAV "Boomerang" are presented in Table 8. Table 3 presents the input data of the optimization algorithm for calculating the objective function. A couple of changes were introduced to match the characteristics of the "Boomerang"; the mass of the equipment is duplicated, and the flight endurance was changed from 1 h to 40 min. Table 4 presents the intervals of the design variables for UAV parameter optimization. The value of the taper ratio was fixed to 1 due to the almost null correlation with either the objective function or the constraints. Table 5 presents the parameters of the SHADE E-PSR evolutionary algorithm for UAV parameter optimization.
The initial population size consists of 60 vectors, in each generation the algorithm reduces the population based on the upper limit of MTOW. Figure 11 shows the evolution of the population over the iterations, and the maximum and minimum values of the penalty function at each iteration. In this study, convergence was achieved after 88 iterations. The design variables corresponding to the optimal vector and their corresponding objective and constraint function values are shown in Table 8. The layout of the resulting optimum UAVs and the UAV subjected to modification are shown in Figure 12. The initial population size consists of 60 vectors, in each generation the algorithm reduces the population based on the upper limit of MTOW. Figure 11 shows the evolution of the population over the iterations, and the maximum and minimum values of the penalty function at each iteration. The geometric, flight performance, weight, mass, and aerodynamic characteristics of the UAV "Boomerang" are presented in Table 8. Table 3 presents the input data of the optimization algorithm for calculating the objective function. A couple of changes were introduced to match the characteristics of the "Boomerang"; the mass of the equipment is duplicated, and the flight endurance was changed from 1 h to 40 min. Table 4 presents the intervals of the design variables for UAV parameter optimization. The value of the taper ratio was fixed to 1 due to the almost null correlation with either the objective function or the constraints. Table 5 presents the parameters of the SHADE E-PSR evolutionary algorithm for UAV parameter optimization.
The initial population size consists of 60 vectors, in each generation the algorithm reduces the population based on the upper limit of MTOW. Figure 11 shows the evolution of the population over the iterations, and the maximum and minimum values of the penalty function at each iteration. In this study, convergence was achieved after 88 iterations. The design variables corresponding to the optimal vector and their corresponding objective and constraint function values are shown in Table 8. The layout of the resulting optimum UAVs and the UAV subjected to modification are shown in Figure 12. In this study, convergence was achieved after 88 iterations. The design variables corresponding to the optimal vector and their corresponding objective and constraint function values are shown in Table 8. The layout of the resulting optimum UAVs and the UAV subjected to modification are shown in Figure 12. The geometric, flight performance, weight, mass, and aerodynamic characteristics of the UAV "Boomerang" are presented in Table 8. Table 3 presents the input data of the optimization algorithm for calculating the objective function. A couple of changes were introduced to match the characteristics of the "Boomerang"; the mass of the equipment is duplicated, and the flight endurance was changed from 1 h to 40 min. Table 4 presents the intervals of the design variables for UAV parameter optimization. The value of the taper ratio was fixed to 1 due to the almost null correlation with either the objective function or the constraints. Table 5 presents the parameters of the SHADE E-PSR evolutionary algorithm for UAV parameter optimization.
The initial population size consists of 60 vectors, in each generation the algorithm reduces the population based on the upper limit of MTOW. Figure 11 shows the evolution of the population over the iterations, and the maximum and minimum values of the penalty function at each iteration. In this study, convergence was achieved after 88 iterations. The design variables corresponding to the optimal vector and their corresponding objective and constraint function values are shown in Table 8. The layout of the resulting optimum UAVs and the UAV subjected to modification are shown in Figure 12.  Notably, as a result of the optimization, we obtained a UAV (Figure 12b,c) of smaller area compared with the original (Figure 12a), which can be explained not only by an increase in the wing loading as a consequence of the reduction in the wing structure weight due to the reduction of its area but also by a reduction in the aircraft weight due to the batteries weight reduction, due to the lowering of the specific energy consumption. Despite some increase in flight speed, the decrease in power consumption was achieved due to an increase in the lift-to-drag ratio caused by an increase in the wing aspect ratio and the cruise flight execution at more optimal angles of attack. The reduction in wing structure weight by 32% led to an improvement in the objective function by about 9.8%. This example also clearly demonstrates the effect of the minimum elevon arm constraint in the longitudinal axis on the optimization result: increasing the wing aspect ratio to improve the lift-to-drag ratio also requires increasing the leading-edge sweep angle of the wing to ensure pitch control. The effect of the constraint is clearly visible in the two competing variants B and C, and it is noticeable that these modifications do not differ much in terms of weight summary. The improvement in take-off weight for the "B" and "C" aircraft relative to the "A" aircraft is shown in Figure 13. Figure 12 illustrates the change in the appearance of the aircraft top view for the variants considered. Moreover, the efficacy of the resulting wings by optimization was demonstrated by comparing the aircraft gross weight factor of the initial and optimized wings; the value of k MGW was reduced by 10.7%.
Drones 2022, 6, x FOR PEER REVIEW 15 relative to the "A" aircraft is shown in Figure 13. Figure 12 illustrates the change in appearance of the aircraft top view for the variants considered. Moreover, the effica the resulting wings by optimization was demonstrated by comparing the aircraft g weight factor of the initial and optimized wings; the value of kMGW was reduced by 10 (a) (b) Figure 13. Weight summary of the UAV modernizing design task: (a) absolute weights in k relative weights.

Effect of Change in Payload and Endurance on MTOW
For the second and third cases of study, the same input data were used, as prese in Table 3. The weight of the propeller is considered part of the reserve weight. Figure 14 shows the weight summary for several payload weights. The ratio change of different payloads to the MTOW and variable weights are presented in Tab It is observed that the change in the payload weight generates a linear increment in variable weights. This is confirmed in Figure 14, in which the dependence of diffe weight components on the payload is more intuitive. Moreover, the battery weigh creases at almost twice the increment of the structure weight.
(a) (b) Figure 14. Weight summary of the optimal vector of the MTOW dependence on the change of load weight: (a) absolute weights in kg; (b) relative weights. Table 9. Component weight ratio with respect to the change in payload weight.

Effect of Change in Payload and Endurance on MTOW
For the second and third cases of study, the same input data were used, as presented in Table 3. The weight of the propeller is considered part of the reserve weight. Figure 14 shows the weight summary for several payload weights. The ratios of change of different payloads to the MTOW and variable weights are presented in Table 9. It is observed that the change in the payload weight generates a linear increment in the variable weights. This is confirmed in Figure 14, in which the dependence of different weight components on the payload is more intuitive. Moreover, the battery weight increases at almost twice the increment of the structure weight.
Drones 2022, 6, x FOR PEER REVIEW 15 relative to the "A" aircraft is shown in Figure 13. Figure 12 illustrates the change in appearance of the aircraft top view for the variants considered. Moreover, the efficac the resulting wings by optimization was demonstrated by comparing the aircraft g weight factor of the initial and optimized wings; the value of kMGW was reduced by 10 (a) (b) Figure 13. Weight summary of the UAV modernizing design task: (a) absolute weights in kg relative weights.

Effect of Change in Payload and Endurance on MTOW
For the second and third cases of study, the same input data were used, as prese in Table 3. The weight of the propeller is considered part of the reserve weight. Figure 14 shows the weight summary for several payload weights. The ratio change of different payloads to the MTOW and variable weights are presented in Tab It is observed that the change in the payload weight generates a linear increment in variable weights. This is confirmed in Figure 14, in which the dependence of diffe weight components on the payload is more intuitive. Moreover, the battery weigh creases at almost twice the increment of the structure weight.    The mass growth factor and its effects are presented in Table 10. If the payload weight is changed by 100%, the mass growth factor changes by −23%. For a 200% conversion, it is −27%. Table 10 and Figure 14 clearly outline the cause of the problem for large electric aircraft: as the payload weight increases, the payload weight efficiency rapidly deteriorates, while the rest of the conditions remain equal, and it is possible to find the minimum value of the mass growth factor of the given design task. For the study of the effect of the endurance on the MTOW, the endurance of 1 hr was increased and the change of the mass growth factor was recorded. The results are summarized in Table 10. Figure 15 shows the weight summary for several endurances. It is observed that the change in the endurance generates a linear increment in the variable weights, as in the previous scenario, except for the structure weight, which falls to the lower interval of the wing loading during the optimization process. The mass growth factor and its effects are presented in Table 10. If the payload we is changed by 100%, the mass growth factor changes by −23%. For a 200% conversion, −27%. Table 10 and Figure 14 clearly outline the cause of the problem for large ele aircraft: as the payload weight increases, the payload weight efficiency rapidly dete rates, while the rest of the conditions remain equal, and it is possible to find the minim value of the mass growth factor of the given design task. For the study of the effect of the endurance on the MTOW, the endurance of 1 hr increased and the change of the mass growth factor was recorded. The results are s marized in Table 10. Figure 15 shows the weight summary for several endurances. observed that the change in the endurance generates a linear increment in the vari weights, as in the previous scenario, except for the structure weight, which falls to lower interval of the wing loading during the optimization process.  Table 11 presents the effect of the endurance change on the MTOW. It was not that if the endurance is changed by 50%, the mass growth factor changes by 48%. If range grows by 100%, the factor grows by 94%. This sensitivity is highly dependen the initial endurance at the maximum fixed weight. In the right column, the effect is sh on the original additional weight. The more endurance grows, the greater the m  Table 11 presents the effect of the endurance change on the MTOW. It was noticed that if the endurance is changed by 50%, the mass growth factor changes by 48%. If the range grows by 100%, the factor grows by 94%. This sensitivity is highly dependent on the initial endurance at the maximum fixed weight. In the right column, the effect is shown on the original additional weight. The more endurance grows, the greater the mass growth factor is. The more the UAV endures, it carries more batteries, which increases its fraction; hence, the mass growth factor is correspondingly high. The greater the sum of the variable weight fractions approaches 1, the more the mass growth factor will grow. By comparing the effect of the change of the payload weight and the change of the endurance on the MTOW, it is observed that the endurance is more sensitive than when increasing the payload. Indeed, the MTOW dependence on endurance is non-linear. Figure 16 shows the main and standard deviations of the angle of attack and the lift-to-drag ratio. The results of the third case of study show that the angle of attack is highly dependent on the static margin. The angle of attack shows little deviation compared to the lift-to-drag ratio, mainly at margins 0 and 0.5. By comparing the effect of the change of the payload weight and the change o endurance on the MTOW, it is observed that the endurance is more sensitive than w increasing the payload. Indeed, the MTOW dependence on endurance is non-linear. Figure 16 shows the main and standard deviations of the angle of attack and the to-drag ratio. The results of the third case of study show that the angle of attack is hi dependent on the static margin. The angle of attack shows little deviation compare the lift-to-drag ratio, mainly at margins 0 and 0.5. An increase in static stability leads to an increase in pitch movement. Trimming UAV, in this case, is ensured by deflecting the elevons upward, or by using airfoils a more upward leading-edge angle. The lift coefficient required for the horizontal f can be maintained by increasing the angle of attack (Figure 16a). The flying wing tri different angles of attack due to changes in the static stability causes changes in the lif drag ratio. By shifting the center of gravity, the static margin, it is possible to trim aircraft at the most favorable angles of attack with the maximum lift-to-drag ratio. H ever, for each degree of the static longitudinal static margin, a choice of optimal para ters was made; for example, for low static margins, it was required to increase the w loading and, in contrast, increasing the static margin led to the need to reduce the w loading to change the angles of attack to a more optimal configuration. Therefore, in An increase in static stability leads to an increase in pitch movement. Trimming the UAV, in this case, is ensured by deflecting the elevons upward, or by using airfoils with a more upward leading-edge angle. The lift coefficient required for the horizontal flight can be maintained by increasing the angle of attack (Figure 16a). The flying wing trim at different angles of attack due to changes in the static stability causes changes in the lift-to-drag ratio. By shifting the center of gravity, the static margin, it is possible to trim the aircraft at the most favorable angles of attack with the maximum lift-to-drag ratio. However, for each degree of the static longitudinal static margin, a choice of optimal parameters was made; for example, for low static margins, it was required to increase the wing loading and, in contrast, increasing the static margin led to the need to reduce the wing loading to change the angles of attack to a more optimal configuration. Therefore, in the range considered, the degree of the aircraft static margin did not significantly affect the lift-to-drag ratio.

Discussion
First, this work presented an MDAO method and two optimization processes for minimizing the MTOW of a flying wing UAV during the conceptual design step by choosing the optimal geometrical and flight conditions. The Gauss-Seidel method and response surface method for a couple of parameters were used to verify the correctness of the obtained results, and to visualize the solution during the differential evolutionary algorithm testing, which also allowed us to trace the correctness of applying restrictions. Similarly, using the response surface method, both the graphical solution of UAV trimming and the solution to the COBYLA algorithm were attained, which were used within the developed MDAO method. The optimization process, based on a differential evolution algorithm and penalty functions, allowed for the generation of not just one, but several design solutions (depending on the minimum threshold specified within the solver) which converged to the same MTOW value with a coefficient of variation of 0.31%. This means that it is highly probable, with the application of the second algorithm, to find the global minimum of the design task. In both optimization processes, the UAV was trimmed by solving an optimization routine within the main optimization process. The twist angle had demonstrated to have a slight influence on the MTOW. It was found that the aspect ratio and the leading-edge sweep angle are highly correlated (not a surprise given that both are linked through the maneuverability constraint) making it possible to reduce the number of design variables, by making one dependent on the other, in future works.
Second, results on the modernization of an existing UAV showed that, with the application of this method, it is possible to optimize the conceptual design; in our study, by a 22% improvement. Results concerning the investigation on the effect of different payloads and endurance on the MTOW showed that by incrementing the payload weight, the mass growth factor decreases, converging to 2.5; this is an indication that by incrementing the payload up to the value corresponding to the mass growth factor, the UAV becomes more efficient. However, by increasing the payload, the MTOW increases to a point (that depends on the structural material), in which the structure weight calculation must account for the stresses. The MTOW dependence on endurance is non-linear, and, by looking at the divergence of the mass growth factor, we observed that by incrementing the endurance, the task becomes unsolvable. When the sum of the weight ratios is more than 1, the computed value lies outside the sizing equation domain. These studies are valid for mass increasing close to the center of gravity; if the mass increment happens, for instance, at the nose of the UAV, the results are expected to be different. This tells us that the work could be further improved by adding a weight distribution model that could not only more precisely calculate the aircraft weight factor when the weight increment is far from the center of gravity, but also account for a more accurate interaction between stability and weight disciplines.
And finally, the results of the investigation on the effect of different static margins on the aerodynamic characteristics showed that one of the main benefits of choosing a less stable, even a neutral or unstable configuration, would be that trimming the UAV is possible at a lower angle of attack; especially important with a flying wing configuration. However, a negative impact could be that the uncertainty regarding lift-to-drag ratio is higher.
Future works on this topic will be aimed toward extending this technique to other aerodynamic schemes of UAVs, introducing into the algorithm the requirements of lateral stability and control, and adding a weight distribution model that could not only more precisely calculate the aircraft weight factor when the weight increment is far from the center of gravity, but also account for a more accurate interaction between stability and weight disciplines.

Conclusions
The developed MDAO method for conceptually designing electric flying-wing UAVs makes it possible to obtain a sufficiently detailed view of the designed aircraft, including its optimal geometric, flight parameters, as well as the corresponding flight performance and weight summary of the components at a sufficiently low computational and time costs. Several used constraints allow us to obtain the optimal solution, consider aircraft trimming with a given static margin, and consider the requirements of pitch control. The method was verified by performing a mesh independence study, comparing optimization solutions for the trimming of the UAV against geometrical solutions, assessing the SHADE E-PSR algorithm, and evaluating the deviation of the objective function.
The introduction of the MTOW as a parameter within the differential evolutionary algorithm and the consequent reduction of its upper interval value by taking the maximum MTOW, in the previous iteration, made it possible to dispense of an extra internal cycle for solving the sizing equation at each optimization step; therefore, increasing the algorithm performance. Another advantage of using this method in design is the possibility to obtain not just one optimal solution, but a few concurrent options due to a low coefficient of variance; the values of the design variables can significantly vary, indicating the possible presence of several minima towards the global minima. The application of traditional optimization methods would lead to only one of the minima, ignoring the others.
The application of this method in UAV conceptual design will make it possible to quickly study the impact of changes in various parameters and flight conditions on flight performance, and technical and economic indicators, which are relevant in the design of new vehicles, and in the modernization of existing ones. Typical examples of solving such design problems have been considered in the work.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, O.U.E.B., upon reasonable request.