0 1 4 10 0 8 2 7 3 0 6 9 5 0

**Figure 1.** An example of a solution representation for a CARP instance with 10 required tasks, where 0 is the *id* of the dummy task. In this example, there are 3 routes. The first route services the tasks with *id*s 1, 4 and 10. The second services the tasks with *id*s 8, 2, 7 and 3. The third services the tasks with *id*s 6, 9 and 5.

#### 3.1.2. Objective and Constraints

The objective of the CARP is to minimize the total cost of the solution *S* subject to some constraints, which are defined in this section. The total cost of a solution *S* (i.e., *TC*(*S*)) is calculated with the following formula (Equations (4)–(6)):

$$TC(S) = \sum\_{k=1}^{|S|} DC(r\_k) + SC(r\_k) \tag{4}$$

$$\text{pDC}(r\_k) = \text{mdc}\left(t\_0, \text{head}(t\_{k,1})\right) + \sum\_{i=1}^{l\_k - 1} \text{mdc}\left(\text{tail}(t\_{k,i}), \text{head}(t\_{k,i+1})\right) + \text{mdc}\left(\text{tail}(t\_{k,l\_k}), t\_0\right) \tag{5}$$

$$SC(r\_k) = \sum\_{i=1}^{l\_k} sc(t\_{k,i}) \tag{6}$$

where *DC*(*rk*) and *SC*(*rk*) are the total dead-heading and service cost of the route plan *rk*.

The solution *S* has to satisfy the following constraints. First, each route plan starts and ends at the depot. Second, each task is executed exactly once. Therefore, the total number of tasks executed on each route plan (excluding the dummy task *t*0) is equal to *n*:

$$\sum\_{k=1}^{|S|} l\_k = n$$

Moreover, a task cannot be executed more than once, neither in the same route nor in another route:

$$(t\_{a,i} \neq t\_{b,j}) \forall (a,i) \neq (b,j) \tag{8}$$

where *ra* and *rb* are route plans within *S*, *ta*,*<sup>i</sup>* is the *i*-th task in the route plan *ra*, and *tb*,*<sup>j</sup>* is the *j*-th task in the route plan *rb*. If a task *t* has an inverse (i.e., ∃*inv*(*t*)), then either *t* or *inv*(*t*) is executed. Both cannot be executed in the same solution. Third, the total demand served each route plan does not exceed the capacity (*q*) of the vehicle:

$$\sum\_{i=1}^{l\_k} dem(t\_{k,i}) \le q\_\prime \qquad \forall k \in \{1, 2, \dots, |S|\} \tag{9}$$

#### *3.2. The Data-driven DCARP*

There are various approaches for DCARP, but in this work, the data-driven version of DCARP is considered, which was recently formulated in [20].

In this problem, instead of one static CARP instance, there is a series of DCARP instances (i.e., a DCARP scenario [21]) that needs to be solved. A DCARP scenario is denoted by I = *I*0, *I*1, ... , *Ii*, ... , *Im*−1, where *m* is the number of DCARP instances within the scenario (i.e., the number of dynamic events that occurred and changed the previous DCARP instance is *m* − 1). Each *Ii* (*Ii* ∈ I) DCARP instance contains all the information about the current problem. The previous DCARP instance *Ii*−1, the execution of the accepted solution for *Ii*−<sup>1</sup> and the occurred event(s) define the next DCARP instance *Ii*, where 0 < *i* < *m*. The initial instance (*I*0) can be viewed as a static (data-driven) CARP instance, since initially every vehicle is in the depot (in good state) and no task has been executed yet.

For a data-driven DCARP instance, information needs to be stored about all the vehicles and route plans. For each vehicle, the current location and state have to be known. Furthermore, identifiers are needed to be used for the vehicles and the route plans, since a vehicle may follow multiple route plans (one after another) and it is important to know for each route plan whether a vehicle already executed it, its execution is still in progress, or a vehicle still needed to be assigned to it to start its execution.

Instead of the number of vehicles (*w*), a set of identifiers of all the vehicles is needed to be defined, which is denoted by *H*. The set of the identifiers of the (currently) free vehicles is denoted by *Hf* (*Hf* ⊆ *H*), which is initially equal to *H*. The identifier of a vehicle is added to *Hf* , if the vehicle finishes the execution of a route plan, and the identifier is removed, when a new route plan is assigned to the vehicle. If the execution of all the tasks is finished and all the vehicles are returned to the depot (i.e., there are no broken down vehicles outside on the roads), then *Hf* = *H*, otherwise *Hf* ⊂ *H*.

The set of identifiers of all the route plans is denoted by *R*, and the set of identifiers of the route plans that currently cannot be modified and not executed by any vehicle is denoted by *Re* (*Re* ⊆ *R*). When a new route plan is created, its identifier is added to *R*, and when the execution of it is finished or suspended (due to vehicle breakdown), its identifier is added to *Re*. If the execution of all the route plans is finished and there are no more tasks to execute, then *Re* = *R*, otherwise *Re* ⊂ *R*. The identifier is removed from *Re* only if the vehicle which is assigned to it was broken, but got fixed and can continue the execution of the plan. The function that defines which vehicle is assigned to a specific route plan is denoted by *rv* : *R* → *H*.

To store the current location of the vehicles in the instance, the virtual task strategy introduced in [21] is used, which replaces the executed tasks in each route plan with "virtual tasks". A "virtual task" is an arc whose *head* is the depot vertex *v*<sup>0</sup> and *tail* is the current location of the vehicle, vertex *v* (*v* ∈ *V*). For the sake of simplicity, it is assumed that when an unexpected event occurs, every vehicle is located exactly at a vertex. Since this task is "virtual", it cannot be traversed, for this reason, it has an infinite traversing cost (i.e., *dc*(*v*) = ∞). Furthermore, since it is a "task", there is a demand and a service cost assigned to it, which are calculated according to the provided data: the service cost is the total cost produced by the vehicle so far (i.e., it is the sum of traversing and serving cost of the arcs that were crossed or served by the vehicle) , and the demand is the total demand served by the vehicle so far. A route plan can have at most one virtual task. Therefore, if a route plan already has a virtual task, then it is updated taking into account the arcs traversed and the tasks executed since then by the corresponding vehicle. The set of all virtual tasks is denoted by *Tv* (*Tv* ⊆ *T*), and the function that defines which virtual task belongs to a specific route plan is denoted by *rt* : *R* → *Tv*.

The set of arc tasks that need to be executed is denoted by *T*. If according to the gathered information a task *t* (*t* ∈ *T*) was executed by the vehicle *h* (*h* ∈ *H*), then in the new DCARP instance *t* needs to be removed from *T*. Furthermore, the virtual task of the route plan of the vehicle (e.g., *tk*,*v*, where *rv*(*k*) = *h* and *tk*,*<sup>v</sup>* ∈ *Tv*) needs to be updated. The new virtual task is generated in such a way that *t* is included in it (along with the other tasks the vehicle executed and arcs the vehicle traversed). If *t* has an inverse (i.e., ∃*inv*(*t*)), then it is removed from *T* as well. Accordingly, the total number of tasks that have to be executed (*n*) is decreased by one or two.

The initial DCARP instance (*I*0) is similar to a static CARP instance. The sets of the route plan identifiers (*R*) and the function *rv* are created and filled only after the solution is found for *I*0. The set *Hf* is initially equal to *H*, then based on *rv*, all the vehicle identifiers that are assigned to a route plan are removed from *Hf* . At this stage, the sets *Re* and *Tv* are empty sets, therefore the function *rt* is an empty function as well. According to these, the initial DCARP instance (*I*0) is defined as follows:

$$\mathbf{l}\_0 = (V, \upsilon\_0, \mathbf{A}, T, n, w, q, H, head, tail, dc, inv, dem, sc, m dc) \tag{10}$$

The subsequent DCARP instances (*Ii*, where 0 < *i* < *m*) are defined as follows:

$$I\_i = (V, v\_0, A, T, T\_\upsilon, n, q, H, H\_f, R\_\tau R\_\varepsilon, \text{rt}, \text{r}\upsilon, \text{head}, \text{tail}, \text{clc}, \text{inv}, \text{dem}, \text{sc}, \text{undc}) \tag{11}$$

#### 3.2.1. Structure of a Scenario

A new DCARP instance is constructed and added to the DCARP scenario, when an unexpected event happens that changes the current problem to such an extent that it has effect on the currently executed solution. In [21] all the possible events (based on realistic assumptions) were collected and analyzed based on their effect.

It is assumed that the roadmap, the number of vehicles, and the maximum capacity of the vehicles cannot change (at least during the execution of the solution). Therefore, *V*, *v*0, *A*, *head*, *tail*, *inv*, *q* and *H* are the same in all the DCARP instances of a DCARP scenario.

It is assumed that roads can become closed/opened (it changes *dc*, thus *mdc*, too), the traffic can decrease/increase (it changes *dc* and in some cases *sc*, thus *mdc*, too), tasks can get cancelled/added (it changes *T*, *n*, *dem*, and *sc*) and vehicles can breakdown/restart (it changes *Re*), which are unexpected events. The expected events are the events that normally occur during the execution of the solution: a task is executed (it changes *T*), a vehicle moves (it changes *Tv*, thus *rt*, too), or a vehicle returns to the depot (it changes *Hf* and in some cases *rv*). The affected components are updated only when a new instance is constructed. If rerouting is performed, then *R* and *rv* may change, but the changes are visible only in the next DCARP instance. Therefore, *T*, *Tv*, *n*, *Hf* , *R*, *Re*, *rt*, *rv*, *dem*, *sc*, *dc*, and *mdc* may be different among the DCARP instances of a DCARP scenario.

Since due to the unexpected events some components of the DCARP instance change, the optimal solution may change, too. It is one's choice to construct a new DCARP instance and reroute when there might be a better solution available, but the current solution is still feasible. However, constructing a new DCARP instance and rerouting is necessary when the current solution is not feasible anymore.

#### 3.2.2. Solution Representation

For each DCARP instance, the solution representation is mainly the same as for static CARP instances. The only difference is that if the route plan has a virtual task assigned to it, then the virtual task is the second task within the route plan (since the first task is always the dummy task *t*0). For instance, if the route plan *rk* = *t*0, *tk*,1, *tk*,2, ... , *tk*,*lk* , *t*0 has an identifier *k* (*k* ∈ *R*) and there is a virtual task *tv*,*<sup>k</sup>* assigned to it (i.e., *rt*(*k*) = *tv*,*k*), then *tk*,1 = *tv*,*k*.

#### 3.2.3. Objective and Constraints

For each DCARP instance, the objective and the constraints are mainly the same, as well as for static CARP instances. The only difference is at the second constraint, which requires that the total number of tasks in the solution *S* (excluding the dummy task *t*0) is equal to the sum of the number of tasks that still need to be executed (*n*) and the total number of virtual tasks (|*Tv*|):

$$\sum\_{k=1}^{|S|} l\_k = n + |T\_v| \tag{12}$$

The attributes of a virtual task guarantee that a solver will always place the virtual task right after the dummy task within a route plan of a (nearly optimal) solution, so there is no need to add a constraint regarding it.

#### 3.2.4. Finding a Solution

The data-driven DCARP framework allows rerouting when a critical event (i.e., an unexpected event that may change the feasibility of the current solution) occurs. These events are the task appearance, the demand increased and the vehicle breakdown.

The data-driven DCARP framework allows the use of static CARP solvers by converting the current data-driven DCARP instance into a static CARP instance. After a (sufficiently good) solution is found by the CARP solver, the solution is converted into a data-driven DCARP solution.

Converting a data-driven DCARP instance into a static CARP instance works as follows: the sets of vehicle and route plan identifiers (i.e., *H*, *Hf* , *R*, and *Re*) are omitted, along with the related functions (*rt* and *rv*). Furthermore, all virtual tasks related to finished and suspended route plans are removed from *T*. For instance, if *rt*(*k*) = *tv*,*<sup>k</sup>* (*tv*,*<sup>k</sup>* ∈ *T*) is the virtual task of the route plan with identifier *k* (*k* ∈ *R*) and *k* ∈ *Re*, then *tv*,*<sup>k</sup>* is removed from *T* (i.e., *T* \ {*tv*,*k*}).

Converting a static CARP solution into a data-driven DCARP solution works as follows: the virtual tasks that are related to finished and suspended route plans are added to the solution in separate route plans to keep track the total cost of the DCARP scenario. Furthermore, if there are any new route plans within the solution, the framework gives them identifiers and also attempts to assign each of them to a free vehicle. For the other route plans, it can be easily determined which route plan identifier belongs to which route plan, based on the virtual task within them.

#### *3.3. The Basic ABC Algorithm*

This section introduces the basic ABC algorithm for combinatorial problems, based on [16]. Just like in the original ABC algorithm [13], the artificial bees are classified into the three groups:


In the ABC algorithm, a food source is corresponded to a solution and the nectar amount of a food source is corresponded to the fitness of a solution.

The ABC algorithm is an iterative process with four phases in total. It begins with the initial phase, then it iterates three bee phases (always in the same order) until a predefined termination criterion is met. In the initial phase, the population is initialized with randomly generated food sources. In the first phase, the employed bee phase, the employed bees are sent to the food sources, where they determine the nectar amounts of the food sources. In the second phase, the onlooker bee phase, the probability value of the sources are calculated based on their nectar amount, then the onlooker bees are sent to the preferred food source to find neighboring food sources and determine their nectar amount. In the third phase, the scout bee phase, the exploitation process of the sources exhausted by the bees are stopped and the scout bees are sent out to randomly discover new food sources within the search area. In each phase, the best food source found so far is memorized. The phases are described in more detail in the subsections below.

#### 3.3.1. Initialization Phase

In the initialization phase, the parameters and the population are initialized. The parameters of the ABC algorithm can be defined as follows:


The population is initialized by randomly generating *sn* number of food sources and assigning one employed bees to each of them. The employed bees evaluate the fitness of these solutions.

#### 3.3.2. Employed Bee Phase

At this phase, each employed bee *xi* generates a new food source *xnew* in the neighborhood of its current position. Once *xnew* is obtained, it will be evaluated and compared to *xi*. If the nectar amount of *xnew* is equal to or higher than that of *xi*, *xnew* replaces *xi* and becomes a new member of the population, otherwise *xi* is retained. In other words, a greedy selection mechanism is employed between the old and the new candidate solutions.

#### 3.3.3. Onlooker Bee Phase

An onlooker bee evaluates the nectar information taken from all the employed bees and selects a food source *xi* depending on its probability value *pi* calculated by the following expression:

$$p\_i = \frac{fit\_i}{\sum\_{j=1}^{su} fit\_j} \tag{13}$$

where *fiti* is the nectar amount (i.e., the fitness value) of the *i*-th food source *xi*. The higher the value of *fiti* is, the higher the probability of that the *i*-th food source is selected.

Once the onlooker has selected her food source *xi*, she produces a modification on *xi* by using a local search operator. The local search operator randomly selects a position in the neighborhood of *xi*. As in the case of the employed bees, if the modified food source has a better or equal nectar amount than *xi*, the modified food source replaces *xi* and becomes a new member in the population.

#### 3.3.4. Scout Bee Phase

If a food source *xi* cannot be further improved through a predetermined number of trials limit, the food source is assumed to be abandoned, and the corresponding employed bee becomes a scout. The scout produces a food source randomly.

In the basic ABC algorithm, in each cycle, at most one scout bee goes outside to search for a new food source.

#### *3.4. Move Operators for the CARP*

In population-based evolutionary algorithms, to enrich the diversity of the population, move operators with different levels of step-size are utilized to generate new, neighboring solutions. These move operators can be divided into two main categories: small step-size

operators and large step-size operators. Small step-size operators can modify the position and/or the direction of the tasks within one or two route plans. In contrast, large step-size operators are able to modify more than two route plans. The most commonly used small step-size operators in the literature, which are used in this work as well, are inversion, (single) insertion, swap, and two-opt operators [6–8]. In this work, a novel small step-size operator is used as well, the sub-route plan operator, which is introduced in this work, in Section 5. The only large step-size operator used in this work is merge-split, which was introduced in [6]. It is called a large-step-size operator, since it is able to modify more than two route plans.

The inversion and the sub-route operators can only change the direction and the order of the tasks within one route plan, so they do not change the feasibility of the solution. In contrast, the insertion, the swap, and the two-opt operators may change the amount of demand that needs to be served in some of the route plans, so the feasibility of the solution may change, too. For this reason, based on the settings, the output solution of these operators could be different. If infeasible solutions are not accepted and the calculated output solution is infeasible, then the operator returns the original, input solution instead (assuming that the input solution is a feasible solution).

#### 3.4.1. Inversion Operator

The inversion operator randomly selects a task *t* ∈ *T* within the input solution. If this task has an inverse (i.e., ∃*inv*(*t*) ∈ *T*), then the operator replaces *t* with *inv*(*t*) within the solution, else it returns the input solution.

#### 3.4.2. Insertion Operator

The insertion operator randomly selects a task *t*<sup>1</sup> ∈ *T*, then replaces (inserts) it before or after another randomly selected task *t*<sup>2</sup> ∈ *T* within the input solution. The selected tasks can be in different route plans or in the same route plan, but they cannot be the same tasks (i.e., *t*<sup>1</sup> = *t*2).

It creates two potential output solutions, based on where *t*<sup>1</sup> is inserted (before or after *t*2). If *t*<sup>1</sup> has an inverse, then the operator creates other potential output solutions, which contains the inverse task of the task (i.e., *inv*(*t*1)) instead of *t*1. It selects the solution as the output solution which has the smallest total cost among the potential output solutions.

#### 3.4.3. Swap Operator

The swap operator randomly selects two tasks (*t*<sup>1</sup> and *t*2, where *t*1, *t*<sup>2</sup> ∈ *T*), then replaces them with each other (i.e., swaps them). Similarly to the insertion operator, the selected tasks can be from the same route plan or different route plans, but they cannot be the same tasks.

It creates potential output solutions, which contain one or two inverse task(s) of the selected tasks instead of the task(s). All the four possible combinations are considered. It selects the solution as the output solution which has the smallest total cost among the potential output solutions.

#### 3.4.4. Two-Opt Operator

The two-opt operator randomly selects two route plans (e.g., *r*<sup>1</sup> and *r*2) of the solution. Based on the selected two route plans, two cases exist for this move operator. If the selected two route plans are the same (i.e., *r*<sup>1</sup> = *r*2), then a sub-route plan (i.e., a part of the route plan) is selected randomly and its direction is reversed. If the selected two route plans are different (i.e., *r*<sup>1</sup> = *r*2), then these two route plans are randomly cut into four sub-route plans, and then two new potential output solutions are generated by reconnecting the four sub-route plans and the best one from them is selected. For example, *r*<sup>1</sup> and *r*<sup>2</sup> are cut into sub-route plans *r*<sup>11</sup> *r*<sup>12</sup> and *r*<sup>21</sup> *r*22, respectively. Two new solutions are generated by connecting them in the following ways: (1) *r*<sup>11</sup> with *r*<sup>22</sup> and *r*<sup>21</sup> with *r*<sup>12</sup> and (2) *r*<sup>11</sup> with reversed *r*<sup>21</sup> and reversed *r*<sup>12</sup> with *r*22.

#### 3.4.5. Merge-split Operator

As it was mentioned before, the merge-split operator can make large changes in the solution (e.g., it can modify the order of all the tasks within one or more route plans), so it is considered a large step-size operator. This operator randomly selects *x* number of different route plans in the input solution, where *x* is a random number (1 ≤ *x* ≤ |*S*|). It obtains an unordered list of tasks by merging the tasks of the selected route plans into one list, and then sorts this unordered list with a path scanning heuristic (e.g., [22], which is used in this work as well). The obtained ordered list is then optimally split into new route plans using Ulusoy's splitting procedure [24].

The ordered list is constructed by the path scanning heuristic the following way. First, an empty path is initialized, then, the affected tasks are added one by one into the current path, until no tasks are left in the unordered list. In each iteration, only those tasks are taken into account that can be added to the current path without breaking the capacity constraint. If there are no any tasks like that, then the depot is added to the current path and a new path is initialized (that becomes the current path). When a task or the depot is added to the current path, the task/depot is connected to the end of the current path with the shortest path between them. If there are multiple tasks that can be added, the one that is closest to the end of the current path is added. If there are multiple tasks that are closest to the end of the current path, then one of the following rules are applied to determine which task should be added next:


In the rules above, *t* (*t* ∈ *T*) is a task and *v*<sup>0</sup> (*v*<sup>0</sup> ∈ *V*) is the depot. In one run, only one of the rules can be used. Therefore, the path scanning heuristic is ran five times, which results in five ordered lists.

The Ulusoy's splitting procedure creates five new candidate output solutions from the five ordered lists by splitting the lists into route plans. How the procedure works is best summarized in [50]. The procedure starts with constructing the Directed Acyclic Graph (DAG) from the ordered list. A DAG is a graph with arcs that represent feasible sub-tours of one giant tour. Next, the shortest path through the graph is calculated, which gives the optimal partition of the giant tour into feasible route plans. As the final step, a new candidate solution is created from the untouched route plans of the input solution and the route plans returned by the procedure. From the five candidate solutions the best one is chosen and returned by the operator.

#### **4. The Proposed ABC Algorithm for CARP (CARP-ABC Algorithm)**

In this chapter, the ABC algorithm developed for CARP (CARP-ABC algorithm) is presented. The notations used for the CARP-ABC algorithm are collected in Table A2 in Appendix A.

The algorithmic description of the main CARP-ABC algorithm can be seen in Algorithm 1. The main algorithm can be divided into four main phases: initialization, employed bee, onlooker bee, and scout bee phases. The algorithm begins with the initialization phase, then enters a cycle, where it repeats the mentioned phases in the respective order until the termination criterion is satisfied (line 4). In the initialization phase (line 2), the colony (*C*), the age of the solutions within the colony (A), the global best solution (*S*∗), and its age (*α*∗) are initialized. In the employed bee phase (line 6), local search is performed around the members of the colony. In the onlooker bee phase (line 7), a more in-depth local search is performed around one solution from the colony. In the scout bee phase (line 8), global search is performed.

#### **Algorithm 1:** Main CARP-ABC algorithm.

```
input: I; ncs, nmi, ngsl, nlsl, nsal ∈ N
1 begin
2 C, A, S∗, α∗ ← initializationP(I, ncs);
3 t ← 0;
4 while t < nmi ∧ α∗ < ngsl do
5 C¯ ← C;
6 C, P ← employedBP(I, C, nlsl);
7 C ← onlookerBP(I, C, P, ncs, nlsl);
8 C, A, S∗, α∗ ← scoutBP(I, C¯, C, A, S∗, α∗, nsal);
9 t ← t + 1;
10 return S∗;
```
The parameters of the algorithm are the followings:


#### *4.1. Initialization Phase*

The algorithmic description of the initialization phase of the CARP-ABC algorithm can be seen in Algorithm 2. The algorithm in this phase first initializes the sets for the solutions *Si* for *i* = 1, 2, ... , *ncs* and their age (lines 2–3). To guarantee an initial population with certain quality and diversity, the solutions are generated randomly by using the Random Solution Generation (RSG) algorithm for CARP (line 6), which is introduced in this work (in Section 4.1.1). Other population based evolutionary algorithms usually use the Randomized Path-Scanning Heuristic (RPSH) [22] to generate initial solutions. However, our experiments showed that, in the case of the CARP-ABC algorithm, it does not improve the convergence speed of the algorithm, so only the RSG is used.

After the initialization of the colony, the algorithm selects the solution *Si* ∈ *C* with the best (highest) fitness value by using the *selectBestSolution* function (line 11). The fitness of a solution *Si* is defined by its total cost *TC*(*Si*) (it needs to be as small as possible). Therefore, the fitness value of a solution *Si* is computed by the following fit function:

$$fit(S\_i) = \frac{LB}{TC(S\_i)}\tag{14}$$

where *LB* is the lower bound of the solution (i.e., the total service cost of all the tasks, involving only one of each tasks which has an inverse). Its value ranges between 0 and 1. Solutions with greater fitness values are preferred, since greater fitness value means that the total cost of the solution is closer to the lower bound.


#### **Algorithm 2:** initializationP (Initialization Phase of the CARP-ABC algorithm).

#### 4.1.1. Random Solution Generation Algorithm

As the first step, the algorithm generates *ncs* random permutations, which contain positive integers from 1 to *n* (i.e., the id of every arc task and one of the two *id*s of every edge task) in random order. As the next step, the algorithm reads the *id*s in the permutation one-by-one from left to right, while summing up the demand of the corresponding tasks. If the task assigned to the currently read *id* would break the capacity constraint of the current route plan, the algorithm inserts a "0" (the *id* of the dummy task *t*0) before the *id* of the task in the sequence (i.e., the task is added to a new route plan). After it is finished with separating the *id*s of the tasks into route plans, the algorithm checks each task in the solution and, if it has an inverse task, then randomly (e.g., with 0.5 probability) replaces the id of the task with the *id* of its inverse. As the final step, the algorithm inserts a "0" as the first and last task of the solution to make it a valid solution.

#### *4.2. Employed Bee Phase*

The algorithmic description of the employed bee phase of the CARP-ABC algorithm can be seen in Algorithm 3. The algorithm in this phase, for each employed bee, generates new candidate solutions in the neighborhood of *Si* with each small step-size operator, then evaluates and selects the best solution (lines 2–11). In this phase, only the inversion operator (line 6) and the sub-route plan operator (line 7) are used, because only these operators guarantee that the new candidate solution will be feasible. It is repeated until the known best local solution *S*∗ *<sup>i</sup>* cannot be improved within the defined number of iterations (i.e., its age reached *nlsl*). If the fitness value of the new candidate solution *S*<sup>∗</sup> *<sup>i</sup>* is greater than or equal to the fitness value of the current solution *Si*, the new solution replaces the current one in the population (lines 14–15).

As the next step, the algorithm calculates the winning probability values *pi* for the solutions *Si* (lines 16–19). The probability values *pi* are calculated with the same function as in the basic ABC algorithm (Equation (13)).

**Algorithm 3:** employedBP (Employed Bee Phase of the CARP-ABC algorithm).


#### *4.3. Onlooker Bee Phase*

The algorithmic description of the onlooker bee phase of the CARP-ABC algorithm can be seen in Algorithm 4. The algorithm in this phase, depending on the *pi* values, selects a solution *Si* with the *selectSolution* function. This function first performs a roulette selection *Int*( <sup>√</sup>*ncs*) times to select *Int*( <sup>√</sup>*ncs*) number of solutions from the colony. Next, it compares the selected solutions to each other and selects the best one from them (i.e., the one with the greatest fitness value).

As a next step, the algorithm generates *ncs* number of new candidate solutions *Si*,*<sup>j</sup>* in the neighborhood of *Si* (i.e., one solution for each onlooker bee) with the merge-split operator (lines 5–6). It generates new candidate solutions in the neighborhood of these solutions with the small step-size operators, until the known best local solution *S*∗ *<sup>i</sup>*,*<sup>j</sup>* cannot be improved within the defined number of iterations (i.e., the age of the solution, *α*∗ *i*,*j* , reaches *nlsl*) (lines 7–21). In this phase, all the small step-size operators (i.e., inversion, insertion, swap, 2-opt, and sub-route plan) are applied to *S*- *i*,*j* , which is the best solution that was found in the previous iteration (lines 11–15). From the resulted solutions, the best one is chosen with the *selectBestSolution* function as the new *S*- *<sup>i</sup>*,*<sup>j</sup>* (line 16). If the new *S*- *<sup>i</sup>*,*<sup>j</sup>* is better than the currently known best solution in the neighborhood of *Si*,*<sup>j</sup>* (i.e., *S*<sup>∗</sup> *i*,*j* ), then it is set as the new best solution *S*∗ *<sup>i</sup>*,*<sup>j</sup>* (lines 17–19). Otherwise, the age of *S*<sup>∗</sup> *<sup>i</sup>*,*<sup>j</sup>* (i.e., *α*<sup>∗</sup> *i*,*j* ) is increased by one (lines 20–21). After the search ends in the neighborhood of *Si*,*j*, it is checked whether the best solution found (i.e., *S*∗ *i*,*j* ) is better than the best solution found in the whole neighborhood of *Si* (i.e., *S*∗ *<sup>i</sup>* ). If *S*<sup>∗</sup> *<sup>i</sup>*,*<sup>j</sup>* has a higher fitness value and is also feasible (its total excess demand is zero), then it will be set as the new *S*∗ *<sup>i</sup>* (lines 22–23).

If the best solution found in this phase (i.e., *S*∗ *<sup>i</sup>* ) is better than the current solution *Si*, then *Si* is replaced by *S*<sup>∗</sup> *<sup>i</sup>* in the colony (lines 25–26).

**Algorithm 4:** onlookerBP (Onlooker Bee Phase of the CARP-ABC algorithm).

**input:** *<sup>I</sup>*; *<sup>C</sup>* <sup>=</sup> {*S*1, *<sup>S</sup>*2,..., *Si*,..., *Sncs*}; *<sup>P</sup>* <sup>=</sup> {*p*1, *<sup>p</sup>*2,..., *pi*,..., *pncs*}; *ncs*, *nlsl* <sup>∈</sup> <sup>N</sup> **<sup>1</sup> begin <sup>2</sup>** *Si* ← *selectSolution*(*C*, *P*); **<sup>3</sup>** *S*∗ *<sup>i</sup>* ← *Si*; **<sup>4</sup>** *j* ← 1; **<sup>5</sup> while** *j* ≤ *ncs* **do <sup>6</sup>** *Si*,*<sup>j</sup>* ← *mergeSplit*(*Si*); **<sup>7</sup>** *S*∗ *<sup>i</sup>*,*<sup>j</sup>* ← *Si*,*j*; **<sup>8</sup>** *S*- *<sup>i</sup>*,*<sup>j</sup>* ← *Si*,*j*; **<sup>9</sup>** *α*∗ *<sup>i</sup>*,*<sup>j</sup>* ← 0; **<sup>10</sup> while** *α*∗ *<sup>i</sup>*,*<sup>j</sup>* < *nlsl* **do <sup>11</sup>** *Si*,*j*,1 ← *inverse*(*I*, *S*- *i*,*j* ); **<sup>12</sup>** *Si*,*j*,2 ← *insert*(*I*, *S*- *i*,*j* ); **<sup>13</sup>** *Si*,*j*,3 ← *swap*(*I*, *S*- *i*,*j* ); **<sup>14</sup>** *Si*,*j*,4 ← *twoOpt*(*I*, *S*- *i*,*j* ); **<sup>15</sup>** *Si*,*j*,5 ← *subRoutePlan*(*I*, *S*- *i*,*j* ); **<sup>16</sup>** *S*- *<sup>i</sup>*,*<sup>j</sup>* ← *selectBestSolution*({*Si*,*j*,1, *Si*,*j*,2, *Si*,*j*,3, *Si*,*j*,4, *Si*,*j*,5}); **<sup>17</sup> if** *fit*(*S*∗ *i*,*j* ) < *fit*(*S*- *i*.*j* ) **then <sup>18</sup>** *S*∗ *<sup>i</sup>*,*<sup>j</sup>* ← *S*- *i*,*j* ; **<sup>19</sup>** *α*∗ *<sup>i</sup>*,*<sup>j</sup>* ← 0; **<sup>20</sup> else <sup>21</sup>** *α*∗ *<sup>i</sup>*,*<sup>j</sup>* ← *α*<sup>∗</sup> *<sup>i</sup>*,*<sup>j</sup>* + 1; **<sup>22</sup> if** *fit*(*S*∗ *<sup>i</sup>* ) < *fit*(*S*<sup>∗</sup> *i*,*j* ) ∧ *totalExcessDem*(*S*<sup>∗</sup> *i*,*j* ) = 0 **then <sup>23</sup>** *S*∗ *<sup>i</sup>* ← *S*<sup>∗</sup> *i*,*j* ; **<sup>24</sup>** *j* ← *j* + 1; **<sup>25</sup> if** *fit*(*Si*) < *fit*(*S*<sup>∗</sup> *<sup>i</sup>* ) **then <sup>26</sup>** *Si* ← *S*<sup>∗</sup> *i* ; **<sup>27</sup> return** *C;*

#### *4.4. Scout Bee Phase*

The algorithmic description of the scout bee phase of the CARP-ABC algorithm can be seen in Algorithm 5. The algorithm in this phase increases the age of unchanged solutions (lines 3–4) and sets the age to zero for new solutions (lines 8–9) within the colony. Furthermore, if there is an abandoned solution (i.e., a solution which could not be improved through a predetermined number of trials, which is called *nsal*), the algorithm replaces it with a new solution (lines 5–7), which is generated by using the RSG algorithm (as in the initialization phase).

In this phase, the algorithm also updates the global best solution, *S*∗. First, the best solution of the new colony is selected with the *selectBestSolution* function as solution *S*- (line 10). If *S* is better than *S*∗, then *S* is set as the new global best solution (lines 11–13). Otherwise, the age of *S*∗ (i.e., *α*∗) is increased by one (lines 14–15).

**Algorithm 5:** scoutBP (Scout Bee Phase of the CARP-ABC algorithm).

**input:** *<sup>I</sup>*; *<sup>C</sup>*¯ <sup>=</sup> {*S*¯ 1, *S*¯ 2,..., *S*¯ *<sup>i</sup>*,..., *S*¯ *ncs*}; *C* = {*S*1, *S*2,..., *Si*,..., *Sncs*}; <sup>A</sup> <sup>=</sup> {*α*1, *<sup>α</sup>*2,..., *<sup>α</sup>i*,..., *<sup>α</sup>ncs*}; *<sup>S</sup>*∗; *<sup>α</sup>*∗, *nsal* <sup>∈</sup> <sup>N</sup> **<sup>1</sup> begin <sup>2</sup> forall** *Si* ∈ *C* **do <sup>3</sup> if** *fit*(*Si*) = *fit*(*S*¯ *<sup>i</sup>*) **then <sup>4</sup>** *α<sup>i</sup>* ← *α<sup>i</sup>* + 1; **<sup>5</sup> if** *α<sup>i</sup>* = *nsal* **then <sup>6</sup>** *Si* ← *RSG*(*I*); **<sup>7</sup>** *α<sup>i</sup>* ← 0; **<sup>8</sup> else <sup>9</sup>** *α<sup>i</sup>* ← 0; **<sup>10</sup>** *S*- ← *selectBestSolution*(*C*); **<sup>11</sup> if** *fit*(*S*∗) < *fit*(*S*- ) **then <sup>12</sup>** *S*<sup>∗</sup> ← *S*- ; **<sup>13</sup>** *α*<sup>∗</sup> ← 0; **<sup>14</sup> else <sup>15</sup>** *α*<sup>∗</sup> ← *α*<sup>∗</sup> + 1; **<sup>16</sup> return** *C*, A, *S*∗, *α*∗*;*

#### *4.5. Computational Complexity Analysis*

In this section, the computational complexity of the proposed CARP-ABC algorithm is discussed. The computational complexity is expressed by using the big-*O* notation. For the sake of simplicity, approximations are used and the constant values are omitted. The computational complexity of the whole algorithm depends on the given parameter values and the complexity of the input CARP instance, mainly on *n* (i.e., the number of tasks that have to be executed).

#### 4.5.1. Initialization Phase

The computational complexity of the initialization phase is *O*(*ncs* ∗ *n* + *ncs*), in which *O*(*n*) is the complexity of the RSG algorithm and *O*(*ncs*) is the complexity of selecting the best solution. *O*(*n*) is multiplied by *ncs*, because RSG is executed *ncs* times to create the initial population.

Within the RSG algorithm, the complexity of generating a random permutation of the task identifiers is *O*(*n*), assuming that the Fisher–Yates shuffle algorithm [51] is used for it. After a permutation is generated, the algorithm iterates over each element, which also has *O*(*n*) as complexity. Therefore, the complexity of the RSG algorithm is around *O*(2 ∗ *n*), which is *O*(*n*) if the constant multiplier is omitted.

#### 4.5.2. Employed Bee Phase

The computational complexity of the employed bee phase is *O*(*ncs* ∗ *nlsl* ∗ (*n* + log *n* + *nmax*) + *ncs*), in which the complexity of the local search is *O*(*nlsl* ∗ (log *n* + *n* + *nmax*)). The complexity of the probability calculation is *O*(*ncs*) (assuming, that the sum of the fitness values is calculated only once). *O*(*nlsl* ∗ (log *n* + *n* + *nmax*)) is multiplied by *ncs*, because the local search is executed for all the *ncs* members of the population.

Within the local search, the complexity of the inversion operator is *O*(*n*) and the complexity of the sub-route plan operator is *O*(log *n* + *n* + *nmax*). Within the sub-route plan operator, the complexity of selecting a route plan is *O*(log *n*), since in the worst case |*S*| = *n* (i.e., every task is on a separate route). After a route plan is selected, one of the methods of the operator is executed. From the methods, the sub-route plan rotation method has the greatest complexity, which is *O*(*n* + *nmax*), since in the worst case *lk* = *n* (i.e., there is only one route plan in the solution).

#### 4.5.3. Onlooker Bee Phase

The computational complexity of the onlooker bee phase is *O*(*ncs* + *ncs* ∗ (log *n* + *n* + *nlsl* ∗ (*n* + log *n* + *nmax*))), in which the complexity of selecting a solution from the colony is *O*(*k* ∗ *ncs*) or *O*(*ncs*) is the constant *k* is omitted. The complexity of the main search is *O*(*ncs* ∗ (log *n* + *n* + *nlsl* ∗ (*n* + log *n* + *nmax*)).

Within the main search, the complexity of the merge-split operator is *O*(log *n* + *n*), because the complexity of selecting the number of route plans is *O*(log *n*) and the complexity of the other components of the operator (i.e., selecting the route plans, collecting the affected tasks, and executing the RPSH) is *O*(*n*). The complexity of RPSH is *O*(*n*), since in the worst case all the *n* tasks are affected in the solution. After the merge-split operator returns a solution, search is performed around this solution. The complexity of this search is *O*(*nlsl* ∗ (*n* + log *n* + *nmax*)), because the complexity of the sub-route plan operator is *O*(log *n* + *n* + *nmax*) and the complexity of the other operators (i.e., inversion, insertion, swap, and two-opt) are *O*(*n*).

#### 4.5.4. Scout Bee Phase

The computational complexity of the scout bee phase is *O*(*ncs* ∗ *n* + *ncs*), because in the worst case, all the solutions have to be replaced in the colony for exceeding *nsal*), so RSG is executed *ncs* times. Next, the best solution is chosen from the colony, which has *Ocs* complexity.

#### 4.5.5. Whole Algorithm

The computational complexity of the whole CARP-ABC algorithm is composed of the complexity of the initialization phase and the multiplication of the other phases by *nmi*, since in the worst case the algorithm runs till the maximum number of iterations is reached. If the duplications are removed, then it is the following: *O*(*ncs* ∗ *n* + *ncs* + *nmi* ∗ (*ncs* ∗ *nlsl* ∗ (*n* + log *n* + *nmax*) + *ncs* + *ncs* ∗ (log *n* + *n* + *nlsl* ∗ (*n* + log *n* + *nmax*)))).

If the parameters of the CARP-ABC algorithm and the sub-route plan operator are set to a fixed value, then the computational complexity is the following: *O*(*n* + log *n*). Therefore, the time complexity of the CARP-ABC algorithm is mostly linear but it contains components with logarithmic time complexity (e.g., when a route plan is selected).

#### **5. Sub-Route Plan Operator**

The sub-route plan operator is based on the GSTM operator for TSP [47]. The main differences between the modified version and the original version are due to the differences between the TSP and the CARP. Therefore, the modified version works with arcs instead of nodes. Furthermore, since the solution for a TSP is always one route plan, while the solution for a CARP (usually) consists of more than one route plan, the modified version takes into account only a part of the solution (one route plan) instead of the whole solution.

The sub-route plan operator is a complex move operator which consists of two different greedy search methods (greedy reconnection and sub-route rotation) and a method that provides distortion. In all three methods, inversion of the affected tasks is considered. Inversion of the tasks has real importance when a sequence of tasks is inverted in the sub-route rotation method because when the execution order of the tasks changes, the direction in which the tasks are executed should be changed too, to keep the traveling cost minimal. The used notations within this chapter are collected in Table A3 in Appendix A.

#### *5.1. The Main Algorithm*

The algorithmic description of the main algorithm of the sub-route plan operator can be seen in Algorithm 6. As input, a CARP instance *I*, a solution *S* of *I*, and the parameters of the algorithm are expected. The parameters are the following:


**Algorithm 6:** subRoutePlan (main algorithm of the sub-route plan operator).


The maximal length of the sub-route plan (*lmax*) is determined after the route plan is selected. In the proposed CARP-ABC algorithm, the parameters of this algorithm are given as constant values, so only *I* and *S* are expected.

In the first step of the algorithm, a route plan *rk* is selected from the solution *S* (line 2), then, if the number of (not dummy) tasks within *rk* is sufficient (i.e., *lk* is greater than or equal to *lmin*, line 3), the algorithm proceeds to the next step. Otherwise, it returns the input solution *S* unchanged (line 21).

In the following step of the algorithm, the parameters are initialized and the (sub- )route plans are generated. The maximum length of the sub-route plan (*lmax*) is determined based on the number of tasks within *rk* (*lk*) and the predefined minimum length of the sub-route plan (*lmin*) (line 4), and the new route plan (*r*- *<sup>k</sup>*) is initialized (line 5). The length of the sub-route plan *l* is determined randomly based on *lmin* and *lmax* (line 6). The position index of the starting task of the sub-route plan (*s*) is randomly selected taking into account *l* (line 7). The position index of the ending task of the sub-route plan (*e*) is determined by *s* and *l* (line 8). The sub-route plan *r*∗ *<sup>k</sup>* is constructed by taking the sub-route plan that is enclosed by the tasks *tk*,*<sup>s</sup>* and *tk*,*<sup>e</sup>* from *rk* (line 9). The route plan without *r*<sup>∗</sup> *<sup>k</sup>* is denoted by *r*# *<sup>k</sup>* (line 10).

As the next step of the algorithm, a random number is generated (*prnd*) between 0 and 1 (line 11), which determines the operation of the operator. If *prnd* is less than or equal to the predefined reconnection probability *prc* (i.e., *prnd* ≤ *prc*, line 12), then the greedy reconnection method is executed (line 13, Section 5.2), otherwise, a new random number is generated (line 15). If the new value of *prnd* is less than or equal to the predefined correction and perturbation probability *pcp* (i.e., *prnd* ≤ *pcp*, line 16), then distortion is added to *rk* (line 17, Section 5.3), otherwise, the sub-route plan rotation method is executed (line 19, Section 5.4). As the final step, the solution *S* is updated by removing the old route plan *rk* and adding the new one, *r*- *<sup>k</sup>* (line 20), then the updated solution is returned (line 21).

#### *5.2. Greedy Reconnection Method*

The greedy reconnection method inserts *r*∗ *<sup>k</sup>* into the position within *<sup>r</sup>*# *<sup>k</sup>* that generates the least amount of increase in the total cost of the route plan.

#### 5.2.1. Algorithm

The algorithmic description of the greedy reconnection method within the sub-route plan operator can be seen in Algorithm 7. As input, the CARP instance *I*, the original route plan *rk* of the solution *S*, the selected sub-route plan *r*<sup>∗</sup> *<sup>k</sup>* , and the truncated route plan *<sup>r</sup>*# *<sup>k</sup>* (i.e., *rk* without *r*<sup>∗</sup> *<sup>k</sup>* ) are expected.

**Algorithm 7:** greedyReconnection (algorithm of the greedy reconnection method within the sub-route plan operator).

**input:** *I*; *rk*; *r*<sup>∗</sup> *<sup>k</sup>* ; *<sup>r</sup>*# *k* **<sup>1</sup> begin <sup>2</sup>** *r*- *<sup>k</sup>* ← *rk*; **<sup>3</sup>** *i* ← 1; **<sup>4</sup> while** *i* ≤ *l* # *<sup>k</sup>* + 1 **do <sup>5</sup>** *rk*,*<sup>i</sup>* ← *insertSubRoutePlan*(*r*<sup>∗</sup> *<sup>k</sup>* ,*r*# *<sup>k</sup>* , *i*); **<sup>6</sup> if** *TC*({*rk*,*i*}) < *TC*({*r*- *<sup>k</sup>*}) **then <sup>7</sup>** *r*- *<sup>k</sup>* ← *rk*,*i*; **<sup>8</sup>** *i* ← *i* + 1; **<sup>9</sup> return** *r*- *k;*

In the first step of the algorithm, the new route plan *r*- *<sup>k</sup>* is initialized with the current route plan *rk* (line 2). The position index *i* that is used to find the best position for insertion of *r*∗ *<sup>k</sup>* into *<sup>r</sup>*# *<sup>k</sup>* is initialized as well (line 3). The value "1" refers to the first (not dummy) task within *r*# *<sup>k</sup>* (i.e., *t* # *<sup>k</sup>*,1). The value "0" would refer to the first dummy task (*t*0) and "*l* # *<sup>k</sup>* + 1" to the last dummy task within *r*# *<sup>k</sup>* (assuming *l* # *<sup>k</sup>* is the number of not dummy tasks within *<sup>r</sup>*# *k* ). In the following step, the algorithm checks each position within *r*# *<sup>k</sup>* to find the best one to insert *r*∗ *<sup>k</sup>* into (lines 4–8). In each iteration, before task *t* # *<sup>k</sup>*,*<sup>i</sup>* within *<sup>r</sup>*# *<sup>k</sup>* , it inserts *r*<sup>∗</sup> *<sup>k</sup>* with the *insertSubRoutePlan* function (line 5). The total cost of the resulting route plan (*rk*,*i*) is then compared with the total cost of *r*- *<sup>k</sup>* (line 6). If *rk*,*<sup>i</sup>* is better than *r*- *<sup>k</sup>* (i.e., it has lower total cost), then it becomes the new value of *r*- *<sup>k</sup>* (line 7).

In the final step of the algorithm, *r*- *<sup>k</sup>* is returned by the function (line 9).

#### 5.2.2. Example

For a better understanding of the method, see the following example. Let the selected route plan be *rk* = *t*0, *tk*,1, *tk*,2, ... , *tk*,13, *t*0 and the length of the sub-route plan be *l* = 3. Based on these, let the selected starting and ending task be *tk*,*<sup>s</sup>* = *tk*,5 and *tk*,*<sup>e</sup>* = *tk*,7, then the selected sub-route plan is *r*∗ *<sup>k</sup>* = *tk*,5, *tk*,6, *tk*,7 and the route plan *rk* without *r*<sup>∗</sup> *<sup>k</sup>* is *r*# *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,8, *tk*,9, *tk*,10, *tk*,11, *tk*,12, *tk*,13, *t*0. Let us assume that inserting *r*<sup>∗</sup> *k* between *tk*,8 and *tk*,9 in *r*# *<sup>k</sup>* results in the least amount of increase in the total cost of the solution, then *r*- *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,8, *tk*,5, *tk*,6, *tk*,7, *tk*,9, *tk*,10, *tk*,11, *tk*,12, *tk*,13, *t*0 will be the new *k*-th route plan in the solution.

The example discussed in the previous paragraph is depicted in Figures 2 and 3. In both figures, the arc tasks served within a route plan are depicted with solid lines, and the other arcs, which are only traversed, are depicted with dashed lines. The original route plan *rk* can be seen in Figure 2. In Figure 3, the sub-route plan *r*<sup>∗</sup> *<sup>k</sup>* and the truncated route plan *r*# *<sup>k</sup>* are shown, highlighted with red and blue colors, respectively. The selected starting and ending tasks (*tk*,*<sup>s</sup>* and *tk*,*e*) are depicted with thicker lines.

**Figure 2.** Example route plan *rk*.

**Figure 3.** Greedy reconnection method: (**a**) sub-route plan *r*∗ *<sup>k</sup>* subtracted from route plan *rk*; (**b**) sub-route plan *r*∗ *<sup>k</sup>* connected to the route plan *<sup>r</sup>*# *k* .

#### *5.3. Distortion Method*

The distortion method takes the tasks in *r*∗ *<sup>k</sup>* and inserts them one-by-one into *<sup>r</sup>*# *k* , starting from the position index *s*, and by rolling or mixing with the predefined linearity probability (*pl*). Rolling means selecting the current last task in *r*<sup>∗</sup> *<sup>k</sup>* and mixing means selecting a random task in *r*∗ *k* .

#### 5.3.1. Algorithm

The algorithmic description of the distortion method within the sub-route plan operator can be seen in Algorithm 8. As input, the CARP instance *I*, the selected sub-route plan *r*∗ *<sup>k</sup>* , the truncated route plan *<sup>r</sup>*# *<sup>k</sup>* (i.e., *rk* without *r*<sup>∗</sup> *<sup>k</sup>* ), the position index *s* of the starting task of *r*∗ *<sup>k</sup>* within the original route plan *rk*, and the linearity probability *pl* are expected.

**Algorithm 8:** Distortion (algorithm of the distortion method within the sub-route plan operator).

```
input: I; r∗
            k ; r#
               k ; s; pl
1 begin
2 i ← s;
3 r-

      k ← r#
           k ;
4 while |r∗
             k | = 0 do
5 prnd ← randFloat(0, 1);
6 t ← t
              ∗
              k,|r∗
                k |
                 ;
7 if prnd ≤ pl then
8 j ← randInt(1, |r∗
                            k |);;
9 t ← t
                 ∗
                 k,j
                   ;
10 r-

         k ← insertTask(t,r-

                           k, i);
11 r∗
         k ← removeTask(t,r∗
                            k );
12 i ← i + 1;
13 return r-

             k;
```
In the first step of the algorithm, the position index *i* is initialized with *s* (line 2) and the new route plan *r*- *<sup>k</sup>* is initialized with the truncated route plan *<sup>r</sup>*# *<sup>k</sup>* (line 3). While there are tasks in *r*∗ *<sup>k</sup>* , the algorithm executes the following steps (lines 4–12). First, a random number is generated (*prnd*) between 0 and 1 (line 5) and the last task is selected from *r*<sup>∗</sup> *<sup>k</sup>* into *t* (line 6). If *prnd* is less than or equal to *pl*, then the value of *t* is changed into a random task from *r*<sup>∗</sup> *k* (lines 7–9). The selected task *t* is then inserted into the position *i* in *r*- *<sup>k</sup>* with the *insertTask* function (line 10) and removed from *r*∗ *<sup>k</sup>* with the *removeTask* function (line 11). The task is always inserted right after the previously inserted task. When *r*∗ *<sup>k</sup>* runs out of tasks (i.e., all the tasks within it have been inserted into *r*- *<sup>k</sup>*), the algorithm returns *r*- *<sup>k</sup>* (line 13).

#### 5.3.2. Example

For a better understanding of the method, see the following example. Let the selected route plan be *rk* = *t*0, *tk*,1, *tk*,2, ... , *tk*,13, *t*0 and the length of the sub-route plan be *l* = 3. Based on these, let the selected starting and ending task be *tk*,*<sup>s</sup>* = *tk*,5 and *tk*,*<sup>e</sup>* = *tk*,7, then the selected sub-route plan is *r*∗ *<sup>k</sup>* = *tk*,5, *tk*,6, *tk*,7 and the route plan *rk* without *r*<sup>∗</sup> *<sup>k</sup>* is *r*# *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,8, *tk*,9, *tk*,10, *tk*,11, *tk*,12, *tk*,13, *t*0. The length of *r*<sup>∗</sup> *<sup>k</sup>* is three, so in this case the algorithm has three iterations. Let the linearity probability be *pl* = 0.2.

In the first iteration, let the random number be *prnd* = 0.1. Since it is smaller than *pl*, a random task is selected from *r*∗ *<sup>k</sup>* . Let the selected task be *tk*,6. In this case, the new route plan is *r*- *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,6, *tk*,8, ... (i.e., *tk*,6 is inserted between *tk*,4 and *tk*,8) and *r*∗ *<sup>k</sup>* = *tk*,5, *tk*,7 (i.e., *tk*,6 is removed).

In the second iteration, let the random number be *prnd* = 0.8. Since it is greater than *pl*, the currently last task is selected from *r*<sup>∗</sup> *<sup>k</sup>* , which is *tk*,7. In this case, the new route plan is *r*- *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,6, *tk*,7, *tk*,8, ... (i.e., *tk*,7 is inserted between *tk*,6 and *tk*,8) and *r*∗ *<sup>k</sup>* = *tk*,5.

In the third iteration, since only one task left in *r*∗ *<sup>k</sup>* , regardless of the value of *prnd*, *tk*,5 is selected. Therefore, *r*- *<sup>k</sup>* = *t*0, *tk*,1, *tk*,2, *tk*,3, *tk*,4, *tk*,6, *tk*,7, *tk*,5, *tk*,8, ... and *r*<sup>∗</sup> *<sup>k</sup>* = *tk*,5. Since *r*∗ *<sup>k</sup>* is now empty, the algorithm returns *r*- *k*.

#### *5.4. Sub-route Plan Rotation Method*

The sub-route plan rotation method selects one neighbor task randomly from the neighbors of *tk*,*<sup>s</sup>* and *tk*,*<sup>e</sup>* (*tk*,*s*<sup>∗</sup> and *tk*,*e*<sup>∗</sup> , respectively), then inverts the sequence of tasks enclosed by *tk*,*<sup>i</sup>* and *tk*,*i*<sup>∗</sup> (including *tk*,*i*<sup>∗</sup> in the sequence), where (*i*, *i* <sup>∗</sup>) ∈ {(*s*,*s*∗),(*e*,*e*∗)}. The inversion of the sequence is performed in such a manner, that *tk*,*<sup>i</sup>* and *tk*,*i*<sup>∗</sup> (or *inv*(*tk*,*i*<sup>∗</sup> )) become direct neighbors in the new route plan *r*- *k*.

5.4.1. Algorithm

The algorithmic description of the sub-route plan rotation method within the sub-route plan operator can be seen in Algorithm 9. As input, the CARP instance *I*, the original route plan *rk*, the position index *s* of the starting task and the position index *e* of the ending task of *r*∗ *<sup>k</sup>* within the original route plan *rk*, and the size of the neighborhood *nmax* are expected.

**Algorithm 9:** subRoutePlanRotation (algorithm of the sub-route plan rotation method within the sub-route plan operator).

**input:** *I*; *rk*; *s*; *e*; *nmax* **<sup>1</sup> begin <sup>2</sup>** *s*<sup>∗</sup> ← *selectNeighborTask*(*rk*,*s*, *nmax*); **<sup>3</sup>** *e*<sup>∗</sup> ← *selectNeighborTask*(*rk*,*e*, *nmax*); **<sup>4</sup>** *r*- *<sup>k</sup>* ← *rk*; **<sup>5</sup> forall** (*i*, *i* <sup>∗</sup>) ∈ {(*s*,*s*∗),(*e*,*e*∗)} **do <sup>6</sup>** *r*- *<sup>k</sup>*,*<sup>i</sup>* ← ; **<sup>7</sup> if** *i* < *i* ∗ **then <sup>8</sup>** *r*- *<sup>k</sup>*,*<sup>i</sup>* ← *t*0, *tk*,1,..., *tk*,*i*·*tk*,*i*<sup>∗</sup> , *tk*,*i*∗−1,..., *tk*,*i*+1·*tk*,*i*∗<sup>+</sup>1,..., *tk*,*lk* , *t*0; **<sup>9</sup>** *j* ← *i* ∗; **<sup>10</sup> while** *i* < *j* **do <sup>11</sup> if** ∃*inv*(*tk*,*j*) **then <sup>12</sup>** *r*- *<sup>k</sup>*,*<sup>i</sup>* ← *replaceTask*(*tk*,*j*, *inv*(*tk*,*j*),*r*- *k*,*i* ); **<sup>13</sup>** *j* ← *j* − 1; **<sup>14</sup> else <sup>15</sup>** *r*- *<sup>k</sup>*,*<sup>i</sup>* ← *t*0, *tk*,1,..., *tk*,*i*∗−1·*tk*,*i*−1, *tk*,*i*−2,..., *tk*,*i*<sup>∗</sup> ·*tk*,*i*,..., *tk*,*lk* , *t*0; **<sup>16</sup>** *j* ← *i* − 1; **<sup>17</sup> while** *i* <sup>∗</sup> ≤ *j* **do <sup>18</sup> if** ∃*inv*(*tk*,*j*) **then <sup>19</sup>** *r*- *<sup>k</sup>*,*<sup>i</sup>* ← *replaceTask*(*tk*,*j*, *inv*(*tk*,*j*),*r*- *k*,*i* ); **<sup>20</sup>** *j* ← *j* − 1; **<sup>21</sup> if** *TC*({*r*- *k*,*i* }) < *TC*({*r*- *<sup>k</sup>*}) **then <sup>22</sup>** *r*- *<sup>k</sup>* ← *r*- *k*,*i* ; **<sup>23</sup> return** *r*- *k;*

In the first step of the algorithm, one position index of the *nmax* closest neighbor tasks is selected randomly for both *tk*,*<sup>s</sup>* and *tk*,*<sup>e</sup>* (*s*<sup>∗</sup> and *e*∗, respectively) with the *selectNeighborTask* function (lines 2–3) and the new route plan *r*- *<sup>k</sup>* is initialized with the original route plan *rk* (line 4). In the next step, for all (*i*, *i* <sup>∗</sup>) ∈ {(*s*,*s*∗),(*e*,*e*∗)}, the following steps are executed (lines 5–22). First, the potential new route plan *r*- *<sup>k</sup>*,*<sup>i</sup>* is initialized with an empty sequence (line 6), then, based on the relationship between *i* and *i* ∗, a sub-route plan is selected and inverted. If *tk*,*<sup>i</sup>* is before *tk*,*i*<sup>∗</sup> (i.e., *i* < *i* <sup>∗</sup>) in *rk*, then *tk*,*<sup>i</sup>* is directly followed by *tk*,*i*<sup>∗</sup> (or *inv*(*tk*,*i*<sup>∗</sup> )) in *r*- *<sup>k</sup>* (lines 7–8). Otherwise, if *tk*,*i*<sup>∗</sup> is before *tk*,*<sup>i</sup>* in *rk*, then *tk*,*i*<sup>∗</sup> (or *inv*(*tk*,*i*<sup>∗</sup> )) is directly followed by *tk*,*<sup>i</sup>* in *r*- *<sup>k</sup>* (lines 14–15). In both cases, each task that has an inverse and is within the inverted sub-sequence, is replaced by its inverse task in the new route plan *r*- *k* with the *replaceTask* function (lines 9–13, lines 16–20). From the sub-route plan rotations, the one that has the lowest total cost is chosen (lines 21–22) and returned by the algorithm (line 23).

#### 5.4.2. Determining the Neighborhood

The neighbors are determined according to the predefined size of the neighborhood (*nmax*). The distance between arc task *tk*,*<sup>i</sup>* and another arc task is calculated based on their order within *rk* and whether the other arc task has an inverse task (i.e., it is from an edge task) or not. Let *tk*,*<sup>j</sup>* be an arc task in *rk* that is not *tk*,*<sup>i</sup>* (i.e., *tk*,*<sup>i</sup>* = *tk*,*j*). If *tk*,*<sup>j</sup>* is before *tk*,*<sup>i</sup>* in *rk*

(i.e., *j* < *i*) and has inverse task (i.e., *inv*(*tk*,*j*) ∈ *T*), then the distance between the two arc tasks is the shortest path between the head vertices of *tk*,*<sup>j</sup>* and *tk*,*i*, so it can be calculated with the expression *mdc head*(*tk*,*j*), *head*(*tk*,*i*) . The shortest path is calculated starting from *head*(*tk*,*j*), because during the sub-route plan rotation *tk*,*<sup>j</sup>* gets reversed (i.e., it gets replaced by its inverse in the route plan) and it is known that the head vertex of the task is the same as the tail vertex of the inverse task (i.e., *tail inv*(*tk*,*j*) = *head*(*tk*,*j*)). If the task does not have inverse, then the shortest path is calculated starting from *tail*(*tk*,*j*). If *tk*,*<sup>j</sup>* is after *tk*,*<sup>i</sup>* in *rk* (i.e., *<sup>j</sup>* > *<sup>i</sup>*) and has inverse task, then *mdc tail*(*tk*,*i*), *tail*(*tk*,*j*) is calculated. If the task does not have inverse, then the shortest path is calculated ending at *head*(*tk*,*j*). The reasoning behind what expression to use in each case is summarized in Table 1.

**Table 1.** Summary table about what expression to use to calculate the distance between an arbitrary arc task *tk*,*<sup>j</sup>* from the route plan *rk* and *tk*,*<sup>i</sup>* (*i* ∈ {*s*,*e*}).


#### 5.4.3. Example

For a better understanding of the method, see the following example. Let the selected route plan be *rk* = *t*0, *tk*,1, *tk*,2, ... , *tk*,13, *t*0, the length of the sub-route plan be *l* = 3, and the size of the neighborhood be *nmax* = 5. (At this method, *l* only defines the distance between the two selected arc tasks, it has no effect on the length of the rotated sub-route plans.) Based on these, let the selected two arcs task be *tk*,*<sup>s</sup>* = *tk*,6 and *tk*,*<sup>e</sup>* = *tk*,8. Let us assume that all the arc tasks in *rk* are from an edge task, so they all have inverse task.

The route plan, the selected arc tasks, and their neighborhood are illustrated on Figure 4. The arc tasks served within the route plan are depicted with solid lines, and the other arcs, which are only traversed, are depicted with dashed lines. The arc tasks *tk*,*<sup>s</sup>* and *tk*,*<sup>e</sup>* and their neighborhood are highlighted with red and blue color, respectively. Only the arc tasks that are completely covered by the ellipses are part of the neighborhood. It must be noted that the ellipses are only for representational purpose. Since identifying the neighborhood is quite complex due to the distance calculation, this yields to that the shape that covers only the neighbor arc tasks varies. Based on this, for *tk*,*s*, the set of the neighbor arc tasks is {*tk*,2, *tk*,3, *tk*,7, *tk*,8, *tk*,11}, and for *tk*,*e*, it is {*tk*,4, *tk*,6, *tk*,7, *tk*,9, *tk*,10}.

For *tk*,*s*, let us assume that the selected neighbor task is *tk*,2 (i.e., *tk*,*s*<sup>∗</sup> = *tk*,2), then the sub-route plan is *r*∗ *<sup>k</sup>* = *tk*,2, *tk*,3, *tk*,4, *tk*,5 (Figure 5a). Since *tk*,*s*<sup>∗</sup> precedes *tk*,*<sup>s</sup>* (i.e., *s*<sup>∗</sup> < *s*), the sub-route plan is reversed in a manner that in the new route plan *inv*(*tk*,*s*<sup>∗</sup> ) is directly followed by *tk*,*s*. The reversed sub-route plan is *inv*(*tk*,5), *inv*(*tk*,4), *inv*(*tk*,3), *inv*(*tk*,2), therefore the new route plan is *r*- *<sup>k</sup>* = *t*0, *tk*,1, *inv*(*tk*,5), *inv*(*tk*,4), *inv*(*tk*,3), *inv*(*tk*,2), *tk*,6, ... (Figure 5b).

For *tk*,*e*, let us assume that the selected neighbor task is *tk*,10 (i.e., *tk*,*e*<sup>∗</sup> = *tk*,10), then the sub-route plan is *r*∗ *<sup>k</sup>* = *tk*,9, *tk*,10 (Figure 6a). Since *tk*,*e*<sup>∗</sup> is after *tk*,*<sup>e</sup>* (i.e., *e* < *e*∗), the sub-route plan is reversed in a manner that in the new route plan *inv*(*tk*,*e*<sup>∗</sup> ) directly follows *tk*,*e*. The reversed sub-route plan is *inv*(*tk*,10), *inv*(*tk*,9), therefore the new route plan is *r*- *<sup>k</sup>* = *t*0, *tk*,1,..., *tk*,8, *inv*(*tk*,10), *inv*(*tk*,9), *tk*,11,..., *tk*,13, *t*0 (Figure 6b).

**Figure 4.** Nearest neighbors of the arc tasks *tk*,*<sup>s</sup>* (*tk*,6) and *tk*,*<sup>e</sup>* (*tk*,8) within the route plan *rk*, when *nmax* = 5.

**Figure 5.** Rotation of the sub-route plan enclosed by *tk*,*s*<sup>∗</sup> and *tk*,*s*: (**a**) An arc task (here *tk*,2) is randomly selected as *tk*,*s*<sup>∗</sup> from the neighbor list of arc task *tk*,*<sup>s</sup>* (*tk*,6), thus a sub-route plan is obtained; (**b**) The sub-route plan is inverted, so *inv*(*tk*,*s*<sup>∗</sup> ) will be directly followed by *tk*,*<sup>s</sup>* in the new route plan *r*- *k*.

**Figure 6.** Rotation of the sub-route plan enclosed by *tk*,*<sup>e</sup>* and *tk*,*e*<sup>∗</sup> : (**a**) An arc task (*tk*,10) is randomly selected as *tk*,*e*<sup>∗</sup> from the neighbor list of arc task *tk*,*<sup>e</sup>* (*tk*,8), thus, a sub-route plan is obtained; (**b**) The sub-route plan is inverted, so *tk*,*<sup>e</sup>* will be directly followed by *inv*(*tk*,*e*<sup>∗</sup> ) in the new route plan *r*- *k*.

#### **6. Experiments**

The proposed CARP-ABC algorithm (ABC algorithm from now on in this section) along with the sub-route plan operator was implemented in Python (3.6) in the Spyder (4.2.1) development environment. To compare the ABC algorithm with other static CARP solvers, the HMA and the ACOPR were also implemented. To test the ABC algorithm and the HMA as complete rerouting algorithms within the DCARP framework and compare them to the minimal rerouting algorithm RR1, the implementations from our previous work [20] were used. Python programming language was chosen for the implementation, because the DCARP framework will be supported in the future with the PM4Py process mining platform, which is written in Python. The experiments were performed on a laptop PC with Windows 10 operation system, equipped with an Intel(R) Core(TM) i5-3320M 2.60 GHz 2-core CPU and 8 GB of RAM.

It must be signified, that the HMA and the ACOPR were implemented based on their algorithmic description (in [7] and in [8], respectively), because the original implemented version of them is not available. Therefore, the implementations used in this work might have errors that decrease the effectiveness of these algorithms.

In this section, first, the setups of the experiments are specified (Section 6.1). Next, the results of the CARP experiments (Section 6.2), the results of the DCARP experiments (Section 6.4), and the results of the operator experiments (Section 6.3) are discussed, in the respective order.

#### *6.1. Experimental Setups*

For the CARP experiments, five CARP instances of different sizes were used. The ABC algorithm, the HMA, and the ACOPR were run 30 times with a time limit set to 10 min and applied to the CARP instances, independently, then the recorded outputs were compared and analyzed.

For the DCARP experiments, one CARP instance of medium size was used. Since the initial DCARP instance is fundamentally a static CARP instance and the HMA is the currently known the most accurate metaheuristic for CARP, the HMA was used to obtain the initial solution. The travel and service logs and the events were generated with the algorithms introduced in [20]. For each event type, 15 events were independently generated, then executed on the initial instance, creating new DCARP instances. The RR1 algorithm, the ABC algorithm, and the HMA were run with a time limit set to 1 min and applied to these instances, independently, then the recorded outputs were compared and analyzed.

For the operator experiments, three CARP instances of different sizes were used. The sub-route plan operator and the other small step-size operators for CARP (i.e., inversion, insertion, swap, and two-opt) were used as local search operators within the employed bee phase of the ABC algorithm. The employed bee phase was chosen instead of the onlooker bee phase so the efficiency of the operators can be measured on solutions of different qualities, not only on good quality solutions. The (modified) ABC algorithm was run 30 times with a time limit set to 10 min and applied to the CARP instances, independently, then the recorded outputs were compared and analyzed.

At the CARP and the DCARP experiments, during the execution of each algorithm, the new global best solution and the time it took for the algorithm to find the new global best solution (i.e., the elapsed time since the beginning of the execution of the algorithm) were recorded and analyzed. At the operator experiments, during the execution of the algorithm, the number of search trials in which the move operators found a better local best solution (i.e., *S*∗ *<sup>i</sup>* ) was recorded and analyzed.

The used instances and the parameter settings of the used algorithms are specified in the subsections (in Section 6.1.1 and in Section 6.1.2, respectively).

#### 6.1.1. Used Instances

To test the CARP solvers, benchmark test sets are often used in the literature. These test sets can be divided into two main categories: synthetic (e.g., containing randomly generated instances) [52–54] and real-life based (containing examples based on real road networks and tasks) [2,28,29].

Since testing an algorithm on all the instances would be time-consuming, only the following five instances were selected for the CARP experiments:


For the DCARP experiments, only the "egl-e1-A" instance was used. For the operator experiments, the "egl-e1-A", "egl-s1-A", and "egl-g1-A" instances were used.

The EGL and the EGL-Large sets originate from the data of a winter gritting application in Lancashire (UK). The EGL set contains 24 instances in which two different graphs are combined with various attribute values. The instances "egl-e1-A" and "egl-s1-A" were selected to represent one of each graph. The EGL-Large set contains 10 instances in which the graph is the same but the number of task edges is 347 in 5 instances and 375 in the other 5 instances. The instances "egl-g1-A" and "egl-g2-A" were selected to represent one from both kinds of instances.

The attributes of all the five selected instances are briefly summarized in Table 2. These CARP instances were selected to represent CARPs of very different sizes, thus requiring very different complexity levels to solve them. It can be seen, that the "kshs1" instance is a small synthetic CARP instance, which means a small search space for the algorithm, so for problem difficulty it can be put into the easy category. The EGL instances are real-life based CARP instances, so they are naturally more complex. Based on their size, the difficulty of the "egl-e1-A" instance is medium, the difficulty of the "egl-s1-A" instance is hard, and the difficulty of the instances "egl-g1-A" and "egl-g2-A" are very hard.


**Table 2.** Attributes of the used CARP instances.

<sup>1</sup> Based on the literature.

#### 6.1.2. Parameter Settings

The ABC algorithm was tested with multiple parameter settings. Based on the results, the following settings provided the best quality results without unnecessarily increasing the running time of the algorithm, thus these are used in the experiments:


According to the investigation in [47], the ideal parameter values for the GSTM operator are the followings: *prc* = 0.5, *pcp* = 0.8, *pl* = 0.2, *lmin* = 2, *lmax* = *Int*( <sup>√</sup>*n*), and *nmax* = 5, where n is the number of cities. In the experiments, the same parameter values are used for the sub-route plan move operator. The only difference is that *n* is the number of tasks within the selected route plan.

For the HMA and the ACOPR, the optimal parameter settings defined by the corresponding works [7,8] were used. For the RR1, no parameters are needed.

#### *6.2. Results of the CARP Experiments*

The charts on Figures 7–9 show the convergence speed of 30 independent runs of the algorithms for the selected instances. As it was mentioned before, for these, the new global best solution and the time it took for the algorithm to find the new global best solution were recorded. The y-axis shows the total cost of the solution, and the x-axis shows the elapsed time since the algorithm started running in seconds. The different colors indicate the outputs of the different algorithms. The colored lines indicate the average convergence speed and the colored areas cover all the values that were recorded (i.e., the areas are enclosed by the minimum and maximum values). The closer the line is to the intersection of the axes, the better the convergence speed of the algorithm.

In the case of the "kshs1" instance (Figure 7), it can be seen that the convergence speed of the ACOPR and the ABC algorithm is around twice as fast as the speed of the HMA. However, even though the speed of the ACOPR and the ABC algorithm is nearly the same, the ACOPR algorithm has failed to find the best solution in 30 runs, thus making the ABC algorithm the best solver for CARPs of small size like the "kshs1" instance.

**Figure 7.** The convergence speed of 30 independent runs of HMA, ACOPR, and ABC algorithms on the "kshs1" instance, plotted on one chart.

**Figure 8.** The convergence speed of 30 independent runs of HMA, ACOPR, and ABC algorithm on the "egl-e1-A" instance, plotted on one chart, with time limit (*t* ≤ 600 s).

In the case of the "egl-e1-A" instance (Figure 8), the differences between the convergence speed of the algorithms start to show. It can be seen that in all cases, the ABC algorithm provides better solutions and faster, than the ACOPR algorithm. The HMA algorithm has a very slow cycle time, thus it has a very slow convergence speed as well. If time is not taken into account, the HMA can generally provide better solutions than the other algorithms.

**Figure 9.** The convergence speed of 30 independent runs of HMA, ACOPR, and ABC algorithm on the "egl-s1-A" instance, plotted on one chart with time limit (*t* ≤ 600 s).

In Table 3, the total cost of the globally best solution within different time limits is examined, based on 30 independent runs of the ABC algorithm and the HMA. The calculated statistics are the following: minimum (Min.), maximum (Max.), average (Avg.), and standard deviation (Std.). It can be seen that within 1 min, the ABC algorithm always provided better solutions. Within 5 min, in some cases, the HMA algorithm found better results (it has smaller Min. value), but in average the ABC algorithm still performed better (it has smaller Max. and Avg. values). Nevertheless, within 10 or more minutes, the HMA algorithm provided better solutions. Regardless of the time limit, the ABC algorithm is slightly more stable than the HMA algorithm, in terms of how much the solution varies for different runs (it has smaller Std. values).


**Table 3.** Statistics of the total cost of the globally best solution of the HMA and the ABC algorithms on the "egl-e1-A" instance, within different time limits.

In the case of the "egl-s1-A" instance (Figure 9 and Table 4), the differences between the convergence speed of the HMA and the ABC algorithm is more complex. It can be seen that before 200 s, the ABC algorithm performs better, between 200 and 400 s they perform around the same, then after 400 s the HMA performs better.


**Table 4.** Statistics of the total cost of the globally best solution of the HMA and the ABC algorithms on the "egl-s1-A" instance, within different time limits.

The results were similar for the "egl-g1-A" and "egl-g2-A" instances (Figure 10 and Tables 5 and 6). In most of the runs, the set time limit was not enough for the HMA to improve its initial solution, so only its initial solution was recorded. That is why the graph for the HMA looks like a straight line in Figure 10. As a result of this, in the measured time period, the ABC algorithm performed better than the HMA after around 100 s.

**Table 5.** Statistics of the total cost of the globally best solution of the HMA and the ABC algorithms on the "egl-g1-A" instance, within different time limits.


**Table 6.** Statistics of the total cost of the globally best solution of the HMA and the ABC algorithms on the "egl-g2-A" instance, within different time limits.


**Figure 10.** The convergence speed of 30 independent runs of HMA, and ABC algorithm on the "egl-g1-A" instance, plotted on one chart with time limit (*t* ≤ 600 s)

Based on the results, it can be concluded that the ABC algorithm can provide a good enough solution within a short amount of time. Since it has a small cycle time, the best global solution can be updated more frequently. The ABC algorithm is better than the ACOPR algorithm in all aspects. The ABC algorithm has faster convergence speed and finds better quality solutions, than the ACOPR algorithm. Furthermore, it is competitive with the HMA, when the running time of the algorithms is set to a short time interval.

#### *6.3. Results of the Operator Experiments*

The results of the operator experiments are summarized in Table 7. In each row, the percentage of the number of trials is shown, in which the operator (specified by the column header) found a better solution, compared to the total number of trials in which a better solution was found for the instance (that is specified in the first column). For the sake of simplicity, let us call this measure efficiency. It can be seen that the sub-route plan operator has the highest efficiency in all the three cases, thus, among the examined operators, it has the highest chance to improve the current solution, regardless of the problem size.

A correlation can be observed between the size and complexity of the CARP problem and the efficiency of the operators. As it was mentioned before, the complexity of the "egl-e1-A" instance is medium, the "egl-s1-A" instance is difficult, and the the "egl-g1-A" instance is the most difficult. By increasing the size of the problem, the efficiency of the inversion and the sub-route plan operator increases compared to the other operators, and by decreasing the size of the problem, the efficiency of the insertion, the swap, and the two-opt operators increases compared to the other operators.



Based on the results, it can be concluded that the sub-route plan operator is more likely to find a better solution than the other operators, especially when a greater modification is needed on the current solution (since it is a randomly generated solution and/or it is a solution of a larger CARP instance).

#### *6.4. Results of the DCARP Experiments*

The results for the task appearance events, the demand increase events, and the vehicle breakdown events for the "egl-e1-A" instance can be seen in Table 8, Table 9, and Table 10, respectively. In all the three tables, the first few columns contain the parameter values that can be used to reconstruct the event data components of the problem:


The best total cost calculated by the RR1 rerouting algorithm, the ABC algorithm, and the HMA are contained by the last three columns. The best total cost of each run is highlighted with bold font.

The results of all the events are summarized in Table 11. It can be seen that in total, the ABC algorithm performed better than the other examined algorithms (RR1 and HMA). The HMA performed better only at the vehicle breakdown events, but the difference is negligible. The RR1 algorithm gave the best results in nearly the same amount of times as the ABC algorithm, in case of task appearance and vehicle breakdown events.

It is not shown in the tables, but the RR1 algorithm has the shortest run time (in the test cases, it was always less than one second). The run time of the other algorithms (ABC and HMA) is approximately the same whether a DCARP or a CARP instance is used as input, since it is the complexity of the problem that mainly defines the convergence speed.

Based on the results, similar conclusions can be made as in the CARP experiments. The ABC algorithm outperforms the HMA for a certain period of time, but then the HMA slowly takes over the lead. If time is the priority, then in the case of task appearance and vehicle breakdown events, the RR1 algorithm should be used. If time and the quality of the solution are equally important, then the ABC algorithm should be used for all events. If the quality of the solution is the priority, then the HMA should be used.

**Table 8.** The best total costs of 15 independent runs, calculated by the RR1, ABC, and HMA within one minute, after the occurrence of a random task appearance event in the "egl-e1-A" instance.



**Table 9.** The best total costs of 15 independent runs, calculated by the RR1, ABC, and HMA within one minute, after the occurrence of a random demand increased event in the "egl-e1-A" instance.

**Table 10.** The best total costs of 15 independent runs, calculated by the RR1, ABC, and HMA within one minute, after the occurrence of a random vehicle breakdown event in the "egl-e1-A" instance.


**Table 11.** The total number of the best outputs (and their percentage compared to the total number of outputs) of the algorithms RR1, ABC, and HMA summarized for the "egl-e1-A" instance, for each event type.


#### **7. Conclusions and Future Work**

In this study, an ABC algorithm for the CARP (CARP-ABC) was developed along with a new move operator, the sub-route plan operator, which is utilized by the proposed CARP-ABC algorithm. The CARP-ABC algorithm was tested both as a CARP and a DCARP solver, then, its performance was compared with other algorithms. The results showed that for both CARP and DCARP instances, the CARP-ABC algorithm excels in finding a relatively good quality solution in a short amount of time. It makes the algorithm highly competitive with the currently most accurate CARP solver, the HMA, when the running time of the algorithms is limited to around one minute.

In the future, the CARP-ABC algorithm will be improved upon, to increase the accuracy of the algorithm without increasing its runtime. The goal is to make the algorithm better than the HMA, even when the running time is unlimited.

**Author Contributions:** Conceptualization, Z.N.; Formal analysis, Z.N.; Investigation, Z.N.; Methodology, Z.N.; Software, Z.N.; Supervision, Á.W.-S.; Validation, Z.N.; Visualization, Z.N.; Writing original draft, Z.N.; Writing—review & editing, Z.N., Á.W.-S. and T.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Széchenyi 2020 under the EFOP-3.6.1-16-2016-00015.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The used CARP datasets are available at https://www.uv.es/belengue/ carp.html (accessed on 21 June 2022). The output data of the experiments discussed in the paper are available at https://drive.google.com/file/d/1LFgct7Z8\_W\_yx\_CppVmN1kAYry3VGM8Q/ (accessed on 21 June 2022).

**Acknowledgments:** The authors acknowledge support from the Slovenian–Hungarian bilateral project "Optimization and fault forecasting in port logistics processes using artificial intelligence, process mining and operations research", grant 2019-2.1.11-TÉT-2020-00113, and from the National Research, Development and Innovation Office–NKFIH under the grant SNN 129364.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A**

This appendix shows all the notations that are used in this work, categorized by the context where they appear and along with their meaning.


**Table A1.** Notations used in the CARP and the DCARP.




#### **References**


## *Article* **A Path Planning Model for Stock Inventory Using a Drone**

**László Radácsi, Miklós Gubán, László Szabó \* and József Udvaros**

Faculty of Finance and Accountancy, Budapest Business School, 1149 Budapest, Hungary **\*** Correspondence: szabo.laszlo4@uni-bge.hu

**Abstract:** In this study, a model and solution are shown for controlling the inventory of a logistics warehouse in which neither satellite positioning nor IoT solutions can be used. Following a review of the literature on path planning, a model is put forward using a drone that can be moved in all directions and is suitable for imaging and transmission. The proposed model involves three steps. In the first step, a traversal path definition provides an optimal solution, which is pre-processing. This is in line with the structure and capabilities of the warehouse. In the second step, the pre-processed path determines the real-time movement of the drone during processing, including camera movements and image capture. The third step is post-processing, i.e., the processing of images for QR code identification, the interpretation of the QR code, and the examination of matches and discrepancies for inventory control. A key benefit for the users of this model is that the result can be achieved without any external orientation tools, relying solely on its own movement and the organization of a pre-planned route. The proposed model can be effective not only for inventory control, but also for exploring the structure of a warehouse shelving system and determining empty cells.

**Keywords:** drone; inventory management; GA model; route planning; warehouse

**MSC:** 90B05

#### **1. Introduction**

In large, lightly structured warehouses in logistics centers, especially those where different products from several companies are stored, it is often difficult to pinpoint the exact location of stored goods. This is primarily the case when storage is not carried out with the help of automatic forklifts, or when errors occur in the registered and actual location of the goods due to mistakes during the picking process. To make matters worse, it is very difficult or in some cases impossible to use GPS-based identifications in these warehouses. If Time of Flight (ToF) cameras are not available, the positioning of the automatic devices is not possible. In this study, a model capable of updating the stock-on-hand registry using a drone is introduced, and the proposed model applies to a specific warehouse.

Storage (loading and unloading) and picking, unit load formation, and labeling take place in the warehouse building. In the warehouse, a double shelving system is used. Between the double shelves, there is a wide corridor in which the drone can travel. A shelf is divided into rows or compartments (cells). Within one compartment there can be one, two, or three locations created. Each unit stored in these locations is identified with a QR code that contains all the important data that are relevant to our solution.

Driving down the aisle, the drone searches for the QR codes placed on the outer cartons, and as soon as the drone camera finds one, a photo is taken and it is sent to a processing device (in our study, it is a large capacity tablet). According to the proposed model, the drone knows its exact location in the warehouse, i.e., the 3D data inside the warehouse are defined. This and the QR code or codes sent (in the case of a multi-cell location) can tell the inventory management model exactly where the product is located, or if the location is empty. Further processing depends solely on the capabilities of the application, which is beyond the scope of this study.

**Citation:** Radácsi, L.; Gubán, M.; Szabó, L.; Udvaros, J. A Path Planning Model for Stock Inventory Using a Drone. *Mathematics* **2022**, *10*, 2899. https://doi.org/10.3390/ math10162899

Academic Editors: Zsolt Tibor Kosztyán, Zoltán Kovács and Ripon Kumar Chakrabortty

Received: 3 July 2022 Accepted: 10 August 2022 Published: 12 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The mid-range drone and its camera can take a good quality picture from the middle of the aisle and identify the stored goods in the picture based on the blue and red colors used on the shelves. This helps the drone to move forward, backward, up, and down in the middle of the aisle, with no right-to-left movement. At the end of the aisle, the dock is manually placed where the battery can be charged and replaced.

After this introduction, an optimal route is described that is optimized primarily for the operating time of the drone, i.e., the drone will continue checking the locations until its battery wears out.

Based on the literature review, we have identified a gap in the research, as we have not encountered an approach such as ours to inventory accounting; therefore, we can say that the drone inventory based on this paper can be considered as a development in the existing discussion in this field.

Our study is founded on the following research question: How can warehouse stocking be automated by drones? Within this question, we pose a sub-question: How can the time of stocking be optimized by drones? In this study, we aimed to develop an algorithm to mathematically answer these questions; thus, we have proposed a model for automated warehouse stocking.

#### **2. Path Planning**

Our analysis suggests that there are very different approaches to path planning for the problem we are considering. In general, these studies focus primarily on a specific problem. Furthermore, the solutions that are used also vary (e.g., [1]).

To begin with, we investigated several solution methods. According to the literature, studies have used the Swarm Intelligence approach to solve, for example, the collision avoidance problem of robots in 3D media such as water [2] using a velocity-matching method. The route planning used in this paper is based on the ant algorithm. Our investigations have shown that this solution does not provide sufficient convergence and that the motions in our problem are different from it, but several elements of the basic model can be incorporated into our model. A similar and similarly unsuccessful attempt based on swarm intelligence is presented in [3], where an evolutionary algorithm was also chosen as the suitable approach. Another approach that can be used concerns automatic parking systems. It essentially provides a similar baseline for the path selection in a planar domain to the vertical-plane path-search algorithm. Of course, the authors focus primarily on the narrow setting, i.e., they focus on maneuvering, and their path-finding algorithm focuses on finding a parking space, although this is solved using a less well-known RRT (Rapidexploring Random Trees) algorithm. The chosen method for the model in [4] provides a useful solution in this direction. In this paper, the route planning of UAVs (Unmanned Aerial Vehicles) under terrain conditions is considered. The path is based on a polynomial model, a solution whereby only discontinuous linear motions need to be simulated. The advantage of their method is that it improves the initial population after construction using an ACO.

Several articles deal with logistics inventory control models and propose a specific solution. In their article, F. Benes et al. describe that, in the case of large outdoor warehouses, general identification methods are lengthy and inadequate. One way to easily and quickly determine the inventory is to deploy a UAV (unmanned aerial vehicle) for product identification. In this case, however, there is a problem in determining the location of the goods. A drone moves at a higher altitude, which can lead to a situation where we cannot accurately determine the location of the goods. This paper deals with the development of the definition of the correct flight level, which is suitable for distinguishing one of the identified elements at a distance of at least 2 m. The evaluation is based on an RSSI (received signal strength indicator) value. The experiment proved that the two objects can be distinguished even at the maximum reading distance of the selected passive UHF RFID tags [5]. This solution does not provide a suitable solution, since in this case the products

are placed in one flat, whereas we need to check the inventory of a warehouse with a shelf system (3D).

Indoor drone or unmanned aerial vehicle (UAV) operations in automated or pilotcontrolled drone use cases are addressed by Kurt Geebelen et al. Automated indoor flights have stricter requirements for stability and localization accuracy compared to classic outdoor use cases, which rely primarily on (RTK) GNSS for localization. In this paper, the effect of multiple sensors on 3D indoor-position accuracy is investigated using the OASE flexible sensor fusion platform. This evaluation is based on real drone flights in an industrial laboratory, with mm-accurate ground truth measurements provided by motion capture cameras, which enable the sensors to be evaluated based on their deviation from ground reality in 2D and 3D. The sensors considered in the research are: IMU, the sonar, SLAM camera, ArUco markers, and Ultra-Wideband (UWB). The article shows that with this setup, the achievable 2D (3D) indoor localization error varies from 4.4 cm to 21 cm depending on the sensor set selected. They also include cost/accuracy tradeoffs to indicate the relative importance of different sensor combinations depending on the (engineering) budget and use case. These laboratory results were validated in a Proof-of-Concept deployment of an inventory-scanning drone with more than 10 flight hours in a 65,000 m<sup>2</sup> warehouse. By combining the laboratory results and real-world deployment experience, the different subsets of the sensors represent a minimum viable solution for three different indoor use cases, considering accuracy and cost: a large drone with low weight and cost constraints, one or more medium-sized drones, and a swarm of weight- and cost-constrained nano drones [6]. Our solution also eliminates these flight inaccuracies, as the drone takes a picture from the approaching position, which is processed by software (which extracts the QR code from it). Therefore, the inaccuracies recommended in the article are negligible in our case.

As part of their research project, A. Rhiat, L. Chalal and A. Saadane developed a prototype named "Smart Shelf" that simulates a smart warehouse where mobile robots with grips managed by the ROS "Robotic Operating System" can autonomously navigate through inverse kinematics and different obstacles between other robots; on the other hand, RFID and iBeancon technology have a Smart Shelf to manage stocks. All items on the shelf are identified by RFID tags. The connection of the robots must be applicable to the various predefined or non-obstacles to optimize their search and items they encounter and to accessing the local network through a predefined map in the database. In addition, the embedding of Bin Packing Optimization techniques helps to improve the utilization of static volumes. Optimization algorithms can take into account some robotic constraints, including accessibility, improving the quality of placements and minimizing damaged goods. Their project aims to minimize human intervention and gain time [7]. We consider the solution very useful; it can also be used in our case. However, a problem is caused by the fact that we are inspecting the stock during the storage process, while material handling is taking place. Therefore, we have to discard this solution as well, since it is not possible to use the Smart Shelf during the material-handling process.

Haishi Liu et al. show that tobacco companies must regularly take inventory of finished products as well as raw and auxiliary materials, and drones with radio frequency identification (RFID) readers are becoming a major application trend for inventory applications. Under the condition of ensuring the accuracy of the inventory, this paper considers the limitations of the drone's physical performance, the limitations of the RFID reader, etc., and presents the power of the drone with respect to modeling and a task-planning model for the UAV. An inventory library equipped with an RFID reader is recommended. Thus, considering the problem whereby the greedy strategy in the traditional differential evolution (DE) algorithm causes the loss of location information preserved by other individuals, we propose a hybrid DE algorithm based on lion swarm optimization. Finally, the proposed algorithm was verified by environmental modeling based on data from the tobacco industry warehouse [8]. This solution is already similar to our model, but our proposed Generation Algorithm (OGA) simplifies the solution.

Our study also concerns the determination of optimal routes. The location search strategies used during picking can also provide a good basis, as the goal is to work as quickly and accurately as possible. Based on this, two basic strategies are distinguished: in the case of the strategy regarding the increasing horizontal location coordinates, the order of visiting each storage location corresponds to the increase in the horizontal coordinates of the storage sites. The zoning strategy is a variation of the strategy of increasing the horizontal location coordinates, in which the scaffolding field is first divided into (even-numbered) superimposed zones, and storage locations within each zone are accessed in ascending order of the horizontal location coordinates [9]. However, regarding warehouse robotization, in many cases, besides route planning, the other task is to balance the size of the robot fleet and the load stabilizing between the robots, which is not the case of this study [10].

Moving to the next aisle is a manual movement, as the aisles may not be processed in a sequential way, there may be a closed aisle, or the next aisle might be in another part of the building. This adds another aim to our proposed model of enabling the drone to recognize its own position.

These types of solutions are already used in several places, mainly for large multinational companies. There are some excellent solutions such as the Eyesee system developed by the Hardis Group, a comprehensive drone inventory solution that includes a drone capable of flying unmanned and equipped with a system for automatically capturing and identifying barcode data and handling the collected data via Amazon Web Services [11].

A similar problem is analyzed in [12], involving the target tracking of a drone moving in a defined closed area, but this goes beyond the scope of our study, as their goal is to recognize the target and choose its speed; however, it is still of relevance. The article chooses the Fuzzy solution, which is a viable option for this study as well, since in many cases the position of the QR code within the given location is not clear and the camera unfortunately does not see it in all cases. For this reason, fuzzy control after the basic and deterministic movement should be incorporated into the model. This can also be applied to differences due to columns. The authors of [13] solve the positioning of a drone in a GPS-free environment using a ToF camera used in robotics and self-guidance. However, this approach demands considerable investment into hardware as well. The camera is installed on the ceiling, and this monitors the positioning of the drone in the x, y direction and changes in height. Important technical solutions are presented to prevent interference from rotors. With the use of Gaussian function filtering, it was feasible to accurately analyze the situation in 3D. The article also provides an exact algorithm as a solution. Another methodology worth considering for closed-area drone control is presented in [14]. The authors use a voxel model and perform two types of route calculations: one for the shortest route and the other for the cheapest route.

Their problem is mainly the use of obstacle-avoidance image analysis, for which they use the distance transformation method for an abnormal image laid down by A. Rosenfeld and J. L. Pfaltz in *Distance Functions on Digital Pictures*, from the journal *Pattern Recognition* [15]. The obstacles, in the form of shelves, are assumed to be fixed, which is problematic because this is not always the case. Their approach keeps the drone safe from the obstacle within an effective range. Another article suggests a search for known trajectory options [16]. In their proposal, a usable genetic algorithm is used, which would be very useful in choosing the optimal crawl route for our study as well. The initial population is selected by a randomly generated greedy algorithm. An interesting solution is a low-complexity, machine-learning-based algorithm for optimizing the location of DBSs (drone base stations) [17] based on minimizing the collective wireless reception signal strength experienced by active terminals. The proposed algorithm reduces propagation loss in the system and provides a lower bit error rate compared to the Euclidean cost comparison. The main result of this study is the creation of a model that can provide input for a genetic Algorithm and can even be effective as a precise tool. This is pre-operation processing, as all processing can start with the knowledge of the optimum. The second result is a clear position-definition that is in-service and real-time.

#### **3. Materials and Methods**

The first step in creating a model is to clarify the notations and their content; therefore, a complete parameter exploration was carried out. The data set was defined with the help of practitioners. The basic data collected were grouped according to their roles in the model. Then, after selecting the parameters needed, it was important to capture their dimensions. In the second step, we created a parametric model of the warehouse and the factory parameters of the drone we wanted to use, improved with the data we had found. We created a manageable simplified model of motion and velocity (the simplification was within reason). The next step was to work out a practical positioning strategy. Then, we created a mathematical model of the operation. After that, the necessary items and environment for the optimization were selected.

During the testing of the GA method, we encountered some inconsistencies, which caused the model to undergo several manual tunings (the main reasons were the discrepancy between the factory and real drone parameters and the not necessarily smooth motion, which could have been due to several reasons, e.g., temperature and stock saturation problems). These inconsistencies, as they required relatively small modifications, were easily eliminated.

#### *3.1. Known Data*

In this section, we provide the mathematical model of the problem and the associated known and unknown data. The known data include the most important parameters of the warehouse and the drone: specifically, they were provided by the warehousing company, and the drone parameters are given in the technical specifications. Tables 1 and 2 summarize the data used for the model and provide the notations. We will use each index consistently for the same data, so we will also provide a table of indexes.

**Table 1.** Summary and notation of known data.



**Table 2.** Table of indexes.

The schematic structure of the shelving system is shown in Figures 1–4.

**Figure 1.** The floor plan of the warehouse with the dimensions. (Source: self-edited figure).

**Figure 2.** Shelf system details. (Source: self-edited figure).

**Figure 3.** Structure of a shelf. (Source: self-edited figure).

**Figure 4.** The shelving system to be tested. Black shelves indicate compartments that have already been checked. The green circle is the start position of the examination, the red square is the center of the corridor (reference point). The red circles are reference points. (Source: self-edited figure).

#### *3.2. Variable Data*

Search for the

$$S = S((k\_1i\_1j\_1), (k\_2i\_2j\_2), \dots, (k\_{2mn}i\_{2mn}j\_{2mn})) \tag{1}$$

series of compartments, where

$$(k\_{l\_1}i\_{l\_1}j\_{l\_1}) \neq (k\_{l\_2}i\_{l\_2}j\_{l\_2}) \text{ ha } l\_1 \neq l\_2. \tag{2}$$

The drone passes over all the compartments of the two shelving systems during operation. The series of compartments *S* describes the order of these compartments where the first parameter is the shelf system number (1 or 2) and the other two are the number of shelves (rows) in the shelf system and the number of compartments (columns) within it.

The drone needs to be charged from time to time so it may not be possible to process the two shelf systems with one charge (especially when the batteries are not fully charged from the start). For this reason, the series of compartments *S* must be divided into parts that can be processed with one charge, i.e., the drone contains compartments processed during operation. (Table 3 shows summary and notation of variable data). Thus, the processing time consists of three parts:



**Table 3.** Summary and notation of variable data.

The following correlations illustrate the abovementioned 3 points. The drone goes to the charger when it arrives at a compartment in which the available operating time for the previous compartment is still greater than Δ*t*, but in the case of the current compartment it is already less, i.e., the operating time falls below a predefined level.

In addition, it must be fulfilled that the sum of each subchain must be equal to the entire chain, since each compartment must be examined in *S* order.

Mark the series of compartment tests carried out simultaneously in a corridor

$$\mathcal{S} = \mathcal{S}((k\_p i\_p j\_p), \dots, (k\_q i\_q j\_q)) \subset \mathcal{S} \ (q \ge p) \quad q, p \in \mathbb{N}^+ \tag{3}$$

if the shelving system is to be changed, it will change either in the direction of *C* or in the direction of *C*.

Let *S*˘ denote the corridor change in *G*, i.e.,

$$\mathbf{G} = \begin{cases} \begin{pmatrix} k\_p, k\_{p+1}, 0 \end{pmatrix}, \text{ if it changes in direction } \mathbf{C} \\ \begin{pmatrix} k\_p, k\_{p+1}, -1 \end{pmatrix}, \text{ if it changes in direction } \overline{\mathbf{C}} \end{cases} \tag{4}$$

Let

$$\mathcal{S} = \mathcal{S}\left(\mathcal{S}\_1, \mathcal{G}\_1, \mathcal{S}\_2, \dots, \mathcal{G}\_{f-1}, \mathcal{S}\_f\right) \tag{5}$$

then, the execution time of the *S*´ series is from the charger to the charger

0 ≤ *Tw* − *Tf*

 *S*´ 

$$T\_f\left(\dot{\mathcal{S}}\right) = t\_{T, R\_{\hat{k}}\left(i\_p, j\_p\right)} + t\_P\left(\dot{\mathcal{S}}\right) + t\_{R\_{\hat{k}}\left(i\_q, j\_q\right), T} \tag{6}$$

≤ Δ*t* (7)

if the following is met

Let

$$\mathcal{L}\_{\omega} \mathcal{L}((l\_1 \dot{\imath}\_1 \dot{\imath}\_2)(k\_2 \dot{\imath}\_2 \dot{\imath}\_2) \qquad (l\_2 \dot{\imath}\_2 \dot{\imath}\_2 \dot{\imath}\_2 \dot{\imath}\_2)) = \sum\_{i=1}^{r} \mathcal{L}\_{\omega} \tag{69}$$

$$S = S((k\_1 i\_1 j\_1), (k\_2 i\_2 j\_2), \dots, (k\_{2mn} i\_{2mn} j\_{2mn})) = \sum\_{u=1} S\_u \tag{8}$$

where every *S*´ *u*

$$0 \le T\_w - T\_f(\mathcal{S}\_u) \le \Delta t \ (u = 1, \dots, r - 1). \tag{9}$$

The last one no longer needs to meet this condition.

The final report must also include a flight to the charger. For this reason, *S* must be supplemented with the **D** drone charger in the appropriate places. Its location will be determined by a later algorithm (Target Function Generation Algorithm). The function *Tf S*´ *u* represents the operating time of the drone from a charger to the next charger. The sum of these yields the total operating time.

Then the execution sequence is as follows:

$$S\_{\upsilon} = S\_{\upsilon}(\mathbf{D}, \boldsymbol{\xi}\_1, \mathbf{D}, \dots, \boldsymbol{\xi}\_r, \mathbf{D}). \tag{10}$$

Then the goal is

$$T = \sum\_{u=1}^{r} T\_f \left( \mathcal{S}\_u \right) \to \min! \tag{11}$$

If *Tf*(*Sv*) is minimal, then *Sv* is the optimal route

#### *3.3. The Speed of the Drone*

The speed of the drone is determined by knowing the speed of horizontal travel and the speed of vertical travel by projecting it on a diagonal path, as shown in Figure 5. The red arrow indicates velocity as a vector and its magnitude is the magnitude of the current velocity.

**Figure 5.** Drone speed (source: self-edited figure).

Ascension:

$$V\_t = \sqrt{v\_h^2 + v\_a^2} \tag{12}$$

$$
\alpha = \arctan \frac{\underline{v}}{\underline{x}}, \beta = \arctan \frac{v\_d}{v\_h} \tag{13}
$$

Descension:

$$V\_{\mathfrak{c}} = \sqrt{\upsilon\_h^2 + \upsilon\_d^2} \tag{14}$$

$$
\beta = \arctan \frac{v\_d}{v\_h} \tag{15}
$$

It follows that the ascension speed of the drone is:

$$V\_1(\mathbf{x}, \mathbf{y}) = V\_{\ell^\*} \cdot \cos(\beta - \pi) = V\_{\ell^\*} \cdot \cos\left(\arctan\frac{\upsilon\_d}{\upsilon\_h} - \arctan\frac{\mathbf{y}}{\pi}\right) \tag{16}$$

It follows that the descension speed of the drone is:

$$V\_2(x, y) = V\_c \cdot \cos\left(\arctan\frac{v\_d}{v\_h} - \arctan\frac{y}{x}\right) \tag{17}$$

Let *iv*, *jv* denote the indices of the compartment corresponding to the location of the drone, and *iv*+1, *jv*+<sup>1</sup> denote the indices of the next compartment to be examined by the drone.

If

$$\mathbf{x} = abs(j\_v - j\_{v+1}) \cdot \mathbf{h}\_h \tag{18}$$

$$y = \operatorname{abs}(i\_{\upsilon} - i\_{\upsilon + 1}) \cdot h\_{\epsilon}. \tag{19}$$

If *x* = 0, then

 $\text{If } i\_{\upsilon} - i\_{\upsilon+1} > 0$ 
$$V\_3(x, y) = V\_d \tag{20}$$

 $\text{If } i\_v - i\_{v+1} < 0$ 
$$V\_4(x, y) = V\_a$$

then it follows that the speed of the drone is:

$$V = V(\mathbf{x}, \mathbf{y}) = \begin{cases} j\_{\text{v}} - j\_{\text{v}+1} \ge 0 \text{ és } \mathbf{x} \ne 0 \text{ akkor } V\_2(\mathbf{x}, \mathbf{y})\\ j\_{\text{v}} - j\_{\text{v}+1} < 0 \text{ és } \mathbf{x} \ne 0 \text{ akkor } V\_1(\mathbf{x}, \mathbf{y})\\ j\_{\text{v}} - j\_{\text{v}+1} \ge 0 \text{ és } \mathbf{x} = 0 \text{ akkor } V\_3(\mathbf{x}, \mathbf{y})\\ j\_{\text{v}} - j\_{\text{v}+1} < 0 \text{ és } \mathbf{x} = 0 \text{ akkor } V\_4(\mathbf{x}, \mathbf{y}) \end{cases} \tag{22}$$

Thus, substituting (*iv* − *iv*+1)·*hh* for *x*, and substituting (*jv* − *jv*+1)·*he* for *y*, we obtain the following

$$\begin{array}{l} V\_1 = V\_1 \left( (j\_\upsilon - j\_{\upsilon+1}) \cdot h\_{l\iota} \left( i\_\upsilon - i\_{\upsilon+1} \right) \cdot h\_\iota \right) = \\ = V\_{\ell^\*} \cdot \cos \left( \arctan \left( \frac{\upsilon\_a}{\upsilon\_{l\iota}} \right) - \arctan \left( \frac{abs \left( i\_\upsilon - i\_{\upsilon+1} \right) \cdot h\_\iota}{abs \left( j\_\upsilon - j\_{\upsilon+1} \right) \cdot h\_\iota} \right) \right) \end{array} \tag{23}$$

$$\begin{split} \boldsymbol{V}\_{2} &= \boldsymbol{V}\_{2} \left( (\boldsymbol{i}\_{\upsilon} - \boldsymbol{i}\_{\upsilon+1}) \cdot \boldsymbol{h}\_{\mathrm{lt}} \, (\boldsymbol{j}\_{\upsilon} - \boldsymbol{j}\_{\upsilon+1}) \cdot \boldsymbol{h}\_{\mathrm{t}} \right) = \\ &= \boldsymbol{V}\_{d} \cdot \cos \left( \arctan \left( \frac{\boldsymbol{v}\_{d}}{\boldsymbol{v}\_{\mathrm{lt}}} \right) - \arctan \left( \frac{\mathrm{abs} \, (\boldsymbol{i}\_{\upsilon} - \boldsymbol{i}\_{\upsilon+1}) \cdot \boldsymbol{h}\_{\mathrm{t}}}{\mathrm{abs} \, (\boldsymbol{j}\_{\upsilon} - \boldsymbol{j}\_{\upsilon+1}) \cdot \boldsymbol{h}\_{\mathrm{lt}}} \right) \right) \end{split} \tag{24}$$

This depends on the horizontal and vertical travel distances.

#### *3.4. The Length and Time of Each Route*

The drone begins all processing by centering the two shelving systems, as shown in the green circle's position in Figure 6. This point denotes the center of the bottom row (which can be adjusted if necessary). The drone can take the best shot from the center of the compartments. The steps of processing include:


**Figure 6.** Route and guides. The green circle is the start position of the examination; the blue line is the route of drone; the red square is the center of the corridor (reference point). The red circles are referece points. (source: self-edited figure).

3.4.1. The Journey and Time from the Dock to the Starting Point It consists of three parts (Figure 6):


Step 1. Time to reach the starting position of the shelves Distance of compartment *Rk* = *Rk*(*i*, *j*) from charger:

$$d\_{T, \mathcal{R}\_k(i, j)} = d\_{T, \mathcal{C}} + d\_{\mathcal{C}, \mathcal{R}\_k}. \tag{25}$$

*Rk* = *Rk*(*i*, *j*) compartment access time from charger:

$$t\_{T,R\_k(i,j)} = t\_{T,C} + t\_{C,R\_{k'}} \tag{26}$$

where *dT*,*<sup>C</sup>* is the path between the docking station and the starting position of the middle guide path between the two shelves to the green square.

$$d\_{T, \mathbb{C}\_k} = \sqrt{\left(\frac{h\_\varepsilon}{2}\right)^2 + (c\_{kx} - d\_x)^2 + \left(c\_{ky} - d\_y\right)^2} \tag{27}$$

*tT*,*<sup>C</sup>* is the time required to reach the starting position of the middle guide path between the Dock and the two shelves:

$$t\_{T,C\_k} = \frac{l\_{l\_k}}{2V\_a} + \frac{\sqrt{(c\_{kx} - d\_x)^2 + \left(c\_{ky} - d\_y\right)^2}}{V\_h}.\tag{28}$$

Step 2. Distance and time from the starting position of the shelves to the starting compartment.

Distance and time from the center to the starting compartment:

$$d\_{\mathbb{C}\_k, \mathbb{R}\_k} = \sqrt{\left((j-1)\cdot h\_h\right)^2 + \left((i-1)\cdot h\_c\right)^2} \tag{29}$$

$$t\_{\mathbb{C}\_k, \mathbb{R}\_k} = \frac{d\_{\mathbb{C}\_k, \mathbb{R}\_k}}{V((j-1) \cdot h\_{h\prime}(i-1) \cdot h\_\varepsilon)} + \frac{90}{V\_r} + T\_{PH} \tag{30}$$

Step 3. The camera

Where <sup>90</sup> *Vr* is the camera rotation time.

The way to the starting position:

$$t\_{T,R\_k(i,j)} = t\_{T,C\_k} + t\_{C\_k,R\_k} = 1$$

$$t = \frac{h\_c}{2V\_a} + \frac{\sqrt{(c\_{kx} - d\_x)^2 + (c\_{ky} - d\_y)^2}}{V\_h} + \frac{\sqrt{((j-1) \cdot h\_h)^2 + ((i-1) \cdot h\_c)^2}}{V((j-1) \cdot h\_h)^2 (i-1) \cdot h\_c)} + \frac{90}{V\_r} + T\_{PH} \tag{31}$$

3.4.2. The Time to Travel to the Dock from a Given Point Can Be Determined Similarly This is the reverse of the previous since


Summarized as follows:

$$t\_{R\_k(ij),T} = t\_{R\_k\mathcal{E}} + t\_{C\_k,T} = $$

$$t = \frac{90}{V\_r} + \frac{\sqrt{\left(\left(j-1\right)\cdot h\_h\right)^2 + \left(\left(i-1\right)\cdot h\_\mathcal{E}\right)^2}}{V\left(\left(j-1\right)\cdot h\_{h\prime}\left(i-1\right)\cdot h\_\prime\right)} + \frac{h\_\mathcal{E}}{2V\_d} + \frac{\sqrt{\left(c\_{kx} - d\_x\right)^2 + \left(c\_{ky} - d\_y\right)^2}}{V\_h} \tag{32}$$

3.4.3. The Length and Time of a Point-To-Point Route

The path and time between the next two compartments are determined as shown in Figure 7. Due to the middle plane, it is sufficient to determine the distance between the two compartments by creating a Pythagorean theorem. It is enough to calculate the speed for the time and by dividing the distance we obtain the time. In addition, the camera can rotate if necessary (180◦). This can be determined by subtracting the position of the shelf system of the two compartments and then taking the absolute value. This is 0 if the 2 cells are in the same shelving system and 1 if the two compartments are in a separate shelving system. If the rotation time is multiplied by this value, it will either rotate or not as per our need. At the end, the time of photography is added.

**Figure 7.** Compartment to compartment. The green circle is the start position of the examination; the blue line is the route of drone; the red square is the center of the corridor (reference point). The red circles are referece points. (source: self-edited figure).

Let us take (*kviv jv*),(*kv*+1*iv*+<sup>1</sup> *jv*+1) as the two points. Their distance is:

$$d\left( (k\_{\upsilon}i\_{\upsilon}j\_{\upsilon})\_\iota (k\_{\upsilon}i\_{\upsilon+1}j\_{\upsilon+1}) \right) = \sqrt{\left( (j\_{\upsilon} - j\_{\upsilon+1}) \cdot h\_h \right)^2 + \left( (i\_{\upsilon} - i\_{\upsilon+1}) \cdot h\_c \right)^2} \tag{33}$$

Time spent on this route by the drone:

$$\begin{split} \frac{t\_{P\_1}((k\_{\overline{v}}i\_{\overline{v}}j\_{\overline{v}}),(k\_{\overline{v}+1}i\_{\overline{v}+1}j\_{\overline{v}+1}))}{\sqrt{((j\_{\overline{v}}-j\_{\overline{v}+1})\cdot h\_{\overline{h}})^2 + ((i\_{\overline{v}}-i\_{\overline{v}+1})\cdot h\_{\overline{c}})^2}} & \quad + \frac{ab\mathbf{s}(k\_{\overline{v}}-k\_{\overline{v}+1})\cdot 180}{V\_r} + T\_{PH} \end{split} \tag{34}$$

where

$$\frac{abs\,(k\_1 - k\_2) \cdot 180}{V\_r} \tag{35}$$

is the rotation of the camera.

#### 3.4.4. From One Compartment to Another Compartment

When considering change from one compartment in one row to another compartment in another row, we used the following steps (See Figure 8):


**Figure 8.** From compartment to compartment to another row. The green circle is the start position of the examination; the blue line is the route of drone; the red square is the center of the corridor (reference point). The red circles are referece points.

The situation is similar if the change is made at the other end of the shelving system. Let

$$a(j) = \begin{cases} j - 1, \text{ if } \mathbb{C} \text{ is the point under test} \\ n - j, \text{ if } \overline{\mathbb{C}} \text{ is the point under test} \end{cases} \tag{36}$$

$$\begin{aligned} T\_{R\_{k\_1}(i\_1, j\_1), R\_{k\_2}(i\_2, j\_2)}(\mathbb{C}) &= t\_{R\_{k\_1}\mathbb{C}\_{k\_1}} + t\_{\mathbb{C}\_{k\_1}\mathbb{C}\_{k\_2}} + t\_{\mathbb{C}\_{k\_1}\mathbb{R}\_{k\_1}} = \\ &= \frac{d\_{R\_{k\_1}\mathbb{C}\_{k\_1}}}{V(a(j\_1) \cdot h\_{h}(i\_1 - 1) \cdot h\_{\varepsilon})} + \frac{90}{V\_r} + \frac{\frac{h\_r}{\rho\_l} + \operatorname{abs}(k\_1 - k\_2) \cdot (h\_w + h\_d) + \frac{h\_c}{\rho\_l}}{v\_h} + \\ &+ \frac{d\_{\mathbb{C}\_{k\_2}, \mathbb{R}\_{k\_2}}}{V(a(j\_1) \cdot h\_{h}(i\_2 - 1) \cdot h\_{\varepsilon})} = \\ &= \frac{d\_{R\_{k\_1}, \mathbb{C}\_{k\_1}}}{V(a(j\_1) \cdot h\_{h\_{\varepsilon}}(i\_1 - 1) \cdot h\_{\varepsilon})} + \frac{h\_c + \operatorname{abs}(k\_1 - k\_2) \cdot (h\_w + h\_d)}{v\_h} + \\ &+ \frac{d\_{\mathbb{C}\_{k\_2}, \mathbb{R}\_{k\_2}}}{V(a(j\_2) \cdot h\_{h\_{\varepsilon}}(i\_2 - 1) \cdot h\_{\varepsilon})} + \frac{90}{V\_r} + T\_{PI} \end{aligned} \tag{37}$$

Similar to the optimal path, calculate the

$$T\_{R\_{k\_1}(i\_1, j\_1), R\_{k\_2}(i\_2, j\_2)}(\\\\\overline{\mathbf{C}}) = t\_{R\_{k\_1}, \mathbf{\mathcal{T}}\_{k\_1}} + t\_{\mathbf{\mathcal{T}}\_{k\_1}, \mathbf{\mathcal{T}}\_{k\_2}} + t\_{\mathbf{\mathcal{T}}\_{k\_1}, R\_{k\_1}} \tag{38}$$

Let

$$t\_{R\_{\hat{k}\_1}(i\_1;j\_1), R\_{\hat{k}\_2}(i\_2;j\_2)} = \min \left( T\_{R\_{\hat{k}\_1}(i\_1;j\_1), R\_{\hat{k}\_2}(i\_2;j\_2)}(\mathbb{C}), T\_{R\_{\hat{k}\_1}(i\_1;j\_1), R\_{\hat{k}\_2}(i\_2;j\_2)}(\overline{\mathbb{C}}) \right) \tag{39}$$

Taking into account that there may be two consecutive compartments to being examined within a row or between different rows, the function *tP* is composed of two parts, of which only the current one is always included. If there are two opposing shelving systems, then if *kv*+<sup>1</sup> = *kv* + 1 and *mod*(*kv*, 2) = 1, or *kv* = *kv*+<sup>1</sup> + 1 and *mod*(*kv*, 2) = 0, or *kv*+<sup>1</sup> = *kv*, then

$$t\_P((k\_\nu i\_\upsilon j\_\upsilon)\_\upsilon(k\_{\upsilon+1} i\_{\upsilon+1} j\_{\upsilon+1})) = t\_{P\_1}((k\_\nu i\_\upsilon j\_\upsilon)\_\prime, (k\_{\upsilon+1} i\_{\upsilon+1} j\_{\upsilon+1})) \tag{40}$$

otherwise

$$t\_P((k\_{\text{ $\boldsymbol{v}$ }}i\_{\text{ $\boldsymbol{v}$ }}),(k\_{\text{ $\boldsymbol{v}$ }+1}i\_{\text{ $\boldsymbol{v}$ }+1}j\_{\text{ $\boldsymbol{v}$ }+1})) = t\_{\mathbb{R}\_{k\_{\text{ $\boldsymbol{v}$ }}}(i\_{\text{ $\boldsymbol{v}$ }}j\_{\text{ $\boldsymbol{v}$ }})} \mathbb{R}\_{k\_{\text{ $\boldsymbol{v}$ }+1}}(i\_{\text{ $\boldsymbol{v}$ }+1}j\_{\text{ $\boldsymbol{v}$ }+1}) \tag{41}$$

Execution time of the *S*´ *u* sub-sequence (without docking):

$$t\_P(\boldsymbol{\dot{\xi}}\_\nu) = \sum\_{\upsilon=p}^{q-1} t\_P((k\_\upsilon i\_\upsilon j\_\upsilon)\_\prime (k\_{\upsilon+1} i\_{\upsilon+1} j\_{\upsilon+1})\_\prime). \tag{42}$$

$$T\_f\left(\dot{\mathcal{S}}\_u\right) = t\_{T, \mathbb{R}\_k(i, j)} + t\_P\left(\dot{\mathcal{S}}\_u\right) + t\_{\mathbb{R}\_k(i, j), T} \tag{43}$$

#### *3.5. The Constraints*

Based on the above, the sum of the compartment variables is exactly the sum of the number of cells from 0 to the number of cells −1.

Two conditions must be provided in the model's conditional framework:


This can be specified by two group conditions. The value of the variable cells must not be negative (since it is an entry sequence); therefore,

$$0 \le \varkappa\_{kij}.\tag{44}$$

The value of *xkij* can be clearly bounded from above:

$$
\omega\_{kij} \stackrel{\ast^\*}{=} r \cdot m \cdot n. \tag{45}
$$

Each number must be different if you are looking at different compartments:

$$k, k\_1 \in [1, \dots, r], \ i, i\_1 \in [1, \dots, m], j\_i, j\_1 \in [1, \dots, n]. \tag{46}$$

If (i.e., two compartments are different), then

$$\begin{aligned} 1000k\_1 + 100i\_1 + j\_1 &\neq 1000k\_2 + 100i\_2 + j\_2, \\ \varkappa\_{k\_1, i\_1, j\_1} &\neq \varkappa\_{k\_2, i\_2, j\_2}. \end{aligned} \tag{47}$$

In summary, the mathematical model

$$\begin{aligned} k, k\_1 &\in [1, \dots, r], \ i, j\_1 \in [1, \dots, m], j\_\prime \not{j\_1} \in [1, \dots, n] \\ &0 \le \mathbf{x}\_{kij} \le r \cdot m \cdot \mathbf{n} \end{aligned} $$

$$\begin{aligned} \mathbf{x}\_{kij} \ne \mathbf{x}\_{k\_2 j\_2 j\_2} &\text{if } 1000\mathbf{k} + 100\mathbf{i} + j \ne 1000\mathbf{k}\_1 + 100\mathbf{i}\_1 + j\_1 \text{.} \\ &T = \sum\_{u=1}^p T\_f(\mathcal{S}\_u) \to \min! \end{aligned} \tag{48}$$

#### **4. The Solution Method**

The computational model can be further simplified by

$$0 \le \varkappa\_{k\bar{i}j} \le \ r \cdot m \cdot n \tag{49}$$

The constraint can be replaced by the

$$0 \le \mathbf{x}\_{kij} \le m \mathbf{t} \mathbf{p} \cdot \mathbf{r} \cdot \mathbf{m} \cdot \mathbf{n} \tag{50}$$

(e.g., may be *mtp* = 100), and use the following routine as the nonlinear constraint. This is important for the genetic algorithm. Then, it is sufficient that all values are different.

*4.1. The Objective Generation Algorithm (OGA)*

Step 1.

Construct a random set of

$$S = S((k\_1i\_1j\_1), (k\_2i\_2j\_2), \dots, (k\_{2mn}i\_{2mn}j\_{2mn})) \tag{51}$$

compartments where

$$(k\_{l\_1}i\_{l\_1}j\_{l\_1}) \neq \ (k\_{l\_2}i\_{l\_2}j\_{l\_2}), \text{ if } l\_1 \neq l\_2. \tag{52}$$

Each compartment should be tested once and only once. Step 2. Let *u* = 0, *Sv* = ∅, and *g* = (*k*1*i*<sup>1</sup> *j*1) the first item in the series (if any). Step 3. If there are no more items in the series, continue from Step 5. If there are, then

$$\begin{array}{c} \mu := \mu + 1 \\ \mathring{S}\_{\mathfrak{U}} = \mathcal{Q} \end{array} \tag{53}$$

Step 4.

Add element *g* of the series to the subsequence.

*S*´ *<sup>u</sup>* := *S*´ *<sup>u</sup>* + {*g*} (54)

Calculate the value of

If

$$T\_w - T\_f \left(\dot{\mathcal{S}}\_u\right) > \Delta t \tag{56}$$

(55)

take the next element *g* and continue from Step 4.

If

$$T\_w - T\_f \left( \mathcal{S}\_u \right) \le 0 \tag{57}$$

i.e., the drone would not return with the addition of the new element, then

*Tf S*´ *u*

$$\begin{aligned} \mathcal{S}\_u &:= \mathcal{S}\_u - \{ \mathcal{g} \} \\ T &:= \; T + T\_f(\mathcal{S}\_u) \\ \mathcal{S}\_v &:= \mathcal{S}\_v + \{ \mathbf{D}\_r \mathcal{S}\_u \} \end{aligned} \tag{58}$$

Then proceed to Step 3.

If

$$0 \le T\_w - T\_f(\mathcal{S}\_u) \le \Delta t \tag{59}$$

Then take the next item *g* and continue from Step 3.

$$\begin{aligned} T &:= \; T + T\_f \left( \dot{\mathcal{S}}\_u \right) \\ \mathcal{S}\_v &:= \mathcal{S}\_v + \left\{ \mathbf{D}\_\prime \mathcal{S}\_u \right\} \end{aligned} \tag{60}$$

Step 5. Let

$$S\_v := S\_v + \{\mathbf{D}\}\tag{61}$$

At this point, the task list was completed with the charges, and the total processing time *T* was obtained.

#### *4.2. Determining the Optimal Route*

We used an evolutionary algorithm for the solution.

Step 1. Set the initial parameters (population number, number of steps, etc.). The formation of the initial population, i.e., each chromosome, will be one

$$S = S((k\_1i\_1j\_1), (k\_2i\_2j\_2), \dots, (k\_{2mn}i\_{2mn}j\_{2mn})) \tag{62}$$

series.

Step 2.

Run the algorithm by determining the fitness value for each individual in the population using the Target Function Generation Algorithm (OGA). This should also be built into factory algorithms.

Step 3.

After stopping the algorithm, we evaluate the result and program the obtained optimal solution (*Sv*) into the drone together with the necessary camera rotations and photography.

#### **5. Results**

The model and procedure described above leads to the following results.

The first result is in an enclosed double-bay warehouse that does not have a satellite positioning system nor an IoT positioning system; the inventory can be checked automatically using a mid-range drone. It follows that the above model generally handles the movement of the drone in the aisles of a double-bay warehouse between the charging station (dock) and the corresponding starting compartment, and the last compartment and the dock. The outlined procedure provides a near-optimal route-planning method that can scan as many compartments as possible relative to the operating time of the drone and ensures a safe return to the starting point. Since the warehouse solution does not require a complete optimal solution, a well-tested OGA device can clearly be used.

The results can be divided into two parts: in the first part, a route-optimizing model was created. The characteristics of the warehouse were considered, such as the unified double-shelf system structure and the resolution of a compartment up to three parts. The drone was able to take photos from the middle of the aisle, so there was no need to move out of the bisecting plane towards the shelving systems. This was a simplifying condition discovered during the study in the logistics decenter. The drone was not responsible for finding the QR codes associated with the compartments, nor was it responsible for evaluating whether there were any goods placed in it. All these tasks can be easily performed by employing software to evaluate the photos, which could further accelerate the process. Another possible option is to take a picture of two opposite compartments in a row by rotating the camera, provided that it happens in a shorter time. This is influenced by specific values of the parameters. The first result obtained is the route optimization; this should be the procedure carried out prior to the launch of the drone.

The second result is on-the-fly control: after launching the drone from the dock, it receives instructions where to go, how to get there, and what activities to perform there (if camera rotation and photography are needed). The drone is instructed in real time for the next task for each position, while (of course) internal energy level monitoring is also performed.

#### **6. Discussion**

During our research of the literature, we came across several solutions related to optimal road planning. Most of them solve the problem with the help of neural networks. In their paper, Andrey V. Gavrilov and Artem Lenskiy propose a model for a new biologically inspired mobile robot navigation system. The novelty of their work is the combination of short-term memory and online neural-network learning using the event history stored in this memory. The neural network is trained with a modified error backpropagation algorithm that uses the principle of reward and punishment when interacting with the environment [18]. The robot navigation mechanism is one of the most challenging research topics in mobile robots, which requires the robot to find the right path and travel from its current position to the target position without encountering obstacles. In their paper, Cheng and Cheng use the intelligent method of reinforcement learning to find a solution to the abovementioned problem. It considers the distances detected by a laser beam and the relative movement angle as the input of the neural network model, and the action posture of the robot as the output. This neural network model is trained by a deep Q-learning network (DQN) algorithm through positive and negative feedback rewards defined by task-specific learning goals. In this sense, the trained model helps the robot determine the appropriate steps to take in each state to safely reach the target without manual intervention. According to the results achieved on the simulation platform, the trained neural network model successfully moves the robot from a random starting point to a random destination, which proves the effectiveness of the DQN algorithm in the field of robotic navigation. [19] These models could also be applied to our problem, but it should be considered that the drone must not move in a plane, but in space. This would greatly complicate the solution. In fact, the solution we provide provides a simpler, more efficient solution. In particular, the application of neural networks would have required many examples, which we did not have at our disposal.

Another task from the research was to create an image processing and inventory analysis and a recording software that can be used for the task. The first task in Image Processing is to recognize the QR code itself and to identify the specific compartment. The identification of the compartment is provided by the routing model, and the route itself must be known to the processing application. Hence, automatically received photos are taken according to the route activities. The shelf system and the compartment where the photo was taken can be assigned to each photo by the application. This is a simple synchronization task.

The processing of the images is not very complicated, as most QR code readers are able to validate, and of course, interpret the data in the code. However, the fact that a compartment may be divided into several sub-compartments in a manner not previously known, so that several QR codes may have to be recognized, complicates the situation. It is also important that the QR code of another compartment should not be included in the photo and should not be handled by the recognizing application either. This can be easily prevented, as the shelf system is uniformly colored in most logistics centers—there were blue columns and red shelves in the warehouse we examined. With the help of these colors, the cells can be delimited, and the division can be determined based on the expected size by linear image processing. In other words, based on the image, one can see how many QR codes to search for, which can be easily carried out with the QR code search procedure. (Focusing was included into the shooting time of the drone).

The final processing application checks whether the actual unit load in the given compartment meets the requirements based on the position and QR code and indicates the compliance or deviations to the warehouse operator accordingly.

#### **7. Conclusions**

In this paper, we presented a model of a drone-driven inventory control solution. As a result of the method and with the help of the implemented application, the time of inventory testing can be greatly reduced, and the task can be performed with the help of an average human resource. A task that had once taken several days to complete can now be completed in about two days with this solution. In addition, the platform solution currently used, which endangers personal safety, can be avoided. It is obvious that in the case of the entire warehouse, the use of several drones can further reduce this time, and erroneous picking due to human errors can be eliminated by frequent inspections, so that the use of drones can also reduce the resulting losses.

Determining location based on QR code recognition is not a viable option, and one of the best known and most widely used methods is the Viola–Jones object detection framework. The mentioned framework provides an efficient way to focus the detection process on specific parts of the image. Based on our solution design, the initial state already shows a possible solution, and the intermediate states will always be possible solutions and a good GA procedure is sufficient to apply a parameter reduction solution. In reality, there are many "good" paths to choose from. However, our tests show that the path depends significantly on the relationship between the underlying data.

At the beginning of our article, a full parameter exploration was performed. The range of data was defined with the help of practitioners, and then the basic collected data were grouped according to their role. Once the required parameters had been selected, it was important to record their dimensions. Next, we created a parametric model of the warehouse and the factory parameters of the drone to be used, improved with the data we had found. We created a manageable, simplified motion and speed. The next step was to develop a practical positioning strategy. Then, we created a mathematical model of the operation and assigned the theorems and determined environment needed for optimization.

**Author Contributions:** Conceptualization, M.G.; methodology, M.G.; writing—review and editing, L.S.; visualization, L.R.; investigation, J.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research project was funded by Budapest Business School Research Fund.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research project was carried out within the framework of Centre of Excellence for Future Value Chains of Budapest Business School.

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

#### **References**

