Two linearized finite difference schemes for time fractional nonlinear diffusion-wave equations with fourth order derivative

1.   Introduction

Interest on the operation of warehouses and distribution centers has grown considerably in recent years, particularly due to the growing importance of e-commerce and the increasing number of deliveries of small size packages ^[1]. Improving the efficiency of the operational process in storing facilities is critical for the overall performance of a firm ^[2,3,4,5]. These processes include all the activities carried out in distribution centers, like receiving, transferring, putting away, sorting, picking-up, classifying, grouping, accumulating and dispatching orders ^[6,7,8]. This requires coordinating different technical teams and designing the lay-out of the storage center ^[9,10].

The receiving activities consist in unloading orders from vehicles, carrying them inside the depot, updating the inventory databases and inspecting the items to detect eventual inconsistencies with the quality, quantity and packaging declared. The putting away tasks amount to moving the items from the unloading dock to their assigned place in the storage area, registering their corresponding placement. This involves handling the goods, verifying the suitability of the placement sites and proceeding to store them until they are requested in an order ^[11]. The order picking process covers the most expensive activities in most warehouses, consisting in picking up the right amount of requested goods from the corresponding storage sites and taking them to the preparation/dispatch area ^[3,12,13]. Classifying and grouping the requested and picked up articles consists in consolidating the orders of the clients in the modular packages (boxes, pallets, containers, etc.) to be dispatched ^[14,15]. In most cases this involves marking, labeling and grouping several goods into a single unit load. Dispatch is the process of verifying the load units to be transported, checking that the orders received are fulfilled, elaborating the documents to accompany them, to finally load them on transports.

Given the high costs in resources and time of the logistics in the storage facilities, leads us to focus on the order picking activities. We are particularly interested in the solution of the integral batching and pick-up problems, by taking the displacement costs and the earliness and tardiness penalties into consideration. We treat this combined problem in a multi-level storage system (with a 3D positioning of items) solved by applying a hybrid evolutionary algorithms.

The rest of this paper is organized as follows: Section 2 describes the problem and reviews the relevant literature. Section 3 presents the formal model, while Section 4 describes the algorithm devoted to solve this question. Section 5 presents datasets with different order characteristics and 3D warehouse environments and section 6 discusses the results of computer experiments Finally, section 7 discusses the results and the prospects for further work.

2.   Problem description and literature review

To pick up orders fastly and at a minimal cost, three planning problems, involving operations in warehouses or distribution centers, have to be addressed. One is the allocation of storage positions for the received articles. Another one is batching-up the orders in lots to be collected. Finally, we have the problem of scheduling the sequence of pick-ups and delivery displacements to the dispatch area ^[16]. In this paper we focus on the integrated treatment of the two last problems, which are critical for the efficiency of the activities in the storage floor. Furthermore, they generate most of the costs associated to the operation of the distribution center since they involve activities that are intensive in labor and equipment ^[17,18,19].

A more precise description of the order picking process is the following. The process starts with incoming multiple customer orders. Each order involves requests of different articles made by a customer, detailing the amount of each article and the due time at which the order has to be available at the dispatch area. The articles must be picked-up from the storage positions at the right time by a pick-up team that has to visit them in a sequence ^[20,21]. Masae et al. ^[22] review systematically the literature on order picker routing policies.

Small orders may reduce the displacement time, by finishing them in single picking tours. This can be taken advantage by order batching procedures that split up larger orders into smaller ones. Alternatively, small orders can be combined in a single large order that can then be picked in a single picking tour ^[23]. Figure 1 depicts this process, indicating that it starts and ends at the dispatch area.

So, the integrated process consists in designing a plan to reduce the cost of selecting, picking-up batches of different orders simultaneously, given precise deadlines for the finishing of each one. This cost is proportional to the displacement times across the storage floor and the punctuality in finishing the requested tasks. Since the final goal is to dispatch the batches, the schedule of pick-ups depends on the schedule of deliveries to customers. If the order were finished with delays it would lead to violations of the contracts with customers. Penalties and related costs ensue. On the other hand, if the orders are finished before the stipulated time the goods may clog the regular flow in the dispatch area, increasing the processing time of urgent requests, increasing the total cost of operations.

In more formal terms, this problem integrates two already known NP-Hard problems, the order batching problem (OBP) and the order picking problem (OPP). OBP consists in determining how to conform the batches that will be collected by the pick-up teams, taking into account the capacities of the teams and the time at which each article should be at the dispatch area. OPP in turn consists in identifying optimal visit sequences to the storage positions in which the goods of a given batch can be found, minimizing the distance covered and the processing time of the batch, visiting each storage position only ones. We denote the integrated problem as OBP/OPP.

De Koster et al. ^[7] present a thorough literature revision on OBP and OPP. Ho and Tseng ^[24] reviewed the heuristics used to solve OBP. Chen and Wu ^[14] use a clustering-based approach to solve OBP taking into account demand patterns instead of the distance covered in each sequence of visits. Henn et al. ^[15] proposed using heuristics like taboo search and attribute-based hill climbing for OBP. Henn and Schmid ^[25] added the iterated local search heuristic to address OBP. Later on, Lam et al. ^[26] proposed and integer programming model for OBP in which the distance covered by each sequence is estimated and the problem is solved using a fuzzy logic-based heuristic.

Van Gils et al. ^[27] reviewed and classified the literature on OPP. Petersen ^[28], De Koster et al. ^[6] as well as Theys et al. ^[29] presents revisions of the heuristics for OPP. Henn et al. ^[30] used ant colony optimization and iterated local search to solve OPP, while Chen and Lin ^[31] used a very efficient two-stage method. Lu et al. ^[32] present a routing algorithm for dynamical OPP.

Tsai et al. ^[2] used a multiple-GA to solve the integrated OBP/OPP on 2D and 3D item positions. The novelty of this work resided in the use of flexible time windows for the delivery of orders, allowing a certain degree of earliness and tardiness. Miguel et al. ^[3], presented an evolutionary algorithm hybridized with a local search method to solve the 2D instances in reference ^[2]. Later on, Miguel et al. ^[4] changed the representation of individuals, improving the results of reference ^[3].

In this paper we seek to further test and improve the performance, using experimentation, of the algorithm presented in reference ^[3] applying it on the 3D instances of reference ^[2].

The main contributions of this paper are the following:

ⅰ) We addressed the order batching and order picking with a more realistic model, by considering a multi-level storage system (with a 3D positioning of items), explicitly incorporating the due times of the requests.

ⅱ) We take into account that the actual procedures in warehouses involve several pick-up teams. This has been barely approached in the literature on the processing and routing of pick-ups in distribution centers.

ⅲ) The improvements are obtained by implementing an algorithm that solves the problem in a single stage, unlike the more complex algorithm presented in ^[2], which requires 2 stages.

ⅳ) The algorithm is tested on sets of orders generated probabilistically using the methodology presented in reference ^[2], showing its robustness.

The literature has not covered sufficiently the resolution of the integrated problem by dividing larger problems into smaller ones to be fulfilled by multiple pick-up teams. While reference ^[33] comes close to this, it does not allow the splitting of orders. This aspect of our work indicates its scientific relevance.

3.   Specification of the problem

We present here a specification of OBP/OPP based on a mathematical programming formulation ^[2]. In this setting we define the decision variables, the objective function as well as the constraints that involve delivery deadlines, number of pick-up teams and the operational restrictions on the use of the storage facility.

3.1.   Parameters

$\mathcal{P}$ denotes the set of different items in storage¹. We assume that there are $nArt$ items. Each unit of an item has a weight, being the set of those weights $\mathcal{W} = \left\{{\mathcal{w}}_{1}, \dots , {\mathcal{w}}_{p}, \dots , {\mathcal{w}}_{nArt}\right\}$. ${\mathcal{P}}_{i}$ is the subset of items requested by customer $i$. We assume that each customer places a single order or request of different articles, with different amounts of them. Thus, the number of customers $nC$, equals the number of requests $nReq$, being the set of customers/orders $\mathcal{I} = \left\{1, \dots , i, \dots , nC\right\}$. Each request has a specified deadline, with the class of those deadlines being $\mathcal{T} = \left\{{t}_{1}, \dots , {t}_{i}, \dots , {t}_{nPed}\right\}$. ${\mathcal{P}}_{r}$ is the set of items in a batch $r$. The items in ${\mathcal{P}}_{r}$ can be requested by a single or several customers. $\mathcal{L} = \left\{{\mathcal{l}}_{0}, {\mathcal{l}}_{1}, \dots , {\mathcal{l}}_{p}, \dots , {\mathcal{l}}_{nArt}\right\}$ represents the set of storage positions of the types of articles plus ${\mathcal{l}}_{0}$ which indicates the dispatch area. For an item $p\in \mathcal{P}$, its storage position ${\mathcal{l}}_{p}$ is given by its coordinates in the store, ${\mathcal{l}}_{p} = \left({x}_{p}, {y}_{p}, {z}_{p}\right)$. $\mathcal{R} = \left\{1, .., r, \dots , \left|\mathcal{R}\right|\right\}$ indicates that set of batches to be picked-up. ${\mathcal{S}}_{r} = \left\langle{{s}_{1}, \dots , {s}_{u}, \dots , {s}_{\left|{\mathcal{S}}_{r}\right|}}\right\rangle$ is the sequence of positions to be visited to conform batch $r$, where ${s}_{u}$ is the $u-th$ position to be visited. ${q}_{i, p}\in Q$ denotes the number of units of item $p$ requested by customer $i,$ with ${\mathcal{Q}}_{i} = \sum _{p\in {\mathcal{P}}_{i}}{q}_{i, p}$ being the total number of goods demanded by $i$ while ${\mathcal{Q}}_{p} = \sum _{i\in \mathcal{I}}{q}_{i, p}$ denotes the total number of requested units of $p$. Finally, $\mathcal{K} = \left\{1, \dots , \left|\mathcal{K}\right|\right\}$ is the set of pick-up teams, each one with total weight capacity $Cap$.

¹ Each item has a different reference code (SKU). So for instance, 200 units of the smartphone EjS22 share the same SKU, which is not the same of the smartphone EjS7, another item.

We define an undirected graph $\mathcal{G} = \left(\mathcal{V}, \mathcal{A}\right)$, where $\mathcal{V}$ are the nodes, each one corresponding to a storage position of an item $p\in \mathcal{P}$, plus two copies ($0$ and $n+1$) of the node corresponding to the dispatch area. $\mathcal{A},$ represents the class of edges connecting pairs of nodes of $\mathcal{V}$. Each $edge\left(h, l\right)\in \mathcal{A}$ has an associated travel time ${t}_{hl}$ given by the distance between positions $h$ and $l$ over the speed of a picking-up team, $v$ (i.e., ${t}_{hl} = {D}_{h, l}/v$) and an operational cost for each unit of time, $\varsigma$. ${t}_{pick}$ is the average time to pick a unit of any item, once the team has reached the corresponding storage site.

3.2.   Binary flow variables

${x}_{hlkr} = 1$ if $h$ is picked-up right before item in storage position $l$ by the pick-up team $k$ in the sequence corresponding to batch $r$, where $h, l\in \mathcal{V}$, $k\in \mathcal{K}$ and $r\in \mathcal{R}$. This means that ${x}_{hlkr} = 1$ if team $k$ has to travel through $edge(h, l)$ to form batch $r.$

3.3.   Binary index variables

${y}_{hkr} = 1$ if $k$ picks up item in storage position $h$ for batch $r$, where $h\in \mathcal{V}$, $k\in \mathcal{K}$ and $r\in \mathcal{R}$.

3.4.   OBP/OPP model

s.t.

The objective function Eq (1) represents the total operational cost of the pick-up process, expressed in monetary units, corresponding to the collection of batches, plus penalties for earliness and tardiness. The travel times in the first term obtain from the analysis of the lay-out of the storage facility.

With respect to the second term, $\alpha$ is the penalty per unit of time for earliness, while $\beta$ is the current penalty per unit of time for tardiness. ${E}_{r}$ is the earliness in making up order $i$, while ${T}_{i}$ is tardiness in making up order $i$. Formally, ${E}_{i} = max\left\{0, {\mathcal{t}}_{i}-{c}_{i}\right\}$ and ${T}_{i} = max\left\{0, {c}_{i}-{\mathcal{t}}_{i}\right\}$, where ${c}_{i}$ is the finishing time of order $i$, defined as the time at which all the units of all items corresponding to order $i$ are picked up and collected to form the order.

Constraints Eq (2) forbid that the total weight of a batch exceed the capacity of the teams and their equipment, ${\mathcal{V}}_{r}$ is the set of storage positions of the items that belong to lot $r$ ($p\in {\mathcal{P}}_{r}$). Constraints Eq (3) indicate that each storage position cannot be visited more than once for each batch $r$. Constraints Eq (4) ensure that all pick-up teams start and end at the dispatch area. In turn Eqs (5) and (6) preserve the flow of pick-up operations. If pick-up team $k$ gets item $l$ for batch $r$, it has to have picked up item $h$ or vice versa. Constraints Eq (7) indicate that all the requests of item $p$ are satisfied, while constraints Eq (8) state that all the requests of customer $i$ have to be satisfied. Constraints Eqs (9) and (10) restrict the range of values of variables.

3.5.   Lay-out of the distribution center

Figure 2 shows the lay-out of the distribution center used in our previous analysis of the OBP/OPP problem. The bottom left corner corresponds to the access to the dispatch area, where all the sequences of pick-up operations start and end. That is, a team leaves the dispatch area following a pre-specified sequence, going from a storage position to another, picking-up the items before moving on to the next storage position, until all the articles in a batch have been collected. After that, the team returns with the articles to the dispatch area. Graphically, we have surrounded with a dashed red curve the storage and dispatch areas, where all the actual OBP/OPP operations are carried out.

This configuration agrees with the one presented originally by Tsai et al. ^[2]. We keep the parameters used by them to test and validate our proposed solution method for the OBP/OPP problem. More precisely, we assume that there are two lateral racks and two double racks in the middle of the storage area. Pick-ups can be carried out along two transversal and three longitudinal aisles. The distance traveled from the storage position of item $l$ to that of $h$, i.e., from ${\mathcal{l}}_{l}$ to ${\mathcal{l}}_{h}$, where ${\mathcal{l}}_{l} = ({x}_{l}, {y}_{l}, {z}_{l})$ and ${\mathcal{l}}_{h} = ({x}_{h}, {y}_{h}, {z}_{h})$, can be expressed as follows:

In this expression ${\mathcal{A}}_{l}$ and ${\mathcal{A}}_{h}$ denote the aisles along which ${\mathcal{l}}_{l}$ and ${\mathcal{l}}_{h}$ can be reached, respectively. $\mathcal{H}$ is the second coordinate of ${\mathcal{l}}_{0}$, the start and end of a pick-up sequence. $\mathcal{Y}$ is the length of the picking aisles.

4.   Hybrid evolutionary algorithm (HEA)

The algorithm applied was introduced in reference ^[3] to solve the 2D instances presented by Tsai et al. ^[2]. It consists of an evolutionary algorithm ^[34] that evolves the solutions to the OBP/OPP problem hybridized with a constructive method based on the $k$ closest neighbors heuristic, to improve the sequences using a local search with $\lambda$-exchanges. It uses the representation of as permutation of integers, usual in the treatment of combinatorial problems. Each chromosome consists of two genomes. A simple example shows how this representation works. Consider three requests and four items. Table 1 indicates that, for instance, request 1 consists of 2 items (one unit of A and two units of C), totaling 3 articles.

Table 2 shows a representation consisting of two rows, one for items and the other for requests. The columns indicate the total number of requested articles (in this case, 12).

Table 3 indicates how the information in Table 2 is used to define two genomes. The first genome contains the information of the items assigned to batches (we consider three batches in this case). So, for instance, one unit of item A in request 1, is assigned to batch 2. On the other hand, genome 2 corresponds to the sequence of visits to the storage positions of the items. For instance, in batch 2, item A is collected from position 2. Table 4 presents the final chromosome to be evolved.

Table 5 shows how the chromosome of Table 4 is decoded, indicating the batches and their corresponding pick up sequences.

We can see that this chromosome is not yet a solution for the problem since, for instance for batch 1, this sequence prescribes going from the dispatch area to the position where item C is, pick up one unit, then go the position of item D, pick up one unit and return to the position of C to pick up another unit. After that the sequence prescribes going back to the position of D and then again back to C. This is clearly inefficient, since it could pick up three units of C and then two of D (or the other way around) with a lower cost. The evolutionary process ends up discarding chromosomes like that of Table 4.

The main advantage of this type of representation that incorporates specific knowledge about the problem is that it allows reaching higher levels of efficiency than Holland's original binary representation ^[35]. The downside is that it requires using operators adapted to this representation instead of general ones ^[36].

The penalty term in the objective function facilitates discarding non-feasible chromosomes. On the other hand, the satisfaction of the family of constraints Eqs (2)–(6) is ensured by the hybridization of the genetic operators with the closest neighbor heuristic. The constraints given by Eqs (7) and (8) are trivially satisfied by the representation.

The initial population is generated randomly, to ensure a larger variety of alternatives. The process terminates according a cost-oriented criterion, which limits the maximal number of iterations.

We use here the tournament selection approach ^[37], according to which the individuals that get selected are those with the higher scores among k individuals chosen at random from the current population. This is repeated until a new population is generated.

Table 6 presents the pseudo-code of the algorithm, as presented in reference ^[3].

5.   Datasets and parameter settings

Tsai et al. ^[2] generated instances for OBP/OPP positioning items in two and three dimensions, as shown in Table 7. We have previously applied the HEA to find very good solutions to some 2D instances (DS0, DS1, DS2 and DS3) ^[3,4]. Here we address 3D instances (DS4, DS5 and DS6). In this, more complex setting, we can test and improve the performance of the algorithm.

We assume that the number of units of item $p$ requested by a customer $i$ is uniformly distributed between 1 and 10, i.e., ${q}_{i, p}\sim U(1, \dots , 10)$. In turn, the number of different ítems requested by a customer follows a normal distribution, with mean 10 and standard deviation 5, i.e., $\left|{\mathcal{P}}_{i}\right|\sim N\left(\mathrm{10, 5}\right)$. The deadline for finishing the request of $i$ follows a uniform distribution over the range of seconds between 10:00 am and 18:00 pm, ${\mathcal{t}}_{i}\sim U(36000, \dots , 64800)$. The unit weight of each item $p$ is uniformly distributed between 8 and 24 kilograms, i.e., ${\mathcal{w}}_{p}\sim U(8, \dots , 24)$. With respect to the parameters of the pick-up teams, we assume an average travel speed of $v = 2\;m/s$, an average pick-up time for each unit of ${t}_{pick} = 15$ seconds, a travel cost per unit of time $\varsigma = {\$}\;0.05$ and a capacity ($Cap$) that depends on the instance under consideration.

For the objective function we consider an earliness penalty per unit of time of $\alpha = 0.5$ while the corresponding penalty for tardiness is $\beta = 1$.

With respect to the strategy to limit the search space and facilitating the comparison with the results of Tsai, Liou and Huang ^[2], we postulate lower and upper bounds on the number of possible batches, ${\left|\mathcal{R}\right|}_{min}$ and ${\left|\mathcal{R}\right|}_{max}$, respectively. These bounds are defined as follows: ${\left|\mathcal{R}\right|}_{min} = ({\varphi }_{1} \cdot \sum _{p\in \mathcal{P}}{\mathcal{w}}_{p})/Cap$ and ${\left|\mathcal{R}\right|}_{max} = ({\varphi }_{2} \cdot \sum _{p\in \mathcal{P}}{\mathcal{w}}_{p})/Cap$, where ${\varphi }_{1}$ and ${\varphi }_{2}$ are constants such that ${\varphi }_{2}\ge {\varphi }_{1}$. If ${\varphi }_{1}$ is too large or ${\varphi }_{2}$ is too small, unfeasible sequences may be generated. The latter case generates longer pick-up routes and the batches may exceed the capacity of the pick-up teams. Similarly, a large ${\varphi }_{1}$ would induce a very large number of batches that may generate large travel costs.

Tables 8–10 summarize the information about instances DS4, DS5 and DS6, respectively. Each row indicates the items and amounts of them requested by customer $i$. The last column shows the deadline for finishing the requests.

A first stage of calibration allows assigning values to other parameters. So, the maximal number of iterations is $MaxGen = 500$, the size of populations is $PopSize = 150$, the size of tournaments in the selection process is $SizeTournament = 2$, the parameters in the bounds on the number of batches are ${\varphi }_{1} = 2$ and ${\varphi }_{2} = 4$, the probability of crossover is ${Pr}_{cross} = 0.9$, the probability of mutations is ${Pr}_{mut} = 0.15$ while the number of individuals in the elite is (5% of the population) $nElite = 0.05\cdot PopSize$, using direct sampling as elite selection rule.

6.   Computational experiment

In this section we present the results of running our HEA on 3D instances (DS4, DS5 and DS6) generated by Tsai et al. ^[2]. We used a PC with an Intel Core i7 3.00 GHz processor with a 8 GB RAM. 100 runs for each instance yielded the results presented in this section. Figures 3–5 depict the $TOC$ when the number of batches $\mathcal{R}$, varies between ${\left|\mathcal{R}\right|}_{min}$ and ${\left|\mathcal{R}\right|}_{max}$ for each instance.

We can see that the number of batches is optimal when $TOC$ reaches its lowest value. In order to find more precise optimal values, we penalized the solutions with excess weight, adding an additional cost to them.

Table 11 presents the best solutions obtained by the HEA for instances DS4, DS5 and DS6. We compare them to results in reference ^[2]. The rows indicate some natural measures of performance, like the total travel distance in the pick-up plan $\left({D}_{total}\right)$, the optimal number of batches, the mean of the distribution of distances travelled on batches $({D}_{med}$), the standard deviation of that distribution (${D}_{std})$, the sum of earliness and tardiness $\left({ET}_{time}\right)$ and the total operational cost $\left(TOC\right).$ Here distances are measured in meters, times in seconds and monetary costs in US dollars.

Table 12 indicates that HEA yields a better TOC on the different instances. This is, in general not the case for the total distance traveled (except for instance DS5).

Figure 6 presents boxplots indicating the variability, degree of asymmetry and the extreme values of the distributions of ${D}_{total}$, ${ET}_{time}$ and $TOC$ obtained with our HEA on instances DS4, DS5 and DS6, using the optimal number of batches in each case.

We can see that in the boxplot of instance DS4 for ${D}_{total}$, we have an atypical value close to the minimal one for 8 batches, but 50% of the values of the distribution of ${D}_{total}$ is way above that value. In the case of ${ET}_{time}$ the central values are below their corresponding benchmark values. This is also the case for $TOC$. In summary, even if HEA does not yield better results in terms of ${D}_{total}$, the quality of solutions measured by $TOC$ and ${ET}_{time}$ improves over the known results.

7.   Conclusions

In this paper we addressed an integrated problem that combines the Order Batching Problem (OBP) and the Order Picking Problem (OPP) using a hybrid evolutionary algorithm (HEA). This algorithm uses a novel representation of the problem using specific knowledge that allows restricting the search space. The crossover and mutation operators are hybridized with a heuristic that allows finding better routes for the pick-up teams. Tsai et al. ^[2] generated 2D instances of this problem (in which items are stored at the same level) and 3D ones (with items placed at different levels). We have previously applied the HEA to 2D instances, with very good results ^[3,4]. In this paper we have addressed the more complex 3D instances. The results, in terms of the objective function, improve over those in reference ^[2], which were obtained using a Multiple-GA. This completes the testing of the HEA.

Future work involves considering this problem but in a multi-objective and cross-dock setting. The reach of the problem can be extended, studying the relation between the placement of items in the different layouts and the distribution of the last mile. Other aspects should also be considered, in particular those impacting on the operational planning of the activities in storage sites, like the volume and size of the packages, the relative demands of the different items, their rotation and similar considerations.

Acknowledgments

The authors are grateful for partial support from the following sources: National Scientific and Technical Research Council CONICET Grant PIP No. 11220150100777, Universidad Nacional del Sur Grant PGI 24/J086 and Universidad Nacional de Río Negro PI-JI 40-A-917.

Conflict of interest

The authors declare they have no financial interests. In addition, the authors have no conflicts of interest to declare that are relevant to the content of this article.