Cost-Minimizing System Design for Surveillance of Large, Inaccessible Agricultural Areas Using Drones of Limited Range

: Drones are used increasingly for agricultural surveillance. The limited ﬂight range of drones poses a problem for surveillance of large, inaccessible areas. One possible solution is to place autonomous, solar-powered charging stations within the area of interest, where the drone can recharge during its mission. This paper designs and implements a software system for planning low-cost drone coverage of large areas. The software produces a feasible, cost-minimizing charging station placement, as well as a drone path speciﬁcation. Multiple optimizations are required, which are formulated as integer linear programs. In extensive simulations, the resulting drone paths achieved 70–90 percent of theoretical optimal performance in terms of minimizing mission time for a given number of charging stations, for a variety of ﬁeld conﬁgurations.


Introduction
Agricultural surveillance is an important part of ensuring optimal crop yields, by detecting problems in the growth process or during early stages so that proactive measures can be taken. Traditionally, surveillance of agricultural areas was conducted by human caretakers who personally surveyed the fields. This approach is error-prone, and limited to small size areas. To obtain better surveillance, some researchers have proposed video surveillance systems with IP cameras in a wireless network for farms or forests [1][2][3][4][5]. Although this technology can provide improved surveillance, several drawbacks remain, including the cost of the system for large area surveillance and the risk of blind spots. In addition, positioning of cameras may be difficult, and image quality may not be sufficient. To overcome those challenges, drones have been proposed as a solution to agricultural surveillance [6][7][8][9]. In the basic scenario, a drone flies over an agricultural area farm and collects image data that will be analyzed later for various applications, including weed detection and mapping [10][11][12][13]; growth monitoring [14][15][16]; crop health monitoring [17][18][19]; irrigation management [20,21]; yield prediction [22][23][24]; crop lodging detection [25][26][27]; and vineyard monitoring [28,29]. Drones may carry cameras sensors that detect various types of radiation, including visible light (RGB), thermal, multispectral, or hyperspectral sensors. Multispectral and hyperspectral sensors yield much more information than RGB but at much higher cost. Both fixed-wind and rotary-wing drones are used, with multirotor drones being most commonly represented in the literature [30]. Until now, most deployments have been located in developed countries.
Agricultural surveillance missions with drones are generally done in automatic mode, meaning that the path the drone will follow to cover the given farm is determined beforehand. This problem is known as the Coverage Path Planning (CPP) [31]. The earliest works on CPP in agriculture were focused on ground-based robots [31,32]. In some works, field topography as well as field geometry have been taken into account [33,34]. CPP usually aims to determine an optimal path that enables a drone to complete a mission in a minimum time. Obviously, the mission completion time depends on the size and the form of the area to cover. The area studied is usually a planar polygon and is decomposed based on the Ground Sampling Distance (GSD) into smaller polygons [35]. The robots' paths are typically based on a back-and-forth pattern referred to as "boustrophedon" or "lawn-mower".
More recently, CPP for unmanned aerial vehicles (UAV) has also been considered [36][37][38][39]. A related active area of research is cooperative patrolling with multiple agents [40,41], in which a variety of techniques have been applied, including reinforcement learning and graph theoretic approaches, among others.
An important challenge in CPP for drones remains the limited capacity of the battery that restricts the drone to short-range missions. Several efforts have been made in the development of energy aware path planning algorithms to extend the size of the area a drone can cover. A common idea in most of the current algorithms for CPP using drones relies on the assumption that the drone requires more energy when making turns than going straight [42,43].
Although energy-aware algorithms allow the efficient use of a drone's battery, the coverage remains limited and the drone requires more than one trip for larger areas. In a simple surveillance mission, DJI Phantom 3 professional requires around ten flights to cover an area of about four square kilometers [44]. To conduct surveillance mission in large area, charging stations (CSs) have recently been introduced to assist drones [45,46]. A CS is a platform on which a drone can land to recharge its battery when it runs out of energy. The charging can be done wirelessly [47], and the station can even move [48]. However, in most scenarios, stations are deployed at predefined locations (as in [49]), where only one drone is used with the aim of minimizing the completion time of the mission by deciding depending on the state of charge when the drone will land to recharge its battery. An alternative approach was used in [50], where a single CS installed at the center of the area was used to service a fleet of drones. In this case, the area of interest is limited to a circle with diameter equal to the traveling distance of a drone, in order to ensure that a drone can make a round trip to the CS.
Planning efficient large scale surveillance missions using drones of limited range will require one or more CSs, optimally placed in the area of interest. Relatively few works attempted to solve the placement issue of CSs in a given area of interest. Hong et al. [51] give a procedure for locating a set of CSs with a view towards maximally satisfying customer demands. In their formulation, the number of stations is known beforehand. This assumption is relaxed in [52,53], where the authors attempt to optimize the overall mission performance in terms of energy, communication, and safety.
Though CSs are necessary for large scale automatic surveillance, their use incurs high costs. For instance, the CS M2 Combo costs about USD 4000 (Price as of 23 September 2020, from https: //www.heishatech.com/product/c500-m2-combo/), and is designed for the DJI Mavic 2 Pro drone whose price is 60% lower (The price for a Mavic 2 Pro equipped with an embedded Hasselblad camera with a 20MP 1" CMOS sensoron www.amazon.com was USD 1599, as of 23 September 2020.). Thus decreasing the number of CSs even by just one unit significantly reduces the cost of the architecture.
This paper addresses the coverage path planning of a single drone using multiple CSs for the surveillance of large agricultural areas. Since there is usually no grid power available in remote agricultural areas, we assume the CSs are being powered by solar energy. This work tackles the path planning from an economic perspective in which the number of CSs should be minimized first, and then the completion time of the mission. Our algorithm can be applied to general surveillance regions, and drones with arbitrary flight and range characteristics: however, to demonstrate feasibility we employ a scenario that is appropriate for applications in developing countries, in which low-cost drones are used and the area under surveillance contains irregular terrain such that only a limited number of locations are suitable for charging station installation.
The rest of the paper is organized as follows: Section 2 presents the material and methods and presents the assumptions as well as mathematical and computational details of the algorithm used to solve the problem. Section 3 presents the results, followed by discussion in Section 4. This paper ends with a conclusion and several ideas for future works.

Model Assumptions
Our CPP algorithm generates a drone mission path under the following assumptions: • It is assumed that minimizing the cost of the system is considered to be more important than all other factors. This will usually be the case in resource-poor situations. However, if the user wants to investigate alternative system configurations that give shorter mission times, (s)he may reduce the CS radius parameter (explained below), which will lead to higher-cost systems with reduced mission times.

•
The area that is to be surveilled (henceforth referred to as "the field") is polygonal and may thus be completely specified by the vertices of the bounding polygon.

•
The drone is assumed to have a single flight mode, in which the drone's camera runs continuously and video information is either stored onboard or is continuously transmitted back as the drone flies.

•
The drone's height, velocity, and power consumption are assumed to be constant throughout the drone's mission (except when recharging).

•
A start location is specified by the user, which is a point (presumably somewhere on the field's edge) from which the drone is originally launched to conduct its mission. When launched from the start location, the drone is fully charged.

•
Possible CS locations are scattered throughout the region and are specified to the algorithm by the user. This assumption is included to accommodate the situation where some locations within the field are more favorable than others for CS installation. If all locations are equally favorable, the user can input a regular or quasirandom lattice of points that fills the entire region.

Mission Path Planning Algorithm Justification and Outline
In developing a cost-and time-minimizing mission path, we first take into account the limited range of drones, which imposes hard constraints on the set of charging stations. In particular, every point in the field must be reachable by a drone that originates its flight from a charging station and subsequently returns to a charging station. In our algorithm, we impose the condition that for every field point there is at least one charging station from which a round trip can be made to that particular field point. Since minimizing the number of CSs is our primary goal, we restrict possible solutions to those sets of CSs that both satisfy the constraint and have minimal size.
Given these possible CS solutions, our goal is to select a solution that minimizes mission time. It would be much too computationally intensive to compute mission times for a large number of possible solutions. Instead, we select a set of CSs that minimizes the maximum distance from field point to nearest charging station, where the maximum is taken over all points within the field. This criterion is motivated by the fact that field points that are far from the nearest charging station tend to produce coverage inefficiency, because reaching them requires flying multiple flyovers of closer field points (this point is explained in more detail in Section 2.9).
Once the CS have been decided based on the above criteria, the next step is to design an efficient path that takes full advantage of the characteristics that have been used to select the CS solution.
As noted above, field points that are reached from a far-away CS are likely to have inefficient paths. For this reason, we assign each field point to nearest charging station. This divides the field into Voronoi regions. Our strategy is to completely cover each Voronoi region before passing on to the next. It follows that we need a spanning walk that begins and ends at the starting location.
Finally, the drone's path for covering each Voronoi region must be specified. Our construction is based on the observation that the shortest covering path for a drone with infinite range is a "boustrophedon" or "lawn-mower" pattern, which provides single coverage throughout the field. Although a finite-range drone cannot execute a boustrophedon, nonetheless a similar pattern of nonoverlapping parallel segments can be constructed to obtain a solution that closely resembles a boustrophedon path, and is nearly as efficient. Our algorithm has a recursive structure, based on the self-similarity of boustrophedon paths within polygons.
The steps in the drone path planning algorithm may be summarized as follows: 1.
Initialization of parameters for field specifications, potential CS locations and drone range specifications; 2.
Create binary grid representation of the field; 3.
Use iterative linear programming to find a minimal set of CSs that minimizes the maximum distance from field points to the nearest CS; 4.
Decompose field into Voronoi cells based on selected CS locations; 5.
Compute a closed spanning walk visiting all CSs using a modified traveling salesman algorithm; 6.
Perform a triangular decomposition of each Voronoi cell; 7.
Construct drone paths within each triangular region.
In the following subsections, we will describe in more detail the steps listed above. The entire algorithm was coded in Python, using several pre-existing libraries (as described below). The code is available on Github (https://github.com/LuisEVT/Optimal_Drone_Field_Coverage).

Initialization
In order to use the system, the following parameters must be supplied by the user:

•
Field of view (FOV) radius: The camera's field of view is assumed to be circular with fixed radius. As the drone flies, it takes continuous video pictures, thus covering a strip whose width is twice the FOV radius.

•
Drone range: This is the maximum distance that a drone can reliably travel on a full charge. • CS radius: The CS radius is the maximum distance between any point in the field and the nearest CS. In other words, a CS located at point C services a region that lies within a circle whose radius is equal to the CS radius with center at C. The CS radius must necessarily be less than or equal to half the drone's range. However, the user has the option of choosing a smaller CS radius. A smaller CS radius will tend to produce solutions that have more CS and shorter surveillance mission times. • Mesh step: For the sake of numerical computations, the coverage region is represented as a rectangular grid, whose mesh size is given by the mesh step. • Field vertices: The user specifies the polygonal field area by specifying the field vertices in consecutive order. In our simulations we employed three basic shapes (rectangle, square, and hexagon, as shown in Figure 1), and each shape had three different sizes corresponding to total areas 25, 50, and 100 km 2 .

•
Start location: This is the drone launching point, as explained in Section 2.1. • CS vector length: according to our assumptions, there will be a number of favorable locations within the field for placing the CSs. In our simulation, these favorable locations are chosen randomly. The CS vector length gives the number of randomly-generated favorable locations within the field. In our simulation, we chose the CS vector length so that there was an average of 1 favorable location per square kilometer. In practice, the user would determine favorable locations by observation and supply these locations to the algorithm.  Table 1 summarizes the values of these parameters that were used in the simulations in this paper, along with symbols used to represent them in our subsequent description of the algorithm.

Coverage Field Specification
Next, x and y mesh grid matrices are generated to discretize a rectangular region containing the field. The x values run from x min ≡ min n x n to x max ≡ max n x n and the y values run from y min ≡ min n y n to y max ≡ max n y n , both with the same mesh step s. These x and y values are used to create a list of (x, y) points within the rectangle, of which some are in the field and some are not, as shown in Figure 2a. The points in this list that lie in the field are extracted by inputting the field vertices into the Path class and using the contains_points() function from the matplotlib.path package in Python. The list of points that lie within the field is further reduced by removing points that are within a distance R of the start location (ι x , ι y ): these points are removed because they can be reached from the start location and will not need another charging station to cover them. The remaining list of (x, y) points form a grid that discretizes the field area that needs to be covered using charging stations, as shown in Figure 2b. Figure 2. (a) Points in the discretized rectangular (x,y) meshgrid that includes field points, for a pentagonal field.The outline of the field is shown as a dashed line. (b) Remaining grid points within the field after applying contain_points() and removing points within R of the start location (0,0).

Charging Station Location Selection
The primary optimization parameter is the number of CSs, because minimizing the number of CSs minimizes the cost of the system. The optimal solution must also ensure that each grid point within the field is within a distance R of the nearest CS. We may formulate the minimization problem as a linear program as follows: where: • K and M are the number of candidate CS and the number of points in the field discretization, respectively.
A is a M × K incidence matrix. The entry A mk is equal to 1 if the mth grid point in the field discretization is within a distance R of CS k, and 0 otherwise. • z is a K × 1 solution vector. Interpreted as a logical index vector: z k = 1 if CS k is in the solution, and z k = 0 otherwise.
In (1), the scalar quantity 1 T z represents the total number of CS in the solution, which is the quantity to be minimized. The condition A z ≥ 1 represents M inequalities that guarantee that each of the M grid points within the field lies within a distance R of at least one CS. Our code uses the glpk.ilp integer linear program solver from the cvxopt convex optimization package. This solver implements the branch-and-bound algorithm. As mentioned in Section 2.2, the value of R is successively reduced from its initial value until either the number of CSs in the solution increases or no solution exists.

Creation of Voronoi Regions
Once we have selected a particular set of CS locations, the next step is to design the drone's path. The first step in this process is to identify the (x, y) locations belonging to the jth charging station: that is, the points in the field which are reached following a recharge at station j.
Any finite set of J points P j , j = 1 . . . J in the (x, y) plane determines a set of J Voronoi cells, where each Voronoi cell is a convex polygon that contains exactly one of the points (see Figure 3a). The Voronoi cell for point P j consists of all points in the plane that are closer to point P j than any of the J − 1 other points. Thus the Voronoi cell of each CS corresponds to the region which is closest to that CS from among all CS. If we intersect the j'th Voronoi cell with the field, the result gives the region to be covered by means of drone trips between which the drone is recharged at CS j.
The Python package spatial in the library SciPy (version 1.4.1) contains the class Voronoi which return vertices of Voronoi cells, when the point (x, y) locations are supplied as inputs. The class only returns vertices of finite regions (see Figure 3a): thus, an additional four proxy sites located outside of the field were used to ensure that all Voronoi cells of interest were finite. The four proxy sites formed a rectangle that contained the entire field, as shown in Figure 3b. Once these finite Voronoi regions were generated, the field boundaries were used to restrict the Voronoi regions so that the resulting cells were all contained within the field. For this purpose, the Python library Shapely was used to determine the polygons obtained by intersecting the polygonal Voronoi regions with the user-specified field polygon, as shown in Figure 3c.

Charging Station Spanning Walk Construction
Once the locations of the CSs have been fixed, the next step is to specify the actual drone mission path. To cover the entire region, the drone will have to pass from Voronoi region to Voronoi region and cover each visited region by repeated trips to the region's CS. In order to characterize the drone's motion from CS to CS, we recast the situation in graph theoretic terms. Each CS is identified as a vertex in a graph, and two vertices are joined by an edge if the Voronoi regions of the two corresponding CSs share a common bounding edge. This path linking the different CSs may be described in graph theoretic terms as a closed spanning walk: that is, a sequence of vertices joined by edges that contains all vertices in the graph, such that the first and last vertex is the same. Note that according to this definition, both vertices and edges can be repeated in the walk: repetition may in fact be necessary to attain the spanning property, as shown in Figure 4a. In our construction of the closed spanning walk, we shall minimize the total number of edges, because extra edges may lead to duplicate coverage.
In order to pose the closed spanning walk construction as a linear programming problem, we associate the walk with a Hamiltonian cycle in a larger graph as follows. With the given vertices, we create a complete weighted graph such that between every two vertices there is an edge. The edge between two vertices is assigned a weight equal to the number of edges in the shortest path within the original graph joining the two vertices. A closed spanning walk in the original graph that minimizes edges corresponds to a Hamiltonian cycle in the new graph which minimizes the sum of edge weights. Finding a weight-minimizing Hamiltonian cycle is equivalent to the classical traveling salesman problem. The steps of the algorithm for finding the closed spanning walk that visits each CS are summarized as follows:  In our system, for step (b) we made use of the implementation of Dijkstra's algorithm in the scipy class scipy.sparse.csgraph. The weight-minimizing Hamiltonian cycle in (d) was found using the following linear program: Minimize w T e subject to M e = 2, where e k ∈ {0, 1} where: • E and N denote the number of vertices and edges in the graph, respectively; • w is an E × 1 vector , where w j is the weight associated with the j'th edge in the complete graph; • M is a N × E vertex-edge incidence matrix. The entry m ij is equal to 1 if node i is an endpoint of edge j; • 2 is an E × 1 matrix of all 2's. • e is an E × 1 solution vector, interpreted as a logical index vector: e j = 1 if edge j is in the solution, and e j = 0 otherwise.
The scalar quantity w T e gives the sum of weights of edges in the Hamiltonian cycle in the complete graph, which is equal to the number of edges in the corresponding closed spanning walk in the original graph. The condition M e ≥ 2 represents the fact that each vertex in the cycle must have an entering and leaving edge. When a solution to (2) is computed, the corresponding tour is checked for connectivity. If the cycle is not connected, then an additional constraint is added to the linear program that prevents this solution-and this process is iterated until an optimal connected cycle is contained. The minimum connecting paths found using Dijkstra's algorithm are then used to construct the closed spanning walk on the original graph.
In the actual drone path, the drone does not fly directly from one CS to another but rather passes by way of a vertex that is shared by the coverage cells of the two CSs (see Figure 5). There are two reasons for this. First, the paths between adjacent CSs cover vertices that are far from any CS, when they do not need to be covered by long round trips to-from a single CS. Second, the straight segments of the paths between CSs lie along the edges of triangles used in the triangulations of the Voronoi regions described in the next section, so that the segments are not double-covered when the triangle interiors are covered as described below. This serves to improve the algorithm's efficiency in terms of reducing total mission time. The algorithm is configured so that no segment of the path linking the CSs is traveled twice in a complete path traversal, as may be seen in Figure 5. CSs are shown in red, while Voronoi regions are outlined in black. Note that drone does not fly directly from CS to CS but rather travels by way of common vertices between Voronoi regions as described in the text.

Drone Trajectories within Individual Voronoi Cells
To simplify the drone's coverage of Voronoi cells, the cells were further divided into triangles, as shown in Figure 6. Suppose that the vertices of the j'th Voronoi region are given by (p 1 , q 1 ), . . . , (p n j , q n j ). Then the region is divided into n j triangles, where the vertices of the i'th triangle are given by {(u j , v j ), (p i , q i ), (p i+1 , q i+1 )}, i = 1, . . . n j (in this expression (p n j +1 , q n j +1 ) is identical to (p 1 , q 1 )).
The algorithm for finding the drone path within each triangle was programmed as a subroutine. Before calling the subroutine for a given triangle, the triangle is rotated so that the longest side is parallel to the y axis and the bottom vertex of the longest side is translated to (0, 0), as shown in Figure 7. If possible, this vertex is chosen as the CS location, and if necessary, the triangle is reflected so that all points in the triangle lie in the first quadrant. The linear transformation applied to (x, y) may be expressed as: where (x 0 , y 0 ) is the vertex that is translated to (0, 0); (x 1 , y 1 ) is the other vertex of the longest side; and the ± corresponds to whether the triangle is or is not flipped (− or + respectively). The drone path algorithm used within the triangle is based on a "boustrophedon" (or "lawn mower") pattern. Some double coverage is necessary, when multiple recharges are required for the drone to cover the entire triangle. The drone path that is constructed within the triangle consists of segments parallel to one the three sides, joined (if necessary) by short connecting segments. The path is constructed segment by segment, and each time a segment is constructed the remaining region to be covered remains triangular in shape but with reduced area. The basic idea of the algorithm is as follows. Let V 0 , V 1 , V 2 be the vertices of the triangle, where V 0 is the vertex in which the CS is located, V 1 is the farthest vertex from CS, and V 2 is the remaining vertex. The algorithm will always start by traveling from V 0 to V 1 . Ideally, we want the drone to cover the area farthest away from the CS to minimize double coverage. One flight path that would accomplish this is: Go back and forth parallel to V 1 V 2 until the entire region is covered.
This path is shown in Figure 8. However, the fact that the drone has limited range means that, in practice, the drone will have to break off from this pattern and return to the CS to recharge. After recharging, the algorithm will restart by having the drone travel to V 0 , as seen in Figure 9. As can seen in the figure, after each drone path segment, the remaining uncovered triangle is similar to the original triangle but successively smaller and smaller. Consequently, we can reapply the same algorithm to each smaller triangle to obtain a path that covers the entire original triangle. A flowchart for the algorithm used to generate path segments in shown in Figure 10. Pseudocodes for the algorithm and its subfunctions are shown in the appendix. Once the path is completed on a the transformed triangle, we apply the inverse transformation to obtain the actual path.    We illustrate the algorithm process in more detail by showing the construction of the initial path segment. We consider the case where the CS location (denoted by V 0 ) is located at (0, 0) as shown in (see Figure 7). In this case, the initial segment is parallel to V 0 V 1 at a perpendicular distance ρ, as shown in Figure 11. The endpoints of this initial segment are denoted by V 0 ≡ ( x 0 , y 0 ) and V 1 ≡ ( x 1 , y 1 ), as shown in the figure. Note that the point V 1 is on the segment V 1 V 2 . The mathematical expressions for V 0 and V 1 are: Since we are beginning at V 1 , we must also include the short segment from V 1 to V 1 in the path. Once the segment from V 0 to V 2 is completed, then the remaining uncovered area is a smaller triangle with vertices {V 0 , V 1 , V 2 } where V 0 and V 1 are given by: ; as shown in Figure 11. At this point, the algorithm has a choice. If there is still enough remaining energy, the drone will pass from V 1 to V 2 by way of point V 1 , where Consequently, the remaining uncovered area is another smaller triangle with vertices {V 0 , V 1 , V 2 } where V 1 and V 2 are given by: On the other hand, if there is not enough energy, the drone will return from V 1 to V 0 along a path that is parallel to V 0 V 1 at a distance ρ. The equations for this return segment are similar to Equation (6) except that x 1 , y 1 , x 2 , y 2 are replaced by x 1 , y 1 , x 2 , y 2 respectively.
In case the drone does travel to V 2 , the drone once again has two options: either return to recharge or travel back towards V 1 . Let us suppose the drone needs to recharge: then it will travel parallel to V 0 V 2 , and then return to the CS at V 0 to recharge. The mathematical expression for V 2 and V 0 are: The remaining uncovered area is another smaller triangle with vertices {V 0 , V 1 , V 2 } where V 0 and V 2 are given by: After the drone is recharged, the algorithm will start again by sending the drone to V 0 and beginning its journey to cover the remaining area of the triangle. Figure 12 shows the final drone path produced by the system for the field configuration shown in Figures 5 and 6.  Figure 5 shows as white lines.

Performance Evaluation
The two performance parameters of practical interested are the number of CSs and the total mission time. The efficiency of the algorithm with respect to these two parameters was measured as follows.
It is known that the most efficient cover of a plane by circles is obtained by a circumscribing the hexagons in a "honeycomb" lattice [54]. This implies that regular hexagonal-shaped cells cover the maximum area with the fewest possible cells, given the cells' radius is constrained to be less than R. To measure the CS efficiency for a configuration with CS coverage radius R, we divided the mean area per Voronoi cell by the area of a hexagon with radius R. This gives a dimensionless measure that can be compared across different field sizes and shapes.
As far as mission time, mission time may be broken up into (travel time ) + (recharge time). The travel time was computed as (travel distance) / cruise velocity. The recharge time is computed as (travel time) × (power consumption) / recharge power. It follows that the tour time is equal to: Tour time = (Tour distance) cruise velocity + power consumption (cruise velocity)(recharge power) .
Equation (11) shows that tour time is proportional to tour distance. It follows that minimizing tour time is equivalent to minimizing tour distance.
We have mentioned in Section 2.2 that a maximally efficient boustrophedon path that passes over all field points exactly once is not possible for drones of limited range. Instead, a drone may need to travel multiple times over the same area, consequently increasing the total travel distance.
We may derive a lower bound on the total travel distance for a complete surveillance mission using the following notations: • D(r) is the total distance that the drone flies while its distance to the nearest CS is greater than r, for 0 ≤ r ≤ R. Note that D(R) = 0, while D(0) is equal to the total distance that the drone flies during its surveillance mission (denoted by D). • A(r) is the field area that lies farther than a distance r of the nearest CS, for 0 ≤ r ≤ R. Note that A(R) = 0, while A(0) is the total field area (denoted by A).
Then we have: The inequality in (12) expresses the fact that on each single charge, the drone can fly a total distance at most d of which at least 2r is within radius r. We also have: In (13) the left-hand side represents the total area which the drone camera's FOV passes over while flying farther than r from the nearest CS, which must be greater than or equal to the total field area that is farther than r from the nearest CS. Combining (12) and (13) gives: A dimensionless and scalable measure of efficiency was obtained by taking the theoretical minimum mission distance derived in (14) and dividing by the achieved mission distance. In this way, a mission efficiency of 1 represents the shortest possible mission for the given CS configuration.

Simulation Specifications
In order to evaluate the effectiveness of the algorithm, we applied it to a variety of scenarios. Parameter values used are listed in the last column of Table 1. These parameters were chosen to match the specifications of the Mavic Air drone with 2375 mAh battery [55], together with the Heisha C300 solar-powered charging pad [56]. Three different field shapes (3:1 rectangle, square, octagon) were used, as well as three different total areas (25, 50, 100 km 2 ). For each of these 9 field configurations, 100 different random CS configurations were generated, and run for 5 different initial CS coverage radii, as given in Table 1. Smaller initial radii led to solutions with more CSs and shorter mission times.
All simulations were done on a MacBook Pro with a single 8-core 2.3 GHz Intel Core i9 processor and 16 GB RAM, running Python 3.7.6 with a Spyder 4.0.1 interface. Figure 13 shows the runtimes achieved, where each single run of the algorithm is shown as a dot on the graph. Runtimes were under 100 s, even for the largest configuration, demonstrating the algorithm's practical applicability even for large scenarios.

Results
Results are graphically summarized in Figures 14 and 15, while Table A1 gives numerical statistics. In calculating the coverage time, the cruise speed was taken as 25 km and a charge time: flight time ratio of 2:1 was used.
The first row of graphs in Figure 14 shows the expected field per CS. For fields of size 25 km 2 , we may expect to use 1 CS for every 4 km 2 when the CS radius is 2 and about 1 CS for every 11 km 2 when the CS radius is 4. For larger areas, a much greater efficiency can be achieved; 1 CS for every 6 km 2 when the CS radius is 2 and about 1 CS for every 20 km 2 when the CS radius is 3.9. This implies that only about 1/3 as many CSs would be required if the larger CS radius is used for the design.
The second row of graphs in Figure 14 measures the CS efficiencies, which ranged from about 35% to above 70%. Greater efficiencies were observed with more regular fields (square and octagon as opposed to rectangle), larger fields, and smaller CS coverage radii. Mission times per field area are shown in the first row of graphs in Figure 15. The similarity between the graphs for different field sizes shows that mission times are roughly proportional to total field area. For small areas the coverage radius had little effect on the total mission time, but for the largest area size total mission times were up to 40% larger for CS coverage radius equal to 4 km compared to CS coverage radius equal to 2 km. Mission efficiencies ranged from about 70% to about 90%. Greater efficiencies were observed with smaller field sizes and smaller CS coverage radii. For the smallest field size, the efficiency was roughly independent of the CS coverage radius.
The information in Figures 14 and 15 is combined in Figure 16 to produce Pareto frontiers showing the tradeoff between mission time and number of charging stations. Information for all of the individual runs is plotted. There is considerable variation in the mission time per area for a given value of covered area per CS, especially for larger coverage areas. The similarity between the three graphs shows that the tradeoff is relatively independent of the field size.

Discussion
The results observed are reasonable from a geometric point of view. For smaller field sizes and large CS coverage radii, it is unavoidable that a large portion of the CS coverage region will lie outside of the field. This is why lower CS efficiencies were observed in Figure 14 for the 25 km 2 field sizes. As the size increases, a smaller proportion of Voronoi cells is near the field boundary, so the wasted area is less. For large cells, the increased CS efficiency comes at the cost or lower mission efficiency. This is because greater CS efficiency means larger cells, which implies more area that is far from the nearest CS which will require repeated trips that cover the same area multiple times.

Conclusions
We have developed a software system that generates a comprehensive drone surveillance plan for large, inaccessible agricultural or forested areas. Our results indicate feasible numbers of CSs and mission times for a variety of field shapes and sizes. We have demonstrated the adaptability of the software. Future possibilities for research include introducing more realistic conditions such as variable winds, drone's multiple operating modes, etc. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.