Task Scheduling Problem of Double-Deep Multi-Tier Shuttle Warehousing Systems

Double-deep multi-tier shuttle warehousing systems (DMSWS) have been increasingly applied for store-and-retrieval stock-keeping unit tasks, with the advantage of a reduced number of aisles and improved space utilization. Scheduling different devices for retrieval tasks to increase system efficiency is an important concern. In this paper, a Pareto optimization model of task operations based on the cycle time and carbon emissions is presented. The impact of the rearrangement operation is considered in this model. The cycle time model is converted into a flow-shop scheduling model with parallel machines by analyzing the retrieval operation process. Moreover, the carbon emissions of the shuttle in the waiting process, the carbon emissions of the lift during the free process, and the carbon emissions of the retrieval operation are considered in the carbon emissions model, which can help us to evaluate the carbon emissions of the equipment more comprehensively during the entire retrieval task process. The elitist non-dominated sorting genetic algorithm II (NSGA-II) is adopted to solve the non-linear multi-objective optimization function. Finally, a real case is adopted to illustrate the findings of this study. The results show that this method can reduce carbon emissions and improve system efficiency. In addition, it also help managers to reduce operational costs and improve the utilization of shuttles.


Introduction
With the rapid development of automation technology and the related economy, the warehousing industry has gradually become a key link in economic activities [1]. Warehousing systems are mainly used to store and retrieve materials and commodities. The double-deep multi-tier shuttle warehousing system (DMSWS) is a new type of intelligent warehousing system, introduced in recent years for store-and-retrieve stock keeping units (SKUs). The basic components of a DMSWS are the double-deep storage rack (SR), buffer position, shuttles, lifts and the I/O point. In a traditional automated storage and retrieval system (AS/RS), each task is executed by an aisle-captive stacker crane capable of horizontal and vertical movement [2,3]. In contrast, the aisle-captive stacker crane is superseded by the cooperation of shuttles and lifts in the DMSWS. The shuttle can move horizontally on each tier for transportation of the required SKU. The vertical movement of the required SKU is controlled by the lift, which is located outside the buffer position of the doubledeep SR. The horizontal movement of the shuttle and the vertical movement of lift are independent of each other. The DMSWS not only has the advantages of high operation efficiency and high space utilization, but also has a high degree of flexibility, making it capable of meeting the fluctuation of orders. Moreover, in terms of economic adaptability and warehousing volume, the DMSWS is one of the main development trends in automated warehousing systems.
In a DMSWS, the double-deep SR has two storage position depths for storing SKUs. In most studies, the double-deep SR increases storage positions by 40-50%, compared with a single-deep SR [4,5]. The double-deep SR requires fewer aisles, resulting in higher space utilization. Many producers of warehousing equipment and logistics companies have begun to apply the DMSWS in practice, such as Knapp, Dematic, Alibaba, and JD. Although these DMSWSs are present in logistics companies for order distribution, to the best of our knowledge, extensive research has not been conducted on the DMSWS according to our investigation. There is a clear necessity to study DMSWS, especially in terms of task scheduling, which significantly affects the system efficiency. This observation motivated us to examine the task scheduling problem of double-deep multi-tier shuttle warehousing systems.
In the AS/RS, when the crane receives some tasks (inbound and outbound), the stacker crane in the aisle reciprocates between the starting position and destination point. The task scheduling problem involves how to complete all warehousing operations in the shortest time; this is highly similar to the travelling salesman problem (TSP). According to our review of the relevant literature, many studies on the task scheduling problem have focused on the TSP and its improved algorithms [6,7]. However, the task scheduling problem of DMSWS is different from TSP, as storage and retrieval tasks can be executed in parallel by multi-shuttles. Thus, new methods need to be discovered to solve task scheduling problem in the DMSWS.
Research on the carbon emissions of warehousing systems has gradually attracted more and more academic and industrial attention. In the literature reviews, most researchers have focused on carbon emissions in warehousing design and path optimization and storage policies [8][9][10][11]. Some scholars have established novel models to study the sustainability of warehousing systems [12,13]. In terms of task scheduling, there have been few studies on carbon emissions. These studies have only considered the carbon emissions of the equipment during operations. To fill this research gap, it is not only required to calculate the carbon emissions of equipment during operations, but also the carbon emissions of the equipment during waiting processes or free time, such that carbon emissions can be investigated more comprehensively. In addition, in this study, an optimization model is established to study the relationship between carbon emissions and total working time. It cannot be simply considered (as in previous studies) that the longer the total working time, the higher the carbon emissions.
The operation of AS/RS and its variant is similar to that of the DMSWS and some of its methods can be applied, but the theoretical study on DMSWS is very scarce. In addition, most previous studies on task scheduling in AS/RS and its variant have focused on TSP and its solution algorithms. However, this case differs from DMSWS because multiple shuttles perform storage and retrieval tasks in parallel. To improve the system efficiency, how the different equipment can be scheduled and made to work together efficiently are essential.
In this paper, a time sequence mathematical model based on rearrangement operation for DMSWS is established, which helps us get accurate total time for retrieval operations. To evaluate carbon emissions more comprehensively, a carbon emissions model that includes working and waiting or idle process of equipment is proposed. Moreover, we provide a multi-objective optimization model for task scheduling by optimizing the total working time, waiting time and carbon emissions. In other words, the purpose of this paper is to find out a better operation execution solution for minimizing the total working time and carbon emissions.
The remainder of this paper is organized as follows. In Section 2, the literature review is presented. The structure of DMSWS and some assumptions are introduced in Section 3. The cycle time model and carbon emissions model are proposed in Section 4, while Section 5 presents the NSGA-II to solve the task scheduling problem. In Section 6, an actual engineering project is adopted to validate the model and its algorithm. Finally, the conclusions of this research and future research topics are given in Section 6.

Literature Review
In this section, studies on AS/RS and its variants (SBS/RS and AVS/RS) are introduced in the first part. Then, the studies relating to the carbon emissions of warehousing systems are presented. Finally, the studies of warehousing systems we review herein are summarized, as shown in Table 1. To date, scholars have carried out a series of studies on the task scheduling problem of AS/RS; however, the research object in most studies is the single-deep SR. Based on the assumption that the single-deep SR is "square-in-time" (SIT), Hausman et al. [2] first formulated the travel-time models for AS/RS with single-deep SR. In addition, the system efficiency and throughput of the three storage policies (full turnover policy, random policy, and class-based policy) have been compared, and the optimization results showed that the class-based storage policy is the most effective. Bozer et al. [3] built a classical analytical model to compare single-and dual-command strategies for AS/RS under the random storage policy. Additionally, various I/O point strategies and alternative dwell-point strategies were examined through actual cases for AS/RS. This research method is still Processes 2021, 9, 41 4 of 21 being used in travel-time research. According to the actual movement characteristics of the automation equipment, Hwang et al. [6] formulated continuous analytical models of the travel-time. Kouvelis et al. [7] and Bortolini et al. [14] established analytical models to evaluate the expected time for both single-and dual-command cycles. Furthermore, Regattieri et al. [15] adopted a simulation model to find the best dwell-point policy under different driving rules.
Multi-deep AS/RS is another research direction of AS/RS, and some scholars have conducted research in this direction. Hu et al. [16] presented a travel-time model of splitplatform SR based on the dwell-point policy and gave the optimal design guidelines of the split-platform SR. Sari et al. [17] studied a multi-deep AS/RS with special structure and proposed closed-form travel-time expressions. To the best of our knowledge, the 3-dimensional (3D) compact AS/RS was first studied by De Koster et al. [18], where a closed-form expression was given to evaluate the expected retrieval travel-time for a singlecommand cycle. As an extension of previous research, Yu et al. [19] obtained the expected single-command cycle time based on the full turnover policy. Double-deep AS/RS is another research highlight. Lerher et al. [20] presented a new analytical travel-time model for double-deep AS/RS, which can calculate the cycle time more accurately. Moreover, the actual movement characteristics of the automation equipment and the rearrangement operations during the retrieval process were considered. Xu et al. [4] used the single-shuttle as the benchmark model to analyze the improvement of the dual-shuttle and developed a "quadruple command" analytical model to evaluate the system efficiency. According to three different rearrangement rules, Xu et al. [21] proposed cycle time models under singleand dual-command. Their study also illustrated the relationship between the fill-grade factor and the total cost. Wang et al. [5] and Lerher [22] studied task scheduling models for double-deep SR. Moreover, research on AS/RS has been widely conducted; Roodbergen et al. [23] and Gu et al. [24] have presented literature reviews in this area.
Apart from AS/RS, some scholars have been conducted research on its variants, including the shuttle-based storage and retrieval system (SBS/RS) and the autonomous vehicle storage and retrieval system (AVS/RS). To the best of our knowledge, Malmborg [25] first proposed the concept of AVS/RS, formulating a system efficiency analysis model based on the features of autonomous vehicles. Fukunari et al. [26] presented a new travel-time model for AVS/RS, comparing the system efficiency with that of the traditional AS/RS. Heragu et al. [27] proposed analytical models based on an open queuing network (OQN) and analyzed the system throughput of the AVS/RS. Marchet et al. [28] presented an analytical model to evaluate the system efficiency of AVS/RS. Ekren and Cai et al. and Zou et al. [29][30][31] studied the system efficiency and throughput, according to different queuing network models. In addition, Wang et al. [32] investigated the task sequencing problem in SBS and proposed a sequence time model of task operation based on the actual movement characteristics of automation equipment. Lerher et al. [33] conducted research on SBS/RS and presented the cycle time model to evaluate the travel-time. Tappia et al. [34] developed novel queuing networking models to calculate the system efficiency of both single-and multi-tier SBS/RS.
According to the global development trend of supply chain and warehousing systems, the operations of systems not only consider the travel-time and total cost but have also come to consider the carbon emissions and environmental aspects of the automation equipment. As the existing AS/RS analysis models are suitable for the already well-known objectives (minimum total cost, minimum expected time, and maximum efficiency), Lerher et al. [1] proposed a new energy efficiency model based on the analysis of the mini-load AS/RS. Their research results showed that this model can effectively reduce the energy consumption, thereby reducing carbon emissions. On the basis of new policies for storage assignment, Bortolini et al. [35] presented a travel-time and energy consumption bi-objective model for aisle-captive AS/RS with single-deep SR. This model was based on the joint minimization of the cycle time and the energy consumption for cranes to store and retrieve the SKU. Ene et al. [11] adopted a genetic algorithm to minimize the energy consumption with a proper storage policy other than service time. For the unit-load multiple rack AS/RS, Nia et al. [8] proposed a dual-command "dynamic sequencing" method. In order to estimate greenhouse gas (GHG) emissions more comprehensively, limitations on the available time of all automation equipment and the total allowable GHG emissions generated by all facilities are calculated in this model. In addition, the penalty cost, tax cost, and discount of the produced GHG emissions are considered. Borovinsek et al. [9] proposed a multi-objective optimization model from the perspective of an environmentally efficient SBS/RS design, in which the total cost, average throughput time, and energy consumption were considered. Torabizadeh et al. [12] described a method to identify and weigh indicators that assess sustainability in a sustainable warehouse management system using structural equation modeling. Hao et al. [13] developed a conceptual model using a technology-organization-environment framework to investigate the factors which influence logistics firms to adopt green technology. To conduct energy evaluation, Emanuele et al. [10] analyzed the energy recovery system of a deep-lane AVS/RS. The results showed that the energy recovery system significantly saved energy in the deep-lane AVS/RS.
The previous studies focusing on single-deep AS/RS, multi-deep AS/RS, and singledeep AVS/RS and SBS/RS are shown in Table 1. However, no study has investigated double-deep multi-tier shuttle storage warehousing systems, even though they have high floor utilization and high equipment utilization. To fill this gap, we study the task sequencing problem in DMSWS and consider the carbon emissions of equipment, including the waiting and idle processes.

System Description and Modelling Preparation
As a new type of automated warehousing system, a double-deep multi-tier shuttle warehousing system has unique structure. The component units of a double-deep multi-tier shuttle warehousing system are the double-deep SR, shuttles, lifts, and control system. A schematic diagram of a DMSWS is shown in Figure 1, where the lift moves back and forth between the I/O point and tier for vertical movement of the required SKU. The shuttle is used for horizontal movement for the storage and retrieval tasks in each storage tier. All storage and retrieval tasks in this warehousing system are performed by vertical movement of the lift and horizontal movement of the shuttle. The shuttle and lift are parked at the buffer position and I/O point, respectively, when they are idle; these processes also generate carbon emissions. Compared with the traditional AS/RS, this combination of shuttle and lift consume more time and generate more carbon emissions for a single task. However, the parallel movement of multi-tier shuttles can improve the operating efficiency and reduce system carbon emissions, to a certain extent.
Processes 2021, 9, x FOR PEER REVIEW 6 of 23 single task. However, the parallel movement of multi-tier shuttles can improve the operating efficiency and reduce system carbon emissions, to a certain extent. According to the structure of a double-deep SR, it has a higher space utilization than a single-deep SR. However, rearrangement operations will increase the retrieval task time, thereby reducing the system efficiency. According to the introduction of the FEM 9.851 guideline, the rearrangement of the blocking SKU will appear during retrieval tasks in DMSWS. The procedure of the rearrangement of the blocking SKU is shown in Figure 2. The SKU (labelled by 'b') in the first lane blocks the SKU to be retrieved in the second lane (labelled by 'a'). Thereby, the shuttle needs to transport the blocking SKU to the nearest empty storage location and send the required SKU to the buffer position. Lastly, the block- According to the structure of a double-deep SR, it has a higher space utilization than a single-deep SR. However, rearrangement operations will increase the retrieval task time, thereby reducing the system efficiency. According to the introduction of the FEM 9.851 guideline, the rearrangement of the blocking SKU will appear during retrieval tasks in Processes 2021, 9, 41 6 of 21 DMSWS. The procedure of the rearrangement of the blocking SKU is shown in Figure 2. The SKU (labelled by 'b') in the first lane blocks the SKU to be retrieved in the second lane (labelled by 'a'). Thereby, the shuttle needs to transport the blocking SKU to the nearest empty storage location and send the required SKU to the buffer position. Lastly, the blocking SKU is transported back to the original location. This operation is called rearrangement of the blocking SKU. According to the structure of a double-deep SR, it has a higher space utilization than a single-deep SR. However, rearrangement operations will increase the retrieval task time, thereby reducing the system efficiency. According to the introduction of the FEM 9.851 guideline, the rearrangement of the blocking SKU will appear during retrieval tasks in DMSWS. The procedure of the rearrangement of the blocking SKU is shown in Figure 2. The SKU (labelled by 'b') in the first lane blocks the SKU to be retrieved in the second lane (labelled by 'a'). Thereby, the shuttle needs to transport the blocking SKU to the nearest empty storage location and send the required SKU to the buffer position. Lastly, the blocking SKU is transported back to the original location. This operation is called rearrangement of the blocking SKU. The fundamental assumptions of a DMSWS are presented, in order to analyze task scheduling and system efficiency in subsequent sections, as follows: 1. The random storage strategy is used. Under this strategy, a SKU will be placed in each empty storage location with the same probability. 2. The automation equipment follow the first-come-first-served (FCFS) rule. There are no emergency orders and every retrieval task shares the same priority level. 3. Shuttles and lifts abide by the dwell-point strategy of point-of-service-completion (POSC). The shuttle remain at the buffer position in each storage tier and the lift dwells at the I/O point. 4. A single-command cycle is considered. The DMSWS will not accept the storage task when the system is performing a retrieval task.
According to the system description and assumptions, retrieval tasks are considered in this study. These tasks represent the most critical operations of a DMSWS, as they directly influence the service level of the system. The operation of multi-tasking in the Step one SKU to retrieval Step two Step three Step four a b The blocking SKU The fundamental assumptions of a DMSWS are presented, in order to analyze task scheduling and system efficiency in subsequent sections, as follows: 1.
The random storage strategy is used. Under this strategy, a SKU will be placed in each empty storage location with the same probability.

2.
The automation equipment follow the first-come-first-served (FCFS) rule. There are no emergency orders and every retrieval task shares the same priority level.

3.
Shuttles and lifts abide by the dwell-point strategy of point-of-service-completion (POSC). The shuttle remain at the buffer position in each storage tier and the lift dwells at the I/O point.

4.
A single-command cycle is considered. The DMSWS will not accept the storage task when the system is performing a retrieval task.
According to the system description and assumptions, retrieval tasks are considered in this study. These tasks represent the most critical operations of a DMSWS, as they directly influence the service level of the system. The operation of multi-tasking in the DMSWS, as shown in Figure 3, relies on collaboration between the horizontal movement of shuttles and the vertical movement of lifts. The operation of retrieval tasks is carried out as follows: Step 1: When the control system releases the retrieval task, the shuttle responds to the control system if it is free at that time.
Step 2: If there is no rearrangement operation, the shuttle moves directly to the retrieval position, otherwise, the shuttle needs to complete the rearrangement operation first.
Step 3: The shuttle picks up the SKU and transports it to the buffer position, where the shuttle and the lift complete the handover of the SKU. Meanwhile, the shuttle sends a retrieval request to the lift through the control system.
Step 4: Based on the FCFS rule, the lift moves from the I/O point to the designated tier when it becomes free.
Step 5: The lift picks up the SKU and it returns to the I/O point to unload the SKU. Meanwhile, the shuttle is released by the control system for the next retrieval task.
In a DMSWS, the control system receives a series of retrieval tasks in a certain time window, some of which are performed simultaneously by the multi-tier shuttle, then all retrieval tasks are completed by the lift in turn. Although the shuttle and the lift abide by the FCFS rule, the performance of a DMSWS can be summarized as the sequential transfer of lifts and parallel retrieval of shuttles. Therefore, the task scheduling problem model of a DMSWS differs significantly from the TSP model, based on the above analysis. first.
Step 3: The shuttle picks up the SKU and transports it to the buffer position, where the shuttle and the lift complete the handover of the SKU. Meanwhile, the shuttle sends a retrieval request to the lift through the control system.
Step 4: Based on the FCFS rule, the lift moves from the I/O point to the designated tier when it becomes free.
Step 5: The lift picks up the SKU and it returns to the I/O point to unload the SKU. Meanwhile, the shuttle is released by the control system for the next retrieval task.  In a DMSWS, the control system receives a series of retrieval tasks in a certain time window, some of which are performed simultaneously by the multi-tier shuttle, then all retrieval tasks are completed by the lift in turn. Although the shuttle and the lift abide by the FCFS rule, the performance of a DMSWS can be summarized as the sequential transfer of lifts and parallel retrieval of shuttles. Therefore, the task scheduling problem model of a DMSWS differs significantly from the TSP model, based on the above analysis.

Models of Double-Deep Multi-Tier Warehousing Systems
In this section, according to the analysis of the retrieval task process, the cycle time model and carbon emissions model are established. Then, a Pareto optimization strategy is designed, based on the three optimization goals.

Cycle Time Model of the DMSWS
The operation of a single retrieval task was analyzed, as shown in Figure 4. The consecutive superscript Q is adopted to depict several critical times in the entire period. Particularly, Q k is defined as the finishing time for the retrieval task k; Q 0 k indicates the time when the control system releases the retrieval task k; Q 1 k presents the time when the shuttle has loaded the SKU and moves to the buffer position; Q 2 k shows the time when the lift responds to the retrieval request; Q 3 k is the moment when the shuttle and the lift complete the delivery of SKU and release the shuttle; and Q 4 k signifies the time when the lift moves to the I/O point. In this section, according to the analysis of the retrieval task process, the cycle time model and carbon emissions model are established. Then, a Pareto optimization strategy is designed, based on the three optimization goals.

Cycle Time Model of the DMSWS
The operation of a single retrieval task was analyzed, as shown in Figure 4. The consecutive superscript Q is adopted to depict several critical times in the entire period. Particularly, k Q is defined as the finishing time for the retrieval task k; 0 k Q indicates the time when the control system releases the retrieval task k;  Based on the FCFS rule, there may be two periods of stagnation during the execution of retrieval tasks, which may be affect the operation time and the system efficiency. Shuttle idling (waiting) is the first kind of stagnation period, which is indicated by k w . This happens when the shuttle has moved to the buffer position and sends a retrieval request to the lift while it is serving the previous task. On the other hand, the lift idle (free) time is induced when lift finishes the previous task and there is no retrieval request from the control system; this is denoted by k f . According to our analysis, the following two situations will occur during each retrieval task period.
1. If the lift has free time in retrieval task k, then 2. If the has shuttle waiting time in retrieval task k, then Based on the FCFS rule, there may be two periods of stagnation during the execution of retrieval tasks, which may be affect the operation time and the system efficiency. Shuttle idling (waiting) is the first kind of stagnation period, which is indicated by w k . This happens when the shuttle has moved to the buffer position and sends a retrieval request to the lift while it is serving the previous task. On the other hand, the lift idle (free) time is induced when lift finishes the previous task and there is no retrieval request from the control system; this is denoted by f k . According to our analysis, the following two situations will occur during each retrieval task period.

1.
If the lift has free time in retrieval task k, then f k = Q 1 k − Q k−1 .

2.
If the has shuttle waiting time in retrieval task k, then w k = Q k−1 − Q 1 k . p k and q k are assumed to be the probabilities of the shuttle and the lift idling, respectively: Thus, the total shuttle waiting time w and the total lift free time f can be calculated as follows: Based on the flow chart for multi-tasking, the shuttle waiting time and lift free time will occur in a single retrieval task, but all retrieval tasks will be completed by the lift in the end. Therefore, according to the vertical movement of the lift, the operation of all retrieval tasks can be divided into three stages, as illustrated in Figure 5: 1.
Stage one. The shuttle moves to the target location to load the SKU, then returns to the buffer position and transmits a retrieval request to the lift. A rearrangement operation may be executed. Based on the flow chart for multi-tasking, the shuttle waiting time and lift free time will occur in a single retrieval task, but all retrieval tasks will be completed by the lift in the end. Therefore, according to the vertical movement of the lift, the operation of all retrieval tasks can be divided into three stages, as illustrated in Figure 5  According to the aforementioned analysis and Figure 5, the entire operation process of retrieval tasks is regarded as a flow-shop scheduling problem in parallel machines mode. Therefore, the total retrieval operation time, total T is equal to the sum of the time taken for stage one of the retrieval task 1 ( 1 1 s ) and the total lift working time ( total s ).
The consumed time for retrieval k in stage one, which is denoted by , can be calculated as follows:  According to the aforementioned analysis and Figure 5, the entire operation process of retrieval tasks is regarded as a flow-shop scheduling problem in parallel machines mode. Therefore, the total retrieval operation time, T total is equal to the sum of the time taken for stage one of the retrieval task 1 (s 1 1 ) and the total lift working time (s total ). The consumed time for retrieval k in stage one, which is denoted by s 1 k , can be calculated as follows: where e 1 k is the time for horizontal movement of the shuttle; e 1 represents the time for loading the required SKU; e d indicates the time for a rearrangement operation; and F ki is a decision variable: If F ki = 1, the shuttle must carry out a rearrangement operation.
In stage two, the operation time for retrieval task k can be obtained by the following formula.
where e 2 k presents the time for the lift to move from the I/O point to the designated tier; and e 2 indicates the time for the SKU handover between the shuttle and the lift.
The operation time for retrieval task k in stage three can be obtained as follows: where e 3 represents the time to unload the required SKU. For all retrieval tasks, the total lift working time can be calculated as follows: According to the analysis, According to the flow-shop scheduling problem with parallel machines mode, the total retrieval operation time T total can be derived as: In the DMSWS, e 1 , e 2 and e 3 are all constant variables, and the motion variables do not consider the effects of acceleration and deceleration on the movement of the shuttle and the lift. v s and v l are the average velocity of the shuttle and the lift. (x, y, z) indicates the original storage position for retrieval task k. The time taken for horizontal movement of the shuttle and for vertical movement of the lift can be derived as follows: Based on [21,22], the expected value of the rearrangement operation distances of the DMSWS can be calculated by Equation (12): where L represents the length of a single storage position and α is the fill-grade factor of the DMSWS. The time of a rearrangement operation can be calculated by Equation (13):

Carbon Emissions Model
The amount of carbon emissions, TE, can be calculated based on the energy consumption of the warehousing system. In this section, the carbon emissions during the retrieval operation are considered, as well as those during the waiting and free processes of the equipment. According to descriptions of practitioners in the field and our investigations, this consideration is very appropriate.
When the operation of the retrieval task k is completed, the total carbon emissions of all outbound equipment in the DMSWS can be calculated by: PF k = f k ·P lw (16) where TE k indicates the total carbon emissions generated by retrieval task k; PT k represents the total energy consumption for the retrieval task k; f e is the conversion factor of greenhouse gas (GHG); PW k indicates the energy consumption of the shuttle in the waiting process; P sw is the power of the shuttle in the waiting process; PF k represents the energy consumption of the lift during the free process; P lw indicates the power of the lift during the free process; PC k is the energy consumption of retrieval operation for retrieval task k; PC sk indicates the energy consumption of the shuttle for the retrieval operation; and PC lk indicates the energy consumption of the lift for the retrieval operation.
In stage one, the energy consumption of the shuttle and the lift for retrieval task k are described by PC 1 sk and PC 1 lk , which can be calculated as follows: where P s f represents the no-load power of the shuttle and P sm indicates the load power of the shuttle. In stage two, the energy consumption of the shuttle and the lift are denoted by PC 2 sk and PC 2 lk , expressed as follows: where P l f indicates the no-load power of the lift. In stage three, the energy consumption of the shuttle and the lift are calculated by: where PC 3 sk represents the energy consumption of the shuttle, PC 3 lk represents the energy consumption of the lift, and P lm indicates the load power of the lift.
Thus, the energy consumption of the retrieval operation for retrieval task k can be obtained according to the following formula: Based on the analysis above, the total carbon emissions TE can be obtained as follows:

Pareto Optimization
When multi-objective functions have an inverse relationship, the multi-objective optimization problem has countless solutions (non-dominated solutions), instead of a unique optimal solution. After obtaining the optimized results, the researcher then selects one or more acceptable solutions among all non-dominated solutions.
The aforementioned analysis of DMSWS indicates that the total retrieval operation time, the total carbon emissions and the total waiting time of shuttles reflect the efficiency of the DMSWS from different perspectives. During the operation of retrieval tasks, the influences of mutual promotion and contradiction are presented. The purpose of this study is to optimize three objective functions-the total retrieval operation time, the total waiting time of the shuttles, and the total carbon emissions-of the DMSWS. Thus, the research objective is to minimize three objective functions through the re-scheduling execution sequence in a certain time window. The modelling objectives are described as follows: where w, T total , and TE can be derived from Equations (3), (9) and (25), respectively, and x i represents the scheduling execution sequence.

Solution Algorithm
The essence of the task scheduling problem is to find the best solution among all permutations and combinations of the operation execution sequence. The elitist nondominated sorting genetic algorithm (NSGA-II) was adopted to solve the task scheduling problem of DMSWS. Since the NSGA-II was first proposed by Deb (2002), it has been widely used for solving multi-objective optimization problems (MOOPs). The key steps of NSGA-II are shown in Figure 6, which presents the detailed program flow. research objective is to minimize three objective functions through the re-scheduling execution sequence in a certain time window. The modelling objectives are described as follows: where w , total T , and TE can be derived from Equations (3), (9), and (25), respectively, and i x represents the scheduling execution sequence.

Solution Algorithm
The essence of the task scheduling problem is to find the best solution among all permutations and combinations of the operation execution sequence. The elitist non-dominated sorting genetic algorithm (NSGA-II) was adopted to solve the task scheduling problem of DMSWS. Since the NSGA-II was first proposed by Deb (2002), it has been widely used for solving multi-objective optimization problems (MOOPs). The key steps of NSGA-II are shown in Figure 6, which presents the detailed program flow. (1) Encoding and decoding First, the encoding and decoding process are introduced, which are very important. As just one problem is involved, the chromosome contains only one important vector (i.e., the task sequence vector). Thus, natural number coding is used for this problem. The task sequence vector has M gene elements, in which each element indicates a retrieval task. The position of each element in the task sequence vector represents the order of execution of the retrieval task. Thus, each chromosome represents an operation execution solution. An example with 14 tasks is depicted in Figure 7. (1) Encoding and decoding First, the encoding and decoding process are introduced, which are very important. As just one problem is involved, the chromosome contains only one important vector (i.e., the task sequence vector). Thus, natural number coding is used for this problem. The task sequence vector has M gene elements, in which each element indicates a retrieval task. The position of each element in the task sequence vector represents the order of execution of the retrieval task. Thus, each chromosome represents an operation execution solution. An example with 14 tasks is depicted in Figure 7. (2) Fitness function To estimate the optimization solution, the objectives to be minimized are transformed to form a maximum optimization problem: (3) Non-dominated sort In each generation, the population is sorted on the basis of the non-dominated Pareto optimality. The initial Pareto front set is not dominated and all chromosomes are transferred from the population. We continue to confirm the next Pareto front set in the remaining population until all chromosomes are classified.
(4) Crowding distance The crowding distance between adjacent individuals in each Pareto front set is calculated, after the non-dominated sort is completed. The crowding distance is equal to the cumulative sum of the fitness values of the chromosomes I − 1 and i + 1: where M presents the number of optimization goals; the fitness values of the optimization goal m for chromosomes I − 1 and i + 1 are indicated by   (2) Fitness function To estimate the optimization solution, the objectives to be minimized are transformed to form a maximum optimization problem: In each generation, the population is sorted on the basis of the non-dominated Pareto optimality. The initial Pareto front set is not dominated and all chromosomes are transferred from the population. We continue to confirm the next Pareto front set in the remaining population until all chromosomes are classified.
(4) Crowding distance The crowding distance between adjacent individuals in each Pareto front set is calculated, after the non-dominated sort is completed. The crowding distance is equal to the cumulative sum of the fitness values of the chromosomes I − 1 and i + 1:  In a population P, the elitist parent set is obtained on the basis of the non-dominated sort crowding distance and the Pareto front rank. H p indicates the elitist set, which has the largest crowding distance and the smallest front rank. The remaining individuals, set by R P , produce the next generation's individuals according to tournament selection. Finally, the genetic operators are the same as in GA, as shown in Figure 8, and the termination condition is usually represented by a fixed number of generations reached.

Case study
In this section, the logistics distribution center of a food company is chosen as an experimental case, in order to validate the proposed optimization model. This distribution center uses a conveyor system and DMSWS. All required SKUs are gathered by the warehouse management system (WMS) within a certain time window, where the WMS obtains the storage positions of all retrieval tasks according to the design of the DMSWS. Then, according to the optimization results, the retrieval tasks are re-scheduled by the warehouse control system (WCS). Meanwhile, the WCS assigns each task to relevant shuttles and the lift, based on the optimized execution sequence. Each required SKU is transferred to the I/O point by the co-operation of shuttles and lift.
As shown in Figure 9, the DMSWS of the logistics distribution center is composed of four parts: each part includes 12 Table 2.

Case study
In this section, the logistics distribution center of a food company is chosen as an experimental case, in order to validate the proposed optimization model. This distribution center uses a conveyor system and DMSWS. All required SKUs are gathered by the warehouse management system (WMS) within a certain time window, where the WMS obtains the storage positions of all retrieval tasks according to the design of the DMSWS. Then, according to the optimization results, the retrieval tasks are re-scheduled by the warehouse control system (WCS). Meanwhile, the WCS assigns each task to relevant shuttles and the lift, based on the optimized execution sequence. Each required SKU is transferred to the I/O point by the co-operation of shuttles and lift.
As shown in Figure 9, the DMSWS of the logistics distribution center is composed of four parts: each part includes 12 aisle, 8 tiers, and 20 storage columns, and there are 16,640 total storage positions. The width of each aisle is W a = 1.5 m, and the length and the width of a single storage position are L = 1.2 m and W = 1.0 m, respectively. The average velocity of a shuttle is v s = 1.5 m/s and the average velocity of a lift is v l = 2.0 m/s. According to the research of [36], the conversion factor for GHG is f e = 974 gCO 2 /kwh. The remaining parameter values are given in Table 2.

Case study
In this section, the logistics distribution center of a food company is chosen as an experimental case, in order to validate the proposed optimization model. This distribution center uses a conveyor system and DMSWS. All required SKUs are gathered by the warehouse management system (WMS) within a certain time window, where the WMS obtains the storage positions of all retrieval tasks according to the design of the DMSWS. Then, according to the optimization results, the retrieval tasks are re-scheduled by the warehouse control system (WCS). Meanwhile, the WCS assigns each task to relevant shuttles and the lift, based on the optimized execution sequence. Each required SKU is transferred to the I/O point by the co-operation of shuttles and lift.
As shown in Figure 9, the DMSWS of the logistics distribution center is composed of four parts: each part includes 12 Table 2.   The WCS gathered 60 retrieval tasks in a certain time window. The original storage positions of all retrieval tasks are listed in Table 3. The operation of retrieval tasks was based on unit load. Thus, the stock was not opened during the stage of retrieval operations.  1  3  14  2  1  31  1  19  5  2  2  2  19  2  2  32  1  8  5  2  3  4  9  2  2  33  4  1  5  1  4  2  3  2  1  34  3  14  5  1  5  2  3  2  2  35  3  2  5  2  6  3  14  2  2  36  2  11  5  1  7  1  8  3  2  37  3  2  6  2  8  1  5  3  2  38  3  14  6  2  9  4  18  3  1  39  4  7  6  1  10  3  17  3  2  40  1  9  6  2  11  2  12  3  2  41  2  15  6  1  12  4  18  3  2  42  3  14  6  1  13  1  8  3  1  43  1  11  7  2  14  3  3  4  1  44  4  20  7  1  15  2  20  4  2  45  1  3  7  1  16  4  13  4  1  46  2  8  7  2  17  4  7  4  2  47  3  17  7  1  18  2  9  4  1  48  3  6  7  1  19  2  4  4  2  49  3  17  7  2  20  3  11  4  2  50  1  3  7  2  21  1  8  4  1  51  4  20  7  2  22  1  16  4  2  52  2  8  7  1  23  3  3  4  2  53  2  9  8  2  24  4  13  4  2  54  3  4  8  1  25  3  11  4  1  55  1  12  8  1  26  2  9  4  2  56  4  18  8  1  27  1  8  4  2  57  2  9  8  1  28  4  1  5  2  58  1  3  8  2  29  3  14  5  2  59  1  12  8  2  30  2  11  5  2  60  3  4  8  2 According to the task scheduling optimization model and NSGA-II, the WCS determined the best execution sequence. The NSGA-II was run using the MATLAB 2016b software, which operated on a control computer with an Intel(R) Core(TM) i7-7700HQ CPU @ 2.80 GHz and 8.00 GB RAM. The parameter values for NSGA-II were selected based on previous analysis and the experience of researchers engaged in the application and development of intelligent algorithms. The parameters directly impact the efficiency of the NSGA-II algorithm. This indicates that the parameters should be careful determined to ensure unprejudiced comparisons are conducted. The factor levels, which are determined by a large number of tests, are shown in Table 4. Because there are two parameters and three factor levels for each parameter, an orthogonal array L 1 (3 2 ) is adopted. The ratio of non-dominated solutions in the final solution set is considered as the response variable. Each case is tested 10 times and the average value is shown in Figure 10. Therefore, the best parameter settings are: crossover probability is 0.90 and mutation probability is 0.1. This parameter setting is used in the next experiments. The size of the initial population was set to 200, and the maximum generation was set to 300. CPU @ 2.80 GHz and 8.00 GB RAM. The parameter values for NSGA-II were selected based on previous analysis and the experience of researchers engaged in the application and development of intelligent algorithms. The parameters directly impact the efficiency of the NSGA-II algorithm. This indicates that the parameters should be careful determined to ensure unprejudiced comparisons are conducted. The factor levels, which are determined by a large number of tests, are shown in Table 4. Because there are two parameters and three factor levels for each parameter, an orthogonal array L1 (3 2 ) is adopted. The ratio of non-dominated solutions in the final solution set is considered as the response variable. Each case is tested 10 times and the average value is shown in Figure 10. Therefore, the best parameter settings are: crossover probability is 0.90 and mutation probability is 0.1. This parameter setting is used in the next experiments. The size of the initial population was set to 200, and the maximum generation was set to 300. The hypervolume convergence curve is shown in Figure 11, and the slope tends to stabilize when the number of generation reaches 200. The convergence curve of optimization objectives is shown in Figure 12. According to the comparison of the minimum and average of optimization objective, the convergence of solution space is achieved. The results of the Pareto optimization are indicated in Figure 13 and Table 5. Figure 13a-c show the dependencies between pairs of objective functions, while Figure  13d shows the final Pareto optimal frontier. According to Figure 13a-c, the total time and the total waiting time were negatively correlated, while the carbon emissions and total The hypervolume convergence curve is shown in Figure 11, and the slope tends to stabilize when the number of generation reaches 200. The convergence curve of optimization objectives is shown in Figure 12. According to the comparison of the minimum and average of optimization objective, the convergence of solution space is achieved. The results of the Pareto optimization are indicated in Figure 13 and Table 5. Figure 13a-c show the dependencies between pairs of objective functions, while Figure 13d shows the final Pareto optimal frontier. According to Figure 13a-c, the total time and the total waiting time were negatively correlated, while the carbon emissions and total time were positively correlated. The diagram in Figure 13d indicates the dependencies between three objective functions (the total time, carbon emissions, and the total waiting time). As the total time continues to increase, carbon emissions gradually increase and waiting time gradually decreases. correlated. The diagram in Figure 13d indicates the dependencies between three objective functions (the total time, carbon emissions, and the total waiting time). As the total time continues to increase, carbon emissions gradually increase and waiting time gradually decreases. Figure 11. The hypervolume convergence curve. The total time (s) Waiting time (s) Figure 11. The hypervolume convergence curve.
Processes 2021, 9, x FOR PEER REVIEW 17 of 23 correlated. The diagram in Figure 13d indicates the dependencies between three objective functions (the total time, carbon emissions, and the total waiting time). As the total time continues to increase, carbon emissions gradually increase and waiting time gradually decreases. Figure 11. The hypervolume convergence curve. According to the analysis, the longer the total waiting time of shuttles, the shorter the free time of the lift. This means that the lift is in working condition for a long time, such that the total time of retrieval requests is shorter. On the contrary, the shorter the total waiting time of shuttles, the longer the free time of the lift, which increases the total time  According to the analysis, the longer the total waiting time of shuttles, the shorter the free time of the lift. This means that the lift is in working condition for a long time, such that the total time of retrieval requests is shorter. On the contrary, the shorter the total waiting time of shuttles, the longer the free time of the lift, which increases the total time of retrieval requests and reduces the utilization of the lift. The carbon emissions are composed of three parts: the carbon emissions of the shuttle in the waiting process, the carbon emissions of the lift during the free process, and the carbon emissions of retrieval operation. The carbon emissions of the retrieval operation is the main component and, so, as the retrieval operation time increases, the total carbon emissions also increase.
The optimization results provided a unique ID to indicate the differences between the non-dominated solutions. Table 4 shows the optimal solutions, sorted by the total time in the increasing direction. The best solution had a total time of 801.53 s, which also had the maximum value of the waiting time. The worst solution was that with the lowest value of the waiting time and total time of 1140.32 s; 338.79 s longer than the best solution. In terms of carbon emissions, the best solution had carbon emissions of 5533.43 gCO2 and total time of 812.91 s, only 11.38 s longer than the minimum total time. There was a positive correlation between carbon emissions and the total time, but the correlation was not strong. The total time is composed of the working time and the free time of the lift. When an increase in the free time of the lift causes the total time to increase, the carbon emissions may not necessarily increase or even decrease. This is a vital finding of this paper.
We choose 3 chromosomes out of Pareto optimal solutions for detailed analysis, namely the chromosome of ID 1, the chromosome of ID 3 and the chromosome of ID 30. The values of ID 1 chromosome are 801.53 s, 243.40 s, and 5617.71 gCO 2 , respectively. It shows that the total time of chromosome 1 is the shortest and the system efficiency is the highest. The chromosome 3 has the smallest carbon emissions that is 5533.45 gCO 2 , which can help the warehousing system achieve sustainable development. The values of ID 30 chromosome are 1140.32 s, 20.55 s, and 6547.53 gCO 2 , respectively, which indicates that waiting time for shuttles is the shortest. It can help warehousing managers improve the utilization of shuttles.
The values of chromosome variables for each non-dominated solution have been enumerated in Table 6. According to the optimization results, a non-dominated solution was selected as an execution without priority. In fact, the operation of retrieval tasks considering the best solution as that with minimal total time for the final operation execution sequence. The solution of ID 1 was chosen as the final operation execution sequence in this real case; that is, the total time, the waiting time, and carbon emissions were 801.53 s, 243.40 s, and 5617.71 gCO 2 , respectively. Additionally, the solution ID 3 was selected as the final operation execution sequence, based on the minimal carbon emissions. Compared to ID 1, the total time was increased by 11.38 s but the carbon emissions were reduced by 84.26 g CO 2 , which is equivalent to the CO 2 absorbed by a broad-leaved forest of about 110 m 2 in one day. Table 6. Detailed execution solutions for non-dominated solution.

ID
Chromosome Variable Value (Operation Execution Solution) Table 6. Cont.

ID Chromosome Variable Value (Operation Execution Solution)
26 [35,7,50,17,11,42,36,13,28,33,58,44,24,23,40,20, The research results show that the Pareto optimization model can reduce carbon emissions and total working time, thereby reducing the operation cost and improving system efficiency. And this method has helped the logistics distribution center of the food company to improve the system efficiency. According to the sustainable development trend of warehousing system, the operation of systems not only considers the travel time but considers the carbon emissions and environment aspect on equipment. It also can help warehouse managers save a lot of operation costs. Meanwhile, for some equipment rental-type warehousing system, managers choose operation execution solution with a high utilization of shuttles, which can save some rental costs. Moreover, the method proposed in this paper can improve system efficiency and reduce carbon emissions without adding additional costs. This study provides a useful and a flexible method for task scheduling problem and this method has also been used in some practical implications. We believe this method will be more widely used in the future.

Conclusions
In this paper, the task scheduling problem was studied for DMSWS, which have high space utilization. A cycle time model based on the rearrangement operation and parallel retrieval of a multi-tier shuttle with sequence transfer to a lift was developed. According to the cycle time composition of a single retrieval task, the task scheduling problem was regarded as flow-shop scheduling model with parallel machines. In addition, a carbon emissions model was proposed, which considers both the storage/retrieval task operations and the idle processes of the equipment; this can help us to evaluate carbon emissions more comprehensively. An optimization model was derived, according to three correlated optimization objectives: the minimum total time, minimum carbon emissions, and minimum waiting time of shuttle. The NSGA-II algorithm was applied to solve the multi-objective optimization model. Finally, a real case was used to validate the effectiveness of the proposed optimization model. The results show that non-dominated solutions can be found effectively by proposed algorithm. Furthermore, ID 1 was selected as the execution sequence for the practical project, improving system efficiency while reducing carbon emissions.
In the future work, the impact of different storage policies (dedicated storage policy or class-based policy) can be considered. The devices acceleration and the procedure of the rearrangement operation are interesting direction for future research. Additionally, the influence of the solution algorithm should be considered. The results of multiple multi-objective optimization algorithms (SMS-EMOA or MOEA/D etc.) will be compared to find an algorithm that is more suitable for this problem.