Routing in Solar-Powered UAV Delivery System

Abstract

As interest grows in unmanned aerial vehicle (UAV) systems, UAVs have been proposed to take on increasingly more tasks that were previously assigned to humans. One such task is the delivery of goods within urban cities using UAVs, which would otherwise be delivered by terrestrial means. However, the limited endurance of UAVs due to limited onboard energy storage makes it challenging to practically employ UAV technology for deliveries across long routes. Furthermore, the relatively high cost of building UAV charging stations prevents the dense deployment of charging facilities. Solar-powered UAVs can ease this problem, as they do not require charging stations and can harvest solar power in the daytime. This paper introduces a solar-powered UAV goods delivery system to plan delivery missions with solar-powered UAVs (SPUs). In this study, when the SPUs run out of power, they charge themselves on landing places provided by customers instead of charging stations. Some advanced path planning algorithms are proposed to minimize the overall mission time in the statically charging efficiency environment. We further consider routing in the dynamically charging efficiency environment and propose some mission arrangement protocols to manage different missions in the system. The simulation results demonstrate that the algorithms proposed in our work perform significantly better than existing UAV path planning algorithms in solar-powered UAV systems.

1. Introduction

As unmanned aerial vehicles (UAVs) become increasingly popular, their applications for drones/UAVs (drones and UAVs can be used interchangeably, although many professionals in the industry believe UAVs need to have autonomous flight capabilities, whereas drones do not; in this study, we only talk about UAVs in the form of quadcopters) are developed and expanded to many civilian or military areas. Soon, UAVs are predicted to play a significant role in the intellectualization of urban cities. However, the limited battery capacity makes it difficult for UAVs to finish long-distance/long-term missions without any power charging on the way.

In some UAV delivery systems, the solution is to build many charging stations along the way, so that UAVs can carry out longer missions. For example, in [1,2,3], charging stations were built for UAV power supply. Unfortunately, it is costly and sometimes illegal to build charging stations on rooftops or personal property in urban cities. Although some hosts are willing to have charging stations on their property, other issues arise, such as the power, the safety of the station and potential UAV-caused accidents. It is, therefore, challenging to build UAV charging infrastructures in both public and personal places.

Solar-powered UAVs can harvest sunlight power anywhere in the daytime given suitable weather. As a result, to achieve more extended distance/time missions by UAVs without building many charging infrastructures, we propose an autonomous solar-powered UAV delivery system to take over the traditional delivery work in urban cities. In our proposed system, places such as roofs of buildings (the study in [4] presented potential solar energy stations on rooftops in modern cities that can be easily discovered by UAVs), high grounds in parks and small places in people’s gardens are identified as possible charging places on city maps. It is much cheaper and more efficient to simply utilize such places as they are rather than building charging infrastructures there instead. In the proposed system, hundreds of solar-powered UAVs are prepared in stores with thousands of charging places determined throughout a city. When the system receives a delivery mission, it automatically finds the proper UAV that takes the least amount of time to finish the task. Since the UAVs can run out of power during the mission, the system helps with UAV path planning and determines how much time they need to stay on each passed landing place for solar power charging. In previous studies, authors paid more attention to path planning for solar-powered UAVs to harvest more energy during the flight (for fixed-wing UAVs), and no path planning methods were proposed for solar-powered UAV delivery systems. Thus, we believe it is necessary to comprehensively analyze the routing problem in solar-powered UAV delivery systems.

This study discusses the path planning problem by considering both the distance and the energy charging efficiency, since SPUs are applied here. The charging rates (efficiency) (the charging rate here represents the charging speed of a solar panel in a landing place in units of energy per minute, depending on the illumination/radiation level of the landing place) of different charging places differ in our scheme. Sequentially, we perform a trade-off between the distance and energy charging efficiency in this problem. The charging efficiency can either be static or dynamic in UAV systems. If the charging efficiency is considered static, it means the energy price in the normal UAV delivery system or the charging rate in the solar-powered UAV system is considered constant for each landing place. At the beginning of our study, both the optimal charging time assignment (CTA) algorithm and the globally optimal path planning algorithm (GOA) with pruning strategies are proposed to solve the path planning problem in the statically charging efficiency environment. As described in the solar radiation model proposed in [5], we also consider a dynamically charging efficiency in our study. It is much more efficient to charge the UAV in the period close to noon than in earlier or later daytime. Hence, the path planning problem became more complicated here. To solve the routing problem in the dynamically charging efficiency environment, we also propose a dynamically greedy CTA algorithm to further reduce the charging time on a known path. In addition, a heuristic path planning algorithm based on the dynamically greedy CTA algorithm is proposed to reduce the path finding time. Finally, some protocols are proposed to optimize the global UAV delivery works in the whole system.

The main contributions in this paper include the following:

• We propose globally optimal CTA and path planning algorithms for UAV delivery problems in a statically charging efficiency environment. Previous studies have not proposed such globally optimal algorithms.
• We propose a greedy CTA algorithm and a heuristic path planning algorithm for the solar-powered UAV delivery problem in a dynamically charging efficiency environment, which has never been proposed before.
• We are the first to propose a SPU delivery system with a prototype, system design, and protocols to operate different missions.

The rest of this paper is organized as follows: Section 2 compares our study with related works. Section 3 proposes the path planning problems we encounter in the SPU delivery system. Section 4 proposes the CTA algorithm and globally optimal path planning algorithm for a single SPU mission in a statically charging efficiency environment. Section 5 proposes a dynamically greedy CTA algorithm and heuristic SPU path planning algorithm to solve the SPU routing problem in the dynamically charging efficiency environment. Section 6 proposes the mission arrangement protocols to solve the system-level problems. Section 7 proves the feasibility of applying SPUs in delivery works. Section 8 presents the simulation results for our proposed algorithms. The last section concludes the paper and suggests future work.

2. Related Studies

The subject of this paper was to solve the path planning problem in solar-powered UAV delivery systems, whereas no other studies explored the same topic as our study. Most papers about solar-powered UAV path planning have focused on fixed-wing solar-powered UAV path planning, which plans the optimal path to create wings with solar panels that can harvest the most power during the flight, such as [6,7]. However, since we aimed to work on the UAV delivery problems, it was not practical to apply fixed-wing solar-powered UAVs in this study. (Fixed-wing UAVs require runways to take off and land on, and it is not practical to build a lot of runways in modern cities. In addition, in our study, the optimization of charging time on a series of landing places was considered, whereas in fixed-wing solar-powered UAV path planning problems, the optimization of continuous charging efficiency on different flight paths is generally considered. Thus, we did not compare our proposed algorithms with any algorithms applied to fixed-wing UAVs.) As a result, we focused on the path planning problem for solar-powered quadcopters, which are more practical for delivery missions.

Some other studies discussed that the path planning problems of UAVs were related to obstacle avoidance, such as [8,9]. In our study, we ignored this kind of problem because we assumed that UAVs would fly to a safe altitude where no obstacles needed to be considered. Consequently, we focused on finding the best route for UAVs to finish a delivery mission with several landing times. Thus, we did not compare our algorithms with such kinds of studies.

Hence, we compared our study with papers related to quadcopter path planning problems and UAV delivery problems. As described in [10], since UAV path planning problems are always treated as being NP-hard (nondeterministic polynomial-time hardness) problems, the time complexity and complexity of the problems increase in complex terrains. They pointed out all the challenges for solar-powered UAV implementations, but did not give a particular solution to the UAV delivery path planning problem. In addition, different from fixed-wing UAVs, solar-powered quadcopter path planning problems are mainly affected by distance and solar power factors. Thus, in our study, we focused on the trade-off between energy charging and flying distance to find the best route. The authors of [11] concluded that in UAV routing problems, which require the finding of a set of locations as paths, authors are always trying to find a path that minimizes the operational cost. Accordingly, in our study, this cost stood in for the mission time, since no money was deducted during the solar-powered UAV delivery missions.

For UAV delivery research, such as [12,13,14], authors proposed routing algorithms for UAV delivery problems with the help of ground vehicles. Solar-powered UAVs are not appropriate for these delivery systems, because it takes more area to charge them than normal UAVs and requires a safe charging environment, while, in [15,16], the authors studied the path planning algorithm for the coverage plan of delivery customers. The UAV delivery system was discussed in some studies before, such as [17,18,19]. For example, in [17], the authors focused on the optimization of mission scheduling, while in our study, we focused on the optimization of path planning. Thus, the algorithms proposed in such kinds of studies cannot be compared with our study as well.

Studies from [1,2,3,20,21] were quite like our study. All these studies talked about path planning algorithms in UAV delivery systems with charging stations deployed. However, authors of [1,3,20] proposed heuristic algorithms to solve the path planning problems in UAV delivery systems, whereas authors of [2,10] proposed Dijkstra-based algorithms to solve the path planning problems. Since we are the first to propose path planning algorithms in solar-powered UAV delivery systems, it is quite difficult to compare our proposed algorithms with existing heuristic algorithms. The existing heuristic algorithms for UAV delivery path planning did not consider the factor of dynamically charging efficiency, which is very important in solar-powered UAV delivery routing problems. Additionally, in our study, we completed a comprehensive study of the optimization of solar-powered UAV routing problems, which has never been performed before. However, it is difficult to compare our proposed algorithms with existing heuristic algorithms, since the weights of factors were different.

On the other hand, when compared to Dijkstra-based algorithms, our proposed algorithm was much better, because of the introduction of charging time assignment algorithms. In Dijkstra-based algorithms, authors assume that UAVs always charge themselves to a full battery state without considering CTA plans (in Dijkstra-based algorithms, the cost between nodes should be constant, and if UAVs are allowed to be partially charged in a charging station/landing place, the cost between nodes becomes a variable when the paths change). As a result, the Dijkstra-based algorithms cannot obtain globally optimal solutions in any scenarios where UAVs can be partially charged. More specifically, we proposed different CTA algorithms that modified the linear greedy algorithm proposed in [21] for vehicle recharging in the UAV delivery system. In a statically charging efficiency environment, our proposed CTA algorithm had an advantage compared with the linear greedy algorithm in SPU problems, because we could prune the path if any node was assigned 0 charging time by the linear greedy algorithm and reduce the total mission time in some cases. In a dynamically charging efficiency environment, our proposed dynamically greedy CTA algorithm further optimized the linear greedy algorithm with a local minimization operation and increased its efficiency when dealing with SPU path planning problems. With the great advantage brought by our proposed CTA algorithms for solar-powered UAV delivery systems, our proposed path planning algorithms showed obvious improvements compared to Dijkstra-based algorithms.

3. The Proposed Model of SPU Delivery System

3.1. SPU Delivery System

To replace the human work in delivery systems in urban cities, we proposed a system with solar-powered UAVs, which allowed for energy charging on landing places during their missions. There are two kinds of solar-powered UAVs, fixed-wing SPUs and multirotor SPUs. We chose multirotor solar-powered UAVs, because it is difficult for fixed-wing SPUs to lift and land in urban cities. In addition, it can be inefficient, and unsafe for multirotor SPUs to fly with their wings expanded during delivery missions when they are expected to carry heavy payloads. Thus, we preferred the transformable multirotor SPUs from [22], because, during the flight, the solar panels can be folded and unfolded again when landing to charge. The feasibility study of such an SPU is discussed in Section 7. Each UAV could only carry one package at a time, since the payload of a single UAV is limited.

Therefore, the solar-powered UAV delivery system could be built with three elements: delivery stores, solar-powered UAVs and landing places. The packages were stored at delivery stores, and these were the start nodes of the delivery missions. The solar-powered UAVs that carried out the delivery missions charged themselves on landing places. The landing places could either be intermediate nodes or the destination of the delivery missions, and in some emergency or the return voyages, they could be the start nodes as well.

In the proposed system, we assumed that the system would operate during the daytime, with optimal weather conditions for flight and solar power charging. As the energy required for flying was larger than the capacity of the battery, the UAVs periodically needed to land and stay grounded until having accumulated enough energy to continue their mission (this is reminiscent of the way that flying bugs travel, such as flies or bees, where these bugs periodically land to rejuvenate); we referred to such times as charging times. The charging rates at different locations were typically different and were assumed to be known (or approximately predictable) to the path planning facility. The goal of the proposed algorithms was to minimize the sum of total flying times and charging times, (see Equation (1) below). Thus, the proposed algorithm determined the path the UAV should follow, and how much time the solar-powered UAV should stay at each passed landing place. Note that for any UAV delivery mission, the path with the shortest distance may not be the path with the shortest total mission time. In other words, we determined a trade-off between flying distance and energy charging efficiency. In practice, the solar-powered UAV path planning problem was more complicated, since our path planning algorithm was proposed for solar-powered UAVs and the charging rate in any landing place varied in the daytime; we also needed to find the best path planning algorithm in a dynamically charging efficiency environment. Finally, after finding the best path planning algorithm for each SPU delivery mission, the arrangement protocol for hundreds of missions was needed.

3.2. SPU Routing Problem Definition (Statically Charging Efficiency)

Before we optimized the whole SPU delivery system, we needed to solve the SPU routing problem for a single UAV. The ultimate objective of this path planning problem, called the single SPU problem, was to find a path that minimized the total mission time. The following diagram shows an example of a single SPU problem in our proposed system.

In Figure 1, the single SPU problem was designed to find a path from node 1 to node 10 with minimized mission time. Path 1-3-8-10 was a candidate optimal solution to this problem. The UAV started at node 1 then landed on nodes 3 and 8 to charge itself, and then arrived at node 10 if the system assigned this routing plan to it.

To simplify the routing problem, we ignored the lifting and landing time and energy costs in this routing problem, since they were negligible when we considered a long-term journey by the UAVs. However, it was still not simply a problem of finding the shortest path, since there was a trade-off between distance and energy charging for any path. For example, even if the length of a path was short, the total mission time would still be long because the charging efficiencies of the nodes along this path were low, and it would take much more time for the UAV to charge along this path. To deal with this issue, we simplified the path planning problem with Assumptions 3.2.1–3.2.3, discussed below.

The path planning problem was equivalent to a mission time minimization problem:

$$\min {\mathrm{T}}_{\mathrm{total}}=\min \left({\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}+{\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}\right)$$

(1)

where ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}$ represents the total flying time and ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}$ represents the total charging time during the mission. ${\mathrm{T}}_{\mathrm{total}}$ is the total mission time that we wanted to minimize. Further, we assumed that the charging efficiency at node x, defined as the fraction of the maximum charging rate of the UAV’s energy storage, was ${\mathsf{\rho }}_{\mathrm{x}}$. If node x had the maximal possible charging rate (i.e., sufficiently large radiation that allowed for the maximal charging rate of the UAV energy storage), ${\mathsf{\rho }}_{\mathrm{x}}=1$, then, in general, $0<{\mathsf{\rho }}_{\mathrm{x}}\le 1$.

We defined the problem formally as follows: given a geographical (local) map with marked n locations (called nodes): ${\mathrm{V}}_{\mathrm{local}}=\left\{{\mathrm{v}}_1,{\mathrm{v}}_2,\dots ,{ \mathrm{v}}_{\mathrm{n}}\right\},$where ${\mathrm{v}}_1={\mathrm{v}}_{\mathrm{S}}$is the start node and ${\mathrm{v}}_{\mathrm{n}}={\mathrm{v}}_{\mathrm{D}}$ is the destination, and ${\mathrm{v}}_2, \dots ,{ \mathrm{v}}_{\mathrm{n}-1}$ are possible landing places on the geographical map. At first, a graph was generated to describe the reachability relation between the nodes $\mathrm{G}=\left({\mathrm{V}}_{\mathrm{local}},{\mathrm{E}}_{\mathrm{reachable}}\right)$. ${\mathrm{E}}_{\mathrm{reachable}}=\left\{{\mathrm{e}}_1,{\mathrm{e}}_2,\dots ,{ \mathrm{e}}_{\mathrm{m}}\right\}$ represents all the edges that follow ${\mathrm{l}}_{{\mathrm{e}}_{\mathrm{i}}}\le {\mathrm{d}}_{\max },$ where ${\mathrm{l}}_{{\mathrm{e}}_{\mathrm{i}}}$ is the length of edge ${\mathrm{e}}_{\mathrm{i}}$, and ${\mathrm{d}}_{\max }$ is the maximum one nonstop flying distance of the UAV with its energy storage maximally charged. Then, the SPU path planning problem could be described as the mission time minimization problem with Equation (1) when graph $\mathrm{G}$ and distance ${\mathrm{d}}_{\max }$ were chosen.

Let ${\mathrm{l}}_{\mathrm{S}-\mathrm{D}}$ denote the straight-line cross node ${\mathrm{v}}_{\mathrm{S}}$ and node ${\mathrm{v}}_{\mathrm{D}}$. Additionally, let ${\mathrm{e}}_{\mathrm{ij}}$ denote the edge between nodes ${\mathrm{v}}_{\mathrm{i}}$ and ${\mathrm{v}}_{\mathrm{j}}$.

Assumption 3.2.1.G is a graph as defined above for all the edges in$E_{reachable}$. Line$l_{S-D}$needs to be the parallel line of the X-coordinate in the local city map. The coordinate of node$v_S$is (0,$Y_S$) and the coordinate of node$v_D$is ($X_D,Y_S$). For any edge in$E_{reachable}$as$e_{ij}$, assume the coordinate of node$v_i$is ($X_i$,$Y_i$) and the coordinate of node$v_j$is ($X_j$,$Y_j$), always keeping$X_j\ge X_i$.

This assumption ensured that the UAV did not fly significantly “back” to nodes away from the destination. Additionally, circulations in paths were also eliminated with this assumption.

For any given path ${\mathrm{P}}_{\mathrm{k}}=\left[{\mathrm{v}}_{\mathrm{S}},{\mathrm{v}}_{{\mathrm{k}}_1},{\mathrm{v}}_{{\mathrm{k}}_2}\dots { \mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}},{ \mathrm{v}}_{\mathrm{D}}\right]$, where ${\mathrm{v}}_{{\mathrm{k}}_1}, \dots ,{ \mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}}$ were possible landing places along the path, ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}$ was fixed, because the speed of the UAV and the distance between any two nodes were fixed in G. As a result, Equation (1) could be written as:

$$\min {\mathrm{T}}_{\mathrm{total}}={\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}+\min \left({\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}\right)$$

(2)

Assuming that the number of nodes in ${\mathrm{P}}_{\mathrm{k}}$ is $\mathrm{L}$, ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}={{\displaystyle \sum }}_{\mathrm{q}=0}^{\mathrm{L}-2}{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{q}}\to {\mathrm{k}}_{\mathrm{q}+1}}}{{\mathrm{d}}_{\max }}$, where ${\mathrm{T}}_{\mathrm{flight}\left(\max \right)}$ is the maximum flying time for a fully charged UAV, ${\mathrm{d}}_{\max }$ is the maximum distance a UAV could fly with a fully charged energy storage, ${\mathrm{d}}_{{\mathrm{k}}_{\mathrm{q}}\to {\mathrm{k}}_{\mathrm{q}+1}}$ represents the distance between node ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{q}}}$ and node ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{q}+1}}$ (${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{q}}}$ is the ${(\mathrm{q}+1)}^{\mathrm{th}}$ node in ${\mathrm{P}}_{\mathrm{k}}$), and the start node ${\mathrm{v}}_{\mathrm{S}}\equiv {\mathrm{v}}_{{\mathrm{k}}_0}$ on path ${\mathrm{P}}_{\mathrm{k}}$. ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}{\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$, where ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ is the charging time in node ${\mathrm{k}}_{\mathrm{i}}$. Subsequently, if the path is known as ${\mathrm{P}}_{\mathrm{k}}$, the total mission time along path ${\mathrm{P}}_{\mathrm{k}}$ can be calculated as:

$${\mathrm{T}}_{\mathrm{total}}={\displaystyle {\displaystyle \sum }_{\mathrm{q}=0}^{\mathrm{L}-2}}{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{q}\to }{\mathrm{k}}_{\mathrm{q}+1}}}{{\mathrm{d}}_{\max }}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}}{{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}}}$$

(3)

where ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ represents the estimated charging time at node ${\mathrm{k}}_{\mathrm{I}}$ when ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}}=1$. Then, the actual charging time for any node ${\mathrm{k}}_{\mathrm{I}}$ could be obtained with ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}=\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}}{{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}}}$. ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}$ depended on the charging time assignment plans along a path ${\mathrm{P}}_{\mathrm{k}}$; e.g., one possible scheme was that we assumed that the UAV was charged only to the level needed to reach the next node on the path. In general, the approximate range of ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}$ could be found with the following formula:

$${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}-{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}{{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}\le {\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}\le {\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\min }}\times \frac{{\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}-{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}{{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}$$

(4)

where ${\mathsf{\rho }}_{\min }$ represents the minimum charging efficiency for all charging places on a local city map, ${\mathsf{\rho }}_{\max }=1$ represents the maximum charging efficiency for all charging places and ${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{\mathsf{\rho }}$ represents the total time to charge the battery from empty to a fully charged state when the charging efficiency was $\mathsf{\rho }.$ ${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}={\mathsf{\rho }}_{\min }{\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\min }}$.

Assuming that ${\mathrm{E}}_{\mathrm{p}}$ was the remaining power in the UAV when the UAV reached the ${\left(\mathrm{p}+1\right)}^{\mathrm{th}}$ node on path ${\mathrm{P}}_{\mathrm{k}}$, ${\mathrm{R}}_{\mathrm{p}}$ represented the charging time to charge the UAV battery to ${\mathrm{E}}_{\mathrm{p}}$ when charging efficiency was ${\mathsf{\rho }}_{\max }$. Assumption 3.2.2 was created to simplify the derivation:

Assumption 3.2.2.$R_1=R_S=T_{charge\left(max\right)}^{{\rho }_{max}}$; if$T_{total-flight}>T_{flight\left(max\right)}$,$R_L=R_D=0$.

Assumption 3.2.2 supposes that we always had a fully charged UAV at the beginning of a mission. If the UAV could not fly from node S to D without recharging, meaning that ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}>{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}$, we assumed that when the UAV reached the destination it had no remaining power. On the other hand, if ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}\le {\mathrm{T}}_{\mathrm{flight}\left(\max \right)}$, we did not need to consider the energy charging time of any path, because we did not need to charge the UAV during its mission.

We also ignored the landing and lifting-off times, as those are typically negligible relative to the flight time between two nodes on a path:

Assumption 3.2.3.

The lifting-off times and the landing times of each mission were assumed to be negligible.

Besides the assumptions above, we needed to claim that, in this study, the energy charging time and distance between nodes were the only two factors we took into consideration. Here, we assumed a fixed flying speed for all solar-powered UAVs in this problem to simplify the path planning problem, since energy and distance are more crucial for solar-powered path planning problems.

In addition, to finding an optimal charging time assignment plan for any path ${\mathrm{P}}_{\mathrm{k}}$, we recorded the sets of reachable nodes for every node ${\mathrm{v}}_{\mathrm{I}}$in ${\mathrm{V}}_{\mathrm{local}}$ as ${\mathrm{U}}_{\mathrm{i}}=\left\{{\mathrm{U}}_1,{ \mathrm{U}}_2, \dots ,{ \mathrm{U}}_{\mathrm{n}}\right\}$. Any node ${\mathrm{v}}_{\mathrm{j}}$ in reachable set ${\mathrm{U}}_{\mathrm{I}}$ should follow ${\mathrm{d}}_{\mathrm{i}\to \mathrm{j}}\le {\mathrm{d}}_{\max }$ and ${\mathrm{X}}_{\mathrm{j}}\ge {\mathrm{X}}_{\mathrm{i}}$. For example, in Figure 2, for the start node ${\mathrm{v}}_{\mathrm{s}}$, which was node one, nodes two, three and four were in the blue circle with the center as node one and the radius as ${\mathrm{d}}_{\max }$. Consequently, the reachable set of node one was ${\mathrm{U}}_1=\left\{2, 3, 5\right\}$.

The routing solution in the statically charging efficiency environment, as described by Equations (1) and (3), was solved in Section 4.

3.3. SPU Routing Problem Definition (Dynamically Charging Efficiency)

In our SPU path planning problem discussed above, we considered the charging efficiency ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }}$of the node ${\mathrm{k}}_{\mathrm{i}}$ as a fixed value. However, in practice, the charging efficiency is a variable and varies during the daytime; thus, it is time dependent: ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }}\left(\mathrm{t}\right)$. Solar radiation/illumination prediction is of great importance for many applications in a modern city [5], and we could see from [5] that the solar radiation in a place can be approximately represented by a curve.

In our study, we adopted the time-dependent solar radiation model called the S. Kaplanis’s model from [5], shown in Figure 3 (the Y-axis label was replaced with charging efficiency, as defined in this paper):

$${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }}\left(\mathrm{t}\right)=({\mathrm{a}}_{1\left({\mathrm{k}}_{\mathrm{i} }\right)}+{\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i} }\right)}\cos \frac{\mathsf{\pi }\left(\mathrm{t}-{\mathrm{t}}_{\mathrm{peak}}\right)}{{\mathrm{t}}_{\mathrm{end}}})$$

(5)

where ${\mathrm{t}}_0=0$ represents the time of 6:00 a.m., ${\mathrm{t}}_{\mathrm{peak}}=360$ is the peak time of radiation, which was noon, and ${\mathrm{t}}_{\mathrm{end}}=720$ represents the end of the daytime, which was 6:00 p.m. (since the unit of time in our study was minutes). Taking the example from [23], where the mean annual solar radiation in the state of Hawaii was shown, we could say that in different places in the same city, radiations are different, so the prediction factors ${\mathrm{a}}_1$ and ${\mathrm{a}}_2$ in Equation (5) were location dependent.

Thus, based on the radiation data collected in any potential landing place, we could generate the radiation prediction formula per Equation (5) with a simulation of the detrending model in [5]. Then, the total mission time along path ${\mathrm{P}}_{\mathrm{k}}$ could be expressed with Equation (6):

$${\mathrm{T}}_{\mathrm{total}}={\displaystyle {\displaystyle \sum }_{\mathrm{q}=0}^{\mathrm{L}-2}}{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{q}\to }{\mathrm{k}}_{\mathrm{q}+1}}}{{\mathrm{d}}_{\max }}+{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{L}-2}}\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}\left({\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}\right)}{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}·\mathrm{dt}}$$

(6)

where ${\mathrm{t}}_{\mathrm{i}0}$ represents the time that the UAV arrived at node ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{t}}_{\mathrm{i}1}$ represents the time the UAV left node ${\mathrm{k}}_{\mathrm{i}}$ at. When ${\mathrm{t}}_{\mathrm{i}0}$ and ${\mathrm{t}}_{\mathrm{i}1}$ for each node were determined, we could find the value of ${\mathrm{T}}_{\mathrm{total}}$ with Equation (6). To solve the total mission time minimization problem described by Equation (6), we studied how these $\mathrm{variables}$ affected the total mission time and tried to find the corresponding algorithm for the dynamically charging efficiency environment in Section 5.

3.4. Mission Arrangement of Problems Definition

The final problem of our study was to plan a UAV delivery system with multiple UAVs operating at the same time. There could be a large number of possible landing places on a city map. We assumed to have a multitude of $\mathrm{C}$ missions $\mathrm{M}=\left\{{\mathrm{m}}_1,{\mathrm{m}}_2, \dots ,{ \mathrm{m}}_{\mathrm{C}}\right\}$to perform every day, and there was a set of UAVs (with B UAVs prepared) $\mathrm{UAV}\left\{{\mathrm{UAV}}_1,{\mathrm{UAV}}_2, \dots ,{ \mathrm{UAV}}_{\mathrm{B}}\right\}$ that could be used to execute these missions.

Considering the coordination of the missions in the whole city, the problem described by Equation (1) was changed to:

$$\min {\mathrm{T}}_{\mathrm{total}-\mathrm{mission}}={\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{c}}}{\mathrm{T}}_{{\mathrm{total}}_{\mathrm{i}}}$$

(7)

where ${\mathrm{T}}_{{\mathrm{total}}_{\mathrm{i}}}$is the total time of the ${\mathrm{i}}^{\mathrm{th}}$ mission, all the missions had to be completed and the total mission time ${\mathrm{T}}_{\mathrm{total}-\mathrm{mission}}$had to be minimized.

In Figure 4, the ten stores were marked with red circles and one hundred possible landing places were marked with blue crosses. The mission destinations were some of the landing places.

In the whole proposed SPU delivery system, we needed to consider the following problems:

• Whether the space in each landing place was limited.
• Whether the local city map needed to be updated from time to time.
• Whether the system should be able to handle emergencies caused by bad weather.

Section 6 further addresses the above problems by our proposed protocols.

4. Proposed Routing Algorithms (Statically Charging Efficiency)

In Section 2, we mentioned that the algorithm in [2,3] did not yield a globally optimal solution. Their CTA plan was not optimal for any given path because, in their assumption, the UAVs always charged their battery to a full battery state. For example, in Figure 5, on the 1-5-9-10 path, let us say that the charging efficiency on node eight was much larger than on node three. When we reached node three, as opposed to fully charging the battery, we only needed to charge enough to reach node eight. This plan was more efficient than plans determined by [2,3]. When we considered both the total flying distance and energy charging efficiency in the path planning problem, even though the path found by the Dijkstra algorithm had a minimized distance (the distance for the Dijkstra algorithm here was from a drone navigation and charging station (DNCS) [2], which was calculated with the summation of flying time and charging time under the scenario of UAVs always charging themselves to a full battery capacity in any station/landing place [3]) of a mission, the path found by this Dijkstra algorithm may not have been the optimal one. For example, in Figure 5, assume that the Dijkstra algorithm found an optimal path 1-5-9-10 for the following graph. However, node nine had a much higher charging efficiency than nodes along path 1-5-9-10. Thus, if we charged the UAV on node nine on path 1-4-5-9-10, the total mission time ${\mathrm{T}}_{\mathrm{total}}$ could be smaller than that on path 1-5-9-10 obtained with Dijkstra’s algorithm. Thus, to find a globally optimal path planning solution for the SPU delivery mission, we needed a principle to traverse all the possible paths and to find the optimal charging time assignment along each such path.

The algorithm proposed in this section found a globally optimal solution of the SPU path planning problem in a statically charging efficiency environment with an acceptable running time complexity (our proposed algorithms worked on local city maps with at most hundreds of nodes).

4.1. Changing Time Assignment Algorithm

The algorithm presented in this section was a simple globally optimal algorithm (GOA) for solving the problem described in Section 3.2. Namely, the solution to this problem was just a combination of the path traverse algorithm with a depth first search (DFS) and the charging time assignment (CTA) algorithm. The CTA algorithm was necessary in the case that a UAV was unable to reach the destination and needed to find the next landing place for charging before completing the mission.

If a UAV could not fly without recharging to the destination from the current node (including the start node), it needed to find the next landing place for energy charging. This was where the CTA algorithm was necessary.

Assuming the UAV at the time was in the ${\mathrm{w}}^{\mathrm{th}}$ node on path ${\mathrm{P}}_{\mathrm{k}}$ ($\mathrm{w} \ge 1$), it could reach the next $\mathrm{h}$ nodes along path ${\mathrm{P}}_{\mathrm{k}}$ (which meant the distances between the ${\mathrm{w}}^{\mathrm{th}}$ node and its next $\mathrm{h}$ nodes were smaller than ${\mathrm{d}}_{\max }$) and the flight time to reach the current ${\mathrm{w}}^{\mathrm{th}}$ node was ${\mathrm{T}}_{\mathrm{flight}-\mathrm{total}\left(\mathrm{w}\right)}$. For any ${\mathrm{w}}^{\mathrm{th}}$node on path ${\mathrm{P}}_{\mathrm{k}}$, assume that the reachable set for the ${\mathrm{w}}^{\mathrm{th}}$node was ${\mathrm{U}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}$, and the set ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)} ({\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}={\mathrm{P}}_{\mathrm{k}}\cap {\mathrm{U}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}$); i.e., ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}$contained only the nodes that were both in ${\mathrm{P}}_{\mathrm{k}}$ and ${\mathrm{U}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}$. Let set ${\mathrm{Z}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\left(\mathrm{w}-1\right)}}$ denote the charging efficiency of the ${\mathrm{w}}^{\mathrm{th}}$ node on path ${\mathrm{P}}_{\mathrm{k}}$ and ${\mathrm{P}}_{\mathrm{X}-\mathrm{Y}}$ denote the part of the path from the ${\mathrm{X}}^{\mathrm{th}}$ node to the ${\mathrm{Y}}^{\mathrm{th}}$ node.

Theorem 4.1.

If there exists$o$  ($o>1)$  such that the  ${(w+o)}^{th}$node on the path$P_k$  has a larger charging efficiency than all the other nodes in  $u_{ P_k\left(w\right)}$, path$P_k$  cannot be an optimal path.

Proof of Theorem 4.1:

Assume the ${(\mathrm{w}+\mathrm{o})}^{\mathrm{th}}$ node $\left(\mathrm{o}>1\right)$has a larger charging efficiency than all the other nodes in ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}$. Additionally, assume two charging plans A and B. Plan A included all the nodes on the path in ${\mathrm{P}}_{\mathrm{k}}$, and in plan B, we skipped the nodes between the ${\mathrm{w}}^{\mathrm{th}}$ node and the ${(\mathrm{w}+\mathrm{o})}^{\mathrm{th}}$ node. Then, ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}\left(\mathrm{planB}\right)}\le {\mathrm{T}}_{\mathrm{total}-\mathrm{flight}\left(\mathrm{planA}\right)}$, because the distance of the straight line between the ${\mathrm{w}}^{\mathrm{th}}$ and ${(\mathrm{w}+\mathrm{o})}^{\mathrm{th}}$ nodes in ${\mathrm{P}}_{\mathrm{k}}$ was never larger than any other path between the nodes that went through other nodes. In fact, ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}\left(\mathrm{planB}\right)}<{\mathrm{T}}_{\mathrm{total}-\mathrm{charge}\left(\mathrm{planA}\right)}$, because we did not waste time to charge the UAV on nodes between the ${\mathrm{w}}^{\mathrm{th}}$ node and the ${(\mathrm{w}+\mathrm{o})}^{\mathrm{th}}$ node, of which the charging efficiencies were lower than the ${(\mathrm{w}+\mathrm{o})}^{\mathrm{th}}$ node. Thus, ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{planB}\right)}<{\mathrm{T}}_{\mathrm{total}\left(\mathrm{planA}\right)}$. Accordingly, plan A with path ${\mathrm{P}}_{\mathrm{k}}$ could not be an optimal path. □

Lemma 4.1.

Assume we are currently in$w^{th}\left(w>1\right)$node in$P_k$. If any$Y$node before the$w^{th}$node along the path$P_k$can reach the$w^{th}$node, and node$Y$has a larger charging efficiency than any other node along the part of$P_k$from$w$to$Y$,$P_k$should be pruned.

Remark 4.1. With the linear greedy algorithm in [21], if we pruned all the paths that satisfied Lemma 4.1, then we would have the following equation to calculate the total charging time on the ${(w+1)}^{th}$ node in $P_k$:

$$\mathrm{If} \max \left\{\mathrm{Z}\left(\mathrm{w}+1\right)\right\}<{\mathsf{\rho }}_{{\mathrm{k}}_{\left(\mathrm{w}\right)}},{\mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{w}+1\right)}={\mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{w}\right)}+\frac{{\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}}{{\mathsf{\rho }}_{{\mathrm{k}}_{\left(\mathrm{w}\right)}}}$$

$$\begin{array}{ll}\mathrm{If} \max \left\{\mathrm{Z}\left(\mathrm{w}+1\right)\right\}& \ge {\mathsf{\rho }}_{{\mathrm{k}}_{\left(\mathrm{w}-1\right)}}\mathrm{or} \mathrm{w}=\mathrm{L}-2,{ \mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{w}+1\right)}\\ & ={\mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{w}\right)}+\frac{\left({\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{w}}\to {\mathrm{k}}_{\mathrm{w}+1}}}{{\mathrm{d}}_{\max }}-{ \mathrm{R}}_{\mathrm{w}}\right)}{{\mathsf{\rho }}_{{\mathrm{k}}_{\left(\mathrm{w}\right)}}}\end{array}$$

(8)

where $\left\{\mathrm{Z}\left(\mathrm{w}+1\right)\right\}$ represents the set of charging efficiencies of nodes in ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}+1\right)}$. If we charged the UAV to full capacity in the ${\mathrm{w}}^{\mathrm{th}}$ node, ${\mathrm{R}}_{\mathrm{w}+1}={\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}-{\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{w}-1}\to {\mathrm{k}}_{\mathrm{w}}}}{{\mathrm{d}}_{\max }}$; otherwise, ${\mathrm{R}}_{\mathrm{w}+1}=0$.

It was shown in [21] that if there existed an optimal solution, the greedy algorithm of Equation (8) could discover it.

We, hereby, present Algorithm 1 to find an optimal CTA on path ${\mathrm{P}}_{\mathrm{k}}$:

Algorithm 1: CTA algorithm.

Input: Path ${\mathrm{P}}_{\mathrm{k}}$, graph G, charging efficiency $\mathrm{Z}$ for all the nodes in ${\mathrm{V}}_{\mathrm{local}}$, reachable set U Output: Minimized total charging time ${\mathrm{T}}_{\mathrm{charge}-\mathrm{total}}$ along ${\mathrm{P}}_{\mathrm{k}}$ and charging time on each node ${\mathrm{k}}_{\mathrm{i}}$, which is ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ Variables: $\mathrm{w}, \mathrm{v}$ 1.    For $\mathrm{w}=1:\mathrm{L}$ − 1 2.    For $\mathrm{v}$ = $\mathrm{w}$-1: 1 3.        If the ${(\mathrm{w}+1)}^{\mathrm{th}}$ node is in ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)} { \mathrm{and} \mathrm{Z}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}+1\right)}$ is larger than any other node           along ${\mathrm{P}}_{\mathrm{v}-\left(\mathrm{w}+1\right)}$, path ${\mathrm{P}}_{\mathrm{k}}$ is dropped (by Lemma 4.1) 4.          Break 5.        End//drop nonoptimal paths 6.    End 7.    If ${\mathrm{P}}_{\mathrm{k}}$ is kept, the next node is the next landing place. Additionally, in this ${(\mathrm{w}+1)}^{\mathrm{th}}$       node, the UAV is charged and allowed to have enough power to reach the ${\left(\mathrm{w}+2\right)}^{\mathrm{th}}$       node. Update ${\mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{w}+1\right)}$ with Equation (8) 8.    End//assign charging time according to Equation (8) 9. End 10.  Return ${\mathrm{T}}_{\mathrm{charge}-\mathrm{total}\left(\mathrm{L}\right)}$ and ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$

Algorithm 1 was applied to find the optimal charging time assignment plan in any statically charging efficiency environment and could be tested by comparing the total charging time it found with what the other algorithms found. Since the CTA algorithm took $\mathrm{O}\left({\mathrm{n}}^2\right)$ and time path generation algorithm took $\mathrm{O}({\mathrm{n}}^2$), the total running time was $\mathrm{O}\left({\mathrm{n}}^4\right)$ if we just processed the CTA algorithm on all the possible paths and then found the optimal path.

Reference [21] proposed a related algorithm for vehicle refueling policies. In their study, cars needed to find a reachable refueling place with the lowest price at each stop. However, in UAV path planning problems, if we wanted to skip a node along a path, it would not be necessary to pass this node, and we could just fly straight to the next node. As a result, in our SPU path planning problem, the path was changed when we decided to skip a node along the path. In addition, our proposed Algorithm 1 was different from the linear greedy algorithm in [21]. For example, assume we had two nodes A and B on path P, and there were some nodes between A and B on path P. In the scenario from [21], to decide if A could reach B, the total distance along A to B on path P should be smaller than ${\mathrm{d}}_{\max }$. However, in our SPU routing problems, A reaching B meant that the cartesian distance between A and B was smaller than ${\mathrm{d}}_{\max }$.

4.2. Globally Optimal Routing Algorithm

As discussed above, nonoptimal paths were dropped during the CTA algorithm. However, it was more efficient to prune these nonoptimal paths during the path generation process instead of the CTA process. As Lemma 4.1 claims, in each branch generated with the DFS algorithm, for each new generated node ${\mathrm{v}}_{\mathrm{i}}$ along path P, if any node along path P before ${\mathrm{v}}_{\mathrm{i}}$ could reach ${\mathrm{v}}_{\mathrm{i}}$ and all the nodes between them had lower charging efficiencies than ${\mathrm{v}}_{\mathrm{i}}$, path P should be pruned. Therefore, we skipped node ${\mathrm{v}}_{\mathrm{i}}$ and found the next node to generate a new path. If there was no path allowing the UAV to finish the delivery mission, the system reported a failure.

Another pruning strategy was to use the result of the Dijkstra algorithm as a guide to evaluate if a generated path could be a candidate for a globally optimal path. Combining this with the strategy based on the CTA, the complete GOA with a pruning strategy for SPU is depicted in Figure 6.

The result of the Dijkstra algorithm represented an upper bound for the total mission time, and at the beginning of the GOA algorithm processing, the result of the Dijkstra algorithm was set as the maximal estimated flying time (MEFT), ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$:${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}={\mathrm{T}}_{\mathrm{total}\left(\mathrm{dijkstra}\right)}$. With this bound in mind, we calculated the total mission time for the current path, ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}$, with the iterative application of Remark 4.1 and general Equation (3).

Algorithm 2 shows the detailed process for GOA with pruning.

Algorithm 2: GOA with pruning.

Input: Graph G, ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}={\mathrm{T}}_{\mathrm{total}\left(\mathrm{dijkstra}\right)}$, charging efficiency Z for all nodes in ${\mathrm{V}}_{\mathrm{local}}$, reachable set U Output: Minimized total mission time ${\mathrm{T}}_{\mathrm{total}}$ and the optimal path ${\mathrm{P}}_{\mathrm{optimal}}$ Variables: $\mathrm{L}, \mathrm{w}, \mathrm{v}$ 1. While all the possible unpruned paths are generated, 2.    Process the DFS path to find the algorithm and initialize the current generated path to        ${\mathrm{P}}_{\mathrm{k}}=\left[{\mathrm{v}}_{\mathrm{S}},{\mathrm{v}}_{{\mathrm{k}}_1},{\mathrm{v}}_{{\mathrm{k}}_2}\dots { \mathrm{v}}_{{\mathrm{k}}_{\mathrm{L}-1}}\right]$, and assume to now be in the ${\mathrm{w}}^{\mathrm{th}}$ node in ${\mathrm{P}}_{\mathrm{k}}$ which is ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}}$and        the next node of ${\mathrm{P}}_{\mathrm{k}}$ is ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}+1}}$. Additionally, the total mission time arriving node ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}}$ is         ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}.$//Path generation 3.    For $\mathrm{v}$ = $\mathrm{w}$ −1: 1, 4.       If the ${(\mathrm{w}+1)}^{\mathrm{th}}$ node is in ${\mathrm{u}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}\right)}{ \mathrm{and} \mathrm{Z}}_{{ \mathrm{P}}_{\mathrm{k}}\left(\mathrm{w}+1\right)}$ is larger than any other node           along ${\mathrm{P}}_{\mathrm{v}-\left(\mathrm{w}+1\right)}$, path ${\mathrm{P}}_{\mathrm{k}}$ is pruned and any path starting with ${\mathrm{P}}_{\mathrm{k}}$ is pruned as well 5.          Break//pruning strategy by Lemma 4.1 6.       End 7.    End 8.    If ${\mathrm{P}}_{\mathrm{k}}$ is unpruned by Lemma 4.1, add ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}+1}}$ to ${\mathrm{P}}_{\mathrm{k}}$ and calculate ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}$ when ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{L}}}$ is added 9.       If ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}>{\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$, 10.          ${\mathrm{P}}_{\mathrm{k}}$ is pruned and any path starting with ${\mathrm{P}}_{\mathrm{k}}$ is pruned as well 11.       Else, if ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}\le {\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$ and ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{L}}}={\mathrm{v}}_{\mathrm{D}}$ 12.          ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}={\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}$ and ${\mathrm{P}}_{\mathrm{optimal}}={\mathrm{P}}_{\mathrm{k}}$ 13.       Else, if ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{candidate}\right)}\le {\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$ and ${\mathrm{v}}_{{\mathrm{k}}_{\mathrm{L}}}\ne {\mathrm{v}}_{\mathrm{D}}$ 14.           Back to step 1 15.       End//find candidate optimal path by comparing 16.    End 17. Return ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$ and ${\mathrm{P}}_{\mathrm{optimal}}$

Algorithm 2 was applied to find the optimal path for the SPU delivery problem in the statically charging efficiency environment and could be tested by comparing the total mission time it found with what other algorithms found. The running time of the unpruned GOA was $\mathrm{O}\left({\mathrm{n}}^4\right)$. However, since we only processed the CTA on unpruned paths, we only needed to calculate the total mission time with Equation (8) on these paths. Consequently, the running time of the pruned GOA could save time exponentially, according to the experiment in Section 8.1. Additionally, the great improvement by the GOA in the experiment also showed that it was necessary to carry this out, because we could significantly save customer time with the GOA. (The GOA could also be applied in other kinds of UAV path planning problems, where the charging efficiency is static. For example, in the path planning problem from [2], when the cost coefficient of time and money was fixed, the charging (cost) efficiency became a fixed value, which was related to the price of energy in different nodes. In such kinds of problems, the GOA could be applied and an optimal solution could be obtained as well.)

5. Proposed Routing Algorithms (Dynamically Charging Efficiency)

As described in Section 3.3, the charging efficiency for any landing place varied with time throughout the day, as shown by Equation (5). Given path ${\mathrm{P}}_{\mathrm{k}}=\left[{\mathrm{v}}_{\mathrm{S}},{\mathrm{v}}_{{\mathrm{k}}_1},{\mathrm{v}}_{{\mathrm{k}}_2}\dots { \mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}},{ \mathrm{v}}_{\mathrm{D}}\right]$, the estimated total mission time could be calculated with Equation (6). However, there were different kinds of variables in Equation (6), such as ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{I}}\right)}$, ${\mathrm{t}}_{\mathrm{i}0}$ and ${\mathrm{t}}_{\mathrm{i}1}$, which caused it to be difficult to find a minimized solution of the total mission time. Thus, in this section, we further analyzed this minimization problem and proposed SPU path planning algorithms regarding the dynamically charging efficiency.

5.1. Advanced Routing Problem Analysis (Dynamical Charging Efficiency)

In any known path ${\mathrm{P}}_{\mathrm{k}}$ = $\left[{\mathrm{v}}_{\mathrm{S}},{\mathrm{v}}_{{\mathrm{k}}_1},{\mathrm{v}}_{{\mathrm{k}}_2}\dots { \mathrm{v}}_{{\mathrm{k}}_{\mathrm{j}}},{ \mathrm{v}}_{\mathrm{D}}\right]$ in graph $\mathrm{G}=\left({\mathrm{V}}_{\mathrm{local}},{\mathrm{E}}_{\mathrm{reachable}}\right)$, represented in Figure 7, the total mission time was calculated in Equation (6), when the charging efficiency was dynamic.

Since, in Equation (5), we could not compare the real-time charging efficiency of two different landing places when we did not know the time $\mathrm{t}$ the UAV would arrive at these charging places, the pruning strategy by Lemma 4.1 in Section 4.1 could not be applied in the SPU path planning problem with dynamically charging efficiency. It was also difficult to compare the total mission time of two paths unless we found a proper charging time assignment algorithm for each path, as well as there being too many variables in Equation (6). Please note that, in Equation (6), the total flying time for a UAV mission, which was ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}}=$ ${{\displaystyle \sum }}_{\mathrm{q}=0}^{\mathrm{L}-2}{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{{\mathrm{k}}_{\mathrm{q}\to }{\mathrm{k}}_{\mathrm{q}+1}}}{{\mathrm{d}}_{\max }}$, was only determined with path ${\mathrm{P}}_{\mathrm{k}}$; the total mission time minimization problem, therefore, was only determined by the total charging time, which was ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}\left({\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}\right)}{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}·\mathrm{dt}}$. Hence, if we found the minimized total charging time with a CTA plan and then compared the total mission time for all possible paths from node S to node D, the mission time minimization problem could be solved.

Firstly, Equation (6) needed to be further analyzed to find a more general description of the total mission time minimization problem. Thus, in this section, we tried to simplify the variables in Equation (6). We started by studying the relationship between ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ and ${\mathrm{t}}_{\mathrm{i}0}$, ${\mathrm{t}}_{\mathrm{i}1}$. In any intermediate node (nodes except for the start node and destination) ${\mathrm{k}}_{\mathrm{i}}$ on path ${\mathrm{P}}_{\mathrm{k}}$, the exact charging time ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ could be calculated with Equation (9):

$${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}=\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}\left({\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}\right)}{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}}$$

(9)

where $\frac{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}}{\left({\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}\right)}$ represents the averaged charging efficiency on node ${\mathrm{k}}_{\mathrm{i}}$ between times ${\mathrm{t}}_{\mathrm{i}0}$ and ${\mathrm{t}}_{\mathrm{i}1}$. Moreover, ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ could also be calculated with ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}={\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}$. Thus, we had the following equation:

$$\frac{{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}\left({\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}\right)}{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}}={\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}$$

$${{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}={\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$$

(10)

If we wanted to charge the UAV from an empty to full battery state, ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}={\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}$. (In addition, with dynamically charging efficiency, we always had ${\mathrm{T}}_{\mathrm{charge}\left({\mathrm{k}}_{\mathrm{i}}\right)}>{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{I}}\right)}$, since ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}\le 1$ and $\frac{{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}}{{\mathrm{t}}_{\mathrm{i}1}-{\mathrm{t}}_{\mathrm{i}0}}<1$).

With Equation (5) and Figure 3, we could find the relations between ${\mathrm{a}}_{1\left({\mathrm{k}}_{\mathrm{i}}\right)}$ and ${\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i}}\right)}$:

(1)

${\mathrm{a}}_{1\left({\mathrm{k}}_{\mathrm{i}}\right)}=0;$ with Figure 3, when t=${\mathrm{t}}_0=0$, we had ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}=0$. Thus, $({\mathrm{a}}_{1\left({\mathrm{k}}_{\mathrm{i} }\right)}+{\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i} }\right)}\cos \frac{\mathsf{\pi }\left(\mathrm{t}-360\right)}{720})={\mathrm{a}}_{1\left({\mathrm{k}}_{\mathrm{i} }\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}=0;$

(2)

Assuming the peak charging efficiency for node ${\mathrm{k}}_{\mathrm{i} }$ was ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$, we always had ${\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i}}\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$. With Figure 3, when t = ${\mathrm{t}}_{\mathrm{peak}}=360$, ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$. Thus, ${\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i}}\right)}\cos \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{peak}}-360\right)}{720})={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$, and then we had ${\mathrm{a}}_{2\left({\mathrm{k}}_{\mathrm{i}}\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$. ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}$ could be estimated with the used data collected in different landing places.

Thus, Equation (10) was changed to Equation (11):

$${{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{i}0}}^{{\mathrm{t}}_{\mathrm{i}1}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}\cos \frac{\mathsf{\pi }\left(\mathrm{t}-360\right)}{720}{\mathrm{d}}_{\mathrm{t}}={\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$$

$${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}0}-360\right)}{720}={\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$$

(11)

If ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ and ${\mathrm{t}}_{\mathrm{i}0}$ were known, we could solve the equation above with Equation (12):

$${\mathrm{t}}_{\mathrm{i}1}=\frac{720}{\mathsf{\pi }}{\sin }^{-1}(\frac{\mathsf{\pi }}{720{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}}({\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}+{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}0}-360\right)}{720}))+360$$

(12)

where ${\mathrm{t}}_{\mathrm{i}0}$ could be obtained with Equation (13):

$${\mathrm{t}}_{\mathrm{i}0}={\mathrm{t}}_{\left(\mathrm{i}-1\right)1}+{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{\left({\mathrm{k}}_{\mathrm{i}-1}\right)\to \left({\mathrm{k}}_{\mathrm{i}}\right)}}{{\mathrm{d}}_{\max }}$$

(13)

where ${\mathrm{t}}_{\left(\mathrm{i}-1\right)1}$ represents the time the UAV left node ${\mathrm{k}}_{\mathrm{i}-1 }$ and ${\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\times \frac{{\mathrm{d}}_{\left({\mathrm{k}}_{\mathrm{i}-1}\right)\to \left({\mathrm{k}}_{\mathrm{i}}\right)}}{{\mathrm{d}}_{\max }}$ represents the total flying time from node ${\mathrm{k}}_{\mathrm{i}-1 }$to node ${\mathrm{k}}_{\mathrm{i} }$. It is obvious that Equation (13) is an iterative equation, where each ${\mathrm{t}}_{\mathrm{i}0}$ and ${\mathrm{t}}_{\mathrm{i}1}$ in node ${\mathrm{k}}_{\mathrm{i} }$ are determined not only by the time the UAV leaves the previous node ${\mathrm{k}}_{\mathrm{i}-1}$, but also by the estimated charging time ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ for node ${\mathrm{k}}_{\mathrm{i} }$. However, if the ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ for any node ${\mathrm{k}}_{\mathrm{i}}$ was known in this problem and ${\mathrm{t}}_{01}=0$, it was possible to calculate the total mission time with Equations (6), (12) and (13), since there were no other variables in Equations (12) and (13).

We defined ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ in Section 3.2 to denote the required charging time on node ${\mathrm{k}}_{\mathrm{i}}$ if ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{t}\right)}\equiv 1$. There was a constraint in regard to the total estimated charging time ${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ in this problem. It claimed that the total energy charged during the mission should be the same as the total energy cost during the mission minus the initial energy the UAV took when the mission started:

$$\begin{array}{ll}{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{L}-2}}{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}& ={\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{L}-2}}({\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}. \frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{peak}\right)}. \frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}0}-360\right)}{720})\\ & ={\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{d}}_{\mathrm{total}}-{\mathrm{d}}_{\max }}{{\mathrm{d}}_{\max }}\end{array}$$

(14)

where ${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{d}}_{\mathrm{total}}-{\mathrm{d}}_{\max }}{{\mathrm{d}}_{\max }}$ shows a constant that shows the total energy needed for charging if a UAV were to fly along path ${\mathrm{P}}_{\mathrm{k}}$. Array [${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_1\right)},{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_2\right)}, \dots ,{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{L}-2}\right)}$] represents the estimated charging time sequence in path ${\mathrm{P}}_{\mathrm{k}}$. Additionally, assume all the missions started from t=0, and when the UAVs arrived at the destination, their batteries were always empty (except in the situation where we did not need to charge the UAV during the voyage, where ${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\times \frac{{\mathrm{d}}_{\mathrm{total}}-{\mathrm{d}}_{\max }}{{\mathrm{d}}_{\max }}\le 0$, where ${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}=0$). The total mission time when the UAV left node ${\mathrm{k}}_{\mathrm{i}}$ was defined as ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i}\right)}$ ($\mathrm{i}\in \left[0,\mathrm{L}-1\right]$). ${\mathrm{T}}_{\mathrm{total}-\mathrm{flight}\left(\mathrm{i}\right)}$ describes the total flying time when the UAV left node ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{T}}_{\mathrm{total}-\mathrm{charge}\left(\mathrm{i}\right)}$ describes the total charging time when the UAV left node ${\mathrm{k}}_{\mathrm{i}}$.

Besides Equation (14), we had another constraint which described the range of the estimated charging time on each node:

$$\max [\left({\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\frac{{\mathrm{d}}_{\left({\mathrm{k}}_{\mathrm{w}}\right)\to \left({\mathrm{k}}_{\mathrm{w}+1}\right)}}{{\mathrm{d}}_{\max }}-{\mathrm{R}}_{\mathrm{w}}\right),0]\le {\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{w}}\right)}\le {\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}$$

(15)

where ${\mathrm{R}}_{\mathrm{w}}$=${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{w}}{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}-\frac{{\mathrm{T}}_{\mathrm{flight}-\mathrm{total}\left(\mathrm{w}\right)}-{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}{{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}{\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}$ shows the remaining energy when the UAV arrived at nodes ${\mathrm{k}}_{\mathrm{w}}$ in ${\mathrm{P}}_{\mathrm{k}}$. ${\mathrm{R}}_{\mathrm{w}}$ was in minutes, which meant it represented how much time we needed to charge the UAV to the current remaining energy if ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i} }\left(\mathrm{t}\right)}\equiv 1$. To make sure we had enough charged energy to reach node ${\mathrm{k}}_{\mathrm{w}+1}$, ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{w}}\right)}$ needed to be larger than Max[$\left({\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\frac{{\mathrm{d}}_{\left({\mathrm{k}}_{\mathrm{w}}\right)\to \left({\mathrm{k}}_{\mathrm{w}+1}\right)}}{{\mathrm{d}}_{\max }}-{\mathrm{R}}_{\mathrm{w}}\right),0$]. ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{w}\right)}$ represents the total mission time when the UAV arrived at node ${\mathrm{k}}_{\mathrm{w}}$.

Let $\left[{\mathrm{X}}_1,{ \mathrm{X}}_2, \dots ,{ \mathrm{X}}_{\mathrm{L}-2}\right]$ denote $\left[{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_1\right)},{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_2\right)}, \dots ,{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{L}-2}\right)}\right]$, $\left[{\mathrm{Y}}_1,{ \mathrm{Y}}_2, \dots ,{ \mathrm{Y}}_{\mathrm{L}-2}\right]$ denote $\left[{\mathrm{T}}_{\mathrm{total}\left(1\right)},{ \mathrm{T}}_{\mathrm{total}\left(2\right)}, \dots ,{ \mathrm{T}}_{\mathrm{total}\left(\mathrm{L}-2\right)}\right]$,$\left[{\mathsf{\rho }}_1,{ \mathsf{\rho }}_2, \dots ,{ \mathsf{\rho }}_{\mathrm{L}-2}\right]$ denote $[{\mathsf{\rho }}_{{\mathrm{k}}_{1 }\left(\mathrm{peak}\right)},{ \mathsf{\rho }}_{{\mathrm{k}}_{2 }\left(\mathrm{peak}\right)}, \dots ,{ \mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{L}-2 }\left(\mathrm{peak}\right)}] ]$ and $\left[{\mathrm{d}}_1,{ \mathrm{d}}_2, \dots ,{ \mathrm{d}}_{\mathrm{L}-1}\right]$ denote $\left[{\mathrm{d}}_{\left({\mathrm{k}}_0\right)\to \left({\mathrm{k}}_1\right)},{ \mathrm{d}}_{\left({\mathrm{k}}_1\right)\to \left({\mathrm{k}}_2\right)}, \dots ,{\mathrm{d}}_{\left({\mathrm{k}}_{\mathrm{L}-2}\right)\to \left({\mathrm{k}}_{\mathrm{L}-1}\right)}\right]$. ${\mathrm{T}}_{\mathrm{cm}}={\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}$, ${\mathrm{D}}_{\mathrm{m}}={\mathrm{d}}_{\max }$, ${\mathrm{F}}_{\mathrm{m}}={\mathrm{T}}_{\mathrm{flight}\left(\max \right)}$ and ${\mathrm{F}}_{\mathrm{i}}={\mathrm{T}}_{\mathrm{flight}\left(\max \right)}\frac{{\mathrm{d}}_{\mathrm{i}}}{{\mathrm{D}}_{\mathrm{m}}}$.

Then, the minimization problem for the total charging time along path ${\mathrm{P}}_{\mathrm{k}}$ could be described as:

$$\begin{array}{ll}\underset{\left[{\mathrm{X}}_1,{ \mathrm{X}}_2, \dots ,{ \mathrm{X}}_{\mathrm{L}-2}\right] \in \left[0,{\mathrm{T}}_{\mathrm{cm}}\right]}{\min }& {\mathrm{Y}}_{\mathrm{L}-1}, \mathrm{s}.\mathrm{t}.{ \mathrm{Y}}_{\mathrm{L}-1}={\mathrm{Y}}_{\mathrm{L}-2}+{\mathrm{F}}_{\mathrm{L}-1};{ \mathrm{Y}}_{\mathrm{L}-2}\\ & =(\frac{720}{\mathsf{\pi }}{\sin }^{-1}(\frac{\mathsf{\pi }}{720{\mathsf{\rho }}_{\mathrm{L}-2}}({\mathrm{x}}_{\mathrm{L}-2}\\ & +{ \mathsf{\rho }}_{\mathrm{L}-2}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{Y}}_{\mathrm{L}-3}+{\mathrm{F}}_{\mathrm{L}-2}-360\right)}{720}))+360)\end{array}$$

(16)

where for any ${\mathrm{X}}_{\mathrm{i}}$ and ${\mathrm{X}}_{\mathrm{j}}$, ${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-2}{\mathrm{X}}_{\mathrm{i}}={\mathrm{T}}_{\mathrm{cm}}\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{L}-1}{\mathrm{d}}_{\mathrm{i}}-{\mathrm{D}}_{\mathrm{m}}}{{\mathrm{D}}_{\mathrm{m}}}$ and $\mathrm{Max}\left[\left({\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{d}}_{\mathrm{j}}}{{\mathrm{D}}_{\mathrm{m}}}-{\mathrm{R}}_{\mathrm{j}}\right),0\right]\le {\mathrm{X}}_{\mathrm{j}}\le {\mathrm{T}}_{\mathrm{cm}}$.

This was a multivariable minimization problem with an iterative expression. Hence, it was quite difficult to find a globally optimal solution for the total mission time minimization problem with dynamic charging efficiency in any ${\mathrm{P}}_{\mathrm{k}}$ path (in [3], such a path planning problem was proved to be an NP-hard problem). Even if we simplified the problem by assigning $\left[{\mathrm{Y}}_1,{ \mathrm{Y}}_2, \dots ,{ \mathrm{Y}}_{\mathrm{L}-2}\right]$ to be integers, the problem became an integer factorization problem, which could not be solved in polynomial time. Thus, in Section 5.2, we found a locally optimal solution for Equation (16) instead; then, we proposed a dynamically greedy algorithm to solve the mission time minimization problem.

5.2. Locally Optimal CTA

If we only focused on one single iteration of Equation (16), the problem would become much simpler (the total time for one iteration of Equation (16) could be presented as ${\mathrm{Y}}_{\mathrm{w}}=\left(\frac{720}{\mathsf{\pi }}{\sin }^{-1}(\frac{\mathsf{\pi }}{720{\mathsf{\rho }}_{\mathrm{w}}}({\mathrm{X}}_{\mathrm{w}}+{ \mathsf{\rho }}_{\mathrm{w}}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{Y}}_{\mathrm{w}-1}+{\mathrm{F}}_{\mathrm{w}}-360\right)}{720}))+360\right)$; in such a problem, the iterative expression was removed). In other words, we considered a locally optimal problem to assign the charging time among two adjacent nodes ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{k}}_{\mathrm{i}+1}$ in path ${\mathrm{P}}_{\mathrm{k}}$. Assume that when the UAV arrived at node ${\mathrm{k}}_{\mathrm{i}}$, it stored ${\mathrm{R}}_{\mathrm{i}}$ energy, while when it left node ${\mathrm{k}}_{\mathrm{i}}$, it stored ${\mathrm{R}}_{\mathrm{i}}^{\prime }$ energy. Let ${\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }$ denote the maximum value that ${\mathrm{R}}_{\mathrm{i}}^{\prime }$ could be, and let ${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }$ denote the minimum value that ${\mathrm{R}}_{\mathrm{i}}^{\prime }$ could be (${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }={\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{d}}_{\mathrm{j}+1}}{{\mathrm{D}}_{\mathrm{m}}}$ when ${\mathrm{R}}_{\mathrm{i}}<{\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{d}}_{\mathrm{j}+1}}{{\mathrm{D}}_{\mathrm{m}}}$ while ${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }={\mathrm{R}}_{\mathrm{i}}$ when ${\mathrm{R}}_{\mathrm{i}}\ge {\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{d}}_{\mathrm{j}+1}}{{\mathrm{D}}_{\mathrm{m}}}$). Let ${\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$ denote the time that the UAV achieves ${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }$ stored energy on node ${\mathrm{k}}_{\mathrm{i}}$, ${\mathrm{t}}_{\mathrm{i}1}$ denote the actual time the UAV leaves node ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$ denote the time that the UAV achieves ${\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }$ stored energy on node ${\mathrm{k}}_{\mathrm{i}}$. Let ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}$ denote the time that the UAV achieves the required additional energy ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }={\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }$ on node ${\mathrm{k}}_{\mathrm{i}+1}$, where ${\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }\in \left[0, {\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\right)$ and ${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\in \left[0, {\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }\right]$.

To calculate the total mission time ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}$ for this local CTA problem, we used ${\mathrm{t}}_{\mathrm{i}1}$ as the variable instead of ${\mathrm{X}}_{\mathrm{i}}$, because it described the locally optimal problem more precisely. Equation (17) shows the way to calculate ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}$ with ${\mathrm{t}}_{\mathrm{i}1}$ (obtained by using Equations (11) and (16)):

$${\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }={\mathsf{\rho }}_{\mathrm{i}}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-{\mathsf{\rho }}_{\mathrm{i}}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}$$

$$\begin{array}{ll}{\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}& =\frac{720}{\mathsf{\pi }}{\sin }^{-1}(\frac{\mathsf{\pi }}{720{\mathsf{\rho }}_{\mathrm{i}+1}}\left({\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime } \right)-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}\\ & -\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}+\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})+360\end{array}$$

(17)

where ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }-\left({\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)$ represents the estimated charging time on node ${\mathrm{k}}_{\mathrm{i}+1}$.

Therefore, the locally mission time minimization problem could be described with the following formula using variable ${\mathrm{t}}_{\mathrm{i}1}$:

$$\begin{array}{ll}\underset{{\mathrm{t}}_{\mathrm{i}1}\in \left[{\mathrm{t}}_{\mathrm{i}1\left(\min \right)},{\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\right]}{\min }& {\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)},\mathrm{ s}.\mathrm{t}.{\mathrm{ T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}=\frac{720}{\pi }{\sin }^{-1}(\frac{\mathrm{\pi }}{720{\mathsf{\rho }}_{\mathrm{i}+1}}\left({\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{stone}}^{\prime }\right)\\ & -\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}\\ & +\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})+360\end{array}$$

(18)

where $\max \left[0,\left({\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{F}}_{\mathrm{i}}}{{\mathrm{F}}_{\mathrm{m}}}-{\mathrm{R}}_{\mathrm{i}}\right)\right]\le {\mathsf{\rho }}_{\mathrm{i}}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-{\mathsf{\rho }}_{\mathrm{i}}\frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}\le \left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right).$ ${\mathrm{t}}_{\mathrm{i}1\left(\max \right)}=\left(\frac{720}{\mathrm{\pi }}{\sin }^{-1}(\frac{\mathrm{\pi }}{720{\mathsf{\rho }}_{\mathrm{i}}}(\left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}\right)+{\mathrm{ \rho }}_{\mathrm{i}}\frac{720}{\mathrm{\pi }}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}0}-360\right)}{720}))+360\right)=\left(\frac{720}{\mathrm{\pi }}{\sin }^{-1}(\frac{\mathrm{\pi }}{720{\mathsf{\rho }}_{\mathrm{i}}}(\left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)+{\mathrm{ \rho }}_{\mathrm{i}}\frac{720}{\mathrm{\pi }}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}))+360\right).$ (When ${\mathrm{R}}_{\mathrm{i}}<{\mathrm{T}}_{\mathrm{cm}}\frac{{\mathrm{F}}_{\mathrm{i}}}{{\mathrm{F}}_{\max }}$, ${\mathrm{t}}_{\mathrm{i}1\left(\min \right)}>{\mathrm{t}}_{\mathrm{i}0}$, otherwise, ${\mathrm{t}}_{\mathrm{i}1\left(\min \right)}={\mathrm{t}}_{\mathrm{i}0}$,).

After processing a derivation on Equation (17), we had:

$$\begin{array}{l}\frac{{\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}}{{\mathrm{d}}_{{\mathrm{t}}_{\mathrm{i}1}}}\\ =\frac{\cos \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720}-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\cos \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}}{\cos (\frac{\mathrm{\pi }}{720{\mathsf{\rho }}_{\mathrm{i}+1}}{\mathrm{X}}_{\mathrm{i}}-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1\left(\min \right)}-360\right)}{720}+\sin \frac{\mathrm{\pi }\left({\mathrm{t}}_{\mathrm{i}1}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})}\end{array}$$

(19)

In Equation (19), we were able to find that three singular points could possibly be the candidate ${\mathrm{t}}_{\mathrm{i}1}$ to minimize the total charging time between nodes ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{k}}_{\mathrm{i}+1}$: ${\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$, ${\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$ and the singular point calculated with the following equation with ${\mathrm{t}}_{\mathrm{i}1}$ as the variable:

$$\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\cos \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}-360\right)}{720}-\cos \frac{\mathsf{\pi }\left({\mathrm{t}}_{\mathrm{i}1}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720}=0$$

(20)

The singular point calculated with Equation (20) was:

$${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}=\frac{720}{\mathsf{\pi }}{\tan }^{-1}\frac{\cos \left(\frac{\mathsf{\pi }\left({\mathrm{F}}_{\mathrm{i}+1}\right)}{720}-\frac{{\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\rho }}_{\mathrm{i}+1}}\right)}{\sin \frac{\mathsf{\pi }\left({\mathrm{F}}_{\mathrm{i}+1}\right)}{720}}+360$$

(21)

Conclusion 5.1.In the locally optimal problem described in Equation(18), when the solution of Equation(21)$t_{i1\left(singular\right)}$was in [$t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$T_{total\left(i,\left(i+1\right)\right)}$had three candidate extreme points, which were$t_{i1\left(singular\right)}$,$t_{i1\left(min\right)}$and$t_{i1\left(max\right)}$; when the solution of Equation(21) was not in [$t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$T_{total\left(i,\left(i+1\right)\right)}$had two extreme points, which were$t_{i1\left(min\right)}$and$t_{i1\left(max\right)}$.

However, we could not make sure which of the candidate extreme points resulted in ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}$ being minimized by the formulas above. Thus, we needed to divide the problem into different situations from 5.1 to 5.3, and further analyze the locally optimal problem.

5.2.1. Situation 5.1 ${\mathsf{\rho }}_{\mathrm{i}}>{\mathsf{\rho }}_{\mathrm{i}+1}$

${\mathrm{x}}_1={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$ is the time the UAV needed to charge to have enough energy to reach node ${\mathrm{k}}_{\mathrm{i}+1}$; ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1}$ is the exact time that the UAV left node ${\mathrm{k}}_{\mathrm{i}}$; ${\mathrm{x}}_3={\mathrm{x}}_2+{\mathrm{F}}_{\mathrm{i}+1}$ is the time the UAV arrived at node ${\mathrm{k}}_{\mathrm{i}+1}$; ${\mathrm{x}}_4$ is the time the UAV charged for to have enough energy for the stored energy of ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }$ to be achieved after a charging process on node ${\mathrm{k}}_{\mathrm{i}+1}$ and ${\mathrm{x}}_4={\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}$. Thus, using Equation (11), the estimated charging time for the UAV to charge itself from ${\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }$ stored energy to ${\mathrm{R}}_{\mathrm{i}}^{\prime }$ energy on node ${\mathrm{k}}_{\mathrm{i}}$ could be calculated with the following equation:

$${\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }={{\displaystyle \int }}_{{\mathrm{x}}_1}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$$

(22)

Additionally, the estimated charging time for the UAV to charge itself with enough stored energy to ${\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }$ on node ${\mathrm{k}}_{\mathrm{i}+1}$ could be calculated with Equation (23):

$${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }-\left({\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)=\left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)-\left({\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)=\left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}^{\prime }\right)={{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$$

(23)

where ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }-\left({\mathrm{R}}_{\mathrm{i}}^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)$ represents how much energy was still needed to charge on node ${\mathrm{k}}_{\mathrm{i}+1}$ to achieve the estimated stored energy ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }$.

Combining Equations (22) and (23), we had:

$${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }=\left({\mathrm{R}}_{\mathrm{i}}{}_{\max }^{\prime }-{\mathrm{R}}_{\mathrm{i}}{}_{\min }^{\prime }\right)={{\displaystyle \int }}_{{\mathrm{x}}_1}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+{{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$$

(24)

Assume in plan A that the UAV leaves node ${\mathrm{k}}_{\mathrm{i}}$at time ${\mathrm{x}}_2$, and then plan B where the UAV leaves node ${\mathrm{k}}_{\mathrm{i}}$at time ${\mathrm{x}}_2{}^{*}$. In plan B, ${\mathrm{x}}_3{}^{*}={\mathrm{x}}_2{}^{*}+{\mathrm{F}}_{\mathrm{i}+1}$ is the time the UAV arrives at node ${\mathrm{k}}_{\mathrm{i}+1}$ and ${\mathrm{x}}_4{}^{*}$ is the time the UAV takes to charge to have enough stored energy of ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }$. ${\mathrm{x}}_2{}^{*}-{\mathrm{x}}_2={\mathrm{x}}_3{}^{*}-{\mathrm{x}}_3$, since we had ${\mathrm{x}}_3={\mathrm{x}}_2+{\mathrm{F}}_{\mathrm{i}+1}$ and ${\mathrm{x}}_3{}^{*}={\mathrm{x}}_2{}^{*}+{\mathrm{F}}_{\mathrm{i}+1}$.

Remark 5.1. In Figure 8, if${\rho }_i>{\rho }_{i+1}$,$x_2<360$and$x_3<360$, and$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$]. We assumed in plan A that$x_2=t_{i2\left(singular\right)}$. Then, in plan B, when$360\ge x_2{}^{*}\ge t_{i1\left(singular\right)}$, we could conclude that$x_4{}^{*}<x_4$.

Proof of Remark 5.1.

$\frac{({\mathsf{\rho }}_{\mathrm{i}}\cos \frac{\mathsf{\pi }\left(\mathrm{x}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}=\frac{{\mathsf{\pi }\mathsf{\rho }}_{\mathrm{i}}}{720}\sin \frac{\mathsf{\pi }\left(\mathrm{x}-360\right)}{720}$, $\frac{({\mathsf{\rho }}_{\mathrm{i}+1}\cos \frac{\mathsf{\pi }\left(\mathrm{x}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}=\frac{{\mathsf{\pi }\mathsf{\rho }}_{\mathrm{i}+1}}{720}\sin \frac{\mathsf{\pi }\left(\mathrm{x}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720}$, $\mathrm{x}+{\mathrm{F}}_{\mathrm{i}+1}<360$and ${\mathsf{\rho }}_{\mathrm{i}}>{\mathsf{\rho }}_{\mathrm{i}+1}$. Then, $\frac{({\mathsf{\rho }}_{\mathrm{i}}\cos \frac{\mathsf{\pi }\left(\mathrm{x}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}>\frac{({\mathsf{\rho }}_{\mathrm{i}+1}\cos \frac{\mathsf{\pi }\left(\mathrm{x}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}$ could be concluded easily. In Figure 8, assume that ${\mathrm{x}}_2{}^{*}$ is quite close to ${\mathrm{x}}_2$; then, ${\mathrm{x}}_3{}^{*}$ is quite close to ${\mathrm{x}}_3{}^{*}$ as well. To make $\frac{{\mathrm{x}}_4}{{\mathrm{d}}_{{\mathrm{x}}_2}}=0$, we had ${\mathrm{x}}_4={\mathrm{x}}_4{}^{*}$. Then, to make ${\mathrm{x}}_4={\mathrm{x}}_4{}^{*}$, the bar area formed by four points (${\mathrm{x}}_2$,0), (${\mathrm{x}}_2{}^{*},0)$, (${\mathrm{x}}_2$, ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2\right)}$), (${\mathrm{x}}_2{}^{*},{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2{}^{*}\right)})$, which can be represented as ${{\displaystyle \int }}_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$ and the bar area formed by four points (${\mathrm{x}}_3$,0), (${\mathrm{x}}_3{}^{*},0)$, (${\mathrm{x}}_3$, ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3\right)}$), (${\mathrm{x}}_3{}^{*},{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3{}^{*}\right)})$, which can be represented as ${{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$, should have the same area. Thus, if ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}$, we had ${{\displaystyle \int }}_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}={{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$ when ${\mathrm{x}}_2{}^{*}$ was quite close to ${\mathrm{x}}_2$. It is obvious that only when ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3\right)}$ did ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}$ hold. Hence, if ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}$, ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3\right)}$. □

The following equation shows the relationship between plan A and plan B:

$$\begin{array}{ll}{\int }_{{\mathrm{x}}_1}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+& {\int }_{{\mathrm{x}}_3}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}\\ & ={\int }_{{\mathrm{x}}_1}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+{\int }_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{\mathrm{*}}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+({\int }_{{\mathrm{x}}_3}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}\\ & -\mathrm{ }{\int }_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{\mathrm{*}}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}})+{\int }_{{\mathrm{x}}_4}^{{\mathrm{x}}_4{}^{\mathrm{*}}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}\end{array}$$

(25)

With Equation (25), we could conclude that if ${{\displaystyle \int }}_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}>{{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$, then ${\mathrm{x}}_4{}^{*}<{\mathrm{x}}_4$.

Then, if we compared ${{\displaystyle \int }}_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$ and ${{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$, since ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2\right)}={\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3\right)}$ and $\frac{({\mathsf{\rho }}_{\mathrm{i}}\cos \frac{\mathsf{\pi }\left({\mathrm{x}}_2{}^{*}-360\right)}{720})}{{\mathrm{d}}_{{\mathrm{x}}_2{}^{*}}}>\frac{({\mathsf{\rho }}_{\mathrm{i}+1}\cos \frac{\mathsf{\pi }\left({\mathrm{x}}_2{}^{*}+{\mathrm{F}}_{\mathrm{i}}-360\right)}{720})}{{\mathrm{d}}_{{\mathrm{x}}_2{}^{*}}}$, we could conclude that ${{\displaystyle \int }}_{{\mathrm{x}}_2}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}>{{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_3{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$ when (${\mathrm{x}}_2{}^{*}-{\mathrm{x}}_2$) was not a quite small value.

Remark 5.2. If${\rho }_i>{\rho }_{i+1}$,$x_2<360$and$x_3<360$, and$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$]. Assume that$x_2=t_{i1\left(singular\right)}$and$x_2{}^{*}<x_2$; similarly,${{\displaystyle \int }}_{x_2{}^{*}}^{x_2}{\rho }_{k_i\left(t\right)}{d}_t>{{\displaystyle \int }}_{x_3{}^{*}}^{x_3}{\rho }_{k_{i+1}\left(t\right)}{d}_t$. Thus,$x_4{}^{*}<x_4$. (The proof is like Remark 5.1.)

Remark 5.3. If${\rho }_i>{\rho }_{i+1}$,$x_2<360$and$x_3<360$, and$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$T_{total\left(i,\left(i+1\right)\right)}$is a convex function in range$[t_{i1\left(min\right)}, t_{i1\left(max\right)}$] due to Remarks 5.1 and 5.2. Thus,$T_{total\left(i,\left(i+1\right)\right)}$should be compared when$t_{i1}=t_{i1\left(min\right)}$and$t_{i1}=t_{i1\left(max\right)}$and find the result with the minimum total charging time.

Proof of Remark 5.3.

If ${\mathsf{\rho }}_{\mathrm{i}}>{\mathsf{\rho }}_{\mathrm{i}+1}$, ${\mathrm{x}}_2<360$ and ${\mathrm{x}}_3<360$, when ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}$ and ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\in [{\mathrm{t}}_{\mathrm{i}1\left(\min \right)},{ \mathrm{t}}_{\mathrm{i}1\left(\max \right)}$], we always had ${\mathrm{x}}_4{}^{*}<{\mathrm{x}}_4$. Thus, the result of ${\mathrm{x}}_4{}^{*}$ obtained a maximum value in plan A. Using Remarks 5.1 and 5.2, $\frac{{\mathrm{x}}_4{}^{*}}{{\mathrm{d}}_{{\mathrm{x}}_2{}^{*}}}>0$, when ${\mathrm{x}}_2{}^{*}<{\mathrm{x}}_2$; $\frac{{\mathrm{x}}_4{}^{*}}{{\mathrm{d}}_{{\mathrm{x}}_2{}^{*}}}<0$, when ${\mathrm{x}}_2{}^{*}>{\mathrm{x}}_2$. (${\mathrm{x}}_4{}^{*}$ had only one singular point, where $\frac{{\mathrm{x}}_4{}^{*}}{{\mathrm{d}}_{{\mathrm{x}}_2{}^{*}}}=0,$and it decreased either when ${\mathrm{x}}_2{}^{*}<{\mathrm{x}}_2$ or ${\mathrm{x}}_2{}^{*}>{\mathrm{x}}_2$). □

Remark 5.4. If${\rho }_i>{\rho }_{i+1}$,$x_2<360$and$x_3<360$, and$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$t_{i1\left(singular\right)}$must be less than$t_{i1\left(min\right)}$or larger than$t_{i1\left(max\right)}$. If$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$, assume$x_2=t_{i1\left(max\right)}$, since$\frac{({\rho }_i\cos \frac{\pi \left(x-360\right)}{720})}{d_x}>\frac{({\rho }_{i+1}\cos \frac{\pi \left(x+F_{i+1}-360\right)}{720})}{d_x}$,${\rho }_{k_i\left(x_2\right)}<{\rho }_{k_{i+1}\left(x_3\right)}$. Additionally, since$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$x_4{}^{*}$is a monotonic function by$x_2{}^{*}$. Then, also by$\frac{({\rho }_i\cos \frac{\pi \left(x-360\right)}{720})}{d_x}>\frac{({\rho }_{i+1}\cos \frac{\pi \left(x+F_{i+1}-360\right)}{720})}{d_x},$ $t_{i1}=t_{i1\left(min\right)}$obtains a minimum$x_4{}^{*}$. Moreover, if$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$, assume$x_2=t_{i1\left(min\right)}$, since$\frac{({\rho }_i\cos \frac{\pi \left(x-360\right)}{720})}{d_x}>\frac{({\rho }_{i+1}\cos \frac{\pi \left(x+F_i-360\right)}{720})}{d_x}$,${\rho }_{k_i\left(x_2\right)}>{\rho }_{k_{i+1}\left(x_3\right)}$. Additionally, since$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$x_4{}^{*}$is a monotonic function by$x_2{}^{*}$. Then, also by$\frac{({\rho }_i\cos \frac{\pi \left(x-360\right)}{720})}{d_x}>\frac{({\rho }_{i+1}\cos \frac{\pi \left(x+F_{i+1}-360\right)}{720})}{d_x},$ $x_2{}^{*}=t_{i1}=t_{i1\left(max\right)}$obtains a minimum$x_4{}^{*}$.

Proof of Remark 5.4.

Additionally, in Figure 8, if we assumed ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}>{\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$ and ${\mathrm{x}}_2={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$, then we had ${\mathrm{x}}_2{}^{*}<{\mathrm{x}}_2$, ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left({\mathrm{x}}_2\right)}<{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left({\mathrm{x}}_3\right)}$ and $\frac{({\mathsf{\rho }}_{\mathrm{i}}\cos \frac{\mathsf{\pi }\left(\mathrm{x}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}>\frac{({\mathsf{\rho }}_{\mathrm{i}+1}\cos \frac{\mathsf{\pi }\left(\mathrm{x}+{\mathrm{F}}_{\mathrm{i}+1}-360\right)}{720})}{{\mathrm{d}}_{\mathrm{x}}}$. Thus, ${{\displaystyle \int }}_{{\mathrm{x}}_2{}^{*}}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}<{{\displaystyle \int }}_{{\mathrm{x}}_3{}^{*}}^{{\mathrm{x}}_3}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$; then, we had ${{\displaystyle \int }}_{{\mathrm{x}}_1}^{{\mathrm{x}}_2}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+{{\displaystyle \int }}_{{\mathrm{x}}_3}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}<{{\displaystyle \int }}_{{\mathrm{x}}_1}^{{\mathrm{x}}_2{}^{*}}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}+{{\displaystyle \int }}_{{\mathrm{x}}_3{}^{*}}^{{\mathrm{x}}_4}{\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}+1}\left(\mathrm{t}\right)}{\mathrm{d}}_{\mathrm{t}}$ (if we still charged the UAV on node ${\mathrm{k}}_{\mathrm{i}+1}$ till ${\mathrm{x}}_4$ and we changed the time the UAV left node ${\mathrm{k}}_{\mathrm{i}}$ from ${\mathrm{x}}_2$ to ${\mathrm{x}}_2{}^{*}$, ${\mathrm{x}}_2{}^{*}<{\mathrm{x}}_2$, the stored energy could extend ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }$. Thus, it would take less time than ${\mathrm{x}}_4$ to charge the UAV to ${\mathrm{R}}_{\mathrm{i}}{}_{\mathrm{store}}^{\prime }$). Hence, ${\mathrm{x}}_4{}^{*}<{\mathrm{x}}_4$. Since ${\mathrm{x}}_4{}^{*}$ is a monotonic function by ${\mathrm{x}}_2{}^{*}$, ${\mathrm{x}}_4{}^{*}$ obtains a minimized value on ${\mathrm{x}}_2{}^{*}={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$. Similarly, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}<{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$, when ${\mathrm{x}}_2{}^{*}={\mathrm{t}}_{\mathrm{i}1}={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$, ${\mathrm{x}}_4{}^{*}$ would obtain a minimized value. □

Remark 5.5. Similarly toFigure 8, if${\rho }_i>{\rho }_{i+1}$,$x_3\ge 360$,$x_2<360$, we would have a similar conclusion to Remarks 5.3 and 5.4: when${\rho }_i>{\rho }_{i+1}$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], the result of should be compared when$t_{i1}=t_{i1\left(min\right)}$and$t_{i1}=t_{i1\left(max\right)}$, and the result with the minimum total charging time should be found, while if$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], when$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$ $t_{i1}=t_{i1\left(min\right)}$has the minimum total charging time, and when$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$,$t_{i1}=t_{i1\left(max\right)}$has the minimum total charging time. (It can be easily proved with the same steps as those in Remarks 5.1.–5.4).

Remark 5.6. In situation 5.1, if$x_2>360$,$t_{i1\left(singular\right)}$could not be larger than 360 when${\rho }_i>{\rho }_{i+1}$. Thus, it is obvious$t_{i1}=t_{i1\left(max\right)}$was the solution of Equation (17) under this situation.

Proof of Remark 5.6${\rho }_{k_i\left(x_2\right)}>{\rho }_{k_{i+1}\left(x_2\right)}$,${\rho }_{k_{i+1}\left(x_2\right)}>{\rho }_{k_{i+1}\left(x_3\right)}$, thus${\rho }_{k_i\left(x_2\right)}>{\rho }_{k_{i+1}\left(x_3\right)}$, similarly,${\rho }_{k_i\left(x_2{}^{*}\right)}>{\rho }_{k_{i+1}\left(x_3{}^{*}\right)}$. Thus,${\int }_{x_2}^{x_2{}^{*}}{\rho }_{k_i\left(t\right)}{d}_t>{\int }_{x_3}^{x_3{}^{*}}{\rho }_{k_{i+1}\left(t\right)}{d}_t$. Hence,$x_4{}^{*}<x_4$, and the more we charged the UAV on node  $k_i$, the less total charging time we obtained.

Operation 5.1. When${\rho }_i>{\rho }_{i+1}$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], we compared$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$and$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$and chose the CTA plan with a smaller$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$, we chose$t_{i1}=t_{i1\left(min\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$, while if$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$, we chose$t_{i1}=t_{i1\left(max\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$.

5.2.2. Situation 5.2 ${\mathsf{\rho }}_{\mathrm{i}}={\mathsf{\rho }}_{\mathrm{i}+1}$

Remark 5.7. When${\rho }_i={\rho }_{i+1}$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$360-t_{i1\left(singular\right)}=\frac{F_{i+1}}{2}$. Assume$x_2=t_{i1\left(singular\right)}$, when$x_2{}^{*}<x_2$, we had$x_4{}^{*}<x_4$. When$x_2{}^{*}>x_2$, we also had$x_4{}^{*}<x_4$. Thus,$T_{total\left(i,\left(i+1\right)\right)}$is a convex function in range$[t_{i1\left(min\right)}, t_{i1\left(max\right)}$], the same as situation 5.1. Hence, if we compared$T_{total\left(i,\left(i+1\right)\right)}\left|t_{i1}=t_{i1\left(min\right)} and T_{total\left(i,\left(i+1\right)\right)}\right|t_{i1}=t_{i1\left(max\right)}$, we found the minimized$T_{total\left(i,\left(i+1\right)\right)}$. (The proof is like Remarks 5.1.–5.3).

Remark 5.8. When${\rho }_i={\rho }_{i+1}$, if$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], and if$t_{i1\left(max\right)}<t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$; if$t_{i1\left(min\right)}>t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$. (The proof is like Remark 5.4.)

Operation 5.2. When${\rho }_i={\rho }_{i+1}$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], we compared$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$and$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$and chose the CTA plan with a smaller$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$, we chose$t_{i1}=t_{i1\left(min\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$, while if$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$, we chose$t_{i1}=t_{i1\left(max\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$.

5.2.3. Situation 5.3 ${\mathsf{\rho }}_{\mathrm{i}}<{\mathsf{\rho }}_{\mathrm{i}+1}$

Remark 5.9. Like Remark 5.6, if${\rho }_i<{\rho }_{i+1}$and$x_3<360$, then$t_{i1\left(singular\right)}\notin [t_{i1\left(min\right)}, t_{i1\left(max\right)}$]. Since${{\displaystyle \int }}_{x_2}^{x_2{}^{*}}{\rho }_{k_i\left(t\right)}{d}_t<{{\displaystyle \int }}_{x_3}^{x_3{}^{*}}{\rho }_{k_{i+1}\left(t\right)}{d}_t$,$t_{i1}=t_{i1\left(min\right)}$resulted in a minimized total charging time$T_{total\left(i,\left(i+1\right)\right)}$.

Remark 5.10. If${\rho }_i<{\rho }_{i+1}$and$x_2>360$, assume$x_2=t_{i1\left(singular\right)}$,$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$],$T_{total\left(i,\left(i+1\right)\right)}$is a convex function in range$[t_{i1\left(min\right)}, t_{i1\left(max\right)}$]. Thus, if we compared$T_{total\left(i,\left(i+1\right)\right)}\left|t_{i1}=t_{i1\left(min\right)} and T_{total\left(i,\left(i+1\right)\right)}\right|t_{i1}=t_{i1\left(max\right)}$, we found the minimized$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(max\right)}<t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$; if$t_{i1\left(min\right)}>t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$.

Remark 5.11. In remark 5.10, if$x_2<360$and$x_3>360$, we obtained the same conclusion as Remark 5.10. In conclusion, when$x_3>360$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], and if we compared$T_{total\left(i,\left(i+1\right)\right)}\left|t_{i1}=t_{i1\left(min\right)} and T_{total\left(i,\left(i+1\right)\right)}\right|t_{i1}=t_{i1\left(max\right)}$, we could find the minimized$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(max\right)}<t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$; if$t_{i1\left(min\right)}>t_{i1\left(singular\right)}$,$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$resulted in a minimized value of$T_{total\left(i,\left(i+1\right)\right)}$.

Operation 5.3.When${\rho }_i<{\rho }_{i+1}$, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], we compared$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$and$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$and chose the CTA plan with a smaller$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$, we chose$t_{i1}=t_{i1\left(min\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$, while if$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$, we chose$t_{i1}=t_{i1\left(max\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$.

Operation 5.4. In different situations, if$t_{i1\left(singular\right)}\in [t_{i1\left(min\right)}, t_{i1\left(max\right)}$], we compared$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(min\right)}$and$T_{total\left(i,\left(i+1\right)\right)}|t_{i1}=t_{i1\left(max\right)}$and chose the CTA plan with a smaller$T_{total\left(i,\left(i+1\right)\right)}$. If$t_{i1\left(singular\right)}>t_{i1\left(max\right)}$, we chose$t_{i1}=t_{i1\left(min\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$, while if$t_{i1\left(singular\right)}<t_{i1\left(min\right)}$, we chose$t_{i1}=t_{i1\left(max\right)}$to minimize$T_{total\left(i,\left(i+1\right)\right)}$.

5.3. Dynamically Greedy CTA Algorithm

Section 5.2 shows all the situations we were met with in the locally optimal CTA problem described by Equation (18). ${\mathrm{t}}_{\mathrm{i}1\left(\min \right)}$ could be calculated in any known path P. However, ${\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$ was not determined in the locally optimal problem. If we applied the CTA algorithm to this problem and used the peak charging efficiency in each node as the statically charging efficiency, we could predict a proper ${\mathrm{t}}_{\mathrm{i}1\left(\max \right)}$.

$\left[{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_1\right)},{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_2\right)}, \dots ,{ \mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{L}-2}\right)}\right]$ (shown in Figure 9) could describe different charging time assignment plans obtained in this problem, as follows in the image shown.

The grey bar means we only charged enough power on the current node to reach the next node, while the black bar means we charged the UAV to a full battery state. If we only modified the two adjacent CTAs ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}$ and ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}+1}\right)}$, and kept ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}+{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}+1}\right)}$ unchanged, the CTA problem among node ${\mathrm{k}}_{\mathrm{i}}$ and ${\mathrm{k}}_{\mathrm{i}+1}$ became a locally optimal problem, as discussed in Section 5.2. Combining Algorithm 1 and Operation 5.4, we proposed a dynamically greedy CTA(DG-CTA) algorithm as an advanced greedy algorithm. The DG-CTA algorithm focuses on CTA problems with dynamically charging efficiency. Please note that if the two adjacent nodes were both assigned charging times to charge the UAV to a full battery state, or both to charge the UAV only enough power to reach the next node, we could not apply any changes, since there was only one CTA plan for the current ${\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}}\right)}+{\mathrm{T}}_{\mathrm{charge}-\mathrm{est}\left({\mathrm{k}}_{\mathrm{i}+1}\right)}$.

Algorithm 3 shows the DG-CTA algorithm which improved Algorithm 1 regarding its dynamically charging efficiency. Assume ${\mathrm{t}}_{\mathrm{i}1}\prime$ is the new operated time that the UAV left node ${\mathrm{k}}_{\mathrm{i}}$ at.

Algorithm 3: DG-CTA algorithm.

Input: Path ${\mathrm{P}}_{\mathrm{k}}$, graph G, peak charging efficiency ${\mathsf{\rho }}_{\mathrm{i}}$ for all the nodes in ${\mathrm{V}}_{\mathrm{local}}$ Output: Total charging time along ${\mathrm{P}}_{\mathrm{k}}$ and the arrival time/leaving time for each node Perform Algorithm 1 and find all stayed nodes in path ${\mathrm{P}}_{\mathrm{k}}$ as a new path ${\mathrm{P}}_{\mathrm{k}}\prime$ (L—length (${\mathrm{P}}_{\mathrm{k}}\prime$)), while all the possible unpruned paths are generated 1. For i = 1:L-2, 2.    If ${\mathsf{\rho }}_{\mathrm{i}}^{\prime }<{\mathsf{\rho }}_{\mathrm{i}+1}^{\prime }$, //Situation 5.2.3 is solved with Operation 5.4 3.         If ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime \in \left[{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime },{ \mathrm{t}}_{\mathrm{i}1\left(\max \right)}^{\prime }\right]$ and          ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}^{\prime }={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }>$ ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime >$ 4.             ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$. 5.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime <{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }{\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$. 6.            ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 7.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}^{\prime }>{\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 8.            ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime$ 9.         End 10.    Else, if ${\mathsf{\rho }}_{\mathrm{i}}^{\prime }>{\mathsf{\rho }}_{\mathrm{i}+1}\prime$, //Situation 5.2.3 is solved with Operation 5.4 11.         If ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime \in [{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime ,{ \mathrm{t}}_{\mathrm{i}1\left(\max \right)}^{\prime }$ and            ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}^{\prime }={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }<{\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 12.            ${\mathrm{t}}_{\mathrm{i}1}^{\prime }={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }$ 13.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime <{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }$ 14.            ${\mathrm{t}}_{\mathrm{i}1}^{\prime }={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 15.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}^{\prime }>{\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 16.            ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime$ 17.         End 18.    Else //Situation 5.2.2 is solved with Operation 5.4 19.         If ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime \in [{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime ,{ \mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ and            ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime >{ \mathrm{T}}_{\mathrm{total}\left(\mathrm{i},\left(\mathrm{i}+1\right)\right)}\prime |{\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 20.            ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 21.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}\prime <{\mathrm{t}}_{\mathrm{i}1\left(\min \right)}^{\prime }$ 22.            ${\mathrm{t}}_{\mathrm{i}1}^{\prime }={\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 23.         Else, if ${\mathrm{t}}_{\mathrm{i}1\left(\sin \mathrm{gular}\right)}^{\prime }>{\mathrm{t}}_{\mathrm{i}1\left(\max \right)}\prime$ 24.            ${\mathrm{t}}_{\mathrm{i}1}\prime ={\mathrm{t}}_{\mathrm{i}1\left(\min \right)}\prime$ 25.         End 26.    End 27. Return $\left[{\mathrm{t}}_{10}^{\prime },{ \mathrm{t}}_{20}^{\prime }, \dots ,{ \mathrm{t}}_{\left(\mathrm{L}-2\right)0}^{\prime }\right], \left[{\mathrm{t}}_{11}^{\prime },{ \mathrm{t}}_{21}^{\prime }, \dots ,{ \mathrm{t}}_{\left(\mathrm{L}-2\right)1}^{\prime }\right]$

Algorithm 3 was applied to find the optimal charging time assignment plan in any dynamically charging efficiency environment and could be tested by comparing the total charging time it found with what the other algorithms found. Algorithm 3 was in $\mathrm{O}\left(\mathrm{n}\right)$. However, since the result of Algorithm 1 was required in Algorithm 3, the actual running time of Algorithm 3 was in $\mathrm{O}\left({\mathrm{n}}^2\right)$.

5.4. Heuristic Algorithm Based on DG-CTA

Although Algorithm 3 did not take too much time compared with the discrete enumerate algorithm (which enumerated the minute assignment plan for the CTA and could not be solved in polynomial time), the traversal of paths using the DFS algorithm combined with the CTA algorithms was still time-consuming, because Lemma 4.1 could not be applied in this scheme. To further reduce the total running time of the algorithm, we proposed a heuristic algorithm for the SPU path planning problem with dynamically charging efficiency. The purpose of this heuristic algorithm was to quickly find candidate optimal paths with smaller distances and larger averaged charging efficiencies in the early path generation process. The first ${\mathrm{NOP}}_{\max }$ paths with the smallest cost among all possible paths were chosen to be put in the candidate path set. The best path we found with the heuristic algorithm was the one with a minimum total mission time among all the paths in the candidate path set.

The cost function of the heuristic algorithm was:

$$\mathrm{h}\left(\mathrm{a}\right)=2({\mathsf{\gamma }}_1{\mathrm{T}}_{\mathrm{total}-\mathrm{flight}\left(\mathrm{est}\right)}+{\mathsf{\gamma }}_2{\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}\frac{{\mathrm{T}}_{\mathrm{total}-\mathrm{flight}\left(\mathrm{est}\right)}-{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}}{{\mathsf{\rho }}_{\mathrm{a}\left(\mathrm{k}\right)}{\mathrm{T}}_{\mathrm{flight}\left(\max \right)}})$$

(26)

where ${\mathsf{\rho }}_{\mathrm{a}\left(\mathrm{k}\right)}$ is the averaged charging efficiency on path ${\mathrm{P}}_{\mathrm{k}}$. ${\mathsf{\gamma }}_1$ and ${\mathsf{\gamma }}_2$ are the weight factors for the heuristic cost function and ${\mathsf{\gamma }}_1+{\mathsf{\gamma }}_2=1$.

Algorithm 4 shows the detail of this heuristic algorithm:

Algorithm 4: DG-CTA-based heuristic path planning algorithm.

Input: Path ${\mathrm{P}}_{\mathrm{k}}$, graph G, peak charging efficiency ${\mathsf{\rho }}_{{\mathrm{k}}_{\mathrm{i}}\left(\mathrm{peak}\right)}$ for all nodes in ${\mathrm{V}}_{\mathrm{local}}$ ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}={\mathrm{T}}_{\mathrm{total}\left(\mathrm{Dijkstra}\right)}$, ${\mathrm{NOP}}_{\max }$ which represents the number of candidate paths Output: Minimized total mission time ${\mathrm{T}}_{\mathrm{total}}$ and optimal path ${\mathrm{P}}_{\mathrm{optimal}}$ Variables: a, ${\mathrm{P}}_{\mathrm{k}}$, i 1. Run DFS path generation algorithm on path ${\mathrm{P}}_{\mathrm{k}}$, calculate the cost for each path using Equation (26) 2. Sort all possible paths by the cost and find the first ${\mathrm{NOP}}_{\max }$ paths as the candidate path set 3. //Build the candidate path set 4. For $\mathrm{i}=1:{\mathrm{NOP}}_{\max }$ 5. Calculate ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{a}\right)}$ for ${\mathrm{i}}^{\mathrm{th}}$ path in the candidate path set using Algorithm 3 6.        If ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{a}\right)}>{\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$, 7.           ${\mathrm{P}}_{\mathrm{k}}$ is pruned. 8.        Else, 9.           ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}={\mathrm{T}}_{\mathrm{total}\left(\mathrm{a}\right)}$ and ${\mathrm{P}}_{\mathrm{optimal}}={\mathrm{P}}_{\mathrm{k}}$ 10.            End 11.         End 12. End//calculate mission time for paths 13. Return ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{MEFT}\right)}$ and ${\mathrm{P}}_{\mathrm{optimal}}$

Algorithm 4 was applied to find the optimal path for the SPU delivery problem in the dynamically charging efficiency environment, and could be tested by comparing the total mission time it found with what other algorithms found. Since Equation (26) required a running time in $\mathrm{O}\left(\mathrm{n}\right)$, the cost and path generation process required a running time of $\mathrm{O}\left({\mathrm{n}}^3\right)$. Thus, let the sorting process take $\mathrm{O}\left({\mathrm{n}}^2\right)$, so the total time cost for the candidate path set building process is in $\mathrm{O}\left({\mathrm{n}}^3\right)$. To find the path with a minimized mission time, it required a running time of $\mathrm{O}\left({\mathrm{NOP}}_{\max }\times {\mathrm{n}}^2\right)$. Since ${\mathrm{NOP}}_{\max }$ was a constant in this problem, the total running time of Algorithm 4 was in $\mathrm{O}\left({\mathrm{n}}^3\right)$. As shown in Section 8.2, the performance of our proposed algorithm was much better than any other existing algorithms (Algorithms 3 and 4 can only be used in solar-powered UAV systems, since the optimization analysis in this section was performed on a solar-powered UAV path planning model).

6. UAV Delivery Mission Arrangement for Urban City

Finally, with the single SPU path planning problem solved, we optimized the problem described in Equation (7). In our study, the city map was divided into different local city maps based on the service area of the stores. Thus, we did not find nearby UAVs as [2] was able to.

However, in practical UAV delivery systems, the landing places on a city map may not always be available for the following reasons:

• The landing spaces may be fully occupied.
• The landing space may not be open at certain times.
• Bad weather may be forecasted.

To manage more than one UAV delivery mission, we needed to consider the factors above. Aiming at the first reason, we needed to give the UAV missions different ranks to determine which UAV was chosen to take the limited space in each crowded landing place. Alongside the problems above, we also needed to update the city map and rerun the path planning algorithm in special cases.

6.1. UAV Mission Arrangement Protocol to Solve Space Limitation

We proposed a rank-based protocol in our SPU delivery system to solve the space limitation problem in landing places. Before introducing the algorithms and protocols, we needed to initialize the landing places’ occupation status $\mathrm{SO}\left(\mathrm{city}\right)=\left\{{\mathrm{SO}}_1,{\mathrm{SO}}_2, \dots ,{ \mathrm{SO}}_{\mathrm{n}}\right\}$ and the UAV mission rank status $\mathrm{SM}\left(\mathrm{city}\right)=\left\{{\mathrm{SM}}_1,{\mathrm{SM}}_2, \dots ,{ \mathrm{SM}}_{\mathrm{m}}\right\}$. n is the number of nodes in the local city map and m is the number of missions for the current daytime. ${\mathrm{SO}}_{\mathrm{i}}\le {\mathrm{L}}_{\mathrm{SO}}$, where ${\mathrm{L}}_{\mathrm{SO}}$ is the maximum number of UAVs that could stay in a landing place (the existing technique from [21] allowed UAVs to recognize different landing points where different specific landing points were drawn in different shapes or colors). The $\mathrm{SO}\left(\mathrm{city}\right)$ is a real-time function sequence, such as $\mathrm{SO}\left(\mathrm{city}\right)=\left\{{\mathrm{SO}}_1\left(\mathrm{t}\right),{\mathrm{SO}}_2\left(\mathrm{t}\right), \dots ,{ \mathrm{SO}}_{\mathrm{n}}\left(\mathrm{t}\right)\right\}.{ \mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)=0$ means that there were no UAV missions that decided to land on it. ${\mathrm{SM}}_{\mathrm{j}}$ represents the current rank of mission j. Assuming we had an initialized estimated arriving time for all missions of ${\mathrm{T}}_{\mathrm{MP}\left(\mathrm{initial}\right)}=\left\{{\mathrm{T}}_{\mathrm{MP}}{}_1,{\mathrm{T}}_{\mathrm{MP}}{}_2, \dots ,{ \mathrm{T}}_{\mathrm{MP}}{}_{\mathrm{m}}\right\}$, and at the end of each stop of any UAV, we updated a new estimated arriving time to be ${\mathrm{T}}_{\mathrm{MP}\left(\mathrm{realtime}\right)}=\left\{{\mathrm{T}}_{\mathrm{MP}}{}_{\left(\mathrm{realtime}\right)1},{\mathrm{T}}_{\mathrm{MP}\left(\mathrm{realtime}\right)}{}_2, \dots ,{ \mathrm{T}}_{\mathrm{MP}\left(\mathrm{realtime}\right)}{}_{\mathrm{m}}\right\}$, ${\mathrm{SM}}_{\mathrm{i}\left(\mathrm{initial}\right)}={\mathrm{T}}_{\mathrm{MP}}{}_1$. The rank of the mission was calculated using Equation (27):

$${\mathrm{SM}}_{\mathrm{i}\left(\mathrm{update}\right)}={\mathsf{\tau }}_1{\mathrm{T}}_{\mathrm{MP}}{}_1+{\mathsf{\tau }}_2\left({\mathrm{T}}_{\mathrm{MP}}{}_1-{\mathrm{T}}_{\mathrm{MP}}{}_{\left(\mathrm{realtime}\right)1}\right)$$

(27)

where ${\mathsf{\tau }}_1$ and ${\mathsf{\tau }}_{2 }$are parameters which balanced the contributions for ranks with the estimated total mission time and the customers’ additional waiting time, ${\mathsf{\tau }}_1+{\mathsf{\tau }}_2=1$. Missions with a higher rank had higher weight to charge on a crowded landing place.

Then, we proposed the following delivery mission arrangement protocol to deal with the space limitation problem in landing places:

Delivery mission arrangement protocol:

1.

First, we applied Algorithm 4 to each mission and built a database of $\mathrm{SO}\left(\mathrm{city}\right)=\left\{{\mathrm{SO}}_1\left(\mathrm{t}\right),{\mathrm{SO}}_2\left(\mathrm{t}\right), \dots ,{ \mathrm{SO}}_{\mathrm{n}}\left(\mathrm{t}\right)\right\}$, ${\mathrm{T}}_{\mathrm{MP}\left(\mathrm{initial}\right)}=\left\{{\mathrm{T}}_{\mathrm{MP}}{}_1,{\mathrm{T}}_{\mathrm{MP}}{}_2, \dots ,{ \mathrm{T}}_{\mathrm{MP}}{}_{\mathrm{m}}\right\}$ and the path information ${\mathrm{P}}_{\mathrm{MP}\left(\mathrm{initial}\right)}=\left\{{\mathrm{P}}_{\mathrm{MP}}{}_1,{\mathrm{P}}_{\mathrm{MP}}{}_2, \dots ,{ \mathrm{P}}_{\mathrm{MP}}{}_{\mathrm{m}}\right\}$.

2.

For any occupation function of ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)$, if in any ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)$ we had ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)>\mathrm{C}$, $\mathrm{C}$ was the capacity of the landing place. Then, we needed to arrange the missions that would arrive at the node i. We then had to find missions $\mathrm{C}$ that had the highest weight and kept their fight mission plan unchanged. Assume ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)=\mathrm{K}$.

3.

For missions with lower weights, other than the first $\mathrm{C}$ missions, we ordered them by their mission weights. Choosing missions from a higher weight to lower weight, we reran Algorithm 4 and pruned the current path ${\mathrm{P}}_{\mathrm{MP}}{}_{\mathrm{j}}$, then found other paths until a path did not pass node I, or when during the landing time of a newly generated path, node i was not fully occupied.

4.

Once the path plan for mission j was updated, we updated ${\mathrm{SM}}_{\mathrm{j}}$, ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right),$ ${\mathrm{T}}_{\mathrm{MP}}{}_{\left(\mathrm{realtime}\right)\mathrm{j}}$. If no possible path could be found for mission j, ${\mathrm{T}}_{\mathrm{MP}}{}_1={\mathrm{T}}_{\mathrm{MP}}{}_1+720$ (minutes), and it was delayed to the next day.

5.

We returned back to step three and changed the mission path plan for the next mission.

6.

After we ensured ${\mathrm{SO}}_{\mathrm{i}}\left(\mathrm{t}\right)\le \mathrm{C}$ at all times of day, we went back to step two and changed the mission path plans for the mission in the next node.

With the rank of the mission introduced, the spaces for landing places could be well scheduled without any conflicts.

6.2. Local City Maps Updating Protocol

The algorithms we proposed before was aimed at local city maps, which are different from the integrated city map. While initializing the integrated city map and the local city maps, we considered each partitioned region (local city map) to be the same as the service region for stores by their delivery companies. However, in the following cases, we needed to update the region partition for stores in a city:

• Some nodes in a local city map could not be reached by the store in charge of that region.
• It always took less time to reach a node from other stores than the store serving the node’s region.
• Some nodes were not available, and the owners reported the situations.

Hence, we proposed the following strategy as a protocol to update the local city map:

• When a mission could not be completed because a landing place could not be reached, we calculated whether stores from the adjacent service regions could reach that landing place, and then updated the local city map of the current store and the store that could reach that node with a minimized time.
• A simulation program wan run to calculate the mission time for each node to fulfil the delivery from a store in the current local city map and other adjacent local city maps. Suppose that node a is currently in a region with service covered by store A, and store B has a shortest mission time to reach the node compared with all other adjacent stores; in that case, this node should be removed from the local city map for store A and added to the local city map for store B.
• When any unavailability of landing places was reported, the local city map temporarily removed it and added it back when it was available again.

6.3. Real-Time Rerouting Protocol

In our proposed SPU delivery system, when a bad weather prediction report was received, the missions related to the corresponding region needed to be rerouted to avoid accidents or to reduce the total mission time.

A bad weather report included the following three conditions:

• The cloudy level exceeds the threshold of the UAV charging process taking a much longer time.
• It is raining in this region and the UAVs are not able to fly.
• The wind speed level exceeds the threshold that UAVs can fly safely with payloads.

When these bad weather reports were received, the related UAVs were rerouted. The landing places these UAVs currently stayed on or the next landing places they would have arrived at could be the start node of a new routing problem. In emergency cases, we applied Algorithm 4 with a smaller $NO{P}_{max}$ to quickly reroute the UAVs.

With the proposed protocols in this section, the delivery missions in the city map were well organized and the problems we were met with in the system could be solved.

7. SPU Delivery Feasibility Study

Before we started the experiments to test the path planning algorithms proposed in this study, the feasibility of solar-powered UAV delivery work needed to be discussed. Since the subject of this paper was the path planning algorithms proposed for solar-powered UAV delivery systems, we did not expand on how to implement solar-powered UAVs for delivery, we just showed the prototype and proof that it is feasible in modern techniques.

As we discussed before, we chose a transformable solar-powered UAV [22] with solar panels because of the following reasons: (1) it is not possible to build runways for fixed-wing UAVs to take off or land on in crowded cities [24] and (2) it is not practical to use solar-powered UAVs with the solar panel expanded [25] to perform delivery works, since the large area taken over by solar panels must highly reduce the flying speed and safety.

However, the transformable solar-powered UAV from [22] was not practical, since it was not proven to be safe if such types of solar-powered UAVs with payloads were processing a transformation in the sky. In addition, it is difficult for a transformable solar-powered UAV in the fixed-wing state to adjust its flying direction. As a result, this kind of UAV is not quite practical for UAV delivery missions, since the purpose of our solar-powered UAV system was to deliver accurately and on time. Hence, we preferred a new type of transformable solar-powered UAV that only harvested solar energy during the landing time with the solar panel expanded. Additionally, when lifting-off and flying in the sky, the solar panel should be folded to guarantee the speed and safety of the flight. This new prototype of solar-powered UAVs was called transformable delivery solar-powered UAVs (TDSPUs). As a result, the energy consuming model of TDSPUs was like that of normal quadcopters [26], and the energy harvesting model was like that of the solar-powered quadcopter from [27].

At the beginning of the feasibility study of TDSPUs, we needed to make sure solar-powered UAVs could take payloads. From [27], we knew that 88 pieces of SunPower C60 solar cells, each with dimensions of $12.5\times 12.5{ \mathrm{cm}}^2$, required an area of $13.75{ \mathrm{m}}^2$ and provided 293.92 W of power in standard test conditions. Such solar panels required a payload of 1.2 kg and an extra 1.13 kg due to the construction constraints.

When looking at the most advanced quadcopter designed to carry payloads this year [28], the total mass of the quadcopter was 10.3 kg, and the payload was 5 kg. With the solar panels installed, the additional payload of this quadcopter was still more than 2.5 kg. Despite innovations in solar panel efficiency in recent years, it would be feasible for TDSPUs to take payloads with a maximum mass of 2.5 kg.

Secondly, we needed to find a proper flight duration for TDSPUs, so that we could perform experiments to test our proposed algorithms. The flight duration of the quadcopter needed to be reduced due to the change of structure, and the full payload it took. Thus, let us assume it was reduced from 55 min to 25 min.

Then, based on [28], the advanced quadcopter used a lithium polymer battery, which has a capacity of $2\times \mathrm{21,000} \mathrm{mAh}=\mathrm{42,000} \mathrm{mAh}$. Assume that the battery charging efficiency is 90%, the MPPT efficiency is also 90% (in [29], the battery charging efficiency and MPPT efficiency for lithium polymer batteries charged by SunPower C60 solar cells were both 90%, whereas the other two efficiency parameters in [29] were only used for fixed-wing types of UAVs, so we did not consider them here), so the power harvested by the battery from the 88 pieces of SunPower C60 solar cells can be calculated as 293.92 $\times$ 90% $\times$ 90% = 238.08 W and the battery voltage is 22.2 V. Then, the total energy that could be generated from this battery was $42\times 22.2\times 3600=\mathrm{3,356,640} \mathrm{J}$. Then, the estimated charging time to charge the UAV from 0 to a full battery state was ${\mathrm{T}}_{\mathrm{charge}\left(\max \right)}^{{\mathsf{\rho }}_{\max }}=$ $\mathrm{3,356,640}÷238.08=234.99 \min$.

Finally, to improve the safety of the delivery mission, we assumed that the UAVs had a speed of 1m/s [4], so that ${\mathrm{d}}_{\max }=1\times 25\times 60=1.5 \mathrm{km}$. Then, we showed an instance of how a UAV finished a delivery mission during the daytime.

Assume our SPU system could assign a UAV with a mission from node one to node four, and the UAV is able to harvest solar power on node two and node three. The distance between these nodes was: ${\mathrm{d}}_{1\to 2}=1;{\mathrm{d}}_{2\to 3}=0.6;{ \mathrm{d}}_{3\to 4}=1.2$. The UAV started its mission at 6:00 am and charged itself on landing places as much as the system assigned.

At time ${\mathrm{t}}_{10}={\mathrm{t}}_0+\frac{1}{1.5}\times 25=16.67$, the UAV would arrive at node two. Then, assume the peak charging efficiency on node two was 0.6, and the system assigned the UAV to charge on node two to obtain only enough energy to reach the next node. Using Equation (12), ${\mathrm{t}}_{11}=\frac{720}{\mathsf{\pi }}{\sin }^{-1}\left(\frac{\mathsf{\pi }}{720\times 0.6}\left(234.99\times \frac{0.1}{1.5}+0.6\times \frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left(16.67-360\right)}{720}\right)\right)+360=111.76$. Then, with Equation (13), ${\mathrm{t}}_{20}={\mathrm{t}}_{11}+\frac{0.6}{1.5}\times 25=121.76$. Assume the peak charging efficiency in node three was 0.9, still using Equation (12), ${\mathrm{t}}_{21}=\frac{720}{\mathsf{\pi }}{\sin }^{-1}\left(\frac{\mathsf{\pi }}{720\times 0.9}\left(234.99\times \frac{1.2}{1.5}+0.9\times \frac{720}{\mathsf{\pi }}\sin \frac{\mathsf{\pi }\left(121.76-360\right)}{720}\right)\right)+360=371.29$. Then, finally, the time the UAV arrives at the destination node four is ${\mathrm{t}}_{30}={\mathrm{t}}_{21}+\frac{1.2}{1.5}\times 25=391.29$. In other words, the UAV could finish its mission at 12:31 pm.

As a result, TDSPUs were feasible under today’s techniques to finish delivery missions. Additionally, we believe that with the development of UAV technology and solar panel technology, such kinds of solar-powered UAVs are likely to become more and more utilitarian in our daily lives.

8. Experimental Results

The experiments in this section were performed on MATLAB. We simulated the process of UAV delivery missions on a local city map built with landing places and stores. For each delivery mission, we recorded the total time cost for any single delivery mission using our proposed algorithms and other existing algorithms. The parameters were chosen from Section 7.

8.1. Simulation of Proposed Algorithms (Statically Charging Efficiency)

As discussed in Section 2, previous studies about UAV delivery path planning problems could be classified as Dijkstra-based algorithms and heuristic algorithms. Since it is difficult to determine the parameters in previously proposed heuristic algorithms for SPU path planning problems in a dynamically charging efficiency environment, we only compared our proposed algorithms with Dijkstra-based UAV path planning algorithms.

The GOA was compared with the following two kinds of Dijkstra-based UAV path planning algorithms according to the total mission time:

1.

Distance-based Dijkstra algorithm which did not apply our optimal CTA algorithm (D-Dij)

(Distance-based Dijkstra means the cost function in the Dijkstra algorithm only concerns flying time, which was mentioned in [8]. A Dijkstra algorithm that did not apply the optimal CTA algorithm meant that the UAVs always charged themselves to a full battery state. In addition, to make the comparison fairer, if the UAV arrived at a node that could reach the destination on path P, the UAV only charged enough power to reach the destination instead of to a full battery state. This operation exhausted the UAV’s battery when it arrived at the destination, just like what we assumed in Section 4).

2.

Weight-based Dijkstra algorithm which did not apply the optimal CTA algorithm (W-Dij)

(Weight-based Dijkstra means the cost function in the Dijkstra algorithm was calculated through the summation of the flying time to reach the next node and the charging time on the next node to compensate for the energy consumed during this flight [8].)

Each time running a case, we generated an $\mathrm{H}\times {\mathrm{H} \mathrm{km}}^2$ local city map. The density of landing places varied from 40 to 60 and 80 per $49{ \mathrm{km}}^2$, with a uniform distribution.

Assuming that ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{GOA}\right)}$ represents the averaged total mission time for the pruned GOA (Algorithm 2), and ${\mathrm{T}}_{\mathrm{total}\left(\mathrm{Dij}\right)}$ represents the Dijkstra algorithm applied in [2], the percentage of time saved by the GOA algorithm was:

$${\mathrm{Ef}}_{\mathrm{GOA}}=\frac{{\mathrm{T}}_{\mathrm{total}\left(\mathrm{Dij}\right)}-{\mathrm{T}}_{\mathrm{total}\left(\mathrm{GOA}\right)}}{{\mathrm{T}}_{\mathrm{total}\left(\mathrm{Dij}\right)}}$$

(28)

The following diagrams describe when H was fixed, the minimum charging efficiency (MCE) in the local city map ranged from 0.1 to 0.9 (nine points total) and the percentage of time saved by the GOA algorithm compared with different kinds of Dijkstra algorithms. The improvement was measured using Equation (28), and we only showed the maximum improvement.

Based on Figure 10, we could conclude that our GOA algorithm showed a remarkable improvement compared with the Dijkstra algorithms. As the MCE increased, the improvement by the GOA compared with D-Dij dropped linearly, while the improvement by the GOA compared with W-Dij would not decrease so quickly when the MCE moved closer to one.

These great improvements happened in special cases, where the following situations occurred:

• None of the paths found by the Dijkstra algorithm were the optimal path.
• The CTA plan by the Dijkstra algorithm was not optimal.
• The variance of charging efficiency along the candidate optimal paths was large.

The improvement by the CTA algorithm contributed more to the experiment results for special cases. In some paths, if we used the CTA plan of the Dijkstra algorithms, UAVs would always charge themselves to a full battery state. However, if the charging efficiency of landing places in the early trip was low, and the charging efficiency of landing places in the later trip was high, our CTA algorithm charged the UAV less early in the trip, and then charged more time later in the trip; then, the total charging efficiency for the mission was improved. If the charging efficiencies for nodes on such paths had a large difference, the optimal CTA plan chosen by our CTA algorithm could obtain a large percentage of improvements when compared with the Dijkstra-based CTA plans. The different path planning results between the GOA and Dijkstra algorithms also contributed to the improvement, since the paths found by the Dijkstra algorithms were sometimes not the optimal path (proved in Section 4). As a result, in the situation that nodes in the candidate optimal paths had large differences (or, in other words, the variance of charging efficiency along the candidate optimal paths was large), the GOA could achieve great improvements when compared with the Dijkstra algorithms. For example, in Figure 10a, the GOA (Algorithm 2) could save 80% mission time compared with Dijkstra-based algorithms when the MCE was 0.1 on a local city map. Let ${\mathrm{T}}_{\mathrm{running}\left(\mathrm{pruned}\right)}$ represent the total running time for the GOA with pruning, and ${\mathrm{T}}_{\mathrm{running}\left(\mathrm{unpruned}\right)}$ represent the GOA without pruning; therefore, the ratio of ${\mathrm{T}}_{\mathrm{running}\left(\mathrm{unpruned}\right)}$ to ${\mathrm{T}}_{\mathrm{running}\left(\mathrm{pruned}\right)}$ would be:

$${\mathrm{ER}}_{\mathrm{pruning}}=\frac{{\mathrm{T}}_{\mathrm{running}\left(\mathrm{unpruned}\right)}}{{\mathrm{T}}_{\mathrm{running}\left(\mathrm{pruned}\right)}}$$

(29)

Then, the ${\mathrm{ER}}_{\mathrm{pruning}}$ using the GOA with pruning compared with the unpruned GOA is shown in the following diagram:

The result from Figure 11 shows that the pruning strategy could exponentially reduce the running time of the GOA when the number of nodes in the local city map increased.

8.2. Simulations of Proposed Algorithms (Dynamically Charging Efficiency)

For the experiment in the dynamically charging efficiency scheme, our proposed DG-CTA algorithm was compared with the CTA plan proposed by the Dijkstra algorithms (always charging the UAV to a full battery state) and the CTA algorithm using Algorithm 1, which we proposed in Section 4. Since this was a comparison between CTA algorithms, the paths needed to be fixed. We chose every ten paths that had the minimized total Dijkstra distance in each topology and applied different CTA algorithms.

The heuristic algorithm combined with the DG-CTA algorithm was compared with W-Dij in this experiment as well.

These special cases from Figure 12 described by the red curve happened when Algorithm 1 had an opposite local CTA plan to Algorithm 3. For example, if the charging efficiency in the first landing place was a little higher than the second node, then Algorithm 1 assigned the UAV to charge to a full battery state. However, since it was early in the morning, the charging efficiency on all the nodes was low. Then, our DG-CTA algorithm assigned the UAV to charge only long enough to have enough energy to reach the next node, because the charging efficiency increased dramatically in the morning, and it was more efficient to charge on the next node. As a result, based on Operation 5.4, in the case of when the UAVs arrived at nodes with higher charging efficiencies in the morning or the late afternoon, the DG-CTA algorithm had different results compared with Algorithm 1; then, the improvement by the DG-CTA was obvious.

Described by the blue curve, the special cases where the DG-CTA had a great improvement compared with the weighted Dijkstra-based CTA plans happened when the DG-CTA had a remarkable improvement compared with Algorithm 1 and the variances of charging efficiency along the paths were large.

The special cases described by the green curve that our heuristic algorithm had a great improvement compared with the weighted Dijkstra algorithm happened when the DG-CTA algorithm already had a great improvement and there was a better path chosen by the heuristic algorithm compared with the Dijkstra algorithms.

In conclusion, although our DG-CTA-based heuristic algorithm was not the globally optimal solution, it could still achieve a large amount of improvement compared with existing CTA algorithms. Moreover, our proposed heuristic algorithm could achieve a great improvement compared to the existing Dijkstra algorithm in dynamically charging efficiency schemes with an acceptable running time. For example, in Figure 12a, Algorithm 4 could save 70% of the mission time compared with the Dijkstra-based algorithms when the MCE was 0.1 on a local city map.

9. Conclusions

In this paper, a solar-powered UAV delivery system was proposed to extend the endurance of UAVs for delivery missions. To find the globally optimal solution to the UAV path planning problem in the statically charging efficiency scheme, we proposed a very efficient pruning-based globally optimal algorithm. When considering the routing solution according to the dynamically charging efficiency environment, we also proposed a DG-CTA-based heuristic algorithm. Finally, we designed mission arrangement protocols to manage the UAV delivery missions and solve system-level issues in the SPU delivery system. Simulation results showed that our proposed algorithms had significant improvements compared with previous algorithms when applied in the SPU schemes.

For future research, it is suggested that solar-powered UAVs should be able to automatically find the landing places during voyages by themselves. In addition, we plan to study the path planning problem and the system model in UAV delivery systems with more dense landing places deployed.