Coverage Path Planning Based on the Optimization Strategy of Multiple Solar Powered Unmanned Aerial Vehicles

Abstract

In some specific conditions, UAVs are required to obtain comprehensive information of an area or to operate in the area in an all-round way. In this case, the coverage path planning (CPP) is required. This paper proposes a solution to solve the problem of short endurance time in the coverage path planning (CPP) problem of multi-solar unmanned aerial vehicles (UAVs). Firstly, the energy flow efficiency based on the energy model is proposed to evaluate the energy utilization efficiency during the operation. Moreover, for the areas with and without obstacles, the coverage path optimization model is proposed based on the undirected graph search method. The constraint equation is defined to restrict the UAV from accessing the undirected graph according to certain rules. A mixed integer linear programming (MILP) model is proposed to determine the flight path of each UAV with the objective of minimizing operation time. Through the simulation experiment, compared with the Boustrophedon Cellular Decomposition method for coverage path planning, it is seen that the completion time is greatly improved. In addition, considering the impact of the attitude angle of the solar powered UAV when turning, the operation time and the total energy flow efficiency are defined as the optimization objective. The bi-objective model equation is established to solve the problem of the CPP. A large number of simulation experiments show that the optimization model in this paper selects different optimization objectives and applies to different shapes of areas to be covered, which has wide applicability and strong feasibility.

1. Introduction

In recent years, with the development of the UAV and its related technologies, the unmanned aerial vehicles (UAV) have been widely used in the fields of national defense military [1,2], post-disaster search and rescue [3], and environmental remote sensing [4,5]. With the improvement of solar cell efficiency, the solar powered UAV has begun to enter people’s vision with the advantages of environmental protection and the ability to achieve long endurance [6,7,8,9,10,11].

In 1974, the first solar powered UAV named “Sunrise I” came out [12,13]. In 2003, the Swiss Federal Institute of Technology Zurich began to develop a small solar powered UAV called “Sky-sailor” and completed a 27-h low-altitude flight in Switzerland in June 2008 [14,15]. In 2005, American AC company successfully tested the electric gliding solar powered UAV called “SoLong”, which was used in the field of civil remote sensing [16]. In 2017, the Swiss Federal Institute of Technology Zurich developed a solar powered UAV called “Atlantik Solar”, which broke the voyage record of light and small solar powered UAVs with a continuous flight of 81.5 h laying a foundation for the practical application of solar powered UAV [17,18].

Path planning is one of the most important problems in UAV operation which can greatly improve the work efficiency. Traditional path planning is point-to-point path planning. That is, after knowing the location of the starting point and the end point, the best path is mapped out to the end point [19]. In some specific cases, such as sensing mapping, agricultural plant protection, and post-disaster search and rescue, UAVs are required to obtain comprehensive information of an area or to operate in the area in an all-round way. In this case, coverage path planning (CPP) is required. CPP refers to the planning of the best route to traverse non-obstacle areas in a certain area. When talking about the method of covering path planning, boustrophedon is a typical method of the coverage path planning [20]. For the obstacle area, Oksanen et al. proposed the trapezoidal segmentation method [21], which is used to divide into sub-areas with or without obstacles through a secant line parallel to the direction of coverage. Then, each sub-area without obstacles is performed. Later, Choset et al. proposed Boustrophedon Cellular Decomposition [22] which can divide an area into fewer cells. After that, Huang et al. proposed a “line-sweep-based decomposition” with variable secant direction [23]. In this respect, researchers considered the turning time of the UAV and determined the coverage direction by minimizing the number of turns [24,25]. In addition, Nedjati et al. divided the entire area into different areas, such as the areas that need to be covered, the areas that do not need to be covered, and the areas that UAVs take off. Then, the entire area is divided into a uniform grid. Finally, the problem is transformed into a coverage problem [26].

For the algorithm of the CPP, different application scenarios may have different characteristics. In recent years, it has mainly included geometric algorithm, grid-based algorithm, column generation algorithm, heuristic algorithm, machine learning algorithm, clustering-based algorithm, and genetic algorithm. Some examples are listed in

Table 1. Different algorithms employed in CPP.

Author	Algorithm	Optimization Goal	Characteristics
Torres et al. [27]	Geometric algorithm: the back-and-forth (BF) algorithm for concave or multiple polygon coverage.	Reduce energy consumption as much as possible.	Decompose a complex region into multiple simple regions.
Cabreira et al. [28]	Geometric algorithm: E-spiral algorithm.	Reduce energy consumption.	Set different speeds at different stages of direct flights.
Valente et al. [29]	Grid-Based algorithm: algorithm for irregular-shaped areas.	Reduce energy consumption.	Discretize the target area into regular squares, establish cost function to reduce rotation.
Cabreira et al. [30]	Grid-Based algorithm: an energy-aware grid-based algorithm.	Minimize the energy consumption.	Combine with pruning technology and improve the speed of algorithm.
Cho et al. [31]	Grid-Based algorithm: a two-phase method of the CPP.	Minimize the completion time.	Can be used for path planning for multiple UAVs with different performance.
Apuroop et al. [32]	Grid-Based algorithm: the actor-critic experience replay (ACER) reinforcement learning algorithm.	Minimize energy consumption.	Cover the maximum area with the least number of shape reconfigurations possible and generates the tiling shapes and the trajectory with minimum overall cost.
Choi et al. [33,34]	Column generation algorithm: algorithm for an aerial imaging mission with multiple UAVs.	Reduce energy consumption.	Decompose the overlay task and track the energy consumption required for each segment.
Modares et al. [35]	Heuristic algorithm: the Lin–Kernighan heuristic algorithm for drones (LKH-D).	Minimize the energy consumption.	Improve the traditional LKH algorithm.
Yu et al. [36]	Heuristic algorithm.	Minimize the time cost.	Take into account the symbiotic relationship between UGV and UAV.
Chen et al. [37]	Heuristic algorithm: ant colony system (ACS)-based algorithm.	Minimize the completion time.	Cooperate multiple heterogeneous UAVs opreation.
Melo et al. [38]	Heuristic algorithm: a distributed execution of the algorithm through fog-edge computing.	Minimize the cost.	Presents a method to perform CPP in 3D environment for UAVs applications.
Biundini et al. [39]	Metaheuristic algorithm: algorithm based on point cloud data to inspect slope and dams structures.	Minimize flight time.	Find the best path to coverage of a determined area respecting the operation’s restrictions.
Theile et al. [40]	Machine learning algorithm: a double deep Q-network.	Maximum efficiency.	Balance limited power budget and coverage task.
Challita et al. [41]	Machine learning algorithm: algorithm based on echo state network	Balance limited power budget.	Solve this problem of the CPP in cellular unmanned aerial vehicle networks.
Chen et al. [42]	Clustering-based algorithm: density-based clustering methods.	Coverage tasks would be carried out correctly and efficiently.	Inspired from density-based clustering methods and obtain approximate optimal point-to-point paths for UAVs.
Chen et al. [43]	Clustering-based algorithm: an adaptive clustering-based algorithm with a symbiotic organisms search-based optimization strategy.	Minimizing the time consumption.	Efficiently settle the path planning problem and generate feasible paths for heterogeneous UAVs.
Tang et al. [44]	Genetic algorithm: an improved Depth-First-Search (DFS) algorithm.	Minimize the number of sub-regions and reduce the sub-regions to be merged.	Reduce the number of turns, increase the coverage rate, and shorten the non-working distance.
Yuan et al. [45]	Genetic algorithm: the good point set algorithm with the help of better crossover operator and mutation operator for the genetic operation.	Reduce energy consumption.	Solve the energy consumption during turning between adjacent tracks, raise the task execution efficiency.

.

However, the above algorithms have some defects: someone ignores the time cost [27,28,30,32,33,34,35,45], the energy consumption is not taken into account incompletely in some conditions [29,31,36,37,39,43,44], and at the same time, some algorithms cannot complete the task when the CPP area has obstacle [38,40,41,42].

For the more complex application fields, due to the limited operating capacity of a single UAV, multiple UAVs are often used to work together to cover an area. Among them, a way to solve the problem is to use different algorithms to transform the CPP problem into point-to-point path planning to find the optimal solution according to different application scenarios and constraints [37,42,43]. In normal cases, traditional electric UAVs are used in the CPP problems, which has the disadvantage of short endurance time. Meanwhile, due to approximating the CPP problem to a point-to-point problem, the optimal solution is only an approximate value which cannot apply in areas with obstacles. Therefore, this paper proposes the coverage path planning problem of applying multiple solar powered UAVs under the condition that the environmental map is known. The innovations are summarized as follows:

(1) Aiming at the areas to be covered with different shapes, an optimization model of coverage path is proposed based on undirected graph search. The shortest total operation time is selected as the optimization objective, and then the optimal flight path of the UAV is solved by defining constraints. This model does not require that there are no obstacles in the working area, which greatly improves the feasibility and practicability.

(2) The shortest total operation time and the highest energy flow efficiency are selected as optimization objectives, and the bi-objective optimization equation is proposed to solve the problem of the CPP. Making use of the characteristics of undirected graph, this method traverses all the schemes, regardless of whether the shape of the map is a convex region or not, which can ensure that the planning path meets the requirements of the lowest energy efficiency and the least operation time without dividing the map. For nonconvex regions, comparing this scheme with the Boustrophedon Cellular Decomposition method, it is obvious that the proposed method in this paper has shorter completion time.

The remainder is as follows. The energy input model, output model, and energy consumption model of solar powered UAV is established in Section 2. For the convex regions, the model is established with the shortest time as the optimization objective, and the effect is verified by simulation in Section 3. In Section 4, when the operation area is broadened to the nonconvex regions, the model is established and verified by the simulation experiment. At the same time, this model is also compared with other coverage path planning schemes to verify the superiority of the scheme. Taking the energy flow efficiency and the shortest operation time as bi-objective optimization, the model is established in the nonconvex regions and the feasibility is verified by simulation experiments in Section 5. Finally, Section 6 presents the conclusions of this work.

2. Energy Model

For fixed-wing UAV, the fuselage is covered with solar panels, which convert the solar energy absorbed by solar radiation into electrical energy. One of the most important problems is the impact of solar powered UAV energy model in the operations.

2.1. Energy Model

In terms of solar radiation models, some researchers have proposed many statistical models of solar radiation based on the solar radiation in different regions of the world [46]. In this section, the ASHRAE clear sky radiation model is selected to build the energy production model of solar powered UAV due to simple calculation and high accuracy.

The solar radiation intensity absorbed by the surface of the solar panel is linearly related to the angle of incidence of solar rays on the wing surface. The solar radiation perpendicular to the horizontal plane can be calculated as:

$$I=I_0(1+0.034\cos \frac{2\pi n_{day}}{365.25})$$

(1)

where $I_0$ denotes the solar constant and the value is $I_0=1367 \mathrm{W}$. $n_{day}$ represents the number of solar days from 1 January.

Then, the direct solar radiation intensity and direct solar scattering intensity absorbed per unit area of the solar panel horizontal plane can be calculated.

$$I_b=I{e}^{-{\tau }_b{m}_r{}^b}$$

(2)

$$I_d=I{e}^{-{\tau }_d{m}_r{}^d}$$

(3)

$$b=1.219-0.043{\tau }_b-1.151{\tau }_d-0.204{\tau }_b{\tau }_d$$

(4)

$$d=0.202+0.852{\tau }_b-0.007{\tau }_d-0.357{\tau }_b{\tau }_d$$

(5)

$$m_r=\frac{1}{\sin {\alpha }_e}$$

(6)

where $b$ represents the air quality index, $d$ represents the scattered air quality index, $m_r$ represents the air mass ratio, ${\alpha }_e$ represents the solar altitude angle, ${\tau }_b$ represents the direct optical depth, and ${\tau }_d$ represents the scattered optical depth. ${\tau }_b$ and ${\tau }_d$ can be obtained from available literature.

The energy absorption model that comprehensively considers dynamics and energetics can be obtained in reference [47,48]. According to the model, the energy production of the solar powered UAV at a given time, place, and attitude can be calculated. Figure 1 shows the positional relationship between the solar rays and the solar powered UAV.

In Figure 1, ${\alpha }_s$ represents the solar azimuth, which is the clockwise angle between the north direction of the target object and the incident direction of the solar rays. ${\alpha }_e$ represents the solar altitude, which is the angle between the ground plane and the incident rays of the sun. From [49], we can calculate ${\alpha }_s$ and ${\alpha }_e$ at a certain location on the Earth.

$$\sin {\alpha }_e=\sin n_{lat}\sin \delta +\cos n_{lat}\cos \delta \cos \omega (t)$$

(7)

$$\sin {\alpha }_s=\frac{\cos \delta \cos \omega (t)}{\cos {\alpha }_e}$$

(8)

$$\delta =0.4093\sin [2\pi (284+n_{day})/365]$$

(9)

$$\omega (t)=0.2618\times (12-t_{local})$$

(10)

where $n_{lat}$ denotes the latitude of the UAV’s flight location, $\omega (t)$ denotes the solar time, and $\delta$ denotes solar declination angle. $t_{local}$ denotes the current time of the flight location,

$$\cos \lambda =S_b\cdot S_{pb}$$

(11)

$$S_b=\left[\begin{array}{c}0\\ 0\\ -1\end{array}\right]$$

(12)

where $S_b$ is the unit vector in the opposite direction of $o_b{z}_b$ axis in the body axes coordinate system, $S_{pb}$ is the unit vector $S_p$ of the sun rays in the body axes coordinate system. Then, it can be calculated as:

$$S_{pb}=R(\varphi )R(\theta )R(\psi )\left[\begin{array}{c}\cos {\alpha }_e\cos {\alpha }_s\\ \cos {\alpha }_e\sin {\alpha }_s\\ -\sin {\alpha }_e\end{array}\right]$$

(13)

Therefore, the cosine of the solar incident angle can be calculated as:

$$\begin{array}{l}\cos (\lambda )=\cos {\alpha }_e\sin {\alpha }_s(\cos \psi \sin \varphi -\cos \varphi \sin \psi \sin \theta )-\cos {\alpha }_e\cos {\alpha }_s(\sin \varphi \sin \psi +\\ \cos \varphi \cos \psi \sin \theta )+\cos \varphi \sin {\alpha }_e\cos \theta \end{array}$$

(14)

According to the ASHRAE clear sky radiation model, the solar radiation intensity absorbed per unit wing area and the energy input power of the solar powered UAV can be calculated as:

$$P_S(\lambda )=I_b\cos \lambda +I_d{\cos }^2(\frac{\varphi }{2})$$

(15)

$$P_{in}={\eta }_{sol}S{P}_S$$

(16)

where $S$ represents the area of the battery plate and ${\eta }_{sol}$ represents its photoelectric conversion efficiency. Then, the total energy absorbed by the solar powered UAV can be calculated as:

$$E_{in}(t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}{P}_{in}dt}={\displaystyle {\int }_{t_1}^{t_2}{\eta }_{sol}S{P}_{S}dt}$$

(17)

2.2. Energy Consumption Model

According to reference [19], the energy consumption model is established for the solar powered UAV. The energy consumption mainly includes the air resistance and airborne equipment which the UAV need to overcome. Therefore, its energy consumption power can be calculated as:

$$P_{out}=\frac{TV}{{\eta }_{prop}{\eta }_{power}}+P_A$$

(18)

where $V$ represents the flight speed, $T$ represents the pulling force of propeller, ${\eta }_{prop}$ represents the propeller efficiency, ${\eta }_{power}$ represents the power system efficiency, and $P_A$ represents the power consumed by airborne operations. The flight speed can be obtained as:

$$L\cos \varphi =mg$$

(19)

$$L=\frac{\rho V^{2}S{C}_L}{2}$$

(20)

$$V=\sqrt{\frac{2mg}{C_L\rho S\cos \varphi }}$$

(21)

where $C_L$ represents the lift coefficient, $\rho$ represents the air density, $\varphi$ represents attack angle of the solar powered UAV. Then, the propeller tension can be calculated:

$$T=D$$

(22)

$$D=\frac{1}{2}\rho V^{2}S{C}_D$$

(23)

$$K=\frac{1}{\epsilon \pi R_a}$$

(24)

where $R_a$ represents the aspect ratio of the solar powered UAV, $C_D$ represents the flight resistance coefficient, $K$ represents the aerodynamic coefficient, and $\epsilon$ represents the Oswald efficiency factor.

Based on the above, the consumed energy of the solar powered UAV can be calculated as:

$$E_{out}(t_1,t_2)={\displaystyle {\int }_{t_1}^{t_2}{P}_{out}dt}={\displaystyle {\int }_{t_1}^{t_2}\frac{TV}{{\eta }_{prop}{\eta }_{power}}dt}$$

(25)

The “Sky-Sailor” solar powered UAV was selected from reference [14] as the UAV’s model in this paper. The main parameters are listed in

Table 2. Parameters of solar UAV model.

Parameters	Value	Unit
Wing area	0.776	m2
Mass	2.6	Kg
Wingspan	3.2	M
Angle of attack	11	°
Aspect ratio	12.9
Solar cell efficiency	0.169
Propeller efficiency	0.85
Oswald factor	0.9

.

2.3. Energy Flow Efficiency

In this paper, an evaluation index of energy flow efficiency is proposed to represent the utilization rate of solar energy of solar powered UAV during flight operations, which can be calculated as:

$${\eta }_{flow}=\frac{E_{out}}{E_{in}}\times 100\%=\frac{P_{out}}{P_{in}}\times 100\%$$

(26)

Therefore, the value of energy flow efficiency is less than 1, and the larger the value, the higher the utilization rate of solar energy during the flight operations. In addition, for the same solar powered UAV, the higher energy flow efficiency means that the battery load is smaller for the operation. When the energy flow efficiency is inefficient, the battery load is a burden of the operation task. Therefore, the energy flow efficiency can be used as the optimization objective to solve the coverage path planning of solar powered UAV.

3. Single-Objective Optimization of Coverage Path in Non-Obstacle Areas

In this paper, the purpose is to cover the known map area based on the boustrophedon method, which avoids the large amount of calculation but has the disadvantage that is not conducive to subsequent task assignment due to the large number of units to be divided. Therefore, this part establishes the basic model and basic constraints of the CPP in non-obstacle areas, and then expands the constraints to improve the strategy of the CPP when obstacles exist in later sections. It is assumed that the flying height of the UAV is constant, and the limitation of battery life is not considered during the flight. Then, we can solve the optimal path of each UAV when the shortest total time is selected as the optimization objective.

3.1. Determination of Converage Flight Direction

Firstly, the coverage direction needs to be determined. In order to ensure the quality of coverage operations, UAVs often need to fly out of the area for turning. The increase in the number of turns often corresponds to the increase in the total path length [50]. In order to reduce the operating time of the UAVs, the turning times must be as small as possible.

Only when the total operation time is considered as the optimization objective, could the direction be chosen that minimizes the number of turns is set as the best coverage flight direction. In the convex polygonal area, the direction perpendicular to the minimum height is the direction with the least number of turns of the polygon.

In order to determine the direction of coverage flight, the minimum width direction of the area to be covered is selected through traversal. The minimum width direction is found by rotating the polygon representing the area to be covered on a two-dimensional coordinate system and measuring its height. Then, the optimal coverage flight direction is determined, that is, the angle from the x axis of the coordinate system to the optimal coverage flight direction. The search process is shown in Figure 2.

The area to be covered can be divided into lines to prepare for the subsequent coverage task assignment when the coverage flight direction is determined. It is assumed that the minimum width of the area to be covered is $h_{\min }$, and the coverage operation width of the UAV is $L$, and then the number of rows $N_l$ of the coverage operation of the UAV can be obtained.

$$N_l=\left[\frac{h_{\min }}{L}\right]$$

(27)

Then, the coverage width is corrected as:

$$d_l=\frac{h_{\min }}{N_l}$$

(28)

3.2. The Optimization Model of the Coverage Path Based on Undirected Graph Search

In this paper, the optimization model of the coverage path is established based on the undirected graph search method. First, a set of parallel lines with a spacing of $d_l$ is intersected with the boundary of the area to be covered to obtain a set of intersection points $P=\{p_i=(x_i,y_i),i=1,2,\cdots ,N-1\}$. Then, the node set $V$ is composed of this set of intersection points and the take-off nodes in the undirected graph $G=(V,E)$ of the UAV, as shown in Figure 3.

Then, each node is numbered in the set, and the number of all nodes is $N$. The take-off node is set to number 1, the nodes related to the first covered row are numbered 2 and 3, the nodes related to the second covered row are numbered 4 and 5, and so on. In the end, each covered row is related to its corresponding even and odd nodes. The edge set $E$ in an undirected graph $G=(V,E)$ is composed of all the line segments connecting $N$ nodes in the graph, thus a completed undirected graph is formed. The undirected graph is shown in Figure 4.

The undirected graph $G$ can be represented by a cost matrix $C$ of $N\times N$. Moreover, the element $C_{ij}$ represents the Euclidean distance between node $i$ and $j$.

In the optimization model of the coverage path based on undirected graph search, the task of several UAVs at the take-off nodes is to traverse the area to be covered divided into rows. The task can be further described as search the edges in the undirected graph using UAVs to ensure that it accesses all nodes in the undirected graph according to the direction of the covered flight. That is, it accesses all edges in the undirected graph that represent the rows to be covered. According to the optimization objective, the objective model is defined. It forces the UAV to fly according to certain rules in the route by defining the constraint conditions, and finally ensures that each row to be covered has a responsible UAV. The optimal order of each UAV is obtained to access the node by solving the equation set representing the optimization problem, that is, the optimal route of each UAV can be finally obtained. The solution process is shown in Figure 5.

Then, according to the objective function and the constraint equation in the optimization model, it can be solved by the mathematical method of mixed integer linear programming. First, the constants and variables involved in the problem are defined. Among them, the element $C_{ij}$ represents the distance between two nodes in the cost matrix $C$, that is, the covered flight cost of the edge $(i,j)$. The binary variable $X_{ij}^k\in \{0,1\}$ represents whether the kth UAV flies from node $i$ to node $j$. $X_{ij}^k=0$ means that the kth UAV cannot fly from node $i$ to node $j$, $X_{ij}^k=1$ means the kth UAV which flies from node $i$ to node $j$. The constant $V_{ij}^k\in ℝ$ is the speed of the kth UAV which fly from node $i$ to node $j$. It is assumed that the speed of each UAV is equal and constant. $m\in \mathbb{N}$ represents the number of the UAVs participating in the coverage task.

Based on the definition of the above variables, the operating time of the kth UAV can be calculated as:

$$T_k={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N\frac{C_{ij}}{V_{ij}^k}}}{X}_{ij}^k, k=1,2,\cdots ,m$$

(29)

In this section, the single objective of the shortest total operation time is taken as the optimization objective, that is, ensuring that the UAV with the longest flight distance according to the optimal path working in the shortest time (the total time to complete the mission). In this problem, an additional variable $V$ is introduced to represent the operating time of the UAV with the longest path, which is shown in Equation (30):

$$\min (V)$$

(30)

Subject to

$${\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N\frac{C_{ij}}{V_{ij}^k}}}{X}_{ij}^k\le V, k=1,2,\cdots ,m$$

(31)

$${\displaystyle {\sum }_{k=1}^m{\displaystyle {\sum }_{i=1}^N{X}_{ij}^k=1,j=2,3,4,\cdots ,N}}$$

(32)

$${\displaystyle {\sum }_{i=1}^N{X}_{ip}^k-{\displaystyle {\sum }_{j=1}^N{X}_{pj}^k}}=0, p=2,3,\cdots ,N, k=1,2,\cdots ,m$$

(33)

$$u_i-u_j+N{\displaystyle {\sum }_{k=1}^M{X}_{ij}^k}\le N-1, i,j=2,3,\cdots ,N$$

(34)

$${\displaystyle {\sum }_{k=1}^M{X}_{i,i+1}^k+}{\displaystyle {\sum }_{k=1}^M{X}_{i+1,i}^k}=1, i=2,4,6,\cdots ,N$$

(35)

$${\displaystyle {\sum }_{k=1}^M{X}_{i,i+1}^k}={\displaystyle {\sum }_{k=1}^M{\displaystyle {\sum }_{j=\{1,3,\cdots \}\backslash \{i+1\}}^N{X}_{i+1,j}^k}},i=2,4,6,\cdots ,N$$

(36)

$${\displaystyle {\sum }_{k=1}^M{X}_{i,i-1}^k}={\displaystyle {\sum }_{k=1}^M{\displaystyle {\sum }_{j=\left\{1\right\}\cup \{2,4,\cdots \}\backslash \{i-1\}}^N{X}_{i-1,j}^k}},i=3,5,7,\cdots ,N$$

(37)

In Equation (31), the left side of the equation represents the operating time of the kth UAV, and the variable $V$ is used to constrain the operating time of each UAV. As shown in Equation (30), the objective function can be expressed as a linear min-max problem. The constraint equation is defined to restrict the UAV from flying according to certain rules.

Constraint 1: The constraint in Equation (32) indicates that each node except the take-off nodes can only be accessed once by one UAV.

Constraint 2: The constraint in Equation (33) ensures that the UAV that arrives at a node is the same as the UAV that leaves the node.

Constraint 3: Each UAV is to take off at take-off node 1 and return to the take-off node. It is ensured that there is no internal sub-circulation path in this path. Equation (34) represents the standard sub-circulation removal constraint where $u_i\in \mathbb{Z}, i=2,3,4,\cdots ,N$.

Constraint 4: The UAV can fly in the best coverage direction and each coverage line can be covered by the UAV. Constraint Equation (35) also needs to be defined.

The constraint of the Equation (35) guarantees that the UAV which accesses one node of a certain coverage row will access the other node of the coverage row. Considering the numbering order of the nodes, the constraint is realized to force the UAV which accesses an even node to access its corresponding next odd node. The UAV which accesses an odd node will access its corresponding previous even node. Therefore, this constraint is crucial to make this method a feasible solution to the coverage path optimization problem.

Constraint 5: The UAV does not pass through the area to be covered along a side that is not parallel to the coverage direction. Two optional constraints need to be added to the optimization model to prevent turning sharply, as shown in Equations (36) and (37).

The solution can be obtained by solving the objective function in Equation (30) and the constraint equations from Equations (31)–(37). That is, a given number of the UAVs can traverse the area to be covered in the undirected graph in the shortest time.

In addition, it should be considered that the object in the objective function (30) is to ensure that the maximum flight distance of the UAV is as short as possible on the premise of completing the mission, so as to reduce the total time consumption. Therefore, if the longest path of the UAV is found, the path of other UAVs will not be minimized. In order to ensure that the path length of all UAVs is minimum, the following iterative solution is designed here. This method uses multiple iterations to ensure that the path of each UAV is optimal. First of all, in the first iteration, the original undirected graph $G$ is used to solve the optimization problem mentioned above. In the second iteration, we remove the UAV assigned to the longest path in the first iteration and all nodes belonging to the path except the take-off node. After updating the parameters, the method in the first iteration is used to solve the problem. It will be stopped until all nodes are traversed. Through the iterative method, it can be ensured that the path of each UAV is optimal. That is, if the objective function is unchanged, the total cost of the coverage task will be reduced to the lowest possible amount.

3.3. Consider the Optimization Model of Take-Off Operation Time

In the previous section, the optimal path was solved under the single objective of the shortest total operation time. The precondition is that all UAVs take off at the same time. However, when the UAV takes off, it needs a certain operation time, including the installation of the UAV battery module, the start-up of the navigation module, and the time for take-off operation. This section considers the impact of take-off time, the number of operators, and the number of the UAVs.

The new constants and variables need to be defined before describing the optimization model of the UAV take-off time. $M\in \mathbb{N}$ represents the total number of available UAVs. $O\in \mathbb{N}$ represents the best number of the UAVs in a coverage operation. The constant $t_s\in ℝ$ represents the take-off operation time of each UAV. The variable $d_k\in ℝ$ represents the preparation time for launching the kth UAV, including the take-off operation time of the kth UAV and the waiting time before the operator operates the UAV. Then, the time $d_k$ of launching the kth UAV can be calculated, as shown in Equation (38).

$$d_k=t_s\left[\frac{k}{O}\right]{\displaystyle {\sum }_{j=1}^N{X}_{1j}^k}, k=1,2,\cdots ,M$$

(38)

The total time $T_k$ of the kth UAV in the coverage operation and the objective function for updating can be calculated as shown in Equations (39) and (40).

$$T_k={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N\frac{C_{ij}}{V_{ij}^k}}}{X}_{ij}^k+d_k,k=1,2,\cdots ,M$$

(39)

$${\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N\frac{C_{ij}}{V_{ij}^k}}}{X}_{ij}^k+d_k\le V,k=1,2,\cdots ,M$$

(40)

The take-off operation time of the kth UAV is described a cost except for covering the flight cost in Equation (38). This part includes the time cost for the take-off operator to prepare for the take-off mission, the take-off operation, and the waiting time. When there is only one operator, the take-off time of each UAV is superimposed. In an operating system with two UAVs, when the operator operates the first UAV to take off, the second UAV is in the waiting state. The operator can operate the second UAV only after the first UAV takes off. When ${\displaystyle {\sum }_{j=1}^N{X}_{1j}^k=0}$, it means $d_k$ is null, and the kth UAV has not left the take-off node 1. That is, the kth UAV is not used in this operation. The following is an example to illustrate the effect of take-off time on the total operation time.

It is supposed that the number of available UAVs is M = 3, the number of operators is O = 1, and the take-off operation time of each UAV is $t_s=10 \min$. In addition, it is assumed that the task is to cover a rectangular area, which is divided into eight coverage lines, and the time required to cover each line is 2.5 min. The additional time cost $d_k$ of each UAV can be calculated as follows: $d_1=10\left[1/1\right]\times 1=10$; $d_2=10\left[2/1\right]\times 1=20$; $d_3=10\left[3/1\right]\times 1=30$. It can be seen that the cumulative take-off operation time of the third UAV is $d_3=30 \min$, which is equal to the total time $(8\times 2.5+10 \min )$ required for one UAV to cover the area. This shows that it is unnecessary to use three UAVs in this mission. In this case, it is not the best solution of using three UAVs. According to the objective function and constraint equations, we can solve the optimal number of the UAVs and the optimal path of each UAV.

3.4. Simulation Experiment

According to the optimization model built above, the simulation experiments are carried out for the two cases with and without considering the take-off operation time after assuming the experimental conditions. In this paper, the simulation experiment was carried out in MATLAB 2018b, which was configured to run memory 4 GBram. The CPU was Intel corei5, and it was solved by Yalmip toolbox.

3.4.1. Simulation Experiment without Considering Take-Off Operation Time

First of all, it was assumed that the operation site was $(125°\mathrm{E},50°\mathrm{N})$, and the operation time was 09:00 a.m. on 1 May 2020. Moreover, it was supposed that the coverage width of the UAV is $L=130 \mathrm{m}$. When the angle of attack of the UAV is $\alpha =2°$, the velocity can be calculated as $V=10.7784 \mathrm{m}/\mathrm{s}$.

According to the shape of the area to be covered, the experiment without considering the take-off operation time was divided into two groups. In the first group of simulation experiments, the coordinates of the boundary of the operation area to be covered were (−2204, −170), (−533, −1227), (92, −453), (11, −283), (72, 76) and (−1691, 390), the coordinates of the UAV take-off node were (−300, −400). In the second group of simulation experiments, the coordinates of the boundary in the operation area to be covered were (1196, −1021), (188, −565), (51, −184), (−238, 877), (1622, 763) and (1504, −523), the coordinates of the UAV take-off node were (−300, −400). Both of the maps of the area to be covered are shown in Figure 6. Different colors represent different paths of the UAVs. The main parameters for evaluating simulation results are listed in

Table 3. Parameters for evaluating simulation results.

Parameters	Symbol	Unit
Optimal coverage direction	$\delta$	rad
Modified coverage width	$d_l$	m
Operation time	$T_k$	min
Total operation completion time	$t$	min
Number of the UAVs	$m$

.

When $m=2,3,4$, the optimal path of each UAV is shown in Table 4.

Table 4. Results of Experiment 1 of single objective optimization.

Map	m	Tk (min)	t (min)	Optimal Path
The first group of simulation experiments.	2	14.43, 14.43	14.43	As shown in Figure 7a.
	3	10.73, 11.07, 10.92	11.07	As shown in Figure 7b.
	4	9.39, 6.08, 6.75, 9.39	9.39	As shown in Figure 7c.
The second group of simulation experiments.	2	16.37, 16.73	16.73	As shown in Figure 8a.
	3	12.92, 12.05, 12.74	12.92	As shown in Figure 8b.
	4	9.73, 10.49, 10.38, 10.71	10.71	As shown in Figure 8c.

Table 4. Results of Experiment 1 of single objective optimization.

Map	m	Tk (min)	t (min)	Optimal Path
The first group of simulation experiments.	2	14.43, 14.43	14.43	As shown in Figure 7a.
	3	10.73, 11.07, 10.92	11.07	As shown in Figure 7b.
	4	9.39, 6.08, 6.75, 9.39	9.39	As shown in Figure 7c.
The second group of simulation experiments.	2	16.37, 16.73	16.73	As shown in Figure 8a.
	3	12.92, 12.05, 12.74	12.92	As shown in Figure 8b.
	4	9.73, 10.49, 10.38, 10.71	10.71	As shown in Figure 8c.

Figure 7. Experiment 1 of single objective optimization without considering $T_k$: (a) The optimal path of the UAV when m = 2; (b) the optimal path of the UAV when m = 3; (c) the optimal path of the UAV when m = 4. Among them, the different colors lines present different optimal path for different UAVs.

Figure 7. Experiment 1 of single objective optimization without considering $T_k$: (a) The optimal path of the UAV when m = 2; (b) the optimal path of the UAV when m = 3; (c) the optimal path of the UAV when m = 4. Among them, the different colors lines present different optimal path for different UAVs.

Figure 8. Experiment 2 of single objective optimization without considering $T_k$: (a) The optimal path of the UAV when m = 2; (b) the optimal path of the UAV when m = 3; (c) the optimal path of the UAV when m = 4. Among them, the different colors lines present different optimal path for different UAVs.

Figure 8. Experiment 2 of single objective optimization without considering $T_k$: (a) The optimal path of the UAV when m = 2; (b) the optimal path of the UAV when m = 3; (c) the optimal path of the UAV when m = 4. Among them, the different colors lines present different optimal path for different UAVs.

In the first group of experiments, it was obtained that the optimal coverage flight direction of the UAV was $\delta =1.75 \mathrm{rad}$, the modified coverage width was $d_l=126.49 \mathrm{m}$. In the second group of experiments, it was obtained that the optimal coverage flight direction of the UAV was $\delta =2.88 \mathrm{rad}$, the modified coverage width was $d_l=126.24 \mathrm{m}$.

Then, according to the equations in Section 2, the energy production power $P_{in}$ and energy consumption power $P_{out}$ were calculated for the solar UAV. Thus, the energy flow efficiency flow ${\eta }_{flow}$ could be obtained.

It was assumed that the power consumption was set as $P_A=60 \mathrm{W}$ of the UAV in the coverage operation, the energy production power $P_{in}$ and energy consumption power $P_{out}$ of each solar powered UAV could be calculated as

$$P_{in}=72.4626 \mathrm{W}, P_{out}=68.3057 \mathrm{W}$$

(41)

3.4.2. Simulation Experiment Considering Take-Off Operation Time

It was assumed that the coverage width of the UAV was set to $L=200 \mathrm{m}$, and other experimental conditions were the same as previous experiments. The main parameters for evaluating simulation results are listed in

Table 5. Parameters for evaluating simulation results.

Parameters	Symbol	Unit
Optimal coverage direction	$\delta$	rad
Modified coverage width	$d_l$	m
Operation time	$T_k$	min
Total operation completion time	$t$	min
Take-off operation time	ts	min
Optional number of the UAVs	$m$
Number of operators	$O$
Available number of the UAVs	$M$

.

The coordinates of the boundary in the operation area to be covered were (1196, −1021), (188, −565), (51, −184), (−238, 877), (1622, 763) and (1504, −523), the coordinates of the UAV take-off nodes were (−300, −400). When $O$ is different from $M$, the results of the comparison experiment considering $T_k$ is shown in Table 6. Among them, different colors represent different paths of the UAVs.

Table 6. Results of comparison experiment considering $T_k$.

	Condition	${\mathit{d}}_{\mathit{k}}$	${\mathit{t}}_{\mathit{k}}$	${\mathit{T}}_{\mathit{k}}$	t	m	Optimal Path
M = 2	O = 1, ts = 0	0, 0	20.78, 20.94	20.78, 20.94	20.94	2	As shown in Figure 9a.
	O = 1, ts = 6	6, 12	23.11, 17.79	29.11, 29.79	29.79	2	As shown in Figure 9b.
	O = 2, ts = 6	6, 6	20.94, 20.78	26.94, 26.78	26.94	2	As shown in Figure 9c.
M = 3	O = 1, ts = 6	6, 12, 18	20.94, 14.39, 9.47	26.94, 26.39, 27.47	27.47	3	As shown in Figure 10a.
	O = 1, ts = 10	10, 20	25.15, 15.29	35.15, 35.29	35.29	2	As shown in Figure 10b.
	O = 1, ts = 6	6, 6, 12	17.66, 17.66, 10.59	23.66, 23.66, 22.59	23.66	3	As shown in Figure 10c.
M = 4	O = 1, ts = 6	6, 12, 18	21.86, 14.44, 9.17	27.86, 26.44, 27.17	27.86	3	As shown in Figure 11a.
	O = 2, ts = 6	6, 6, 12, 12	16.56, 16.43, 9.47, 10.10	22.56, 22.43, 21.47, 22.10	22.56	4	As shown in Figure 11b.
M = 5	O = 2, ts = 6	6, 6, 12, 12	16.49, 16.56, 10.10, 10.59	22.49, 22.56, 22.10, 22.59	22.59	4	As shown in Figure 12.

Table 6. Results of comparison experiment considering $T_k$.

	Condition	${\mathit{d}}_{\mathit{k}}$	${\mathit{t}}_{\mathit{k}}$	${\mathit{T}}_{\mathit{k}}$	t	m	Optimal Path
M = 2	O = 1, ts = 0	0, 0	20.78, 20.94	20.78, 20.94	20.94	2	As shown in Figure 9a.
	O = 1, ts = 6	6, 12	23.11, 17.79	29.11, 29.79	29.79	2	As shown in Figure 9b.
	O = 2, ts = 6	6, 6	20.94, 20.78	26.94, 26.78	26.94	2	As shown in Figure 9c.
M = 3	O = 1, ts = 6	6, 12, 18	20.94, 14.39, 9.47	26.94, 26.39, 27.47	27.47	3	As shown in Figure 10a.
	O = 1, ts = 10	10, 20	25.15, 15.29	35.15, 35.29	35.29	2	As shown in Figure 10b.
	O = 1, ts = 6	6, 6, 12	17.66, 17.66, 10.59	23.66, 23.66, 22.59	23.66	3	As shown in Figure 10c.
M = 4	O = 1, ts = 6	6, 12, 18	21.86, 14.44, 9.17	27.86, 26.44, 27.17	27.86	3	As shown in Figure 11a.
	O = 2, ts = 6	6, 6, 12, 12	16.56, 16.43, 9.47, 10.10	22.56, 22.43, 21.47, 22.10	22.56	4	As shown in Figure 11b.
M = 5	O = 2, ts = 6	6, 6, 12, 12	16.49, 16.56, 10.10, 10.59	22.49, 22.56, 22.10, 22.59	22.59	4	As shown in Figure 12.

Figure 9. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 2: (a) The optimal path of the UAV when O = 1, t_s = 0; (b) the optimal path of the UAV when O = 1, t_s = 6; (c) the optimal path of the UAV when O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 9. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 2: (a) The optimal path of the UAV when O = 1, t_s = 0; (b) the optimal path of the UAV when O = 1, t_s = 6; (c) the optimal path of the UAV when O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 10. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 3: (a) The optimal path of the UAV when O = 1, t_s = 6; (b) the optimal path of the UAV when O = 1, t_s = 10; (c) the optimal path of the UAV when O = 1, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 10. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 3: (a) The optimal path of the UAV when O = 1, t_s = 6; (b) the optimal path of the UAV when O = 1, t_s = 10; (c) the optimal path of the UAV when O = 1, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 11. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 4: (a) The optimal path of the UAV when O = 1, t_s = 6; (b) the optimal path of the UAV when O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 11. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 4: (a) The optimal path of the UAV when O = 1, t_s = 6; (b) the optimal path of the UAV when O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 12. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 5, O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

Figure 12. Considering take-off operation time experiment of single objective optimization without considering $T_k$ when M = 5, O = 2, t_s = 6. Among them, the different colors lines present different optimal path for different UAVs.

It can be seen that when the take-off operation time of the UAV is considered, the operation time of each UAV increases, and the total operation completion time also increases. Therefore, it affects the optimal number of the UAVs. That is, when the take-off operation time is enough, the number of the UAVs needed is reduced. It is shown that if the number of operators increases, the number of available UAVs will increase and the total operation time will reduce. However, it is not the case that the more UAVs are used, the more reasonable the results.

4. Single Objective Optimization of Coverage Path with Obstacles

The corresponding optimization model needs to be re-established when the area to be covered is concave polygon or contains obstacles inside. With the foundation of the optimization model of the coverage path of the area to be covered by the convex polygon in the previous chapter, this section takes the shortest total operation completion time as the optimization objective and establishes the optimization model of the coverage path of the concave polygon or the area with obstacles to be covered. Another approach to planning in nonconvex regions is to divide the regions into some convex part, and the Boustrophedon Cellular Decomposition method is one of those. Therefore, several groups of simulation experiments were conducted to compare the completion time of the two CPP methods when the coverage area was concave polygon or contained obstacles inside. The results show that the proposed optimization model significantly reduces the completion time and greatly simplifies the CPP problem.

4.1. Determination of Converage Flight Direction

First, a set of coordinate points was defined to indicate the boundary of the area to be covered, assuming that the set of coordinate points was $B^{\prime }=\{{b_k}^{\prime }=(x_k,y_k),k=1,2,\cdots ,n\}$. It could be further obtained that the edge set of the area to be covered was $L=\{l_k=[(y-y_k)(x_k-x_{k+1})$$=(x-x_k)(y_k-y_{k+1}),x\in (x_k,x_{k+1}),k=1,2,\cdots ,n]\}$. When the area to be covered contained concave edges or internal obstacles, the edge set L was divided into a concave edge or obstacle boundary set $L_1=\{l_k,k=1,2,\cdots ,h\}$ and a convex edge or non-obstacle boundary set $L_2=\{l_k,k=1,2,\cdots ,n-h\}$, as shown in Figure 13.

When the area to be covered contains concave edges or internal obstacles, the UAV needs to turn. Therefore, the number of turns is not only related to the width of the area to be covered, but also to the position and shape of the concave edges or obstacles. In this case, the selection principle of the optimal coverage direction was as follows: first, a group of lines with the width of coverage were intersected with the area to be covered at different angles when the coverage width was determined. Then, the angle range was obtained with the minimum number of intersection nodes. Secondly, the vertical direction with the minimum width direction of the area to be covered was selected from this set of angle ranges as the optimal coverage flight direction of the UAV.

4.2. Optimal Model of Coverage Path

The values of $N_l$ and $d_l$ were calculated after the direction of coverage flight was determined. A set of parallel lines with a distance of $d_l$ was used to intersect the boundary line of the area to be covered. As shown in Figure 14, a set of intersection points was obtained, which can be expressed as $P^{\prime }=\{{p_i}^{\prime }=(x_i,y_i),i=1,\cdots ,N-1\}$.

The new node set $V^{\prime }$ in the new undirected graph $G^{\prime }=(V^{\prime },E^{\prime })$ was composed of all the intersection points in $P^{\prime }$ and the take-off nodes of the UAV. Similarly, each node in the set was numbered, and the number of all nodes is $N$. In this case, the number of the UAV take-off node is 1, the node numbers related to the first coverage row are 2 and 3, the node numbers related to the second coverage row are 4 and 5, and so on. The edge set $E^{\prime }$ of the undirected graph $G^{\prime }=(V^{\prime },E^{\prime })$ consists of line segments of $N$ nodes connected to each other, and then a complete undirected graph is formed.

Similarly, the undirected graph $G^{\prime }$ can be represented by the cost matrix C of $N\times N$, where element $C_{ij}$ represents the Euclidean distance between node $i$ and node $j$. The edge set in the undirected graph $G^{\prime }$ is:

$$E^{\prime }=\{l_{ij}=[(y-y_i)(x_i-x_j)=(x-x_i)(y_i-y_j),x\in (x_i,x_j)],i,j=1,2,\cdots ,N\}$$

(42)

If the area with concave edges and obstacles cannot be traversed, such as mountains, etc., then the path which intersects the concave edge and the obstacle boundary are not selected. To facilitate the subsequent definition of the constraint equation, the elements in the cost matrix that represent the line segment that intersects the concave edge are set to zero, as shown in Equation (43).

$$\begin{array}{c}\mathrm{if} \exists \forall (x_a,y_a), x_a\in (x_k,x_{k+1})\cup (x_i,x_j)\\ (y_a-y_k)(x_k-x_{k+1})=(x_a-x_k)(y_k-y_{k+1}),x_a\in (x_k,x_{k+1}),k=1,\cdots ,n\\ (y_a-y_i)(x_i-x_j)=(x_a-x_i)(y_i-y_j),x_a\in (x_i,x_j),i,j=1,2,\cdots ,N\\ C_{ij}=0\end{array}$$

(43)

After removing the line segment that intersects the concave edge and the boundary of the obstacle area in the undirected graph, the undirected graph representing the optimized model for the concave polygon or the area containing obstacles can be updated, as shown in Figure 15.

After determining the undirected graph for this problem, the optimization model of coverage path is also established based on undirected graph search. Then, the optimal flight path of each UAV can be solved by defining constraint equations and objective functions, as shown in Figure 16.

The objective function and constraint equation in this optimization problem are also defined. It is solved by mixed integer linear programming. The definitions of variables $C_{ij}$, $X_{ij}^k\in \{0,1\}$, $m\in \mathbb{N}$, $T_k$, and $V$ remain unchanged. Moreover, to constrain the UAV to fly according to certain rules, the definition of the objective function and constraints 1–4 are unchanged.

It is different from the optimization model for the area without obstacle. When the area to be covered is concave polygon or contains obstacles, constraint 5 does not need to be added in order to make the problem as solvable as possible, which means that the path will not be trapped at the edge of the concave edge or obstacle. As shown in Figure 17, when constraint 5 is added, the UAV can only fly along the edges $l_{35}$, $l_{37}$, and $l_{39}$ that intersect the concave edge or obstacle boundary, which will have no solution. When this constraint is not added, the UAV can fly from node 3 along the edges $l_{36}$ and $l_{38}$ to node 6 or 8, then, flying out of the concave edge or obstacle boundary area. Therefore, no solution can only occur when the number of rows to be covered is very small. However, when the number of rows to be covered is large, the problem may be solved. In this case, if the constraint 5 is added, the path length of each UAV may be uneven, which will cause the total operation completion time to become longer, especially when the number of the UAVs is large. This will be illustrated by comparative experiments in the next section.

Constraint 6: In the problem model for concave polygons, in order to prevent the UAV from flying along the concave edge, a new constraint condition should be added, as shown in Equation (44).

$${\displaystyle {\sum }_{k=1}^m{X}_{ij}^k=0, \mathrm{if} C_{ij}=0}$$

(44)

By solving the objective function and the constraint equations defined by constraints 1–4 and constraint 6, the optimal path for a given number of the UAVs to traverse the area to be modeled as undirected graph $G^{\prime }$ in the shortest time can be found. Similarly, in order to ensure that each UAV has the shortest path, it will be solved by an iterative method.

4.3. Contrast Simulation Experiment

In order to verify the superiority of the algorithm in this paper, the Boustrophedon Cellular Decomposition algorithm is employed to experiment on the same area compared with the proposed method in this paper. The cellular decomposition algorithm is a common method of the CPP for the area with obstacles or nonconvex parcels, whose strategy is decomposing a complete map into several simple small blocks [51]. To solve the problem that the cellular decomposition method cannot cover the complex shape area with too many obstacles, and for multiple drones to cooperate, Choset et al. proposed the Boustrophedon Cellular Decomposition algorithm [22]. Based on this improved algorithm, the intersection of the sweeping involute and the obstacle is marked. Among them, the swept area is marked as the “old” area, and the new area in front is marked as the “new” area.

This method will divide the area in multiple directions, and a certain decomposition direction formula is given in Equation (45).

$$w={\displaystyle \sum _{i=1}^m{y}_{\max ,i}}-y_{\min ,i}$$

(45)

where $m$ represents the number of monotone units, $y_{\max ,i}$ represents the y coordinates of the uppermost vertex, and $y_{\min ,j}$ represents the y coordinates of the lowest vertex in the monotone polygon unit. After the completion of the decomposition algorithm, the rotating calipers algorithm is used to determine the coverage path planning direction, and then the comparative path simulation diagram and the corresponding completion time are given to compared with the proposed method in this paper.

In the contrast experiment, we mainly discuss the influence of the shape of the area to be covered on the path planning result of multiple UAVs, so the influence of take-off operation time is ignored. As in Section 3, the flying speed of the UAV was $V=10.7784 \mathrm{m}/\mathrm{s}$ and the coverage width was $L=130 \mathrm{m}$. It was assumed that the number of available UAVs was different and then different shapes of areas were selected to be covered. The main parameters for evaluating simulation results are listed in

Table 7. Parameters for evaluating simulation results.

Parameters	Symbol	Unit
Optimal coverage direction	$\delta$	rad
Modified coverage width	$d_l$	m
Operation time	$T_k$	min
Total operation completion time	$t$	min
Optional number of the UAVs	$m$

.

4.3.1. Contrast Simulation Experiment for Concave Polygonal Area

In the first group of experiment, the coordinates of the boundary in the operation area to be covered were (−2555, −986), (−2770, −99), (341, −49), (345, −961), (−833, −1073), (−1754, −683), and (2074, −895), the UAV take-off node was (−3335, −289). The map is shown in Figure 18.

When $m=2,3$, the contract results and the optimized path of each UAV are shown in Table 8.

Table 8. Results of the first contrast experiment for concave polygonal area.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.55	125.96	2	20.63, 21.72	21.72	As shown in Figure 19a.
			3	17.03, 16.29, 16.89	17.03	As shown in Figure 19c.
Boustrophedon Cellular Decomposition			2	22.75, 23.46	23.46	As shown in Figure 19b.
			3	18.21, 15.26,17.58	18.21	As shown in Figure 19d.

Table 8. Results of the first contrast experiment for concave polygonal area.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.55	125.96	2	20.63, 21.72	21.72	As shown in Figure 19a.
			3	17.03, 16.29, 16.89	17.03	As shown in Figure 19c.
Boustrophedon Cellular Decomposition			2	22.75, 23.46	23.46	As shown in Figure 19b.
			3	18.21, 15.26,17.58	18.21	As shown in Figure 19d.

Figure 19. The optimized path of the first contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Figure 19. The optimized path of the first contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Through the first contrast experiment, the completion time of the proposed method in this paper was 21.72 (when m = 2), and 17.03 (when m = 3). Meanwhile, the completion time of the Boustrophedon Cellular Decomposition method was 23.46 (when m = 2) and 18.21 (when m = 3). It is obviously seen that the proposed method in this paper has a shorter completion time than the Boustrophedon Cellular Decomposition method for a concave polygonal area.

In the second set of experiment, the coordinates of the boundary in the operation area to be covered were (−2707, −1152), (−2702, −38), (440, 119), (553, −932), (−279, −1111), (−969, −1140), (−986, −791), (−1592, −433), (−2071, −475), and (−2058, −1152), the UAV take-off node was (−3334, −47). The map is shown in Figure 20.

When $m=2,3$, the contract results and the optimized path of each UAV are shown in Table 9.

Table 9. Results of the second contrast experiment for concave polygonal area.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.52	119.83	2	26.74, 26.64	26.74	As shown in Figure 21a.
			3	22.91, 22.38, 22.17	22.913	As shown in Figure 21c.
Boustrophedon Cellular Decomposition			2	29.78, 31.28	31.28	As shown in Figure 21b.
			3	24.21, 23.54, 25.87	25.87	As shown in Figure 21d.

Table 9. Results of the second contrast experiment for concave polygonal area.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.52	119.83	2	26.74, 26.64	26.74	As shown in Figure 21a.
			3	22.91, 22.38, 22.17	22.913	As shown in Figure 21c.
Boustrophedon Cellular Decomposition			2	29.78, 31.28	31.28	As shown in Figure 21b.
			3	24.21, 23.54, 25.87	25.87	As shown in Figure 21d.

Figure 21. The optimized path of the second contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3.

Figure 21. The optimized path of the second contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3.

Through the second contrast experiment, the completion time of the proposed method in this paper were 26.74 (when m = 2) and 22.913 (when m = 3). Meanwhile, the completion time of the Boustrophedon Cellular Decomposition method were 31.28 (when m = 2) and 25.87 (when m = 3). It is obviously seen that the proposed method in this paper has a shorter completion time than the Boustrophedon Cellular Decomposition method for a concave polygonal area.

4.3.2. Contrast Simulation Experiment for Area with Obstacles

In the first group of experiment, the coordinates of the outer boundary in the work area to be covered were defined as (−2780, −1495), (−2780, −305), (−550, −305), and (−550, −1495), the coordinates of the boundary of the obstacle in the area to be covered were defined as (−1118, −642), (−1260, −626), (−1352, −497), (−1304, −364), (−961, −299), and (−921, −396), the UAV take-off node was (−3424, −300). The map is shown in Figure 22.

When $m=2,3$, the contract results and the optimized path of each UAV are shown in Table 10.

Table 10. Results of the first contrast experiment for area with obstacles.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.57	119	2	22.92, 22.85	22.92	As shown in Figure 23a.
			3	16.43, 16.54, 16.96	16.96	As shown in Figure 23c.
Boustrophedon Cellular Decomposition			2	23.48, 27.56	27.56	As shown in Figure 23b.
			3	17.25, 19.32, 20.52	20.52	As shown in Figure 23d.

Table 10. Results of the first contrast experiment for area with obstacles.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.57	119	2	22.92, 22.85	22.92	As shown in Figure 23a.
			3	16.43, 16.54, 16.96	16.96	As shown in Figure 23c.
Boustrophedon Cellular Decomposition			2	23.48, 27.56	27.56	As shown in Figure 23b.
			3	17.25, 19.32, 20.52	20.52	As shown in Figure 23d.

Figure 23. The optimized path of the first contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Figure 23. The optimized path of the first contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Through the first contrast experiment, the completion time of the proposed method in this paper were 22.92 (when m = 2) and 16.96 (when m = 3). Meanwhile, the completion time of the Boustrophedon Cellular Decomposition method were 27.56 (when m = 2) and 20.52 (when m = 3). It is obviously seen that the proposed method in this paper has the shorter completion time than the Boustrophedon Cellular Decomposition method for an area with obstacles.

In the second group of experiment, the coordinates of the outer boundary of the work area to be covered were defined as (−2204, −170), (−533, −1227), (92, −453), (11, −283), (72, 76), and (−1691, 390), the coordinates of the boundary of the obstacle in the area to be covered were defined as (−2283, −958), (−2210, −525), (−1838, −482), and (−1799, −910), the UAV take-off node was (641, −5). The map is shown in Figure 24.

When $m=2,3$, the contract results and the optimized path of each UAV are shown in Table 11.

Table 11. Results of the second contrast experiment for area with obstacles.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.75	124.38	2	18.56, 18.35	18.56	As shown in Figure 25a.
			3	12.86, 13.62, 13.01	13.62	As shown in Figure 25c.
Boustrophedon Cellular Decomposition			2	21.94, 19.56	21.94	As shown in Figure 25b.
			3	13.75, 15.57, 11.08	15.57	As shown in Figure 25d.

Table 11. Results of the second contrast experiment for area with obstacles.

Method	δ (rad)	dl (m)	m	Tk (min)	t (min)	Optimized Path
Proposed method in this paper	1.75	124.38	2	18.56, 18.35	18.56	As shown in Figure 25a.
			3	12.86, 13.62, 13.01	13.62	As shown in Figure 25c.
Boustrophedon Cellular Decomposition			2	21.94, 19.56	21.94	As shown in Figure 25b.
			3	13.75, 15.57, 11.08	15.57	As shown in Figure 25d.

Figure 25. The optimized path of the second contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Figure 25. The optimized path of the second contrast experiment: (a) The optimized path of the proposed method in this paper when m = 2; (b) the optimized path of the Boustrophedon Cellular Decomposition method when m = 2; (c) the optimized path of the proposed method in this paper when m = 3; (d) the optimized path of the Boustrophedon Cellular Decomposition method when m = 3. Among them, the different colors lines present different optimized path for different UAVs.

Through the second contrast experiment, the completion time of the proposed method in this paper were 18.56 (when m = 2) and 13.62 (when m = 3). Meanwhile, the completion time of the Boustrophedon Cellular Decomposition method were 21.94 (when m = 2) and 15.57 (when m = 3). It is obviously seen that the proposed method in this paper has the shorter completion time than the Boustrophedon Cellular Decomposition method for an area with obstacles.

4.3.3. Comparative Experiment of Adding Constraint 5

When constraint 5 is added, the UAV will not be allowed to fly along nodes that are not parallel to the coverage row. Here, we selected the test conditions in the fourth group of experiment in this section. The comparison is discussed when adding constraint 5 and not adding constraint 5 when the number of available UAVs is 3. Figure 26 shows the experimental result of adding constraint 5.

Table 12. Results of comparative experiment of adding constraint 5.

Constraint 5	Tk (min)	t (min)
Not added	16.43, 16.54	16.96
Added	21.81, 24.05, 8.96	24.05

shows the operation time and total operation completion time of each UAV in two cases.

Through the comparison between Figure 23c and Figure 26, it can be seen that there are three paths which are not parallel to the row to be covered in the UAV path without the constraint condition. However, in the experiment with constraints, the UAV flies in the direction parallel to the coverage row. According to the experimental results in Table 12, when this constraint is not added, the flight operation time of each UAV is relatively average, and the total operation completion time is relatively small. However, when this constraint is added, in order to make the problem solvable, the path length of each UAV is not average, and the total operation completion time is relatively long.

5. Bi-Objective Optimization of Coverage Path Considering Energy Flow Efficiency in Areas Containing Obstacles

The UAV’s attitude change will affect its energy production and energy consumption when turning. Then, it will affect the energy flow efficiency. Therefore, in this section, for the area with obstacles, the attitude change of the UAV is considered during flight. The shortest total operation time and the maximum total energy flow efficiency are taken as the optimization objectives, so as to determine the best coverage flight direction and the optimal path of each UAV.

5.1. Turning Path Refinement Strategy

The attitude angle of the UAV remains unchanged when it is in the horizontal flight state. The roll angle needs to be adjusted to enable the UAV to achieve a turning flight according to a certain turning radius while in the turning state. It is necessary to determine the turning strategy of the UAV in different situations when the influence of the UAV’s attitude change is considered during turning.

There are two situations that need to be considered in this problem: the turning path to switch between the two covered rows of the UAV and the turning path of flying from the take-off node to the first node or before returning from the last node to the take-off node. The principle of the shortest total path length is used to determine the best turning path when the turning radius is determined. First, the turning path is refined when the UAV changes between two covered lines.

The path before the UAV turns and refines is assumed as shown in Figure 27. The distance between the two nodes is $d_{AB}$. The distance between the two covered rows is $d$, the turning radius of the UAV is $R$, and the direction of the arrow in the figure indicates its flight direction.

It can be seen from the reference [52] that the optimal turning strategy of the UAV needs to be classified and discussed according to its turning radius $R$, the distance $d$ between the two lines, and the angle $\beta$ between the polyline AB at the turn and the covered line $L_1$. The optimum in each case path is shown in Figure 28.

Then, the refinement strategy is considered of the turning path when the UAV flies from the take-off node to the first node or returns from the last node to the take-off node. Two cases where the angle is acute or obtuse between the take-off node and the first/last node and the line to be covered are considered. In order to ensure the coverage rate of the coverage operation and make the path as short as possible, the best flight path is designed as shown in Figure 29.

The minimum turning radius $R_{\min }$ of the UAV can be calculated, as shown in Equation (46).

$$R_{\min }=\frac{v^2}{g\tan {\varphi }_{\max }}$$

(46)

It can be seen from the reference [53] that the safe range of the fixed-wing UAV roll angle is $-\pi /4~\pi /4 \left(\mathrm{rad}\right)$, that is, the maximum roll angle is ${\varphi }_{\max }=\pi /4 \left(\mathrm{rad}\right)$. Therefore, $R_{\min }=7.15 \mathrm{m}$ from Equation (46).

5.2. Coverage Path Optimization Model

The roll angle will change when the UAV turns during the flight, and this angle will affect the energy flow efficiency of the solar powered UAV. Different turn path proportions affect the energy flow efficiency when the energy flow efficiency is also used as the optimization target. Different coverage flight directions will affect the turn path proportions. In this section, we take the energy flow efficiency and the shortest total operation completion time as the bi-objective optimization to determine the best coverage flight direction and the optimal path for each UAV. The bi-objective optimization equation is obtained as Equation (47).

$$\min (y)=k\frac{t_x}{t_{\min }}+(1-k)\frac{{\eta }_{flow\max }}{{\eta }_{flowx}}$$

(47)

where $k$ and $1-k$ represent the scale factor of the total time and the total energy flow in the optimization equation.

Due to the different attitudes of the UAV in level flight and turning, there will be different energy flow efficiency. It is assumed that the energy production power and energy consumption power are $P_{in1}$ and $P_{out1}$ of the UAV in the turning state. The energy production and energy consumption power are $P_{in2}$ and $P_{out2}$ of the UAV in the horizontal flight state. The time of the UAV used for turning and non-turning are $t_1$ and $t_2$, then the energy flow efficiency ${\eta }_{flow}$ when considering attitude change of the UAV can be calculated as shown in Equation (48).

$${\eta }_{flow}=\frac{P_{out1}\times t_1+P_{out2}\times t_2}{P_{in1}\times t_1+P_{in2}\times t_2}$$

(48)

In the following simulation experiment, $10°$ is selected as a unit. Then, the coverage flight direction is divided into 18 cases, and the operation time and energy flow efficiency of each case is calculated. Different values of $k$ are set in each case. According to Equation (48), the optimal coverage flight plan is obtained in different values.

5.3. Simulation Experiment

The simulation experiments in each of the following cases were repeated 10 times, and the results listed are the average of each result. By comparing the total energy flow efficiency in different flight directions, the best coverage flight direction was determined. The most suitable weighting factor $k$ was obtained by comparing the effects of different $k$ values.

It was assumed that the width of the coverage operation was $L=130 \mathrm{m}$ and the speed of the UAV was set to $V=10.7784 \mathrm{m}/\mathrm{s}$. According to Equation (42), the minimum turning radius can be calculated as $R_{\min }=11.85 \mathrm{m}$.

Moreover, the turning radius of the UAV was set to $R=60 \mathrm{m}$, and the rolling angle was $\varphi =11.2°$ based on Equation (42) when the UAV turns. Then, the energy production power and energy consumption power in the two states could be calculated as follows:

$$P_{in1}=83.5183 \mathrm{W}, P_{in2}=72.4627 \mathrm{W}$$

(49)

$$P_{out1}=85.0733 \mathrm{W}, P_{out2}=68.3057 \mathrm{W}$$

(50)

It was assumed that the location of the coverage operation was $(125°\mathrm{E},50°\mathrm{N})$ for 10:00 am on 1 May 2022. The boundary of the area to be covered and the take-off node conditions were the same as the second group of experimental conditions of the coverage path with obstacles. The number of available UAVs was set to $m=2$ and the take-off operation time was ignored. The best UAV path when covering different flight directions is shown in Figure 30.

According to the turning path refinement strategy, the flight path can be obtained after the turning path refinement in each case, as shown in Figure 31.

According to the results of simulation experiments, the operation time $t_k$ and total operation time $t$ of a single UAV can be obtained when covering different flight directions. Then, the total turning path time $t_1$ and non-turning path time $t_2$ can be obtained according to the refined path. Therefore, the total energy flow efficiency ${\eta }_{flow}$ can be obtained in each case. The experimental results are listed in

Table 13. Experimental results with different direction of coverage. Among them, the shortest total operation time and the biggest total energy flow efficiency are shown in bold and highlight.

$\mathit{\delta }(°)$	$\mathit{t}(\mathbf{min})$	${\mathit{\eta }}_{\mathit{f}\mathit{l}\mathit{o}\mathit{w}}(\mathit{\%})$	$\mathit{\delta }(°)$	$\mathit{t}(\mathbf{min})$	${\mathit{\eta }}_{\mathit{f}\mathit{l}\mathit{o}\mathit{w}}(\mathit{\%})$
0	19.87	95.24%	90	23.48	95.57%
10	19.36	95.26%	100	23.43	95.55%
20	19.42	95.46%	110	21.86	95.84%
30	19.84	95.80%	120	22.19	95.74%
40	20.73	95.32%	130	21.00	95.60%
50	22.25	95.54%	140	21.38	95.36%
60	23.41	95.71%	150	21.04	95.51%
70	22.56	95.35%	160	21.41	95.46%
80	24.75	95.35%	170	21.02	95.48%

.

It can be seen from Table 13 that the total energy flow efficiency (${\eta }_{flow\max }=95.84\%$) is the largest when the flight direction $\delta =110°$, and the total operating time ($t_{\min }=19.36 \min$) is the smallest when the flight direction $\delta =10°$. From Equation (44), we can get the optimization objective $y$ when the value $k$ is different under different coverage flight directions.

In the case of the weighting factor $k=0.1$ and the flight direction $\delta =30°$, the total objective value $y$ of the bi-objective optimization is the smallest, which is $\min (y)=1.0029$. In the cases of the weighting factor $k=0.2, 0.3, 0.4$ and the flight direction $\delta =20°$, the total objective value $y$ of the bi-objective optimization is the smallest, which is $\min (y)=1.0038, 1.0037, 1.0036$. In the cases of the weighting factor $k_{}=0.5, 0.6, 0.7, 0.8, 0.9$, and the coverage flight direction $\delta =30°$, the total objective value $y$ of the bi-objective optimization is the smallest, which is $\min (y)=1.0030, 1.0024, 1.0018, 1.0012, 1.0006$.

6. Conclusions

In order to solve the problem of short endurance of traditional electric UAV during coverage operations, this paper proposed to use solar powered UAVs for coverage operations. We proposed an optimization model of coverage path based on undirected graph search. The main works are summarized as follows:

(1) We established the energy production model and energy consumption model of solar powered UAV during the covering operations and defined the evaluation index of energy flow efficiency to evaluate the utilization rate of solar power absorbed by solar powered UAV during operations.

(2) For the area with and without obstacles, the direction of coverage flight is determined based on the principle of the minimum number of turns, and the optimization objective is based on the shortest total operating time. An optimization model of coverage path based on undirected graph search is established to solve the optimal flight path of the UAVs. Compared with the Boustrophedon Cellular Decomposition method, the completion time of proposed method in this paper is obviously improved.

(3) We considered the impact of the attitude change when turning of the UAVs on the energy flow efficiency. Moreover, the shortest total operation time and the highest total energy flow efficiency were used as the optimization index. Then, a bi-objective optimization equation was established to solve the optimal coverage flight direction and the optimal flight path of each solar powered UAV.

In this paper, the coverage path optimization model proposed has great applicability, and the optimal objective of the UAV coverage path planning can be obtained under different complex areas to be covered. Therefore, the method is feasible.