Design of Low Altitude Long Endurance Solar-Powered UAV Using Genetic Algorithm

: This paper presents a novel framework for the design of a low altitude long endurance solar-powered UAV for multiple-day ﬂight. The genetic algorithm is used to optimize wing airfoil using CST parameterization, along with wing, horizontal and vertical tail geometry. The mass estimation model presented in this paper is based on structural layout, design and available materials used in the fabrication of similar UAVs. This model also caters for additional weight due to the change in wing airfoil. The conﬁguration is optimized for a user-deﬁned static margin, thereby incorporating static stability in the optimization. Longitudinal and lateral control systems are developed for the optimized conﬁguration using the inner–outer loop strategy with an LQR and PID controller, respectively. A six degree-of-freedom nonlinear simulation is performed for the validation of the proposed control scheme. The results of nonlinear simulations are in good agreement with static analysis, validating the complete design process.


Introduction
At present, aviation focuses on economic, fuel-efficient, and long-endurance systems that incur low operating costs. The growing demand for airborne operations has resulted in the adoption of unmanned aircraft that can remain in the air for significantly longer times than aircraft with pilots on board. The UAV design has been developed extensively to attain these objectives. At present, UAVs can have an endurance of more than a day [1,2].
Engines that consume fuel are not the solution to the ever-growing demands for airborne time. The endurance of an aircraft must satisfy the mission requirement, which can be a few hours (broadcasting a soccer match), a few days (border patrolling, search and rescue missions), or even years (continuous surveillance, communication relay, and meteorological investigations). Exceptional endurance of months or years is possible with a solar-powered aircraft. Photovoltaic cells mounted on wings can be used to capture solar energy during the daytime. A part of this energy can be used directly to run the propulsion system and onboard electronics, and the remaining part can be stored in onboard batteries for use during the night. If sufficient solar energy can be stored in onboard batteries during the day to last the succeeding night, we can design an aircraft that can fly for years. Theoretically, these solar-powered aircraft can fly forever, and their endurance is limited only by the reliability of the subsystems. Solar airplanes can be divided into two categories: high altitude long endurance and low altitude long endurance aircraft. Zephyr [3] and Solara [4] are examples of high altitude long endurance solar aircraft. Designed by QinetiQ, Zephyr achieved three world records in July 2010. The UAV was launched for flight trials on 9 July 2010, and stayed aloft for 14 nights (336 h 22 min) at an altitude of 70,740 ft (21,561 m) above the US Army's Yuma strategy is used with an LQR and PID controller, respectively. Finally, 6-DOF nonlinear simulation is performed in Matlab Simulink to validate the design process and control system.
The remainder of this paper is organized as follows: The design methodology is presented in Section 2, including mass estimation, airfoil parameterization, and stability. The optimization framework is discussed in Section 3. The design of the solar UAV and the optimization results are discussed in Section 4. The dynamic analysis, linear control system design and nonlinear 6-DOF results of the optimized configuration are presented in Section 5. Section 6 presents a few conclusions and discussions, followed by the References section.

Power Required for Level Flight
During steady level flight, the lift and propeller thrust is equal to the weight and drag, respectively: and where W is weight, T is thrust, C L is lift coefficient, C D is drag coefficient, S is surface area, V is the velocity, m is mass and g is gravitational constant. We can determine the velocity from Equation (1): Using Equations (2) and (3), the power required for level flight, P Level , is given by: We replace the surface area in Equation (4) with the aspect ratio, AR, to obtain the following equation:

Daily Energy Requirement
The efficiency of the motor η mot , controller η clrt , gearbox η grb , and propeller η plr must be considered to calculate the power consumption. The power required for the payload P pld and avionics P av is generally known. If the voltage is to be reduced for the payload and avionics, the efficiency of the step-down converter η bec must also be considered. Therefore, the total power consumption P tot is calculated using the equation given in Noth [5].
η clrt η mot η grb η plr P Level + 1 η bec P av + P pld (6) Considering the efficiency of charging η chrg and discharging η dchrg of the batteries, the daily energy consumption E elec tot is given by: where T day and T night are duration of day and night, respectively.

Mass Prediction Model
It is highly important to accurately estimate the mass of a solar UAV because the power required for level flight is directly proportional to it. The mass prediction model in this study is adopted from Noth [5]. The exception is that rather than estimating the structural mass statistically, a new model based on structural layout, design and available materials is proposed.

Fixed Masses
The fixed mass includes the masses of the payload and avionics. These masses are known and do not vary with the aircraft's dimensions: m f ixed = m av + m pld (8) where m f ixed is fixed mass, m pld is payload mass and m av is avionics mass.

Battery Mass
Depending on the mission requirements, the battery adds significant mass to the system. The battery mass is directly proportional to the product of the power required and the duration of the night and is inversely proportional to its energy density [5]. m bat = T night η dchrg k bat P tot (9) where m bat is battery mass and k bat is energy density of battery.

Structural Mass
The most significant challenge in solar UAV design is the estimation of the structural weight. An overestimation of the weight would result in selecting a heavier battery, which would take a long time to charge and increase the weight. An underestimation of the weight would result in selecting a lighter battery, which would be insufficient to last a night. The weight estimation methods discussed in Section 1 are based on statistical data of different manned sailplanes, unmanned radio-controlled models, NASA high-altitude prototypes, and twin-boom configurations. These models are derived from their geometric parameters such as wing area, wing aspect ratio and wingspan, rather than the structural layout, design and materials used. The structural weights of the four solar UAVs estimated using different models are presented in Table 1. It is important to mention here that, for solar UAV design studies discussed in the preceding section, structural weight was quoted for very few prototypes. The most common wing structural layouts are double circular spare [16], single circular spare [7,18] and D-box [5,20,21]. Even for a specific structural layout, there can be many variations in the numbers of ribs, stringers, the diameter and thickness of spares and the thickness of ribs and skin. However, for a given span and aspect ratio, with different structural layouts, any mass prediction models presented in Table 1 will give the same prediction. In the proposed design process, weight estimation is based on the specific The wing structural layout selected for the current study is double circular spares with ribs and stringers. The reasons for this selection are the simplicity and practical experience in fabricating this layout. The selection of the material is also very important for solar UAV design. The materials should be lightweight but adequately strong enough to withstand flight loads. For the current study, the materials considered are 3 k carbon fiber, light density balsa wood and plywood. Several prototypes such as MARAAL [16], Sky-Sailor [5], AtlantikSolar [7], SunSailor [33] and LEEUAV [20] also used carbon fiber and balsa wood as primary materials.
There is no study conducted for structure design. The structure design is adopted from a reference UAV available at the UAV design lab, Beihang University, Beijing. The reference UAV is also a solar-powered UAV with a total weight of 7.5 kg. This UAV has been extensively tested in various experimental flights with a cruise speed up to 15 m/s. The reference UAV has a front and rear spare with a thickness of 1 mm and diameter of 20 mm and 12 mm, respectively. Ribs are made from balsa wood with a thickness of 2 mm. Four stringers run along the span. The stringer at the trailing edge is made from balsa wood. The other three stringers are made from plywood. Wing surface is also balsa wood of 1 mm thickness. The reference UAV also has a third shorter spare adjacent to the front spare. The purpose of this spare is to provide extra stiffness and install batteries if required. The horizontal and vertical tail design is similar to the wing, except they have only one spare with a thickness and diameter of 1 mm and 12 mm, respectively. The tail boom is also a carbon fiber tube with a thickness and diameter of 1 mm and 20 mm, respectively. The fuselage is made from a carbon fiber sheet of 1 mm thickness. To estimate the length of the fuselage, similar solar UAVs are studied (Table 2). From the literature discussed in Section 1, the fuselage length of only three prototypes was quoted. From the geometric data of these three prototypes, it is concluded that the fuselage length could be assumed to be 21% of the wingspan. The fuselage mass is determined based on its length and the average diameter of the reference solar UAV. A heat sink plastic film is used for the wing and tail surfaces in conjunction with balsa wood to save weight. The user can input the percentage coverage of balsa wood (100% implies that the entire wing and tail surface would be of balsa wood, i.e., plastic film would not be used).
The shape of the rib is identical to that of the airfoil and varies as the airfoil undergoes modifications during the optimization process. The user can change the number of holes and their respective locations. The diameters of these holes are input as percentages of the airfoil thickness at the respective locations. This prevents misleading results when the input diameter value is larger than the airfoil thickness at that particular location. The wing airfoil of the reference UAV is NACA 6409, and its tail airfoil is NACA 0009.
The reference UAV is disassembled for measurements. Mass is calculated using geometric dimension, structural layout, design and material densities. Material densities considered for carbon fiber, balsa wood and plywood are 1600 kg/m 3 , 90 kg/m 3 and 260 kg/m 3 , respectively. The structural components of the reference UAV are shown in Figure 1. The reference UAV is disassembled for measurements. Mass is calculated using geometric dimension, structural layout, design and material densities. Material densities considered for carbon fiber, balsa wood and plywood are 1600 kg/m 3 , 90 kg/m 3 and 260 kg/m 3 , respectively. The structural components of the reference UAV are shown in Figure 1. The calculated and actual masses are compared in Table 3. The total estimated mass of the wing also contains the mass of stringers.  The calculated and actual masses are compared in Table 3. The total estimated mass of the wing also contains the mass of stringers. The individual components and total structural mass of the reference UAV and the calculated mass are very close. The advantage of using this type of method is that we can accurately predict the mass given the available materials, structural layout and design. This method is very flexible and can be extended to D-box structural layout. The design of the wing and tail ribs is shown in Figure 2. The modeled airfoils exactly replicate the actual airfoils. The individual components and total structural mass of the reference UAV and the calculated mass are very close. The advantage of using this type of method is that we can accurately predict the mass given the available materials, structural layout and design. This method is very flexible and can be extended to D-box structural layout. The design of the wing and tail ribs is shown in Figure 2. The modeled airfoils exactly replicate the actual airfoils.

Solar Cells Mass
To estimate the mass of the solar cell, the required solar cell area is computed using the following equation: where is the area of the solar cells, is the efficiency of the solar cells, is the efficiency of the camber, is the efficiency of the Maximum Power Point Tracker (MPPT), is the margin factor for clouds and is the maximum irradiation. The mass of the solar cells can be calculated by considering the encapsulation and mass density of the solar cell.
where is the mass of solar cells, and and are solar cells and encapsulation mass density, respectively.

Maximum Power Point Tracker Mass
To extract the maximum power from the solar cell, we must track the optimal working point on its current-to-voltage curve using MPPT. This point shifts continuously because of the irradiance conditions. The mass of the MPPT is calculated using: where is the mass of MPPT and is the mass-to-power ratio of MPPT.

Solar Cells Mass
To estimate the mass of the solar cell, the required solar cell area is computed using the following equation: where A sc is the area of the solar cells, η sc is the efficiency of the solar cells, η cbr is the efficiency of the camber, η mppt is the efficiency of the Maximum Power Point Tracker (MPPT), η wthr is the margin factor for clouds and I max is the maximum irradiation. The mass of the solar cells can be calculated by considering the encapsulation and mass density of the solar cell. m sc = A sc (k sc + k enc ) (11) where m sc is the mass of solar cells, and k sc and k enc are solar cells and encapsulation mass density, respectively.

Maximum Power Point Tracker Mass
To extract the maximum power from the solar cell, we must track the optimal working point on its current-to-voltage curve using MPPT. This point shifts continuously because of the irradiance conditions. The mass of the MPPT is calculated using: m mppt = k mppt I max η sc η cbr η mppt A sc (12) where m mppt is the mass of MPPT and k mppt is the mass-to-power ratio of MPPT.

Propulsion System Mass
The propulsion system consists of an electric motor, controller, electronics, gearbox, and propeller. Based on the statistical analysis from Noth [5], the mass of the entire propulsion system is given by: where m prop and k prop are the mass and mass-to-power ratio of the propulsion system, respectively.

Stability
Stability is generally defined as the capability of an aircraft to return to its equilibrium state after any imbalance, gust, and control input. There are two types of stability: static and dynamic. For an airplane to be statically stable, the sign of the C Mα curve must be negative. A statically stable plane tends to return to its equilibrium position. A statically unstable plane continues to increase the orientation after disturbance. Meanwhile, a statically neutral plane regains its position, implying that the net force or moment acting on the aircraft in the new orientation is zero. Static stability is generally measured in terms of the static margin, which is defined as follows: L re f is the reference length, X CG is the center of gravity and X NP is the neutral point. In the proposed framework, the solar UAV is optimized for a user-specified static margin at the equilibrium position. As the geometry of the UAV is changed during optimization, the neutral point will also change. One methodology is to fix X CG and calculate the static margin using X NP and L re f . In this case, an additional constraint is required to fulfill the static margin requirement. In the current study, the required X CG is calculated using X NP , L re f and the static margin. The aerodynamic moment is then shifted to the required X CG for analysis. For the design of a solar-powered UAV, a static margin of 15% is selected.

Airfoil Parameterization
The airfoil must be parameterized to modify it during optimization. Several methods for airfoil parameterization are available in the literature (Hicks-Henne bump function [34], CST (class function/shape function transformation) [35], and PARSEC [36]). In the present study, CST is used for airfoil parametrization. This method can define a wide range of airfoils and can be extended to other shapes such as squares and circles. It is defined as: The class function C N1 N2 (x) is defined as: where ∆Z defines the trailing edge thickness and x ∈ [0, 1]. For airfoil shapes, N1 = 0.5 and N2 = 1. S(x) is defined as a linear combination of Bernstein polynomials and is given by: a i is the Bernstein coefficient and n is the degree of the polynomial. The quality of the CST airfoil is significantly dependent on the degree of the polynomial. In Figure 3, CST airfoils for polynomials of order zero, two, and four for NACA 63-137 are plotted in conjunction with the original airfoil. It is evident that the polynomial of order four provides a very good representation of the original airfoil. Hence, a fourth-order polynomial is used in the present optimization framework. It provides ten airfoil parameters: five each for the lower and upper portions.
is the Bernstein coefficient and is the degree of the polynomial. The quality of the CST airfoil is significantly dependent on the degree of the polynomial. In Figure 3, CST airfoils for polynomials of order zero, two, and four for NACA 63-137 are plotted in conjunction with the original airfoil. It is evident that the polynomial of order four provides a very good representation of the original airfoil. Hence, a fourth-order polynomial is used in the present optimization framework. It provides ten airfoil parameters: five each for the lower and upper portions.

Optimization Framework
The optimization framework is modeled in Matlab using the genetic algorithm (GA) with selection, mutation and crossover functions. Seventeen design variables are considered. The first 10 design variables are the CST parameters for the airfoil.  This results in very wide design space and the possibility of the production of nonfeasible and unsmooth airfoils. The airfoils produced during optimization are passed through a series of assessments (maximum thickness, camber, curvature of the lower portion, and quality of the trailing edge) to ensure a feasible and smooth airfoil. The details of the design variables are listed in Table 4. There is no geometric and aerodynamic twist in the wing. The wing planform is a rectangle without taper and sweep. The dihedral angle of 7° is used from 45% of the wing semi-span for this study. The tail is a simple Tshape design with the horizontal tail mounted at the top of the vertical tail. The fuselage is not considered in the optimization, as recommended by the developer of Xflr5 (Xflr5

Optimization Framework
The optimization framework is modeled in Matlab using the genetic algorithm (GA) with selection, mutation and crossover functions. Seventeen design variables are considered. The first 10 design variables are the CST parameters for the airfoil. is the Bernstein coefficient and is the degree of the polynomial. The quality of the CST airfoil is significantly dependent on the degree of the polynomial. In Figure 3, CST airfoils for polynomials of order zero, two, and four for NACA 63-137 are plotted in conjunction with the original airfoil. It is evident that the polynomial of order four provides a very good representation of the original airfoil. Hence, a fourth-order polynomial is used in the present optimization framework. It provides ten airfoil parameters: five each for the lower and upper portions.

Optimization Framework
The optimization framework is modeled in Matlab using the genetic algorithm (GA) with selection, mutation and crossover functions. Seventeen design variables are considered. The first 10 design variables are the CST parameters for the airfoil.  This results in very wide design space and the possibility of the production of nonfeasible and unsmooth airfoils. The airfoils produced during optimization are passed through a series of assessments (maximum thickness, camber, curvature of the lower portion, and quality of the trailing edge) to ensure a feasible and smooth airfoil. The details of the design variables are listed in Table 4. There is no geometric and aerodynamic twist in the wing. The wing planform is a rectangle without taper and sweep. The dihedral angle of 7° is used from 45% of the wing semi-span for this study. The tail is a simple Tshape design with the horizontal tail mounted at the top of the vertical tail. The fuselage is not considered in the optimization, as recommended by the developer of Xflr5 (Xflr5  This results in very wide design space and the possibility of the production of nonfeasible and unsmooth airfoils. The airfoils produced during optimization are passed through a series of assessments (maximum thickness, camber, curvature of the lower portion, and quality of the trailing edge) to ensure a feasible and smooth airfoil. The details of the design variables are listed in Table 4. There is no geometric and aerodynamic twist in the wing. The wing planform is a rectangle without taper and sweep. The dihedral angle of 7 • is used from 45% of the wing semi-span for this study. The tail is a simple T-shape design with the horizontal tail mounted at the top of the vertical tail. The fuselage is not considered in the optimization, as recommended by the developer of Xflr5 (Xflr5 GUI). However, the fuselage and tail boom's skin friction drag and form drag are considered in the optimization and calculated using W.H. Mason's FRICTman code [37].
During optimization, Xflr5 is used to obtain aerodynamic coefficients. Xflr5 is an analysis tool for airfoils based on Xfoil's airfoil analysis and also includes wing design and analysis capabilities using the lifting line theory, vortex lattice method and 3D panel model. The validation of Xfoil with experimental data of E387 airfoil at Reynolds number 2.0 × 10 5 is shown in Figure 5. The results are in good agreement in the linear region. As solar UAVs fly at a lower angle of attack, Xflr5 is a good choice for preliminary design.  During optimization, Xflr5 is used to obtain aerodynamic coefficients. Xflr5 is an analysis tool for airfoils based on Xfoil's airfoil analysis and also includes wing design and analysis capabilities using the lifting line theory, vortex lattice method and 3D panel model. The validation of Xfoil with experimental data of E387 airfoil at Reynolds number 2.0 × 10 5 is shown in Figure 5. The results are in good agreement in the linear region. As solar UAVs fly at a lower angle of attack, Xflr5 is a good choice for preliminary design. The analysis is performed in the following sequence: 1. Matlab GA provides 17 design variables. 2. The first 10 variables are used for airfoil generation using CST. The quality of the CST airfoil is assessed. When the airfoil is not smooth, the iteration is terminated, the objective function is assigned a value of 200, and the next set of design variables is considered. When the airfoil is smooth, we proceed to Step 3. 3. The CST airfoil and other geometry parameters are used to perform aerodynamic analysis in Xflr5. A VB script is written that calls the Xflr5.exe, writes all the variables in the respective fields, performs analysis, and writes the output file. 4. Matlab reads the Xflr5 output file. This file contains all the aerodynamic forces and moment coefficients. The pitching moment coefficient is transferred to the required C.G. location calculated from the static and neutral points using Equation (14). 5. The trim (zero pitching moment) angle of attack and the corresponding lift and drag coefficients are calculated. 6. The structural mass of the solar UAV is calculated (Section 2.3.3). 7. The total mass is calculated by solving a cubic equation as presented in [5]. 8. Equations (1), (6) and (9)- (13). If the required solar cell area is greater than the wing area, the iteration is terminated, the objective function is assigned a value of 200, and the next set of design variables is considered. The analysis is performed in the following sequence: 1.
The first 10 variables are used for airfoil generation using CST. The quality of the CST airfoil is assessed. When the airfoil is not smooth, the iteration is terminated, the objective function is assigned a value of 200, and the next set of design variables is considered. When the airfoil is smooth, we proceed to Step 3.

3.
The CST airfoil and other geometry parameters are used to perform aerodynamic analysis in Xflr5. A VB script is written that calls the Xflr5.exe, writes all the variables in the respective fields, performs analysis, and writes the output file.

4.
Matlab reads the Xflr5 output file. This file contains all the aerodynamic forces and moment coefficients. The pitching moment coefficient is transferred to the required C.G. location calculated from the static and neutral points using Equation (14).

5.
The trim (zero pitching moment) angle of attack and the corresponding lift and drag coefficients are calculated. 6.
The structural mass of the solar UAV is calculated (Section 2.3.3). 7.
The total mass is calculated by solving a cubic equation as presented in [5]. 8.
Equations (1), (6) and (9)- (13). If the required solar cell area is greater than the wing area, the iteration is terminated, the objective function is assigned a value of 200, and the next set of design variables is considered. 9.
If the required solar cell area is less than the wing area, the objective function is calculated as follows: The first term ensures that the lift produced is equal to the weight, whereas the second term helps to decrease mass and increase C 3/2 L /C D . The number of generations is set to 100 with a population size of 25.

Design of Solar UAV
The solar UAV is designed to fly over a specific location on a specific day of the year. Given these specifications, t night and t day are retrieved. The designed UAV is oriented to fly for t night hours. The next step is to add robustness to the design [7].
The first consideration is to have multiple-day endurance, t Date exc . This implies a window of a specific number of days or months for flight operation. To achieve this, the additional night duration, owing to date change, must be considered: Meteorological factors such as clouds and water fog in the early morning and evening, t Weather exc , and aggressive flight conditions during the night, t Plevel exc , such as gusts, may result in additional power consumption. The total demand surplus time is: This design is not fully robust. This is because the energy stored in the batteries starts to drain when the available solar power is less than the power required for level flight, which can occur more than 1 h before sunset. To illustrate this, a solar UAV is designed to perform a full-night operation in Beijing or Tianjin (40 • N). The optimum time to perform a full night operation is 21 June, when t day = 14.8 h and t night = 9.2 h. In this design (referred to as Design-1), t req exc = zero. This implies that we assume ideal conditions during flight. The input parameters for the design process are listed in Table 5. The energy simulation for Design-1 is shown in Figure 6a. It is evident that although the battery is designed for t night = 9.2 h, it is insufficient for a continuous 24 h flight and is drained completely before sunrise. This is because the battery starts to supply the energy when the available solar power is less than the power required to maintain level flight. In this case, it is approximately 1.4 h earlier than sunset. In the next design (referred to as Design-2), an additional time of 1.4 h (denoted as t Enight ) is included. The energy simulation for Design-2 is presented in Figure 6b. This design demonstrates the capability for continuous 24 h flight. The battery shows a minimum capacity of 10.8 Wh at sunrise, which corresponds to less than half an hour of flight. Although Design-2 displays potential, it is an ideal design. To add robustness, Equation (22) is incorporated (referred to as Design-3). The associated parameters are as follows: night k cc f is the cloud or fog thickness factor [38] and is considered to be 0.2.
the battery is designed for = 9.2 h, it is insufficient for a continuous 24 h flight and is drained completely before sunrise. This is because the battery starts to supply the energy when the available solar power is less than the power required to maintain level flight. In this case, it is approximately 1.4 h earlier than sunset. In the next design (referred to as Design-2), an additional time of 1.4 h (denoted as ) is included. The energy simulation for Design-2 is presented in Figure 6b. This design demonstrates the capability for continuous 24 h flight. The battery shows a minimum capacity of 10.8 Wh at sunrise, which corresponds to less than half an hour of flight. Although Design-2 displays potential, it is an ideal design. To add robustness, Equation (22) is incorporated (referred to as Design-3). The associated parameters are as follows: is the cloud or fog thickness factor [38] and is considered to be 0.2.
(a) (b) An important parameter of a battery is the state of charge, SOC. It is defined as the current battery capacity to the rated battery capacity. The over-discharging of a battery can attenuate its usable capacity. In addition, the minimum SOC may be limited to 5-10% to increase the battery life cycle [39]. A minimum SOC of 10% is imposed in Design-3. The output parameters for these three designs are listed in Table 6.
An important parameter of a battery is the state of charge, SOC. It is defined as the current battery capacity to the rated battery capacity. The over-discharging of a battery can attenuate its usable capacity. In addition, the minimum SOC may be limited to 5-10% to increase the battery life cycle [39]. A minimum SOC of 10% is imposed in Design-3. The output parameters for these three designs are listed in Table 6. The increase in the mass of the battery and total mass of the UAV owing to the additional time of flight is highly significant. Design-1 has a battery mass of 1.2 kg and a total mass of 4.22 kg. With the addition of t Enight of 1.4 h, Design-2 has a battery mass of 1.52 kg and total mass of 5.04 kg. The further addition of a t req exc of 5.56 h and minimum SOC of 10% increases the battery mass to 3.4 kg and total mass to 7.1 kg. The convergence of the objective function is presented in Figure 7. The mean and best values of the objective function for the last 20 generations are almost equal. This implies that all the configurations in the generation are identical. One design iteration requires 45 s on a personal i5 laptop with 16 GB RAM and 2.43 GHz processor. A total of 100 generations with 25 population sizes may require more than one day. The increase in the mass of the battery and total mass of the UAV owing to the additional time of flight is highly significant. Design-1 has a battery mass of 1.2 kg and a total mass of 4.22 kg. With the addition of of 1.4 h, Design-2 has a battery mass of 1.52 kg and total mass of 5.04 kg. The further addition of a of 5.56 h and minimum SOC of 10% increases the battery mass to 3.4 kg and total mass to 7.1 kg. The convergence of the objective function is presented in Figure 7. The mean and best values of the objective function for the last 20 generations are almost equal. This implies that all the configurations in the generation are identical. One design iteration requires 45 s on a personal i5 laptop with 16 GB RAM and 2.43 GHz processor. A total of 100 generations with 25 population sizes may require more than one day. The lift, drag, and pitching moment coefficients and the lift-to-drag ratio of the optimized configuration are shown in Figure 8a. The trim angle of attack is 3.1°, where the pitching moment is zero. For this trim angle of attack, the lift coefficient is 0.96, and the lift-to-drag ratio is 28.4 (which is nearly the maximum for this design). The Xflr5 model showing panels' density and control surfaces is also shown in Figure 8b. The lift, drag, and pitching moment coefficients and the lift-to-drag ratio of the optimized configuration are shown in Figure 8a. The trim angle of attack is 3.1 • , where the pitching moment is zero. For this trim angle of attack, the lift coefficient is 0.96, and the lift-to-drag ratio is 28.4 (which is nearly the maximum for this design). The Xflr5 model showing panels' density and control surfaces is also shown in Figure 8b. The optimized airfoil, shown in Figure 9, has a thickness of 9.7% at 25% and a maximum camber of 4.6% at 43.5% of the chord. The airfoil is highly smooth, and the trailing edge quality is good. For Design-3, the wing incidence is set at 2°. As the steady-state alpha for Design-3 is The optimized airfoil, shown in Figure 9, has a thickness of 9.7% at 25% and a maximum camber of 4.6% at 43.5% of the chord. The airfoil is highly smooth, and the trailing edge quality is good.
For Design-3, the wing incidence is set at 2 • . As the steady-state alpha for Design-3 is 3.1 • , the wing airfoil encounters a local incidence of 5.1 • with a lift coefficient of 1.07. The Reynolds number for a cruise speed of 8.5 m/s is 1.65 × 10 5 . Fifteen different airfoils suggested for solar UAVs in the literature are analyzed in Xflr5 at a Reynolds number of 1.65 × 10 5 . cl 3/2 /cd for these airfoils for a lift coefficient of approximately 1.07 are plotted in Figure 10. The performance of optimized airfoil is higher than most of the airfoils and very close to Sky-Sailor [5] airfoil. This plot clearly shows the importance of airfoil selection with regard to the design lift coefficient and drag coefficient. These airfoils can be used for this lift coefficient, but a few airfoils evidently operate at higher cl 3/2 /cd. The optimized airfoil, shown in Figure 9, has a thickness of 9.7% at 25% and a maximum camber of 4.6% at 43.5% of the chord. The airfoil is highly smooth, and the trailing edge quality is good. For Design-3, the wing incidence is set at 2°. As the steady-state alpha for Design-3 is 3.1°, the wing airfoil encounters a local incidence of 5.1° with a lift coefficient of 1.07. The Reynolds number for a cruise speed of 8.5 m/s is 1.65 × 10 5 . Fifteen different airfoils suggested for solar UAVs in the literature are analyzed in Xflr5 at a Reynolds number of 1.65 × 10 5 . ⁄ ⁄ for these airfoils for a lift coefficient of approximately 1.07 are plotted in Figure 10. The performance of optimized airfoil is higher than most of the airfoils and very close to Sky-Sailor [5] airfoil. This plot clearly shows the importance of airfoil selection with regard to the design lift coefficient and drag coefficient. These airfoils can be used for this lift coefficient, but a few airfoils evidently operate at higher ⁄ ⁄ .   For Design-3, the wing incidence is set at 2°. As the steady-state alpha for Design-3 is 3.1°, the wing airfoil encounters a local incidence of 5.1° with a lift coefficient of 1.07. The Reynolds number for a cruise speed of 8.5 m/s is 1.65 × 10 5 . Fifteen different airfoils suggested for solar UAVs in the literature are analyzed in Xflr5 at a Reynolds number of 1.65 × 10 5 . ⁄ ⁄ for these airfoils for a lift coefficient of approximately 1.07 are plotted in Figure 10. The performance of optimized airfoil is higher than most of the airfoils and very close to Sky-Sailor [5] airfoil. This plot clearly shows the importance of airfoil selection with regard to the design lift coefficient and drag coefficient. These airfoils can be used for this lift coefficient, but a few airfoils evidently operate at higher ⁄ ⁄ .  The energy simulation for 21 June is presented in Figure 11a. The SOC is also shown in black color. The solar UAV takes off at 07:00 with a completely drained battery (SOC = 0). The available solar power is higher than the required power. The batteries become completely charged (SOC = 1) within 6 h. The available solar energy can be used to attain an altitude or increase speed. At approximately 19:00, the available solar power is less than that required, and the batteries start to supply energy. At 20:00, solar power becomes unavailable, and the UAV operates completely on batteries. The next morning, when the solar power is higher than required, the batteries show a remaining capacity of 335 Wh. This corresponds to an SOC of 0.41. The energy simulation for 21 June is presented in Figure 11a. The SOC is also shown in black color. The solar UAV takes off at 07:00 with a completely drained battery ( = 0). The available solar power is higher than the required power. The batteries become completely charged ( = 1) within 6 h. The available solar energy can be used to attain an altitude or increase speed. At approximately 19:00, the available solar power is less than that required, and the batteries start to supply energy. At 20:00, solar power becomes unavailable, and the UAV operates completely on batteries. The next morning, when the solar power is higher than required, the batteries show a remaining capacity of 335 Wh. This corresponds to an SOC of 0.41.
(a) (b) A similar simulation for 1 May is also presented in Figure 11b. Owing to the lower solar irradiation on 1 May, the battery takes a longer time to charge fully. More battery is consumed to perform full-night flight because of the increased night hours. The next morning, when the available solar power is equal to the required power, the battery shows A similar simulation for 1 May is also presented in Figure 11b. Owing to the lower solar irradiation on 1 May, the battery takes a longer time to charge fully. More battery is consumed to perform full-night flight because of the increased night hours. The next morning, when the available solar power is equal to the required power, the battery shows a minimum SOC of 0.37. It is concluded that Design-3 has the potential to perform continuous flight operations. However, it is important to note that the charge margin time [7] is reduced owing to the date change.

Dynamic Stability and Control System Design
In this section, the response of the optimized solar-powered UAV (Design-3) is discussed by applying the linearized equation of motion. These equations are based on small perturbations from the trim conditions. Additionally, these equations are decoupled. That is, the perturbations in longitudinal forces and moments depend neither on the lateral/directional perturbations nor the lateral/directional control inputs, and the perturbations in lateral/directional forces and moments depend neither on the longitudinal perturbations nor the longitudinal control inputs. Altitude hold and bank to turn control (BTT) systems are developed. The proposed control laws are then implemented in the nonlinear 6-DOF Matlab block for validation.
The first step is to estimate the moments of inertia. As the proposed optimization framework provides only the C.G. location for the desired static margin, a reasonable mass distribution that provides the required C.G. location is assumed. This distribution is passed on to Xflr5 to calculate the moments of inertia [40]. For Design-3, the neutral point and C.G. location are 0.72 m and 0.675 m, respectively (origin at fuselage nose). This yields a static margin of 15%. The estimated mass distribution is shown in Figure 12.

Attitude and Altitude Control
The linearized small disturbance longitudinal rigid body equations of motion in the state-space form are [41]: and = ⎣ ⎢ ⎢ ⎡ The longitudinal state variable vector, , is given by: where , , , and are the forward velocity, vertical velocity, pitch rate, and pitch an-

Attitude and Altitude Control
The linearized small disturbance longitudinal rigid body equations of motion in the state-space form are [41]: where and The longitudinal state variable vector, x long , is given by: where u, w, q, and θ are the forward velocity, vertical velocity, pitch rate, and pitch angle, respectively. A longitudinal control vector, η long , is given by: δ e and δ T are perturbations from the trim in the elevator and throttle settings, respectively. Xflr5 is used to calculate the stability and control derivatives. Xflr5 predictions can be used for preliminary design analysis. These predictions can be improved using high fidelity tools or wind tunnel testing. The definitions and relationships between the dimensional stability derivatives for longitudinal dynamics and dimensionless derivatives of aerodynamic coefficients are presented in [40,41]. The longitudinal derivatives for Design-3 are listed in Table 7. XFLR5 does not compute derivatives with respect to . w. Hence, these are assumed to be zero. Table 7. Longitudinal derivatives of aerodynamic coefficients.

Derivatives
Value The resulting plant matrix is given by: The eigenvalues of the longitudinal plant matrix A are −11.13 + 4.68i, −11.13 − 4.68i, −0.068 + 0.91i, and −0.068 − 0.91i, which corresponds to a short period and the phugoid mode. The short-period mode has a relatively short time period and is generally damped substantially. The phugoid mode has a significantly longer time period and is damped marginally. The properties of the short period and phugoid mode of the present system are given in Table 8. In the attitude and altitude control system design, it is assumed that the forward velocity is controlled by a separate controller and maintained constant at 8.5 m/s. The engine control and actuator dynamics are omitted for simplicity. Control surfaces are modeled as simple trailing edge flaps with x and y hinge position at 70% of chord and 50% of thickness, respectively. The control surface sizing can be refined using high fidelity tools. The resulting state-space model is expressed as: Figure 13a shows the longitudinal open-loop step response. The vertical velocity and pitch rate have non-zero values. The pitch angle continues to decrease dramatically. Therefore, a controller must be designed for the desired pitch angle response. In this study, a linear quadratic regulator (LQR) controller is used to determine the gains for the feedback loop. The goal is to identify a control vector η(t) that drives the state from a specified initial state x(t) to a desired final state x d t f such that a specified performance index (J) is minimized: where g is specified as where Q is a positive semi-definite weight matrix and R is a positive-definite weight matrix. Q and R are selected to provide the optimal response as desired by the designer. For this design: The corresponding gain vector is given by: The feedback response for the attitude-hold control system in response to a 1° step input is shown in Figure 13b. The rise and settling time of the pitch angle are 0.784 s and 1.48 s, respectively, with zero steady-state error. As the designed solar UAV must demonstrate exceptional endurance at a constant altitude, an altitude-hold control system is required. To design this system, we introduce a vertical height equation [41] in the current state-space model.
The resulting state space model is given as:  The feedback response for the attitude-hold control system in response to a 1 • step input is shown in Figure 13b. The rise and settling time of the pitch angle are 0.784 s and 1.48 s, respectively, with zero steady-state error. As the designed solar UAV must demonstrate exceptional endurance at a constant altitude, an altitude-hold control system is required. To design this system, we introduce a vertical height equation [41] in the current state-space model. The resulting state space model is given as: The corresponding LQR gains are [0.2308 − 0.0482 − 1.000], which correspond to the state alpha, pitch rate, and pitch angle, respectively. This system is the same as Equation (30), except vertical velocity is replaced by alpha. A PID controller is suggested in the outer loop that would generate a theta command resulting from the difference between the target height and current height. The parameters of the PID controller are P = 0.07025, I = 0.0003309, D = −0.0300, and N = 0.5372. The response to the 1 m step input is shown in Figure 14. The linear pitch angle accurately matches the pitch angle command, and the target height is achieved with zero steady-state error. The maximum elevator deflection to produce desired pitch angle is 1 • .

Heading Control
The linearized small disturbance lateral rigid body equations of motion in state-space form are given by [41]: Lateral state vector and input control vector are given by: where , , and are the lateral velocity, roll rate, roll angle and yaw rate, respectively. and are the perturbations from trim in the rudder and aileron. It is also convenient to replace lateral velocity with sideslip angle, , using the following relationship.
The resulting plant matrix is given by [42]:

Heading Control
The linearized small disturbance lateral rigid body equations of motion in state-space form are given by [41]: where Lateral state vector x lat and input control vector η lat are given by: where v, p, φ and r are the lateral velocity, roll rate, roll angle and yaw rate, respectively. δ r and δ a are the perturbations from trim in the rudder and aileron. It is also convenient to replace lateral velocity with sideslip angle, β, using the following relationship.
The resulting plant matrix is given by [42]: The dynamics of the resulting system are the same as those of the original system except that the state vector is now β p φ r T . For heading angle control, both skid to turn (STT) and bank to turn can be implemented. In STT, rudders are the primary control surface. Rudders are used to create sideslip and aircraft turns in the direction of rudder deflection. Ailerons can be used to encounter roll produced due to yaw and to keep the wings level if desired. In BTT, ailerons are the primary control surface. Ailerons are deflected to bank the aircraft into a turn. In the current study, a BTT maneuver is modeled. Lateral aerodynamic derivatives estimated using Xflr5 are listed in Table 9. Table 9. Lateral derivatives of aerodynamic coefficients.

Derivatives Value
The lateral dynamics are complex as compared to longitudinal dynamics, involving one force and two moments. Eigenvalues of lateral plant matrix A lat are −26.09, −1.47 + 2.12i, −1.47 − 2.12i and 0.116, which corresponds to roll subsidence, Dutch roll and spiral mode, respectively. The positive value of the spiral mode indicates that it is unstable. Before designing a linear BTT control system, we must first add the BTT equation to our system. The BTT equation is given by: It states that any non-zero roll angle will induce a yaw rate. This equation can be incorporated into state-space using . ψ = r. After augmenting the BTT equation in Equation (41), the plant matrix and control matrix are given by: 28i. All the roots are now negative; hence, the system is now stable. The response to a 110 • turn is shown in Figure 15. The yaw rate is limited to 2 • /s.

Six-DOF Non-Linear Simulation
The proposed control strategy is validated in the 6-DOF Matlab Simulink model. The aerodynamic model used for this simulation is presented below. to an altitude of 350 m at the rate of 2 m/s, after being hand-launched from an initial altitude of 0 m. After climbing to 350 m, there is a small cruise segment. At t = 250 s, the UAV starts a climbing turn, climbing with a rate of 1 m/s and turning with a yaw rate of 2°/s. After climbing to 700 m, the UAV holds this altitude and continues to fly in a circular path with the same yaw rate. As noted in Figure 15, the BTT maneuver also produces a sideslip

Six-DOF Non-Linear Simulation
The proposed control strategy is validated in the 6-DOF Matlab Simulink model. The aerodynamic model used for this simulation is presented below.
C l r r + C l δa δa + C l δr δr N A = C n qSb where C n = C n 0 + C n β β + b 2V C n p p + b 2V C n r r + C n δa δa + C n δr δr where L, D, M A , L A , N A , Y, q and c are lift, drag, pitching moment, rolling moment, yawing moment, side force, dynamic pressure and mean aerodynamic chord, respectively. C L 0 , C D 0 , C m 0 , C l 0 , C n 0 and C Y 0 are lift, drag, pitching, rolling, yawing and side force coefficient at α = δe = β = δa = δr = 0, respectively. The final trajectory is to climb to an altitude of 350 m at the rate of 2 m/s, after being hand-launched from an initial altitude of 0 m. After climbing to 350 m, there is a small cruise segment. At t = 250 s, the UAV starts a climbing turn, climbing with a rate of 1 m/s and turning with a yaw rate of 2 • /s. After climbing to 700 m, the UAV holds this altitude and continues to fly in a circular path with the same yaw rate. As noted in Figure 15, the BTT maneuver also produces a sideslip angle. In the 6-DOF simulation, the rudder is used to encounter the sideslip angle using a proportional controller, often known as a coordinated turn. Throughout the trajectory, the velocity of 8.5 m/s is maintained. To account for any expected sensor noises, the Matlab Simulink built-in band limited white noise model is used (Figure 16). The angle of attack and sideslip angle are shown in Figure 17a. At the final cruise altitude of 700 m, alpha is around 3.2°, which is consistent with static analysis and design. Sideslip is zero throughout the trajectory. In Figure 17b, x and y coordinates are plotted to visualize the circular motion. For the given yaw rate of 2°/s, the UAV flies in a circular path with 480 m diameter. The height profile of the solar UAV with axial location is also shown in Figure 18a. After reaching the target height of 700 m, the UAV successfully maintains this altitude. Euler angles are also presented in Figure 18b. During climbing turn, the pitch angle is around 10°. After achieving a target height, the pitch angle is reduced to 3.2°, which is equal to the angle of attack. Hence, the solar UAV is flying with zero flight path angle. For the yaw angle, the range of the plot is −180° to 180° (one complete circle). To elaborate, the trajectory of the solar UAV is shown in a 3D plot ( Figure 19). Different segments of the trajectory are shown with different colors. The angle of attack and sideslip angle are shown in Figure 17a. At the final cruise altitude of 700 m, alpha is around 3.2°, which is consistent with static analysis and design. Sideslip is zero throughout the trajectory. In Figure 17b, x and y coordinates are plotted to visualize the circular motion. For the given yaw rate of 2°/s, the UAV flies in a circular path with 480 m diameter. The height profile of the solar UAV with axial location is also shown in Figure 18a. After reaching the target height of 700 m, the UAV successfully maintains this altitude. Euler angles are also presented in Figure 18b. During climbing turn, the pitch angle is around 10°. After achieving a target height, the pitch angle is reduced to 3.2°, which is equal to the angle of attack. Hence, the solar UAV is flying with zero flight path angle. For the yaw angle, the range of the plot is −180° to 180° (one complete circle). To elaborate, the trajectory of the solar UAV is shown in a 3D plot ( Figure 19). Different segments of the trajectory are shown with different colors. The height profile of the solar UAV with axial location is also shown in Figure 18a. After reaching the target height of 700 m, the UAV successfully maintains this altitude. Euler angles are also presented in Figure 18b. During climbing turn, the pitch angle is around 10 • . After achieving a target height, the pitch angle is reduced to 3.2 • , which is equal to the angle of attack. Hence, the solar UAV is flying with zero flight path angle. For the yaw angle, the range of the plot is −180 • to 180 • (one complete circle). To elaborate, the trajectory of the solar UAV is shown in a 3D plot ( Figure 19). Different segments of the trajectory are shown with different colors.

Conclusions
In this paper, the genetic algorithm is used to optimize a solar-powered UAV for a specific height, cruise speed and static margin. The wing airfoil is also considered as a design variable and a new structural mass estimation model is proposed. This model also caters for the change in structural mass with the change in wing airfoil during optimization. The objective is to design an airfoil for specified flight conditions, considering overall performance and weight estimation. To add robustness in the design, an extra battery is added as the system starts to consume the battery as soon as the required power is lower than the available solar power. The optimized configuration has a total mass of 7.08 kg with a battery mass of 3.4 kg. The optimized airfoil is very smooth, with a thickness and camber of 9.7% and 4.6%, respectively. To ensure control and handling, dynamic stability and control systems are also discussed using linearized equations of motion. The control laws are developed using an LQR and PID controller. The designed altitude and BTT controllers are implemented in a closed-loop nonlinear 6-DOF simulation. The results of the 6-DOF simulation are in good agreement with linear static analysis, validating the complete design process.

Conclusions
In this paper, the genetic algorithm is used to optimize a solar-powered UAV for a specific height, cruise speed and static margin. The wing airfoil is also considered as a design variable and a new structural mass estimation model is proposed. This model also caters for the change in structural mass with the change in wing airfoil during optimization. The objective is to design an airfoil for specified flight conditions, considering overall performance and weight estimation. To add robustness in the design, an extra battery is added as the system starts to consume the battery as soon as the required power is lower than the available solar power. The optimized configuration has a total mass of 7.08 kg with a battery mass of 3.4 kg. The optimized airfoil is very smooth, with a thickness and camber of 9.7% and 4.6%, respectively. To ensure control and handling, dynamic stability and control systems are also discussed using linearized equations of motion. The control laws are developed using an LQR and PID controller. The designed altitude and BTT controllers are implemented in a closed-loop nonlinear 6-DOF simulation. The results of the 6-DOF simulation are in good agreement with linear static analysis, validating the complete design process.

Conclusions
In this paper, the genetic algorithm is used to optimize a solar-powered UAV for a specific height, cruise speed and static margin. The wing airfoil is also considered as a design variable and a new structural mass estimation model is proposed. This model also caters for the change in structural mass with the change in wing airfoil during optimization. The objective is to design an airfoil for specified flight conditions, considering overall performance and weight estimation. To add robustness in the design, an extra battery is added as the system starts to consume the battery as soon as the required power is lower than the available solar power. The optimized configuration has a total mass of 7.08 kg with a battery mass of 3.4 kg. The optimized airfoil is very smooth, with a thickness and camber of 9.7% and 4.6%, respectively. To ensure control and handling, dynamic stability and control systems are also discussed using linearized equations of motion. The control laws are developed using an LQR and PID controller. The designed altitude and BTT controllers are implemented in a closed-loop nonlinear 6-DOF simulation. The results of the 6-DOF simulation are in good agreement with linear static analysis, validating the complete design process.