Atmospheric Aircraft Conceptual Design Based on Multidisciplinary Optimization with Differential Evolution Algorithm and Neural Networks

1. Introduction

The growing importance of aviation technology in various fields of human activity today requires the emergence of new, more effective manned and unmanned aerial vehicles (UAVs). A key role in ensuring the competitiveness of newly created models or modifications of existing aircraft is played by the initial stages of design, during which up to three-quarters of the main technical decisions are made. The success of the project as a whole depends on the correctness of these decisions [1,2,3,4,5,6,7].

Traditional methods of aircraft design [1,2,3,4,5,6,7,8] involve a sequential consideration of the problems of aerodynamics, flight dynamics, strength, weight calculations and a number of other disciplines. This determines the iterative character of the initial stages of design, which requires repeated execution of computational approximations in each discipline. Many solutions are based on the experience of the designer and the established traditions in the development team.

The development of high-performance computing technology has made possible the emergence of a new design paradigm of "concurrent design" [9,10,11], which involves parallel consideration of weight, energy and aerodynamic efficiency indicators at the early stages of design through the use of multidisciplinary optimization (MDO) methods and numerical mathematical modeling. This approach is already being introduced into the practice of designing manned [12,13,14,15,16,17,18] and unmanned aerial vehicles [19,20,21,22,23,24,25,26,27]. The use of MDO at the initial stages of design for UAVs is of particular importance, since the unmanned aerial vehicles being created today cover almost the entire variety of existing and possible unusual aerodynamic configurations [28] and require consideration of a larger number of design variables across various types. The key requirements of the paradigm “concurrent design” for the design methods used are the high accuracy of mathematical models and the quick impact of optimization algorithms that allow considering a large number of design variables in the limited time at the initial design stages.

A monograph [29] should be considered a certain stage in the development of automation of aircraft design processes, in which it is proposed to use the formulation and solution of design problems in terms of nonlinear mathematical programming (NLP), including the problems of choosing the appearance of aircraft of various types and configurations [30].

In a number of works on optimizing the appearance of aircraft, various search methods are used, in particular, gradient [24,25,31], but they require the calculation of partial derivatives of the objective function, including the refinement of the calculation of take-off weight according to the sizing equation at each iteration. An additional difficulty in solving problems of choosing a design is the heterogeneity and discreteness of design variables. For example, a hybrid power plants may include a heat engine and an electric motor in various combinations. Therefore, to solve such problems, metaheuristic search methods are being developed, in particular, evolutionary methods borrowed from the natural environment are considered the main idea. In the beginning, these were relatively simple algorithms and programs that solved combinatorial optimization problems from homogeneous elements, for example, the choice of orientation and number of layers of composite material in different parts of the wing skin [32] or the optimization of the structural and geometric parameters of complex composite parts [33,34]. At present, algorithms and programs in this area are intensively developed and implemented in the field of preliminary design of aviation technology, [22,35,36]. Two main components of any evolutionary algorithm are: the selection of the best solutions and randomization. The selection of the best ones ensures that the solutions will converge to the optimal values, while randomization prevents the solutions from getting trapped in the local optima and, at the same time, increases the diversity of the solutions. The good combination of these two components will usually ensure that the overall optimization can be achieved. One of the most used evolutionary algorithms in engineering problems is the differential evolution (DE) algorithm. DE optimizes a problem by maintaining a population of candidate solutions and creating new candidate solutions by combining existing ones according to their simple operators, and then maintaining the candidate solution that has the best score or aptitude in the optimization problem in question. In this way, the optimization problem is treated as a black box that simply provides a measure of quality given a candidate solution. In addition, it presents a great flexibility to incorporate adaptive methods, use different forms of vector coding and to work with different surrogate methods (especially with methods based on deep learning models) [37,38,39,40]. However, the potential of this optimization method in the field of conceptual design of aircraft has so far been demonstrated in a limited number of publications.

The purpose of this work is to show the possibilities of the DE algorithm with certain adaptability in solving multidisciplinary optimization problems of aircraft-type aerial vehicles with many design variables, including structural ones. In addition, the use of the DE algorithm and neural networks for the optimal selection of the wing airfoil of the aircraft is shown.

2. Methods and Models

2.1. Formulation of the Conceptual Design Problem

The problem of choosing the optimal parameters of the appearance of an aircraft-type UAV will be considered in terms of NLP, following [7,29,41] and most recent works in this area. In canonical form, the problem of NLP is formulated as follows [42]:

$$F\left({\mathit{x}}^{\mathit{o}\mathit{p}\mathit{t}}\right)\le F\left(\mathit{x}\right)\forall \mathit{x}\in \mathsf{\Omega },$$

$$\mathsf{\Omega }=\left\{\mathit{x}:q_j\left(\mathit{x}\right)\le 0,h_k\left(\mathit{x}\right)=0\right\}$$

(1)

where:

F($\mathit{x}$) - objective function;

$\mathit{x}=\left\{x_1,x_2,...,x_n\right\}$ - vector of design variables;

${\mathit{x}}^{\mathit{o}\mathit{p}\mathit{t}}$- optimal solution of the problem;

Ω - area of permissible design variables;

q_j($\mathit{x}$) $\le$0, j = 1, …, p- constraints in the form of inequalities;

h_k($\mathit{x}$)=0, k = 1, …, k- constraints in the form of equations;

n, m, l - the number of variables and constraints.

In the proposed methodology, various technical characteristics of UAVs that are of interest to developers can be declared as a objective function. However, due to the complexity of formulating and solving the problem in a fairly general form, the features of the methodology are further considered through the problem of improving technical characteristics - minimizing weight take-off - of two existing heavy drones with different aerodynamic configurations, the U-40 [43] and MQ-1 [44] (see Figure 1).

2.2. Digital UAV Model

When solving problems using NLP methods, the mathematical model of the optimization object must answer two questions:

- What is the value of the objective function at the considered point in the space of the design variables, defined by the vector x?

- Does this point belong to the area of permissible design variables or not? And if “no,” then to what area should the value of the objective function be degraded – increased in minimization problems using a penalty?

2.2.1. Objective Function

In the considered tasks, the take-off weight f($\mathit{x}$) = $W_{TO}$ is chosen as the objective function. If another indicator (energy or transport efficiency, etc.) is selected as the objective function, a sufficiently accurate calculation of the take-off weight is required in any case, since on the basis of this value, the aircraft is balanced, the weight and energy balance is calculated, and absolute geometric parameters are determined. The take-off weight of an aircraft depends on almost all design parameters. To determine it, an iterative solution of the sizing equation is used – the balance of weight [45,46].Interventionary studies involving animals or humans, and other studies that require ethical approval, must list the authority that provided approval and the corresponding ethical approval code.

2.2.2. Constraints

In solving NLP problems, it is convenient to divide restrictions into three groups:

The first is restrictions on the maximum and minimum values of the design variables of the form a_i ≤ x_i ≤ b_i. These constraints determine the range of values of the corresponding variable and are called geometric regardless of their nature. They can be linear, angular, specific wing load, etc.

The second group in the form of inequalities q_j(x) ≤ 0 and the third group in the form of equations h_k(x) = 0 determine the technical requirements for the operation of the aircraft. As a rule, they cannot be expressed explicitly through design variables. To assess the fulfillment of functional constraints, special algorithms are used, in some cases quite complex, such as the vortex-lattice method (VLM) in aerodynamics, statistical data or criteria “Force coefficient” [46] in weight calculation, etc.

Among the second and third groups of restrictions, the following are used in the work:

- UAV equilibrium condition in the vertical plane with a given static stability margin:

$$h_1\left(\mathit{x}\right)=C_m\left(\mathit{x}\right)=0$$

(2)

$$h_2\left(\mathit{x}\right)=C_{L\_\mathit{bal}}(\mathit{x}\left)\text{–}{C}_L\right(\mathit{x})=0$$

(3)

$$q_1\left(\mathit{x}\right)=C_L\left(\mathit{x}\right)\le C_L^{*};$$

(4)

- Constraint on the value of the horizontal tail volume coefficient to ensure the required characteristics of UAV controllability influencing the size of the horizontal tail arm and its relative area:

$$q_2\left(\mathit{x}\right)=c_{HT}(\mathit{x}\left)\text{∈}\right[c_{HT\_min},c_{HT\_max}]$$

(5)

Here: C_m – the coefficient of pitch moment relative to the center mass; C_L – the lift coefficient; C_{L_bal} – lift the coefficient in the equilibrium condition; $C_L^{*}$ - permissible lift coefficient in a given flight mode; $c_{HT}$ – the horizontal tail volume coefficient.

2.2.3. Design Variables

The vector of design variables x includes the geometric and specific energy parameters of the aircraft appearance, as well as the kinematic parameters of the considered flight modes.

A feature of this work is the addition of the input value (preliminary approximation) of the take-off mass $W_{TO }^{in}$into the vector x. On the basis $W_{TO}^{in}$, the output (refined) value of the take-off weight $W_{TO }^{out}$is calculated and then the value of the objective function $W_{TO }^{+}$ = $W_{TO }^{out}$ + ψ, where ψ is the value of the penalty function that implements the requirements of the constraints. Convergence of the optimization process is ensured when $W_{TO}^{in}$≈ $W_{TO}^{out}$= $W_{TO}^{+}$. For the variable $W_{TO}^{in}$, the first index in the vector x should be used, so in this work the proposed vector x consists of 14 design variables: x = [$W_{TO}^{in}$, AR₁, AR₂, ${\mathsf{\Lambda }}_1$, ${\mathsf{\Lambda }}_2$, ${\lambda }_1$, ${\lambda }_2$, δ₁, ${\overline{L}}_2$, ${\bar{S}}_2$, V, W/S, δ₂, α], where $W_{TO}^{in}$- the input value of the take-off weight, [kg]; AR₁, AR₂ - the aspect ratio of the forward lift surface and the aftward lift surface, respectively; ${\mathsf{\Lambda }}_1$, ${\mathsf{\Lambda }}_2$– sweep angle the leading edge of the forward lift surface and the aftward lift surface, respectively, [^o]; ${\lambda }_1$, ${\lambda }_2$ – taper ratio of the forward lift surface and the aftward lift surface, respectively; δ₁ – incidence angle of the forward lift surface, [^o]; ${\overline{L}}_2=\frac{L_2}{\overline{c}}$ – the relative distance between the lift surfaces; ${\overline{S}}_2=\frac{S_2}{S_1}$ – the relative area of ​​the aftward lift surface compared to the forward lifting surface; V – the flight speed, [m/s] and W/S – lift system loading, [kg/m²]. The selection of the values incidence angle δ₂ of the aftward lift surface and the angle of attack α of the aircraft is carried out by the balancing algorithm within the general optimization cycle.

A set of tools for calculating the value of the take-off weight of an aircraft at certain values of design variables and assessing the level of compliance with constraints is an information model of the optimization object in NLP problems.

2.3. Methodology for Selecting the Optimal Parameters of an Aircraft-Type UAV

The proposed method is based on the use of the differential evolutionary optimization algorithm “Successful-History based Adaptive Differential Evolution” (SHADE) [47,48,49] with the use of penalty functions [50,51,52], method of population size reduction [53,54,55] and numerical mathematical modeling. The optimization method used (SHADE) involves the transformation of information contained in populations of individuals.

In the considered problems, an individual is one of the possible variants of the project – vector x_s, the components of which are certain values of the design variables of the designed aircraft x_i:

$${\mathit{x}}_{\mathit{s}}={[x_1, x_2,\text{…} x_i, \text{…} x_n]}_s,s=1,\dots w$$

${\mathit{x}}_{\mathit{s}}$ – the vector design variables of the individual s in the population P_g;

s – index of the individual;

w – number of individuals in the population;

i and n – index and number of design variables in the individual s.

A population P_g is a set of individuals s at an iteration of optimization (generation) g, includes vectors of individuals ${\mathit{x}}_{\mathit{s}}$:

$${\mathit{P}}_{\mathit{g}}={\left\{{\mathit{x}}_1,{\mathit{x}}_2,\dots {\mathit{x}}_{\mathit{s}},\dots ,{\mathit{x}}_{\mathit{w}}\right\}}_{\mathit{g}}, g=1, \dots , m.$$

Where: g – the index of population.

The methodology of this work is based on an algorithm of [56], but have new approach with significant improvements that allow:

- consider various aerodynamic configurations of UAVs with one or two lift surfaces (flying-wing, normal, duck, tandem);

- design UAVs of different dimensions;

- use various types of power plants on UAVs;

- increase the performance of calculations through the use of an analytical tools and parallel calculations.

The methodology is implemented on the Python platform with the connection of the AVL open source code for aerodynamic calculations [57].

A block-scheme of the methodology is presented in Figure 2.

In general, the methodology includes a number of basic blocks.

• Block 1

the process of entering initial data, design constants, entering the settings of the optimization method, choosing the size of the first NP population, assigning the value of the stop criterion ε, setting the ranges of values for the design variables [x_i(min)… x_i(max)], the values of constraints q_j ≤ 0 as well as specifying the range of possible values of the input take-off weight of the individuals of the initial approximation at the first iteration [$W_{TOg=1}^{in}$_(min)...$W_{TOg=1}^{in}$_(max)].

• Block 2

the process of initialization of the first populations P_g=1, consisting of vectors ${\mathit{x}}_{\mathit{s}}$of individuals. The selection of the values of the design variables x_g=1,i in the vectors of each individual x_g=1,s is carried out randomly from the user-defined range of values design variables [x_i(min)… x_i(max)], using the Latin Hypercube Sampling (LHS) method [58], including the initial take-off weight approximation $W_{TO }^{in}$_g=1,s taking into account a given range [$W_{TOg=1}^{in}$_(min)...$W_{TOg=1}^{in}$_(max)]. The size of the first population should be at least NP_g=1 = 10⋅n, where n is the number of design variables.

• Block 3

Block 3 is designed to calculate the values of the objective function $W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$ of individuals ${\mathit{x}}_{\mathit{s}}$of the population P_g based on the values of their design variables x_i, including the input take-off weight of the initial approximation $W_{TO }^{in}$_g,s.The value of the objective function of each individual is the value of the output (refined) take-off weight $W_{TO }^{out}$(x_g,s), calculated using the sizing equation (see Figure 3, blocks 3.1.s) and supplemented by the value of the penaltyψ (see Figure 3, block 3.2):

$$W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)=\left\{\begin{array}{c}{W}_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if \psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)=0 \\ R\psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)+U^{*} if \psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)>0 \bigwedge W_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)\le U^{*}\\ R\psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)+W_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if \psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)>0 \bigwedge W_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)>U^{*}\end{array}\right.$$

(6)

Where is:

$$\psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)={\psi }_1\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)+{\psi }_2\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$$

(7)

$${\psi }_1\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)=\left\{\begin{array}{c}0 if C_L\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)\le C_L^{*} \\ C_L\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)-C_L^{*} if C_L\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)>C_L^{* } \end{array}\right.$$

$${\psi }_2\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)=\left\{\begin{array}{c}\begin{array}{c}0 if c_{HT}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) \left[c_{HT\_min}, c_{HT\_max}\right] \\ c_{HT}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)-c_{HT\_max} if c_{HT}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)>c_{HT\_max} \end{array}\\ {c_{HT\_min}-c}_{HT}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if c_{HT}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)<c_{HT\_min} \end{array}\right.$$

$\psi \left(\mathit{x}\right)$ – penalty function; $U^{*}$- upper limit of the possible take-off weight; R – the penalty amplification parameter that matches the dimensions and orders of the penalty function and take-off weight values. The value of this parameter is selected expertly taking into account the take-off weight.

The result of block 3 is a vector of values of objective functions $W_{TO}^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$for each individual x_s in the population P_g (see Figure 3).

The use of the differential evolution method allows the calculation process $W_{TO }^{out}$(x_g,s) to be performed using parallel calculations (see Figure 3, blocks 3.1.s) using the Joblib library [59] in order to improve performance, since the vectors of the design variables of each individual do not depend on each other.

Blocks 4-7 are responsible for the process of generating new populations.

A new population P_g+1 of individuals x_g+1,s at each subsequent step of optimization is formed on the basis of the selection of the best individuals x_g,s from the previous populations P_g = {x_1, x _2,… x _{s, …} x _w}_g and crossover populations ${\mathit{P}}_{\mathit{g}}^{\mathit{c}\mathit{r}\mathit{o}\mathit{s}\mathit{s}}$ = {u_1, u_{2, …} u_{s, …} u_w}_gwhich is obtained by crossing a population of P_g and a mutant population ${\mathit{P}}_{\mathit{g}}^{\mathit{m}\mathit{u}\mathit{t}}$ = {v_1, v_2,…v_{s, …} v_w}_g.

• Block 4 – Mutation

A population ${\mathit{P}}_{\mathit{g}}^{\mathit{m}\mathit{u}\mathit{t}}$ includes mutated vectors of individuals ${\mathit{P}}_{\mathit{g}}^{\mathit{m}\mathit{u}\mathit{t}}$ = {v_1, v_2,…v_{s, …} v_w}_g, each of which is calculated by the formula:

$${\mathit{v}}_{\mathit{s}}={\mathit{x}}_{\mathit{s}}+F_s\left({\mathit{x}}_{\mathit{p}\mathit{b}\mathit{e}\mathit{s}\mathit{t}}-{\mathit{x}}_{\mathit{s}}\right)+F_s\left({\mathit{x}}_{\mathit{r}1}-{\mathit{x}}_{\mathit{r}2}\right)$$

(8)

Where:

${\mathit{x}}_{\mathit{p}\mathit{b}\mathit{e}\mathit{s}\mathit{t}}$ – a randomly selected vector from the group of the best vectors of the population P_g. The group of best vectors is formed according to the principle of the minimum value of the objective function $W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$, and the size of the group of best vectors is determined by the settings of the algorithm [47].

${\mathit{x}}_{\mathit{r}1}$ – a random vector from the population P_g, except ${\mathit{x}}_{\mathit{s}}$;

${\mathit{x}}_{\mathit{r}2}$ – a random vector from the combined population P_g and the archive of worst solutions A [47].

The value of the scaling factor F_s is determined by the algorithm settings [47].

• Block 5 – Crossover

Vectors u_s = (u_1, u_2,… u_{i, …} u_n)_s of crossover population ${\mathit{P}}_{\mathit{g}}^{\mathit{c}\mathit{r}\mathit{o}\mathit{s}\mathit{s}}$ = {u_1, u_{2, …} u_{s, …} u_w}_gis formed by crossing the values of mutant vectors v_s =(v_1, v_2,…v_{i, …} v_n)_s of mutant populations ${\mathit{P}}_{\mathit{g}}^{\mathit{m}\mathit{u}\mathit{t}}$ = {v_1, v_2,…v_{s, …} v_w}_gand vectors x_s = (x_1, x _2,… x _{i, …} x_n)_s of the current population P_g = {x_1, x _2,… x _{s, …} x _w}_g.

Crossing is carried out randomly according to the condition:

$$u_{i,s}=\left\{\begin{array}{c}{v}_{i,s} if U\left[\mathrm{0,1}\right]\le {CR}_s or i=i_{rand}\\ x_{i,s } if otherwise \end{array}\right.$$

(9)

${CR}_s$ – the value of the crossover speed;

$U\left[\mathrm{0,1}\right]$ – random number range.

In the event that after the mutation and crossover operations, the new values of the design variables of any individual go beyond the range of variables, they will be reassigned according to the condition:

$$u_i=\left\{\begin{array}{c}\min \left(x_{i }\right) if u_i<\min \left(x_{i }\right)\\ \mathrm{ma}x\left(x_{i }\right) if u_i>\mathrm{ma}x\left(x_{i }\right)\end{array}\right.$$

(10)

• Block 6 – Selection of individuals and the formation of new population

The values of the objective functions $W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)$of the crossover vectors u_s are calculated using block 3 (see Figure 2) and their values are compared with the values of the objective functions ${W_{TO }^{+}(\mathit{x}}_{\mathit{g},\mathit{s}})$of the vectors x_s of the current population P_g. The selection of individuals of the new population P_g+1 for the next optimization step is carried out according to the condition:

$${\mathit{x}}_{\mathit{g}+1,\mathit{s}}=\left\{\begin{array}{c}{\mathit{u}}_{\mathit{g},\mathit{s}} if W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)\le W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)\\ {\mathit{x}}_{\mathit{g},\mathit{s}} if otherwise \end{array}\right.$$

(11)

Accordingly, the values of the objective functions, penalty functions, and input values of the take-off weight of each individual will be obtained for the new population according to the following conditions:

$$W_{TO }^{out}\left({\mathit{x}}_{\mathit{g}+1,\mathit{s}}\right)=\left\{\begin{array}{c}{W}_{TO }^{out}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right) if W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)\le W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)\\ W_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if otherwise \end{array}\right.$$

(12)

$$W_{TO }^{+}\left({\mathit{x}}_{\mathit{g}+1,\mathit{s}}\right)=\left\{\begin{array}{c}{W}_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right) if W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)\le W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)\\ W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if otherwise \end{array}\right.$$

(13)

$$\psi \left({\mathit{x}}_{\mathit{g}+1,\mathit{s}}\right)=\left\{\begin{array}{c}\psi \left({\mathit{u}}_{\mathit{g},\mathit{s}}\right) if W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)\le W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) \\ \psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right) if otherwise\hfill \end{array}\right.$$

(14)

At this stage, an archive of the worst individuals is also formed, which are used for mutation in subsequent optimization iterations:

$${\mathit{A}}_{\mathit{g}+1}\mathit{A} +{\mathit{x}}_{\mathit{g},\mathit{s}}if W_{TO }^{+}\left({\mathit{u}}_{\mathit{g},\mathit{s}}\right)\le W_{TO }^{+}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$$

In this case, a set of vectors of individuals is formed, the value of the penalty function of which is equal to zero ψ = 0 according to the following condition:

$${\mathit{Mxf}}_{\mathit{g}}=\mathit{Ø}+{\mathit{x}}_{\mathit{g}+1,\mathit{s}} if \psi \left({\mathit{x}}_{\mathit{g}+1,\mathit{s}}\right)=\mathbf{0}$$

• Block 7 – Population reduction and convergence of the sizing equation

The convergence of the sizing equation is ensured by narrowing the range of values [$W_{TO}^{in}$_(min)...$W_{TO}^{in}$_(max)] at each subsequent optimization step according to the condition:

$$W_{TOg+1}^{in}\min \left({\mathit{Mxf}}_{\mathit{g}}\right)$$

$$W_{TOg+1}^{in}\max \left({\mathit{Mxf}}_{\mathit{g}}\right)$$

That is, new boundaries [$W_{TOg+1}^{in}$_(min)...$W_{TOg+1}^{in}$_(max)] for the next generation are determined by the most successful individuals (ψ = 0).

If the set of individuals with ψ = 0 at the previous optimization step is empty Mxf_g = ∅ (that is, all $\psi \left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$>0), then the new range of values [$W_{TO}^{in}$_(min)...$W_{TO}^{in}$_(max)] remains the same as in the previous generation:

$${W_{TOg+1}^{in}}_{\left(\min \right)}=W_{TOg}^{in}\min$$

$${W_{TOg+1}^{in}}_{\left(\max \right)}=W_{TOg}^{in}\max$$

At each optimization step, the population size is reduced by excluding the worst individuals with the highest values of the objective function $W_{TO}^{+}$ (ψ >> 0). The population size of the next generation is defined by an exponential population size reduction method:

$${NP}_{g+1}=round\left[{NP}_0{\left(\frac{{NP}_{min}}{{NP}_0}\right)}^{\frac{NFE}{{NFE}_{max}}}\right]$$

(15)

where: NP_min – the minimum population, NP₀ – the initial population, NFE_max – the maximum number of evaluated functions, NFE is the current number of evaluated functions. The NP_g+1 individuals with the best fitness are selected..

If the size of the new population is greater than the number of individuals with ψ = 0, then the population will be supplemented by individuals with the lowest values of the objective function of those with ψ ≠ 0.

The input value of the take-off weight of the next approximation $W_{TOg+1, s}^{in}$for any individual is a product of the mutation and crossing process of population at the previous step. If these actions result in the input take-off weight value $W_{TOg+1, s}^{in}$of an individual of the new generation going beyond the boundary of the new value range [$W_{TOg+1}^{in}$_(min)…$W_{TOg+1}^{in}$_(max)], then this individual will be reassigned a boundary value from this range.

Therefore, individuals with ψ ≠ 0 will be gradually eliminated in accordance with the process of population reduction, and the boundaries [$W_{TO}^{in}$_(min)...$W_{TO}^{in}$_(max)] will narrow until the convergence of the sizing equation $W_{TO}^{in}$ _≈ $W_{TO}^{out}$ = $W_{TO}^{+}$together with the convergence of the overall optimization process max$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)$ - min$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)\le \epsilon$.

The convergence of objective function when solving the sizing equation is shown schematically in Figure 4.

• Block 8

Evaluation the convergence condition max$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)$ - min$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)\le \epsilon$. If this condition is not met, blocks of algorithm 3-7 are executed until it converges.

• Block 9

Selects the best individual from the final population ${\mathit{x}}^{\mathit{o}\mathit{p}\mathit{t}}$and the corresponding value of the objective function $W_{TO}^{+}$(${\mathit{x}}^{\mathit{o}\mathit{p}\mathit{t}}$) (as a result, $W_{TO}^{+}$(${\mathit{x}}^{\mathit{o}\mathit{p}\mathit{t}}$) = ${W_{TO}^{out}}^{opt}$). In addition, the result of the algorithm is other output values that are important for comparing the result with alternative solutions and for carrying out further design stages, including the weights of UAV parts, the value of its take-off weight, flight technical, energy, aerodynamic characteristics and appearance (see Figure 2, blocks 9.1-9.2).

The pseudocode of the algorithm (see Figure 2) is presented in Appendix A.

The algorithm for calculating the output (refined) take-off weight $W_{TO}^{out}$(x_g,s) (see Figure 3, blocks 3.1.s) is one of the cycles of the general optimization algorithm and is presented in more detail in Figure 5:

Blockm1 – receives as input the initial data from Block 1 and the vector of the design variables of the individual х_g,s on the generation g, which also includes the initial approximation (input value) of the take-off weight $W_{TOg,s}^{in}$.

Blockm2 is used to calculate the absolute geometric characteristics of the aircraft based on the input value of the take-off weight $W_{TOg,s}^{in}$and the specific load on the lift surface (W/S)_g,s, the value of which is in the vector of the design variables of the considered individual х_g,s.

Geometric characteristics are used in the algorithm to:

- automated generation of three-dimensional geometric models of UAVs;

- generation of numerical models for calculating the aerodynamic characteristics of UAVs;

- application of engineering formulas for aerodynamics, taking into account compressibility and viscous friction;

- calculation of the weight of the UAV airframe structure.

Blockm3 – calculates the aerodynamic characteristics of the UAV in order to determine the lift-to-drag ratio.

Unitm4 – checks the balancing condition of the UAV in the vertical plane.

Blockm5 – implements the process of balancing the UAV by selecting the angle of attack α and the incidence angle δ₂ of the balancing surface, taking into account the specified static margin of longitudinal stability.

Blockm6 – calculates the required power characteristics of the power plants and the required amount of energy carrier at all stages of the flight.

Blockm7 – calculates the weights of the main parts of the UAV.

Blockm8 – calculates the output take-off weight $W_{TO}^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)$ using the sizing equation (16) based on the values of the parameters x_i and the input take-off weight value${ W}_{TO g,s}^{in}$.

• Calculation of the weights of aircraft parts

Taking into account the features of the optimization algorithm, the take-off weight is calculated using formula (16) for each individual x_g,s of the current generation g.

$$W_{TO }^{out}\left({\mathit{x}}_{\mathit{g},\mathit{s}}\right)=\frac{W_{pay}}{1-\left(\frac{W_{pow}}{W_{TO}}\right)-\left(\frac{W_{fuel}}{W_{TO}}\right)-\left(\frac{W_{struct}}{W_{TO}}\right)-\left(\frac{W_{equip}}{W_{TO}}\right)}$$

(16)

where: $W_{pay}$– payload; $\left(\frac{W_{pow}}{W_{TO}}\right)$– relative weight of the power plants; $\left(\frac{W_{fuel}}{W_{TO}}\right)$– relative weight of fuel (or $\left(\frac{W_{bat}}{W_{TO}}\right)$ – relative weight of batteries in the case of an electric power plants); $\left(\frac{W_{struct}}{W_{TO}}\right)$– relative weight of structure; $\left(\frac{W_{equip}}{W_{TO}}\right)$– relative weight of equipment and control.

1.

Weight of energy carrier

In the case of the use of an internal combustion engine UAV (or other type of fuel engine), the relative weight of the fuel is determined by the formula:

$$\left(\frac{W_{fuel}}{W_{TO}}\right)=\left(\frac{P}{W_{TO}}\right)Ct$$

(17)

where: $\left(\frac{P}{W_{TO}}\right)$ – power-to-weight ratio [kW/daN]; C – specific fuel consumption [kg/(kW.h)]; t – the flight endurance [h].

In the case of using an electric motor UAV, the relative weight of the batteries will be equal to:

$$\frac{W_{bat}}{W_{TO}}=\frac{g\left(\frac{P}{W_{TO}}\right)t}{E{\eta }_{pow}Soc}$$

(18)

where: ${\eta }_{pow}$ – efficiency of the power plants; SoC – the state of charge of the battery; E – the specific energy capacity of the battery [kg/(kW.h)]; g – the acceleration of gravity [m/s²].

2.

The weight of the power plants

The weight of the power plants includes the weight of the propellers and the engine with performance support systems – a fuel system for the internal combustion engine and controllers of electric motors in the case of their use as a power unit. The relative weight of the engine is calculated according to the formula:

$$\left(\frac{W_{pow}}{W_{TO}}\right)=k\left(\frac{P}{W_{TO}}\right)\gamma$$

(19)

where: k – coefficient takes into account the increase in the weight of the power plants due to the systems; $\gamma$ – the specific weight of the engine.

3.

Weight of structure and on-board equipment

The weight of the UAV airframe structure, which includes the weights of the wing, tail, fuselage, landing gear, as well as the weight of the onboard equipment, is determined using weight formulas presented in the specialized literature [2,4,6,7,8].

The values included in the denominator in (16) directly or indirectly depend on the take-off weight and vice versa. Equation (16) is solved by successive approximations: it requires the value of the take-off weight of the initial approximation $W_{TO}^{in}$ as an input and gives a refined value $W_{TO}^{out}$ at the output. The convergence of the solution (16) is achieved with a given accuracy $W_{TO}^{in}$ _≈ $W_{TO}^{out}$. The process of convergence of the sizing equation is organized in the general optimization cycle and is described in detail above.

The pseudocode of the take-off weight calculation algorithm is presented in Appendix B.

• Geometric characteristics (see Figure 5, Block m2)

The geometric characteristics of each individual of the current generation x_g,s (see Figure 5, Block m2) are calculated on the basis of the input value of the take-off weight $W_{TOg,s}^{in}$, the lift system loading (W/S)_g,s and other relative geometric parameters from the vector of the individual x_g,s generated by the optimization cycle.

The total area of the lift surfaces of the UAV S_Σ is determined by the formula:

$$S_{\sum }=\frac{W_{TO}}{(W/S)}$$

(20)

where: W_TO – take-off weight [kg]; $W/S$ – the lift system loading [daN/m²].

The remaining absolute geometric characteristics of each lift surface are found by the relationships:

$$S_{\sum }=S_1+S_2$$

(21)

$$b_j=\sqrt{{AR}_j{S}_j}$$

(22)

$$c_{tj}=\frac{2S_j}{b_j\left(1+{\lambda }_j\right)}$$

(23)

$$c_{rj}=\frac{c_{tj}}{{\lambda }_j}$$

(24)

$$\overline{c}=\frac{2}{3}{c}_r\frac{1+\lambda +{\lambda }^2}{1+\lambda }$$

(25)

where: S₁ – the area of the forward lift surface [m²]; S₂ – the area of the aftward lift surface [m²]; b – the span [m]; $AR$ – aspect ratio; j – the index (j = 1 or 2); c_r – root chord [m]; c_t – tip chord [m]; $\lambda$ – taper ratio; $\overline{c}$– the mean aerodynamic chord [m]. The geometric twist $\tau$ is determined by the law of distribution of the incidence angles of the flow sections of the wing.

The geometry of the fuselage is described by relative and absolute parameters: fuselage aspect ratio ${AR}_{\mathrm{f}}$; aspect ratio of the nose and tail fuselage ${AR}_{\mathrm{n}}$, ${AR}_{\mathrm{t}}$; equivalent mid-section diameter $d_{mf}^{\mathrm{э}}$, loading on mid-section (W/S)_mf.

The scientific novelty of the proposed method is the use of geometric parameters $S_{\sum }$and $\overline{S_2}$ = $S_2/S_1$, which allow considering aerodynamic configurations of one or two lift surfaces without dividing them into a normal, a duck, a tandem or a flying-wing. The value $\overline{S_2}$ is a variable parameter during optimization. In the case $\overline{S_2}$< 1, UAV has a normal configuration, in the case $\overline{S_2}$> 1 – duck configuration, if $\overline{S_2}$ ~ 1 – tandem and $\overline{S_2}$ = 0 – this is the flying-wing or tailless configuration.

In this case, the main lift surface is distinguished: the forward if $S_1$> $S_2$and otherwise, the aftward if $S_1<$ $S_2$. Normalization of aerodynamic coefficients and relative geometric characteristics is carried out relative to the mean aerodynamic chord of the main surface.

• Calculation of the aerodynamic characteristics of the UAV (see Figure 5, Block m3)

The calculation of the aerodynamic characteristics of the UAV (see Figure 5, Block m3) is carried out for the aerodynamic configuration of the UAV obtained for each individual x_g,s. Determination of the main properties and the inductive component of drag is carried out using the method VLM [60,61,62] using open software AVL [57] in conjunction with Python, and taking into account the compressibility and calculating the viscous friction forces is carried out on the basis of engineering methods.

Pseudocode for creating a calculation file for AVL and determining aerodynamic characteristics is presented in Appendix C.

The main problem of the aerodynamics model used is to determine the lift-to-drag ratio at each of the stages of flight in order to assess the required energy characteristics of the power plants and the amount of energy carrier, determined by the required power-to-weight ratio at various stages of flight (26).

The lift-to-drag ratio, which determines the energy costs at each stage of flight, is determined by the ratio of the lift coefficient C_L and drag coefficient C_D corresponding to this flight mode.

The energy balance is provided at each of the stages of flight by the expression:

$$\left(\frac{P}{W_{TO}}\right)=\frac{V}{{\eta }_{vint}}\frac{\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$D$}\right.sin\gamma +cos\gamma }{sin\alpha +\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$D$}\right.cos\alpha }$$

(26)

where: $\gamma$ – flight path angle [^o]; α - angle of attack [^o]; $\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$D$}\right.$ – lift-to-drag ratio; V – flight speed [m/s]; ${\eta }_{vint}$ – efficiency of the propeller.

• UAV balancing (see Figure 5, Block m4)

Determination of the lift-to-drag ratio and required energy characteristics of the UAV is carried out when the UAV equilibrium is ensured in the vertical plane, that is, the condition is met:

$$\left\{\begin{array}{c}M=0 \left(a\right)\\ L=L_{bal } \left(b\right)\\ D=Tcos+Wsin\gamma \left(c\right)\hfill \end{array}\right.$$

(27)

where: M – the pitch moment relative to the centre mass; L – lift force; $L_{bal}$ – lift required for equilibrium along the axis z_w of the wind coordinate frame; D – drag force and T – the thrust of the power plants; W – the weight of the UAV.

Condition (27) is fulfilled at the expense of (26). The values of forces and moments in equations (27a) and (27b) are more conveniently converted into coefficients.

The value of the lift coefficient $C_{L\_bal}$ required to ensure equilibrium on the axis z_w is determined by the formula:

$$C_{L\_bal}=\frac{g\left(W/S\right)cos\gamma }{\frac{1}{2}\rho V^2}$$

(28)

where: g – the acceleration of gravity [m/s²]; $\raisebox{1ex}{$W$}\!\left/ \!\raisebox{-1ex}{$S$}\right.$ – the lift system loading [kg/m²]; $\gamma$ – flight path angle [^o]; $\rho$ – air density [kg/m³]; V – the flight speed [m/s].

The lift coefficient $C_L$ of the aircraft depends on the angle of attack and the angular location of the lift surfaces relative to each other, and the coefficient of pitch moment relative to the center mass also depends on the selected static margin of longitudinal stability.

The equilibrium problem (27) is an optimization problem to ensure equations (27a) and (27b) with two variables:

$$\left\{\begin{array}{c}{C}_m(\alpha ,{\delta }_2)=0 \left(а\right)\\ C_L(\alpha ,{\delta }_2)-C_{L\_bal}=0 \left(b\right)\end{array}\right.$$

(29)

where: α - angle of attack [^o]; ${\delta }_2$ – the incidence angle of the balancing aerodynamic surface [^o].

The aerodynamic models used determine the linear characteristic of the dependencies in (29), so it is proposed to use the following algorithm with the wide use of analytical methods in order to increase the speed of solving this problem:

- Calculation of the characteristics $C_L$and $C_{mA}$by numerical method VLM in the AVL software regarding the leading edge of mean aerolynamic chord of the main lift surface for two different combinations α and ${\delta }_2$;

- Based on the data obtained, a linear approximation of the analytical dependencies ${\mathrm{С}}_L(\alpha ,{\delta }_2)$ = $\frac{\partial C_L}{\partial \alpha }(\alpha -{\alpha }_0+{\delta }_2)$ and $C_{mA}(\alpha ,{\delta }_2)$ = $C_{m0}+{\left(\frac{\partial C_m}{\partial C_L}\right)}_A{\mathrm{С}}_L(\alpha ,{\delta }_2)$ is performed;

- The position of the aerodynamic focus is calculated by the angle of attack relative to the leading edge of mean aerolynamic chord of the main lift surface $x_{AC}^{\alpha }=-\frac{\partial C_m}{\partial C_L}$ and the required position of the center mass is calculated, which will provide a given static margin of longitudinal stability: ${\overline{x}}_{AC}+\left({\overline{x}}_{CG}-{\overline{x}}_{AC}\right)$, where $\left({\overline{x}}_{CG}-{\overline{x}}_{AC}\right)$= Δ - the required static margin of longitudinal stability;

- the dependence is recalculated $C_{m\_CG}(\alpha ,{\delta }_2)$ relative to the required position of the center mass (see Figure 6):

$$C_{m_{CG}}=C_{m0}+C_L=C_{m0}+\frac{C_{mA} – C_{m0}}{{\left(\frac{\partial C_m}{\partial C_L}\right)}_A}·\Delta$$

(30)

- analytical expressions of two lines obtained by the intersection of plane $C_L=C_L(\alpha ,{\delta }_2)$ with plane $C_L=C_{L\_bal}$ and plane $C_{m\_CG}=C_{m\_CG}(\alpha ,{\delta }_2)$ with plane $C_{m\_CG}$ = 0 are found, solutions (29a) and (29b) are found independently of each other (see Figure 7);

- the solution of the problem (27) is analytically calculated, which is the point B = (${\alpha }^{*},{\delta }_2^{*}$) of the intersection of two lines $C_L(\alpha ,{\delta }_2)=C_{L\_bal}$ and $C_{m\_CG}(\alpha ,{\delta }_2)=0$ (see Figure 7).

- Values $C_L\left({\alpha }^{*},{\delta }_2^{*}\right)$ and $C_{m\_CG}({\alpha }^{*},{\delta }_2^{*})$ are also calculated based on the models of the method VLM in the AVL software to verify the balancing conditions.

The pseudocode of the algorithm is presented in Appendix D.

and ${\delta }_2$; (b) Dependence of C_L from $\alpha$and ${\delta }_2$; (c) Solution of the condition UAV balancing.

2.4. Visualization of Optimization Results in the Form of a Planned Projection of the Optimal Option and Its Three-Dimensional Model

After the optimization cycle is completed, a file *.xlsx containing all the geometric parameters of the last generation is saved.

At the end of the main optimization program, an additional module is added to run the FreeCAD program [63]. Through a Python subprocess package.

3D models are built using macro files written in Python. This macro file creates a 3D model based on the geometric parameters contained in the *.xlsx file, saves the model image in .png format, closes the FreeCAD program after the task is completed, and places it in the FreeCAD installation folder. With the “run macro on startup” feature, FreeCAD automatically runs the macro file and the whole process happens in the background on the system.

2.5. Selection of the Wing Airfoil with Multilayer Perceptron

Returning to the $C_L^{*}$ condition used for the optimization of the prototypes. The selection of the wing airfoil of the aircraft is mathematically defined as:

$$f\left({\mathit{x}\mathit{\text{'}}}^{\mathit{o}\mathit{p}\mathit{t}}\right)\ge f\left({\mathit{x}}^{\mathit{\text{'}}}\right)=\frac{C_l^{1.5}}{C_d} \forall \mathit{x}\mathit{\text{'}} \in {\mathsf{\Omega }}_1,$$

$${\mathsf{\Omega }}_1=\left\{{\mathit{x}}^{\mathit{\text{'}}}: C_l^{*}\left({\mathit{x}}^{\mathit{\text{'}}}\right)-C_l\ge 0,y_t\left({\mathit{x}}^{\mathit{\text{'}}}\right)-y_{t,min}\ge 0\right\}$$

(31)

where C_l and C_d are the lift and drag coefficients of the profile; $C_l^{*}$ is the maximum permissible lift coefficient for the given flight condition ($C_l^{*}$= $C_L^{*}$); y_t is the maximum thickness of the profile; y_t,min is the minimum allowable value of y_t; ${\mathit{x}}^{\mathit{\text{'}}}$ is the vector of design parameters. The objective function of maximizing the $\frac{C_l^{1.5}}{C_d}$ parameter was selected in order to improve the endurance of the aircraft [64].

For the resolution of this optimization problem, the same optimization algorithm described in the previous section was used, only adapting it for this optimization task. The design parameters are determined by the Bezier-PARSEC (BP) parameters [65]. The parameters used to create the profiles are the following: r_le – radius of the leading edge; α_te – angle of the camber line at the trailing edge; β_te – angle of the thickness line at the trailing edge; z_te – vertical displacement of the trailing edge; γ_le – angle of the camber line at the leading edge; (x_c, y_c) – position of the maximum value of the camber line; k_c – curvature at the maximum point of the thickness line; (x_t, y_t) – position of the maximum value of the thickness line; k_t – curvature at the maximum point of the thickness line; dz_te – half the thickness of the trailing edge; α is added to these parameters, which is the angle of attack of the airfoil.

Figure 8. BP parameters and Bezier reference points for: (a) thickness line; (b) camber line.

Figure 8. BP parameters and Bezier reference points for: (a) thickness line; (b) camber line.

To determine the aerodynamic coefficients of the profile, a multilayer perceptron (MLP) was used. The MLP is in charge of performing a non-linear mapping of the vector that indicates the aerodynamic coefficients and the design parameters (the BP parameters and the angle of attack) [66]. The MLP used has the following architecture: an input layer with 13 neurons, an output layer with two neurons, and 6 hidden layers (256, 128, 64, 32, 16, 8 neurons, respectively) [40]. The input and output values were normalized with values between 0 and 1, for a better training. The hidden layers have as activation function the RELU function, while the output layer has the sigmoid function. For the training the ADAM optimizer was used and the loss function was the mean square error (MSE). The database used for MLP training was created from tests carried out in XFOIL.

3. Assessment of the Reliability of Models and Methods

3.1. Validation of Mathematical Models of Aerodynamics

The aerodynamic characteristics of the UAV according to the specified parameters are calculated based on the method VLM with the connection of AVL software in combination with the formulas of engineering methods on the Python platform. Validation of the obtained calculation results was carried out on the basis of comparison with experimental data obtained in the T-3 wind tunnel of Samara University by the weight method. A description of the experimental setup is presented in [67,68]. The object of the study is a model with geometric characteristics presented in Figure 9.

For validation, the integral aerodynamic coefficients C_D, C_L, C_m depending on the angle of attack were considered.

Figure 10 present the results of a comparison of experimental data and calculation results.

The validation results showed a good correspondence between the calculation results using the mathematical models and experimental data. The coefficient of variation in the deviation of data on the coefficients C_D, C_L and C_m was 6.8%, 6.1%, and 7.3%, respectively.

3.2. Validation of Models and Algorithm for Calculating the Objective Function

To assess the reliability of the calculation of the objective function and the characteristics obtained from the results of the take-off weight calculation cycle, the calculation of the main characteristics of existing UAVs with known data was carried out using the proposed models: U-40 [43] and MQ-1 [44]. Both UAVs are equipped with Rotax 914 engines. Some characteristics of these devices (geometric, kinematic, specific energy) used in the calculations are given in Table 1.

The convergence of the solution of the sizing equation in the proposed conceptual design methodology is provided in the general optimization cycle. Therefore, to solve the sizing equation in a direct one-time calculation of the objective function with the initial data from Table 1, Table 2 and Table 4, a regular iterative cycle is used. The stages of the flight are considered: take-off, climb, cruise, descent and landing. Taking into account the balancing of the UAV in the vertical plane in all flight modes.

Table 1 shows the results of calculating the objective function - take-off weight, as well as empty UAV weight, fuel weight, engine weight and their maximum power. Similar characteristics of existing considered UAVs [43,44] are also given in Table 1, 4 for comparison with the results of calculations.

The difference in the results of calculating the take-off weight, as well as the weights of the main parts of the two existing UAVs according to the proposed methodology, did not exceed 3-4%. The estimated required power of the power plants and the actual maximum power of the UAV power plants are in the range of 4-5%. The comparison results from Table 1 allow to assess the accuracy of the calculation the main characteristics of the UAV using the proposed methodology, as sufficient for design.

4. Solving Applied Optimization Problems

Using the proposed optimization methodology, the solution of problems for possible improvement of the characteristics of two existing UAVs U-40 and MQ-1 by optimizing a number of their main parameters is considered.

These examples use the general statement of problem (1) presented in Section 2 and use functional constraints (2)-(5).

The objective function is the take-off weight $W_{TO}$.

A number of values used as initial data are taken from Table 1, Table 2 and Table 4. Table 3 shows some of the optimization algorithm settings.

Table 3. Some optimization algorithm settings.

Parameters	Value
Lift coefficient constraint (CL)	< 0.6
Constraint of the horizontal tail volume coefficient (cHT)	[0.2…0.6]
Initial value of penalty (U*) [kg]	60000
Optimization stopping criterion: max$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)$ - min$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)\le \epsilon$	1
Number of design variables (D)	12
Number of initial population (NP0)	10D
Minimum number of generations	12

Table 3. Some optimization algorithm settings.

Parameters	Value
Lift coefficient constraint (CL)	< 0.6
Constraint of the horizontal tail volume coefficient (cHT)	[0.2…0.6]
Initial value of penalty (U*) [kg]	60000
Optimization stopping criterion: max$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)$ - min$\left({\mathit{W}}_{\mathit{T}\mathit{O}}^{+}\left({\mathit{P}}_{\mathit{g}+1}\right)\right)\le \epsilon$	1
Number of design variables (D)	12
Number of initial population (NP0)	10D
Minimum number of generations	12

Table 4. Results of optimization of UAV U-40 and MQ-1 parameters.

	Parameters	Optimization
		Min value	Max value	Optimal value		Initial values
				U-40	MQ-1	U-40	MQ-1
1	${AR}_1$	4	20	14	7.6	20	19
2	${\mathsf{\Lambda }}_1$ [o]	0	25	2.7	0	1	7.22
3	${\lambda }_1$	0.3	1	0.53	0.34	0.33	0.36
4	${\delta }_1$[o]	0	5	2.5	0	2,5	2.5
5	${AR}_2$	4	20	4	18.7	20	6.75
6	${\mathsf{\Lambda }}_2$ [o]	0	25	0	3	1	0
7	${\lambda }_2$	0.3	1	0.63	0.62	0.33	1
8	${\delta }_2$ [o]	-10	+10	-2.2	-2.0	2.2	-2.62
9	${\overline{L}}_2$	3	8	5.1	4.6	5.55	4.37
10	${\overline{S}}_2$	0.2	5	0.2	4.9	1	0.26
11	α [o]	-10	+10	2.8	5.2	-	-
12	V [m/s]	40	90	50	45	55	47
13	W/S [kg/m2]	20	110	92	75	90	73.2
14	WTO [kg]	500	3000	1687	914	2000	1020

Table 4. Results of optimization of UAV U-40 and MQ-1 parameters.

	Parameters	Optimization
		Min value	Max value	Optimal value		Initial values
				U-40	MQ-1	U-40	MQ-1
1	${AR}_1$	4	20	14	7.6	20	19
2	${\mathsf{\Lambda }}_1$ [o]	0	25	2.7	0	1	7.22
3	${\lambda }_1$	0.3	1	0.53	0.34	0.33	0.36
4	${\delta }_1$[o]	0	5	2.5	0	2,5	2.5
5	${AR}_2$	4	20	4	18.7	20	6.75
6	${\mathsf{\Lambda }}_2$ [o]	0	25	0	3	1	0
7	${\lambda }_2$	0.3	1	0.63	0.62	0.33	1
8	${\delta }_2$ [o]	-10	+10	-2.2	-2.0	2.2	-2.62
9	${\overline{L}}_2$	3	8	5.1	4.6	5.55	4.37
10	${\overline{S}}_2$	0.2	5	0.2	4.9	1	0.26
11	α [o]	-10	+10	2.8	5.2	-	-
12	V [m/s]	40	90	50	45	55	47
13	W/S [kg/m2]	20	110	92	75	90	73.2
14	WTO [kg]	500	3000	1687	914	2000	1020

The number of individuals of the initial population was 120.

The convergence of the optimization was achieved after 123 generations and spent 62 minutes on a computer with Windows 10, an Intel(R) Core i7-6700 processor @ 3.40 GHz, 64 GB of RAM. The number of the last generation was 22 individuals. Figure 11 shows the convergence of take-off weight calculation during the optimization process using the algorithm used.

Table 4 shows the results of the parameters optimization. For a visual assessment of the improvements in the characteristics obtained as a result of optimization, the initial characteristics of the considered aircrafts are entered in Table 4.

Table 5. Results of calculating the weight of components and aerodynamic characteristics of the UAV as a result of optimization.

Parameters		U-40					MQ-1
		Initial values	Optimal value		$\mathit{\xi }$, %		Initial values		Optimal value		$\mathit{\xi }$, %
Take-off weight [kg]		2000	1667		16.7		1020		914		10.4
Payload [kg]		600	600		0		204		204		0
Equipment weight [kg]		n/a	137		-		n/a		84.1		-
Weights [kg]	forward lift surface	n/a		181.7		-		n/a		22.4		-
	aftward lift surface	n/a		25.2		-		n/a		123.4		-
	vertical tail	n/a		13.5		-		n/a		19.6		-
	fuselage	n/a		171.6		-		n/a		77.6		-
	landing gear	n/a		84.1		-		n/a		36.5		-
Engine weight [kg]		152	132		13.2		76		76		0
Fuel weight [kg]		450	313.5		30.3		302		266.6		11.7
Propeller weight [kg]		n/a	8.4		-		n/a		3.8		-

Table 5. Results of calculating the weight of components and aerodynamic characteristics of the UAV as a result of optimization.

Parameters		U-40					MQ-1
		Initial values	Optimal value		$\mathit{\xi }$, %		Initial values		Optimal value		$\mathit{\xi }$, %
Take-off weight [kg]		2000	1667		16.7		1020		914		10.4
Payload [kg]		600	600		0		204		204		0
Equipment weight [kg]		n/a	137		-		n/a		84.1		-
Weights [kg]	forward lift surface	n/a		181.7		-		n/a		22.4		-
	aftward lift surface	n/a		25.2		-		n/a		123.4		-
	vertical tail	n/a		13.5		-		n/a		19.6		-
	fuselage	n/a		171.6		-		n/a		77.6		-
	landing gear	n/a		84.1		-		n/a		36.5		-
Engine weight [kg]		152	132		13.2		76		76		0
Fuel weight [kg]		450	313.5		30.3		302		266.6		11.7
Propeller weight [kg]		n/a	8.4		-		n/a		3.8		-

Table 5. cont.

Parameters	U-40			MQ-1
	Initial values	Optimal value	$\mathit{\xi }$, %	Initial values	Optimal value	$\mathit{\xi }$, %
Engine power [kW]	2х84.5	2х62.9	25.6	1х84.5	1х80.6	4.6
Relative position of the center mass$\left({\overline{x}}_{CG}=\frac{x_{CG}}{MAC}\right)$	n/a	0.54	-	n/a	-0.651	-
Cruise lift coefficient	n/a	0.59	-	n/a	0.59	-
Cruise lift-to-drag ratio	n/a	26	-	n/a	24.7	-

Figure 12 and Figure 13 show histograms of the weight distribution of the UAV main parts and a comparison of the required engine power before and after optimization.

Changes in the geometric parameters of the considered UAVs appearance as a result of optimization are shown in Figure 14, the appearance of the prototype is shown in red, and the appearance after optimization is shown in blue.

The figures are obtained automatically using the FreeCAD application, which is used in Batch mode to automatically render the general view of the last population of optimization individuals.

The following design constraints were taken to make the airfoil selection: C_l = 0.59, y_t,min = 11%. For the optimization process, NP₀ = 10*D (D = 11 because z_te and dz_te are kept constant), 100 generations and U* = 0.2 were considered. The ranges of each design parameter are shown in Table 6.

The optimization test was performed 2 times. The results are shown in Table 7 and Figure 15.

Each test required 72 seconds of computing time.

Table 5 shows that the result of optimization allowed to reduce the take-off weight of the considered UAVs and reduce the required power of the power plants.

The UAV U-40 was developed according to the aerodynamic configuration tandem. The results of optimizing its appearance showed that the proposed methodology makes it possible to perform parametric optimization not only of a specific aerodynamic layout, but also to partially solve the problems of circuit synthesis. The tandem configuration with an equal distribution of areas between the wings, all other things being equal, will have the lowest lift-to-drag ratio and the largest weight of the structure. This circumstance is due to the fact that the aftward wing in cruising flight mode (the main mode for long-range UAVs) is "underloaded": in the example considered at ${\overline{S}}_2$= 1, the forward wing creates 64% of the lift, and the aftward wing - 36%. Attempts to load the rear fender result in a normal or a duck configuration. The increase in the weight of the U-40 structure is also associated with the greater aspect ratio of both wings. The proposed algorithm redistributed the air load more rationally and reduced the solution to a normal configuration with the relative areas ${\overline{S}}_2$= 0.2. This solution allowed to slightly reduce the aspect ratio of the main wing and thereby reduce the weight of the wing structure and the UAV as a whole. Due to the increase in lift-to-drag ratio by 3.84% and a decrease in the take-off weight of the UAV by 16.7%, the required power of the power plants was reduced by 25.6%. This would allow the UAV U-40 to replace the existing engines of the power plants (Rotax 914) with less powerful ones, such as the Rotax 912, which would further reduce the weight of the power plants and the UAV as a whole.

The aerodynamic layout of the UAV MQ-1 was chosen by the developers more successfully, but the placement of the heavy power plants in the rear fuselage shifted the center mass back, which required the wing to be shifted back closer to the tail. In order to maintain the required coefficient of the horizontal tail volume coefficient, due to the reduction of the horizontal tail arm, its area was forced to be increased and the dynamic characteristics of the UAV in the longitudinal channel were worsened.

The proposed algorithm made it possible to obtain a different layout - a forward horizontal tail with a smaller area and a longer arm. The forward horizontal tail creates a positive lift and thus somewhat unloads the wing. Thus, the wing area has been reduced and the lift-to-drag ratio of the UAV has been increased. Taken together, this resulted in a 10.4% reduction in take-off weight and 4.6% in the required power of the power plants.

5. Conclusion

A methodology for optimizing the appearance of UAVs is proposed. The problem statement is presented in terms of nonlinear mathematical programming. Certain features of the methodology related to the conceptual design of the aircraft are considered on two real examples of heavy-class UAVs with fundamentally different aerodynamic configurations.

The main attention in the technique is paid to ensuring high accuracy in predicting the take-off weight of an aircraft. For this purpose, the method VLM used in aerodynamic calculations has been tested by our own full-scale experiment in a wind tunnel. The peculiarity of the proposed approach also lies in taking into account and ensuring longitudinal balancing in algorithm for determining the take-off weight. Integrally, the accuracy of the mathematical models, algorithms and programs used in the methodology was tested by a single calculation of the take-off weight, power-to-weight ratio and fuel weight of each of the considered UAVs according to their geometric characteristics. The discrepancy between the published main parameters and those obtained by the calculation is within the range of 3-4%.

Estimation of the accuracy of optimization algorithms in complex multidisciplinary technical problems deserves special consideration and, apparently, does not have any general methods, except for test functions, test models and specially designed problems with a known optimal solution.

The solution of the optimization take-off weight problems of two UAVs using 14 design variables with a given payload and flight endurance, characteristic of the U-40 and MQ-1, helping to significant reduce the take-off weight of these vehicles by making certain changes in appearance and geometric parameters. The obtained results can be considered as a test of the reliability and effectiveness of the methodology as a whole.

When developing the methodology, certain attention was paid to its speed. For this purpose, the optimization algorithm provides the possibility of parallel calculation of the take-off weight of each individual in the population and a one-time calculation of $W_{TO}$ and several other measures. Particularly noteworthy are the results of a computational experiment, in which the current value of the objective function - take-off weight - was added to the number of design variables. According to preliminary studies, this measure speeds up the optimization process by approximately three times, obtaining an equivalent result when compared with the traditional method of convergence of the sizing equation at each optimization step. The analysis of the reasons for such acceleration is of certain scientific interest. Using these accelerations, the time to obtain a solution with 12 design variables is about one hour on a personal computer with an Intel(R) Core i7-6700 @ 3.40 GHz processor, 64 GB of RAM, which, taking into account the complexity of the multidisciplinary combinatorial problem, seems quite acceptable.

The development of a methodology and software for the conceptual design of UAVs with the achieved parameters of accuracy and speed can be considered as a certain step towards concurrent design technology.

The proposed methodology, based on the considered examples, showed the possibility of obtaining new competitive models of unmanned aviation technical or improving the characteristics of existing modifications in terms of the required energy of the power plants and the take-off weight of the device as a whole.

At the moment only the wing profile selection was made using neural networks in the future, it is intended that the proposed methodology makes use of neural networks to predict the aerodynamic characteristics of the aircraft in the balance problem and determine the lift-to-drag ratio instead of the mathematical models used to significantly increase the performance of the algorithm.