Flight Strategy Optimization for High-Altitude Solar-Powered Aircraft Based on Gravity Energy Reserving and Mission Altitude

Abstract

High-altitude long-duration (HALE) flight capability is one of the ultimate goals pursued by human aviation technology, and the high-altitude solar-powered aircraft (SPA) is the most promising technical approach to achieve this target as well as wide application prospects. Due to the particularity of the energy system, the flight strategy optimization through the storage of gravity potential energy and other methods is a significant way to enhance the flight and application abilities for the SPA. In this study, a flight strategy optimization model has been proposed for the aim of HALE flight capability, which is based on the gravity energy reserving and mission altitude in practical engineering applications. This integrated model contains the five flight path phase model, the three-dimensional kinematic model, aerodynamic model, solar irradiation model and energy store and loss model. To solve the optimization problem of three-dimensional flight strategy, the Gauss pseudo-spectral Method (GPM) was employed to establish and calculate the optimal target as its advantages in treating process constraints and terminal constraints for the multiphase optimization problem. At last, the flight trajectory optimization with minimal battery mass for Zephyr 7 was studied by the GPOPS with some function files in MATLAB. The results indicate that the Zephyr 7 can reduce the battery mass from 16 kg to 12.61 kg for the day and night cycle flight and missions, which equals to increasing the battery specific energy by 23.1%. Meanwhile, the optimization results also show that the attitude angel may contribute to increasing the energy gained by photovoltaic cells. In addition, this optimized flight strategy brings the possibility of monthly or annual continuous flight for SPA as the simulation date is set to the autumnal day.

1. Introduction

Aircraft that use solar radiation as propulsive energy are called solar-powered aircraft (SPA). As a new energy aircraft, it has irreplaceable advantages than the conventional aircraft [1,2]. The SPA can flight higher, with a long time cruise, covering a wide area. In August 2001, Helios achieved the first goal with an unofficial world-record altitude of a 29524 m, 40 min flight above 29261m near Hawaii [3]. In June 2018, the “Zephyr-S” flew a record 25 d 23 h 57 min without interrupted [4]. Therefore, high-altitude SPA with annual continuous flight capability in area residency can be widely applied in the natural monitoring, planetary exploration, border patrol, field investigations, communications relay and electronic warfare. As a result, the high-altitude SPA is definitely a research hotspot for world aviation countries in science and technology [5,6,7,8].

Over the past decade, the study of HALE SPA is mainly focused on the overall parameter design methods, the solar cell technology and applications, efficient aerodynamic layout, the design of internal energy storage systems and energy manage system, the structure design for high aspect ratio wing, et al. [8,9,10,11,12,13,14], and these researches have made significant contributions to the development of SPA. Simultaneously, some HALE flight strategy studies for SPA are also been studied relatively as the special nature of energy systems. As shown in these works, the research of flight strategy can be carried out in two major areas, and one of which is the SPA was considered in particle. Then the SPA can use gravitational potential, wind shear and other methods [2] as the flight strategy to reduce the weight of the energy storage battery. Brandt et al. [15] first proposed the conception of gravity energy storage to reduce the battery mass and did not give a detailed analytical model or quantitative conclusions. G Sachs et al. [16] carried out the optimized flight trajectory for SPA so that it could achieve the day and night flight theoretically without battery for the first time, and the SPA need to glide from 20 km to 2 km. Langelann and Roy [17] describe the method of dynamic gliding that can improve the flight time of unmanned aerial vehicles and expand application in Science journal. The same authors also design a method for autonomously sensing the gradient wind field intensity in a dynamic process [18], which is of great significance for SPA to use wind shear potential energy by gliding. Besides, Gao et al. [2,19] who based on the basic problems of HALE flight and dynamics simulations have proposed the method of energy storage by wind shear, and this idea and method have great significance for further research. Gao et al. also explored the equivalence of energy storage battery and gravity energy storage as well as the energy management strategy, showing that energy storage in altitude is most efficient when the initial altitude is below 10km [2,12].

On the other hand, the SPA can increase solar radiation in the normal direction of wing for photovoltaic cells by adjusting the attitude. Andrew et al. [11], with the objective of maximizing surplus energy, have studied the relationship between the hovering and the max remaining energy for the SPA at fixed altitude. The final result shows that the remaining electric energy obtained by the optimized hovering flight trajectory is much larger than the horizontal straight flight trajectory. Sara et al. [20] have further expanded the flight space and researched the three-dimensional cylindrical space for the relations between the remaining energy and flight trajectory, which has given the optimization method for designing the flight path and some similar qualitative conclusions. Huang et al. [9] also study cylindrical space path planning as the same kinetic equation model as Sara and indicate that the net energy can increase by 70% of the horizontal circular flight in a recent research. Therefore, flight strategy optimizations in these two major areas have made significant contributions to the development for SPA.

Currently, researches on flight strategy optimization that combined the altitude and attitude optimizations have began to attract the interest of research scholars. The first most similar completed study on this topic was carried out by Hosseini et al. [21] in 2016. The author establishes a system for HALE flight that uses the altitude and attitude change to store and account more energy, and the final state of charge of the SPA battery has been maximized. R. Abraham Martin et al. [22] have proposed a simulation to maximize the total remaining energy in the system while the aircraft stay within a 3 km mission radius as well as meeting other system constraints, showing an 8.2% increase in end-of-day battery energy for the most limiting flight condition of the winter solstice. WANG et al. [23] also investigate the three-dimensional flight strategy optimization without the mission area constraint, compared with constant-altitude constant-velocity flight strategy, finding that the optimization flight strategy has more profits when the flight conditions are in lower sunlight intensity at a shorter daytime. However, none of these studies have a higher altitude level flight process in engineering applications, which is generally required for the energy storage battery smooth charging and task execution for the advantage of wind field environment in higher altitude. Besides, the SPA need to cruise in a lower altitude at night based on the gravity energy reserving and lower power consumption. Therefore, two different level flight altitudes including the higher mission altitude at daytime and the lower cruising altitude at night should be designed in the HALE flight strategy for SPA, and this flight strategy is hardly involved in published literatures.

This paper was inspired by those considerations and an integrated flight strategy optimization model was proposed to solve the problem of the three-dimensional flight strategy optimization for the HALE SPA, which is based on gravity energy reserving and mission altitude in engineering application. In this study, the integrated model includes: the flight path phase model, the three-dimensional kinematic model, aerodynamic model, solar irradiation model and energy store and loss model. To solve this nonlinear strategy optimization control problem with dynamic constraints and boundary constraints, the Gauss pseudo-spectra method (GPM) is employed to obtain the minimal energy consumed flight strategy with minimum battery in given configuration of SPA. The contributions of this paper can be summarized as follows.

• A flight path phase model based on gravity energy reserving and mission altitude is proposed for engineering applications and day and night cycle flight.
• In order to calculate the three-dimensional solar energy collection, an energy collection model which contains geolocation, height, aircraft attitude angle is established.
• Due to the mission altitude level flight of the high-attitude SPA, this study shows that the mission altitude may be more beneficial to engineering applications for multiple task execution, like circle hover and Regional residency for natural monitoring, field investigations, communications relay and so on.
• A unique family of solution model of three-dimensional flight strategy optimization for SPA with minimal battery mass is proposed by GPM.
• Due to the particularity of the solar aircraft energy harvesting system, the study results also indicate that a three-dimensional flight strategy optimization may be more beneficial for gaining more solar energy.

2. Models and Methods

The SPA energy system [23,24] is shown in Figure 1. During the daytime, the aircraft generates electricity through photovoltaic cells, and then directly supplies the motor through the energy management system to drive the propellers. Meanwhile, the excess solar energy is stored in battery. At night, the battery powers the motor and on-board electronics to ensure that the SPA can fly uninterrupted at day and night.

Based on this energy system, the models and methods used in this study are established and presented below.

2.1. The Flight Path Phase Model

Based on the gravity energy reserving and mission altitude at daytime for engineering applications, the SPA HALE flight path is divided into five phases, as shown in Figure 2 below. It includes two level flight phases, two ascent phases and a gravity gliding phase. The two level flight phases in this study are defined as cruising altitude at night and mission altitude at daytime, which are marked $h_1$ and $h_2$.

Firstly, the establishment of mission altitude $h_2$ is beneficial to practical engineering applications, which is more conducive to the smooth charging of the energy storage battery. Besides, the mission altitude also provides possibilities for the SPA to execute a flight mission in a stable process, in which the wind speed is lower. Generally, the flight mission could be straight back, circle hover and ra egional mission for monitoring and communication, as shown in Figure 2. Secondly, as we known, the SPA cruising at the altitude $h_1$ in night means lower power consumption. Meanwhile, some tasks such as detection can be also completed at this altitude. Therefore, the above-proposed SPA HALE flight path is very suitable for a day-high and night-low flight strategy and engineering task application, which is different from the three-stage gravity energy storage flight strategy. The details of the phases are as follows.

Phase 1: Maximum power ascending

This flight phase starts at time $t_0$ when the solar radiation energy acquired by the HALE SPA at level flight altitude $h_1$ is sufficient to support the climb flight. The phase ends at time $t_1$ when the flight mission altitude $h_2$ is reached. The flight power supply mode of this phase is: all the solar energy obtained is used for the climbing movement. Therefore, the electric power converted by the energy management system is all available for the aircraft.

Inequality (including greater than or equal to $\ge$ and less than or equal to $\le$) path constraints (IPC), in the optimal control problem usually refers to the constraints of certain functions in a certain path process (“a phase”). Generally speaking, in addition to control variables, time and static parameters are also within the scope of these constraint functions. As a result, the SPA flight power $P$ and the change rate $q$ of the energy storage battery in the Phase 1 are as follows:

$$\{\begin{array}{l}P=I·S_{sc}{\eta }_{MPPT}{\eta }_{mot}{\eta }_{plr}{\eta }_{SC}=T·v\\ \frac{dQ}{dt}=q=0\end{array}$$

(1)

where $P$ is the SPA flight power, $q$ the rate of the battery, $T$ the thrust of propulsion system, $v$ the aircraft speed, $Q$ the battery capacity, $S_{SC}$ the photovoltaic cell area, ${\eta }_{MPPT}$ the efficiency of Maximum Power Point Tracking (MPPT) device, ${\eta }_{mot}$ the efficiency of aircraft motor, ${\eta }_{plr}$ the efficiency of aircraft propeller, ${\eta }_{SC}$ the efficiency of photovoltaic cell. In addition, the photovoltaic cell area $S_{SC}$ can be obtained by the wing area as follows:

$$S_{SC}=k_W·S_W$$

(2)

where $S_W$ is the wing area, $k_W$ the pavement ratio of photovoltaic cell.

Phase 2: Level flight for mission and charging battery

The phase begins at time $t_1$ when the predetermined mission altitude $h_2$ is reached, and ends at time $t_2$ when the energy storage battery is fully charged. In the process, the SPA flies at the altitude $h_2$ and charges the energy storage battery. The flight power supply mode of this phase is: part of the electrical energy converted by photovoltaic cells is used to maintain the aircraft’s flight, and the rest is used to charge the energy storage battery. In this phase, the IPC are as follows:

$$\frac{dQ}{dt}=q=\left(I·S_{SC}{\eta }_{SC}{\eta }_{MPPT}-\frac{P}{{\eta }_{mot}{\eta }_{plr}}\right)·{\eta }_{batin}$$

(3)

$$\{\begin{array}{l}\frac{d\gamma }{dt}=0\\ \frac{dh}{dt}=0\end{array}$$

(4)

where $I$ is the solar radiation intensity, ${\eta }_{batin}$ the charging efficiency of energy storage battery, $\gamma$ the aircraft flight path angle, $h$ the aircraft flight altitude.

Phase 3: Maximum power ascending for gravity energy storage

This flight phase starts at time $t_2$ when the energy storage battery is brimming, and ends at time $t_3$ when the solar irradiation is zero. This phase includes maximum power climbing and powered gravity gliding and depends on whether the solar radiation intensity can meet the power consumption required for flight at its altitude. Before the solar radiation disappears in this process, the maximum flying altitude that the SPA can reach relies on the obtained solar radiation energy. As a result, the flight phase can make full use of the excess solar energy for gravity energy storage. The flight power supply mode of this phase is: all the electrical energy converted by solar cells is used to power the aircraft. Therefore, the IPC of the aircraft in this phase are the same as phase 1, as follows:

$$\{\begin{array}{l}P=I·S_{sc}{\eta }_{MPPT}{\eta }_{mot}{\eta }_{plr}{\eta }_{SC}=T·v\\ \frac{dQ}{dt}=q=0\end{array}$$

(5)

Phase 4: Longest gliding endurance path

The phase begins at time $t_3$ when there is no irradiance, and ends at time $t_4$ when the SPA glides to the cruising altitude $h_1$ at night. Throughout the whole stage, the SPA reduces altitude by gliding without thrust. The SPA energy system does not provide any electricity to the aircraft during this flight phase. In the phase, the IPC of SPA satisfy:

$$\{\begin{array}{l}T=0\\ \frac{dQ}{dt}=q=0\end{array}$$

(6)

Phase 5: Low-power level flight by discharging battery at night

This flight phase starts at time $t_4$ when the SPA reaches the cruising altitude $h_1$, and ends at time $t_5$ when the solar radiation energy acquired by the SPA at altitude $h_1$ is sufficient to support the climbing flight to the next day phase 1. Basically, $t_0$ is equal to $t_5$. During this stage, the SPA only relies on energy storage batteries for cruising. The power supply mode is: the energy storage battery provides all power consumption during the night. Thence, the IPC of the aircraft are as follows:

$$\frac{dQ}{dt}=q=-\frac{P}{{\eta }_{mot}{\eta }_{plr}{\eta }_{batout}}$$

(7)

$$\{\begin{array}{l}\frac{d\gamma }{dt}=0\\ \frac{dh}{dt}=0\end{array}$$

(8)

where ${\eta }_{batout}$ is the discharging efficiency of energy storage battery.

2.2. Aerodynamic Model

Generally, the aerodynamic forces acting on an aircraft can be calculated by the following formula:

$$\{\begin{array}{l}L=\frac{1}{2}\rho v^2·S_W{C}_L\\ D=\frac{1}{2}\rho v^2·S_W{C}_D\end{array}$$

(9)

where $L$ is the lift of aircraft, the $D$ drag force of aircraft, $C_L$ and $C_D$ are the lift and drag coefficients, $\rho$ is the atmospheric density. In this paper, the $\rho$ value on altitude use the relevant section of the 1976 Standard Atmosphere, as described in reference [25,26].

As we known, the HALE SPA usually has the characteristics of large aspect ratio and its wing are larger than the fuselage. Therefore, $C_L$ and $C_D$ of this type aircraft are closely related to airfoil profile, angle of attack and Reynolds number. In order to simplified calculation, the $C_L$ and $C_D$ of SPA can have the following formulas:

$$\{\begin{array}{l}{C}_L=0.9C_{L0}\\ C_D=C_{D0}+C_{Dp}+C_{Di}\\ C_{Di}=\frac{C_L^2}{e\pi AR}\\ e=1.78(1-0.045A{R}^{0.68})-0.64\end{array}$$

(10)

where $C_{L0}$ is the lift coefficient of airfoil, $C_{D0}$, $C_{Dp}$ and $C_{Di}$ are the coefficients of the airfoil zero-lift drag, the parasitic drag and the lift-induced drag, $e$ is the Oswald ’s coefficient, $AR$ is the aspect ratio of SPA.

In this study, the Zephyr 7 is the research object. Hence, the airfoil coefficient of aerodynamic characteristic refers to FX 63-137, 13.7% smoothed, Re = 20000 and the parasitic resistance coefficient $C_{Dp}$ is 0.005. The detailed data is shown in the

Table 1. The airfoil coefficient of lift and drag.

$\alpha$(°)	−6	−5	−4	−3	−2	−1	0	1	2	3
$C_{L0}$	0.0789	0.226	0.3815	0.5343	0.6464	0.7596	0.8763	0.9906	1.0971	1.2013
$C_{D0}$	0.0344	0.0204	0.0162	0.0132	0.014	0.0145	0.0146	0.0148	0.0149	0.0152
$\alpha$(°)	4	5	6	7	8	9	10	11	12	13
$C_{L0}$	1.3039	1.3961	1.4721	1.5644	1.6309	1.6707	1.6853	1.6725	1.6659	1.6776
$C_{D0}$	0.0156	0.016	0.0162	0.0177	0.0196	0.0224	0.0272	0.0357	0.0463	0.0569

.

2.3. Aircraft Kinematic Model in Three Dimensional

For the research of flight strategy optimization, the focus of SPA during the flight is on the macroscopic characteristics of the trajectory. Therefore, a simplified particle model can be used for research [16]. As shown in Figure 3, the three-dimensional dynamic equation of SPA can be expressed as follows:

$$\{\begin{array}{c}\begin{array}{l}\frac{dv}{dt}=\frac{T}{m}-\frac{D}{m}-g·\sin \gamma \\ \frac{d\gamma }{dt}=\frac{L}{mv}\cos \psi -\frac{g}{v}\cos \gamma \end{array}\\ \frac{d\phi }{dt}=\frac{L}{mv\cos \gamma }\sin \psi \\ \frac{dx}{dt}=v\cos \gamma \cos \phi \\ \frac{dy}{dt}=v\cos \gamma \sin \phi \\ \frac{dh}{dt}=v\sin \gamma \end{array}$$

(11)

where $x$, $y$ and $h$ are the SPA position coordinates in flat earth-fixed frame, $m$ the total mass of SPA, $g$ the acceleration of gravity, $\psi$ the aircraft roll angle, $\phi$ the aircraft yaw angle.

2.4. Solar Irradiation Model and Mission Conditions

During the flight, with the change of flight altitude and attitude, the intensity of solar radiation received by photovoltaic cells on the wing surface of the aircraft will change accordingly. As the main power source of SPA, the amount of solar energy obtained will directly affect its flight performance. Therefore, a relatively accurate solar flux model needs to be established to calculate the energy collection for SPA. In this study, based on the solar radiation intensity calculation model in the earth-fixed frame, the integration of coordinate conversion principle and the flight attitude angle is used to obtain the solar radiation intensity for SPA in the body coordinate system. The details are as follows.

2.4.1. Calculation Model of Solar Radiation Intensity in Ground Coordinate System

This paper defines the earth-fixed frame as $O_g{X}_g{Y}_g{Z}_g$, where $O_g{X}_g$, $O_g{Y}_g$, $O_g{Z}_g$ are expressed out the north, west, and vertical upward through the ground, as illustrated in Figure 3. Then in the earth-fixed frame, the solar radiation intensities $I_{Xg}$, $I_{Yg}$, and $I_{Zg}$ (positive and negative represent directions) in each direction are the following formulas:

$$\{\begin{array}{l}{I}_{Xg}=I_{SC}·{\xi }_0·P_{CO}·\cos H_S·\cos A_S\\ I_{Yg}=I_{SC}·{\xi }_0·P_{CO}·\cos H_S·\sin A_S\\ I_{Zg}=I_{SC}·{\xi }_0·P_{CO}·\sin H_S\end{array}$$

(12)

where $I_{SC}$ is the solar irradiance constant 1367 W/m², ${\xi }_0$ the eccentricity correction coefficient related to the distance between the sun and the earth, $P_{CO}$ the atmospheric transparency coefficient, $H_S$ the sun altitude angle, $A_S$ the sun azimuth. Those parameters can be calculated by the formulas in literature [27,28], which are not be described in detail here.

2.4.2. Calculation Model of Solar Radiation Intensity in the Body Coordinate System

This paper defines the aircraft body-fixed frame as $O_b{X}_b{Y}_b{Z}_b$, where $O_b{X}_b$ is located in the symmetry of the aircraft and points forward parallel to the fuselage axis, $O_b{Y}_b$ points to the left perpendicular to the symmetry plane, $O_b{Z}_b$ is expressed out to point upwards and perpendicular to $O_b{X}_b$ and $O_b{Y}_b$, as shown in Figure 3. Therefore, the solar radiation intensities $I_{Xb}$, $I_{Yb}$, and $I_{Zb}$ received in aircraft body-fixed frame can be obtained by the attitude angle (yaw angle $\phi$, pitch angle $\gamma$, roll angle $\psi$ ) through coordinate transformation, as shown in Equation (12):

$$\left[\begin{array}{c}{I}_{Xb}\\ I_{Yb}\\ I_{Zb}\end{array}\right]=\left[\begin{array}{ccc}\cos \gamma \cos \phi & \cos \gamma \sin \phi & -\sin \gamma \\ \sin \gamma \sin \psi \cos \phi -\cos \psi \sin \phi & \sin \gamma \sin \psi \sin \phi +\cos \psi \cos \phi & \sin \psi \cos \gamma \\ \sin \gamma \cos \psi \cos \phi +\sin \psi \sin \phi -\sin \psi \cos \phi & \sin \gamma \cos \psi \sin \phi & \cos \psi \cos \gamma \end{array}\right]\left[\begin{array}{c}{I}_{Xg}\\ I_{Yg}\\ I_{Zg}\end{array}\right]$$

(13)

Consequently, a solar irradiation model is established in a MATLAB file for later model integration, which can calculate solar radiation intensity with the altitude, date, longitude, latitude, moment, roll angle, pitch angle and yaw angle.

2.5. Problem Optimization Frame

Obviously, the goal of three-dimensional flight strategy optimization for the HALE SPA based on gravity energy reserving and mission altitude is a multi-phase optimization problem, and can be disguised as a nonlinear programming problem. Figure 4 shows the connection of this multi-phase optimization problem. Therefore, the optimization objective, boundary conditions and dynamic equations for each segment need to be analyzed.

The boundary conditions and the optimization objective of HALE SPA for each phase are summarized as bellows:

Phase 1: Maximum power ascending($t_0$-$t_1$)

As described in Section 2.1, the main task of the aircraft at this phase is to climb to the mission altitude with the maximum available power, so the performance indicators of the aircraft can be set as follows:

$$Min{J}_1=t_1-t_0$$

(14)

The three-dimensional dynamic equation of SPA can be expressed as same as Equation (9), and the applied boundary conditions can be as follows:

$$\{\begin{array}{cc}{v}_{\left\{1\right\}}(t_0)=v_0& {\gamma }_{\left\{1\right\}}(t_f)=0\\ {\phi }_{\left\{1\right\}}(t_0)={\phi }_0& h_{\left\{1\right\}}(t_f)=h_2\\ {\gamma }_{\left\{1\right\}}(t_0)=0& Q_{\left\{1\right\}}(t_f)=Q_{\min }\\ x_{\left\{1\right\}}(t_0)=x_0& \\ \begin{array}{l}{y}_{\left\{1\right\}}(t_0)=y_0\\ h_{\left\{1\right\}}(t_0)=h_1\\ Q_{\left\{1\right\}}(t_0)=Q_{\min }\end{array}& \end{array}$$

(15)

Phase 2: Level flight for mission and charging battery($t_1$-$t_2$)

The aim of phase 2 is the task execution in predetermined mission altitude as well as charge battery, so the power consumption of the level flight is concerned at this phase. Hence, the target optimization function can be expressed as per the following formula:

$$Min{J}_2={\displaystyle {\int }_{t_2}^{t_1}{P}_{level}dt=}{\displaystyle {\int }_{t_2}^{t_1}{T}_{h_2}vdt}$$

(16)

The applied boundary conditions satisfy:

$$\{\begin{array}{cc}{v}_{\left\{2\right\}}(t_0)=v_{\left\{1\right\}}(t_f)& {\gamma }_{\left\{2\right\}}(t_f)=0\\ {\phi }_{\left\{2\right\}}(t_0)={\phi }_{\left\{1\right\}}(t_f)& h_{\left\{2\right\}}(t_f)=h_2\\ {\gamma }_{\left\{2\right\}}(t_0)={\gamma }_{\left\{1\right\}}(t_f)=0& Q_{\left\{2\right\}}(t_f)=Q_{\max }\\ \begin{array}{c}{x}_{\left\{2\right\}}(t_0)=x_{\left\{1\right\}}(t_f)\\ y_{\left\{2\right\}}(t_0)=y_{\left\{1\right\}}(t_f)\end{array}& \\ \begin{array}{c}{h}_{\left\{2\right\}}(t_0)=h_{\left\{1\right\}}(t_f)=h_2\\ Q_{\left\{2\right\}}(t_0)=Q_{\min }\end{array}& \end{array}$$

(17)

Phase 3: Maximum power ascending for gravity energy storage($t_2$-$t_3$)

In the phase 3, the final altitude is the key to the optimization. The higher altitude the SPA gains, the greater advantage for the phase 4. Therefore, the performance indicators can be set as follows:

$$Max{J}_3=h_{\left\{3\right\}}(t_f)$$

(18)

The applied boundary conditions in this phase can be expressed as follows:

$$\{\begin{array}{cc}{v}_{\left\{3\right\}}(t_0)=v_{\left\{2\right\}}(t_f)& \\ {\phi }_{\left\{3\right\}}(t_0)={\phi }_{\left\{2\right\}}(t_f)& h_{\left\{3\right\}}(t_f)\ge h_2\\ {\gamma }_{\left\{3\right\}}(t_0)={\gamma }_{\left\{2\right\}}(t_f)=0& Q_{\left\{3\right\}}(t_f)=Q_{\max }\\ x_{\left\{3\right\}}(t_0)=x_{\left\{2\right\}}(t_f)& \\ \begin{array}{c}{y}_{\left\{3\right\}}(t_0)=y_{\left\{2\right\}}(t_f)\\ h_{\left\{3\right\}}(t_0)=h_{\left\{2\right\}}(t_f)=h_2\\ Q_{\left\{3\right\}}(t_0)=Q_{\max }\end{array}& \end{array}$$

(19)

Phase 4: Longest gliding endurance path($t_3$-$t_4$)

The aim of the aircraft at this phase is to achieve gliding for longer time, which is helpful for delaying the moment of battery discharge and the night flight for SPA. Hence, the target optimization function can be expressed as follow:

$$Max{J}_4=t_4-t_3$$

(20)

The applied boundary conditions in phase 4 satisfy:

$$\{\begin{array}{cc}{v}_{\left\{4\right\}}(t_0)=v_{\left\{3\right\}}(t_f)& {\gamma }_{\left\{4\right\}}(t_f)=0\\ {\phi }_{\left\{4\right\}}(t_0)={\phi }_{\left\{3\right\}}(t_f)& h_{\left\{4\right\}}(t_f)=h_1\\ {\gamma }_{\left\{4\right\}}(t_0)={\gamma }_{\left\{3\right\}}(t_f)& Q_{\left\{4\right\}}(t_f)=Q_{\max }\\ x_{\left\{4\right\}}(t_0)=x_{\left\{3\right\}}(t_f)& \\ \begin{array}{l}{y}_{\left\{4\right\}}(t_0)=y_{\left\{3\right\}}(t_f)\\ h_{\left\{4\right\}}(t_0)=h_{\left\{4\right\}}(t_f)\\ Q_{\left\{4\right\}}(t_0)=Q_{\max }\end{array}& \end{array}$$

(21)

Phase 5: Low-power level flight by discharging battery at night($t_4$-$t_5$)

The aim for this phase is to keep the SPA aloft in the cruising altitude $h_1$ with minimum power consumed. Hence, the performance indicators can be expressed as per the following formula:

$$Min{J}_5={\displaystyle {\int }_{t_5}^{t_4}{P}_{level}dt=}{\displaystyle {\int }_{t_5}^{t_4}{T}_{h_1}vdt}$$

(22)

The applied boundary conditions in this phase satisfy:

$$\{\begin{array}{cc}{v}_{\left\{5\right\}}(t_0)=v_{\left\{4\right\}}(t_f)& {\gamma }_{\left\{5\right\}}(t_f)=0\\ {\phi }_{\left\{5\right\}}(t_0)={\phi }_{\left\{4\right\}}(t_f)& h_{\left\{5\right\}}(t_f)=h_1\\ {\gamma }_{\left\{5\right\}}(t_0)={\gamma }_{\left\{4\right\}}(t_f)=0& Q_{\left\{5\right\}}(t_f)=Q_{\min }\\ x_{\left\{5\right\}}(t_0)=x_{\left\{4\right\}}(t_f)& \\ \begin{array}{c}{y}_{\left\{5\right\}}(t_0)=y_{\left\{4\right\}}(t_f)\\ h_{\left\{5\right\}}(t_0)=h_{\left\{4\right\}}(t_f)=h_1\\ Q_{\left\{5\right\}}(t_0)=Q_{\max }\end{array}& \end{array}$$

(23)

To summarize, the total target optimization function can be expressed as follow:

$$\begin{array}{l}MinJ=J_1+J_2-J_3-J_4+J_5\\ =(-t_0+t_1+t_3-t_4)+{\displaystyle {\int }_{t_1}^{t_2}{P}_{level}dt}+{\displaystyle {\int }_{t_5}^{t_4}{P}_{level}dt}-h_{\left\{3\right\}}(t_f)\end{array}$$

(24)

In addition, the day and night cycle boundary conditions also satisfy the following formulas in order that the HALE SPA can perform monthly flight mission.

$$\{\begin{array}{l}{v}_{\left\{5\right\}}(t_f)=v_{\left\{1\right\}}(t_0)\\ {\phi }_{\left\{5\right\}}(t_f)={\phi }_{\left\{1\right\}}(t_0)\\ {\gamma }_{\left\{5\right\}}(t_f)={\gamma }_{\left\{1\right\}}(t_0)\\ x_{\left\{5\right\}}(t_f)=x_{\left\{1\right\}}(t_0)\\ y_{\left\{5\right\}}(t_f)=y_{\left\{1\right\}}(t_0)\\ h_{\left\{5\right\}}(t_f)=h_{\left\{1\right\}}(t_0)\\ Q_{\left\{5\right\}}(t_f)=Q_{\max }\\ t_5-t_0=24\mathit{hours}\end{array}$$

(25)

2.6. Optimization Approach and Solving Process

Generally speaking, the trajectory optimization problem can be transformed into a nonlinear numerical optimal solution process by establishing a mathematical model. Some of the original conditions can be transformed into corresponding states and control constraints after transformation. The GPM is a direct global method, which has advantages in calculation efficiency and accuracy. Therefore, the GPM is use to solving this multi-phase optimization problem in this paper and the above optimization control problem can be disguised as a nonlinear programming problem. The basic idea of GPM is the polynomial interpolation of state and control variables in the whole time region [29] so that the continuous optimal control problem can be formally transformed. Hence, the problem is cleverly to become a traditional solution to solve process of nonlinear algebraic equations constrained by IPC, and the GPM solution steps are roughly as follows:

Firstly, the time domain transformation for each phase is modified to transform the time period $t^n\in \left[t_0{}^n-t_f{}^n\right]$ into time interval ${\tau }^n\in \left[-1,1\right]$. Assuming that there are $M$ numbers LG allocation points on each interval ${\tau }^n\in \left[-1,1\right]$, and then those distribution points can be characterized as characteristic roots of Legendre polynomials in $M$ numbers order [30]. As result, the state and control variables can be approximation by $M+1$ order Legendre polynomials. Hence, the dynamic constraint equations can be converted to algebraic constraint equations by the above-approximation. Meanwhile, the boundary constraints also can be transcribed by ${\tau }^n{}_k$, which is the $k$ numbers of the LG points. At the same time, because the performance indicator function does not contain state variables, the total target optimization function still can be expressed by Equation (24).

In summary, the original multi-phase continuous optimization problem for SPA can be transformed into nonlinear programming problem, which is defined by target optimization function in Equation (23) and dynamic constraints in above-mentioned. Besides, the optimization parameters can be obtained by a MATLAB software named Gauss pseudo-spectral Optimization Software (GPOPS) [31] with some reasonable optimality and feasibility tolerances.

In this study, the angle of attack $\alpha$, the roll angel $\psi$ and the thrust $T$ are chosen as three control variable to change the attitude and flight altitude of the SPA. The paper combines the GPOPS and function files of other models in Section 2 to solve the nonlinear programming problem for SPA in MATLAB. The optimization solving process is shown as follows.

Step 1: Initialize the SPA parameters, geographical parameters and mission parameters.

Step 2: Define the state and control variables and set up the limits on variables.

Step 3: Provide boundary and initial values in each phase of the problem, and set up the event constraints which is linked in phases.

Step 4: Provide an initial Mesh for each phase, the phase number is 5 and the LG allocation points number is 40 for each phase in this study.

Step 5: Assemble all information into setup structure, which contains the aircraft kinematic model, aerodynamic model, solar irradiation and day-night cycle boundary conditions, etc., function files.

Step 6: Load the target optimization function and the variable battery mass $m_{bat}$ of SPA for satisfying the $t_5-t_0=24$ h .m function files;

Step 7: Solve the problem by using GPOPS in MATLAB;

Step 8: Summarize the information of each node and phase to obtain the optimal flight trajectory for SPA base on energy reserving and mission altitude.

3. Simulation and Results Discussion

3.1. Simulation Condition

Firstly, the study has tried to analyze the environmental characteristics of the flight area by the above model in Section 2.4 for further parameter settings. Figure 5 display the variation of the maximum solar irradiance intensity with altitude in the south china sea (105°E, 4°N, data = 266). Besides, a realistic wind speed data as a function of the altitude from high altitude detection balloon is also been shown. It can be concluded from the analysis that the maximum solar irradiance intensity will not change drastically with the altitude after reaching 20 km. Meanwhile, the wind speed in the environment from 17.5 km to 27.5 km is relatively low, which is beneficial to mission execution.

After the above preparation works, a series of numerical simulations are computed to verify the proposed method. It will be very efficient to solve the flight trajectory optimization problem into a nonlinear programming problem. The parameters of initial configuration are listed in

Table 2. The parameters of initial configuration for model.

Parameter	Value	Description
$m_{struct}$ (kg)	37	Mass of structure
$m_{bat}$ (kg)	16	Battery mass
$E_{bat}$ (Wh/kg)	350	Energy density of battery
$b$ (m)	22.5	Span length
$S_W$ (m2)	25.3	Wing area
$S_{SC}$ (m2)	20.24	Area of photovoltaic cells
${\eta }_{SC}$ (-)	0.2	efficiency of photovoltaic cell
${\eta }_{MPPT}$ (-)	0.9	Efficiency of MPPT
${\eta }_{mot}$ (-)	0.85	efficiency of aircraft motor
${\eta }_{plr}$ (-)	0.8	efficiency of aircraft propeller
${\eta }_{batin}$ (-)	0.9	Efficiency of battery for charging
${\eta }_{batout}$ (-)	0.9	Efficiency of battery for discharging

, which is referred from Zephyr 7 [2,4,8], as the Zephyr 7 is been introduced comprehensive and detailed in all the literatures among the high-altitude SPA. In addition, the accurate parameters of Zephyr 7 is more conducive to the advantage and reliability of the simulation results by this strategy.

The set of simulation conditions and constraints of variables are listed in

Table 3. The parameters of simulation conditions.

Parameter	Value	Description
$h_1$ (km)	15	Cruising altitude
$h_2$ (km)	23	Battery mass
${\mathit{SOC}}_0$	0.05	Initial battery percentage
$d_n$ (-)	266	Autumn equinox day number
$\mathit{Lati}$ (°)	4	Latitude of the starting point
$\mathit{Long}$ (°)	105	Longitude of starting point
$t_{sr}$ (h)	6.8611	The sunrise time
$t_{noon}$ (h)	12.8564	The standard noon time
$t_{ss}$ (h)	18.8517	The sunset time

and

Table 4. The constraints of variables.

Parameter	Minimum Value	Maximum Value
$h$ (km)	15	30
$v$ (m/s)	0	50
$SO{C}_0$	0.05	0.99
$Lati$ (°)	4	15
$Long$ (°)	105	11
$\alpha$ (°)	5	10
$\gamma$ (°)	-15	15
$\psi$ (°)	-1	1
$T$ (N)	0	1000

, respectively.

3.2. Result and Discussion

In this study, the SPA selects the South China Sea (105°E-116°E, 4°N-15°N) as the target flight area. Besides, the two level flight altitude are 15km and 23km, which is obtained according to the environmental characteristics of the flight area. As described in Section 2.6 above, the study firstly use the GPM method to solve the optimal flight path, and then look for the minimum battery mass that meets long flight time constraints. As results, the final optimized battery mass $m_{bat}$ of Zephyr 7 for satisfying the $t_5-t_0=24$ h in Equation (25) is 12.61 kg, which is shown in MATLAB workspace. Compared with the original initial parameter of battery mass, it is equivalent to a 23.1% reduction for Zephyr 7. It also means how much potential energy can be gained by this flight strategy. Meanwhile, The the three-dimensional trajectory optimization, the speed with altitude, the quantity of electricity, the thrust, the flight path pitch angle with attack angle, the yaw angle and the roll angle with minimal battery mass are shown in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10, respectively.

From the Figure 7a, the time of each phase for this flight strategy is 2.73 h, 5.03 h, 4.09 h, 2.06 h and 10.09 h, respectively. Meanwhile, the speed has good following to the altitude, and the highest flight altitude of Zephyr 7 can reach is about 28.47km at time 18:03. In addition, it can be seen that the flight trajectory of phase 2 from Figure 6 is estimated to have a symmetrical trend at noon time as the same as the trend of yaw angle in Figure 10, which can conclude that the SPA will adjust the trajectory with the solar altitude angle changing symmetrically at the standard noon time. Besides, in order to better analyze the changes of battery charge percentage at noontime, we have performed a simple polynomial interpolation on the data points of the battery percentage at 11:00 to 14:00 by EXCEL. As shown in Figure 8, the battery charge rate has the same trend to the intensity of solar radiation, which is different from the same charge rate in some other studies. More importantly, the unpowered gliding time of the phase 4 is 2.06 h (Figure 8), which did not consume the energy of the battery. This optimization result can be seen as the result of gravity energy reserving, and thus the quality of the battery can be reduced by this strategy.

Another thing has to note is that there is a phenomenon of oscillation at the connection points over the flight time. The literature [32] has specifically studied the oscillation of this numerical result, and the research results show that the calculation results of GPM are stable. Besides, the calculation results of the literature [33] also indicate that the optimal gravity glide track has the characteristics of oscillating with time by GPM. Therefore, the calculation results are reasonable and credible in this paper.

4. Conclusions

In this paper, research on HALE flight trajectory with minimal battery mass for SPA based on gravity energy reserving and mission altitude is conduced to achieve the following.

• A flight path phase model based on gravity energy reserving and mission altitude has been carried out for analysis and research background for engineering applications is taken into consideration for the day and night cycle flight and tasks.
• GPM is feasible for the research on flight trajectory optimization for high-altitude SPA due to its advantages in calculation efficiency and accuracy as well as the integration between GPOPS and function files in MATLAB.
• This flight trajectory optimization results shows that Zephyr 7 can reduce the battery mass from 16 kg to 12.61 kg for the day and night cycle flight and missions, which equals to reducing its total mass by 6.4%.
• According the optimization results, the yaw angel may have more influence on high-altitude SPA 3D flight strategy in energy harvesting than the roll angle and the battery charge rate has good following to the solar irradiation intensity, which is instructive for SPA to use solar radiation reasonably at noon phase in mission execution.

This study focuses on flight strategy optimization of the HALE SPA in three-dimensional, and gives the optimal flight trajectory with minimal battery mass based on gravity energy reserving and mission altitude. It will be very interesting to have the result of anther case in different simulation conditions based on this flight strategy. As a result for the future works, there are two aspects of the related work. Firstly, we will further design mission path planning frame-work for the mission altitude in the integrated model, which is more complicated than the proposed framework. Secondly, the energy management for this flight strategy and the influences of other environmental factors for SPA such as the wind will also be taken into our considerations in future research.