Path Planning with Drones at CSP plants

Path Planning with Drones at CSP plants Adrián Gutiérrez Camacho Path Planning with Drones at CSP plants Adrián Gutiérrez Camacho Trabajo de fin de máster que forma parte de los requisitos para la obtención del título de Máster Universitario en Matemáticas por la Universidad de Sevilla. Dirigido por José Miguel Díaz Báñez Agradecimientos: En primer lugar, me gustaría agradecer la inestimable ayuda de mi tutor, José Miguel Díaz-Báñez, que ha sabido encauzar el trabajo y darme los hilos de los cuáles ir tirando a lo largo de su desarrollo. También a mi jefe, el profesor Juan Sebastián Valverde por ofrecernos un problema tan interesante y sobre todo realista. La aplicación de las matemáticas es algo que siempre me ha fascinado. También querría agradecerle a mi compañero José Morales, experto ya en el arte de las inspecciones en plantas CSP, por darme toda la información de la que disponía para poder hacer el trabajo de forma más realista. De igual manera, le agradezco a mi compañero Miguel Ángel Pérez-Cutiño las aportaciones de ideas y correcciones que ha hecho del trabajo. Por último, pero no menos importante, me gustaría agradecerle a todas las personas que forman parte de mi vida personal y que son las que hacen posible mi dedicación al trabajo de investigación: a mi pareja Miriam, que es el motor que me impulsa a seguir siempre adelante, a mis amigos y amigas, que son capaces de sacarme de la cabeza todas las cuentas y ecuaciones y dejarme solo una sonrisa, y a mi familia, sin la cual nunca me hubiera propuesto estos objetivos y por supuesto nunca los habría alcanzado. Abstract The goal of this work is to apply mathematics knowledge and skills to efficiently solve a practical problem posed by the industry. We study an actual problem related to the inspection of Concentrated Solar Power (CSP) plants. Due to the big extension of solar fields, Unmanned Aerial Vehicles (UAV), commonly called drones, are used to inspect all the tubes of the CSP plant. We introduce a new problem, named the drone CSP inspection problem, that aims the computation of the tours to be performed by the drone in order to cover the CSP plant so that some penalization function is minimized. Specifically, we take into account two objective functions: the total time or the number of refills. First, we model the energy consumption of the UAV and the individual time inspection costs in a realistic fashion and use them as inputs for the procedures described. We also propose several formulations adapting classical optimization problems. In addition, we prove that this particular problem is NP-complete and develop some heuristics. An extensive comparison against the current approach adopted by the industry shows best performance of our algorithms, saving a considerable amount of time for inspection. Resumen El objetivo de este trabajo es aplicar conocimiento y habilidades matemáticas para resolver eficientemente un problema práctico propuesto por la industria. Estudiaremos un problema real relacionado con la inspección de plantas de concentración solar de potencia (CSP). Debido a la gran extensión de los campos solares se utilizan vehículos aéreos no pilotados (UAV), comúnmente llamados drones, para inspeccionar todos los tubos de la planta CSP. Introduciremos un nuevo problema, el problema de inspección CSP con drones, donde se propone calcular las trayectorias a realizar por el dron de manera que se cubra la planta CSP mientras se minimiza una cierta función de penalización. Concretamente, tendremos en cuenta dos funciones objetivo: el tiempo total de inspección y el número de recargas que el dron necesita. Primero, modelaremos el consumo de energía del UAV y los tiempos individuales de inspección de forma realista y los usaremos como entrada de los procedimientos descritos. Propondremos varias formulaciones adaptando problemas de optimización clásicos. Además, probaremos que este problema particular es NP-completo y desarrollaremos algunos heurísticos. Comparando éstos con procedimiento actual adoptado por la industria, probamos que nuestros algoritmos tienen un mayor rendimiento, ahorrando una considerable cantidad de tiempo total de inspección. 4path planning with drones at csp plants Figure 1.3: Absorber element of a parabolic trough collector (borrowed from [1]) tems by using long rows of flat or slightly curved mirrors to reflect the sun rays onto a downward facing linear receiver. The receiver is a fixed structure mounted over a tower above and along the linear reflectors. The reflectors are mirrors that can follow the sun on a single or dual axis regime. The main advantage of LFR systems is that their simple design of flexibly bent mirrros and fixed receivers requires lower investment costs and facilitates direct steam generation, thereby eliminating the need of heat transfer fluids and heat exchangers. LFR plants are however less efficient than PTC and SPT in converting solar energy to electricity. Moreover, it is more difficult to incorporate storage energy into their design. Parabolic dish collectors (PDC) concentrate the sun rays at a focal point supported above the center of the dish. The entire system tracks the sun, with the dish and receiver moving in tandem. This design eliminates the need for a HTF and for cooling water. PDCs offer the highest transformation efficiency of any CSP system. PDCs are expensive and have a low compatibility with respect of thermal storage and hybridization. Promoters claim that mass production will allow dishes to compete with larger solar thermal systems. Each parabolic dish has a low power capacity (typically tens of kW or smaller), and each dish produces electricity independently, which means that hundreds or thousands of them are required to install a large scale plant like built with other CSP technologies. Within the commercial CSP technologies, parabolic trough collector (PTC) plants 1. introducing the problem. state of the art 5 are the most developed of all commercially operating plants. In terms of land occupancy, PTC requires more land than SPT and LFR to produce a given output. PDC has the smallest land requirement among CSP technologies. Water requirements are of high importance for those locations with water scarcity, for example, in most of the desserts. CSP requires water for cooling and condensing processes. Dry cooling (using air instead of water) can be an effective alternative. Figure 1.4: Thermal energy storage system in parabolic trough collector plant (borrowed from [1]) 1.1 The drone CSP inspection problem In this work, the goal is to study a drone path planning problem in CSP plants. We will focus on PTC ones, but this problem can be, at least partially, extended to other CSP plant types or even to another optimization problems. CSP plants are composed of regions or solar fields with absorber tubes or HCE (heat collector element) which warm up some fluid to very high temperatures. The HCEs are covered with glasses which are usually broken due to vibrations and other phenomena, and this leads to heat losses (a pipe whose glass is broken will have leak of temperature in the fluid inside) which, at the end, provoke a reduction in electrical availability. Due to big extension of solar fields, it can be very hard to identify broken pipes in order to fix them and recover the good performance of the plant. To overcome this issue, a drone is used to inspect all the absorber tubes of the CSP plant, one-by- one. This drone flies over the plant and takes thermal pics by means of an special 6path planning with drones at csp plants camera. Then, the drone CSP inspection problem or CSP problem, for short, asks for the computation of the trajectories followed by the drone along the CSP plant so that some penalization function is minimized: the needed time or the number of refills. The constraint of the problem is the autonomy of the drone that can be set to 30 min. Therefore, our goal is to study the complexity of the CSP problem and propose optimal or sub-optimal solutions. To accomplish that, we explore various related optimization problems in the literature, adapt them to our CSP problem and develop heuristic algorithms based on the known heuristics for these problems. We use the Python programming language to develop tools for solving, representing and computing the solutions associated to each approach. Moreover, we are specially interested in the practical value it might have for companies engaged in drone inspections. In fact, this work has been developed in the framework of Virtualmechanics1, a competitive company in the field of plants inspections. They put a real problem on the table and also have been reporting to us all the information of particularities of the problem so that we could work with a realistic model. One of the purposes of this thesis is to improve the performance of the current drone inspections that this company is conducting in CSP plants. Figure 1.5: Typical CSP plant top view Firstly, in order to state the model, we need to introduce the targets or points to be visited by the drone. In a CSP plant the absorber elements are lied in couples: 1https://virtualmech.com/ 1. introducing the problem. state of the art 7 the upstream tube and the downstream tube, one aside the other. Another constraint of our problem is that these couples have to be analyzed consecutively by the drone (although they can be inspected in any order). Because of this fact, we will model the upstream and downstream pipes as one point to visit, instead of two. Then the input of the problem is just a cloud of npoints Sand we want to describe a set of tours starting and ending at a given base station so that all the points are visited. In addition, this couples are placed making big extensions called batches. Figure 1.6 shows two faced batches with 41 elements to be covered. So, actually, we can identify our CSP plant with a cloud of npoints in R2which are also separated in bbatches. Unfortunately, due to the limited battery life of the used small drone, the size of the problem becomes very large if we consider the entire CSP plant as an instance. Moreover, the CSP plants are delicate environments where accuracy is required and also laws usually do not allow to fly drones far away from the pilot. There are similar limitations associated with the RTK system and the drone antenna, which give the drone special accuracy in those environments. Thus, we will focus in how to efficiently cover a pair of faced CSP batches. Figure 1.6: Picture with two faced CSP batches. The upper batch has 20 couples of pipes and the lower batch has 21. We identify each couple with an inspection node The problem is how specifically the drone can visit all the inspection points so that the total time is minimized. More formally, given the complete digraph whose nodes are the npoints and the edges cost is given by the time to travel between the nodes, which is the set of tours or cycles to visit all nodes of the graph in the minimum time?. The constraints for each tour are the following: 8path planning with drones at csp plants All node which is different from the depot node (base station) must appear exactly in one tour. On the other hand, all the tours must start and end at the depot node. A set of tours verifying this and such that its union visits all the nodes in the given graph is called a covering set. The total energy consumption effectuated by the drone during a tour can not be grater than the total battery level of the drone. In fact, we will use a 75% to 90% of the maximum battery level because of security reasons. A tour with this characteristic is called feasible tour. The CSP instances are defined by: the complete digraph Ggenerated from the n inspection points and the depot node, the energy costs wij, the inspection time costs tij and the maximum level of battery W. In this scenario we consider two criteria giving rise to the following CSP problems: 1. CSP problem 1: How could we design a covering set of feasible tours for a given CSP instance minimizing the total inspection time? 2. CSP problem 2: How could we design a covering set of feasible tours for a given CSP instance minimizing the total number of tours or refills the drone makes? In this work, we usually use the term CSP problem for the two problems defined above. However, the problems are different, as we will see in future sections, because a solution for one could not be a solution for the other. However, due to the symmetry of the instances, we will see that the two solutions can be quite similar. Moreover, we could approximate the minimum time cost covering set by the covering set with the lesser number of tours. We could also wonder about where to place the battery station but, in fact, the places where it can be put are limited to the roads around the CSP plant, with some of them crossing it. This is an interesting facility location problem that is out of the scope in this work. We will suppose that the battery station is placed in the road below or above any of the two faces batches, in a position whose xcoordinate is the mean of the xcoordinates of the inspection nodes in the nearest batch. 1.2 Related classical problems In this section we will refer to some classical path planning problems that appear in the literature that come to mind to relate with the CSP problem. We also will 1. introducing the problem. state of the art 9 see some recent related work. The notations and concepts will be useful in the CSP modelling that we are going to develop in the next chapters but, we do not use these problems to solve the CSP instances. Suppose we have a directed graph G= (V, E) where Vcontains the nnodes or points to be visited by the drone and Eis the set of all the possible arcs or paths between two different nodes in V, i.e, Gis a complete digraph. The costs are asymmetric because of the wind. Let us define the decision variables xij ∈ {0,1}so that xij = 1 if and only if we use the arc (i, j)in our solution. For sake of simplicity, a tour in Gmay be given as an ordered subset of V,{i1, i2, . . . , in}which means we travel from ijto ij+1,j= 1, . . . , n −1and then from inagain to i1. The first known optimization problem to be related with the CSP problem is the Travelling Salesman Problem (TSP) whose linear formulation can be given ([2]) as: (TSP) min n X i,j=1 cijxij (1.1a) s.t.: X i=1,i=j xij = 1, j = 1, . . . , n (1.1b) X j=1,j=i xij = 1, i = 1, . . . , n (1.1c) ui−uj+nxij ≤n−1,2≤i=j≤n(1.1d) ui∈ {1,2, . . . , n −1}, i = 2, . . . , n (1.1e) xij ∈ {0,1}(1.1f) This is known as the MTZ formulation of TSP. The objective function (1.1a) minimizes the total travel cost. The constraints (1.1b) and (1.1c) ensure that each node of Gis reached by an arc exactly once and that only another one leaves from it. (1.1d) are the tour elimination constraints: using the dummy variables u2, . . . , un(excluding 1 as it is a depot node) which are forced to grow up in value with the order that the tour has. The last constraints (1.1e) and (1.1f) are the domains of the decision variables. However, if we want to be more realistic, we should deal with the fact that the drone usually can not inspect all the nodes in graph without returning to the battery station (BS). In fact, the limit of duration in past real CSP inspections has been usually about 9of these points (around 30 minutes of battery). We have a certain kind of capacity limitation in our problem. Because of that, the duration of the battery will be an additional parameter to consider in the model. Also, we must include the BS in our graph as the depot or 0-node. Capacity and depot are key words that remind to the Vehicle Routing Problem (VRP), which we have formulated here in a Dantzig, Fulkerson and 10 path planning with drones at csp plants Johnson way: (V RP) min n X i,j=1 cijxij (1.2a) s.t.: X i∈V,i=j xij = 1, j ∈V\{0}(1.2b) X j∈V,j=i xij = 1, i ∈V\{0}(1.2c) X i∈V xi0=K(1.2d) X j∈V x0j=K(1.2e) X i∈SX j∈S xij ≥r(S),∀S⊆V\{0}, S =∅(1.2f) xij ∈ {0,1},∀i, j ∈V(1.2g) In this formulation, Kis the number of available vehicles (the number of tours or charges needed by the drone to visit all the inspection points) and r(S)the minimum number of vehicles needed to serve a set S. The VRP optimal solution is the minimal low-cost set of tours, which start and end in a depot node, that reach all the nodes of the graph while ensuring enough capacity in each of them. In our problem we have a vehicle (drone) which has to visit nnodes, starting from the BS depot, while watching its battery duration. When the drone reaches a minimum amount of battery it has to come back to the base station and start a new one tour. The main goal is to plan the tours in order to minimize the time we spent in the CSP inspection. Although this model has similarities with our problem, this formulation is still too simple. A more advanced version of VRP is the Capacitated Vehicle Route Problem (CVRP) [3]. Its formulation is pretty similar to the previous one, but we must add some new components. In this case, let G= (V, H, c)a complete directed graph with V= {0,1,2, . . . , n}as the nodes and H={(i, j): i, j ∈V, i =j}as the arcs, where node 0is the depot for a fleet of pvehicles with the same capacity Q. Each node i∈V\{0}has a certain positive demand di≤Q. We also consider a travel cost cij associated with each arc (i, j)∈H. The minimum number of vehicles to needed to serve all customers is ⌈Pn i=1 di Q⌉(note we could also add a battery efficiency coefficient to consider uncertainty). The binary decision variables xrij ∈ {0,1}are defined to 1. introducing the problem. state of the art 11 indicate if vehicle r,r∈ {1,2, . . . , p}, traverses an arc (i, j)∈Hin an optimal solution. Then, the linear programming model of the CVRP can be written as: (CV RP) min p X r=1 n X i=0 n X j=0,i=j cijxrij (1.3a) s.t.: p X r=1 n X i=0,i=j xrij = 1, j ∈V\{0}(1.3b) n X j=1 xr0j= 1, r ∈ {1, . . . , p}(1.3c) n X i=0,i=j xrij = n X i=0 xrji, j ∈V, r ∈ {1, . . . , p}(1.3d) n X i=0 n X j=1,i=j djxrij ≤Q, r ∈ {1, . . . , p}(1.3e) p X r=1 X i∈SX j∈S,i=j xrij ≤ |S| − 1,∀S⊆V\{0}, S =∅(1.3f) xrij ∈ {0,1},∀r∈ {1, . . . , p}, i, j ∈V, i =j(1.3g) The objective function (1.3a) minimizes the total travel cost of all the vehicles in the fleet. The degree constraints (1.3b) ensure that each node is visited by exactly one vehicle. The flow constraints (1.3c) and (1.3d) guarantee that each vehicle can leave the depot only once, and the number of the vehicles arriving at each customer and entering the depot is equal to the number of the vehicles leaving. The capacity constraints (1.3e) make sure that the sum of the demands of the customers visited in a route is less than or equal to the capacity of the vehicle performing the service. The tour elimination constraints (1.3f) ensure that the solutions contains no cycles disconnected from the depot. The last constraints (1.3g) specify the definition domains of the variables. This model is known as a three-index vehicle flow formulation and has the drawback that the number of inequalities in (1.3f) grows exponentially with the number of nodes. A more related problem is studied in [4]. Although the costs and topologies in their graphs are different to the ones in a CSP problem, the goal is the same: a drone has to visit a set of points S, while recharging its battery in some locations of a discrete set C, minimizing the total trip time. They studied how to estimate, through a 12 path planning with drones at csp plants regression model, the battery consumption of determined drone flight actions. The regression uses the basic kinematics of the drone as prediction variables. We want to take advantage of this work using their model to compute realistic energy consumption weights. The authors also develop some heuristics to give approximations to the problem. Mainly, their contributions are: 1. A model the energy consumption of drones, considering various flight scenarios. 2. A study on the joint problem of flight tour planning with recharging optimization for drones with the goal of completing a tour mission for a set of locations of interest in the shortest time. 3. A real implementation of their algorithms in an intelligent drone management system, which supports real-time flight path tracking and re-computation in dynamic environments. Since we borrow some notation of that paper, we include here more details on their model and frormulation. Consider a set of cities Sthat a drone has to visit, and a set of charging stations Cwhere a drone can charge its battery. The base location or depot of drones is denoted by {v0}. The problem is to find a sequence of locations or flight plan in S∪C, such that the drone can start its tour at v0, visit all sites in Sand then returns to the depot, with the objective of minimizing the total trip time, while keeping the SoC (state-of-charge) within the operational range. Given a pair of locations (u, v), the flight path between them is noted by l(u, v), and the flight time by τ(u, v). Let f(l(u, v), τ(u, v)) be the energy consumption related to a drone travelling along l(u, v)within flight time τ(u, v). Also, let V= S∪C∪ {v0}. For a drone flying between two sites u, v ∈V, there is a battery consumption by an amount of ndf(u, v). If the drone returns to a charging station u∈C, it can recharge its battery by an amount of energy ncb(u). The coefficients nd and ncstand for discharging and charging efficiency performance. In addition, when recharging its battery there will be a incurred charging time denoted as τc(b(u)). Let T⊆S∪C∪{v0}be a flight plan and denote its k-th location by Tk. In order to find a flight plan which minimizes the total trip time (travel and charging time) it is needed to impose that T1=T|T|=v0. Let xkbe the SoC when arriving at Tk, assuming that xk∈[B0, B1], that is, the SoC lies within a feasible range. The lower bound B0ensures sufficient energy for the drone to return to the depot. Then, the drone flight plan optimization problem (DFP) is formulated as follows: 1. introducing the problem. state of the art 13 (DFP) min T,b(·),x |T|−1 X k=1 τ(Tk, Tk+1) + |T| X k=1: Tk∈C τc(b(Tk)) (1.4a) s.t.: T1=T|T|=v0(1.4b) S⊆T⊆S∪C∪ {v0}(1.4c) xk=           xk−1−ndf(l(Tk, Tk+1), τ(Tk, Tk+1)) if Tk∈S xk−1+ncb(Tk+1)−ndf(l(Tk, Tk+1), τ(Tk, Tk+1)) if Tk∈C (1.4d) B0≤xk≤B1(1.4e) This formulation accounts for the situation of recharging the drone without filling up the entire battery. Our CSP inspection problem is not as complex in that aspect: at Cthe drone replaces its battery obtaining xk=B1. Also, for us, C={v0}, that is, our depot is the only charging station. The authors then developed some heuristic algorithms to apply its DFP problem and found lemmas to support its work with some theory. In [5], the authors study the Drone Arc Routing Problems (DARP) and its relation with the Postman Arc Routing Problems. The main difference is that drones can travel directly between any two points in the plane without following the edges of the network. Thus the DARP is a continuous optimization problem with infinite feasible solutions. However, they solve it as a discrete optimization problem by approximating curves as polygonal chains. In their research, they consider a set of drones that must jointly service a set of edges of a network. Drone start and end at a specified vertex (depot) and have a route length limit that applies to the total route (capacity). Arc routing problems consist of, given a network, and given a set of lines that are required to be covered, each one with an associated cost, finding a tour covering all the required lines with total minimum cost. Typical examples are the Chinese Postman Problem , the Rural Postman Problem. They give same results relating optimal drone tours with optimal postman tours. Specifically, the Drone Arc Routing Problem is: Given a set of lines or curves on the plane, each one with an associated service cost, and a point called the depot, and assuming that the cost of deadheading between any two points of the plane is the Euclidean distance, find the minimum cost tour starting and ending at the depot that services all the given lines and curves. They develop 20 path planning with drones at csp plants Turning: optionally, but very often, a drone effectuates a turn to head its objective (in general, a node to inspect) and then, at the arrival, to head parallel to the pipe. In addition, during the inspection phase, inside a node, the drone needs to turn at least 3 times, in U shape, to visit the entire couple defining a point of interest. It is important to remark that one of the most expensive displacements a drone can do are the ones which imply more turnings, in the sense of increasing the time and energy costs. We will refer to turns as Ag. The notations for each type of action will help us to differentiate states of the energy consumption parameters, as their values change with the action performed. Let [0, T] = [0, t1]∪[t1, t2]∪[t2, T]be the time interval during the performance of a given action. Also, each displacement type (except turnings) will be decomposed in three parts: 1. Acceleration part: the drone increases its speed in the direction of movement (horizontal, vertical or negative vertical) by using propellers. This is an uniformly accelerated rectilinear motion (linear acceleration is negative). This part lasts from t= 0 to t=t1. 2. Uniform part: the drone moves with constant speed in the desired direction. This is an uniformly rectilinear motion (linear acceleration is zero). It lasts from t=t1to t=t2. 3. Braking part: the drone decreases its speed in the direction of movement, until it is 0. This is an uniformly decelerated rectilinear motion (linear acceleration is positive). It lasts from t=t2to t=T. Turning is considered to be uniform motion, so it has neither acceleration nor braking parts. So now, we can distinguish an action or displacement by its type, A∈ A ={At, Ah, Al, Ag}, and the time interval when it occurs, t∈[0, t1],[t1, t2] or [t2, T]. In addition, note that the total acceleration of an object, aT, can be decomposed as linear acceleration, a, plus the gravity acceleration, g. We will define the gravity vector1as g = [0,0,−9.81]m/s2. For example, in order to positively accelerate while ascending, the drone has to generate a positive vertical linear acceleration Z component. Moreover, we must use some formula to compute the energy consumption from drone kinematics. For this we will use [4] and also [10]. In the articles, authors de- 1We will assume that Z axis is perpendicular to ground and positive while ascending. 2. some formulations for the csp problem 21 fine a standard regression model which allows to compute the instantaneously power consumption, ˆ P, as function of s= [vxy,axy,vz,az, ωxy, m]where: vxy and axy are the horizontal speed and acceleration vectors. vzand azare the vertical speed and acceleration vectors. ωxy is the horizontal wind speed vector. mis the constant mass of the drone payload. Then, ˆ Pcan be computed2as follows: ˆ P(s) =   β1 β2 β3   T  ∥vxy∥ ∥axy∥ ∥vxy∥∥axy∥ +  β4 β5 β6   T  ∥vz∥ ∥az∥ ∥vz∥∥az∥ +  β7 β8 β9   T  m vxy ·ωxy 1 =βTˆ s where ˆ s= [∥vxy∥,∥axy∥,∥vxy∥∥axy∥,...,1]. If during a period Tthe parameters needed for computation remain unaltered, the energy consumption can be computed as ˆ E=ˆ P·T. However, the reader would likely observe that 2/3 of each displacement type which we have defined is not constant in speed. To tackle this, we will generalize the energy consumption computation in a more analytic fashion. Let sbe now a time dependent function modelling the kinematics state of the drone, which also depends on the action type performed. Fixed A∈ A, we could define the instantaneous energy consumption during an infinitesimal period dt (where the state sremains constant) as: dˆ E(s(t)) = ˆ P(s(t)) ·dt =βTˆ s·dt and then compute ˆ E(s(t)) integrating along [0, T]. In the Appendix Athere is a more detailed explanation of how we can compute kinematics states for each action in order to obtain ˆ E. Now, we have all the ingredients needed to model the energy costs for visiting a point of interest in our CSP plant. We will distinguish between trips where the drone is coming from the base station node (beginnings), trips where the drone is going to the base station node (endings) and trips where the drone is travelling between points of interest (normal trips). In beginnings, the drone starts with a takeoff, heads (turning) to the next node to visit and then does an horizontal displacement. The endings start heading the depot node, horizontal movement is performed and then the landing happens. The normal trips are just turnings to the following node and 2Obviously we will need the values of β= [β1, . . . , β9]but we will discuss this question later. 22 path planning with drones at csp plants horizontal displacements between nodes. There is an important detail in how the energy costs are modelled, they depend on current and previous visited node: let wt ij be the energy consumption of the drone when it is travelling from node ito node j, i, j ∈ {0,1, . . . , n}, i =j, that is, wt ij =ˆ E(s(i, j)) where, here, s(i, j)models the parameter state the drone will present if it travels from ito j. However, we usually don’t just visit the nodes but we also inspect them (except for the depot). The inspection cost associated to the node i= 0,we i, can be computed creating fictitious nodes along pipes and performing normal trips between them. Recall that when we reach a point of interest and the inspection begins, the drone flies until the end of pipe, jumps to the other part of the couple, and finishes at the beginning of it. So, actually, the start position at node iand the end position after inspection doesn’t match. But we can just create fictitious nodes i1, i2, i3, i4and compute the energy costs by parts, we i=ˆ E(s(i1, i2)) + ˆ E(s(i2, i3)) + ˆ E(s(i3, i4)), where i1stands for the start position coordinates at i,i2represents the coordinates at the end of the first pipes, i3plays as the coordinates at the end of the other part of the loop and i4defines the end position coordinates (see Figure 2.3). Clearly, we will use different iand jnodes coordinates when computing wt ij (output coordinates for i and input coordinates for j) that when getting wt ji (output coordinates for jand input coordinates for i). The depot node is the only one with the same input and output coordinates. Figure 2.3: Fictitious nodes scheme for inspection cost computation at nodes i, j with trajectories of inspection (in blue) and between nodes (in red). Besides the energy consumption costs, the time costs tt ij model the spent time of going from ito j. To compute them, we could use basic kinematics formulas and the 2. some formulations for the csp problem 23 distance between the nodes, d(i, j)(if the speed is constant, the spent time equals to the distance divided by the speed). This is deepened in the Appendix A. There is also an inspection time cost, te i, which can be computed creating fictitious nodes as we i. To sum up, we define tij =tt ij +te iand wij =wt ij +we ias travel-inspection time costs and energy weights between nodes iand j, respectively. In order to model the costs of coming from the base station and returning to it, we create a fictitious node 0′which is at the same XY coordinates as 0but at some altitude above it, z0. Then, the costs wt 0jcan decompose in a takeoff from 0to 0′and an horizontal displacement from 0′to j(similar for wt i0, which implies a landing from 0′to 0). Evidently, wi0=wt i0because in this node there is no inspection. However, to be more realistic, the time costs at depot node will be ti0=tt i0+tr, where tris the recharging or battery replacement time cost. 2.2 Multiple Knapsack Problem approach In this section we look at the CSP problem as an instance of the MKP. Suppose that now the npoints to visit are items from a collection. We want to pack these items in k bins (tours or cycles) with kin {1, . . . , n}. Each of the items has a weight (associated energy consumption) and a profit, which is related with the time cost. At the same time, each bin has a maximum weight capacity (total battery energy when it is full of charge), which we denote by W. Then, we would like to choose the items that will be in each bin in the most efficient way, i.e, maximizing the profits while ensuring the weight of formed bins is less than W. The Multiple Knapsack Problem is known in the literature as an immediate extension of Knapsack Problem, and which is NP-hard in the strong sense [11]. The basic formulation for Multiple Knapsack problems can also be found in [11], and it is as follows: let nthe number of items, mthe number of knapsacks (k≤n) and let pj=the profit of item j,wjthe weight of item jand cithe capacity of knapsack i. We want to select mdisjoint subsets of items so that the total profit of selected items is a maximum, and each subset can be assigned to a different knapsack whose capacity is no less than the total weight of items in the subset: 24 path planning with drones at csp plants (MKP) max z= m X i=1 n X j=1 pjxij (2.1a) s.t.: n X j=1 wjxij ≤ci, i ∈M={1, . . . , m}(2.1b) m X i=1 xij ≤1, j ∈N={1, . . . , n}(2.1c) xij = 0 or 1, i ∈M, j ∈N(2.1d) Where xij = 1 if item jis assigned to knapsack iand xij = 0 otherwise. When m= 1 the MKP reduces to the 0-1 (single) knapsack problem. However, we cannot use this formulation directly. Firstly, we want to pack all the items in any of the bins, but the constraint (2.1c) admits to not assign a item to any bin. In fact, another problem which is related to MKP and CSP problem is the Generalized Assignment Problem, which ensures that all the tasks or items are assigned to exactly one task processor or bin. Also, the profits and weights that model our problem depend on the previous visited point or, in knapsack terms, the last added item. These complications will be tackle in the following. The main idea is to consider the nodes to visit as items and the tours as knapsacks with some capacity, i.e, the energy battery level when the drone is full of charge. In the following, we will develop a Mixed Integer Programming formulation which generalizes the MKP so that we can use it to solve instances of our particular problem. 2.2.1 Generalized MKP Firstly, let V={0,1, . . . , n},V={1, . . . , n}and K={1, . . . , m}be the items and knapsacks sets, respectively, where m≤nis fixed. If we recalled MKP formulation from Section 2, we would like to have just one wjand pjby node. However, the profits and the weights usually depend on previous visited node. To overcome this, let: yij =(1if item i∈Vgoes before item j∈V 0otherwise 2. some formulations for the csp problem 25 Then, we can compute the weight of node jas wj=P i∈V wijyij. To compute the profit of node j pj, suppose that we have pij as the profit for going from node ito j, then pj=P i∈V pijyij. On the other hand, let: xjk =(1if item j∈Vis assigned to knapsack k∈K 0otherwise Actually, xjk are the main decision variables in our problem that assign nodes to knapsacks. The extra yij variables are needed for our CSP problem. Then, the MIP formulation for the GMKP could state as follows. (GMKP) max z=X i=j∈VX k∈K pijyijxikxjk (2.2a) s.t.: X i=j∈V wijyijxikxjk ≤W, ∀k∈K(2.2b) X k∈K xjk = 1,∀j∈V(2.2c) X j∈V xjk =X i=j∈V yijxikxjk,∀k∈K(2.2d) 2yij ≤max k∈Kxik +xjk,∀i=j∈V(2.2e) X k∈K x0k=m, X i∈V yi0=m, X j∈V y0j=m(2.2f) X i∈V yij = 1,∀j∈V , X j∈V yij = 1,∀i∈V(2.2g) Li+wijyij −W(1 −yij)≤Lj(2.2h) yijwij ≤Lj≤W, ∀i∈V, ∀j∈V , i =j, (2.2i) Lj∈R,∀j∈V , xjk, yij ∈ {0,1} ∀i=j∈V, k ∈K(2.2j) The objective function (2.2a) maximizes the total profit in all the knapsacks. The terms of the sum, pijyijxikxjk, mean that we are only adding a profit when the items i and jare assigned, in this order, to the same bin, i.e, we visit node iand then node jin the same tour k. The first two constraints, (2.2b) and (2.2c), correspond to the ones at MKP formulation, with slight changes: for example, we impose that every item which is not the depot must be assigned to some bin. Also, the capacities between knapsacks 26 path planning with drones at csp plants remain constant and equal to W. The capacity constraints are formulated using the same idea as objective function, we are only adding a given weight wij when item jis assigned to bin kpreceded by item i, which is also assigned to bin k. Constraint (2.2d) ensures that the number of items assigned to a given bin kmatches the sum of item pairs (i, j)that are assigned to bin kbeing jpreceded by i. Constraint (2.2e) ensures that if jis preceded by i, then they are assigned to the same bin: xik, xjk ∈ {0,1}so maxk∈Kxik, xjk ∈ {1,2}(it can be zero because we have to assign iand jto some bin). This maximum will be 2 if and only if exists k∈Kwith xik =xjk = 1. So, if yij = 1 (iis assigned before j) the inequality implies that i, j are assigned to the same bin. Constraints at (2.2f) mean that the depot node is assigned to all the minitialized bins and also that it goes before and goes after other items the same number of times. This means that we start and end tours at this node. Constraints at (2.2g) force all the items j∈Vto go before and after another item, that is, all the nodes are reached exactly once in a tour. Constraints (2.2h) are useful to ensure the continuity of the route and to eliminate tours that doesn’t include the depot node. The continuous variables Lj, j ∈Vshow the drone battery consumption after visiting node j. We use the constraints at (2.2i) to define lower and upper bounds for Lj. The load or accumulated battery consumption Ljmust be greater or equal than the weight wij only if igoes before jin the same bin. On the other hand, the load is not allowed to surpass the battery capacity W. Finally, constraints at (2.2j) define the domain of load and decision variables and indexes. 2.2.2 Solver for GMKP In order to obtain exact solutions to our CSP problem instances via the GMKP formulation we will need some sort of solver. In this work we have decided to use Gurobi. The Gurobi Optimizer is a commercial optimization solver for linear programming (LP), quadratic programming (QP), quadratically constrained programming (QCP), mixed integer linear programming (MILP), mixed-integer quadratic programming (MIQP), and mixed-integer quadratically constrained programming (MIQCP). The Gurobi Optimizer supports a variety of programming and modeling languages including: C++, Java, MATLAB, R and Python, and also links to standard modeling languages as AMPL. The selected programming language in this work was Python, so this solver fits enough our needs. Also, we acquired an academic license in order to access a bigger computation power. The Gurobi Optimizer also includes a number of features to support the building of optimization models including support for: 2. some formulations for the csp problem 27 A powerful Python modeling API3that integrates with vector and matrix objects from NumPy and SciPy, two of the main libraries in Python. Multiple objective functions with flexibility in how they are prioritized. General constraints as MIN, MAX, ABS, AND, OR, and indicator constraints help avoid having to turn commonly occurring constraints into linear constraints. Client-server computing, to offload optimization or tuning runs to a compute server. Cloud computing, to perform optimization or tuning jobs in the cloud. Distributed computation, to harness the power of multiple machines working together on a single optimization or tuning job. The use of Gurobi in Python is very simple, as we can see in the next example. Suppose we want to optimize the next model: max x+y+ 2z(2.3a) s.t.: x+ 2y+ 3z≤4(2.3b) x+y≥1(2.3c) x, y, z ∈ {0,1}(2.3d) We can program this model with library Gurobipy in Python as follows: 1import gurobipy as gp 2from gurobipy import GRB 3# Create a new model 4m = gp . Model (" mip1 ") 5# Create variables 6x = m. addVar ( vtype = GRB . BINARY , name ="x") 7y = m. addVar ( vtype = GRB . BINARY , name ="y") 8z = m. addVar ( vtype = GRB . BINARY , name ="z") 9# Set objective 10 m. setObjective (x + y + 2 * z, GRB . MAXIMIZE ) 11 # Add constraint : x + 2 y + 3 z <= 4 12 m. addConstr (x + 2 * y + 3 * z <= 4, "c0") 13 # Add constraint : x + y >= 1 14 m. addConstr (x + y >= 1, "c1") 15 # Optimize model 16 m. optimize () 17 for vin m. getVars (): 3Application Programming Interface is a software intermediary that allows two applications to talk to each other 28 path planning with drones at csp plants 18 print (’%s %g’ % (v. varName , v.x)) 19 print (’Obj: %g’ % m. objVal ) The two first lines are imports: we import gurobipy library to have access to its functions and classes through the prefix gp and also import class GRB which contains useful constants we use to define the variables domain. The next line is for creating the model. Lines 6 to 8 define variables x, y, z as binary ones. Line 10 set the objective function of the model. After that we add the constraints of our model. Line 16 starts the process of optimization of the model. The last lines print the optimal solution and value. This was a simple model that also had a simple code implementation. However, Gurobi is flexible enough to allow us implement the GMKP formulation, although that implementation won’t be as simple. 2.3 Bin Packing Problem approach Now we put an instance of the CSP problem as an instance of another well-known problem in the literature, the Bin Packing Problem (BPP). It is similar to the Multiple Knapsack Problem but, in this case, we would like to minimize the number of bins rather than maximize the total profit among them. The BPP is NP-hard in the strong sense [11]. The BPP can be described, using the terminology of knapsack problems, as follows. Given nitems and nknapsacks (or bins), with wj=weight of item j, c=capacity of each bin, assign each item to one bin so that the total weight of the items in each bin does not exceed cand the number of bins used is a minimum. One possible mathematical formulation of the problem is: (BPP) min z= n X i=1 yi(2.4a) s.t.: n X j=1 wjxij ≤cyi,∀i∈N={1, . . . , n}(2.4b) n X i=1 xij = 1,∀j∈N(2.4c) xij ∈ {0,1} ∋ yi, i, j ∈N(2.4d) 2. some formulations for the csp problem 29 where yi=(1if bin iis used 0otherwise xij =(1if item jis assigned to bin i 0otherwise A formulation that generalizes BPP one could give us the minimum number of tours (or returns to the battery station) that the drone can do in a given inspection, m∗. We denote this formulation by GBPP. If we combine GBPP with the GMKP formulation, we could find the best m∗tours to perform the inspection, in the sense of minimizing the total time cost. However, we must highlight that finding the minimum number of tours needed to properly4cover the inspection using GBPP and then solve the GMKP may lead us to a solution worse than if we use more tours, as it is shown in the Figure 2.4. If we consider the instance which is described in the figure, is obvious that we have enough battery to visit all the nodes and then return to base. If we do this, i.e, if we use one tour, the minimum total time cost we can achieve is 8, because to cover both sides we need to spend 3units of time costs, to travel to the other side of the graph costs 4units and return to base costs 1. Now, consider we use two tours, one for each side of the graph, that is, for example, {0,1,2}and {0,3,4}, then the total profit is 3 + 3 = 6, which is less than using one tour. To overcome this, we can solve multiple GMKP instances for m=m∗, m∗+ 1, . . . and observe the behaviour of the obtained solutions. Nevertheless, solutions for GMKP using fewer number of tours usually give the optimal in total time cost. In addition, in order to use the BPP for solving the CSP problem, we must adapt this formulation to one which accounts for the restrictions of our situation. The main differences are related to weights dependence on previous assigned items. Another difference, as we have commented in GMKP, is that we have a depot node acting as a required item to assign to each bin, the 0-item. So, let us now consider the Generalized Bin Packing Problem (GBPP), based on BPP and GMKP, where kvaries in K={1, . . . , n}and we are supposing that we have enough capacity to, at least, visit any of the nodes (adding it to some bin) and returning to base station (adding the 0-item to the same bin), that is, W≥w0jx0jk +wj0xj0k,∀j, k ∈ {1, . . . , n}5. Also, 4Not violating the problem constraints. 5This supposition will be done in all the instances because if they do not verify this, they are infeasible instances. 36 path planning with drones at csp plants Figure 2.5: A 6x10 pentomino puzzle example that: A={1,4,7}, B ={1,4}, C ={4,5,7}, D ={3,5,6}, E ={2,3,6,7}, F ={2,7}. Then the subcollection S∗={B, D, F}is an exact cover of Xusing S, since each element in Xis contained in exactly one of the subsets. Moreover, the given S∗is the only exact cover, as the following argument demonstrates: because Aand Bare the only subsets containing 1, an exact cover must contain Aor B, but not both. If an exact cover contains A, then it doesn’t contain B, C, E or F, as each of these subsets overlap with A. Then Dis the only remaining subset, but the collection {A, D}doesn’t cover the element 2. So, there is no exact cover containing A. On the other hand, if an exact cover contains B, then it doesn’t contain Aor C, as each of these subsets has an element on common with B. Because Dis the only remaining subset containing 5,Dmust be part of the exact cover. If an exact cover contains D, then it doesn’t contain E, as Eoverlaps D. Then Fis the only remaining subset, and the collection {B, D, F}is indeed an exact cover. Clearly, the hardness of this problem grows when there are a lot of subsets in Sand Xis big. Now, we want to explain how we can relate the Exact Cover Problem with the CSP one. Suppose we have a directed graph G= (V, E)where Vis the set of nodes where inspections take place and Eare the possible paths between two of the nodes. What we are looking for is a collection of cycles, C, such that it covers the entire set Vand we don’t repeat a node vin two of the cycles, that is, a ECP where pieces are cycles. Obviously, there might be more than one solution. Since we have a drone battery restriction, many of the cycles that may appear in some solutions could be infeasible, so we will remove them from potential solutions to the CSP problem. Finally, when we have found the feasible solutions, we will keep that with the minimum time consumption or size, in the sense of less tours performed. In order to solve CSP problem by ECP, we must associate weights to cycles: a 2. some formulations for the csp problem 37 cycle will be feasible if it correspond to a feasible tour, that is, the total weight of the cycle doesn’t exceed a given battery consumption bound, W. Also, we can associate profits or costs to each cycle and maximize or minimize them, respectively. Note that each cycle will start and end in the depot node (because the drone has to recharge its battery), so Sis the collection of all the cycles (subsets of V) which visit the depot node, and whose total weight is less or equal than W. Actually, we are not interested in all the cycles, but in the elementary ones (see [12]), which are those where no vertex appears twice (except the first one which matches the last one), and also two elementary cycles are different if they are not cyclic permutations of each other. There are exactly: n−1 X i=1 n n−i+ 1(n−i)! elementary cycles in a complete digraph with nvertices. Computing the total weight of a cycle can be done as follows: let c={u, v, v′, u} the visited nodes in a given one (uis the depot node), then we will compute w(c) = w(u, v)+w(v, v′)+w(v′, u), where w(i, j)are computed like in the previous sections, given a pair of nodes i, j. The ECP has been studied for a long time as it is derived from the famous Set Cover Problem, but nowadays there are few efficient algorithms for solving it. In our research we have found two main references that will help us with the task of solving CSP problem instances. First of all, we need to obtain the collection S, and this is not current because, at the beginning, the feasible cycles in Gare unknown, although we know that there will be a large number of them if there are many vertices and edges. For this first task we will use the Johnson’s Algorithm for finding all elementary cycles in a given digraph, deeply explained in [12]. We must keep only the elementary cycles that visit the depot node. 2.4.1 The Johnson’s algorithm The Johnson’s algorithm finds all the elementary circuits of a digraph in time bounded by O((n+e)(c+ 1)) and space bounded by O(n+e), where there are n vertices, eedges and celementary circuits in the graph. The algorithm is faster than previous ones because it considers each edge at most twice between any one circuit and the next in the output sequence. We can also consider elementary paths, which are paths where no vertex appears twice. 38 path planning with drones at csp plants In this document, we will be using a nonrecursive, iterator/generator version of Johnson’s algorithm which is implemented in the Python library networkx. In fact, we have found that for regular instances of CSP problem, where N≈40, it is very difficult to found all elementary circuits with arbitrary length. To overcome this disadvantage, we use an adapted version of networkx’s implementation which only consider paths whose length is dominated by a given limit L. Also, we have modified it to only search for paths starting at 0node, which are actually the ones we want to obtain. Another possible relaxation we can assume is to only use edges between neighboring nodes, instead of considering a complete digraph. This algorithm allows us to obtain all the elementary circuits with maximal length L≤13 of particular graphs with N≈40 nodes and |E|edges (those related to neighboring nodes in two faced batches associated graphs) relatively quickly. Since the drone usually inspects 9to 10 interests points before returning for charging, it seems like this adapted Johnson’s algorithm is good enough for our purposes. 2.4.2 The Knuth’s Algorithm X To solve ECP in Python, we have used some libraries which implement the dancing links technique of Donald Knuth for his Algorithm X, which solves ECP instances. In the following we will briefly explain how the algorithm X works and also detail the dancing links technique which leads to DLX algorithm and give a small example in order to understand better how it operates. Once he have obtained all elementary circuits Sin graph Gwe can transform the instance of the CSP problem in one for Exact Cover Problem. Sis a collection of subsets of V={0,1, . . . , N}all of them starting at 0 (or at least containing it, we can reorder them). To reduce the computational costs of solving the ECP, we also can discard every cycle in Ssuch that its total associated weight exceeds W, that is, infeasible cycles, obtaining S∗. We will solve the ECP instance (S∗, V )using a Python implementation of the Donald Knuth’s X algorithm, described in [13]. We probably will obtain several possible exact covers Cfor our instance, due to typical big size of S∗. In order to solve the CSP problems 1 and 2 associated with the instance, we’ll select the best cover solution C∗in terms of total inspection time cost and number of cycles used, respectively. The Knuth’s Algorithm X is an algorithm that finds all solutions to an exact cover problem. DLX is the name given to Algorithm X when it is implemented efficiently 2. some formulations for the csp problem 39 using Donald Knuth’s Dancing Links technique on a computer. It is a straightforward recursive, nondeterministic, depth-first7, backtracking algorithm used by Donald Knuth to demonstrate an efficient implementation called DLX, which uses the dancing links technique. Dancing links is a technique for reverting the operation of deleting a node from a circular doubly linked list8. The idea of DLX is based on the observation that in circular doubly linked list of nodes, if we let L[x]and R[x]point to the predecessor and successor of that element, L[R[x]] ←L[x], R[L[x]] ←R[x] remove xfrom the list, while L[R[x]] ←x, R[L[x]] ←x will put xback into the list again, assuming that L[x]and R[x]haven’t been modified. Knuth observed that a naive implementation of his Algorithm X would spend an inordinate amount of time searching for 1’s. To improve this search time from complexity O(n)to O(1), Knuth implemented a sparse matrix where only 1’s are stored. At all times, each node in the matrix will point to the adjacent nodes to the left and right (1’s in the same row), above and below (1’s in the same column), and the header for its column (a special node each column will have, they form the control row consisting of all the columns which still exist in the matrix9). Each row and column in the matrix will consist of a circular doubly-linked list of nodes. The reciprocal procedures for deleting and re-adding an element xto a doubly linked list are extremely useful when updates to a data structure are not indented to be permanent. For example, this occurs in backtrack programs, which enumerate all solutions to a given set of constraints. Backtracking is also called depth-first search. These programs can require to maintain a very large stack and might take too much time, if we use a naive implementation. The Dancing Links technique improve the performance in this algorithms. For example, the exact cover problem, which we can represent by a matrix Aconsisting of 0’s and 1’s such that we want to select a subset of the rows where the digit 1 appears in each column exactly once, is a natural candidate 7Algorithm which searches in tree or graph data structures by starting at a root node and exploring as far as possible along each branch before backtracking. 8A linked data structure that consists of a set of sequentially linked records called nodes. Each node contains three fields: references to the previous and to the next node in the sequence and one data field. 9Also, each column header may optionally track the number of nodes in its column, so that locating a column with the lowest number of nodes is of complexity O(n)rather than O(mn). Selecting a column with a low node count is a heuristic which improves performance in some cases. 40 path planning with drones at csp plants for backtracking. The algorithm X, that we will use for solving it, is detailed in the following: Algorithm 1: Algorithm X input : A output: Solutions to exact cover problem represented by A 1. If Ais empty, the problem is solved; terminate successfully. 2. Otherwise, choose a column, c (deterministically, we can select it by the lowest node count, for example). 3. Chose a row, r, such that A[r, c]=1(non-deterministically). 4. Include rin the partial solution. 5. For each jsuch that A[r, j]=1, 6. delete column jfrom matrix A; 7. for each isuch that A[i, j]=1, 8. delete row ifrom matrix A. 9. Repeat this algorithm recursively on the reduced matrix A. The non-deterministic choice of rmeans that the algorithm generates independent subalgorithms. Each subalgorithm inherits the current matrix A, but reduces it with respect to a different row r. If a column cis entirely zero, there are no subalgorithms and the process terminates unsuccessfully. The subalgorithms form a search tree in a natural way, with the original problem at the root and with level kcontaining each subalgorithm that corresponds to kchosen rows. Any systematic rule for choosing column cin this procedure will find all solutions, but some rules work much better than others. Golomb and Baumert suggested choosing, at each stage of a backtrack procedure, a subproblem that leads to the fewest branches, whenever this can be done efficiently. In the case of an exact cover problem, this means that we want to choose at each stage a column with fewest 1’s in the current matrix A. We will see that dancing links allow us to do this quite nicely. A way to implement algorithm Xis to represent each 1 in the matrix Aas a data object xwith five fields L[x],R[x],U[x],D[x],C[x]. Rows of the matrix are doubly linked as circular lists via the Land Rfields (“left” and “right”) and columns are doubly linked as circular lists via the Uand Dfields (“up” and “down”). Each column list also includes a special data object called its list header. The list headers are part of a larger object called a column object. Each column object ycontains the fields L[y],R[y],U[y],D[y]and C[y]of a data object and two additional fields, S[y](“size”) and N[y](“name”). The size is the number of 1’s in the column, and the name is a 2. some formulations for the csp problem 41 symbolic identifier for printing the answers. The Cfield of each object points to the column object at the head of the relevant column. Land Rfields of the list headers link together all columns that still need to be covered. This circular list also includes a special column object called the root, h, which serves as a master header for all the active headers. Fields U[h],D[h],C[h],S[h]and N[h]aren’t used. To deepen this with more detail, consider the following 0-1 matrix: A=          0010110 1001001 0110010 1001000 0100001 0001101          It represents an exact cover problem with an exact cover formed by rows {1,4,5}. This matrix would be represented by the objects shown in Figure 2.6, if we name the columns A, B, C, D, E, F and G. The diagram wraps around toroidally at the top, bottom, left, and right. The Clinks are not shown because they would mess up the pic. Every Cfield points to the topmost element in its column. Figure 2.6: Four-way-linked representation of the exact cover problem (borrowed from [13]) The non-deterministic algorithm to find all exact covers can now be cast in the following explicit, deterministic form as a recursive procedure search(k), which is 42 path planning with drones at csp plants invoked initially with k= 0. But, before to see the pseudo-code of the procedure, let’s define some concepts: Operation of printing the current solution. We successively print the rows containing O0, O1, . . . , Ok−1, where the row containing data object Ois printed by printing N[C[O]],N[C[R[O]]],N[C[R[R[O]]]], etc. Operation of choose a column object c. We could simply set c←R[h](leftmost uncovered column). Or if we want to minimize the branching factor, we could select the column with the smallest number of 1’s (we can obtain it by using S fields and iterating jin R[h],R[R[H]],..., while j=h). Operation of covering column c. Removes cfrom the header list and removes all rows in cown list from the other column lists they are in: Algorithm 2: COV ER(c) 1. Set L[R[c]] ←L[c]and R[L[c]] ←R[c]. 2. For each i←D[c], D[D[c]],...,while i=c, 3. for each j←R[i], R[R[i]],...,while j=i, 4. set U[D[j]] ←U[j]and D[U[j]] ←D[j], 5. set S[C[j]] ←S[C[j]] −1. Operation of uncovering column c. Notice that uncovering takes place in precisely the reverse order of covering operation, using the fact that re-adding an element undoes deleting it from the list. We must be careful to unremove the rows from bottom to top, because we removed them from top to bottom. Similarly, it is important to uncover the columns of row rfrom right to left, because we covered them from left to right. Algorithm 3: UNCOV ER(c) 1. For each i←U[c], U[U[c]], . . . , while i=c, 2. for each j←L[i], L[L[i]],...,while j=i, 3. set S[C[j]] ←S[C[j]] + 1. 4. set U[D[j]] ←jand D[U[j]] ←j, 5. Set L[R[c]] ←cand R[L[c]] ←c. We can note this operations as PRINT, CHOOSE, COVER and UNCOVER. Once we have defined them, we can state the search(k)algorithm: 2. some formulations for the csp problem 43 Algorithm 4: search(k) 1. If R[h] = h,PRINT the current solution. 2. Otherwise, CHOOSE a column object c. 3. COVER column c. 4. For each r←D[c], D[D[c]], . . . , while r=c, 5. set Ok←r; 6. for each j←R[r], R[R[r]], . . . , while j=r, 7. COVER column j; 8. search(k+ 1); 9. for each j←L[r], L[L[r]], . . . , while j=r, 10. UNCOVER column j. 11. UNCOVER column cand return. The reader may observe that this pseudo-code is pretty similar to Algorithm X’s one, although it explore all the possible reductions of matrix respect to rows rin a recursive and more efficient way. To understand the printed solutions, we will explain another possible representation for matrix A. Recall that we named its seven columns A to G, so what if we just name rows by the sequence of letters associated with their non-zero columns? If we do this we can get: A=          0010110 1001001 0110010 1001000 0100001 0001101          ⇐⇒          CEF ADG BCF AD BG DEG          . Consider, for example, what happens when search(0) is applied to the data of the matrix Awhich we defined before, as it is represented by the four-way-linked objects in Figure 2.6. Firstly, column A is covered by removing both of its rows from their other columns, the structure now takes the form of Figure 2.7. Continuing search(0), when rpoints to the A element of row (A,D,G), we also cover columns D and G. Figure 2.8 shows the status as we enter search(1), this data structure represents the reduced matrix   B C E F 0 1 1 1 1 1 0 1 ⇐⇒ CEF BCF. 44 path planning with drones at csp plants Figure 2.7: The links after column A in Figure 2.6 has been covered (borrowed from [13]) Now search(1) will cover column B, and there will be no 1’s left in column E. So search(2) will find nothing. Then search(1) will return, having found no solutions, and the state of Figure 2.8 will be restored. The outer level routine, search(0), will proceed to convert Figure 2.8 back to Figure 2.7, and it will advance rto the A element of row (A,D). Now the reduced matrix would be:     B C E F G 0 1 1 1 0 1 1 0 1 0 1 0 0 0 1     ⇐⇒   CEF BCF BG  . After some few steps more, the solution will be found. It will be printed as   AD BG CEF . The solution obtained matches with the subset of rows {1,4,5}that we gave before explaining the search(k)procedure. The main problem we will face when using this approach is that even being the DLX algorithm such an efficient one, our CSP graphs have a lot of nodes and then the collection Sis too large. To reduce at most the size of 2. some formulations for the csp problem 45 subcollection S, we will only consider some certain edges in G, those which are drawn in Figure 2.9. That is, the edges will be those between neighboring nodes, in the sense of being placed immediately next to another node (vertical, horizontal or diagonally). This doesn’t apply to depot node: It will have arcs leaving to and arriving from every other node. Doing this we are relaxing the processing of the ECP approach. However, since n≈40 in the most faced batches pairs, we have find that usually the algorithm can’t solve the ECP in a satisfactory amount of time, because there are many (>1010) different solution covers and it is computationally hard to sort them by its number of tours or total inspection time. Even n≈20 is a high size for solving a CSP instance using ECP, but we should be able to, at least, develop a divide and conquer simple algorithm to chop the faced batches pairs in computationally feasible instances for ECP method. Then we can merge the sub-solutions in order to approximate the global one. Figure 2.8: The links after columns D and G in Figure 2.7 have been covered (borrowed from [13]) 52 path planning with drones at csp plants once. Otherwise, Sis not a solution. This is clearly a simple task which can be done in polynomial time. 3. In addition, we must compute the total weight in each s∈Sand check that is less or equal than W. Otherwise, Sis not a solution. This task is also performed in polynomial time. We define w(s) = w({0, j1, j2, . . . , jr}) = w0j1+wj1j2+· · · +wjr0 as the total energy weight of s∈S. If Spasses all the checkings, then the instance is a YES instance for CSP decision problem. It is easy to see that the asymptotic running time of this certifier algorithm is polynomial. Namely, CSP problem is a NP one. Let’s now see which of the NP-complete problems is easier to polynomial transform to the CSP problem. The classical NP-complete problem that we have chosen here is the Bin Packing decision problem: deciding if Nitems can fit into a specified number kof bins, given an array wjof weights associated to each item j∈ {1, . . . , N}and capacity value cfor the bins. We are following the steps of the recipe now: 1. We have already shown that CSP ∈NP. 2. The NP-complete problem that we use in the proof will be BPP. 3. Now we are going to describe an algorithm to define the function fwhich transforms an instance of BPP into an instance of CSPP. Consider xof BPP, x= (N, k, (wj), c). Let fbe the function which maps x→f(x), which is an instance for CSP decision problem, then f(x) = f(N, k, (w1, . . . , wN), c) =        N+ 1, k,        0w1w2· · · wN w00w2· · · wN w0w10· · · wN . . .. . .. . ..... . . w0w1w2· · · 0        , c +w0        , where w0can be chosen as 1 N N P j=1 wj, for example. Note that the matrix (wij) is formed using w0and the (wj)weights array. This matrix represents the arc costs in a complete digraph G= (V, E, C)where V={0,1, . . . , N}is the node set, E={(i, j): u=v∈V}is the edge set, and C={cij =wj: (i, j)∈ E}is the edge costs set. The computation of the algorithm is almost direct, and is clear that it runs in polynomial time (the most complicated operations 3. complexity and heuristics for the csp problem 53 are computing a mean and reordering (w0, w1, . . . , wN)in a specific way in a (N+ 1) ×(N+ 1) matrix). In Figure 3.2 it is shown how this function works. Figure 3.2: Transforming an example of BPP instance into a CSP one by f 4. The fourth step is about prove that xis a YES instance for BPP if and only if f(x)is a YES instance for CSPP. Suppose, first, that xis a YES instance for BPP, so there exists a collection Bof bins bl={i1, . . . , ir}such that |B| ≤ kand w(bl):=wi1+· · · +wir≤cfor each bl∈B. Now consider f(x) = (N+1, k, (wij), c +w0). For each blwe can construct Cl={0, i1, . . . , ir}, which is a cycle in Gthat starts from zero node and don’t repeat intermediate nodes (as all items in blare different). It’s pretty evident, by Figure 3.3 that w(Cl) = wi1+· · · +wir+w0≤c+w0. Then, each feasible bin produces a feasible cycle. Besides, there the same number of bins that cycles, as we have defined a one-to-one correspondence here, |{Cl}| ≤ k. So f(x)is a YES instance for CSPP. On the other hand, suppose that f(x)=(N+ 1, k, (wij), c +w0)is a YES instance for CSPP. Then there exists a collection of cycles Csuch that |C| ≤ k and for each Cl={0, i1, . . . , ir} ∈ C,w(Cl) = wi1+· · · +wir+w0≤c+w0. Let’s define bl={i1, . . . , ir}. Since w0+w(bl) = w(Cl)≤c+w0, then w(bl)≤cand then blis a feasible bin. As before, we have the same number of bins that cycles, hence |{bl}| ≤ k. Thus, xis a YES instance for BPP. Further explanation can be found in Figure 3.3. 54 path planning with drones at csp plants Figure 3.3: We have that xis a YES instance for BPP if and only if f(x)is a YES instance for CSP. Where xis the example instance from Figure 3.2 5. We have noted yet that the algorithm which defines fruns in polynomial time. This ends the proof. In the next section we will explain how transform classical algorithms and procedures for solving known NP problems to solve CSP instances. 3.2 Algorithms for solving CSP problem In this section we will show the modifications we made to classical algorithms for solving the CSP problems 1 and 2. 3.2.1 Adapting BPP heuristics In the Bin Packing Problem section we described four heuristics that could be adapted to the CSP problem 2. Although our problem is similar to BPP, there is a main difference: CSP weights depend on previous visited node. This affects in two ways: we have to record our last location and also we did not know at first how to sort the items if we would want to apply the sorting versions of these heuristics. For example, let us analyze the following adaptation of Best-Fit algorithm: 3. complexity and heuristics for the csp problem 55 1def bestfitCSP (N,w,W): 2items = list(range (1 ,N +1)) 3bins = [[]]* N 4last_item = 0 5res = 0# Initialize result ( Count of bins ) 6# Create an array to store remaining space in bins 7# there can be at most n bins 8bin_rem = [W]*N 9# Place items one by one 10 while items : 11 # Find the first bin that can accommodate weight [i] 12 j=0 13 # Initialize minimum space left and index of best bin 14 mini = W + 1 15 bi = 0 16 item = items . pop (0) 17 for jin range ( res): 18 last_item = bins [j ][ -1] 19 travel = w[ last_item ][ item ]+w[ item ][0] 20 c1 = bin_rem [j] >= travel 21 c2 = bin_rem [j] - travel < mini 22 if c1 and c2: 23 bi = j 24 mini = bin_rem [j] - travel 25 last_item_bi = last_item 26 # If no bin could accommodate weight [i], create a new bin 27 if (mini == W + 1): 28 bin_rem [ res ] -= w [0][ item ] 29 bins [ res ] = [0 , item ] 30 res += 1 31 last_item = item 32 else:# Assign the item to best bin 33 bin_rem [bi ] -= w[ last_item_bi ][ item ] 34 bins [ bi ]. append ( item ) 35 solution = [b for bin bins if len (b) >0] 36 return solution ,res We first initialize some lists for the N+ 1 items (inspection and depot nodes), the Npossible bins, and the remaining space or energy in each bin. We also initialize the last item variable at 0 (we start from the depot node). Then we enter a loop until items list is empty. At each iteration, we get the index of the next item to assign from the remaining unassigned items. Once we have the next item, we search for the bin with minimum space left and its index. Note that weights depend on last item which is assigned to each bin. If the item can fit in any of the bins, then a new 56 path planning with drones at csp plants one is filled. Otherwise, we assign it to the best bin, those with minimum space left. Also we have added some variables to save and return the solution assignment of items. The adaptations are pretty similar for all the mentioned heuristics and their sorting versions, then we omit here the other codes. In the sorting versions, the main difference is that we select the next item by a given ordering procedure. In this work, we select the next item which returns the smaller weight, given the last added item to each bin. We denote these versions by NFS, FFS, BFS and WFS. 3.2.2 Implementing GBPP with Gurobi In this section we will analyze the implementation of the Generalized Bin Packing Problem formulation with Python and Gurobi. To accomplish this task, we have defined a Python function, SolverBPP_MIP, which takes an input instance and returns the exact solution for GBPP as a MIP problem. 1def SolverBPP_MIP ( instance ): 2# Getting instance relevant parameters 3_,n,w,W,_,_ = instance 4points = [i for iin range (1,n+1)] # points to inspect 5nodes = [0]+ points # nodes set 6 7# Creating Gurobi model 8model = Model (" Generalized_BPP ") 9 10 # Defining variables 11 ind_y_vars = [i for iin nodes ] 12 y = model . addVars ( ind_y_vars , vtype = GRB . BINARY , name = ’y’) 13 ind_a_vars = [(i,j) for iin nodes for jin nodes if i!=j] 14 a = model . addVars ( ind_a_vars , vtype = GRB . BINARY , name = ’a’) 15 ind_x_vars = [(j,k) for jin nodes for kin points] 16 x = model . addVars ( ind_x_vars , vtype = GRB . BINARY , name = ’x’) 17 ind_s_vars = [(i,j,k) for iin nodes for jin nodes for kin points if j!=i] 18 s = model . addVars ( ind_s_vars , vtype = GRB . INTEGER , lb =0 ,ub =2 , name =’sum ’) 19 z = model . addVars ( ind_a_vars , vtype = GRB . INTEGER , lb =1 ,ub =2 , name =’summax’) 20 b = model . addVars ( ind_s_vars , vtype = GRB . BINARY , name = ’bmult ’) 21 ind_L_vars = [i for iin points] 22 L = model . addVars ( ind_L_vars , vtype = GRB . CONTINUOUS , lb =0 , ub =W , name=’L’) 23 3. complexity and heuristics for the csp problem 57 24 # Setting objective function 25 model . setObjective ( quicksum (y[i] for iin nodes ) ,GRB . MINIMIZE ) 26 27 # Defining Constraints 28 model . addConstrs ( quicksum (a[i,j]* w[i,j]* b[i,j,k] for iin nodes for jin nodes if j!=i)<=W*y[k] for kin points) 29 model . addConstrs (b[i,j,k] == x[j,k]* x[i ,k] for iin nodes for jin nodes for kin points if j!=i) 30 model . addConstrs ( quicksum (x[j,k] for kin points ) ==1 for jin points) 31 model . addConstrs ( quicksum (x[j,k] for jin nodes ) == quicksum (a[ i,j]*b[i,j,k] for iin nodes for jin nodes if j!=i) for kin points) 32 model . addConstrs ( quicksum (a[i,j] for iin nodes if j!=i)==1 for jin points) 33 model . addConstrs ( quicksum (a[i,j] for jin nodes if j!=i)==1 for iin points) 34 model . addConstrs (2* a[i ,j] <=z[i ,j] for iin nodes for jin nodes if j!=i) 35 model . addConstr ( quicksum (x[0 ,k] for kin points ) == quicksum (y[ i] for iin points )) 36 model . addConstr ( quicksum (a[i ,0] for iin points ) == quicksum (y[ i] for iin points )) 37 model . addConstr ( quicksum (a[0 ,j] for jin points ) == quicksum (y[ i] for iin points )) 38 model . addConstrs (w[i,j]* a[i ,j] <= L[j] for iin nodes for j in points if j!=i) 39 model . addConstrs (L[i ]+ w[i,j ]* a[i ,j]-W *(1 - a[i,j ]) <= L[j] for iin points for jin points if j!=i) 40 model . addConstrs (s[i,j,k] == x[i,k]+ x[j ,k] for iin nodes for jin nodes for kin points if j!=i) 41 model . addConstrs (z[i,j] == max_ (s[i ,j,k] for kin points) for iin nodes for jin nodes if j!=i) 42 43 # Optimization 44 model . params . outputflag = 0 # Silent mode : 0 45 model . optimize () 46 47 # Extracting solution in a more fancy structure 48 b = model . ObjVal 49 x_sol = np . zeros (( n+1 ,n)) 50 x_list = [] 51 a_sol = np . zeros (( n+1 ,n +1) ) 52 a_list = [] 53 for vin model . getVars () : 58 path planning with drones at csp plants 54 if v. varName [0] in [’s’,’L’,’y’,’b’]: 55 continue 56 if v.xn > 0.9: 57 row ,col = [ int(i) for iin v. VarName [2: -1]. split ( ’,’) ] 58 if v. VarName [0] == ’x’: 59 x_sol [row , col -1] = float (v.xn) 60 x_list . append (( row , col , float (v.xn))) 61 else: 62 a_sol [row , col ] = float (v.xn ) 63 a_list . append (( row , col )) 64 a_list_copy = a_list . copy () 65 sequences = [] 66 while a_list_copy: 67 seq = [] 68 ai = a_list_copy . pop (0) 69 last = ai [1] 70 seq . append ( ai ) 71 while last != 0: 72 al = [aj for aj in a_list_copy if aj [0]== last] 73 if len (al) == 0: 74 print ("FAILED") 75 break 76 ai = al [0] 77 last = ai [1] 78 seq . append ( ai ) 79 a_list_copy = [a for ain a_list_copy if a != ai ] 80 sequences . append ( seq ) 81 solution = [] 82 for seq in sequences : 83 sol = [t [0] for tin seq] 84 solution . append ( sol ) 85 86 return solution We first extract the weights wij, total weight limit Wand the number of inspection nodes n. We use nto define points and nodes sets, which correspond to Vor Kand Vin GBPP formulation. Then we generate the Gurobi model and define the variables xjk, aij , yk, Ljand two auxiliar ones, sijk =xik +xjk ∈ {0,1,2} and zij = maxk∈Kxik +xjk ∈ {1,2}, which we use to implement the constraints 2yij ≤maxk∈Kxik +xjk for i=j= 0, . . . , n. We also define variables bijk which will stand for products xikxjk. After that, we add the objective function and the constraints to our problem. The objective function minimizes the number of initialized bins. The constraints are those described in Chapter 2, in the GBPP section. The 3. complexity and heuristics for the csp problem 59 only difference is the definition of auxiliary variables to compute the products and the maximum sum of xjk pairs. Finally, we launch the optimization and then extract the solution in a more convenient format: we would like to return a covering of feasible tours with minimum size. We express it as a list of tours, for example solution = [[0,1,2,3],[0,4,5,6]] means that we cover the inspection with two tours which are walked in that order. As long as we have seen, Gurobi Optimizer can not solve GBPP instances for N= 40, at least with an academic license. However, we managed to solve N= 20 instances, which is half a faced batches pair. Therefore, we can halve the pair and solve each half as a CSP instance, joining the solutions as a global one. 3.2.3 Implementing GMKP with Gurobi The implementation of the Generalized Multiple Knapsack Problem formulation with Python and Gurobi is pretty similar to the given for GBPP. The main differences with GBPP implementation are the following: It solves CSP problem 1. That is, it maximizes total inspection profit. We define the profit of inspections as negative inspection time costs so that maximizing profits is minimizing time costs. Total profits closer to 0 mean that we spend little total inspection time. 1model . setObjective ( quicksum (y[i,j ]*p[i,j ]*b[i,j ,k] for iin nodes for jin nodes for kin Kif j!= i),GRB . MAXIMIZE ) GMKP takes the number of initialized bins mas a constant. For example: 1model . addConstr ( quicksum (x[0 ,k] for kin K)==m) Since we do not know which mis the best for solving CSP problem 1, we can find the minimum number of feasible tours that is needed for covering the inspection, m∗, using GBPP. Then, we will have that GMKP is feasible for m=m∗, m∗+ 1, . . . , n. Finally, we can perform a search of best the covering in the multiple defined GMKP problems. We usually will find that is enough to only solve instances for m=m∗, m∗+ 1, m∗+ 2 in order to obtain satisfactory solutions. Since we cannot find solutions for GBPP with N= 40 in a satisfactory amount of time, we are not able to solve GMKP for these large instances. However, we could 60 path planning with drones at csp plants solve the problem for N= 20, which is half faced batches pair. In addition, we could consider relaxations for the GMKP, as the neighboring nodes one, which states that the only allowed edges are the associated with neighboring nodes and the associated with the depot node. This relaxation allows us to solve instances until N= 40 nodes. 3.2.4 Adapting ECP algorithms In this section we will explain how we combine Johnson’s algorithm and Algorithm X to solve the CSP problems (versions 1 and 2). The code is the following: 1import numpy as np 2from algorithm_x import AlgorithmX 3import graph_utils as gu 4from graph_utils import tic ,toc 5def SolverECP ( instance , cycle_limit_length =13 , 6maxit = -1, timeLimit = -1 , stop_when_imp =0, plot = True ): 7# Getting instance 8G,N,w,W,p, pos2 = instance 9print (’\n’+’-’*30+’ ECP SOLVER ’+’-’*30) 10 11 # Getting cycles . Applying Johnson ’s algorithm 12 limit = cycle_limit_length 13 C = list( gu . limit_simple_cycles (G, limit , start =0) ) 14 Cw = [c for cin Cif gu . GetCycleWeight (c,w) <=W] 15 16 # Solving ECP. Applying AlgorithmX with DLX implementation 17 C0 = [c [1:] for cin Cw ] 18 C0. sort (key = lambda x:-len (x)) 19 solver = AlgorithmX (N) 20 for i,c0 in enumerate (C0): 21 solver . appendRow (np . array ( c0 ) -1, i+1) 22 solutions = solver . solve () 23 24 # Getting optimal covers 25 tStart = tic () 26 i0 = 0 27 C0_arr = np . array (C0 , dtype =’object’) 28 min_len_sols = [] 29 min_len = N 30 max_profit_sols = [] 31 max_profit = -1e10 32 n_covers = 0 33 tsol = gu. GetcurrentSolution (N, W, w) 3. complexity and heuristics for the csp problem 61 34 tsol_profit = gu. GetCoverProfit ( tsol , p) 35 tStartGetSols = tic () 36 print (" Starting getting covers process ... ( feasible cycles : {}) "\ 37 .format(len (Cw))) 38 for solution in solutions : 39 cover_rows = [s -1 for sin solution] 40 cover = C0_arr [ cover_rows ] 41 current_profit = gu.GetCoverProfit(cover,p) 42 current_len = len( cover ) 43 n_covers += 1 44 if current_len < min_len : 45 min_len = current_len 46 min_len_sols = [ cover ] 47 elif current_len == min_len : 48 min_len_sols . append ( cover ) 49 if current_profit > max_profit : 50 max_profit = current_profit 51 max_profit_sols = [ cover ] 52 elif current_profit == max_profit : 53 max_profit_sols . append ( cover ) 54 # Terminate because of max iteration exceeded 55 i0 += 1 56 if maxit > 0: 57 if i0 >= maxit : 58 print ("\ nMAXIT_Warning : there are more than {} covers!"\ 59 .format( maxit )) 60 break 61 # Terminate because of time limit exceeded 62 if timeLimit >0: 63 if toc ( tStartGetSols )>timeLimit : 64 print ("\ nTIMELIM_Warning : processing has exceeded {}s of processing "\ 65 .format( timeLimit )) 66 break 67 # Terminate because of max iteration exceeded 68 if stop_when_imp > 0: 69 current_imp = (tsol_profit - current_profit )/ tsol_profit 70 if current_imp >= stop_when_imp : 71 print ("\ nReached an improvement of {:.2 f}% over current solution "\ 72 .format( current_imp *100) ) 73 break 74 # Stats 68 path planning with drones at csp plants Figure 4.2: Total profit of current company solution with N= 20:4004s s ≈67 min Figure 4.3: Total profit of current company solution with N= 40:7978 s≈132 min 4.1 BPP heuristics results The results of the BPP-based heuristics are shown in this sub-section. This is the fastest approach for solving CSP instances due to its low computational cost. For example, an instance with N= 40 can be resolved in 0.01 swith the heuristic BPP 4. results 69 solver that we developed, which applies the four heuristics presented in this work and the sorting versions. However, we must keep in mind that this procedure solves CSP problem 2, although it may return good feasible solutions for CSP problem 1. For obtaining this results, we have implemented the adapted algorithms in Python and also have included them into a more fancy style. A typical output for the N= 10 instance presented before is the following: 1-------------- BPP ( Heuristic ) SOLVER ------------- 2Ellapsed ( global ): 0.0 3Current solution : Profit : -2007.00 , N. bins : 2 4NF: Profit : -2007.00 , Improvement : -0.00% , N. bins : 2 5NFS : Profit : -1565.00 , Improvement : 22.02% , N. bins : 1 6FF: Profit : -2007.00 , Improvement : -0.00% , N. bins : 2 7FFS : Profit : -1565.00 , Improvement : 22.02% , N. bins : 1 8BF: Profit : -2007.00 , Improvement : -0.00% , N. bins : 2 9BFS : Profit : -1565.00 , Improvement : 22.02% , N. bins : 1 10 WF: Profit : -2007.00 , Improvement : -0.00% , N. bins : 2 11 WFS : Profit : -1565.00 , Improvement : 22.02% , N. bins : 1 12 13 ---------------- END OF SOLVER ----------------- The total elapsed time during the processing of the 4 heuristics in its two versions (normal and sorting) is displayed in line 1. Line 2 shows the profit and number of used bins (tours) for the current company solution. Next lines show the total profits (negative total inspection times), the improvement over the current company solution and number of used bins obtained, by each heuristic: Next-Fit (sorting), First- Fit (sorting), Best-Fit (sorting) and Worst-Fit (sorting). As we have highlighted, BPP heuristics solve CSP problem 2, that is, the minimization total tours problem. For that, we expected to rarely find good enough solutions for CSP problem 1 (minimizing total inspection time cost). However, as we can observe in the output for the N= 10 instance, we managed to improve the current solution in a 22.02% with the sorting versions of these heuristics. This saves 442 seconds, which is more than 7 minutes! CSP inspections are developed in all pipes of the plant (≈160−180 inspection nodes, in small CSP plants). So, if we divided the CSP plant in groups of ten nodes and apply the two strategies, we would obtain an overall improvement of approximately two hours! In Figures 4.4 and 4.5 we represent the two solutions obtained by these approaches (the FFS, BFS and WFS heuristics return the same covering). Also, we have surprisingly found that good enough solutions for the N= 40 instance are given by these procedures. Figures 4.7 and 4.6 show the solutions obtained (in this case the BFS, FFS and WFS also concur). 70 path planning with drones at csp plants Figure 4.4: Total profit of NF, FF, BF and WF solutions with N= 10:2007 s Figure 4.5: Total profit of FFS, BFS and WFS solutions with N= 10:1565 s 4. results 71 Figure 4.6: Total profit of FF and BF solutions with N= 40:7944 s Figure 4.7: Total profit of NFS, FFS, BFS and WFS solutions with N= 40:6466 s 4.2 GBPP results The Generalized Bin Packing Problem that we have developed in this document aims to solve the CSP problem 2. As its purpose is not minimizing the total inspection time, the obtained solutions include many edge crossings that likely delay the total inspection time costs. In fact, we use this algorithm to feed the Generalized Multiple 72 path planning with drones at csp plants Knapsack Problem, which needs to specify the number of used bins. If we combine the two procedures, we can obtain very good solutions for CSPP1, even exact ones for N≤20 instances. However, when Nis bigger, for example 30 or 40, the problem becomes computationally infeasible, spending lot of time without finding solutions, even if we apply relaxations like the neighboring nodes ones or if we configure the Gurobi Optimizer to find at least approximate solutions. A typical output for a N= 10 instance is the following: 1--------------- BPP (MIP ) SOLVER ------------ 2Ellapsed ( global ): 2.6068475246429443 3RESULTS SUMMARY ( min length criteria , nsols : 2) : 4( Retrieving solution with maximum associated profit ) 5Items in each bin: 6Bin 1: [0, 3, 5, 9, 10 , 6, 1, 7, 4, 8, 2] 7 8Current company solution : Profit : -2007.00 , N. bins : 2 9BPP - MIP : Profit : -1761.00 , Improvement : 12.26% , N. bins : 1 10 11 --------------- END OF SOLVER --------------- Figure 4.8: Total profit of GBPP exact solution with N= 10:1761 s Line 2 shows elapsed time solving the CSP instance. The next lines summary the results: the solver finds 2 solutions with minimum number of tours, and it selects that with minimum inspection time cost (maximum profit). After that it comes the solution and the summary of total profit, improvement over current company solution and number of bins. We observe that this instance returns a solution worse than the 4. results 73 given ones by the previous heuristics. Figures 4.8 and 4.9 represent the exact solution of GBPP for the N= 10 and N= 20 instances. Figure 4.9: Total profit of GBPP exact solution. with N= 20:3517 s 4.3 GMKP results The GMKP problem is responsible for solving the CSP problem 1, that is, the total inspection time minimization problem. In order to relax this formulation when the size of the instance is big (N= 30 or N= 40), we can consider the hypothesis on the neighboring nodes and reduce the number of variables and constraints. Until N≈20 we can obtain exact solutions by this procedure. Recall that this formulation requires as input the number of used bins m. Our idea is to solve the CSP problem 2 by GBPP, obtaining the minimum number of feasible tours m∗needed for cover the graph and use it as input. Although technically (we saw an example) the best solution in terms of CSP problem 1 can use more tours than best solution for CSP problem 2, we have observed that this rarely occurs, probably due to the symmetry of the graphs. A typical output of this procedure for the N= 10 CSP instance is shown below. 1------------------ MKP SOLVER ----------------------- 2Ellapsed ( global ): 0.1405172348022461 3RESULTS SUMMARY ( max profit criteria , nsols : 9) : 4Items in each bin: 5Bin 1: [0, 8, 3, 7, 2, 6, 1, 10 , 5, 9, 4] 6 74 path planning with drones at csp plants 7Current company solution : Profit : -2007.00 , N. bins : 2 8GMKP : Profit : -1564.00 , Improvement : 22.07% , N. bins : 1 9 10 ------------------- END OF SOLVER ------------------ Line 2 shows the elapsed time in the processing. The next lines show the results summary: the solver has found 9 optimal solutions, then it is shown the selected solution and the comparison between company solution profits and GMKP ones. The improvement ascend to 22.07%, the inspection of ten nodes lasts 1564 seconds. Figures 4.10 and 4.11 represent the optimal solutions of CSP problem 1 for instances with N= 10 and N= 20. We see that the solutions present a lot of symmetry and that mainly use edges associated with neighboring nodes. We also observe that the solutions avoid horizontal edges, which imply a longer turning phase than other ones, and focuses in use vertical and diagonal arcs. These patterns show that maybe a possible heuristic to develop could be one which smartly chooses between diagonal and vertical edges composing a good covering of feasible tours. As we commented before, the GBPP approach can not solve instances for N > 20 and therefore we do not know exactly the value of m∗, in order to use it as input for GMKP. For large instances we use the BPP heuristics, which usually give one or two more bins that optimal solutions, at least for the used instances. That information gives an estimation of a good mfor feeding the GMKP, and we can also try use one or two less bins that the returned by the heuristics to see if it is feasible. Another interesting point in the GMKP is the neighboring relaxation: as commented above, if we only consider edges associated with neighboring nodes, we can relax the formulation and solve instances for N≥30. We call this the relaxed GMKP. This relaxation can also apply for the GBPP but we have found that it does not increase the solver performance very much. The exact solution for the relaxed GMKP with N= 30 is shown in Figure 4.12. We also have even found an exact solution for a N= 40 CSP instance using relaxed GMKP displayed in Figure 4.13. The gap between the best lower bound found and the objective value was of <2% after 49 minutes of processing. After approximately 24 hours and 30 minutes, the gap was of 0.18%, with a best lower bound of 6448.5s. The processing stopped with an elapsed time of 25 hours and 46 minutes. So, the final exact solution of 6460 sis the optimal for the relaxed GMKP in the N= 40 instance and it would return a saving over the current company solution of 25 minutes and 20 seconds. 4. results 75 Figure 4.10: Total profit of GMKP exact solution with N= 10:1564 s Figure 4.11: Total profit of GMKP exact solution with N= 20:3158 s 76 path planning with drones at csp plants Figure 4.12: Total profit of relaxed GMKP exact solution with N= 30:4837 s Figure 4.13: Total profit of relaxed GMKP approximate solution with N= 40:6460 s 4.4 ECP results Among all the procedures described in this document, Exact Cover Problem one is the most expensive from the computational point of view. We even can not compute the exact solution for the N= 10 instance. However, we have managed to obtain good enough solutions aggregating different relaxations: 4. results 77 Instead of considering the complete digraph, we only keep the edges associated with neighboring edges and the ones who joint the depot node with the inspection nodes. This obviously leads to sub-optimal solutions but they are usually good enough solutions. We add some termination criteria to the algorithm, in the part that searches for optimal solutions for CSP problems 1 and 2. These criteria are three: maximum number of iterations of the search, time spent in the search and improvement over the current company solution. Combining the relaxation and the termination criteria we could obtain solutions even for a N= 40 CSP instance. The output for the N= 10 instance is shown below: 1------------------- ECP SOLVER --------------------- 2Starting getting covers process ... 3ECP (current ,best ) - iter : 4290000 , profit : ( -1942.00 , -1593.00) , len: (4 ,1) 5Ellapsed ( global ): 155.21034479141235 , CycleLengthLimit : 13 , FeasibleCycles : 15798 , Covers : 299613 6RESULTS SUMMARY ( min length criteria , nsols : 1506) : 7Items in each bin: 8Bin 1: [0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5] Ordering : [0 , 6, 1, 7, 2, 8, 3, 9, 4, 10, 5] 9 10 RESULTS SUMMARY ( max profit criteria , nsols : 2) : 11 Items in each bin: 12 Bin 1: [0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5] Ordering : [0 , 6, 1, 7, 2, 8, 3, 9, 4, 10, 5] 13 14 Current company solution : Profit : -2007.00 , N. bins : 2 15 ECP ( min length criteria ): Profit : -1593.00 , Improvement : 20.63% , N.bins: 1 16 ECP ( max profit criteria ): Profit : -1593.00 , Improvement : 20.63% , N.bins: 1 17 18 ----------------- END OF SOLVER ----------------- First line shows a message that helps to know when the searching part of the algorithm starts. After that we added a dynamic printing that shows the best solution found and the current solution at the iteration (this message updates after some iterations). The next lines show the counting of feasible cycles and covers found, and also the selected cycle length limit. As we can see, the cycle limit length is loose enough because there are no feasible cycles with that length. During the searching process 84 path planning with drones at csp plants the pilot in critical environments. We considered two problems related with the inspections of CSP plants, namely the CSP problem 1 and the CSP problem 2. The former asks for a minimization in the inspection total time while the latter proposes minimizing the number of tours. We have developed a series of procedures, which are based in classical optimization problems, we have adapted and implemented them in Python programming in order to solve CSP instances, and compared the results obtained with the current solution adopted by the company Virtualmech. The CSP instances where we evaluated the procedures were configured to be as realistic as possible in such a way the solutions obtained will be suitable estimations for real scenarios. The results show improvements over the current solution of almost 20% for instances with 40 inspection points. This supports the use of mathematics instead of only good intuition. In addition, during the work we also have proven that the CSP problem 2 is NP- complete. We think that the CSP problem 1 could also be proven to be NP-complete, if we use the MKP in the reduction. Finally, although we have made some advances in the study of path planning with drones at CSP plants, there are still some loose ends, which could give ideas for future research in this area: It is necessary to model the real energy consumption weights of the currently used drone. The company uses a DJI Matrice 300RTK1drone model. We have tried several times to extract the kinematics information from log files generated during the inspections but we have not succeed. It would be very useful for the continuity of the work to manage extract this information and then use it to compute the beta coefficients of the standard regression model explained in the Chapter 2. The formulations developed to solve the CSP problems 1 and 2, by exact procedures, GMKP and GBPP, respectively, are very complex formulations that use many variables and constraints. This makes difficult to solve big sized instances exactly. Developing new smaller formulations may be helpful to reduce the computational cost seems possible. Another research line could be to improve our control over Gurobi Optimizer. It includes several features and functionalities. A well-chosen configuration could increase the performance of the exact procedures and reduce the computational cost. 1https://www.dji.com/es/matrice-300/specs 5. conclusions 85 We can also experimented with the neighboring relaxations. For example, we have restricted the edges of the graph to be neighbors no further than 1 in the xcoordinates, but we could increase that threshold and compare both the improvement and the elapsed processing time with the original neighboring relaxation. We observed that the obtained showed some particular patterns. These patterns could be divided into: using a zigzag strategy and using a diagonal strategy. That is, we start from a point and move only using vertical or diagonal edges associated with neighboring nodes. We could use this information as an idea for developing a smart heuristic algorithm which sequentially decides between using some types of edges to cover the graph of the CSP instance. Finally, in real situations, the edge weights are changing during a working day, for example, with wind changes. Therefore, there is a clear need in the inspection of CSP plants for online route planning algorithms that use incomplete information that is acquired in a dynamic manner. A Kinematics for energy consumption weights and time costs In this appendix, the aim is to compute realistic values for energy consumption weights wij. We use basic classical kinematics equations to model the kinematics state, sof a drone travelling during the CSP inspection. The given inputs are: vmax xy m/s, the maximum horizontal speed (in module) reached by the drone. vmax zm/s, maximum speed module while ascending. v−max zm/s, maximum speed module while descending. aa xy,h m/s2, horizontal acceleration module of the drone when it’s moving horizontally, during the acceleration part. ab xy,h m/s2, horizontal acceleration module of the drone when it’s moving horizontally, during the braking part. az,h m/s2, vertical acceleration component of the drone when it’s moving horizontally. It remains constant in the three parts, keeping the drone at the same z coordinate. aa z,t m/s2, vertical acceleration component of the drone when it’s taking off, during acceleration part. ab z,t m/s2, vertical acceleration component of the drone when it’s taking off, during braking part. aa z,l m/s2, vertical acceleration component of the drone when it’s landing, during acceleration part. ab z,l m/s2, vertical acceleration component of the drone when it’s landing, during braking part. d(i, j)m, distance between node iand node j. vmaxr xy m/s, maximum speed module when returning to base. amaxr xy m/s2, maximum acceleration module when returning to base. 88 path planning with drones at csp plants ztarget m, constant vertical distance between the take-off zone and the inspection points. g= [0,0,−9.80665] m/s2, gravity acceleration vector. ωmax rad/s, angular velocity during turnings. It remains constant. m kg, constant payload. This value is senseless for us because we don’t have a load. ωxy m/s, horizontal proyection of wind speed. A.1 Takeoffs Takeoff kinematics can be modelled using equations for uniform and accelerated rectilinear motion. Suppose the drone is on the ground and we want to elevate it ztarget meters. It firstly starts the vertical motion by accelerating at aa z,t, until it reaches vmax zand begins to move uniformly. At the end it brakes at ab z,t and arrives at ztarget with vz= 0. In Figure A.1 we can observe the behaviour of z components for acceleration, speed and position vectors of a drone during a takeoff with inputs: ztarget = 20 m, aa z,t = 10 m/s2, ab z,t =−10 m/s2, vmax z= 4.5m/s. Figure A.1: Drone kinematics during a takeoff. The magnitudes that don’t appear are supposed to remain constant. a. kinematics for energy consumption weights and time costs 89 A.2 Landings Suppose the drone is on the air and we want it to land from zstart meters. It firstly starts the motion by descending at aa z,l, until it reaches v−max zand begins to move uniformly. At the end it brakes at ab z,l and arrives to ground with vz= 0. In Figure A.2 we can observe the behaviour of z components for acceleration, speed and position vectors of a drone during a landing with inputs: ztarget = 20 m, aa z,l = 8 m/s2, ab z,l = −8m/s2, vmin z= 3.5m/s. Figure A.2: Drone kinematics during a landing. The magnitudes that don’t appear are supposed to remain constant. A.3 Horizontal displacements Suppose the drone is on the air and we want it to travel in straight line between two inspection nodes, iand j, or more generally, between two locations with coordinates [x0, y0, z0]and [x1, y1, z1]. At first the drone is hovering at (x, y, z) = (x0, y0, z0). We suppose it is already heading its target. Then it accelerates at aa xy,h until reach vmax xy and performs a uniform motion. Before it arrives at (x, y, z) = (x1, y1, z1)it performs the braking part and decelerates at ab xy,h. At target, vxy = 0. During this displacement, the z components of acceleration, speed and coordinate remain constant. In Figure A.3 we can observe the behaviour of x components for acceleration, speed and position vectors of a drone during a horizontal displacement 90 path planning with drones at csp plants between locations (0,0,20) mand (10,10,20) mwith inputs: vmax xy = 4 m/s, amax xy = 10 m/s2. Figure A.3: Drone kinematics for x components during a horizontal displacement. The kinematics for y components is pretty similar and the z components don’t vary. A.4 Turnings Turnings are modelled as simple uniform rotation at a given constant angular velocity ωmax. We suppose that the drone always turns at the same velocity, it doesn’t perform neither acceleration nor braking part. We also consider that it always choose the orientation of turning with the smallest arc length, respect to the initial and target angles. The angles are defined by convention or using the director vector associated with the displacement from ito j. For example, we suppose that the drone starts at base station with an initial angle of π/2rad respect the OX axis. The kinematics of this part is computed using simple equations: we know the initial and target angle and the angular velocity, so we know how much time the process lasts. A.5 Weights and costs computation The process is then simple, knowing the kinematics we know how much the action last and which are the components for accelerations and velocities. Then we can a. kinematics for energy consumption weights and time costs 91 use the firsts to compute inspection time costs: tij is the last value of the time array associated with kinematics of going from ito j. On the other hand, the acceleration and velocities components are used in the standard regression model showed at Chapter 2to obtain the wij values, i.e, the energy consumption weights of going from ito j. Finally, in order to evaluate the standard regression model, we need some beta coefficients values for that. In [4] and [10] the authors present 3 vectors standing for beta coefficients obtained from 3 different drone models. In this work we have been working with these values, specifically with 3DR-Solo ones, in order to evaluate the model and obtain realistic weights. This model was chosen because weights behave qualitatively pretty similar to the energy consumption behaviour of the drone used at real CSP inspections by the company. Figure A.4: Realistic energy consumption weights computed for an inspection graph with 10 nodes. Wind horizontal vector: ωxy = (1,0). Beta coefficients model: 3DR- Solo 92 path planning with drones at csp plants Figure A.5: Realistic inspection time costs computed for an inspection graph with 10 nodes Bibliography [1] H.L. Zhang, J. Baeyens, J. Degrève, G. Cacères. (2013). Concentrated solar power plants: Review and design methodology. Renewable and Sustainable Energy Reviews, 22, 2013, Pages 466-481 [2] C. E. Miller, A. W. Tucker, and R. A. Zemlin. 1960. Integer Programming Formulation of Traveling Salesman Problems. J. ACM 7, 4 (Oct. 1960), 326–329. DOI:https://doi.org/10.1145/321043.321046 [3] Borcinova, Zuzana. (2017). Two models of the capacitated vehicle routing problem. Croatian Operational Research Review. 8. 463-469. 10.17535/crorr.2017.0029. [4] Tseng, C. M., Chau, C. K., Elbassioni, K. M., & Khonji, M. (2017). Flight tour planning with recharging optimization for battery-operated autonomous drones. CoRR, abs/1703.10049. [5] Campbell, J. F., Corberán, Á., Plana, I., & Sanchis, J. M. (2018). Drone arc routing problems. Networks, 72(4), 543-559. [6] Mesas-Carrascosa, F. J., Verdú Santano, D., Pérez Porras, F., Meroño-Larriva, J. E., & García-Ferrer, A. (2017). The development of an open hardware and software system onboard unmanned aerial vehicles to monitor concentrated solar power plants. Sensors, 17(6), 1329. [7] Jose, K., & Pratihar, D. K. (2016). Task allocation and collision-free path planning of centralized multi-robots system for industrial plant inspection using heuristic methods. Robotics and Autonomous Systems, 80, 34-42. [8] Amarat, S. B., & Zong, P. (2019). 3D path planning, routing algorithms and routing protocols for unmanned air vehicles: a review. Aircraft Engineering and Aerospace Technology.