Modelling Analysis of a Four-Way Shuttle-Based Storage and Retrieval System on the Basis of Operation Strategy

Abstract

In the context of sustainable development, this paper rationalises the outbound process of a four-way shuttle system with a focus on their modelling, performance evaluation and configuration using a parallel operation strategy to reduce resource waste, thus achieving sustainable development. The parallelism of the hoist and shuttle is innovatively incorporated into the four-way shuttle system, so the modelling content is divided into parallel and serial types. In the parallel operation strategy model, a separation–aggregation queueing network model is constructed, and the open-loop queueing network is innovatively solved using the maximum entropy method. In the serial operation strategy model, a semi-open-loop queuing network is constructed and solved using the geometric matrix method. By varying different parameters, the accuracy of the model is verified by Arena simulation with an error range of 10% or less, and the error of the system performance index calculation is reduced by 20% compared with the existing methods. Setting up 18 different sizes of shuttle systems provided a better performance than a single serial-operation strategy through the addition of parallel strategies, with an average reduction of 12.6% in the system response time and a minimum reduction of 1.8%. The conclusions of this paper were verified on the basis of an arithmetic case analysis.

1. Introduction

With the rapid development of information technology and e-commerce, green transportation, green consumption and other areas of high-quality development, green low-carbon lifestyle and production techniques are accelerating the adoption of the concept of sustainable development among the people. As a third source of profit, the logistics industry has gained unprecedented momentum with respect to development. Warehouse systems are the material and information transit centers of enterprise logistics systems and have a decisive influence on the operational efficiency and cost of the whole supply chain [1]. Four-way shuttle systems are the most advanced system currently available, and represent a highly automated, highly flexible and dense storage system, the appearance of which has improved the performance of storage systems, avoiding the defects of traditional automated warehouse equipment, which typically cover a large area, and possess low operational efficiency and poor flexibility. However, the lack of a good operating strategy will result in operational inefficiency, wasted resources and other problems, which are in opposition to the concept of sustainable development. For this reason, it is of great theoretical and practical significance to establish a theoretical model for a four-way shuttle system that is in line with reality, and to conduct research on the operation strategy and system performance. The existing scholarly research on queuing networks in storage systems is shown in

Table 1. Summary of research on AVS/RS systems.

Author	Warehousing System	Model Type	Modeling	Solving Method
Kuo et al. [2]	Multilayer AVS/RS	System performance analysis	NQN	Iterative method
Fukunari and Malmborg [3]	Multilayer AVS/RS	System performance analysis	OQN	Meanvalue analysis
Fukunari and Malmborg (2009) [4], Kuo et al. [5]	Multilayer AVS/RS	System performance analysis	CQN	Meanvalue analysis
Heragu et al. [6]	Monolayer AVS/RS	System performance and configuration analysis	OQN	Manufacturing system performance analyser
Ekren et al. [7]	Multilayer AVS/RS	Statistic analysis	SOQN	Simulation method
Ekren et al. [8]	Multilayer AVS/RS	System performance analysis	Simulation model	Simulation method
Guerrazzi et al. [9]	Multilayer AVS/RS	System performance analysis	Energy Model	Simulation method
Roy D et al. [10,11,12,13,14]	Monolayer AVS/RS	System performance analysis	SOQN	Decomposition method Iterative method, Decomposition method
Antonio et al. [15]	Multilayer AVS/RS	System performance analysis	Evaluation System Model	Decomposition method
Ekren et al. [16]	Multilayer AVS/RS	System performance analysis	Simulation Regression model	Decomposition method

.

For the study of dense storage systems, the establishment of a system model is an essential optimisation step. Scholars have mostly used the queuing network model; however, for complex shuttle systems, a single open-loop or closed-loop queuing network cannot meet the modelling needs under multilayer, multishuttle and multilift conditions. Therefore, it is also necessary to study the system’s performance in order to establish a model that reflects the actual operation mode of the system. At present, studies of storage system shuttle and hoist operation strategies mostly use a single serial strategy; however, parallel operation strategies have been employed in real scenarios, and have a significant positive impact on the system performance. AVS/RS systems, AS/RS systems and SBS/RS systems have been developed to analyse the impact of the operating strategy on the operational efficiency of the whole system [17,18,19]. It can be observed from the available literature that dense storage systems have been widely studied by scholars, and the research has mostly been centred on AVS/RS. The research directions of some scholars have mostly focused on system scheduling, storage strategy and sorting order optimisation, while other scholars have focused on system modelling, system performance evaluation and system configuration. In addition to there being less research on four-way shuttle systems, the research has the following limitations: (1) it mostly presents single operation strategies for the shuttle and hoist in the system; (2) it mostly presents single-level and single-aisle analyses; and (3) generality is ignored when building queuing networks [20,21,22].

Therefore, this study focuses on a four-way shuttle system with a parallel operation strategy; adopts system performance indices, such as equipment utilisation, system response time and external queue waiting time; establishes a separation–aggregation queuing network model and a SOQN model; and designs an algorithm based on the Coxian distribution and the MGM solution. The performance variation of the system is analysed under different operation strategies [23,24,25,26].

2. Problem Description

This paper studies a single-depth four-way shuttle system. The hoist can assist the four-way shuttle in achieving layer changes, but it can also carry goods to complete the entry and exit operations. Figure 1 shows a simple diagram, in which the shuttle car and layer change hoist can be configured according to the specific situation. Usually, one hoist serves multiple lanes and multiple shuttles within the system, and the shuttles move within the corresponding area. The possible movements of the four-way shuttle car are: straight ahead, lane change and layer change; straight-ahead and lane-change movements are performed by themselves, while layer-change movements are performed with the assistance of the layer change hoist.

The basic operation modes of the four-way shuttle system can be divided into inbound and outbound operations. As the inbound and outbound operations are logically similar, and the steps involved in the outbound operations include those of the inbound operations, this paper focuses on the outbound operations with the following flow:

In the outbound operation process, the FCFS strategy is used, and the WMS locates the layer, column and lane of the SKU according to the outbound tasks [27,28,29,30].

3. System Outbound Process Modelling

3.1. System Outbound Process Modelling Assumptions

To establish a reasonable model, the assumptions and unknown parameters first need to be established. First, the object of study in this paper is a four-way shuttle system with one hoist; the aisles it serves, as well as the shuttles it serves, consist of an overall multilayer, multi-aisle unit for mapping the performance of the whole system. The system uses one tote as the access unit, and only one tote can be accessed per bay, with the same bay size. For the sake of convenience, and without loss of generality, the following assumptions are made [31]:

• The shuttles and hoists obey the FCFS strategy.
• The main reason for studying the outbound process is that both the inbound and outbound operations are logically coordinated with the shuttle and hoist, while the outbound operation involves steps that encompass the inbound operation.
• The goods follow a random storage strategy, and the tasks are assigned using a random assignment strategy in which each unused space has the same probability of being used to store goods, and each shuttle and hoist has the same probability of being assigned to a task.
• The hoist and the shuttle follow the POSC (endpoint stay) strategy, which means that when the shuttle and the hoist complete a certain task, they stay at the end position of the task until the next task arrives. In this paper, the hoist will stay at the I/O location (system level 1 location) when it completes the pickup task, and the shuttle will stay in the cache area of the level on which the previous pickup task was located.
• The operation time of the shuttle and hoist obeys a general distribution.
• A single-cycle operation is adopted.

With assumption 5, it is established in this paper that the service time of the shuttle and hoist obey a general distribution, and the queuing network composed of the general distribution takes the form of a nonproduct solution; therefore, the general model for solving queuing networks is not suitable for use in this paper. However, some scholars have studied the approximate relationship between the general distribution and the Coxian-type distribution in storage systems (see [32]), which can be used to simplify the calculation. The Coxian distribution belongs to a special form of PH constructed from a series of exponential distributions, which can be applied in the approximate solution of a service time while obeying the general distribution [33].

3.2. Average Service Time Model for a Four-Way Shuttle System Server

To describe the four-way shuttle system based on the operation time model in detail, the parameters used in the model and their meanings are expressed as shown in

Table 2. Meaning of unknown parameters in the model.

Parameters	Meaning	Parameters	Meaning
T	Storage levels of shelves	$W_A$	Roadway width (m)
C	storage columns for shelves	$H$	Single storage rack height (m)
N	Number of tunnels served by a single hoist	$Vs$	The average speed of shuttle (m/s)
A	Number of longitudinal aisles of single-level shelves	${\epsilon }_s$	Average loading/unloading time for shuttles (t)
S	Total storage system capacity	${\gamma }_s$	Time delay caused by shuttle acceleration/deceleration (t)
L	Length of individual cargo space (m)	$Vl$	The average speed of hoist (m/s)
W	Width of individual cargo space (m)	${\epsilon }_l$	Average loading/unloading time of hoist (t)
M	Number of shuttles served by a single hoist	${\gamma }_l$	Time delay due to hoist acceleration/deceleration (t)

.

Set the pickup task at the t tier as a customer of class t. From assumption 3, the probability of serving a customer of class t is:

$$P\left(t\right)=\frac{1}{T}$$

(1)

Set the goods to be picked up in column$M_{ac}$ of aisle t, where$a=1,2,\dots ,N,c=1,2,\dots ,$ C; then, the probability of picking up the goods is:

$$P\left(M_{ac}\right)=\frac{1}{2·N·C}$$

(2)

The calculation of the equipment service time can be divided into two cases, depending on the availability of the idle four-way shuttles in the layer of goods to be picked up. These are the parallel operation strategy, where there are free four-way shuttles in the pickup layer, and the serial operation strategy, where there are no free four-way shuttles in the pickup layer.

3.2.1. Server Service Time Model Using Parallel Operation Strategy

First, the shuttle and the hoist need to carry out horizontal pickup operations and vertical target-level seating operations, respectively, constituting a parallel operation strategy in which the shuttle and the hoist complete their tasks at the same time. Subsequently, the shuttle carries the goods to the hoist, the hoist accepts the goods to complete the convergence, and finally, the shuttle carries the goods to level 1. The probability in this case P₁ is:

$$P_1=\frac{C_{T-1}^{M-1}}{C_T^M}$$

(3)

Therefore, the average service time of the server in vertical jobs, horizontal jobs and remaining jobs are calculated as follows:

7.

Operation service time of vertical operations

In accordance with assumption 4, this paper adopts the endpoint stay strategy. The hoist stays at the endpoint of the previous task, which is the initial position of system level 1; therefore, the movement distance of the vertical job is the vertical distance from the target level to level 1. Thus, the service time of vertical operation $T_{l1}\left(t\right)$is:

$$T_{l1}\left(t\right)=\frac{\left(t-1\right)H}{V_t}+2\left({\epsilon }_l+{\gamma }_l\right)$$

(4)

8.

Operating service time of horizontal operations

When the four-way shuttle is at the target t level, the movement distance can be divided into (1) moving from the hoistway entrance (buffer zone) to the target cargo location; and (2) returning to the hoistway entrance. The distance of movement includes the horizontal distance across the lane and the vertical distance along with the lane movement, which involves the cross-aisle movement of the shuttle. According to expert consultation, the number of cross-aisles $x$is determined by the number of aisles served by hoist N. If N is an even number, the location of the hoist is the $\frac{N}{2}$th aisle from left to right; if N is an odd number, the location of the hoist is the $\frac{N+1}{2}$th aisle from left to right, then the number of aisles $x$ spanned takes a range of values$0\le x\le \left|\frac{N}{2}\right|$ and the probability of crossing aisle $x$ is:

$$P\left(x\right)=\{\begin{array}{c}\frac{1}{N},x=0\\ \frac{1}{N},Nisevenandx=\frac{2}{N}\\ \frac{2}{N},Nisevenand0<x<\lfloor \frac{N}{2}\rfloor \\ \frac{2}{N},Nisoddand0<x<\lfloor \frac{N}{2}\rfloor \end{array}$$

(5)

The value of x is:

$$x=\left|\lceil \frac{N}{2}\rceil -a\right|$$

(6)

Then, the service time of the shuttle during this phase$T_S\left(M_{ac}\right)$ is:

$$T_s\left(M_{ac}\right)=\frac{2\left(x+\lceil \frac{c}{4}\rceil \right)W_A+2Lx+cW}{V_s}+{\epsilon }_s+4{\gamma }_s$$

(7)

9.

Operating service time of the remaining operations

The hoist moves from the target t level with cargo to level 1, and the operating service time $T_{ls}\left(t\right)$ is:

$$T_{ls}\left(t\right)=\frac{\left(t-1\right)H}{V_l}+2\left({\epsilon }_l+{\gamma }_l\right)$$

(8)

In this case, there is a parallel operation between the hoist movement in 1 and the shuttle movement in 2. While the shuttle receives the outbound task, it applies for the hoist to wait at the target level; if one side arrives first, it will wait for the other side until it completes the rendezvous at the hoistway entrance picking up the goods and unloading them with the hoist to level 1.

3.2.2. Service Time Model for the Server Using a Serial Operation Strategy

In this case, the hoist needs to carry out three operations: the first is to move to the layer with the idle shuttle, to pick up the trolley; the second is for the idle four-way shuttle to move to the target cargo layer when the hoist that is still in the task queue is not released; and the third is for the hoist to carry the goods to the first layer to complete the pickup task; at this time, the shuttle and the hoist are released. The probability in this case $P_2$ is:

$$P_2=1-P_1$$

(9)

Then, the average service times of the server for vertical jobs, horizontal jobs and remaining jobs are calculated as follows:

• Operation service time for vertical operations

The hoist is initially at level 1 under assumption 4, and the idle four-way shuttle is located at the position it was in at the end of the last job, which is in the cache of level t^′, on which the last job was located. Therefore, the service time of the movement of the hoist to the idle four-way shuttle$T_l\left(I{t}^{\prime }\right)$ is:

$$T_l\left(I{t}^{\prime }\right)=\frac{\left(t-1\right)H}{V_l}+2{\gamma }_l$$

(10)

The pickup target level is t, and the operating service time for moving the idle shuttle on the hoist to the target level $T_{ls}\left(t^{\prime }t\right)$is:

$$T_{ls}\left(t^{\prime }t\right)=\frac{\left|t^{\prime }-t\right|H}{V_l}+2\left({\epsilon }_l+{\gamma }_l\right)$$

(11)

The total vertical operation time $T_l\left(t^{\prime }\right)$is:

$$T_l\left(t^{\prime }\right)=\frac{\left|t^{\prime }-t\right|+\left(t-1\right)H}{V_l}+2\left(2{\gamma }_l+{\epsilon }_l\right)$$

(12)

2.

Operating service time of horizontal operations

Similar to the case described in Section 4.2.1, the cross-aisle probability is:

$$P_1\left(x\right)=P\left(x\right)$$

(13)

Therefore, the service time of the shuttle car to complete the horizontal operation is:

$$T_{s2}\left(M_{ac}\right)=T_s\left(M_{ac}\right)$$

(14)

3.

The operating service time of the remaining operations is:

$$T_{ls}\left(t^{\prime }\right)=T_{ls}\left(t\right)$$

(15)

Due to assumptions 4 and 6, the shuttle and hoist are not released until they complete a job; thus, no parallelism occurs in this case.

3.3. Modelling the Outbound Process of a Four-Way Shuttle System Based on Two Operation Strategies

The parallel and serial queuing network model for the four-way shuttle system outbound process are shown in Figure 2.

3.3.1. Parallel Part Modelling

Construct the separation–aggregation queuing network model and analyse the parallel operation strategy by decomposing the complex network method in the following steps:

• Constructing the separation–aggregation queuing network model

In this section, we have shown that the whole parallel operation strategy is divided into vertical operations, horizontal operations and remaining operations. First, the pickup order arrives at the “separation” point at arrival rate ${\lambda }_r$, at which time the pickup task is divided into a vertical operation dominated by a hoist and a horizontal operation dominated by a shuttle, where the T shelves are divided into T different types of customers, arriving at arrival rate ${\lambda }_r^1,{\lambda }_r^2,\dots ,{\lambda }_r^T$. Next, the hoist and the shuttle wait in the vertical queue $Q_l$and the horizontal queue $Q_s$, respectively. After both the vertical and horizontal operations are completed, they converge at the “convergence” point; finally, the hoist carries the goods to complete the last remaining operations. As assumption 6 obeys the single-command cycle, the lifters are not released until the task is completed; therefore, the lifters under node ${\mu }_{l2}$ can provide an unlimited service capacity. The separation aggregation queuing network model is shown in Figure 3.

2.

Constructing a closed-loop queueing network

As the shuttle and hoist obey a general distribution, the parallel operation strategy described in Section 3.3.1 is a nonproduct solution and cannot be solved directly; thus, the Coxian distribution solution model, which approximates the general distribution, is applied in this paper. Figure 4 shows a schematic diagram of a closed-loop queuing network under a parallel strategy, where there are ${\mathit{K}}_{\mathit{r}}$ pick-up operations in the closed-loop network and the arrival rate of each tier of jobs is ${\mathit{\lambda }}_{\mathit{r}}^{\mathit{i}}=\frac{{\mathit{\lambda }}_{\mathit{r}}}{\mathit{T}}$. After entering this network, the shuttles are already bound to the order, the maximum number of accesses that can be accepted in the system cannot exceed the number of shuttles and the time used in the jobs to receive service at the shuttles and hoists is equal to their own service time.

The importance of building this closed-loop queuing network model is to create an equivalent service node to the open-loop queuing network model, which can simplify the complex queuing network.

3.

Constructing an open-loop queuing network

First, the service point formed by the closed-loop queuing network is denoted as node 1, and the service rate of this service point obeys a negative exponential distribution, which leads to an open-loop queuing network with two servers in series, where the service rate of service point 1 is replaced by the throughput rate sought in 2, which is load dependent. In the open-loop queuing network, the pickup task arrives with an arrival rate of ${\mathit{\lambda }}_{\mathit{r}}$, node 1 enters node 2 for the remaining operations after equivalence and as node 1 is composed of a closed-loop queuing network, it compensates for the lack of trolley waiting time for the orders in the open-loop network. Finally, the maximum entropy method is used to solve the metrics of this queuing network model. The open-loop queuing network formed under the parallel strategy is shown in Figure 5.

3.3.2. Serial Part Modelling

If there are no idle shuttles at the target level, the serial operation strategy will be employed, constituting a semi-open-loop queuing network, as shown in Figure 6. The “synchronisation node” can be understood as the system receiving the pickup command. The pickup command, after entering the system and the shuttle for binding, is entered as a task in the queue of the waiting hoist node S1; if there is an idle hoist, it can complete the task, whereby a hoist picks up an idle shuttle. After completing the task at shuttle node S2, which is waiting to pick up the goods, the hoist picks up the shuttle and moves it to the target layer to pick up the goods; throughout the whole process, the shuttle is still in the task and is not released; therefore, the shuttle node is able to directly carry out the task of picking up the goods until the hoist receives the goods once they have been released. Finally, we consider hoist node S3, where the remaining operations are performed. Similar to node S2, the hoist is not released during the task, and therefore, any remaining operations do not need to wait.

Similar to the parallel operation strategy, the Coxian distribution is used to approximate the general distribution of the service times instead. In the process of solving the semi-open-loop queuing network model, the use of more than two service nodes will cause a sharp increase in the state space, meaning that the model cannot be computed. Therefore, two nodes need to be aggregated into one, and as hoist node S3 does not need to be queued for direct service, hoist nodes S1 and S3 are combined into one, as shown in Figure 7.

The whole system model is composed based on the probabilities of the two cases under the completed parallel operation strategy and the serial operation strategy.

4. Model Solution

4.1. Solving the Separation-Aggregation Queueing Network Model Using the Parallel Strategy

As the operating times for the horizontal operations, vertical operations and remaining operations obey the general distribution, and the Coxian distribution is used as an approximation of this in this paper, the mean value of the horizontal operating time can be derived from Equations (2) and (7):

$${\tau }_s={\displaystyle \sum }_{a=1}^N{\displaystyle \sum }_{c=1}^{C}P\left(M_{ac}\right)T_s\left(M_{ac}\right)$$

(16)

where $\mathit{a}=1,2\dots \mathrm{N}$, $\mathit{c}=1,2,\dots ,\mathrm{C}$ and the coefficient of the variation of the square of the horizontal operating time is$\mathrm{c}{v}_s^2$. Similarly, the mean values of the vertical operation and remaining operation time ${\tau }_{l1}$,${\tau }_{l2}$, can be expressed according to Equations (1), (4) and (8), as follows:

$${\tau }_{l1}={\displaystyle \sum }_{t=1}^{T}P\left(t\right)T_{l1}\left(t\right)$$

(17)

$${\tau }_{l2}={\displaystyle \sum }_{t=1}^{T}P\left(t\right)T_{ls}\left(t\right)$$

(18)

Additionally, set the squared coefficients of the variation for the vertical operation and the remaining operation times are set as $\mathrm{c}{v}_{l1}^2$, $\mathrm{c}{v}_{l2}^2$.

Thus, the service rates of the horizontal operation of the shuttle, the vertical operation of the hoist and the remaining operation of the hoist, ${\mu }_s$,${\mu }_{l1}$ and ${\mu }_{l2}$can be expressed by their mean values, as follows:

$${\mu }_s=\frac{1}{{\tau }_s},t=1,2,\dots ,T$$

(19a)

$${\mu }_{l1}=\frac{1}{{\tau }_{l1}}$$

(19b)

$${\mu }_{l2}=\frac{1}{{\tau }_{l2}}$$

(19c)

4.1.1. Constructing a Closed-Loop Queuing Network to Compute State-Based Service Rates

A closed-loop queuing network model based on separation-aggregation is shown in Figure 2, where there are $K_r$ pickups in the closed-loop network, $Q_l$is the number of customers waiting to be served in $Q_l^q$ and $Q_s$D is the number of customers waiting to be served in $Q_s^q$. The state variable $s_k$of this closed-loop queueing network can be given as:

$$s_k=\left(Q_s,Q_L\right),k=Q_l+Q_s+m{Q}_l=Q_l\left(m+1\right)+Q_s$$

(20)

where $\mathrm{m}=\min (K_r,$ M) is the maximum number of customers waiting to be served in $Q_s^q$ and$Q_l+Q_s<m+1,Q_l=0,1,Q_s=0,1,\dots ,m$.

If the system completes an aggregation using the parallel operation strategy, there are two scenarios:

1.

At $Q_l=1$, the hoist completes its operation and reaches the convergence point. This process requires the shuttle to complete the pickup operation and move to the entrance of the aisle; then, the service rate of this process is ${\mu }_s$.

2.

At $Q_s>1$, the shuttle completes the pickup operation and arrives at the hoistway entrance, and this process requires the hoistway to move to the entrance of the target level. Then, the service rate of this process is ${\mu }_{l1}$.

Setting $\pi \left(s_k\right)$ as the state probability in state $s_k$, the throughput rate can be expressed by the two cases of aggregation as:

$${\lambda }_f\left(K_r\right)={\displaystyle \sum }_{Q_s=0}^m{\mu }_{l1}\pi \left(Q_s,0\right)+{\displaystyle \sum }_{Q_s=0}^{m-1}{\mu }_s\pi \left(Q_s,1\right)$$

(21)

The service time of shuttle service point ${\mu }_s$is independently and identically distributed, denoted as $\left(\mathsf{\alpha },\mathrm{T}\right)$, where the initial probability vector$\mathsf{\alpha }={\left[\begin{array}{ccccc}1& 0& 0& 0& \dots \end{array}\right]}_{1\times p}$ is the initial probability vector of the p-order Coxian distribution and T is the p-order square matrix. The service time of hoist service point ${\mu }_{l1}$is independently and identically distributed, denoted as $\left(\beta ,S\right)$, and the initial probability vector$\beta ={\left[\begin{array}{ccccc}1& 0& 0& 0& \dots \end{array}\right]}_{1\times p}$ is the initial probability vector of the q-order Coxian distribution and S is the q-order square matrix.

The state variables of the separation-aggregation closed-loop queuing network model are set to $\left(s_k,i,j\right)$, where $i,j$ are the states in which the shuttle and hoist are located, respectively, where $i=1,2,3\dots p;andj=1,2,3\dots q$. Let ${\pi }_t$ be the steady-state probability distribution of $s_k$. As the two cases complete the parallel operation strategy complete the aggregation, ${\pi }_t=\left[{\pi }_0,{\pi }_1\right]$.

Therefore, the state transfer rate matrix Q of the state variable $s_k$ is:

$$Q=\left[\begin{array}{cc}{B}_{00}& B_{01}\\ B_{10}& B\end{array}\right]$$

(22)

The transfer diagram of the generation and extinction process is shown in Figure 8.

$B_{01}$ represents the rate matrix for the transfer of the ${\pi }_0=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state to the ${\pi }_1=[\pi \left(s_{m+1)},\pi \left(s_{m+2}\right),\pi \left(s_{m+1}\right),\pi \left(s_{m+3}\right),\dots ,\pi \left(s_{2m}\right)\right]$ state.

$$B_{01}=\left[\begin{array}{cccccc}\left(Q_s,Q_l\right)& \left(0,1\right)& \left(1,1\right)& \left(2,1\right)& \dots & \left(m-1,1\right)\\ \left(0,0\right)& I_T\otimes S^{\circ }& & & & \\ \left(1,0\right)& & I_T\otimes S^{\circ }& & & \\ \left(2,0\right)& & & I_T\otimes S^{\circ }& & \\ ⋮& & & & \ddots & \\ \left(m-1,0\right)& & & & & I_T\otimes S^{\circ }\\ \left(m,0\right)& & & & & 0\end{array}\right]$$

(23)

$B_{00}$ represents the rate matrix for the transfer of the ${\pi }_0=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state to the ${\pi }_0=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state.

$$B_{00}=\left[\begin{array}{cccccccc}\left(Q_s,Q_l\right)& \left(0,0\right)& \left(1,0\right)& \left(2,0\right)& \left(2,1\right)& \dots & \left(m-1,0\right)& \left(m,0\right)\\ \left(0,0\right)& T\oplus S& T\alpha \otimes I_S& & & & & \\ \left(1,0\right)& I_T\otimes S\beta & T\oplus S-I_T\otimes {\tilde{S}}^{\circ }& T^{\circ }\alpha \otimes I_S& & & & \\ \left(2,0\right)& & I_T\otimes S\beta & T\oplus S-I_T\otimes {\tilde{S}}^{\circ }& T^{\circ }\alpha \otimes I_S& & & \\ ⋮& & & \ddots & & & & \\ \left(m-1,0\right)& & & & I_T\otimes S\beta & & T\oplus S-I_T\otimes {\tilde{S}}^{\circ }& T^{\circ }\otimes I_S\\ \left(m,0\right)& & & & & & S^{\circ }\beta \otimes \alpha & S\end{array}\right]$$

(24)

where $I_T$ and $I_S$ are unit matrices that have the same order as T and S, respectively, $\otimes$ is the Kronecker product, $\oplus$is the Kronecker sum and$T\oplus S=T\otimes I_S+I_T\otimes S$. In the state transfer rate matrix, the sum of all the elements in each row is 0, so that ${\tilde{S}}^{\circ }$ is found to be:

$${\tilde{S}}^{\circ }=\left[\begin{array}{c}{\mu }_2^{\beta }\left(1-b_1^{\beta }\right)\\ {\mu }_1^{\beta }\left(1-b_2^{\beta }\right)\\ ⋮\\ {\mu }_{q-1}^{\beta }\left(1-b_{q-1}^{\beta }\right)\\ {\mu }_q^{\beta }\end{array}\right]$$

(25)

$B_{10}$ is the rate matrix of the transfer from the ${\pi }_1=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state to the ${\pi }_1=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state.

$$B_{01}=\left[\begin{array}{ccccccc}\left(Q_s,Q_l\right)& \left(0,0\right)& \left(1,1\right)& \left(2,1\right)& \dots & \left(m-1,0\right)& \left(m,0\right)\\ \left(0,1\right)& T^{\circ }\alpha & & & & & \\ \left(1,1\right)& & T^{\circ }\alpha \otimes \beta & & & & \\ \left(2,1\right)& & & T^{\circ }\alpha \otimes \beta & & & \\ ⋮& & & & \ddots & & \\ \left(m-1,1\right)& & & & & T^{\circ }\alpha \otimes \beta & 0\end{array}\right]$$

(26)

B is the rate matrix of the transfer from the ${\pi }_1=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state to the ${\pi }_1=\left[\pi \left(s_0\right),\pi \left(s_1\right),\pi \left(s_2\right)\dots \pi \left(s_m\right)\right]$ state.

$$B=\left[\begin{array}{cccccc}\left(Q_s,Q_l\right)& \left(0,1\right)& \left(1,1\right)& \left(2,1\right)& \dots & \left(m-1,1\right)\\ \left(0,1\right)& T-{\tilde{T}}^{\circ }& T^{\circ }\alpha & & & \\ \left(1,1\right)& & T-{\tilde{T}}^{\circ }& T^{\circ }\alpha & & \\ \left(2,1\right)& & & T-{\tilde{T}}^{\circ }& & \\ ⋮& & & & \ddots & T^{\circ }\alpha \\ \left(m-1,1\right)& & & & & T-{\tilde{T}}^{\circ }\end{array}\right]$$

(27)

The ${\tilde{T}}^{\circ }$ is as follows:

$${\tilde{T}}^{\circ }=\left[\begin{array}{c}{\mu }_1^{\alpha }\left(1-b_1^{\alpha }\right)\\ {\mu }_2^{\alpha }\left(1-b_2^{\alpha }\right)\\ ⋮\\ {\mu }_{P-1}^{\beta }\left(1-b_{p-1}^{\beta }\right)\\ {\mu }_p^{\alpha }\end{array}\right]$$

(28)

3.

According to the smooth equation $\{\begin{array}{c}\overrightarrow{\pi }\cdot Q=\overrightarrow{0}\\ \overrightarrow{\pi }\cdot e=1\end{array}$, the equilibrium equation of the model is obtained as follows:

$$\begin{gathered} \left[\begin{array}{cc}{\pi }_0& {\pi }_1\end{array}\right]\cdot Q=\overrightarrow{0} \\ \left[\begin{array}{cc}{\pi }_0& {\pi }_1\end{array}\right]\cdot e=1 \end{gathered}$$

(29)

4.1.2. An Open-Loop Queueing Network Is Constructed to Calculate the Separation-Aggregation Queueing Network Model Metrics

After the calculations in the previous subsection, the system can reduce the subnetwork in the separation-aggregation queuing network to an equivalent simultaneous flow of servers, with the hoist forming a serial open-loop queuing network, which assumes a negative exponential distribution for node 1, whose service rate is determined by the throughput rate and the load in the system; thus, it is set to ${\mu }_{FES}\left(x\right)={\lambda }_f\left(x\right)$.

Setting the state probability in the open-loop queuing network as $\pi \left(k_1,k_2\right)$, $k_{1}and{k}_2$ as the number of customers served by node 1 and node 2 and$\pi \left(k_1\right)$,$\pi \left(k_2\right)$ as the approximate probabilities of node 1 and node 2, respectively. The state probability of this network can be given by the number of customers served by each node.

$$\pi \left(k_1,k_2\right)={\pi }_1\left(k_1\right)·{\pi }_2\left(k_2\right)$$

(30)

(1)

For service point 1, the average service rate and the average number of customers served are${\pi }_1\left(k_1\right)$; this node can be found on the basis of the approximate probability ${\overline{K}}_1$. The equilibrium equation is shown in

Table 3. Equilibrium equation transfer table.

Number of Customers (Status)	Equilibrium Equation	Simplified Equations
0	${\pi }_1\left(0\right){\lambda }_r={\pi }_1\left(1\right){\mu }_{FES}\left(1\right)$	${\pi }_1\left(1\right)={\pi }_1\left(0\right)\frac{{\lambda }_r}{{\mu }_{FES}\left(1\right)}$
1	${\pi }_1\left(1\right)\left({\mu }_{FES}\left(1\right)+{\lambda }_r\right)={\pi }_1\left(0\right){\lambda }_r+{\pi }_1\left(2\right){\mu }_{FES}\left(2\right)$	${\pi }_1\left(2\right)={\pi }_1\left(0\right)\frac{{\lambda }_r^2}{{\mu }_{FES}\left(1\right)·{\mu }_{FES}\left(2\right)}$
2	${\pi }_1\left(2\right)\left({\mu }_{FES}\left(2\right)+{\lambda }_r\right)={\pi }_1\left(1\right){\lambda }_r+{\pi }_1\left(3\right){\mu }_{FES}\left(3\right)$	${\pi }_1\left(3\right)={\pi }_1\left(0\right)\frac{{\lambda }_r^2}{{\mu }_{FES}\left(1\right)·{\mu }_{FES}\left(2\right)·}{\mu }_{FES}\left(1\right)$
…	…	…

, and the birth and death process is shown in Figure 9.

The approximate probability of finding node 1 is:

$${\pi }_1\left(0\right)=\frac{1}{1+{{\displaystyle \sum }}_{k_1}^{+\infty }{{\displaystyle \prod }}_{i=1}^{k_1}\frac{{\lambda }_r}{{\mu }_{FES}\left(i\right)}}$$

(31a)

$${\pi }_1\left(k_1\right)={\pi }_1\left(0\right){\displaystyle \prod }_{i=1}^{k_1}\frac{{\lambda }_r}{{\mu }_{FES}\left(i\right)}$$

(31b)

As the service rate of node 1 is equivalently represented by the throughput rate calculated in Section 4.1.1 and the throughput rate is related to the load of node 1, the average service rate of node 1 is:

$${\mu }_1={\displaystyle \sum }_{k_1}^{+\infty }\pi \left(k_1\right)\cdot {\mu }_{FES}\left(k_1\right)$$

(32)

The average number of customers served at node 1, ${\overline{K}}_1$, can be calculated from the approximate probability:

$${\overline{K}}_1={\displaystyle \sum }_{k_1}^{+\infty }{k}_1·{\pi }_1\left(k_1\right)$$

(33)

(2)

Solving node 2

As node 1 obeys the negative exponential distribution and node 2 obeys the general distribution, the queueing network formed is a nonproduct solution and the queueing network formed is not easy to solve; therefore, in this paper, the MEM is applied to the solution of the open-loop queueing network, where the entropy is the measure of uncertainty that can be predicted for something. For node 2, the entropy function is introduced according to the maximum entropy principle in the form of a universal ME solution ${\pi }_2\left(n\right)$ by maximising the entropy generalisation function H(p) under the following constraints [34,35]:

$$\mathrm{H}\left(\pi \right)=-{\displaystyle \sum }_{n=1}^{k_2}{\pi }_2\left(n\right)\ln \left[{\pi }_2\left(n\right)\right],k_2\ge 1$$

(34)

The optimisation model of the maximum is:

$$\max H\left(\pi \right)=-{\displaystyle \sum }_{n=1}^{k_2}{\pi }_2\left(n\right)\ln \left[{\pi }_2\left(n\right)\right]$$

(35)

$$\mathrm{s}.\mathrm{t}.\{\begin{array}{c}{\displaystyle {\displaystyle \sum }_{n=0}^{k_2}}p\left(n\right)=1\\ {\displaystyle {\displaystyle \sum }_{n=1}^{k_2}}h\left(n\right){\pi }_2\left(n\right)=\rho \\ {\displaystyle {\displaystyle \sum }_{k=1}^{k_2}}n{\pi }_2\left(n\right)=L\end{array}$$

(36)

The approximate probability of node 2 is:

$${\pi }_2\left(n\right)=\{\begin{array}{cc}\left(1-\rho \right)& n=0\\ \left(1-\rho \right)g{x}^n& n\ge 1\end{array}$$

(37)

where $\mathrm{g}=\frac{{\rho }_2^2}{\left(1-{\rho }_2\right)\left(L-{\rho }_2\right)},x=\frac{L-{\rho }_2}{L}and{\rho }_2=\frac{{\lambda }_r}{{\mu }_{l2}}$, applying the generalised Laplace equation and the z-transform equation [36,37], showing that for any queue, the average queue length is given by the same equation:

$$\text{L}=\frac{{\rho }_2}{2}\left(1+\frac{C_{a_2}^2+{\rho }_2{C}_s^2}{1-{\rho }_2}\right)$$

(38)

where $C_{a_2}^2$ is the squared coefficient of the variation of the customer arrival time at node 2 and $C_s^2$ is the squared coefficient of the variation of the service time at node 2 [38].

Setting the average number of customers at node 2 to ${\overline{K}}_2$, the approximate probability can be derived from:

$${\overline{K}}_2={\displaystyle \sum }_{k_2=1}^{+\infty }{k}_2\cdot {\pi }_2\left(k_2\right)$$

(39)

(3)

Solving for system performance metrics

1. The system response time, which is the time required to complete a pickup operation, including the waiting time and the service time, can be derived using Little’s theorem.

$$E\left[T\right]=\frac{{\overline{K}}_1}{{\lambda }_r}+\frac{1}{{\mu }_{l2}}$$

(40)

2. The waiting time in the external queue for pickup is:

$$W_r=\frac{{\overline{K}}_1}{{\lambda }_r}-\frac{1}{{\overline{\mu }}_1}$$

(41)

3. The external queue waiting for pickup is:

$$L_{eq}=W_r\cdot {\lambda }_r$$

(42)

4. As this paper studies multiple shuttles in a cell, the arrival rate of each shuttle is $\frac{{\lambda }_r}{M}$ and the utilisation rate of the shuttles can be calculated as:

$${\rho }_s=\frac{{\lambda }_r}{M\cdot {\mu }_s}$$

(43)

4.2. Solving the Semi-Open-Loop Queuing Network Model under a Serial Strategy

4.2.1. Solving the Transfer Rate Matrix Q

Define the system state as $\left(n_1,n_2,i,j,k\right)$, where $i\ge 0$is the number of external queue pickup orders, $0\le j\le M$ is the number of hoist node pickup orders and $0\le k\le M$ is the number of shuttle node pickup orders. Set $n_1=i+j$ as the sum of the number of the external queues and first-node pickup orders for the sake of the solution convenience of solving and combine it with$n_2=k$ to form the birth and death process of the two states.

If the sum of the external queue, the first node and the second node pickup order is greater than M, then there must be idle four-way shuttles in the system, and the external queue waiting for pickup orders will be 0. If the sum of the external queue, the first node and the second node pickup order number is more than M, then there are no idle shuttles in the system, and the whole system is busy, so the external queue is defined as [39]:

$$\mathrm{i}=\{\begin{array}{cc}0& n_1+n_2\le M\\ n_1+n_2-M& n_1+n_2>M\end{array}$$

(44)

When the performance metrics are calculated (e.g., average queue length), only $n_1$and$n_2$ are relevant. Therefore, for the sake of brevity, this semi-open-loop intraqueueing network state variable is modified to $\left(n_1,n_2\right)$. Defining the steady-state probability vector as ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)$,${\pi }_0=\left({\pi }_{\left(i,0\right)},{\pi }_{\left(i,1\right)},{\pi }_{\left(i,2\right)}\dots \right){\pi }_i=\left({\pi }_{\left(i,0\right)},{\pi }_{\left(i,1\right)},{\pi }_{\left(i,2\right)}\dots \right)$, on the basis of the Coxian distribution, where $\left(\eta ,H\right)$ denotes the arrival interval time distribution, $\left(\omega ,W\right)$ denotes the service time of node 1 and $\left(\gamma .,U\right)$ denotes the service time of node 2, the transfer rate matrix can be expressed as:

$$Q=\left[\begin{array}{ccccc}{B}_0& C_0& & & \\ A_1& B& C& & \\ & A& \mathrm{B}& C& \\ & \ddots & \ddots & \ddots & \end{array}\right]$$

(45)

where $B_0$is the rate matrix for the transfer from the ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)$ state to the ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)$ state, which is an $\left(M+1\right)\times \left(M+1\right)$-squared matrix, M is the number of shuttles served by the hoist, $I^H$is a unit matrix of the same order as H, $U^{\circ }$ is the absorption probability matrix of U, $\otimes$is the Kronecker product and$\oplus$is the Kronecker sum.

A transfer diagram of the generation and extinction process is shown in Figure 10.

$$B_0=\left[\begin{array}{cccccc}\left(n_1,n_2\right)& \left(0,0\right)& \left(0,1\right)& \left(0,2\right)& \dots & \left(0,M\right)\\ \left(0,0\right)& H& & & & \\ \left(0,1\right)& I^H\otimes U^{\circ }& H\oplus U& & & \\ \left(0,2\right)& & I^H\otimes U^{\circ }\gamma & H\oplus U& & \\ ⋮& & \ddots & & \ddots & \\ \left(0,M\right)& & & & I^H\otimes U^{\circ }\gamma & H\oplus U\end{array}\right]$$

(46)

$C_0$ is the rate matrix of the transfer from the ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)\mathrm{state}$ tothe ${\pi }_1=\left({\pi }_{\left(1,0\right)},{\pi }_{\left(1,1\right)},{\pi }_{\left(1,2\right)}\dots \right)$ state, $I^U$ is the unit matrix of the same order as U, and $H^{\circ }$ is the absorption probability matrix of H.

$$C_0=\left[\begin{array}{cccccc}\left(n_1,n_2\right)& \left(1,0\right)& \left(1,1\right)& \left(1,2\right)& \dots & \left(1,M\right)\\ \left(0,0\right)& H^{\circ }\eta \otimes \omega & & & & \\ \left(0,1\right)& & H^{\circ }\eta \otimes \omega \otimes I^U& & & \\ \left(0,2\right)& & & H^{\circ }\eta \otimes \omega \otimes I^U& & \\ ⋮& & & & \ddots & \\ \left(0,M\right)& & & & & H^{\circ }\eta \otimes I^U\end{array}\right]$$

(47)

$A_1$is the rate matrix of the transfer from that state to the ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)$ state and $W^{\circ }$ is the absorption probability matrix of the same order as W. As the upper limit of the capacity of node 2 is M, the transfer rate is 0 in the transfer from the state (1, M) to (0, M).

$$A_1=\left[\begin{array}{cccccc}\left(n_1,n_2\right)& \left(0,0\right)& \left(0,1\right)& \left(0,2\right)& \dots & \left(0,M\right)\\ \left(1,0\right)& & I^H\otimes W^{\circ }\otimes \gamma & & & \\ \left(1,1\right)& & & I^H\otimes W^{\circ }\otimes I^U& & \\ ⋮& & & & \ddots & \\ \left(1,M-1\right)& & & & & I^H\otimes W^{\circ }\otimes I^U\\ \left(1,M\right)& & & & & 0\end{array}\right]$$

(48)

B is the rate matrix of the ${\pi }_1=\left({\pi }_{\left(1,0\right)},{\pi }_{\left(1,1\right)},{\pi }_{\left(1,2\right)}\dots \right)$ to ${\pi }_1=\left({\pi }_{\left(1,0\right)},{\pi }_{\left(1,1\right)},{\pi }_{\left(1,2\right)}\dots \right)$ state transfer, and $I^W$ represents the unit matrix of the same order as W.

$$B=\left[\begin{array}{cccccc}\left(n_1,n_2\right)& \left(1,0\right)& \left(1,1\right)& \dots & \left(1,M-1\right)& \left(1,M\right)\\ \left(1,0\right)& H\oplus W& & & & \\ \left(1,1\right)& I^H\otimes I^W\otimes U^{\circ }& H\oplus W\otimes U& & & \\ ⋮& & & \ddots & & \\ \left(1,M-1\right)& & & I^H\otimes I^W\otimes U^{\circ }\gamma & H\oplus W\oplus U& \\ \left(1,M\right)& & & & I^H\otimes \omega \otimes U^{\circ }\gamma & H\oplus U\end{array}\right]$$

(49)

C is the rate matrix for the transfer from the ${\pi }_1=\left({\pi }_{\left(1,0\right)},{\pi }_{\left(1,1\right)},{\pi }_{\left(1,2\right)}\dots \right)$ to the ${\pi }_2=\left({\pi }_{\left(2,0\right)},{\pi }_{\left(2,1\right)},{\pi }_{\left(2,2\right)}\dots \right)$ state and A is the rate matrix for the transfer from the ${\pi }_2=\left({\pi }_{\left(2,0\right)},{\pi }_{\left(2,1\right)},{\pi }_{\left(2,2\right)}\dots \right)$ to the ${\pi }_1=\left({\pi }_{\left(1,0\right)},{\pi }_{\left(1,1\right)},{\pi }_{\left(1,2\right)}\dots \right)$ state.

$$C=\left[\begin{array}{ccccccc}\left(n_1,n_2\right)& \left(1,0\right)& \left(2,1\right)& \dots & & \left(2,M-1\right)& \left(2,M\right)\\ \left(1,0\right)& H^{\circ }\eta \otimes I^W& & & & & \\ \left(1,1\right)& & H^{\circ }\eta \otimes I^W\otimes I^U& & & & \\ ⋮& & & & \ddots & & \\ \left(1,M-1\right)& & & & & H^{\circ }\eta \otimes I^W\otimes I^U& \\ \left(1,M\right)& & & & & & H^{\circ }\eta \otimes I^U\end{array}\right]$$

(50)

$$A=\left[\begin{array}{ccccccc}\left(n_1,n_2\right)& \left(1,0\right)& \left(1,1\right)& \left(1,2\right)& & \dots & \left(2,M\right)\\ \left(2,0\right)& & I^H\otimes W^{\circ }\otimes \gamma & & & & \\ \left(2,1\right)& & & I^H\otimes W^{\circ }\otimes I^U& & & \\ ⋮& & & & \ddots & & \\ \left(2,M-1\right)& & & & & & I^H\otimes W^{\circ }\otimes I^U\\ \left(2,M\right)& & & & & & 0\end{array}\right]$$

(51)

4.2.2. Solving for the Steady-State Probability Distribution

If a Markov state exists as a steady-state, then $\pi Q=0$, using the steady-state probability vector ${\pi }_0=\left({\pi }_{\left(0,0\right)},{\pi }_{\left(0,1\right)},{\pi }_{\left(0,2\right)}\dots \right)$,${\pi }_i=\left({\pi }_{\left(i,0\right)},{\pi }_{\left(i,1\right)},{\pi }_{\left(i,2\right)}\dots \right)$, rewriting πQ = 0 as a matrix equation yields [39,40]:

$$\begin{gathered} {\pi }_{0}B+{\pi }_1{A}_1=0 \\ {\pi }_0{C}_0+{\pi }_{1}B+{\pi }_{2}A=0 \\ {\pi }_{1}C+{\pi }_{2}B+{\pi }_{3}A=0 \\ ⋮ \\ {\pi }_{i-1}C+{\pi }_{i}B+{\pi }_{i+1}A=0 \end{gathered}$$

(52)

where ${\pi }_0,{\pi }_i$is a row vector of 1 × (M + 1), ${\pi }_i$ is denoted as the steady-state probability vector of $i$in each state and 0 is a row vector with element 0. The class of the birth and death state can be classed as a constant return state, i.e., it returns to itself with a probability of 1 in a finite number of steps, so the relationship between the steady-state probabilities can be written as:

$${\pi }_i={\pi }_0{R}^i$$

(53a)

$${\pi }_{i+1}={\pi }_{i}R$$

(53b)

where R is the probability matrix as a square matrix of (M + 1) × (M + 1); therefore, the formula can be simplified as:

$$C+RB+R^{2}A=0\Rightarrow R=-\left(C+R^{2}A\right)B^{-1}$$

(54)

4.2.3. Solving the Performance Index

After solving the steady-state probability vectors, the captain of each node can be derived from the expected value of each state. Determine the state space of shuttle node S₂ as $S_2=\left[0,1,2,3,\dots ,M,0,1,2,3,\dots .M,0,1,\dots \right]$; then, the captain $L_{S_2}$ of shuttle node S₂ is:

$$L_{S_2}={\pi }_0\left[\begin{array}{c}0\\ 1\\ 2\\ ⋮\\ M\end{array}\right]+{\pi }_1\left[\begin{array}{c}0\\ 1\\ 2\\ ⋮\\ M\end{array}\right]+\dots +{\pi }_i\left[\begin{array}{c}0\\ 1\\ 2\\ ⋮\\ M\end{array}\right]+\dots$$

(55)

Then, the shuttle utilisation rate ${\rho }_s$ is:

$${\rho }_s=\frac{L_{S_2}}{M}$$

(56)

The total captain $L_S$ of the external queue and elevator node S₁ is:

$$L_S={\displaystyle \sum }_{i=1}^{\infty }{\pi }_1{R}^{i-1}ie={\pi }_1{\left(I-R\right)}^{-1}e$$

(57)

The external queue captain $L_{eq}$ is [41]:

$$L_{eq}={\displaystyle \sum }_{i=1}^{\infty }\left|{\pi }_i\right|$$

(58)

where $\left|{\pi }_i\right|$ is the sum of the elements in ${\pi }_i$ and the external queue waiting time $W_r$ is:

$$W_r=\frac{L_{eq}}{{\lambda }_r}$$

(59)

Then, the captain $L_{S1}$ of the hoist node is:

$$L_{S1}=L_S-L_{eq}$$

(60)

The system response time $E\left(T\right)$ for the nodal captain of the hoist and the shuttle as is shown in Equation (61).

$$E\left(T\right)=\frac{L_{S1}+L_{S2}}{{\lambda }_r}$$

(61)

5. Simulation Test and Data Analysis

5.1. Arena Simulation Validation

As the scale of the four-way shuttle system varies according to the size of the business volume, this paper lists seven different four-way shuttle system scales with different numbers of layers and columns and sets five different arrival rates. The other assumptions of the system are consistent with the model established in this paper, i.e., only pickup operations are considered; the first come, first served principle is adopted; and the random storage strategy is adopted to verify the reliability and accuracy of the model at different scales and arrival rates. The reliability and accuracy of the model are verified for different sizes and arrival rates. The initial system parameters set in the simulation are shown in

Table 4. Simulation parameter values.

Parameters	Take Value	Parameters	Take Value
Vs (m/s)	2	${\gamma }_l$ (s)	2
${\epsilon }_s$ (s)	3	L (m)	1.2
(s)	2	W (m)	1
Vl (m/s)	3	WA (m)	1
${\epsilon }_l$ (s)	1	H (m)	1.5
N	3	M	7

, the system size is shown in

Table 5. Simulation scale.

Size	Small Scale				Large Scale
Serial number	1	2	3	4	5	6	7
layers	5	6	7	8	9	10	11
columns	36	40	44	48	64	72	80
Capacities	360	480	616	768	1152	1440	1760

and the arrival rate is set to λ_r = 50, 100, 150, 200, 250 per hour.

Before commencing the seven simulation scenarios, a warm-up period of 100 h was implemented to eliminate errors and stabilise the system, followed by a confidence interval of 95% and a simulation half-step of less than 2% of the mean value, which was run ten times, averaging over 2000 h each time, in order to obtain the final results. The metrics calculated for the simulation are the system response time, the external queue waiting for pickup and the shuttle utilisation. The error value between the simulation results and the model solution results was calculated using Equation (62), where $R_1$ represents the simulation results and $R_2$ represents the model solution results:

$$\mathsf{\xi }=\frac{\left|R_1-R_2\right|}{R_1}\times 100\%$$

(62)

From the results in

Table 6. Comparison of simulation results and theoretical values.

Number	1			2			3			4			5			6
${\lambda }_r=50$	$R_1$	$R_2$	$\zeta$	R1	$R_2$	ζ	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$
$E\left(T\right)$	86.26	88.46	2.55	94.11	98.25	4.40	106.07	107.69	1.53	110.44	113.90	3.14	133.58	139.57	4.48	148.24	153.02	3.22
$W_r$	37.03	37.40	0.99	62.42	60.18	3.59	67.98	64.70	4.82	69.71	66.25	4.96	90.12	87.11	3.34	100.01	95.91	4.10
$L_{eq}$	0.51	0.52	0.99	0.87	0.84	3.59	0.94	0.90	4.82	0.97	0.92	4.96	1.25	1.21	3.34	1.39	1.33	4.10
${\rho }_s$	0.15	0.15	0.21	0.16	0.16	1.00	0.14	0.15	5.82	0.14	0.13	3.53	0.15	0.15	0.13	0.15	0.15	0.39
${\lambda }_r=100$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$
$E\left(T\right)$	95.27	96.78	1.58	102.78	106.94	4.05	110.25	114.73	4.06	114.18	110.59	3.15	140.14	145.26	3.65	154.77	166.82	7.79
$W_r$	64.13	63.71	0.65	69.01	68.85	0.23	74.18	73.76	0.56	76.43	77.56	1.48	97.24	97.78	0.56	108.81	109.68	0.80
$L_{eq}$	1.78	1.77	0.65	1.92	1.91	0.23	2.06	2.05	0.56	2.12	2.15	1.48	2.70	2.72	0.56	3.02	3.05	0.80
${\rho }_s$	0.34	0.35	2.94	0.32	0.32	0.94	0.29	0.30	4.21	0.28	0.27	2.07	0.32	0.30	5.17	0.30	0.30	1.53
${\lambda }_r=150$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	ζ	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	ζ
$E\left(T\right)$	110.37	105.89	4.06	117.37	116.50	0.74	123.70	124.53	0.67	126.45	124.69	1.39	153.86	163.77	6.44	174.07	182.07	4.60
$W_r$	76.89	78.77	2.45	80.94	85.56	5.71	84.63	89.92	6.25	85.75	88.43	3.13	111.26	120.10	7.95	123.62	130.86	5.86
$L_{eq}$	3.20	3.28	0.10	3.37	3.57	5.71	3.53	3.75	0.26	3.57	3.68	0.13	4.64	5.00	7.95	5.15	5.45	5.86
${\rho }_s$	0.51	0.53	2.35	0.48	0.48	0.39	0.45	0.45	0.69	0.41	0.42	1.46	0.45	0.46	0.44	0.45	0.46	0.88
${\lambda }_r=200$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	R1	$R_2$	$\zeta$	R1	$R_2$	$\zeta$
$E\left(T\right)$	124.14	119.83	3.47	144.07	133.11	7.61	145.05	136.89	5.63	145.37	140.10	3.63	188.47	177.51	5.81	208.30	197.58	5.15
$W_r$	92.91	90.63	2.45	107.31	107.39	0.07	107.20	106.40	0.75	102.74	105.45	2.64	139.47	141.07	1.15	155.13	153.68	0.94
$L_{eq}$	5.16	5.04	0.14	5.96	5.97	0.07	5.96	5.91	0.04	5.71	5.86	0.15	7.75	7.84	1.15	8.62	8.54	0.94
${\rho }_s$	0.70	0.70	0.46	0.78	0.77	1.41	0.60	0.60	0.72	0.55	0.56	0.91	0.61	0.61	0.20	0.61	0.61	0.59
${\lambda }_r=250$	$R_1$	$R_2$	$\zeta$	$R_1$	R2	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$	$R_1$	$R_2$	$\zeta$
$E\left(T\right)$	128.23	127.45	0.61	172.63	153.42	11.13	164.92	146.27	11.31	170.04	152.16	10.51	215.88	190.03	11.98	254.21	211.70	16.72
$W_r$	98.23	99.92	1.72	155.41	124.10	20.15	130.98	118.18	9.77	139.65	121.48	13.01	184.04	156.66	14.88	198.91	171.99	13.53
$L_{eq}$	5.46	5.55	1.72	10.79	8.62	1.40	9.10	8.21	0.68	9.70	8.44	0.90	12.78	10.88	14.88	13.81	11.94	13.53
${\rho }_s$	0.87	0.88	0.58	0.80	0.81	0.10	0.75	0.76	0.67	0.69	0.69	0.57	0.76	0.76	0.17	0.91	0.76	16.45

, it can be concluded that the ${\mathsf{\varsigma }}_{E\left(T\right)}$ error average is 5.01%, the ${\varsigma }_{L_{eq}}$ error average is 4.58% and the ${\mathsf{\varsigma }}_{{\rho }_S}$ error average is 2.11%. With the increasing arrival rate and system scale, the error percentage also increases, reaching a value of 16.72% because, while the number of shuttles in the simulation remains constant at five, with the increasing system scale and arrival rate, the shuttles are not able to meet the requirements of high-intensity system operation, thus generating higher error rates. However, as a whole, 74.32% of the data have an error range below 5%, which is sufficient to prove that the model in this paper is able to interpret the four-way shuttle system well; the error range is shown in Figure 11.

5.2. Analysis of the Improved Solution Method

This subsection compares the solution method employing the parallel strategy with the open-loop queueing network solution method developed by Dallery [42] with respect to the general service times. This paper improves the model developed by Dallery in the following ways:

• It replaces the general service time with the Coxian distribution.
• The method of decomposing complex queuing networks is applied with the parallel operation strategy.
• The maximum entropy value method is used to calculate the open-loop queuing network.

Dallery’s paper stated that the number of customers present in the network at the same time must not exceed a given value, i.e., a limited load, which is the same as the limit on the number of shuttles described in this paper. Although it does not involve parallel queuing networks, the decomposition method is similar for complex networks, and given the same parameters, can be used for comparison and validation. The method of solving the parallel operation strategy applied in this paper was compared with Dallery’s method in nine environments with different mean values $\tau$ and squared coefficients of variation $\mathrm{c}{v}^2$. The specific environmental parameters are shown in

Table 7. Model comparison environment.

Serial Number	Configuration	Server Node 1	Server Node 2
1	$\mathrm{M}=3,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
2	$\mathrm{M}=5,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
3	$\mathrm{M}=7,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
4	$\mathrm{M}=3,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(2,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
5	$\mathrm{M}=5,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(2,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
6	$\mathrm{M}=7,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(2,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(3,5\right)$
7	$\mathrm{M}=3,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(1,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(4,5\right)$
8	$\mathrm{M}=5,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(1,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(4,5\right)$
9	$\mathrm{M}=7,\mathsf{\lambda }=0.1$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(1,5\right)$	$\left(\tau ,\mathrm{c}{v}^2\right)=\left(4,5\right)$

.

Using the system response time $E\left(T\right)$ and the external queue waiting for pickup $L_{eq}$ as parameters for calculating the simulation results and the error of the method in this paper and that of Dallery, a comparison of the data is presented in

Table 8. Comparison of model data and error rates.

Response Time						External Waiting Captain
Serial Number	Simulation	This Paper	Error (%)	Dallery	Error (%)	Simulation	This Paper	Error (%)	Dallery	Error (%)
1	22.42	22.80	0.02	16.44	0.27	0.11	0.11	0.03	0.17	0.60
2	23.12	23.90	0.03	18.23	0.21	0.08	0.09	0.07	0.15	0.82
3	22.25	24.77	0.11	19.98	0.10	0.06	0.07	0.18	0.14	1.33
4	22.92	22.62	0.01	21.76	0.05	0.15	0.12	0.17	0.21	0.45
5	23.53	23.99	0.02	23.03	0.02	0.14	0.13	0.09	0.19	0.37
6	24.30	24.86	0.02	24.76	0.02	0.12	0.13	0.05	0.18	0.50
7	25.74	24.81	0.04	26.13	0.02	0.24	0.23	0.02	0.18	0.25
8	27.17	27.80	0.02	28.83	0.06	0.19	0.22	0.13	0.17	0.12
9	28.34	29.81	0.05	30.63	0.08	0.17	0.19	0.09	0.13	0.24

(the data in the table are retained to two decimal places). This paper calculates the system response time more accurately, with an average error of 3.6%, while Dallery reported an average error of 9.2%. For the external queue waiting for pickup, the method described in this paper is better overall than Dallery, with an average error of 9.46%, while the average error of Dallery is 51.9%.

The average error for the nine environments is shown in Figure 12. Most of the error values in this paper are within 10%, with a minimum error of 1.31% and a maximum error of 18.3%, representing an overall reduction in error of 20% compared with Dallery. In summary, the methods applied both by Dallery and in this paper can be applied to solve the system performance index of the response time. However, for the external queue waiting for pickup, a comparison of the error rates between the two methods with respect to these two parameters is shown in Figure 13. The solution method described in this paper outperformed Dallery, reducing the error by approximately 40%.

5.3. Operation Strategy Analysis

In this paper, we mainly consider the combination of parallel and serial strategies for the four-way shuttle. Previous papers have mostly used serial operation strategies to represent the whole system, but in order to better reflect the actual system operation, parallel and serial operation strategies should be considered in combination. This subsection focuses on comparing the model built in this paper, incorporating the parallel operation strategy, with a model that uses the serial operation strategy to verify the effect of the operation strategy on the system.

To eliminate errors, 18 scales of the system were selected for validation and the number of shuttle cars in the system was set to 5 [43]; the selected system scale is shown in

Table 9. System size parameters.

Number of Layers	4			5			6
Number of columns	16	32	40	24	32	48	40	48	56
Specification	128	256	320	240	320	480	480	576	672
Number of layers	7			8			9
Number of columns	48	56	64	56	64	72	64	72	80
Specification	672	784	896	896	1024	1152	1152	1296	1440

.

The overall response time of the serial operation strategy is longer than that of the model developed in this paper for large and small size systems with different numbers of layers and columns, increasing from arrival rate ${\lambda }_r=50$ to arrival rate ${\lambda }_r=250$ in 50 steps. The theoretical serial-parallel combination model established in this paper tends to increase smoothly in response time, while the single serial operation strategy has a faster growth due to the increasing arrival rate. This is because the system service time of the parallel operation strategy is the maximum value of the shuttle service time and the hoist service time, while the system service time under the serial operation strategy is the sum of the two equipment service times.

In the case of a low arrival rate, the two servers, shuttle and hoist produce queues are not too long and the system response time does not result in too much time waiting for the service, so the difference between the two strategies is small. With the increasing arrival rate, the load on the hoist and shuttle increases, and there is an increase in the time that the shuttle needs to wait for the hoist when using the parallel strategy and the serial strategy. However, because the shuttle is independent of the hoist pickup when using the parallel strategy, part of the time waiting for the hoist can be spared, and the service time saved is greater than the time waiting for the hoist, while the serial strategy would cause the overall system response time to become increasingly long because of the accumulation of both the waiting time and service time.

It can be seen from Figure 14 that, at T = 4 and T = 5, there is an average of one shuttle per layer in the system, representing a fully parallel strategy; thus, the percentage reduction in the system response time is greater than that with the serial operation strategy. At T = 6, the model described in this paper, with the addition of the parallel strategy, reduces the response time by an average of 16.8% compared with the serial model, but the percentage reduction in the response time decreases with the increasing numbers of layers. At T = 9, the percentage of the performance reduction is 8.92% on average. At T = 6, T = 7, T = 8 and T = 9, a combination of 12 serial and parallel systems can reduce the system response time by 12.6% on average.

For both small-scale and large-scale systems, the model provided in this paper can reduce the response time by percentages as small as ${\lambda }_r=50$, or 1.8% at T = 9. Thus, the model in this paper can improve the flexibility of the system by reducing the response time and more accurately reflecting the operation of the actual system.

5.4. Case Analysis

This subsection is based on a real warehouse case to verify the conclusions described in this section. The case is a four-way shuttle system for auto parts run by company A, with ten layers of shelves, 64 columns per layer, a single layer deep. The total capacity of the storage system is 1280 totes, and the current system’s average daily outgoing volume is approximately 1800 totes, working 12 h a day. The other parameters of the system are shown in

Table 10. System parameters.

Parameters	Take Value	Parameters	Take Value
Vs (m/s)	2	${\gamma }_l($s)	3
${\epsilon }_s($s)	2	L (m)	0.65
${\gamma }_s($s)	5	W (m)	0.45
$Vl($m/s)	2.5	WA (m)	0.45
${\epsilon }_l($s)	1	H (m)	0.67

.

The comparison between the serial operation strategy in this system and the parallel operation strategy incorporated in this paper is shown in Figure 15. Once the arrival rate${\mathrm{reaches}\mathsf{\lambda }}_r=160$cases per hour, the system will increase in response to the arrival rate, which makes the number of customers served in the system tend towards the maximum load value, while the shuttles will all be in a busy state; thus, the increase in the waiting time causes the system response time to speed up. However, compared with a single serial strategy, the four-way shuttle system outbound flow model provided in this paper with the addition of a parallel strategy alleviates, to a certain extent, the excessive increase in the system response time. It can be seen from the data that when the average arrival rate of this system is 150 cases per hour, the addition of a parallel operation strategy can reduce the response time of the pickup operation by 4.17%. Between the arrival rate of 50 cases per hour and the arrival rate of 200 cases per hour, adding the parallel operation strategy reduces the response time of the picking operation by 4.36% on average.

6. Conclusions

The main research theme of this paper is to account for parallel strategies in a four-item shuttle system, perform system modeling, system performance evaluation and system configuration studies.

In theory, the innovation in this study is the addition of the parallelism of the hoist and shuttle in the four-shuttle system, whereas scholars have typically focused on serial operation strategy, where the response time of the system is the sum of the shuttle service time and hoist service time, thus producing large errors in the calculation and evaluation of the system performance indicators. The parallel operation strategy circumvents this drawback to a certain extent.

In terms of methodology, it is considered that the hoist and shuttle service times obey a general distribution. In the modeling process of storage systems, most scholars will assume that the server service time obeys an exponential distribution, uniform distribution, etc. However, in the storage system, the service time of customers within each server is independent of each other but obeys the same distribution, and this concept is closer to the general distribution. As the general distribution does not have a fixed distribution, the queuing network formed during the solution is not easy to calculate; therefore, this paper uses the coxian distribution to approximate the solution. Arena simulation is used to verify the accuracy of the model, and the results show that the error range of simulation results is within 10% and the error of the system performance index calculation is reduced by 20% compared with the existing methods. The parallel part is integrated into a closed-loop queueing network, which preserves the characteristics of the servers within the closed-loop queueing network to a greater extent than the method used by Dallery, thus reducing the error, and the solution method used in this paper is verified by numerical tests to reduce the error to less than 10%.

In practice, by conducting experiments and data analysis, the result is that the incorporation of a parallel operation strategy outperforms the serial operation strategy in terms of the system performance for both large-scale and small-scale systems, and this advantage rises significantly as the arrival rate increases, reducing the system response time by a minimum of 1.8% and an average of 12.6% in the numerical experiments. For the selection of the number of shuttles to be studied. the four-way shuttle system will stabilise at 80% shuttle utilisation and will not exceed the maximum system load (number of shuttles). After the shuttle utilisation is greater than 80%, the system response time and external queue waiting time will increase significantly, thus allowing 80% shuttle utilisation to be used as the optimal number of shuttles in the system at different arrival rates for the node configuration.

Shortcomings and Prospects

In this paper, the serial and parallel operation strategies are probabilistically distributed, and the “decide” module is used in the simulation. However, in real situations, there are many factors that influence the operation strategy, which are difficult to be covered by the probabilistic approach alone. The current research on automatic storage systems assumes a random storage strategy, but the existing sorting and positioning storage strategies are conducive to improving the system throughput and operational efficiency. Due to the complex distribution of aisles in the four-way shuttle system, the shuttles inevitably encounter conflicts between the shuttles in the access process. Therefore, these factors will be in the focus of future research.