Enhancing Manual Order Picking through a New Metaheuristic, Based on Particle Swarm Optimization

Abstract

This paper proposes a new metaheuristic algorithm called Particle Swarm-based picking time minimization (Pkt_PSO), ideated for picking time minimization in manual warehouses. As the name suggests, Pkt_PSO is inspired by Particle Swarm Optimization (PSO), and it is specifically designed to minimize the picking time in order case picking contexts. To assess the quality and the robustness of Pkt_PSO, it is compared to five alternative algorithms used as benchmarks. The comparisons are made in nine different scenarios obtained by changing the layout of the warehouse and the length of the picking list. The results of the analysis show that Pkt_PSO has a slower convergence rate and suffers less of early stagnation in local minima; this ensures a more extensive and accurate exploration of the solution space. In fact, the solutions provided by Pkt_PSO are always better (or at least comparable) to the ones found by the benchmarks, both in terms of quality (closeness to the overall best) and reliability (frequency with which the best solution is found). Clearly, as more solutions are explored, the computational time of Pkt_PSO is longer, but it remains compatible with the operational needs of most practical applications.

1. Introduction

Warehousing is of paramount importance in supply chains and logistics, as a well-designed warehousing policy leads to benefits within the company and across the entire supply chain [1,2]. Most of the time, high operating performance is achieved using highly efficient Automated Storage and Retrieval Systems (AS/RS) [3,4]. However, in case of small-item storage and/or spare-parts handling, manual warehouses are still more efficient and represent the standard solution [5,6]. In this case, order picking (the problem of collecting a set of items in a warehouse) is an essential issue of order fulfilment, with an immediate impact on customer satisfaction, business reputation and profitability. Picking is generally a major bottleneck in logistics and one of the most expensive and labour-intensive logistic activities, which is why improving this activity is considered a top priority for many companies [7].

When the main operations concern pallet decomposition, short-period storage, pallet recomposition, and shrink wrapping of new pallets, a picker-to-goods strategy is generally adopted [8,9]. Accordingly, the picker moves in the warehouse to collect the required number of cases of each Stock Keeping Unit (SKU) included in the picking list. As the travelling time represents the largest part of the time required to complete an order, many picking strategies have been proposed to reduce the distance covered by the picker. Among them, batch picking, zone picking, wave picking, voice picking and vision picking are certainly the most studied [10,11]. However, all the above-mentioned strategies need specific operating conditions and/or costly equipment. Batch and wave picking assume that items can be sorted after retrieval [12], zone picking has proven to be efficient only if a very high number of SKUs are stored in the warehouse, while vision and voice picking require a computerized control system and expensive augmented reality glasses or Bluetooth headphones [13].

Due to these reasons, the basic picker to goods strategy, is still the best choice in many industrial contexts [14] and it is the most widespread policy in the context of small and spare-parts handling [15]. According to this policy, the picker is in charge of completing an entire order and, to do so, he/she can move freely in each area of the warehouse to collect the cases (of the SKUs indicated in the picking list) needed to form a pallet. Its implementation is simple and cheap, as it does not require any expensive equipment or additional operations such as batch composition or sorting. Operationally, the only decision concerns the routing that the picker must follow, that is, the order in which storage locations must be visited to complete the picking list. This sequence, in fact, determines the distance covered by the picker and the total travelling time. Optimizing the routing assigned to the picker is, thus, of paramount importance to improve productivity and to cut costs. Not surprisingly, literature in this field is rather extensive, and many interesting and efficient metaheuristics have already been introduced to minimize the total travelling time.

To the authors’ knowledge, the last and most efficient framework is the one proposed by De Santis et al. [9], which is also characterized by extraordinary speed of convergence. Although this feature is desirable, it also suggests that most of the provided solutions could be points of local minimum. Therefore, we believe that there is still room for improvement, and to test this hypothesis we propose a novel metaheuristic called Pkt_PSO. As the name suggests (where Pkt stands for picking time), Pkt_PSO is an evolutionary metaheuristic inspired by the well-known Particle Swarm Optimization (PSO). However, to exploit the peculiar structure of the considered problem, many features of the standard PSO have been readapted and modified. Making major modifications to the PSO is indeed a common practice, and it is generally recommended to exploit the potentialities of this approach (see for example [16,17]) This readaptation is, therefore, the main contribution of the paper. The paper is also novel because, to the best of our knowledge, PSO has never been applied to the order-picking domain, although some recent works have shown its effectiveness and flexibility in solving very difficult-to-tackle tasks, such as Bayesian network training [18,19] and different instances of the Traveling Salesman Problem (TSP) [20,21]. The latter is a classic operating research problem very similar to the one herein considered. Assessing the operational benefits that can be achieved using a PSO-based framework is the second main contribution of the paper.

The rest of the paper is organized as follows. Section 2 presents a brief introduction on the state-of-the-art approaches for travelling time minimization in manual warehouses. In the same section, literature gaps are highlighted and the research questions addressed in this paper are clearly defined. For the sake of clarity, Section 3 pinpoints the main features of the standard PSO and Section 4 formalizes the problem and its possible formulation as a TSP. The proposed algorithm (Pkt_PSO) is introduced in Section 5, where all the new procedures and operators developed to exploit the peculiar structure of the considered problem are clearly explained. In view of an experimental section used to validate Pkt_PSO, its operating parameters are fine tuned in Section 6. Next, in Section 7 its performance is tested in nine different scenarios (obtained changing the layout of the warehouse and the length of the picking list) and compared to that of five benchmarks. One of the benchmarks [9] is a top-of-the-class metaheuristic specifically designed for the picking problem in manual warehouses, two [22,23] are outstanding PSO-based approaches for the TSP, while the last two are the well-known Full Depth 2-Opt and the nearest neighbour constructive procedure. Finally, conclusions and directions for future works are drawn in Section 8.

For the sake of clarity, all the abbreviations used in the paper are listed at the end of the paper in the “Abbreviation Section”.

2. Picker-to-Goods Strategy, a Brief Introduction and Literature Review

2.1. The Basic Picker-to-Goods Strategy

The functioning of the basic picker-to-goods strategy is very simple. Each picker receives one picking list at a time and he or she must collect all the SKUs specified therein. Specifically, the picking list includes all the SKUs needed to form an entire pallet (as required by the end customer) and indicates, for each one of them, the number of cases to be picked and their exact locations in the warehouse. The picker, generally onboard of an electric order picker, is free to move in all the aisles of the warehouse and collects cases stored at the ground level. Top racks, in fact, are used to store full pallets that will be lowered to the ground (with forklifts) whenever necessary.

The implementation of this strategy is straightforward: given a set of storage locations, the only decision concerns the order in which they must be visited. This sequence, in fact, determines the routing followed by the picker, which is the only element that can be optimized to increase productivity. Whereas the travelling time highly depends on the routing, all the other times (e.g., pick up time of a case, loading and unloading of the order picker, etc.) are constant and, therefore, irrelevant. It is, therefore, an NP-hard combinatorial problem belonging to the family of the Travelling Salesman Problem (TSP) [24,25]. Even in case of order picking with additional constraints, such as one-way aisles [26], limited handling capacity of the pickers [27], partial overlapping of picking and refilling operations, narrow aisles that can be crossed by a single picker at most [28] or allocation of the same items in multiple locations [29], the task can be reduced to a routing optimization problem with the objective of makespan minimization.

2.2. Solution Approaches Proposed in the Literature to Improve Productivity

Solutions proposed in the literature are manyfold and very differentiated but, as for many other operational research problems, the solution approaches can be classified into two broad categories, namely, exact and heuristic methods. Due to the NP-hard nature of the problem, exact methods exist only for specific or simple warehouse configurations. Just to give an idea, the exact model proposed by Ratliff et al. [30], probably the first one published on the subject, is tailored to a one-block warehouse with 50 aisles and a picking list of 10 locations. Although the problem had such as small size, some minutes of CPU time were needed to find the global optimal solution on a mainframe computer of that time. Another good example of exact models can be found in [31], where two models (based on mixed-integer programming and dynamic programming, respectively), are used to solve the order-picking problem for a two or three cross-aisles warehouse.

Anyhow, even using modern processors, for most real applications the computational time needed to execute an exhaustive approach is too high to be of any practical help. Hence, over the years, exact methods have gradually been ignored, to be replaced by heuristics and metaheuristics capable of providing very good and feasible solutions (or pseudo-optimal solutions) in a reasonably acceptable computational time. A clear discussion about the superiority of heuristic methods over exact ones can be found in [32], where it is shown that heuristics are more flexible and versatile and that can be easily readapted to different layouts and/or to incorporate additional objectives and constraints. In virtue of this, to give an idea of the different trends in heuristics and metaheuristics for the picking problem in manual warehouse, a selected set of interesting algorithms is presented next. For a more comprehensive review, the interested readers are referred to the work by Masae et al. [33].

Some authors studied the effect of different layouts on routing optimization, as in [34], where three heuristics for optimal routing in warehouses with more than two cross aisles are compared, and in [35], where a fishbone layout is considered. In addition to different layouts, in [36], economic and ergonomics issues are also considered as further objectives. More recently, heuristics have been integrated with discrete event simulations, as in [37], that introduce an algorithm for dynamic order-picking problems, in which the picker might receive instructions for new picking requests while other ones are still in progress. A concurrent simulation-based design of experiments is used in [38] to optimize picking activities and to assess the effect of a combination of operating parameters (layout, throughput, size, order pickers, etc.) on overall performance.

Most of these studies assume a random storage assignment policy, but some authors have also considered dedicated assignment policies. Among them, Caron et al. [39] evaluated two simple routing heuristics in a manual single cross-aisle warehouse, with an allocation strategy based on the Cube per Order Index. Petersen [40] carried out a detailed analysis of variance to examine possible interdependencies among routing policies, warehouse layout and depot locations under different operating conditions. The same author also addressed the issue of optimizing the routing strategy in a class-based storage system [41]. In [42], the correlation between the number of aisles, picking list size, and path length in a class-based storage environment is studied, and in [43], the optimal size of ABC zone is determined for different layouts.

In conclusion, the present work fits into this line of research, as it explores the application of metaheuristics to solve the order-picking problem. The reader who is interested in a complete review of the literature with gap analysis and future research directions can refer to over 70 papers published in Scopus; the last two [44,45] published during 2023 are recommended.

2.3. Gap Analysis and Research Questions

As discussed, the use of metaheuristics allows for the extension of non-exact applications that are operationally more attractive. Indeed, pushing the application to more realistic contexts, the TSP that models an order-picking problem becomes intractable in exact terms. As mentioned above, the existing literature on the subject is rather extensive and many interesting and efficient metaheuristics have already been introduced. To the authors’ knowledge, the last and most efficient framework is the one proposed by De Santis et al. [9], which exploits:

• the famous Floyd–Warshall algorithm to solve the all-pair shortest path problem [46,47],
• the Ant Colony Optimization (ACO) to provide an effective solution in a very short time.

Obtained solutions have proved to be superior, or at least equal, to those provided by previous algorithms. However, its extraordinary speed of convergence suggests that most of the times the obtained solution could be a point of local minimum. Therefore, we believe that there is still room for improvement, and to test this hypothesis, we investigate the use of a new metaheuristics inspired by the classic PSO. The rational of this choice can be traced to the fact that, rather recently, some works have demonstrated their good applicability to the TSP (see among others [20,21]). Hence, our main goal is to investigate possible advantages and/or disadvantages of the PSO when applied to this specific instance of the TSP and, most importantly, to modify its underlying structure to exploit the specific structure of the problem at hand.

We make clear from the outset that high-performance improvements are unlikely to be achieved, mainly because of the very constrained structure of the TSP underlying a warehouse organized in blocks and divided by aisles, such as the one considered in this study. Because of this binding layout structure and considering that many SKUs could be found in the same aisle, it is reasonable to expect the existence of many solutions very close in terms of makespan (i.e., point of local minimum). Testing this hypothesis, and thus understanding how sensible it is to further explore optimization techniques, is the second objective of this work.

3. Some Background on Particle Swarm Optimization

PSO is a famous metaheuristic optimization technique inspired by flocks of birds. Its ideation dates to 1995 when Eberhart and Kennedy proposed PSO as an optimization technique for continuous functions [48]. Since then, as stated in the notable works by Banks et al. [49,50] and in the comprehensive literature review by Jain et al. [51], PSO has been widely applied in many fields of operational research both for continuous and discrete problems and in closed-loop supply chain network optimisation [52]. PSO is an evolutionary and population-based metaheuristic that exploits a population of N candidate solutions, which are evolved at each step of the algorithm. For the sake of clarity, the basic notation used to describe PSO is shown in

Table 1. PSO Parameters.

Symbol	Definition
$i$	A generic particle of the swarm.
$C$	A positive constant.
$gbest$(t)	The overall best solution found until iteration $t$.
$n$	The dimension of the solution space.
$N$	The number of particles of the swarm.
$p_i\left(t\right)$	A n-dimensional vector that codifies the position of particle i at iteration t, i.e., a feasible solution of the problem.
$pbes{t}_i\left(t\right)$	The best solution found by particle $i$ until iteration $t$.
$t$	The current iteration or epoch of the algorithm.
$u$	A generic random number uniformly distributed in the range [0, 1].
$v_i\left(t\right)$	A n-dimensional vector that codifies the velocity of particle i at iteration t.
$\omega$	The inertia factor, a random number uniformly distributed in the range [0, 1].
$\chi$	The constraining factor, a constant with a value between $\left(0, 1\right]$.
$f\left(p\right)$	The objective function of the problem.

.

The whole population is referred to as a swarm, while each solution ($i=1,\dots ,N$) of the swarm is referred to as a particle. At each step (or iteration) $t$, each particle moves in the solution space to generate a new and potentially improved solution. Specifically, its movement is formally defined by Equation (1):

$$p_i\left(t+1\right)=p_i\left(t\right)+ v_i\left(t\right)$$

(1)

where $v_i\left(t\right)$ and $p_i\left(t\right)$ are, respectively, the velocity and the position of particle $i$ at iteration $t$. Both $v_i\left(t\right)$ and $p_i\left(t\right)$ are vectors with several components equal to the dimensions of the solution space: the position $p_i\left(t\right)$ corresponds to the current solution coded by particle $i$, and the velocity $v_i\left(t\right)$ defines how far it will move at the next step.

At each iteration, the speed of each particle is updated as in Equation (2).

$$v_i\left(t+1\right)=v_i\left(t\right)+C_1{u}_{1,i}^t\left(p_i\left(t\right)-pbes{t}_i\left(t\right)\right)+C_2{u}_{2,i}^t\left(p_i\left(t\right)-gbest\left(t\right)\right)$$

(2)

where:

• $gbest\left(t\right)$, the best solution found up to iteration $t$ by the whole swarm.
• $pbes{t}_i\left(t\right)$, the best solution found up to iteration $t$ by particle $i$.
• $u_{1,i}^t$ and $u_{2,i}^t$ are two random numbers uniformly distributed in the range [0, 1]. At each iteration $t$ they are generated for each particle $i$.
• $C_1$ and $C_2$ are fixed parameters.

According to Equation (2), each particle tries to keep moving following the direction coded by $v_i\left(t\right)$ that, for this reason, is also referred to as the particle’s intention. However, the trajectory of the particle is subjected to small deviations that are related to its experience and to the knowledge (or culture) of the whole storm. Experience and culture are measured as the distance from the current position $p_i\left(t\right)$ and to the position of $gbest\left(t\right)$ and $pbes{t}_i\left(t\right)$, respectively. The entity of both deviations is random but limited by the two constants C. In this way, the algorithm balances the exploration of new solutions and the exploitation of the best solutions found so far. When $C_1=C_2=0$, all particles fly keeping their current speed and direction until they finally hit the boundaries of the search space. If $C_1>0 \mathrm{and} C_2=0$, particles behave as if they were independent individuals, whereas if $C_1=0 \mathrm{and} C_2=0$, they are all attracted by the same point. A c choice common choice is to set $C_1=C_2=2$. In this way, all particles are attracted towards the midpoint between $gbest\left(t\right)$ and $pbes{t}_i\left(t\right)$ and the strength of this attraction is not negligible (as it would be with $C_1=C_2\ll 1)$.

To summarize, the standard PSO works as follows:

• At every iteration $t,$ each particle $i$ moves from $p_i\left(t\right)$ to $p_i\left(t+1\right)$ according to Equation (1).
• If the new solution $p_i\left(t+1\right)$ is better than $pbes{t}_i\left(t\right)$, the current best solution found by $i$ is updated accordingly.
• When all the swarm has moved, $gbest\left(t\right)$ is updated and the velocity of each particle $i$ is updated using Equation (2).
• The process is repeated until a stopping criterion is met.

For the sake of clarity, the pseudocode is shown in Figure 1, where it is assumed that the PSO ends when the number of iterations without improvement exceeds a predefined threshold value. The pseudocode is written using a Python 3.9© style, and so, indentation is not only for readability but it also defines different blocks of code.

At present, the original PSO has been modified and extended in many ways. A significant contribution is that of [53], where a new parameter called inertia is introduced, and the updating procedure is modified as in Equation (3).

$$v_i\left(t+1\right)={\omega }_i^t·v_i\left(t\right)+C_{1}r\left(p_i\left(t\right)-pbes{t}_i\left(t\right)\right)+C_{2}r\left(p_i\left(t\right)-gbest\left(t\right)\right)$$

(3)

where the inertia ${\omega }_i^t$ is a random number uniformly distributed on the interval [0, 1] generated for each particle $i$ at each iteration $t$. Hence, the scope of $\omega$ is to mitigate the effect of particles’ intention, and in many applications, its value is progressively reduced as the number of iterations increases [53]. The resulting effect is to force the exploration of a wider area of the solution space during the first runs of the algorithm and to exploit the local optima (sort of neighbour search) when the procedure is almost at the end.

As an alternative, [54] proposed updating the velocity, as in Equation (4).

$$v_i\left(t+1\right)=\chi \left[v_i\left(t\right)+C_{1}r\left(p_i\left(t\right)-pbes{t}_i\left(t\right)\right)+C_{2}r\left(p_i\left(t\right)-gbest\left(t\right)\right)\right]$$

(4)

where $\chi$ is a constant value (in the range [0, 1]) called constriction factor that should ensure better reliability and greater stability of the results. As for the inertia, to assure convergence, some authors have also proposed a time dependent constriction factor, as in [55] where, as the particle becomes closer to the global minimum, a lower value of the constriction factor is used which helps in stabilizing the algorithm with fast convergence. We anticipate that neither the inertia nor the constriction factor will be used in the metaheuristic here proposed.

All the above-mentioned approaches work well when the objective function is continuous, and it is possible to introduce a velocity vector $v_i$ on the solution domain. When, as in most industrial settings, problems are binary and/or discrete, their industrial relevance becomes rather scarce. For this reason, many scholars have proposed discretized versions of the PSO. For instance, [56] proposed using a sigmoid function to convert continuous variables into binary ones, whereas [57] suggested working with continuous variables and rounding off the obtained solution to the nearest integer.

Relatively to the discrete domain, PSO has recently been proposed as an effective way to tackle NP-hard routing problems such as the TSP. Among the first contributions, we can cite the relevant works by [58,59]. However, the real burst of interest originated from [20,21], who showed the superiority of PSO in tackling these NP-hard problems. The point to note is that, when operating over a discrete domain, to outperform other metaheuristics, PSO must be subject to major changes that can make it distinctive from the standard and continuous version. Hence, the main question is how the concept of speed can be redefined so that the distance between two discrete solutions can be measured. In this regard, two remarkable contributions were proposed by [22,60]. Specifically, letting $S_1$ and $S_2$ be two sequences representing the order with which $n$ nodes must be visited, both authors expressed speed as a function of the number of swaps needed to convert $S_1$ to $S_2$. Although the results are rather good, the performances cannot be compared with those obtainable with PSO implementations where the concept of speed is completely overturned, as in the works by [23,61]. For a more in-depth description of the latter approach, please see Appendix A, which describes the algorithm by [22,23] and other approaches used as benchmarks.

4. Problem Definition and Mathematical Formulation

4.1. Problem Assumptions

In general terms, the order-picking problem can be defined as “the problem of collecting a set of SKUs in a manual warehouse in a minimum amount of time” [31]. In the present paper, we contextualize the problem using the following assumptions, which hold in most industrial settings.

• The set of SKUs to be collected is indicated in a picking list, which also reports the quantity, expressed in the number of cases, to be picked up.
• Each SKU is stored at the ground level and in a dedicated location in the warehouse. Hence, the picking list univocally defines the locations that must be visited to fulfil an order.
• The aisles of the warehouses are bidirectional and are large enough to avoid any traffic problems.
• There is a single Input Output (I/O) point, from which each tour starts and ends.
• The picker receives a picking list at a time and must collect all the cases indicated in it with a single tour (i.e., the handling capacity of the picker is assumed to be limitless).
• When the picker moves from one location to the next one, he or she proceeds at a constant speed, moving approximately at the centre of an aisle.
• When the picker reaches a picking location, he or she moves from the centre of the aisle to the proximity of the shelves where SKUs are stored. The time needed for this traversal movement is constant.
• The time needed to pick up a case is constant, and it does not depend on the shape, dimension, or weight of the case.

4.2. Problem Formulation and Graph Representation

It is easy to see that the total travelled time is the only relevant quantity of the problem; all the other times are constant (see assumptions 6 and 7) and can be excluded from the optimization procedure. Also, since the travelling time depends on the sequence with which the picker collects the SKUs, the goal is to find the optimal sequence (or optimal picking tour). It is, therefore, a permutational problem that can be reduced to a standard TSP.

To this aim:

• The layout must be schematized as a bidirectional graph, with nodes representing the I/O point, the storage locations, and the access points to the aisles. Considering that each aisle serves two counterposed shelves, the total number of nodes is $\mathcal{N}=\left(0.5\mathcal{L}+2\mathcal{A}+1\right)$, where $\mathcal{L}$ is the number of storage locations and $\mathcal{A}$ is the number of aisles. Furthermore, due to the constrained layout of the warehouse, (based on horizontal and vertical aisles), each node is linked with two or at most with three other nodes. The only exception is the I/O point that is linked with only one node. For the sake of clarity, an example of graph representation is shown in Figure 2, which is relative to a warehouse with 24 locations (from A to X) and three horizontal aisles. The total number of nodes is therefore (24/2 + 2 × 3 + 1) = 19. Also, note that node zero represents the I/O point and that nodes inside the aisles serve two locations each (for instance node two serves both location S and location T).
• The shortest path connecting each pair of nodes (including the I/O point) must be computed and saved in a triangular distance matrix $D$. This can easily be carried out in a polynomial time, using the well-known Floyd–Warshall algorithm [47,48].
• Given a picking list with $M$ SKUs, a bidirectional and fully connected network with $\left(M+1\right)$ nodes and $M·\left(M+1\right)$ edges must be created. The nodes correspond to the I/O point and to the $M$ locations that must be visited, where the edges codify the distances among them. These distances are read in the distance matrix $D$, i.e., the length of the edge connecting node $i$ to node $j$ equals $D\left[i,j\right] = d_{ij}$. For instance, if the picking list was {S, R, V}, the network of Figure 3 should be built, where it is supposed that the distance between two adjacent nodes equals 2 and that between node zero and node one equals 1.
• The optimal sequence is finally found in solving a TSP (as in Equation (5)) formulated for the network generated by the picking list.
• OF:

$$\min \left(\mathrm{z}\right)={{\displaystyle \sum }}_{i=1}^{M+1}{{\displaystyle \sum }}_{j=1}^{M+1}{x}_{ij}{d}_{ij}$$

(5a)

• ST:

$${{\displaystyle \sum }}_{i=1}^{M+1}{x}_{ij}=1 \forall j\in \left\{1, 2, \dots , \left(M+1\right)\right\}$$

(5b)

$${{\displaystyle \sum }}_{j=1}^{M+1}{x}_{ij}=1 \forall i\in \left\{1, 2, \dots , \left(M+1\right)\right\}$$

(5c)

$${{\displaystyle \sum }}_{i,j\in S}{x}_{ij}\le \left|S\right|-1 \forall S$$

(5d)

$$x_{ij}\in \left\{0,1\right\} \forall i\in \left\{1, 2, \dots , \left(M+1\right)\right\}, \forall j\in \left\{1, 2, \dots , \left(M+1\right)\right\}$$

(5e)

where

• $x_{ij}$ is a binary decision variable that equals one if the picker moves from location i to location j, and it is zero otherwise.
• $d_{ij}$ is the distance between location i to location j.
• $M$ is the number of SKUs in the picking list.
• $\left(M+1\right)$ is the number of nodes in the network ($M$ locations for $M$ SKUs and one I/O point).
• $S$ is a generic set of nodes, and |S| denotes its cardinality. Since the cardinality |S| ranges from 2 to $M$ there are $2^M$ possible set S.

The objective (Equation (5a)) consists of finding the minimum travelling distance, and the constraints make sure that each location of the picking list is visited only once (i.e., constraints (5b,c)) and that subtours are avoided (i.e., constraint (5d)). For an additional insight of this integer linear programming model, the interested reader is referred to the extensive literature review on the TSP (see for example [24]).

Continuing with the example, solving the optimization problem, the optimal tour {I/O → P → V → S → I/O} is obtained, for a total distance of 18. We also note that the actual path followed by the picker can be reconstructed, starting from the optimal tour, using the Floyd–Warhsall algorithm in reverse. In this way, the following path is finally obtained: {0, 1, 12, 11, 10, 11, 12, 1, 2, 1, 0}.

5. Pkt_PSO, the Proposed Metaheuristic

Since the TPS is a well-known NP-hard problem, a novel metaheuristic is herein proposed to solve the order-picking problem. The metaheuristic takes inspiration from the PSO and so, from here on, it will be referd to as the Pkt_PSO.

5.1. Notation

As discussed in Section 4, the problem input is a picking list $\mathcal{P}$. We denote the picking list as $\mathcal{P}=\left\{A,B,C,D,E,F\right\}$, where the uppercase letters are the $M$ locations to be visited. Locations are not sorted or can be sorted in lexicographical order just for convenience.

Similarly, a feasible solution is a sorted array of $\left(M+2\right)$ elements, which defines a picking tour. Since a tour always starts and ends at the I/O point, the first and last elements can be removed, and the tour can be fully represented with an array of only $M$ elements. We refer to this array as the picking tour and we denote it as $p=〈h_1,h_2,\dots ,h_M〉$, where each element $p\left[j\right]= h_j$ identifies the $j$-th location that must be visited by the picker. Consider, for instance, a picking list with six codes $\mathcal{P}=\left\{A,B, C,D,E,F\right\},$ and a picking tour $p=〈B,A,C,F,D,E〉$. In this case, the first stop is at location $B$, the second one is at location $A,$ and so on until location $E$ is reached and the picker completes the tour (i.e., $p\left[1\right]= h_1=B, p\left[2\right]= h_2=A, \dots ,p\left[6\right]= h_6=E)$.

The objective function evaluated in a specific tour is therefore computed as in Equation (6):

$$f\left(p\right)=d_{\left(IO, h_1\right)}+{\displaystyle \sum }_{j=1}^M{d}_{\left(h_j, h_{j+1}\right)}+d_{\left(h_M, IO\right)}$$

(6)

where $d_{\left(x,y\right)}$ is the minimal distance between location x and location y and IO is the input–output point.

For the sake of clarity, before delving into the details of the Pkt_PSO, the notation of the problem is fully introduced in

Table 2. Variables related to the picking problem and to the layout of the warehouse.

Symbol	Definition
$\mathcal{P}$	The picking list $\mathcal{P}=\left\{A,B,C,\dots \right\}$ containing all the storage locations that must be visited by the picker to complete an order. Locations in $\mathcal{P}$ are indicated with an uppercase letter and they are not sorted; they may be sorted in lexicographical order just for convenience.
$p$	The picking tour $p=〈h_1,h_2,\dots ,h_M〉$ containing all the locations of $\mathcal{P}$ sorted in the order that the picker must follow during the tour. Please note that locations in $p$ are denoted with the letter $h$.
$j$	The positional index $j=1,\dots ,M$ that specifies at which point of the tour a location is visited.
$h_j$	The j-th location to be visited, i.e., $p\left[j\right] = h_j$.
$d_{x,y}$	The minimum distance between location 𝑥 and location $y$, with $x,y\in \mathcal{P}$. The minimum distance is always obtained using the Floyd–Warshall algorithm.
${\mathcal{S}}_{x,y}$	The set of all the storage locations $l$ located within the shortest path linking $x$ and $y$. That is, to move from $x$ to $y$, the picker also passes in front of all locations contained in ${\mathcal{S}}_{x,y}$.
$s$	The cardinality of the set ${\mathcal{S}}_{x,y}$, i.e., the number of locations $l$ in ${\mathcal{S}}_{x,y}$

and

Table 3. Variables related to the parameters of the algorithm.

Symbol	Definition
$p_i\left(t\right)$	A completed picking tour (solution) generated for particle $i$ at iteration $t$. With “complete” we mean that $p_i\left(t\right)$ is made of $M$ locations.
$p_i^m\left(t\right)$	A partial picking tour generated for particle $i$ at iteration $t.$ The tour is partial because it is made of $m$ locations with $m<M$.
$T$	The maximum number of iterations before the algorithm stops.
$greedy$	A greedy solution shared by all the $N$ particles of the swarm.
$rnd\_intentio{n}_i\left(t\right)$	A random solution generated for particle $i$ at iteration $t$.
$\mathcal{H}\subseteq \mathcal{P}$	A subset of $\mathcal{P}$ containing four candidate locations that can be selected to generate a new tour.
${\mathcal{h}}_k$	One of the four candidate locations contained in $\mathcal{H}.$
$d_{(h_j,{\mathcal{h}}_k)}$	The distance between the $j$-th and last location of $p_i^j\left(t\right)$ and the candidate location ${\mathcal{h}}_k$.
$w_i\left(t\right)$	A four-elements vector $w_i\left(t\right)=〈w_1,w_2,w_3,w_4〉$ that determines the tendency of particle $i$ at iteration $t$ to change its position, moving toward one of reference positions.
$\mathcal{D}\left(t\right)$	A rescaled distance vector generated at each iteration $t$, which is used to create the new solutions $p_i\left(t\right) \forall i$.
$\overrightarrow{\mathcal{H}}$	The array of the candidate solution sorted in ascending order of their rescaled distance.
$g$	The greediness factor used to quantify the propensity of a particle to exploit the greedy solution.
$\alpha$	The shape factor of the quasi-geometric distribution used to select a candidate location ${\mathcal{h}}_k$ among the four included in $\overrightarrow{\mathcal{H}}$.
$\gamma$	The local search probability, i.e., the probability of performing a 2-Opt-based neighbour search on a starting solution.
$\beta$	The probability of exploiting the aisle-based generation mechanism; in a certain way, this probability codifies the propensity of a particle to follow a predetermined path inside an aisle of the warehouse.

. This notation integrates that of Table 1, valid for a generic PSO.

5.2. Basic Framework of the Algorithm

The iterative procedure of the Pkt_PSO follows a standard framework, like the one described in Section 3. Specifically, at every iteration $t\in \left\{1, 2, \dots , T\right\}$:

• Each particle $i$ moves from its previous position $p_i\left(t-1\right)$ to the new position $p_i\left(t\right)$. This part is the core of the optimization procedure, which is explained in full detail in Section 5.4 and Section 5.5.
• Next, with probability $\gamma$, a neighbour search is made on $p_i\left(t\right)$ using a 2-Opt-based approach, as explained in Section 5.6.
• If the solution has improved, the personal best is updated of particle $i$ is updated, i.e., $pbes{t}_i\left(t\right) \leftarrow p_i\left(t\right), \mathrm{if} f(p_i\left(t\right))<f\left(pbes{t}_i\left(t-1\right)\right)$.
• When all particles have moved to a new position, the minimum personal best is computed and compared to the global best. In case of improvement, the global best is updated accordingly, i.e., $gbest\left(t\right)\leftarrow \underset{i}{\min }pbes{t}_i\left(t\right), \mathrm{if} \underset{i}{\min }f\left(pbes{t}_i\left(t\right)\right)<f\left(gbest\left(t-1\right)\right)$.

Lastly, after $T$ iterations (or when a termination criterion is met), the algorithm ends and the best solution $gbest\left(T\right)$ found so far is returned as the output. Since optimality cannot be guarantee, this solution is referred to as pseudo-optimal solution.

5.3. New Solution Generation

Anytime particle $i$ moves, a new solution $p_i\left(t\right)$ is generated in an element-by-element way, by sequentially adding to an initially empty tour $p_i^0\left(t\right)$ one location at a time. Also, the locations used to form $p_i\left(t\right)$ are sequentially extracted from four reference solutions, that are: $pbes{t}_i\left(t\right), gbest\left(t\right), greedy$ and $rnd\_intentio{n}_i\left(t\right)$. As revealed by the notation, $greedy$ is a solution that depends neither on iteration $t$ nor on particle $i$. As better described in Section 5.5, this solution is shared by all the particles of the swarm, and it is generated with a greedy algorithm at the very beginning of Pkt_PSO. From this point onward, it is no longer changed. Conversely, $rnd\_intentio{n}_i\left(t\right)$ is a random solution that is regenerated at every iteration $t$ for each particle $i$.

Concerning how locations are extracted from the four reference solutions, the following procedure is used. Suppose that a total of $j$ locations have already been inserted in the partial tour $p_i^j\left(t\right)$, i.e., the last inserted location is $h_j$ and $p_i^j\left(t\right) =h_1,\dots ,h_j$. Now, we need to define, among the remaining $\left(M-j\right)$ locations of the picking list, which one should be inserted at position $\left(j+1\right)$. To this aim, a set of possible options $\mathcal{H}\subseteq \mathcal{P}$ is created by taking all locations that come immediately after $h_j$ in the four reference solutions. More precisely, let $k_{gb},k_{pb},$ $k_{gr}$ and $k_{rn}$ be the position occupied by location $h_j$ in $gbest\left(t\right)$, $pbes{t}_i\left(t\right)$, $greedy$ and $rnd\_intentio{n}_i$, respectively. So, we have that $\mathcal{H}=\left\{{\mathcal{h}}_{gb},{\mathcal{h}}_{pb},{\mathcal{h}}_{gr},{\mathcal{h}}_{rn}\right\}$, where:

• ${\mathcal{h}}_{gb}=gbes{t}_i\left[k_{gb}+1\right]$,
• ${\mathcal{h}}_{pb}=pbes{t}_i\left[k_{pb}+1\right]$,
• ${\mathcal{h}}_{gr}= greedy\left[k_{gr}+1\right]$,
• ${\mathcal{h}}_{rn}=rnd\_intentio{n}_i\left[k_{rn}+1\right]$.

Note that, if the partial tour $p_i^0\left(t\right)$ is empty, the k-indexes cannot be defined and so they are conventionally set to zero, i.e., $k_{gb}= k_{pb}=$ $k_{gr}=$ $k_{rn}$ = 0.

To better clarify this concept, the extraction procedure is exemplified in Figure 4, for a partial tour made of two locations $p_i^2\left(t\right)=〈h_1=E,h_2=A〉$. Since location $A$ is the last one inserted, to generate the set of the candidate locations $\mathcal{H}$, we need to consider all the locations that come after $A$ in the reference solutions. For instance, in $gbes{t}_i\left(t\right),$ location $A$ is at position four (i.e., $k_{gb}=4$) and so ${\mathcal{h}}_{\mathrm{gb}}=gbes{t}_i\left[k_{gb}+1\right]=gbes{t}_i\left[5\right]=C$. Considering the other reference solutions, we finally have that: $\mathcal{H}=\left\{{\mathcal{h}}_{\mathrm{gb}}=C,{\mathcal{h}}_{\mathrm{pb}}=D,{\mathcal{h}}_{\mathrm{gr}}=\mathrm{B},{\mathcal{h}}_{\mathrm{rn}}=F\right\}$. Therefore, the next and third location of $p_i^3\left(t\right)$ (indicated with three question marks in Figure 1) will be chosen among $\left\{C,D,B,F\right\}$.

Once the four candidate solutions of $\mathcal{H}$ have been identified, one of them must be selected and appended to the partial tour. As the objective is to minimize the overall travelling distance, a natural choice is to link the selection probability to the distance $d_{(h_j,{\mathcal{h}}_{\mathrm{k}})}$ between the last location $h_j$ of the partial tour $p_i^j\left(t\right)$ and each location ${\mathcal{h}}_k$ of $\mathcal{H}$. In other words, the greater is $d_{(h_j,{\mathcal{h}}_{\mathrm{k}})}$, the lower is the probability of selecting ${\mathcal{h}}_k$ and vice versa. Also, to give a different importance to the reference solutions, the distance $d_{(h_j,{\mathcal{h}}_{\mathrm{k}})}$ should be rescaled, taking the Hadamard (or elementwise) product with the weighting vector $w_i\left(t\right)=w_1,w_2,w_3,w_4$.

This is shown in Equation (7), where $\mathcal{D}$ is the rescaled distance vector:

$$\mathcal{D}\left(t\right)=〈d_{(h_j,{\mathcal{h}}_{gb})},d_{(h_j,{\mathcal{h}}_{pb})},d_{(h_j,{\mathcal{h}}_{gr})},d_{(h_j,{\mathcal{h}}_{rn})}〉\circ 〈w_1,w_2,w_3,w_4〉$$

(7)

In general terms, neither $d_{(h_j,{\mathcal{h}}_{gb})}$ nor $d_{(h_j,{\mathcal{h}}_{pb})}$ should be modified, whereas $d_{(h_j,{\mathcal{h}}_{gr})}$ and $d_{(h_j,{\mathcal{h}}_{rn})}$ should be increased. The first two distances correspond to $gbes{t}_i\left(t\right)$ and $pbes{t}_i\left(t\right)$, which are the most important reference solutions, whereas the latter two correspond to $greedy$ and $rnd\_intentio{n}_i\left(t\right)$ that are the less important ones. Therefore, a natural choice is to use the following weight vector $w_i\left(t\right)=〈1, 1, \left(1/g\right),1/\left(g-1\right)〉$, for each particle $i$ and for each iteration $t$, where $g\in \left[0, 1\right]$ is a greedy factor, that measures how much the greedy solution outweighs the random one. A neutral level is obtained for $g=0.5$ and rescales the distances by doubling the original ones. We anticipate that, under some peculiar conditions, the weight vector $w_i\left(t\right)$ could be modified for a limited number of particles, to trigger an antistagnation mechanism, as described in Section 5.5.

To implement this logic, first the four candidate solutions ${\mathcal{h}}_k$ are sorted in ascending order of their rescaled distance, to form the sorted array $\overrightarrow{\mathcal{H}}.$ Next, an element of $\overrightarrow{\mathcal{H}}$ is randomly extracted using the quasi-geometric distribution defined as in Equation (8) [62].

$$f\left(x\right)={\left(1-\alpha \right)}^x$$

(8)

where $\alpha$ is a shape parameter in the interval [0, 1), and $x$ is the position occupied by a candidate solution in the sorted array $\overrightarrow{\mathcal{H}}$. Therefore, the selection probability is maximum for the element at position one and minimal for the element at position four.

It is also important to stress that the use of a quasi-geometric distribution has a twofold purpose, as clearly explained by Juan et al. [63]. First, the adoption of a theoretical distribution provides implementation advantages concerning the use of an empirical distribution. Second, the quasi-geometric distribution has a single parameter, which greatly simplifies the fine-tuning step.

We conclude by noting that two issues could occur during the selection process. Indeed:

• If a candidate location ${\mathcal{h}}_k$ is already included in the partial tour $p_i^j\left(t\right),$ ${\mathcal{h}}_k$ is unusable and cannot be included in $\mathcal{H}$.
• If the location $h_j$ (the last one of the partial tours $p_i^j\left(t\right)$) matches the last location of a reference solution, its candidate location ${\mathcal{h}}_k$ cannot be defined, since position [$M$ + 1] does not exist.

When this happens, the cardinality of $\mathcal{H}$ is less than four; anyhow, provided there is at least one element in $\mathcal{H}$, the selection procedure remains unaltered. Should $\mathcal{H}$ be empty, the next location $h_{\left(j+1\right)}$ is randomly extracted among all locations that have not yet been included in $p_i^j\left(t\right)$.

5.4. Exploiting the Structure of the Problem: Aisle-Based Construction and Greedy Algorithm

To exploit the underlying structure of the picking problem, an additional solution-generation scheme is considered. As described in Section 4.2, the picking problem is characterized by a very constrained layout that makes the path connecting two locations virtually unique [64]. In other words, whenever the picker goes from location $x$ to $y$, he or she necessarily passes in front of other locations. This is typical when $x$ and $y$ are at opposite sides of the same aisle, so that moving from $x$ to $y$, all the intermediate locations are necessarily encountered. For instance, referring to Figure 2, if the picker moves from location V to location D, he or she also passes in front of locations {O, J, C, U, P, I}.

Let $S_{\left(x, y\right)}\subseteq \mathcal{P}$ be the set of all the locations found on the minimal path connecting $x$ and $y.$ To take advantage of this fact, when a new candidate location ${\mathcal{h}}_k$ is selected, all locations $l\in S_{\left(h_j,{\mathcal{h}}_k\right)}$ will also be included in the partial tour, with probability $\beta$. Therefore, the following two partial tours may be generated:

• $p_i^j\left(t\right)=〈h_1,\dots ,h_j,h_{\left(j+1\right)}={\mathcal{h}}_k〉$, with probability $\left(1-\beta \right)$.
• $p_i^{j+s}\left(t\right)=〈h_1,\dots ,h_j,h_{\left(j+1\right)}=l_1,\dots ,h_{\left(j+s\right)}=l_m,h_{\left(j+s+1\right)}={\mathcal{h}}_k〉$, with probability $\beta$

where $s$ is the cardinality of $S_{\left(h_j,{\mathcal{h}}_k\right)}$. Please note that the first one is the standard tour, and the second is the one generated using the aisle-based construction algorithm.

Due to the very constrained graph of routes, we also suggest generating the greedy solution (shared by all the particles of the swarm) using the well-known nearest neighbour algorithm. According to this algorithm, the picker starts the tour from the I/O point and repeatedly visits the nearest location until all locations in $\mathcal{P}$ have been visited. Indeed, due to the constraints induced by the layout on the feasible paths, this solution is known to be good, and it is thus advantageous to include it within the four reference solutions used to build a new tour. Also note that this solution never changes and remains the same from iteration $t$ = 1 to iteration $t$ = $T$. This choice may seem rather odd, but it can be justified as follows. At the single-item level, the nearest neighbour move is certainly optimal; it is therefore useful to ensure that this move is always available when a new tour is generated in an element-by-element way. This is exactly why the greedy solution must be included in the reference solutions.

For clarity, Figure 5 provides a Python 3.9©-based pseudocode, showing the generation of a new solution $p_i\left(t\right)$.

5.5. Local Search via 2-Opt

Any time a new solution $p_i\left(t\right)$ is generated, with probability $\gamma$, a local search is also performed. The local search is based on the 2-Opt, a simple local search algorithm originally introduced by Croes in 1958 [65], which is characterized by a special swapping mechanism.

Specifically, the proposed local search can be executed in a complete and full-depth mode, which differ in the number of 2-Opt exchange steps that they perform. If the complete 2-Opt is used, a new sequence is generated by inverting the order in which the locations between $h_j$ and $h_k$ in $p_i\left(t\right)$ are visited. This procedure is called complete because the inversion is made for each possible couple of nodes $\left(j,k\right)$. To make a simple numerical example, let us consider a four-item sequence $p_i\left(t\right)=〈A,B,C,D〉$. In this case, a complete 2-Opt search generates the following six tours: ${〈B,A,C,D〉}_{1,2}$, ${〈C, B,A,D〉}_{1,3}$, ${〈D, C, B,A〉}_{1,4}$, ${〈D, C, B,A〉}_{1,4}$, ${〈A, C, B,D〉}_{2,3}$ and ${〈A, C, B,D〉}_{2,3}$, where the subscripts refer to node $j$ and $k$, respectively. Instead, when the full-depth 2-Opt is used, the complete 2-Opt is recursively repeated on each improved solution that is found. The two approaches are performed with probability ${\gamma }_c$ and ${\gamma }_f$, respectively, with $\left({\gamma }_c+{\gamma }_f\right)=\gamma$. Also, to balance accuracy and computational time, we suggest using ${\gamma }_f=0.5{\gamma }^2$ so that ${\gamma }_c=\left(\gamma -0.5{\gamma }^2\right)\cong$ $\gamma$.

For clarity, Figure 6 shows the Python 3.9 © pseudocode of the local search via 2-Opt.

5.6. Antistagnation Mechanism

Let $i^{*}$ be a particle for which the full-depth 2-Opt, executed at iteration $t$, has failed (i.e., $pbes{t}_{i^{*}}\left(t\right)$ has not changed). To avoid stagnation, particle $i^{*}$ activates a procedure that is meant to escape possible local minima. The proposed procedure acts on the weight vector $w_{i^{*}}\left(t\right)$, as explained next.

-

Only for step $\left(t+1\right)$ and limited to particle $i^{*}$, the weight vector is changed to: $w_{i^{*}}\left(t+1\right)=〈1, \infty , \left(1/g\right),1/\left(g-1\right)〉$.

-

Additionally, if $pbes{t}_{i^{*}}\left(t\right)=gbest\left(t\right)$, the weight vector becomes: $w_{i^{*}}\left(t+1\right)=〈\infty , \infty , \left(1/g\right),1/\left(g-1\right)〉$.

-

At the end of the generating process:

▪

The current personal best is replaced by the solution just created, regardless of its objective function, i.e., $pbes{t}_{i^{*}}\left(t+1\right)= p_{i^{*}}\left(t\right)$.

▪

The weigh vector is restored to its initial values, i.e., $w_{i^{*}}\left(t+2\right)=〈1, 1, \left(1/g\right),1/\left(g-1\right)〉$.

In this way, all locations $\mathcal{h}\in pbes{t}_{i^{*}}\left(t+1\right)$ become the least likely to be selected and the particle tends to move away (or escape) from its current best. Metaphorically, limited to iteration $\left(t+1\right),$ particle $i^{*}$ temporarily forgets its personal best (and eventually the global best) and explores the solution space without making use of this past information.

5.7. Time Complexity of Pkt_PSO

For the sake of clarity, the whole algorithm is summarized in the flow diagram of Figure 7. From the diagram, it should be clear that in each epoch $t$, each particle i generates at least a new solution using the element-by-element generation scheme, possibly enhanced by the aisle-based generation approach. Next, if a neighbour search is performed, an additional number of solutions is generated. This extra number can be easily determined if the complete 2-Opt is used. In this case $C\left(M,n\right)$ solutions will be considered, where $C\left(M,n\right)$ is the binomial coefficient. Conversely, determining an exact number of different sequences generated by the full-depth 2-Opt before convergence is challenging, as this number will vary depending on the problem size and characteristics. Anyhow, since the complexity of the full-depth 2-Opt is generally considered to be $O\left(M^2\right)$, we can confidently assume that the number of explored solutions is << $M^3$, especially because, in relative terms, the number of iterations required by the 2-Opt to converge tends to decrease as the problem size increases. or a further discussion on this point see [66].

Because of this, considering the number of particles of the storm $N$ and the probability to perform a neighbour search $\gamma$, the total number of generated solutions during $T$ epochs can be estimated as in Equation (9).

$$N\_Solutions=T·\left(N·\left(1+\left(\gamma -0.5{\gamma }^2\right)·\left(\begin{array}{c}M\\ 2\end{array}\right)+0.5{\gamma }^2·\phi \left(M\right)\right)\right)$$

(9)

where $\phi \left(M\right)\ll M^3$ is the number of solutions generated by the full-depth 2-OPT. For instance, using problem instances as the ones discussed next in Section 7, we experimentally found that $\phi \left(20\right)\cong 1400$, $\phi \left(30\right)\cong 3750$ and $\phi \left(30\right)\cong 8400$.

5.8. Concluding Remarks on the Proposed Algorithm

In a standard PSO, particles constantly move according to a velocity profile that codifies the knowledge of the good solutions visited in the past. Indeed, any time a particle moves, a new solution $p_i\left(t+1\right)$ is generated based on its current position $p_i\left(t\right)$ and current velocity $v_i\left(t\right)$. The latter, in turn, depends both on a particle’s personal best and on the swarm’s global best. Metaphorically, each particle represents a bird that explores an area looking for food and changes its trajectory based on its intuition (cognitive behaviour) and/or to imitate the paths followed by the other birds of the storm (social behaviour).

Akin to the standard PSO, we also construct new solutions exploiting $pbes{t}_i\left(t\right)$ and $gbest\left(t\right)$; however, aside from this similarity, a careful reader could spot some evident divergences. Specifically:

-

Due to the discrete nature of the problem, an analytic formula (as in Equation (1)) is not implemented to recombine the reference solutions; an element-by-element approach is used instead.

-

Two additional solutions, namely $greedy$ and $rnd\_intentio{n}_i\left(t\right)$, take part in the generation process.

-

Conversely, the current solution does not explicitly appear in this process.

-

The concept of speed does not seem to be explicitly implemented in the framework.

However, especially concerning the last two points, at a conceptual level, these differences are only apparent. Our algorithm can be reinterpreted in the following way. Each particle always moves to start from its personal best, which coincides with its current position. When the particle moves to a new point, it stops there if the solution has improved; if not, it returns to the starting point. Furthermore, the direction and modulus of the movement depend on four reference solutions and the weight vector $w_i\left(t\right)$. In this sense, the reference solutions can be seen as four cardinal points toward which the particles are attracted. Moreover, the weight vector can be seen as a velocity profile that determines how much a particle deviates from its precedent trajectory, heading toward one of the above-mentioned cardinal points.

Also note that the velocity vector is dynamic, and it is also specific for each particle of the swarm. Most of the time, the velocity vector is shared by all particles, but if a full-depth search is made on particle $i$ at epoch $t$, and the full-depth search fails, for that particle, the velocity vector immediately changes, helping the particle to escape the local minimum, where it presumably was trapped.

To conclude, we note that in the case of discrete problems, the existing literature has two main strands of thought. A minority of authors propose readapting PSO with few marginal changes, while most of the scholars disrupt the original version, giving it a new shape. Given the description of the procedure, we claim that our approach stays in the middle. We reinterpreted the concept of speed, taking inspiration from the recent work by [23], trying at the same time to preserve most of the underlying features of the original PSO, such as the intention of the particles.

6. Parameter Tuning

The algorithm has four operating parameters: (i) the shape factor $\alpha$ of the quasi-geometric distribution, (ii) the probability $\beta$ of activating the aisle-based construction procedure, (iii) the probability $\gamma$ to carry out a local search via 2-Opt, and (iv) the greediness factor $g$, which defines how much the exploration is influenced by the greedy solution. The number of particles of the swarm is another operating parameter that should be fine tuned. However, we verified that the results do not appreciably change if the number of particles remains within the range of 30–40, which is usually used in similar works proposed in technical literature. Therefore, to assure a fair comparison with the other algorithms used as a benchmark, the number of particles was arbitrarily set to 30. Even the probability of triggering a full-depth search should be optimized. However, this event is so rare that performing a detailed optimization would be useless. Hence, also this parameter was arbitrarily fixed to ${\gamma }^2$, as mentioned in Section 5.

For proper fine-tuning, as in a two-level factorial experiment design, we defined reasonable low and high values for each parameter. Specifically:

-

For the probability $\alpha$, the values 0.7 and 0.2 were chosen, as in Juan et al. [63]. When $\alpha =0.7$, the selection of a potential location ${\mathcal{h}}_k$ is strongly influenced by its rescaled distance. For $\alpha$ approaching one, the selection tends to be deterministic and the candidate solution with the lowest distance is almost always returned. Conversely, when $\alpha =0.2$, the selection is almost random (with pure randomness that would be achieved with $\alpha \approx 0$).

-

For the probability $\gamma$, the values 0.2 and 0.05 were chosen. Both values are rather low because the execution of a 2-Opt local search should be an extreme condition with a very low occurrence rate. For values higher than 0.2, the time complexity (see Equation (8)) would grow tremendously, making its practical implementation impossible in a real industrial context.

-

For the greediness factor $g$, we are “pioneers”, as this is a completely new parameter. Lacking scientific references (or experimental evidence), two extreme values 0.8 and 0.1 were tested. When $g=0.8$ the greedy exploration prevails on the random behaviour, and vice versa.

-

Also, probability $\beta$ is a new parameter and so, in line with the previous point, we tested two opposite situations using an extreme value of 0.8 and a low value of 0.1.

All the resulting $2^4$ combinations were tested using five different picking lists with an average length of 30 locations (i.e., $M = 30$). A common layout characterized by three blocks divided by two cross aisles and 10 racks per block was also considered. A more accurate description of this layout is presented in Section 7.2, where the numerical assessment of the algorithm is discussed. For each combination and each picking list, the algorithm was executed five times for a total of 400 repetitions.

Table 4. Parameters’ tuning obtained results (in red the experiments where Pkt_PSO did not find the global optimum, in green the best parameters combination).

Parameters	Picking List 1		Picking List 2		Picking List 3		Picking List 4		Picking List 5
	Mean	Std. Dev.	Mean	Std. Dev.	Mean	Std. Dev.	Mean	Std. Dev.	Mean	Std. Dev.
$g=0.8, \alpha =0.70,$ $\beta =0.8, \gamma =0.20$	533	0	619	0	502	0	432	0	389	0
$g=0.8, \alpha =0.70,$ $\beta =0.8, \gamma =0.05$	533	0	619	0	502	0	433	0	389	0
$g=0.8, \alpha =0.70,$ $\beta =0.2, \gamma =0.20$	512	0	619	0	508	3	432	0	389	0
$g=0.8, \alpha =0.70,$ $\beta =0.2, \gamma =0.05$	512	0	619	0	502	0	432	0	389	0
$g=0.8, \alpha =0.20,$ $\beta =0.8, \gamma =0.20$	534	0	619	0	508	2	432	0	389	0
$g=0.8, \alpha =0.20,$ $\beta =0.8, \gamma =0.05$	534	0	619	0	508	3	432	0	389	0
$g=0.8, \alpha =0.20,$ $\beta =0.2, \gamma =0.20$	513	0	619	0	502	0	432	0	389	0
$g=0.8, \alpha =0.20,$ $\beta =0.2, \gamma =0.05$	513	0	619	0	502	0	432	0	389	0
$g=0.1, \alpha =0.70,$ $\beta =0.8, \gamma =0.20$	512	0	619	0	502	0	432	0	389	0
$g=0.1, \alpha =0.70,$ $\beta =0.8, \gamma =0.05$	512	0	619	0	503	0	432	0	389	0
$g=0.1, \alpha =0.70,$ $\beta =0.2, \gamma =0.20$	512	0	619	0	502	0	432	0	389	0
$g=0.1, \alpha =0.70,$ $\beta =0.2, \gamma =0.05$	512	0	619	0	502	0	432	0	389	0
$g=0.1, \alpha =0.20,$ $\beta =0.8, \gamma =0.20$	512	0	619	0	502	0	438	2	389	0
$g=0.1, \alpha =0.20,$ $\beta =0.8, \gamma =0.05$	512	0	619	0	503	0	432	0	389	0
$g=0.1, \alpha =0.20,$ $\beta =0.2, \gamma =0.20$	512	0	619	0	502	0	432	0	389	0
$g=0.1, \alpha =0.20,$ $\beta =0.2, \gamma =0.05$	520	3	619	0	502	0	432	0	389	0
Global Optimum	512	[-]	619	[-]	502	[-]	432	[-]	389	[-]

reports the results, expressed as the sample mean and standard deviation of the total travelling time of each parameter’s combination and picking list.

The values shown in

Table 5. Comparison of the results, with Avg. Dev expressed in meters; best results are highlighted in bold.

Layout	Pick. List Length	NeNe		Pkt_PSO		PSO_1 [22]		PSO_2 [23]		ACO [9]		2-Opt
		Avg. Dev.	N. Wins	Avg. Dev.	N Wins	Avg. Dev.	N. Wins	Avg. Dev.	N Wins	Avg. Dev.	N Wins	Avg. Dev.	N Wins
L1	20	75.2	0/10	0	10/10	537.1	0/10	0.3	8/10	8.0	3/10	34.9	0/10
	30	88.1	0/10	0	10/10	1095.2	0/10	0.3	8/10	19.1	0/10	59.7	0/10
	40	90.8	0/10	0	10/10	1711.3	0/10	0.6	6/10	27.2	0/10	67.8	0/10
L2	20	74.0	0/10	0	10/10	390.5	0/10	0.0	10/10	5.5	0/10	66.8	0/10
	30	108.0	0/10	0	10/10	797.2	0/10	0.0	10/10	19.9	0/10	76.6	0/10
	40	118.7	0/10	0	10/10	1160.9	0/10	0.2	9/10	28.9	0/10	113.9	0/10
L3	20	162.3	0/10	0	10/10	720.7	0/10	0.0	10/10	16.1	0/10	112.8	0/10
	30	129.4	0/10	0	10/10	1262.8	0/10	0.0	10/10	19.4	0/10	124.1	0/10
	40	200.2	0/10	0	10/10	2089.1	0/10	0.0	10/10	41.8	0/10	186.8	0/10

demonstrate that Pkt_PSO is very consistent and reliable. Indeed, for each parameter’s combination, the standard deviation is always very low, and the mean value does not change significantly as the parameters change. Results also show that the only parameter that has a significant impact on the results is the greediness factor, which should be kept at its low level (i.e., $g=0.1$).

For clarity, we also highlighted in red the experiments where Pkt_PSO failed, at least one time, to find the global optimum shown in the last row of the Table 5. The global optimum was computed with the famous Concorde TSP solver (using the Python version available in the PyTSP library) which, according to Mulder and Wunsch [67], is one of the fastest TSP solvers. For instance, the tuple (533, 0), relative to the first parameters’ combination and first picking list), means that the algorithm always converged to a value of 533 higher than the global optimum. Similarly, the tuple (508, 3), relative to the third parameters’ combination and third picking list, means that the algorithm has not always converged to the global optimum. The mean is higher than the global optimum and the standard deviation is not. Similarly, we highlighted in green the combination of parameters that never failed in finding the global optimum. This is the combination we opted for, also because, having a low value of the local search probability $\gamma =0.05$ it is also the combination requiring the lower computational time. From Equation (8) and using $\phi \left(30\right)\cong 3750$, the expected number of generated solutions for each epoch are approximately 950.

7. Case Study

7.1. The Benchmarks

To test Pkt_PSO, a set of competitive benchmarks was chosen. Since our algorithms are inspired by the PSO and since PSO often outperforms other metaheuristics when applied to the TSP (see for example [20]), the very recent and highly performant PSO by Zhou et al. [22] (PSO_1) and by Zhong et al. [23] (PSO_2) were considered. We also included in the benchmarks the ACO by De Santis et al. [9]. This is one of the most recent solutions specifically designed for the picking problem in manual warehouses, and it has been shown to perform better than other high-performance algorithms. Lastly, as simple benchmarks, we also considered the basic solution provided by the nearest neighbour constructive algorithm and by the full-depth 2-Opt, the same algorithms used by Pkt_PSO to generate the greedy solution and to perform the neighbour search. For a more in-depth description of the selected benchmarks, please see Appendix A.

7.2. Layouts and Picking Lists

Three different layouts that differ in terms of the disposition and number of blocks and rack were considered. For the sake of clarity, Figure 8 shows what we mean by blocks and racks.

The first layout (L1) is a very long warehouse with no cross aisles, made of a single block with 30 racks divided by 32 aisles. The second one (L2) has the same storage capacity as the first, but the racks are arranged differently. Specifically, it is made of three blocks divided by two cross aisles, and, in each block, there are 10 racks divided by 12 aisles. The third one (L3) is a high-capacity warehouse made of three blocks divided by two cross aisles and 30 racks per block, each one divided by 32 aisles. For all considered configurations:

-

Each aisle has 11 storage locations on both sides; the only exceptions are the first and the last aisle that, being flanked by the wall, have only 11 storage locations on a single side

-

The I/O position is in the bottom-right corner

-

The storage locations have a standard size of 2 × 2 m, while the aisles and cross aisles are, respectively, 2 and 4 m wide.

For each layout, we tested several problems of increasing complexity. To the best of our knowledge, and according to [60], given a routing problem, the number of locations to visit is a good proxy of its complexity. For this reason, on each layout, we tested 10 picking lists of length 20, 10 of length 30, and 10 of length 40. We also verified that with a few locations $M\le 19$ all the algorithms always converged to the global optimum computed with the Concorde TSP solver. On the other one hand, a number $M>40$ would be very unusual in a manual warehouse and, due to capacity constraints, certainly unrealistic with a single tour. Lastly, we note that for each generated picking list, we randomly selected the storage locations, to simulate a purely random allocation policy.

7.3. Results

All the algorithms were coded in Python 3.9© using the standard CPython interpreter (to increase speed) and run on a standard PC powered by Intel i7 with an 8-core processor 3.2 GHz and 16 Gb RAM. We also open-sourced the source code and the Jupyter-Notebook© used to run the tests at the following link: <GitHub-dmezzogori/op-pso>. To determine the reliability of the compared approaches, for each layout configuration and picking list, each algorithm was launched five times, for a total of $\left(3\times 10\times 3\right)\times 5=450$ tests performed for each algorithm. Clearly, since the nearest neighbour is a constructive and therefore deterministic algorithm, it was tested only one time for each layout configuration and picking list for a total of $\left(3\times 10\times 3\right)\times 1=90$ tests. For each algorithm and for each run, the total travelled distance of the generated picking tour (i.e., the objective function) and the computational time were recorded. We clarify that we did not calculate the global optimum of each scenario, as this would have taken too long (especially for the 40-locations configuration) and also because: (i) the capability to converge to the real global optimum was previously tested with values of $M\le 19$ (and even with $M=30$ for Pkt_PSO) and (ii) all the benchmarks had already been compared with exhaustive procedures and/or with optimization solvers in their original papers, to which we refer for more details. In view of this, comparisons will be made relatively to the overall best solution found by the five algorithms; from here on this solution will be referred to as the pseudo-optimal solution.

Obtained results are reported in Table 5 in terms of two performance indicators:

-

Number of wins (N. Wins column), the number of times an algorithm found the pseudo-optimal solution (i.e., its solution was better or equal to that provided by the other algorithms).

-

Average deviation (Avg. Dev. column) relative to the pseudo-optimal solution, computed as in Equation (10) for layout $l$, picking list length $p$ and algorithm $a$ $Avg\_De{v}_{\left(l,p,a\right)} = \frac{{{\displaystyle \sum }}_{i=1}^{10}(resul{t}_{\left(l,p,a,i\right)}- pseudo\_optima{l}_{\left(l,p\right)})}{10}$

where the index $i\in \left[1, 10\right]$ indicates the $i$-th instance of the 10 randomly generated picking lists of equal length $p$.

For the sake of clarity, and due to space constraints, the starting data from which the performance indicators were computed (i.e., the length of the picking tours generated by each algorithm) have been omitted from the main body of the article. However, they are included in Appendix B.

For clarity, in

Table 6. Comparison of computational times in seconds.

Layout	Picking List Length	Pkt_PSO	PSO_1 [22]	PSO_2 [23]	ACO [5]	2-Opt
		Avg.	Avg.	Avg.	Avg.	Avg.
L1	20	32.57	0.49	22.56	20.55	0.02
	30	139.44	0.99	44.66	39.02	0.09
	40	410.84	2.03	88.44	53.76	0.36
L2	20	28.47	0.50	24.70	21.15	0.02
	30	113.20	1.17	42.75	34.46	0.11
	40	371.87	2.07	77.05	51.68	0.33
L3	20	32.35	0.53	22.51	22.05	0.02
	30	136.15	1.24	44.52	35.24	0.11
	40	425.26	2.39	74.59	49.94	0.38

we highlighted in bold the best result for a given layout and picking list combination, that is, the best result of each row.

A complementary analysis is provided in Table 6 that shows the average computation times in seconds. Note that the nearest neighbour algorithm has not been included in the table, as its running time is negligible, given its simplicity and determinism.

7.4. Discussion

Results synthesized in Table 6 and

Table 7. 3-way ANOVA to compare performance of ACO and Pkt_PSO, *** denotes a very high significant factor.

	Sum_Sq (Type III)	df	F	p-Value
Intercept ***	5.15 × 107	1	238,180.2	9.78 × 10−141
ALG	7.84 × 101	1	0.3623	5.49 × 10−1
LT ***	1.97 × 106	2	4542.17	4.36 × 10−83
PL ***	1.01 × 106	2	2351.34	8.64 × 10−72
(LO:PL) ***	5.32 × 104	4	61.54	1.18 × 10−23
Residuals	1.73 × 104	80	NaN	NaN

clearly show that Pkt_PSO is more accurate and precise than the benchmarks. Indeed, as Table 6 reveals, Pkt_PSO always converges to the pseudo-optimal at every run and so the average deviation (of the solution provided by Pkt_PSO with respect to the pseudo-optimal) is always null. This fact incontrovertibly proves that Pkt_PSO is more stable and consistent than the benchmarks and that the added procedures (i.e., aisle-based generation, 2-Opt, and deep search), tailored to the task at hand, tackle the problem structure effectively. Moreover, a consistent convergence to a superior solution and a null average deviation is observed for all the investigated combinations of layout and picking lists. This behaviour is extremely important, especially from an operational viewpoint, as interested practitioners should not bother to fine tune the hyperparameters of the algorithm for the specific layout at hand.

The only algorithm that shows similar results, although it does not match them, is the PSO_2. This PSO-based algorithm was able to consistently track the optimum in five out of nine configurations. In all the remaining configurations, the deviation from the best solution (found by Pkt_PSO) remained very contained, although not null. Besides being less precise, this algorithm is also less accurate, as it consistently converged to the same (supposedly best) solution at each attempt, only for the third investigated layout (L3). Instead, in all other cases, it found different solutions at every run, as certified by a non-null value of the average deviation. As for the remaining methods, none of them can track the optimum, and the average deviation is always very high.

Concerning the computational times shown in Table 7, we can say that the best operational trade-off is achieved by the PSO_2. This algorithm, in fact, has a low computational time comparable to that of the ACO, which is the fastest, but it provides much better results in terms of total travelled distance. Pkt_PSO, instead, while showing the best results, needs considerably longer computational times, probably due to the application of the 2-Opt procedure with deep search. These computational times are operationally acceptable and could be easily improved by optimizing the code and/or using a low-level programming language (such as C). Nonetheless, we checked the quality of the solution provided by our algorithm if it stopped after an amount of time equal to that required by the ACO to reach convergence. We performed this comparison for all scenarios described above, with five repetitions each. As shown in Figure 9, which displays the trend of the best solution found by the two algorithms (for problems of increasing complexity), the solutions are practically indistinguishable, with a maximum difference (or error) lower than 1%.

These results were also confirmed by the three-way ANOVA of Table 7, where “layout type” (LT), “picking list length” (PL), and “algorithm type” (ALG) are used as blocking factors. In the ANOVA, we also included the interaction between factor LT and PL; the other ones were excluded because they were non-significant in a preliminary analysis.

As the ANOVA shows, ALG is the only factor that is not significant (with a very high p-value), and so we can conclude that the performances of the two algorithms is indistinguishable. We also note that all the other factors are highly significant, and this is a further demonstration of the correctness of the choices we made to define the complexity of the problem.

8. Conclusions and Future Works

In this paper, we focused on the picking problem in manual warehouses, a very relevant problem with many managerial implications. Considering that a plethora of metaheuristics have been proposed to solve this problem, the paper aims to check if, and how much, the solution can be improved. Furthermore, since the problem is very constrained, it could be also characterized by many points of local minimum with a similar value. If so, at least operationally, a good local minimum could be enough, and it would be unnecessary to use a very sophisticated approach to find the absolute minimum.

To shed light on these research questions, we developed a novel metaheuristic (referred to as Pkt_PSO) and we compared it with a set of state-of-the-art benchmarks. Outcomes of the experimental campaign showed that Pkt_PSO provides the best results, in terms of both accuracy and precision. However, as we initially supposed, the obtained improvements are marginal, and to reach the optimum (or the pseudo-optimum), it is necessary to use an approach that favours exploration over convergence and that requires more computational time. If the execution of Pkt_PSO is interrupted as soon as the benchmarks reach convergence, the obtained results are almost indistinguishable. This finding shows that for this peculiar type of problem, research has practically come to an end, and the use of a general-purpose metaheuristic can be considered enough for many industrial applications. Nonetheless, the use of more accurate approaches could be useful when additional features and/or constraints are considered, such as dedicated storage, multilevel shelves, one-way aisles, and transport capacity of the picker. The application of Pkt_PSO to these alternative scenarios will be a topic of future works.

An interesting and promising field of application for future research could be the development of a methodology, based on the PSO, to solve the three-dimensional case picking problem (3D-CPP) in an order-picking strategy. 3D-CPP combines two important issues in logistics: the problem of pallet loading and the routing of pickers in manual warehouses. The proposed method could find an extension in the combined optimization of pickers and 3D bin-packing problem.