Research on route planning for solar UAV based on the intelligent optimization algorithm

Abstract

The solar unmanned aerial vehicle (UAV) route planning needs to comprehensively consider the conversion efficiency of solar cells under the influence of solar ground reflection radiation and sky scattering radiation. On the one hand, it is necessary to consider the cost of radar threat, mileage energy consumption, mountain impact and other costs. On the other hand, it is also necessary to consider the influence of high threat, mountain shadow occlusion cost, and cloud shading cost on solar photovoltaic conversion efficiency. The above problem was solved through using the ant colony intelligent optimization algorithm. By constructing ant colony paths rationally, models of mountain impact cost, high threat, mountain shadow shelter cost, and cloud shading cost were established. The constraints such as the maximum action distance, solar irradiation angle and effective action distance of various costs were introduced into the cost model and exploration factor calculation, and the comprehensive optimization problem of solar UAV route was solved. Finally, the simulation results show that the algorithm path structure is reasonable; the target node can be found independently; the convergence speed can meet the requirements of route planning; the generated route cost is small; the algorithm is reasonable and effective.

Introduction

Solar drone, as the name implies, is a drone that uses solar energy as the power source for flight and operation, which has the main characteristics of ultra-long endurance, long-term empty space, ultra-high flight, and extensive operation, and its value is comparable to “quasi-satellite,” and it has the advantages of flexible deployment and good economy. In recent years, for the purposes of major natural disaster warning, high-altitude sea area patrol and supervision, emergency rescue and disaster relief, counter-terrorism and stability maintenance, mobile communications and other purposes, countries are investing in, accelerating research and development, and actively deploying the development of solar drones, hoping to seize a new highland in the future drone industry.^1–7 UAV route planning usually only considers the impact of the cost of conventional threats such as radar, missiles, and anti-aircraft artillery, and performs combined optimization. Due to the advantages of small energy consumption, long flight range and long air time, solar drones are stepping up research on this problem and have achieved great development. Therefore, route planning for solar-powered drones is very urgent. Solar UAVs not only need to consider the above threat costs, but also need to combine their own characteristics, comprehensively consider the conversion efficiency of solar cells under the influence of solar ground reflection radiation and sky scattered radiation and other factors, in addition to radar and terrain threats, solar angle, terrain shadow, cloud shading and other factors affecting solar radiation into the calculation model, and finally in order to facilitate analysis, the UAV wing solar cell installation surface is simplified to a horizontal plane, and the solar UAV flies near the earth's equatorial plane. The project team carried out the route planning research of ant colony algorithm for ordinary UAVs in the early stage, and the algorithm has good convergence effect and is relatively easy to project and implement. I am responsible for researching the application of intelligent optimization algorithms to solar route planning problems under specific constraints. To this end, I introduced the ant colony optimization algorithm into the engineering application of solar UAV route planning problem, on the one hand, by reasonably constructing the ant colony path, and establishing the mountain impact cost, high threat, mountain shadow shelter cost, cloud shelter cost model; On the other hand, the constraints such as the maximum action distance, solar irradiation angle and effective action distance of various costs are introduced into the cost model and exploration factor calculation. Finally, the solution of the comprehensive optimization problem of solar UAV route is realized, and the effect is very good.^8–10

Problem description^11–20

Spatial grid generation

3D divides the grid into X, Y, Z directions, where the grid sizes are represented as $X_{\mathrm{Gird}}$ , $Y_{\mathrm{Gird}}$ , $Z_{\mathrm{Gird}}$ , and each grid space is a cube of the same size. In the grid, each node has 26 adjacent nodes. When selecting the path, the next node based on the current node must be selected from 26 adjacent nodes centered on the current node. In addition, the grid size depends on the threat point location distribution and other factors, should not be too small and should not be too large.

Route synthesis optimization

In addition to the cost of radar threat, mileage energy consumption and mountain impact threat, altitude threat, mountain shadow cover cost and cloud cover cost should be considered. The goal of route planning is to minimize the cost of the whole, and according to radar, mountain terrain, cloud distribution and other model constraints, there are some constraints, such as the maximum operating distance, the solar projection angle, the distribution of clouds and the effective killing distance. Among them, the radar has the possibility to detect the position of the UAV, the mountain terrain has the possibility to cause the UAV to be hit, and the cloud cover, the terrain cover will affect the solar cell photoelectric conversion efficiency, it is also about the energy security of drones, so they can all be considered a threat. All of these threats must be taken into account in order to calculate course performance. The threat cost function can be expressed as the weighted sum of the threat cost of radar, altitude, terrain shadow, cloud cover, terrain, energy consumption, etc.^9–15

$$\min W=\min {\int }_0^L[{\delta }_{\mathrm{R}}{w}_{\mathrm{R}}(s)+{\delta }_{\mathrm{H}}{w}_{\mathrm{H}}(s)+{\delta }_{\mathrm{D}}{w}_{\mathrm{D}}(s)+{\delta }_{\mathrm{Y}}{w}_{\mathrm{Y}}(s)+{\delta }_{\mathrm{T}}{w}_{\mathrm{T}}(s)+{\delta }_{\mathrm{l}}{w}_{\mathrm{l}}(s)]ds$$

(1)

In Formula (1), L is the route; W is the optimal objective function; $w_{\mathrm{R}}(s)$ is the route’s radar threat cost; $w_{\mathrm{H}}(s)$ is the high threat cost; $w_{\mathrm{D}}(s)$ is the terrain shadow threat cost; $w_{\mathrm{Y}}(s)$ is the cloud cover threat cost; $w_{\mathrm{T}}(s)$ is the route's terrain threat cost; $w_{\mathrm{l}}(s)$ for mileage. Coefficients ${\delta }_{\mathrm{R}}$ , ${\delta }_{\mathrm{H}}$ , ${\delta }_{\mathrm{D}}$ , ${\delta }_{\mathrm{C}}$ , ${\delta }_{\mathrm{T}}$ , ${\delta }_{\mathrm{l}}$ are the cost weight of each route, and the value is equal to 1. The mileage cost depends on the length of the voyage, the longer the cost is, the smaller the cost. Other threat costs are related to factors such as the distance or angle of the threats.

Establishment of solar radiation-related cost model

Solar energy UAV energy system usually includes two parts: solar battery and storage battery. Therefore, the photovoltaic conversion efficiency of solar cells, shading and other factors will directly affect the energy storage of UAV, and then affect the flight range and even flight safety. The photoelectric conversion efficiency is determined by the cell type. Shadow shadowing, on the other hand, requires the drone to be in a cloud or mountain shadow. Solar radiation comes from direct and scattering, respectively. Scattering can be divided into scattering by ground reflection and sky scattering. Since solar panels are mounted on the upper surface of the wing and no solar cells are mounted on the lower surface, the effect of ground reflected radiation on energy can be ignored.^6–11

$$I_{\mathrm{Gs}}=I_{\mathrm{Ds}}+I_{\mathrm{Is}}$$

(2)

In Formula (2), the $I_{\mathrm{Gs}}$ is the solar radiation received by the solar cell mounted on the wing surface of the solar-powered UAV, and the $I_{\mathrm{Ds}}$ is the solar direct radiation received by the solar cell on the wing of the solar-powered UAV. The $I_{\mathrm{Is}}$ is the solar scattering radiation received by a solar cell on the wing of a solar-powered UAV.^13–19

Direct radiation calculations

Direct radiation $I_{\mathrm{Ds}}$ can be calculated as follows:

$$I_{\mathrm{Ds}}=I_0\cdot \exp \left({\displaystyle \frac{-k}{\sin h}}\right)\cdot \cos i$$

(3)

In formula (3), If $I_0$ is the solar radiation constant, then the normal solar radiation intensity can be calculated as $I_0\cdot \exp (-k/\sin h)$ , i as the solar incidence angle, k as the fitting coefficient, and h as the solar altitude angle. The solar altitude angle h can be calculated according to the following equation:

$$\sin h=\sin \phi \cdot \sin \sigma +\cos \phi \cdot \cos \sigma \cdot \cos \omega$$

(4)

Where $\omega$ , $\sigma$ and $\phi$ represent the solar hour angle, declination angle and latitude, respectively. According to the flight date, position longitude, standard time calculation $\omega$ and $\sigma$ .

In addition, $\cos i$ is calculated according to the following equation:

$$\cos i=\cos \theta \cdot \sin h+\sin \theta \cdot \cos h\cdot \cos \epsilon$$

(5)

Where $\theta$ is the angle of the wing surface and $\epsilon$ is the solar azimuth angle of the wing surface.

By synthesizing the equations, we can calculate $I_{\mathrm{Ds}}$ . In this paper, we consider the solar cell laying on the surface of UAV to be simplified as a horizontal calculation. When flying, the $\theta$ angle is approximately 0. So $\cos i=\sin h$ , therefore:

$$I_{\mathrm{Ds}}=I_0\cdot \exp \left({\displaystyle \frac{-k}{\sin h}}\right)\cdot \sin h$$

(6)

Since the UAV is flying at the equator and the flying time is the time when the sun rises in the east, the altitude of $\sigma =0$ , $\phi =0$ and the sun is only related to the sun time angle, then $\sin h=\cos \omega$ .

$$I_{\mathrm{Ds}}=I_0\cdot \exp \left({\displaystyle \frac{-k}{\sin h}}\right)\cdot \sin h$$

(7)

Calculation of sky scattering radiation

The wing panels of a solar-powered unmanned aerial vehicle (UAV) receive two parts of scattered radiation, namely the refracted and reflected sky scattered radiation energy. ⁶ This part of the Energy II is simplified as follows:

$$I_{\mathrm{Is}}=C\cdot I_0\cdot \exp \left({\displaystyle \frac{-k}{\sin h}}\right)$$

(8)

Where C is the scattering coefficient. If there is a mountain shelter, the solar panels cannot receive the direct radiation of the sun, can only receive scattered radiation from the sky.

Both of these components are affected by $I_0$ , while $I_0$ solar radiation constant is measured regardless of cloud cover and mountain cover. For low-altitude aircraft, the effect of cloud cover on the sun should be considered.

When flying in the equatorial plane, only the solar hour angle and the influence of clouds and mountains should be considered. Therefore, the influence of cloud cover and mountain shadow on the solar radiation energy can be analyzed as part of the main cost. Here’s how.

The shadow is seen as part of the cost of the solar energy conversion that is affected by the Sun shading of the UAV into the shadow area of the mountain. In this paper, we consider the time of the Sun rising in the east at the equator, the Sun illumination and the angle of the ground are h. If the central coordinates of the mountains are ( $x_0,y_0$ ) and the height function $z=f(x,y)$ , the maximum height is $f(x_0,y_0)$ . The coordinates of the nodes are $x_1,y_1,z_1$ , as shown.

The shadow of the terrain is calculated as follows:

when $x_1>x_0$ the aircraft is flying on the sunny side of the mountain, unshaded, $P_{\mathrm{D}}(d_{\mathrm{D}})=0$ .

When $x_1<x_0$ , further analysis is needed:

If $z_1>f(x_0,y_1)$ , then the height is higher than the current height of the mountain, not covered, $P_{\mathrm{D}}(d_{\mathrm{D}})=0$ .

If $z<f(x_0,y_1)$ , when $\mathrm{arctng}((f(x_0,y_1)-z_1)/(x_0-x_1))<h$ is not covered by the mountain, $P_{\mathrm{D}}(d_{\mathrm{D}})=0$ .

When $\mathrm{arctng}((f(x_0,y_1)-z_1)/(x_0-x_1))>h$ , is covered by the mountain shadow. The greater the angle, the longer the masking time may be, at which point the calculation is:

$$P_{\mathrm{D}}(d_{\mathrm{D}})={\displaystyle \frac{\mathrm{arctng}((f(x_0,y_1)-z_1)/(x_0-x_1))}{h}}$$

(9)

When the UAV enters the cloud, it will affect the direct radiation of the sun due to the existence of certain shade. The central part of the cloud is the most effective. As it spreads out, the cloud becomes thinner and more effective. The effect of the cloud on solar UAVs can be approximated as follows:

$$P_{\mathrm{Y}}(d_{\mathrm{Y}})=\{\begin{array}{cc}0& (d_{\mathrm{Y}}>d_{\mathrm{Y}\max })\\ 1/d_{\mathrm{Y}}^2& (d_{\mathrm{Y}}\le d_{\mathrm{Y}\max })\end{array}$$

(10)

In formula (10), $P_{\mathrm{Y}}(d_{\mathrm{Y}})$ indicates the extent of the cloud cover effect, and $d_{\mathrm{Y}}$ indicates the distance between the UAV and the cloud; since only direct solar radiation is affected, even if the total cover conversion energy is 0, there is still an emergency battery available to power the UAV, there is not virtually death threat. That is to say it is not existing for $P_{\mathrm{Y}}(d_{\mathrm{Y}})=1$ . When the UAV flies out of the cloud and is no longer affected by the cloud cover, the solar cell conversion energy can not only meet the flight use of the UAV platform power equipment, but also the spare energy can recharge the storage battery and recharge the power in time. $P_{\mathrm{Y}}(d_{\mathrm{Y}})=0$ means that when the UAV is outside the cloud radius, the impact degree of UAV is 0, then the threat cost of UAV is 0; when $d_{\mathrm{Y}}$ is within the cloud radius $d_{\mathrm{Y}\max }$ , the threat level of solar UAVs can be expressed as $1/d_{\mathrm{Y}}^2$ .

Establishment of alternative cost models

There are some constraints in the threat model, such as the maximum operating distance and the effective killing distance. In this paper, the constraints are introduced into the objective function by establishing a reasonable objective cost function, in this paper, the models of different threats are defined as follows: the detection degree of radar to UAV can be approximately expressed as:

$$P_{\mathrm{R}}(d_{\mathrm{R}})=\{\begin{array}{cc}0& (d_{\mathrm{R}}>d_{\mathrm{R}\max })\\ 1/d_{\mathrm{R}}^4& (d_{\mathrm{R}\min }\le d_{\mathrm{R}}\le d_{\mathrm{R}\max })\\ 1& (d_{\mathrm{R}}<d_{\mathrm{R}\min })\end{array}$$

(11)

In formula (11), $p_{\mathrm{R}}(d_{\mathrm{R}})$ is the radar threat level and $d_{\mathrm{R}}$ is the distance between the UAV and the radar; $d_{\mathrm{R}\max }$ is the maximum radius of the radar detection area. $d_{\mathrm{R}\min }$ represents the effective radar detection radius. The threat cost of $P_{\mathrm{R}}(d_{\mathrm{R}})=1$ is Infinite; the threat cost of $P_{\mathrm{R}}(d_{\mathrm{R}})=0$ is 0. When $d_{\mathrm{M}}$ is in between, the detection probability of a solar-powered UAV can be reduced to.

Terrain threat refers to the threat of impact caused by mountain obstacles when the UAV is flying at low altitude. The cone approximation is shown in Figure 2(a). The terrain radius is $R_{\mathrm{T}}$ , the altitude is H, the slope is $\theta$ , the symmetrical wheelbase between UAV and terrain is $d_{\mathrm{T}}$ , the UAV flying altitude is h ( $h\le H$ ; when $h>H$ , there is no threat).

Figure 1.

Schematic diagram of mountain terrain occlusion cost.

Figure 2.

Schematic diagram of terrain threat cost.

The threat probability of terrain to UAVs varies with distance. The farther away from the mountains, the less likely the impact, and vice versa. $d_{\mathrm{T}\min }$ represents the closest distance allowed by the terrain. The range is too small to avoid, and the impact threat is infinite; $d_{\mathrm{T}\min }$ represents the maximum range of terrain effects, and as shown in Figure 2(b), the terrain threat probability of the UAV can be expressed as follows:

$$p_{\mathrm{T}}(d_{\mathrm{T}},h)=\{\begin{array}{cc}0& ((d_{\mathrm{T}}>R_{\mathrm{T}}+d_{\mathrm{T}\max })(h>H))\\ (d_{\mathrm{T}\max }+R_{\mathrm{T}}(h)-d_{\mathrm{T}})/(d_{\mathrm{T}\max }-d_{\mathrm{T}\min })& (h\le H{R}_{\mathrm{T}}+d_{\mathrm{T}\min }\le d_{\mathrm{T}}\le R_{\mathrm{T}}+d_{\mathrm{T}\max })\\ 1& (h\le H{d}_T<R_{\mathrm{T}}+d_{\mathrm{T}\min })\end{array}$$

(12)

$$\theta =\arcsin \left({\displaystyle \frac{H}{\sqrt{R_{\mathrm{m}}^2+H^2}}}\right)$$

(13)

$$R_{\mathrm{T}}(h)={\displaystyle \frac{(H-h)}{\tan \theta }}$$

(14)

In order to prevent UAVs from exposing their targets, it is necessary for UAVs to make use of the terrain as much as possible to cover their targets:

$$P_{\mathrm{H}}(\xi )=\xi$$

(15)

In formula (15), $P_{\mathrm{H}}(\xi )$ indicates the level of threat the drone is exposed to, and $P_{\mathrm{H}}(\xi )$ indicates the height of the drone; the higher the drone flies without any cover, the greater the likelihood that it will be detected and the greater the threat.

When the optimal objective function is determined, the weight of each edge in the network graph is also determined. That is, the cost weight of the first d edge in the network graph is calculated by the following cost function in the process of searching nodes^5–10:

$$w_{ij}={\delta }_{\mathrm{R}}{w}_{\mathrm{R}ij}+{\delta }_{\mathrm{H}}{w}_{\mathrm{H}ij}+{\delta }_{\mathrm{D}}{w}_{\mathrm{D}ij}+{\delta }_{\mathrm{Y}}{w}_{\mathrm{Y}ij}+{\delta }_{\mathrm{T}}{w}_{\mathrm{T}ij}+{\delta }_{\mathrm{l}}{w}_{\mathrm{l}ij}$$

(16)

In Formula (16), $w_{ij}$ is the price from point i to point j. $w_{\mathrm{R}ij}$ , $w_{\mathrm{H}ij}$ , $w_{\mathrm{D}ij}$ , $w_{\mathrm{Y}ij}$ , $w_{\mathrm{T}ij}$ and $w_{\mathrm{l}ij}$ represent other corresponding costs for the edge, respectively.

Therefore, the fundamental problem to be solved in solar route planning is to obtain a route with the lowest comprehensive cost, that is, the cumulative cost of each segment of the route is the smallest, that is, the solution is the least cost $\min \sum w_{ij}$ . Which makes up the waypoints of the route.

Since the essence of this problem is to solve the combinatorial optimization problem with the best value under specific constraints, the algorithm can be constructed and solved by intelligent optimization algorithms. I carried out the research on the route planning of ant colony algorithm for ordinary UAVs in the early stage, and the algorithm has good convergence effect and is relatively easy to engineer. Therefore, the algorithm will continue to be used for the route planning problem of solar UAV, the cost model of the combinatorial optimization problem established in this section, and the specific construction and implementation of the ant colony-based algorithm will be detailed in the next section, including exploration factor design, state transition strategy, and pheromone update strategy.

Ant colony algorithm for route planning

Ant colony algorithm flow for route planning is shown in Figure 3.

Figure 3.

Algorithm flowchart.

Exploratory factor design

The spatial coordinates representing the start and destination of the route, respectively. Represents the granularity of planning spatial meshing. Let $N_x$ , $N_y$ , and $N_z$ have the number of meshes in the x, y, and z axes, respectively, then:

$$N_x={\displaystyle \frac{(x_{\mathrm{end}}-x_{\mathrm{st}})}{x_{\mathrm{Gird}}}}$$

$$N_y={\displaystyle \frac{(y_{\mathrm{end}}-y_{\mathrm{st}})}{y_{\mathrm{Gird}}}}$$

$$N_z={\displaystyle \frac{(z_{\mathrm{end}}-z_{\mathrm{st}})}{z_{\mathrm{Gird}}}}$$

Then the mesh number $N_{\mathrm{g}}$ of the entire planning space is calculated as follows:

$$N_{\mathrm{g}}=N_x\times N_y\times N_z$$

(17)

The pheromone matrix $\tau$ represents the pheromones between different nodes, so $\tau$ is ( $N_{\mathrm{g}}\times N_{\mathrm{g}}$ ) matrix, then the planning space is a three-dimensional space starting and ending point of the UAV, and the number of storage space grids occupied by it is:

$$\begin{array}{c}{M}_{\mathrm{em}}=N_{\mathrm{g}}^2\\ ={(N_x\times N_y\times N_z)}^2\\ ={\left({\displaystyle \frac{x_{\mathrm{end}}-x_{\mathrm{st}}}{x_{\mathrm{Gird}}}\times {\displaystyle \frac{y_{\mathrm{end}}-y_{\mathrm{st}}}{y_{\mathrm{Gird}}}\times {\displaystyle \frac{z_{\mathrm{end}}-z_{\mathrm{st}}}{z_{\mathrm{Gird}}}}}}\right)}^2\end{array}$$

(18)

Exploratory factor design

The discovery factor was designed as a countdown to the threat cost. In the case of node j, the threat cost can be expressed as the sum of the costs of the six parts:

$${\epsilon }_j={\delta }_{\mathrm{R}}{w}_{\mathrm{R}j}+{\delta }_{\mathrm{H}}{w}_{\mathrm{H}j}+{\delta }_{\mathrm{D}}{w}_{\mathrm{D}j}+{\delta }_{\mathrm{Y}}{w}_{\mathrm{Y}j}+{\delta }_{\mathrm{T}}{w}_{\mathrm{T}j}+{\delta }_{\mathrm{l}}{w}_{\mathrm{l}j}$$

(19)

In Formula (19), $w_{\mathrm{R}j}$ , $w_{\mathrm{H}j}$ , $w_{\mathrm{D}j}$ , $w_{\mathrm{Y}j}$ , $w_{\mathrm{T}j}$ and $w_{\mathrm{l}j}$ are the costs of radar, altitude, ground shadow, cloud shadow, terrain and energy consumption, respectively.

According to the formula above, the exploring factor of node j is calculated, which is the reciprocal of the threat cost of the node, so ${\eta }_j=1/{\epsilon }_j$ , the smaller the node j cost is, the larger the exploring factor is and the higher the resolution is, whereas, the exploring factor is small, the differentiation is not big and the resolution is low. The exploratory factor can accelerate the convergence of the algorithm and improve the resolution of nodes, so it can be used as an important supplement of pheromone.

State transition strategies

Let ant scale be $gm$ and ${\tau }_{ij}(N)$ denote the pheromone concentration of route from node i to node j in n iterations. At the beginning of iteration, the pheromone concentration of each route is the same, and ${\tau }_{ij}(0)=\nu$ ( $\nu$ is constant). According to the pheromone concentration in each path and the exploring factor, the ant k ( $k=1,2,\dots ,gm$ ) decides the direction of migration. $P_{ij}^k(N)$ indicates the probability of the N generation ant k transferring from node i to node j.

$$P_{ij}^k(N)=\{\begin{array}{cc}{\displaystyle \frac{{\tau }_{ij}{(N)}^{\alpha }{\eta }_{ij}{(N)}^{\beta }}{\sum _{s\in \mathrm{Allowe}{\mathrm{d}}_k}{\tau }_{is}{(N)}^{\alpha }{\eta }_{is}{(N)}^{\beta }},}& j,s\in \mathrm{Allowe}{\mathrm{d}}_k\\ 0,& \mathrm{otherwise}\end{array}$$

(20)

In Formula (20), $\mathrm{Allowe}{\mathrm{d}}_k=\mathrm{Neigh}{\mathrm{b}}_i\cap \overline{\mathrm{Tab}{\mathrm{u}}_k}$ and $\mathrm{Tab}{\mathrm{u}}_k$ are used to record the set of nodes that ant k has experienced in contemporary computing. $\mathrm{Tab}{\mathrm{u}}_k$ increases dynamically as the ant selects the next node. $\overline{\mathrm{Tab}{\mathrm{u}}_k}$ represents nodes other than $\mathrm{Tab}{\mathrm{u}}_k$ , $\mathrm{Neigh}{\mathrm{b}}_i$ represents the current node I adjacent node set, $\mathrm{Allowe}{\mathrm{d}}_k$ represents ants, and the next step allows the set of nodes to be selected. $\mathrm{Allowe}{\mathrm{d}}_k$ is not directly equal to $\overline{\mathrm{Tab}{\mathrm{u}}_k}$ , but can only be selected from 26 adjacent nodes, and exclude calendar nodes. $\alpha$ and $\beta$ denote the importance of pheromones and exploration factors.

Pheromone update strategy

The parameter $(1-\rho )$ indicates the degree of volatilization of the PHEROMONE. Every time N increases, the pheromone volatilizes. The pheromone concentrations of all routes experienced in the iteration cycle need to be updated and calculated as follows:

$$\{\begin{array}{c}{\tau }_{ij}(N+1)=\rho *{\tau }_{ij}(N)+(1-\rho )*\sum _{k=1}^{gm}\mathrm{\Delta }{\tau }_{ij}{(N)}^k,\rho \in (0,1)\\ \mathrm{\Delta }{\tau }_{ij}=\sum _{k=1}^{gm}\mathrm{\Delta }{\tau }_{ij}{(N)}^k\end{array}$$

(21)

$\mathrm{\Delta }{\tau }_{ij}^{}{(N)}^k$ represents pheromone enhancement for all routes that are traversed by the N sub-optimization, where k represents the kth ant.

$$\mathrm{\Delta }{\tau }_{ij}{(N)}^k=\{\begin{array}{cc}Q/W_N^k,& \mathrm{if}:ij\in \mathrm{tab}{\mathrm{u}}_k(N)\\ 0,& \mathrm{otherwise}\end{array}$$

(22)

Among them, the Q table pheromone increases intensity coefficient, $W_N^k$ represents the smallest cost in the N th iteration, and its corresponding route is the optimal route. Q, $\rho$ determines its value according to the scale of solution. When the number of iterations reached a predetermined maximum number of iterations, the algorithm stopped. In addition, the best route of the previous generation is retained as a route in the contemporary route to achieve elite strategy.

Since the standard ACO cannot find the target node within the given number of ant colonies and the number of iterations when dealing with 3D route planning, it is necessary to improve the algorithm and introduce the target node factor to guide the ant colony to search for the target node.

TNF-ACO algorithm solves route planning problems

TNF-ACO algorithm refers to the construction of target node factors on the basis of the standard ACO algorithm, so that the algorithm.

Target node factor design

To solve the above two problems, this topic introduces the target node factor and sets the maximum number of track nodes. The conventional ant colony algorithm state transfer strategy is selected according to the probabilistic selection of pheromones and heuristic factors, without the target node factor, and the heuristic factor values of the 26 adjacent nodes in the state transfer strategy are not much different, so the local predictability of the target node is not strong, especially in the early stage of algorithm iteration, the pheromones between nodes are not much different, and the state transfer is easy to fall into blind selection, and it is difficult to quickly reach the target node. According to the needs of track planning problems, this paper introduces the target node as part of the state transition strategy, taking node i as an example, and the target node factor Q is expressed as:

$${\lambda }_i={\displaystyle \frac{1}{d_{iD}}}$$

(23)

In Equation (23), $d_{iD}$ represents the distance from node i to the target node $L_{max}$ . It can be seen from the formula ${\lambda }_i=1/d_{iD}$ that the farther away from the target node, the larger the target node factor; conversely, the smaller. Therefore, in the state transfer strategy, the introduction of target node factor can effectively reduce the blindness of ant search, so that ants search for tracks in the direction of target nodes.

In addition, the UAV has a maximum range parameter, that is, the flight range of the UAV throughout the flight, which is limited by the ration of aircraft fuel and flight time. If the maximum track length is $L_{max}$ , the distance $l_i$ of each flight segment should satisfy:

$$\sum _i{l}_i\le L_{max}$$

(24)

According to Equation (22), the number of track nodes cannot exceed a certain range, which is limited by the UAV’s own conditions. In addition, it can be seen from Equation (1) that the track cost includes the cost of fuel consumption, so the track with high fuel consumption may lead to an increase in the cost of the track. Therefore, the cost of fuel consumption must also meet a certain range, beyond which it is considered unacceptable. Therefore, although the track node is not fixed, the maximum number of nodes can be set, and when it is exceeded, the track is considered unfeasible and abandoned.

State transition policy improvements

$P_{ij}^k(N)$ represents the probability of the N generation ant k being transferred from node i to node j, and its calculation formula is shown in (3.9).

$$P_{ij}^k(N)=\{\begin{array}{cc}{\displaystyle \frac{{\tau }_{ij}{(N)}^{\alpha }{\eta }_{ij}{(N)}^{\beta }{\lambda }_{ij}{(N)}^{\gamma }}{\sum _{s\in \mathrm{Allowe}{\mathrm{d}}_k}{\tau }_{is}{(N)}^{\alpha }{\eta }_{is}{(N)}^{\beta }{\lambda }_{is}{(N)}^{\gamma }},}& j,s\in \mathrm{Allowe}{\mathrm{d}}_k\\ 0,& \mathrm{otherwise}\end{array}$$

(25)

In Equation (25), ${\lambda }_i(N)$ represents the Nth generation target node factor of node i, and the other parameters have the same meaning as in Equation (21). To ensure that the ants eventually reach the target node, the ants need to search for tracks in the direction of the target node.

Different from the standard ACO algorithm, this paper introduces the target node factor ${\lambda }_i(N)$ into the state transition probability. In addition, $\gamma$ represents the degree of importance of the target node factor and is a constant.

The values of $\beta$ and $\gamma$ should not be too large, otherwise it is easy to make the algorithm fall into local optimum, resulting in stagnation.

Considering the constraints such as different threat ranges, this paper also considers the conditional constraints for the heuristic factor ${\eta }_j$ , that is, as long as the distance from a threat point to a selected node j is less than the effective distance of the threat type, the threat cost of node j can be determined to be infinite without considering other threat points.

Then ${\eta }_j=0$ , then the probability of selecting an infinitesimal node is 0, so as to exclude the nodes that do not meet the constraints in $\mathrm{Allowe}{\mathrm{d}}_k$ , and ensure that the algorithm will not appear through the effective killing distance of the threat point.

The ES-TNF-ACO algorithm solves the route planning problem

ES-TNF-ACO algorithm

The ES-TTF-ACO algorithm refers to an algorithm that introduces an elite strategy on the basis of the TNF-ACO algorithm and retains the excellent population of the previous generation as a route for the next generation. Specifically, the following three improvements have been made:

Retain the best routes of the previous generation as one of the contemporary routes and achieve an elite strategy. Achieve relatively fast convergence and reduce iteration of invalid paths.

Discard pheromone enhancement on invalid routes where all ant colonies do not reach the target point.

Among the 26 adjacent nodes of the current node, exclude nodes whose search direction does not deviate from the target node, so that the state transition probability of the current node transferring to this node is 0.

Specific implementation methods

• (1) The path node of the first ant is fixed to the path node of the previous generation with the most ants, so as to retain the elite and realize the convergence and search optimization of the algorithm.

The specific implementation code is as follows:

if CN>=2

Ant_x(1,:)=P_Best_Ant_x(CN-1,:) ;

Ant_y(1,:)=P_Best_Ant_y(CN-1,:) ;

Ant_z(1,:)=P_Best_Ant_z(CN-1,:) ;

Ant_lie(1,:)= P_Best_Ant_lie(CN-1,:) ;

Ant_ha(1,:) =P_Best_Ant_ha(CN-1,:) ;

Ant_gao(1,:) =P_Best_Ant_gao(CN-1,:) ;

Ant_NetPoint(1,:)=P_Best_Ant_NetPoint(CN-1,:) ;

end

where, CN represents the current number of iterations, Ant_x(1,:) represents the coordinate x of the first ant, Ant_y(1,:) represents the coordinate y of the first ant, Ant_z (1,:) represents the coordinate z of the first ant, Ant_lie (1,:) represents the column sequence number in the spatial grid of the first ant, Ant_ha(1,:) represents the row serial number of the first ant in the planning spatial grid, and Ant_gao(1,:) Indicates the high ordinal number of the first ant in the planning space grid, and the Ant_NetPoint (1,:) indicates the node ordinal number of the first ant in the planning space grid.

• (2) Pheromone update strategy improvements

If there is an ant k route path length equal to the maximum path length $L_{max}$ allowed for colony search in the (n-1)th generation, the pheromones of all paths of that path are not enhanced, i.e. $\mathrm{\Delta }{\tau }_{ij}^{}{(N)}^k=0$ .

The specific implementation code is as follows:

Delta_Tau = zeros(n,n);

for i = 1:gm

if leght(i)==Lmax

else

for j = 1:(leght(i)-1)

Delta_Tau(Ant_NetPoint(i,j),Ant_NetPoint(i,j + 1))=Delta_Tau(Ant_NetPoint(i,j),Ant_NetPoint (i,j + 1))+Q/Ant_F(i);

end

end

end

Among them: Delta_Tau represents the pheromone increase matrix, gm represents the number of ants, and letht calculates the route length of the ith ant.

(3) Set the direction divergence node shift probability to 0

The specific implementation code is as follows:

if (Xend-Xst)/(Xo-Xi) < 0 || (Yend-Yst)/(Yo-Yi)|| (Zend-Zst)/(Zo-Zi)

P(i) = 0;

else

(note: Shift probability is calculated according to the normal calculation)

End

Example simulation of route planning

In this paper, the route planning of solar energy UAV is simulated by intelligent optimization algorithm. Simulation set the specific radar, clouds, flight starting and ending coordinates, as shown in Table 1. The parameters of radar model and weight are set to $d_{\mathrm{R}\min }=4$ ; $d_{\mathrm{R}\max }=120$ ; $d_{\mathrm{Y}\max }=50$ ; $d_{\mathrm{T}\min }=2$ ; $d_{\mathrm{T}\max }=10$ ; $h=5$ , respectively. Four radar threat points and eight cloud shadows are set up, and the location coordinates are shown in Table 1. In addition, the coordinates of the starting point and target location point of solar navigation planning are (5, 90, 1.4) and (95, 20, 2.8), respectively.

Table 1.

Coordinates of some threat points.

Threat type	Serial number	Coordinates
Radar	1#	61, 72, 2.11
	2#	73, 12, 1.41
	3#	80, 65, 1.72
	4#	22, 33, 1.12
Cloud mass	1#	65, 10, 5.10
	2#	30, 65, 4.01
	3#	72, 57, 6.11
	4#	50, 40, 7.02
	5#	50, 90, 2.94
	6#	25, 72, 0.12
	7#	44, 20, 0.65
	8#	70, 45, 0.14
Target point (95,20,2.8)
Starting point (5,90,1.4)

The terrain model is composed of six mountains with different centers and heights. The figure of the natural index of the mountains with e as the base is expressed as follows.

$$z_i(x,y)=h_i\times e^{-({(x-a_i)}^2/10)-({(y-b_i)}^2/10)}$$

(26)

$$z(x,y)=\sum _{i=1}^{\mu }{z}_i(x,y)=\sum _{i=1}^{\mu }(h_i\times e^{{}^{+({(x-a_i)}^2/10)-({(y-b_i)}^2/10)}})$$

(27)

The mountain parameters are: a = [20 40 50 70 65 85], b = [20 65 40 55 10 85], c = [3 4 7 10 5 7]. In formula (21) and formula (22): X, Y are the coordinates in the digital map (x,y), the planning range is [0, 100]. ( $a,b$ ) indicates the coordinates of the central axis of the mountain, with the c being the highest height of the mountain.

Where, ${\delta }_{\mathrm{l}}$ = 0.0015; ${\delta }_{\mathrm{R}}$ = 0.768; ${\delta }_{\mathrm{H}}$ = 0.0035; ${\delta }_{\mathrm{D}}$ = 0.002; ${\delta }_{\mathrm{Y}}$ = 0.2; ${\delta }_{\mathrm{T}}$ = 0.025. MCN = 200; $gm$ = 35; $\alpha$ = 0.5; $\beta$ =0.25; $\rho$ = 0.1; Q = 10 (&5,1); $\mathrm{Qif}{\mathrm{a}}_{\mathrm{Limit}}$ = 4; $X_{\mathrm{Gird}}$ = 5; $Y_{\mathrm{Gird}}$ = 5; $Z_{\mathrm{Gird}}$ = 0.7; $L_{max}$ = 105.

Figure 4(a) and (b) are 3D plots and contour plots of 3D track planning of standard ACO algorithm, respectively. Figure 5(a), (b) and (c), respectively, simulate the 3D trajectory planning 3D map, contour plot and convergence curve obtained by TNF-ACO algorithm under the condition of Q = 10. Figure 6(a), (b) and (c), respectively, simulate the 3D trajectory planning 3D map, contour map and convergence curve obtained by TNF-ACO algorithm under Q = 5. Figure 7(a), (b) and (c), respectively, simulate the 3D trajectory planning 3D map, contour plot and convergence curve obtained by TNF-ACO algorithm under Q = 1 (Figure 8(a)–(c)). The 3D trajectory planning 3D diagram, contour plot and convergence curve obtained by the ES-TNF-ACO algorithm under the condition of Q = 10 are, respectively, Figure 9(a), (b) and (c), respectively, simulate the 3D trajectory planning 3D plot, contour plot and convergence curve obtained by the ES-TNF-ACO algorithm under the condition of Q = 5. Figure 10(a), (b) and (c), respectively, simulate the 3D trajectory planning 3D map, contour plot and convergence curve obtained by the ES-TNF-ACO algorithm under the condition of Q = 1. In the figure, the square on the left is the starting point of the route, the circle on the right is the end point, the pentagram connected between the two is the route node, and the other space nodes represent the coordinates of the threat radar and the center of the cloud.

Figure 4.

3D trajectory planning, contour plot and convergence curve of standard ACO algorithm.

Figure 5.

3D track planning, contour plot and convergence curve of TNF-ACO algorithm (optimal value 0.3189, Q = 10).

Figure 6.

3D track planning, contour plot and convergence curve of TNF-ACO algorithm (optimal value 0.3299, Q = 5).

Figure 7.

3D trajectory planning, contour plot and convergence curve of TNF-ACO algorithm (optimal value 0.3312, Q = 1).

Figure 8.

ES-TNF-ACO algorithm 3D track planning, contour plot and convergence curve (optimal value Q = 10,0.3130).

Figure 9.

ES-TNF-ACO algorithm 3D track planning, contour plot and convergence curve (optimal value Q = 5,0.3111).

Figure 10.

ES-TNF-ACO algorithm 3D track planning 3D diagram, contour plot and convergence curve (optimal value 0.3061, Q = 1).

It can be seen from Figure 4 that due to the large search range of three-dimensional space, if the target node factor is not introduced and the standard ant colony algorithm is directly adopted, the success rate of the algorithm reaching the target node is very low, that is, the predetermined task of track planning cannot be completed, so the target node factor must be introduced.

It can be seen from Figures 4, 5 and 6 that after the introduction of the target node factor, compared with the standard ant colony algorithm, the success rate of reaching the target node is very low, and the algorithm can reach the target node smoothly, but there is still the problem of the algorithm not convergence, and the randomness is strong.

It can be seen from Figures 8, 9 and 10 that on the basis of the TNF-ACO algorithm, the elite strategy is introduced, and the ES-TNF-ACO algorithm can converge every time and obtain better optimization results under the same Q parameters.

Table 2 shows the statistical analysis of the TNF-ACO algorithm and the improved algorithm ES-TNF-ACO with elite subgroups under the condition that different Q parameters are set. Q indicates the pheromone increase strength factor. Since the standard ACO algorithm does not reach the track target node, its statistics are of little practical significance. In the TNF-ACO algorithm and the improved ES-TNF-ACO algorithm, other parameter configurations are the same. According to Table 2, in Q = 10, Q = 5 and Q = 1, the ES-TNF-ACO algorithm has obvious advantages over the TNF-ACO algorithm, with 100 statistical data minimum, maximum, average, and variance, and 200 iterations in each loop, the former of which has better optimization effect than the latter. The results of 200 iterations of the optimal loop out of 100 cycles are shown in Figures 5, 6, 7, 8, 9 and 10.

Table 2.

3D spatial simulation data description statistics comparison of improved ACO algorithm.

Parameter	Statistical values	Track cost statistics comparison
		TNF-ACO	ES-TNF-ACO
$Q=10$	Min	0.3189	0.3130
	Max	0.4070	0.3999
	Average	0.3779	0.3565
	Variance	2.9706e-004	3.5888e-004
$Q=5$
	Min	0.3299	0.3111
	Max	0.4189	0.3992
	Average	0.3778	0.3556
	Variance	3.4271e-004	2.9686e-004
$Q=1$
	Min	0.3312	0.3061
	Max	0.4063	0.3988
	Average	0.3762	0.3547
	Variance	2.0041e-004	3.1147e-004

It can be seen from Table 2 that the ES-TNF-ACO algorithm obtains better results than other parameters when Q = 1, and its optimization route and convergence curve are shown in Figure 10(c), with the continuous adjustment of iterations, the route cost gradually converges, and after 157 iterations, the cost converges to the minimum cost of the global route 0.3061, indicating that the algorithm is indeed effective and feasible. The radar, altitude, terrain shadow, cloud shading, terrain threat and fuel consumption costs of the route with the lowest cost were 0.0040, 46.9000, 11.1867, 0.2154, 4.7539 and 26, respectively, and the weighted costs were 0.0027, 0.1544, 0.0148, 0.0365, 0.0557 and 0.0420, respectively, for a total cost of 0.3061.

Conclusion

In order to solve the specific route planning problem of solar UAV, this paper proposes the TNF-ACO algorithm, and further improves the ES-TNF-ACO algorithm. In this paper, solar UAV flight at the equator position, at a certain point in the process of rising from the east, and simplifying the laying of solar cell wings to a flat installation, without considering the influence of the flight attitude of the aircraft itself, by establishing the model of mountain impact cost, height threat, mountain shadow shading cost, cloud shading cost, based on ant colony intelligent optimization algorithm, the ant colony path is reasonably constructed, and the maximum action distance of various costs. The constraints such as solar irradiation angle and effective action distance are introduced into the cost model and heuristic factor calculation, which solves the comprehensive optimization problem of solar UAV route planning. The simulation results show that the algorithm path structure is reasonable, the target node can be found independently, and the convergence speed can meet the requirements of route planning, and the generated route cost is small, which can be used for static planning design.

However, the algorithm only considers the laying of solar cells on the upper surface of the drone, and does not consider the laying of the lower surface. And when laying, it is simplified to horizontal plane calculation, and the level flight attitude is maintained when flying, and the situation of adjusting the flight attitude of UAV maneuvering is not taken into account. Therefore, the model still has certain limitations, and we will conduct in-depth research from the above mentioned areas later.

Finally, the potential application areas of the algorithm can be extended to solar airships, overpressure balloons and other fields.