Comparison of MOEAs in an Optimization-Decision Methodology for a Joint Order Batching and Picking System

: This paper investigates the performance of a two-stage multi-criteria decision-making procedure for order scheduling problems. These problems are represented by a novel nonlinear mixed integer program. Hybridizations of three Multi-Objective Evolutionary Algorithms (MOEAs) based on dominance relations are studied and compared to solve small, medium, and large instances of the joint order batching and picking problem in storage systems with multiple blocks of two and three dimensions. The performance of these methods is compared using a set of well-known metrics and running an extensive battery of simulations based on a methodology widely used in the literature. The main contributions of this paper are (1) the hybridization of MOEAs to deal efficiently with the combination of orders in one or several picking tours, scheduling them for each picker, and (2) a multi-criteria approach to scheduling multiple picking teams for each wave of orders. Based on the experimental results obtained, it can be stated that, in environments with a large number of different items and orders with high variability in volume, the proposed approach can significantly reduce operating costs while allowing the decision-maker to anticipate the positioning of orders in the dispatch area.


Introduction
The interest in the operation of distribution centers has grown with the increasing e-commerce deliveries of small-size packages in a shipping supply chain [1].The efficiency of the operational process in storing facilities is critical for the overall performance of a firm [2][3][4][5].The preparation of batches amounts to between 50% and 70% of the operational costs in distribution centers, being labor-intensive activities in manual systems and capitalintensive activities in automatized ones.
The preparation of orders requires picking up goods from certain storage areas in response to specific requests made by customers [6][7][8][9][10].These processes cover the most expensive tasks in most warehouses, namely picking-up the requests from their sites of storage and conveying them to the preparation/dispatch area [11][12][13].Those articles are classified and consolidated in packages (boxes, pallets, containers, etc.) for dispatch [14,15].
Single unit loads usually include several goods that must be marked and labeled.These loads are examined to verify that they fulfill the orders and the corresponding documents are prepared.They are, finally, dispatched by loading them on the transportation units [11].
The costs of these activities are proportional to the processing time of the orders.The strict fulfillment of due times increases the processing times and underutilizes the capacities of pickers [11].This indicates the existence of a trade-off between the time devoted to preparing the batches and the possibility of either anticipating or delaying the positioning of some goods because the times can be achieved by partitioning the pick-up tasks and assigning them to different pickers, but this increases the distances covered and the total work time of pikers.Thus, the goal of minimizing picking times may not necessarily be desirable in all circumstances [16].On the other hand, anticipating the positioning of goods before the corresponding due time may shorten the pick-up distances and times [17].
This indicates that no single criterion of optimality exists in this setting.It is more adequate to find a Pareto set from which the planner may select the best solution given the circumstances.If the main concern is the congestion in the dispatch area, the strict satisfaction of due times will be the right goal.On the contrary, if the planner aims to reduce the costs of the pick-up process, the goal will be to relax the strict compliance with deadlines.This paper seeks to solve the order batching and the order picking problems (both are NP-Hard problems [6,7]) in an integrated and multi-criteria approach.The criteria considered are the minimization of the operational cost of the process of picking up the goods, which is in correspondence with the length of the pick-up path, and the minimization of the total earliness picking time, under a zero-delay-of-delivery policy.In the first stage, a multiobjective optimization (MOO) yields the Pareto optimal solutions.For this purpose, three hybrid evolutionary algorithms based on a dominance relation (NSGA-II [18], SPEA2 [19] and PESA-II [20]) are considered.In the second stage, a procedure orders the alternatives and selects the best one.For this, the Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) [21] is used.The performance of this integrated approach is evaluated by comparing the solutions obtained with the benchmarks of Tsai et al. (2008) [2].
Our main contributions are the following: 1.
The introduction of an integrated multi-criteria treatment of Order Batching and Order Picking problems.The decision criteria are the minimization of both operational cost and earliness in picking times.

2.
A more realistic analysis of the Order Batching and Order Picking problems using a multi-objective optimization model of nonlinear mixed-integer programming, by considering a multi-level storage system with an explicit inclusion of the due times of the requests and a zero tardiness policy.

3.
The incorporation into the model of the actual scheduling procedures for different pick-up teams with capacity restrictions used in real-world warehouses.4.
The consideration of storage systems that combine multiple two and three-dimensional blocks.

5.
Finding that the fronts of efficient solutions that reached the MOO stage dominate the results presented in [2,11], in which a single-objective approach is applied on smaller instances and simpler layouts without zero-delay policies.
This joint order-batching and multiple-pickers warehouse-picking problem has barely been addressed in the literature [22].As described in the next section, only two papers analyze it, but not in a Multiple Criteria Decision-Making (MCDM) approach taking into account economic factors and without considering how to anticipate with no delay the availability of orders for the next stage of the distribution process.
The structure of this paper is as follows: Section 2 presents the problem and examines the relevant literature.Section 3 and 4 develop the model.Section 5 presents the solution method to be applied.Section 6 presents the computer experiments to solve the problem and Section 7 discusses the results and potential further work.

Problem Description and Literature Review
The problem of picking up orders at a minimal cost, both in terms of time and resources, can be subdivided into three planning problems.One is how storage positions must be allocated to the received articles.The second one is how orders must be batched-up in lots to facilitate the collection process.Finally, there is the question of how to schedule the pick-up sequence to deliver the articles to the dispatch area [23].This paper is centered on developing an integral treatment of the last two problems.A solution to them will enhance the efficiency of the processes on the storage floor.These activities generate a large proportion of the costs of the distribution center since they are intensive in equipment and labor [24][25][26].
The order picking process can be described as follows.It starts with the reception of the orders submitted by multiple customers.Each order details the amounts of different articles requested as well as the due time of availability at the dispatch area.The goods will be grabbed from the corresponding storage positions by a pick-up team composed of several pickers.Each of them picks up a batch of items that could belong to different orders [17,27].Masae et al. [28] present a systematic survey of the literature on the specification of the picking routes.
Small orders can be finished in a single tour, reducing the length of displacements.This feature indicates that larger orders may be split up into smaller ones to reduce the total time of the order batching procedure.Alternatively, a bottom-up procedure is to combine small orders into a single large order that can be gathered in a single tour [29].
In either case, an integrated process is required to address the problem of reducing simultaneously the total cost of selecting and picking up the batches of several orders, satisfying the constraints imposed by the deadlines on deliverance.This cost increases monotonically with the time length of the displacements across the storage floor and how strictly the deadlines are complied with.Delays in the fulfillment of orders would lead to breaches of contracts, with ensuing penalties and other costs.If, instead, orders are delivered to the dispatch area before the deadline, they may congest it, delaying urgent requests and inducing larger total costs of operation.
Formally, this problem combines two NP-Hard problems, the Order Batching Problem (OBP) and the Order Picking Problem (OPP).The former amounts to determining the optimal configuration of the batches, under constraints given by the capacity of pick-up teams and the delivery times.OPP, instead, amounts to identifying the optimal sequence of pick-ups, as to minimize the length of displacements and the corresponding processing times, under the constraint that each storage site must be visited only once.This combined problem will be called OBP-OPP in the rest of the paper.
De Koster et al. [6] review thoroughly the literature on both OBP and OPP.Chen and Wu [14] apply a clustering method to solve OBP to satisfy demand patterns instead of minimizing the length of the tours.Ho and Tseng [30] examine the heuristic procedures applied to solve OBP.Henn et al. [15] apply Taboo Search and Attribute-Based Hill Climbing for OBP.Henn and Schmid [31] add the Iterated Local Search heuristic to this toolbox.Lam et al. [32] state OBP as an integer programming problem that is solved using a heuristic based on Fuzzy Logic.
Van Gils et al. [33] present an overview of the published contributions on OPP.Several authors, like Petersen [34], de Koster et al. [8] and Theys et al. [35], analyze the application of heuristics to the solution of OPP.So, for instance, Henn et al. [36] use Ant Colony Optimization and Iterated Local Search, while Chen and Lin [37] apply an efficient twostage method.Lu et al. [38] solve a dynamical OPP using a routing algorithm.
Diefenbach et al. [39] examine exact solutions for U-shaped layouts, resulting from the application of Benders' combinatorial decomposition, which yields accurate solutions to small instances.For larger instances, this research group applies a sweep algorithm as a heuristic approach.
Tsai et al. [2] solve the OBP-OPP problem formulated for storage positions in two and three dimensions by using a multiple genetic algorithm.Some degree of earliness and tardiness is achieved by defining flexible delivery time windows.Miguel et al. [3] develop a hybrid algorithm to solve the two-dimensional instances presented in [2].Miguel et al. [4] improve on the results of [3] by slightly modifying the representation of the problem.Miguel et al. [11] present improvements on the results obtained by [2] for multilevel instances.
Concerning the criteria of choice, Molnar and Lipovszki [40] consider both the distance covered and the weight of the articles.Pan et al. [41] and Battini et al. [42] study a genetic algorithm with Pareto-optimization and niche technique.Congestion and waiting times are considered in [41], while [42] include preparation times, human energy expenditure, and fatigue during the picking process.Elbert et al. [43] analyze the relative efficiency of routing deviation policies.Ardjmand et al. [16] minimize makespan and total travel time in a pick-up system but do not consider due dates for individual orders nor the possibility of being earlier or later with respect to them.
The model proposed by Vázquez et al. [44] minimizes e-fulfillment costs under the requirements of preservation of perishable products.Another single-objective model, proposed by Kocaman et al. [45], focuses on finding the best layouts of a unit-load warehouse for single-command operations.Physical distance between pickers, as a COVID-19 mitigation strategy, is included as a criterion in [46], while an ergonomic evaluation of the effort of pickers is introduced in [47].
Concerning the concrete problem under study here, only two authors have presented similar specifications.Cano et al. (2020) [48] develop several models for OBP-OPP with multi-picker sequencing, solving it for the weighted sum of objectives, one of which is tardiness.Cals et al. (2021) [49] also study this joint problem, minimizing the number of orders with delays, using Deep Reinforcement Learning (i.e., Reinforcement Learning in Deep Neural Networks).A review of the most relevant literature as well as a complete taxonomy of this problem can be found in Pardo et al. (2023) [22].

The Model
The proposed bi-objective OBP-OPP is based on a nonlinear mixed-integer programming formulation that uses the following sets, parameters, and variables.

Sets and Parameters
-P = 1, . . ., p, . . ., P is the set of indexes of different articles required items in storage.P i represents the subset of different articles requested by customer i.P b represents the subset of articles grouped in batch b, and which can be requested by different customers.-I = 1, . . ., i, . . ., Ī is the set of the index of customers and orders in the current wave.Ī is the number of customers and orders, each customer places a single order with different articles.-B = 1, . . ., b, . . ., B is the set of the index of batches to be picked up.B i represents the subset of batches containing items of order i. -W = w 1 , . . ., w p , . . ., w P is the set of weights, each article has a unit weight.-D = {d 1 , . . . ,d i , . . . ,d Ī } is the set of customer order deadlines.-L = ℓ 0 , ℓ 1 , . . ., ℓ p , . . ., ℓ P is the set of storage positions of the requested articles.ℓ 0 represents the dispatch area.ℓ p is the coordinates of the storage position of article p is given by x p , y p , z p .-K = 1, . . ., k, . . ., K is the set of the index of pick-up teams.K the total number of available teams.-Cap is the total capacity of the pick-up teams.-S b = ⟨s 1 , . . ., s u , . . ., s Sb ⟩ is the sequence of positions to be visited to conform batch b. -Q = 1, . . ., q i,p , . . ., Q is the number of units of item p requested by customer i.Q i = ∑ p∈P i q i,p is the total demand of articles by customer i.Q p = ∑ i∈I q i,p denotes the total number of requested units of p.
G = (V, A) defines an undirected graph where V represents the vertices.Each vertex corresponds to a storage location of an article p ∈ P, along with two additional copies (0 and n + 1) of the vertex representing the dispatch area.Consider A to be the set of edges joining pairs of nodes within V. ς represents the operational cost for each unit of time.Each edge (h, l) ∈ A has an associated transit time, t hl , calculated as the ratio of the distance between two positions (h, l) and the speed of the picking equipment, v (i.e., t hl = D h,l v).The average time to pick a unit of any item is represented by t pick once the piking equipment is positioned at the corresponding storage position.

Decision Variables
The model includes the following binary variables: - x hlkb , which equals 1 if h is grabbed immediately before the article at the storage site l by the pick-up team k according to the sequencing of batch b, where h, l ∈ V, k ∈ K and b ∈ B. This means that x hlkb = 1 if k has to go through edge(h, l) to pick up the goods in batch b. -y hkb equals 1 if k picks up the item in storage position h for batch b, where h ∈ V, k ∈ K and b ∈ B.
The model also includes the following continuous variables:

Bi-Objective OBP-OPP Model
In Expression (1) we have two objectives: f 1 represents the total cost of the pick-up process (in monetary units) corresponding to the sum of all the batches, while f 2 represents the total earliness in picking times.The travel times in f 1 are obtained by analyzing the layout of the facility.Constraints (2) indicate that the weight of a batch cannot exceed the capacity of the corresponding picking team.P b are the positions of the items in batch b.Constraints (3) indicate that, for each batch b, items in any given storage position can be picked up only once b.Constraints (4) indicate that the tour of each pick-up team starts at the dispatch area and returns to it.Restrictions ( 5) and ( 6) represent a property of conservation of the flow of pick-up operations.If team k obtains item l for batch b, it must pick up item h or the other way around.Constraints (7) mean that the requests of article p must be satisfied, while Constraints (8) indicate that all the demands of customer i must be satisfied.Constraints (9) guarantee that the picking process of batch b ends at the sum of its starting time and the total time devoted to the task.Constraints (10) say that the finishing time of an order i is the longest completion time of the batches for this order.
Restrictions (11) tell us that the finishing time of order i, corresponding to the time at which all the units have been collected, must be less than the deadline of the order.Constraints ( 12)-( 14) are restrictions on the ranges of the variables.Concerning earliness or tardiness in the orders' due dates, the earliness in making up the order i can be formally expressed as , where c i is the finishing time of order i.The constraints ( 9) allow E i ≥ 0 ∀i ∈ I but force T i = 0 ∀i ∈ I.
The first objective function expresses in monetary units the cost of the collection of all batches.It is proportional to the pick-up time, thus depending on the distances covered and the amount of articles in the batches.For the definition of the total earliness in picking times in the second objective function, notice that it depends on the specification of the pick-up tours because of the capacity limit of each pick-up team k.If two batches must be picked up by the same team, one batch must wait until the other is finished.If the orders in a batch have similar deadlines, the variety of items within the batch will increase, forcing a longer picking time.If instead a single item can be put in the first batch, the efficiency increases as well as the earliness in picking times.

Distribution Center Layout
Figure 1 shows an example of a small distribution center's layout.On the left at the bottom, this blueprint shows the access to the dispatch area.A picking team departs the dispatch area according to a prearranged itinerary, traveling from one position to another, picking up articles before moving on to the next position until all the goods in a batch are collected.Afterwards, the team returns to the dispatch area.The areas in which all the OBP-OPP operations are carried out are enclosed by a dashed red curve.
This type of configuration agrees with the one originally presented by Roodbergen and de Koster [50] and Tsai et al. [2].To check the validity of our solution method for the OBP-OPP problem, the parameters used by them are maintained.More precisely, it is assumed that there are storage racks on the side walls and double racks in the middle of each block of the storage area.Multiple blocks are also considered. Figure 1 shows a representation of a layout with three blocks.Goods are picked out along four vertical and four horizontal shelves.Orders are released in waves, each wave being made up of multiple orders that must be picked within a specific time window (which includes due dates) before the next set (wave) of orders is released.
The length of the path from the site of item l to that of h, i.e., from ℓ l to ℓ h , where ℓ l = (x l , y l , z l ) and ℓ h = (x h , y h , z h ) is given by: where A l and A h denote the vertical shelves or aisles along which ℓ l and ℓ h can be reached, respectively.Let B l and B h be the blocks to which each of them belongs.

Methodology
A two-stage procedure is proposed.The first stage involves a MOO in which a MOEA yields an approximate Pareto front.The second stage consists of the application of an MCDM procedure that assumes an L 1 metric applies the TOPSIS method to automatically obtain an approximate classification of solutions and selects the best-ranked one according to its similarity ratio.The quality of the solutions is tested against those in benchmark instances.

Solution Methods at the MOO Stage
Population-based MOEAs are rather popular in the MOO scientific community, exhibiting a remarkable performance when solving hard optimization problems [51].These algorithms do not guarantee the determination of the exact Pareto-Optimal front, but they come very close to it.The Non-dominated Sorting Genetic Algorithm-II (NSGA-II) [18], the improved Strength Pareto Evolutionary Algorithm (SPEA2) [19] and the improved Pareto Envelope-based Selection Algorithm (PESA-II) [20] are recognized algorithms in the multi-objective literature that are used in this work.
The parameters used to describe the algorithms are as follows: N, population size; P c , crossover probability; P m , mutation probability; N E , external population size; N I , internal population size; T, maximum number of iterations and D, additional parametric information to enable operations and evaluation during evolution, for example, weights D w , orders D I , due dates D dd , different items per order D p , required quantities of each type of article D q , etc.Each algorithm returns a set A * of non-dominated solutions.

NSGA-II
The Non-dominated Sorting Genetic Algorithm [52] classifies the solutions in layers, based on a non-dominance criterion.Each layer is assigned a rank proportional to the fitness of its individuals.It selects individuals for reproduction according to a stochastic remainder proportionate selection procedure.Its improved version NSGA-II [18,53] uses an elitism-based non-dominated sorting method for ranking and sorting each individual, and uses a crowding distance approach in its selection operator for keeping the diversity among the Pareto optimal solutions (see Table 1).[18]).

SPEA2
The Strength Pareto Evolutionary Algorithm (SPEA) was proposed by Zitzler, Laumanns, and Thiele in 2001.Each individual is assigned a strength value proportional to the number of solutions that it dominates.The fitness assignment process of SPEA yields solutions that are closest to the actual Pareto front.Its improved version, SPEA2 [19], uses a fine-grained fitness assignment strategy, taking into account, for each individual, the number of individuals that dominate it as well as the number of those that it dominates.A nearest neighbor density estimation technique enhances the search procedure.A binary tournament yields the selected individuals (see Table 2).[54].It defines a hyper-grid division of the phenotype space to ensure the diversity of solutions according to a crowding measure.In its improved version PESA-II [20], the selection is region-based and the subject of the selection is now a hyperbox, not just an individual, to reduce the computational cost associated with obtaining the Pareto front (see Table 3).end-while 17: end-while 18: return: P * These three MOEAs apply selection techniques to ensure the uniform dispersion of solutions on their Pareto front.PESA-II uses hyperboxes, while NSGA-II uses a measure of nearness among individuals.SPEA2 achieves an adequate degree of dispersion using an indirect procedure that measures the degree of isolation of the individuals.These differences impose distinctions in the ensuing results.

The Decision Method at MCDM Stage
The TOPSIS method was introduced by Hwang and Yoon [21].It is based on the axiom of Zeleny [55]: the rational choice is to select an action closest to the ideal or farthest from the anti-ideal.TOPSIS chooses the alternatives closest to the ideal (I + ) and farthest from the anti-ideal or nadir (I − ).In undefined situations, it uses a notion of similarity with the ideal [56].
The method goes through six steps: 1.
A decision matrix is defined in which each element x ij = f j (a i ) corresponds to the evaluation of the alternative i, according to criterion j .

3.
The initial values are normalized according to a procedure in [57]: where v 1j , v 1j , . . ., v nj is the normalized vector for criterion j.

4.
The ideal and anti-ideal solutions are identified.The ideal solution I + is is the solution with the best possible values for each of the objective functions I + = I + 1 , . . ., I + j , . . ., I + m .The anti-ideal solution I − is is the solution with the worst possible values for each of the objective functions I − = I − 1 , . . ., I − j , . . ., I − m .

5.
The distances of each action to its ideal and anti-ideal are measured by the metric L p .
Here, v ij is the normalized value of action i for criterion j, and m is the number of criteria.

6.
The similarity ratios S(v) are computed, yielding an ordering of actions: In [58] it was shown that even if I + and I − are unknown, v + and v − can be obtained from the approximate Pareto fronts using the L 1 distance.Since this is the case of the problem under consideration, this approach is followed.

Measures of Performance
where n is the number of non-dominated solutions and c i = f 2 1i + f 2 2i being f 1i and f 2i , respectively, the values of the first and second objective functions for the i th non-dominated solution.A lower value of MID indicates a better solution.
• Spread of Non-dominance Solution (SNS): It is a measure of the diversity of the Pareto front solutions.It is given by: • Hypervolume (HV): It measures the size of the region dominated by the Pareto front (P) and is limited by a point of reference dominated by the front.It takes into account both the convergence towards the Pareto front and the distribution of solutions: where x i is a solution in P, n is the number of non-dominated solutions in P, and A(x i ) is the rectangular area confined between the points x i and a reference point.If the Pareto front P a is a better approximation of the real front than the Pareto front P b , it will follow that HV(P a ) > HV(P b ).
All measures of performance are multiplied by the factor λ = 1.0 × 10 −3 to present a clearer representation of numerical results.

Characterization of the MOEAs
MOEAs will be used to find solutions in the instances presented by Tsai et al. [2] for single and multi-level storage systems.A representation as permutations of integers usual in the treatment of combinatorial problems is used.
Each chromosome has two genomes.The first holds information about the articles in each batch, while the second contains the pick-up sequence for each batch.
To illustrate the representation, a greatly simplified example is presented in Table 4, which shows the roster of orders slated for scheduling in the wave.The encoding of the order list in the example wave is presented in Tables 5-7.The first four rows of Table 5 and the first three rows of Table 6 provide additional information that remains unchanged throughout evolution.The last row of each table represents genomes 1 and 2 that compose the chromosome, respectively.
The chromosome composed of the two genomes is shown in Table 7.Each element of the first genome represents the batch to which each item is assigned, while the second genome represents the index in the visit sequence of each batch.The dimension of genome 1 corresponds to the number of orders, and that of genome 2 corresponds to the number of different items in the order list.The advantage of this representation is that it yields higher levels of efficiency than Holland's original binary representation [61], thanks to the incorporation of specific knowledge about the problem.The downside is that, instead of general operators, it requires problem-specific ones [62].
The crossover operator (Table 8) implements a constructive hybridized method based on the k-closest neighbor heuristic, to improve the sequences with λ-exchanges.
The mutation operator (Table 9) applies the well-known insertion operator and local search with λ-exchanges to improve the sequences.
The satisfaction of the family of Constraints ( 2)-( 6) and ( 9) is guaranteed by the hybrid nature of the genetic operators.The representation makes the Constraints ( 7) and ( 8) trivially satisfied.To initiate a satisfactory level of diversity, the population is initially randomized.The process terminates according to a criterion of limitation of costs, which restricts the maximum number of iterations.
We apply the approach of tournament selection [63], which prescribes selecting the individuals with the highest scores between k ones randomly chosen from the current population.Repeated applications generate a new population.Generate random integer, between 1 and B 5: Insert the integer into a random position between 1 and Ī in the first genome of the parent and store in y 31 6: Identify batches with changes in genome 1 7: Apply a lambda swap procedure to enhance the local sequence within the modified batches of genome 2 and save it in y 32 9: Termination: Join the pairs of genomes 1 and 2 of the mutated child and save it in Y 3

Data Sets and Parameters Settings
The OBP-OPP instances generated by Tsai et al. [2] consider the positioning of the element in two and three dimensions, as can be seen in Table 10.We have previously used a single-objective hybrid evolutionary algorithm to obtain satisfactory solutions to some 2D and 3D instances in [3,4,11].These cases are addressed here from a multi-objective approach.
To study this more realistic but more complex setting, we evaluate the performance of the hybridizations of three MOEAs. 1 Small; M: Medium; L: Large.
Each customer i has a request with a due date uniformly distributed within the time range from 10:00 a.m. to 6:00 p.m., denoted as t i ∼ U(36,000, 64,800) s.The uniform distribution q i,p ∼ U(1, 10) describes the required quantities q of an item p by a customer i.The probability distribution of the number of distinct items per order follows a normal distribution, with a mean of 10 distinct items per order and a standard deviation of five distinct items per order, i.e., |P i | ∼ N(10, 5).Each item p has a unit weight evenly distributed between 8 and 24 kg, i.e., w p ∼ U (8,24).We also suppose an average speed of v = 2 m/s, an average time for pick-up of any item of t pick = 15 s, a cost of traveling per unit of time ς = 0.05 and an instance-dependant capacity (Cap).
A strategy used in the literature to limit the search space is to consider the upper and lower limits to the number of possible lots based on the level of capacity utilization.To facilitate a comparison with the results of Tsai et al. [2], these limits are assigned as follows: Bmin = φ 1 ∑ p∈P w p /Cap and Bmax = φ 2 ∑ p∈P w p /Cap, where φ 1 and φ 2 are constants such that φ 2 ≥ φ 1 .If φ 1 is too big or φ 2 is too small, it may generate non-feasible sequences.In that case, longer tours and heavier batches can overtake the capacity of the picking equipment.Analogously, a bigger φ 1 may generate numerous batches incurring significant travel expenses.
The first stage prescribes the assignment of values to other parameters using standard procedures of hyperparameter tuning.Performance evaluations of the algorithms were conducted in terms of the described metrics, followed by the selection of configurations that yielded the most favorable results.So, the maximal number of generations is T = 500, the population size is N = 40, the internal population size is N I = 40, the external population size is N E = 40, the tournament size in the process of selection is st = 2, the parameters in the number of batches bounds are φ 1 = 2 and φ 2 = 4, the probability of mutations is P m = 0.02, the probability of crossover is P c = 0.8, using direct sampling as an elite selection rule.Each algorithm is run 200 times, each run independent of the others, for each test instance selected from [2] under the same initial conditions, starting with a randomly generated population.
Subsequently, the second stage of the decision-making process is carried out and the best-ranked alternative is determined using the TOPSIS [21] methodology, assigning equal weights to the objective functions.

Numerical Experiments
First, the proposed methodology is applied to the joint OBP-OPP bi-objective problem in a multi-level storage system.The objectives of the problem are the minimization of the total operational cost of the pick-up process and the minimization of the total earliness in picking times.The runs were performed on a PC with a 3.00 GHz processor and 8 GB of RAM.The results presented in this section correspond to 200 independent runs for each instance.
Table 11 presents the mean values for each of the four performance measures on the Pareto fronts generated over the 200 independent runs on the seven sizes of the problem.Figures 2a-d present the corresponding box plots.Figure 2a shows that the median number of solutions in the Pareto fronts (NPS) is practically the same for all three algorithms.Figure 2b shows that for large instances the distance to the ideal solution (MID) is shorter for NSGA-II and SPEA2.The measure of the dispersion of non-dominated solutions on the Pareto front (SNS), shown in Figure 2c, yields, in general, better results with PESA-II than with NSGA-II and SPEA2.In terms of hypervolume (see Figure 2d), in large instances, NSGA-II shows a better performance than PESA-II and SPEA2.
Pairwise Tukey tests are performed to provide a statistical validation of the results.This test reduces the risk of incurring Type I errors.The hypotheses to be tested are as follows: where u i − u j is the difference between the values of a pair of MOEAs under a given metric.The significance level (α) is set at 0.05.Table 12 presents the results of the hypothesis tests.A p-value < α indicates that there is sufficient evidence to reject the null hypothesis; therefore, it can be concluded that significant differences are found between the means compared.Taking this into account, it can be seen that, for the NPS metric, in general, no significant differences are observed seen that the solutions of TOPSIS, in the smaller instances, do not differ much from those obtained by Tsai et al. (2008) [2] for the M1 model but, in larger instances, the solutions of TOPSIS are better in terms of TOC.

Conclusions
This paper presents an integrated procedure to solve the joint OBP-OPP problem by combining an MOO and an MCDM.The first one consists of a MOEA that produces an approximate Pareto front.In the second one, TOPSIS is applied, which ranks the solutions of this approximate Pareto front according to their similarity ratios and selects the best one.
The objective of the procedure is the minimization of both the total operational cost of the pick-up process and the total earliness picking time.This multi-criteria approach considers the trade-off between minimization of operational costs and inefficiency due to earliness in the positioning of items to be included in the batches, assuming that no delays are allowed.
Three multiobjective hybrid evolutionary algorithms based on the dominance relation among solutions (NSGA-II, SPEA2, and PESA-II) are applied and compared on small, medium-sized, and large instances of mono and multilevel storage systems by Tsai et al. (2008) [2] Four performance metrics are used to compare the performance of the proposed multi-objective algorithms according to different measures: NPS, MID, SNS, and HV.
The results were statistically validated using Tukey's test of difference of means for each of the metrics.It is found that on medium-sized and large instances, the hybridized NSGA-II algorithm achieves Pareto fronts with mean hypervolumes up to 10% above those of the other two algorithms, while the distances to the ideal solution are around 17% shorter, on average, for the large DS3 and DS6 instances.On the other hand, NSGA-II achieves a similar or even higher mean dispersion on large instances, although in smaller instances the other two MOEAs show better results.
In small-sized instances, TOPSIS does not yield better top-ranked results than those of the M1 model of Tsai et al. (2008) [2] .On the contrary, in large instances, the best-ranked results of the hybridized TOPSIS procedure proposed in this paper dominate those of the benchmark model.
Taking into account the growing importance of e-commerce and the high costs of picking activities in logistic centers, the proposed methodology can contribute to reducing those costs and gaining efficiency in their procedures.
Future work includes finding optimal placings for the articles, simultaneously taking into account aspects such as the size of the goods, the volume of their packaging, their relative demands, as well as operational considerations.Funding: "Consejería de Economía, Industria, Comercio y Conocimiento" of the Government of the Canary Islands through the direct grant awarded to the ULPGC called "Support for R+D+i activity.Campus of International Excellence CEI CANARIAS-ULPGC".

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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

-
A time variable g b defines the starting time of batch b. -A time variable f b represents the end of the picking process of batch b.

HFigure 1 .
Figure 1.Illustrations of a small layout representation with a position marker for a particular batch.(a) Orders batch (highlighted) in warehouse layout.(b) A graph-based representation of the underlying topology.(c) A picking tour for a particular batch of orders.

Four
metrics are applied to compare the performance of the multi-objective algorithms [59,60]: • Number of Pareto Solutions (NPS): This metric counts the number of non-dominated solutions found by each algorithm.• Mean Ideal Distance (MID): It measures the closeness between the Pareto solution and the ideal point (0, 0) as:

Figure 3 .
Figure 3. Non-dominated solutions medium and large problem size.

Figure 4
Figure 4 depicts these differences.It can be seen that the benchmark solutions are dominated by the best-ranked solution by TOPSIS on the front obtained by the hybridized NSGA-II algorithm.

Figure 4 .
Figure 4. Comparison of the TOPSIS-based results with M1 in all instances.

F
.T. and B.G.; Formal analysis, F.M.M., M.F., M.M., F.T. and B.G.; Investigation, F.M.M., M.F., M.M., F.T. and B.G.; Funding acquisition, F.M.M., M.F., M.M., F.T. and B.G.All the authors have contributed equally to the realization of this work.All authors have read and agreed to the published version of the manuscript.

Table 4 .
List of orders.

Table 8 .
Crossover Operation.1:Input:X 1 , X 2 and D. Parent solutions and data parametric information 2: Output: Y 1 and Y 2 .Offspring solutions generated by the crossover operation 3: Initialization: Obtain the genomes of each parent.x11,x12 , x 21 , and x 22 4:Generate k random integers, between 1 and Ī 5:Swap k elements belonging to the first genome in the two parents and store in y 11 and y 21 6:Identify batches with changes in genome 1 of each offspring 7:Apply the closest neighbor heuristic and store the modified batches in the second genomes of each offspring y 12 and y 22 9: Termination: Join the pairs of genomes 1 and 2 of each offspring and save them in Y 1 and Y 2

Table 9 .
Mutation Operation.1: Input: X 3 and D. Parent solution and data parametric information 2: Output: Y 3 .Offspring solution generated by the mutation operation 3: Initialization: Obtain the genomes of parent.x 31 , and x 32 4:

Table 10 .
Order batching and picking problem instances.

Table 11 .
Numerical results of mean performance measures.

Table 14 .
Comparisons of mean criteria ratios.