A permutation-combination heuristics for crane-based automated storage and retrieval systems considering order fulfillment time and energy consumption

1.   Introduction

Since the invention of the automated storage and retrieval system (AS/RS) in the 1950s, it has been widely used worldwide. The AS/RS is a key component of enterprise supply chain logistics, and is advantageous considering its high space utilization, reduced labor costs, short retrieval times and better inventory control ^[1]. The usage of AS/RS has grown rapidly with the growth of e-commerce and the increasing demand for manufacturing. In order to provide AS/RS users with greater economic benefits, AS/RS manufacturers usually focus on continuously improving the throughput efficiency of the system. Furthermore, AS/RS users such as manufacturing enterprises and e-commerce companies care about the economic benefits it brings. As a consequence, more and more high-speed AS/RS have been developed and put into use, and some AS/RS stacker cranes can even reach an astonishing speed of 3.5m/s ^[2]. Higher AS/RS speeds often mean higher power motors, more energy consumption and more carbon emission when the system operates. Therefore, it is necessary to consider the energy consumption of the AS/RS in the design and operation process. According to the classification standards of storage and retrieval machines (S/R machines), common AS/RS can be roughly divided into shuttle-based storage and retrieval system (SBS/RS) and crane-based automated storage and retrieval system (CB-AS/RS). Among them, the SBS/RS mainstream handling machines are shuttles and elevators. Elevators can transport the goods to the corresponding layer, while single-tier shuttles can receive the goods at the corresponding level and move them to the designated storage location. There is also another type of multi-tier shuttles, which can move vertically between two or more aisles and complete the storage and retrieval of goods ^[3]. The CB-AS/RS is a traditional AS/RS, whose S/R machine is only a crane. This type of crane can move freely in the storage rack aisle and use the loaded fork to achieve the handling task of the goods on both sides of the storage rack. Usually, each aisle of the rack would be equipped with only one crane, but that single crane can operate with multiple items simultaneously to improve operational efficiency.

As society develops, people's awareness of protecting the environment and saving energy is increasing, and how to reduce AS/RS's energy consumption during operation has become a hot spot for more and more researchers. Lerher et al. ^[2] have conducted research earlier in the energy consumption problem and presented a model for miniload AS/RS based on warehouse throughput, mechanical analysis of S/R machines and carbon dioxide emissions. They also studied the impact of different warehouse volumes and different speeds and accelerations of the S/R machines on warehouse performance. They thought it would be unwise to blindly increase the speed and acceleration of the S/R machine for some large storage shelves, whereas they considered it advisable to increase the number of the S/R machines. The performance metrics of the AS/RS are defined primarily by the unit-time throughput, or the time consumed to complete a certain amount of delivery tasks. Therefore, while reducing energy consumption, it is necessary to ensure that the performance of the warehouse does not suffer too much losses. Borovinšek et al. ^[4] established an optimization objective model for the SBS/RS with the goal of minimizing average throughput time, energy consumption and total investment costs. They considered the number of channels, layers, columns, shuttle speeds and accelerations, as well as elevator speeds and accelerations in the optimization process, and used the NSGA2 algorithm to solve and optimize the model. They provided an optimization idea for the early stage of warehouse design. Liu et al. ^[5] studied the impact of speeds and accelerations on the energy consumption of elevators and shuttles of the SBS/RS in different configurations under determined throughput requirements. The results showed that the speeds and accelerations of vertical motions have the most significant impact on energy consumption. Also, they concluded that there is a high correlation between throughput and energy consumption. They established a decision-making basis for energy consumption optimization for warehouse designers and SBS/RS managers.

In the AS/RS, it is both the pain point as well as the focus for researchers to establish a mathematical model that fully conforms to reality. The SBS/RS carries relatively light loads and belongs to the miniload type. The shuttles used by the SBS/RS have lower power and are more energy-saving than those of the CB-AS/RS. Thus, there is much more researches on the mathematical models of the SBS/RS than the CB-AS/RS. However, for a heavier category of items, companies are more inclined to use a CB-AS/RS. Therefore, we establish a location optimization model for the existing CB-AS/RS to minimize energy consumption and order fulfillment time. Calculations of the order completion time and crane energy consumption are done from different perspectives, which specifically speaking are specifical input and output instructions. The GA is used to solve the model by spatial mapping coding. A permutation-combination heuristics is proposed, which follows the coding method and cross-operation of the genetic algorithm and conducts both arrange-operation and change-operation.

2.   Literature review

The AS/RS combines modern intelligent technology with traditional warehouse management, and upgrades the warehouse structure and management in an all-round way, which is the development direction of intelligent warehousing. From a structural point of view, the AS/RS is in general composed of aisles with storage racks (SR) on both sides, storage and retrieval machines (S/R machines), input and output (I/O) locations and pick locations ^[6].

With the increase in enterprise demand and the development of e-commerce trade, the AS/RS has achieved rapid development in recent years. Research mainly aims at increasing the performance of the AS/RS as well as increasing the efficiency of operation and reducing the costs of enterprises. The location allocation and storage policy are the hotspots in these attempts. Jiang et al. ^[7] proposed a scattered storage strategy considering the correlation between products. They also developed a mixed genetic algorithm and particle swarm algorithm to solve the storage model, The results proved that the hybrid algorithm can greatly improve the quality of model. Li et al. ^[8] used the product affinity heuristic algorithm of data mining to calculate the correlation among the items to solve the dynamic storage allocation problem. They proposed a greedy genetic algorithm to calculate the model. The results showed that the proposed model based on data mining improves the average picking distance by 7–104% compared with the traditional ABC classification model. Liu and Poh ^[9] proposed a scattered storage space allocation method suitable for Internet distribution warehouse operations in the field of e-commerce retail. The storage model was described as a mixed integer programming (MIP) model, and a heuristic algorithm including a feasibility check and two-stage optimal location allocation is developed for this model. Through the algorithm performance test and the application model efficiency test, it was proven that the proposed model and algorithm improve the picking efficiency. He et al. ^[10] studied three-rack strategies of an automated stereoscopic warehouse from the perspective of space utilization and proposed a real matrix genetic algorithm, including first-fitting and best-fitting location allocation. The experimental results showed that the adaptive shelf strategy and the best-fitting shared shelf strategy can better meet the needs of small inventory and high-density warehouse operations. Yang et al. ^[11] established the optimization function from the perspective of shelf stability and inbound and outbound efficiency of an automated stereoscopic warehouse. They used the genetic algorithms to prove that optimizing the cargo space and the running route of the stacker crane can improve the overall operation efficiency. Chen et al. ^[12] proposed several heuristic repositioning strategies inspired by traditional storage location allocation strategies for the warehousing process of e-commerce. They set up a dynamic scheme to relocate the storage position of non-empty pallets. They also proposed a migration strategy that approximated dynamic programming to improve the efficiency of retrieval operations based on the experiment results. In summary, most researchers use the heuristic algorithms to solve the location allocation problem.

With the strengthening of people's concept of green development, more and more researchers have begun to pay attention to the energy consumption of AS/RS. Meneghetti et al. ^[13] studied the combination of racks of different heights with the cranes of different specifications to optimize time and energy consumption simultaneously. An overall optimization model based on constraint programming and large neighborhood search was proposed, which combined the best storage allocation strategy and the optimal ranking method for joint application for the first time. The results show that the medium-height shelf-shape rack can achieve the best performance in energy efficiency. Lower shelf-shape racks perform better in time performance. Ekren ^[14] took the number of layers, columns, maximum speed and acceleration that can be achieved by the elevators and shuttles as design variables, and studies their impacts on the average cycle time, energy consumption and energy regeneration of the SBS/RS's transactions. She also considered the influence of interactions between design variables on the results. Through the simulation with ARENA, the level of influence of each design variable on warehouse performance was obtained. Hsu et al. ^[15] proposed a framework for the close connection between simulation and optimization of double-deep AS/RS. Taking the minimum energy consumption as the optimization goal, they used a variety of common heuristic algorithms to solve the problem. Moreover, they proposed an improved whale optimization algorithm (IWOA) with adaptive motion and mutation characteristics, and proved that the IWOA combined with dynamic programming can better solve the crane scheduling problem in the double-deep AS/RS.

Lerher et al. ^[2] did not calculate energy consumption with actual input and output orders, but obtained the average power by taking the square root of the motor power at acceleration, deceleration and constant speed, through the process of which the theoretical energy consumption was calculated. Borovinšek et al. ^[4] also obtained the average power by taking the square root. He achieved the theoretical solution based on probability theory and reached at the expected energy consumption of the running time. Liu et al. ^[5] used the average time of acceleration, deceleration and constant speed to obtain energy consumption. It is hard to judge which calculation method is better using different models. In this paper, we calculate energy consumption based on specific orders and obtain energy consumption for each input and output instruction. In order to get closer to the real warehouse operating state, the outbound and inbound instructions of specific orders are randomly generated by computer simulation. Furthermore, to solve the problem that there are few research papers on the CB-AS/RS's energy consumption, we establish a location allocation model based on the CB-AS/RS to minimize order fulfillment time and energy consumption. The GA is used to solve the model by spatial mapping coding, and a PCH algorithm is proposed to solve the model. The influence of different storage utilization rates and different input and output instruction quantities in a batch of orders on the results is also considered in the simulation process.

3.   Model for designing crane-based storage and retrieval systems

The following assumptions are proposed to develop an order fulfillment time and energy consumption model for the CB-AS/RS.

1) Ignore the time consumed when the forks load and unload the items

2) Ignore the energy consumption of the forks when loading and unloading the items

3) Ignore the resistance of rotating masses when the crane accelerates and decelerates

4) Ignore the circuit heating, noise and other consumption of the crane

For assumptions 1), 2) and 4), most researches ^[2,4,5] ignore circuit heating, noise and the energy and time consumption of forks. Their explanation is that the time and energy consumption consumed to load and unload the same order can be considered roughly the same, and it is also irreducible. The consumption of heat and noise is theoretically difficult to calculate, as it only accounts for a small portion of the overall consumption of the crane. Therefore, in order to simplify the calculation, this article also ignores the above consumption.

For assumption 3), Ekren et al. ^[16] considered resistance of rotating masses with variable speed in the SBS/RS model. They introduced factor (fr = 1.15) into the resistance calculation process and multiply the factor by the overall mass. Finally, they achieved results with high accuracy. However, the CB-AS/RS proposed in this article is different from the SBS/RS. The acceleration and deceleration of the crane model studied in this article are much lower than that of the SBS/RS. Second, if this article directly refers to their method, which is multiplying the overall mass by a factor during calculation. This may yield results that may not be applicable to this model. Finally, and most importantly, it is difficult to obtain the precise mass of rotating components without dismantling the crane. Thus, in summary, we will not consider the resistance of rotating masses in the calculation process and adopt random storage strategy. The decision is to highlight the differences in model calculations between the GA and the PCH.

To facilitate analysis, the symbols are introduced to represent parameters (See Appendix). The parameters consist of three parts: basic parameters, auxiliary variables and decision variables.

3.1.   Order fulfillment time model

In this paper, the CB-AS/RS warehouse racks are designated to be single-deep racks, and the crane can handle one location unit per operation. It is necessary to take into account whether the crane is operating in single command (SC) or double command (DC) mode. Take SC mode as an example: The crane operation can be decomposed into vertical movement and horizontal movement, and they do not interfere with each other. The crane, from the starting point to the target point, has undergone acceleration and deceleration from standstill, and if the maximum running speed is reached in the process, there would be a constant speed process in the middle. Therefore, the movement in both directions needs to consider whether the maximum speed can be reached. The running time of the crane between any point i in the racks and point O is calculated as follows, where point O is placed in the I/O position. Point i and j are the locations in the racks and (i_x, i_y) represents the i_x column and i_y tier in the racks.

If $D{H_{oi}} = \left| {o - {i_x}} \right|*Ws \leqslant \frac{{Vx{{\ max }^2}}}{{{a_x}}}$, the crane has not reached its max horizontal velocity, as shown by the solid line a in Figure 1, Ws represents the width of a storage location, then

If $D{H_{oi}} = \left| {o - {i_x}} \right|*Ws > \frac{{Vx{{\ max }^2}}}{{{a_x}}}$, the crane reaches its max horizontal velocity, as indicated by the dotted line b in Figure 1, then

The same process can be applied to calculate the vertical movement time.

If $D{V_{oi}} = \left| {o - {i_y}} \right|*He \leqslant \frac{{Vy{{\ max }^2}}}{{{a_y}}}$, He represents the height of a tier, then

If $D{V_{oi}} = \left| {o - {i_y}} \right|*He > \frac{{Vy{{\ max }^2}}}{{{a_y}}}$, then

3.1.1.   Order fulfillment time calculation in SC mode

In SC mode, the crane first loads the item at the I/O location, then moves to the input position to unload the item. Finally, it returns to the I/O location. The running route is shown in the blue line in Figure 2, so the SC cycle time can be obtained.

3.1.2.   Order fulfillment time calculation in DC mode

The DC mode operation route is shown in the red line in Figure 2. In this mode, the crane first loads the item at the I/O location and moves to the input location (i_x, i_y) to unload the items. Then, the crane runs to the (j_x, j_y) position to load the output item, and finally returns to the I/O point to unload the item. Thus, the DC cycle time can be obtained as follows.

The running time required for one crane with a batch of orders can be finally obtained as follows. a represents the number of the aisles, whose amount is equal to those of the cranes.

The order fulfillment time depends on the longest running time of all the stacker cranes involved in the order. When all cranes work together:

3.2.   Energy consumption model

The calculation of the energy consumption caused by crane operations can be divided into horizontal energy consumption and vertical energy consumption. Taking the SC as an example, the traction force produced in the process of accelerating, decelerating and constant speed movement in horizontal direction can be obtained through mechanical analysis: F_Ha, F_Hd, F_Hc. The power is calculated respectively: P_Ha, P_Hd, P_Hc.

The same process can be applied to the analysis of the traction force produced in the vertical direction.

when the crane is fully loaded, m = m_c + m_s; when the crane is empty loaded, m = m_c. m_c represents the overall mass of the crane, m_v represents the mass of the lifting device. According to the schematic diagram of the stacker crane in Figure 3, it is clear that m_c includes upper crossbeam, lower crossbeam, left column, right column, carriage platform, forks, vertical travel motor, horizontal travel motor and electrical control cabinet, while m_v includes carriage platform and forks.

3.2.1.   Calculation of SC cycle energy consumption

For one SC cycle, the crane if fully loaded in the former half of the whole process and empty loaded in the other way. The energy consumption can be calculated according to the formula W = Fvt.

The product of all vt in the above formula is the distance. It should be noted that the velocity is non-constant during acceleration and deceleration. The e is the energy conversion efficiency.

If $D{H_{oi}} = \left| {o - {i_x}} \right|*Ws \leqslant \frac{{Vx{{\ max }^2}}}{{{a_x}}}$, then v_xt_Haoi and v_xt_Hdoi are the horizontal distances of acceleration and deceleration, and ${v_x}{t_{Hdoi}} = {v_x}{t_{Hdoi}} = \frac{1}{2}\left| {o - {i_x}} \right|Ws$.

If $D{H_{oi}} = \left| {o - {i_x}} \right|*Ws > \frac{{Vx{{\ max }^2}}}{{{a_x}}}$, then ${v_{xmax}}{t_{Hcoi}} = \left| {o - {i_x}} \right|*Ws - \frac{{Vx{{\ max }^2}}}{{{a_x}}}$.

If $D{V_{oi}} = \left| {o - {i_y}} \right|*He \leqslant \frac{{Vy{{\ max }^2}}}{{{a_y}}}$, then v_yt_Vaoi and v_yt_Vdoi are the vertical distances of acceleration and deceleration, and ${v_{ {y}}}{t_{Vdoi}} = {v_y}{t_{Vdoi}} = \frac{1}{2}\left| {o - {i_y}} \right|*He$

If $D{V_{oi}} = \left| {o - {i_y}} \right|*He > \frac{{Vy{{\ max }^2}}}{{{a_y}}}$, then${v_{ymax}}{t_{Vcoi}} = \left| {o - {i_y}} \right|*He - \frac{{Vy{{\ max }^2}}}{{{a_y}}}$

The EC_HEoi and EC_VEoi can be calculated in the same way. Thus, no further explanation is required.

3.2.2.   Calculation of the DC cycle energy consumption

In DC mode, the crane first loads the item at the I/O location and moves to the input location (i_x, i_y) to unload the item. Then, the crane runs to the (j_x, j_y) position to load the output item, and finally returns to the I/O point to unload the item. After obtaining the operation status of each step, the Eq (30) can be reached. Because each step of the calculation is the same as that of the SC mode, no further elaborated explanation would be presented.

After calculating the energy consumption of each crane, the energy consumption of the crane with one batch of orders is finally summed up as shown in Eqs (31) and (32).

3.3.   Objective function

The model proposed in this paper considers both order fulfillment time and energy consumption of the cranes. By introducing weight coefficients, the multi-objective optimization problem is transformed into a single-objective problem. The specific steps are:

a. Eliminate the dimensions and align the order of magnitude of the two optimization goals.

b. Define the importance of the goals by weight coefficient.

The objective function is as follows: min(EC) and min(T_F) are the minimum values obtained by the genetic algorithm after multiple calculation. The former takes energy consumption as the single-objective function, and the latter takes order fulfillment time as the single-objective function.

4.   Solution procedure

The model proposed in this paper includes input and output operations. For output, the core issue lies in the arrangement and combination of output and input instructions. For input, the core issue lies in the location allocation, which has been proven to be an NP-Hard problem. Therefore, the popular genetic algorithm is used to solve the model proposed in this paper. Furthermore, based on the basic logic of the heuristic algorithm, a permutation and combination heuristic algorithm is proposed to solve this model.

4.1.   Procedure of the genetic algorithm

The genetic algorithm is a heuristic optimization algorithm formed by imitating the evolution process of organisms in nature. From a macroscopic point of view, evolution ensures that the living organisms can adapt to their living environment. Individuals with strong survival abilities in the process of population reproduction will have a greater probability of being retained by the environment. Microscopically speaking, evolution means a process in which biological chromosomes continue to copy, cross and mutate to produce new individuals, whereas new individuals are filtered by the environment with the fittest ones surviving. In the genetic algorithm, the digital information of each individual of the population is retained through the initial population definition and coding method. Each iteration of the algorithm corresponds to a filtering of the individual by the environment. The digital information is copied, crossed and mutated again and again to increase the fitness function value, corresponding to the improvement of the adaptability of the biological group to the environment. The basic process of the genetic algorithm includes encoding, population initialization, copy operation, crossover operation and mutation operation, as shown in Figure 4. These aspects will be further described below.

4.1.1.   Spatial mapping encoding

In this article, we use spatial mapping coding to encode the location of the warehouse in real numbers. Each actual storage location (x, y, z) of the warehouse has only one rack number N that corresponds to it. This method can transform the three-dimensional location coordinate information into one-dimensional information and at the same time reduce the length of a feasible solution. Taking A, M, C as an example, because in each aisle a crane can handle the items on two rows of the racks, the warehouse rack information is column C, row 2*A, tier M, (x, y, z) is represented as the x column, y row and z tier in the rack, and the specific spatial mapping rule is:

Digital information should not only exhibit the input instructions but also the output instructions. The one-dimensional information of only one bit cannot fully represent the input and output operation instructions. Therefore, an additional identifier bit is introduced to represent the distinction between input instructions and output instructions. For example, (identifier, N) is an operation instruction. The identifier belongs to 0 or 1, where 1 represents an input instruction and 0 represents an output instruction. (1, 1120) represents that the next input item will be stored on the rack number 1120, and (0, 1886) means that the item will be taken out of the rack number 1886. For one batch of orders, the known output instructions and random input instructions constitute a feasible solution. In the feasible solution, after being decoded, two instructions on the same channel (instructions tracing back to the same crane) can form a DC cycle instruction. The former would be an input instruction and the latter an output instruction.

4.1.2.   Initialize the population

The input and output instructions of a batch of orders are randomly arranged to form an individual, which can also be called a feasible solution. Multiple individuals form a population. Usually, a feasible solution is not the state with the largest number of DC instructions. It is necessary to collect the instructions belonging to the same stacker in the feasible solution, form as many DC instructions as possible according to the order of first input and then output, and combine them into sub-feasible solutions in crane units. To a certain extent, the efficiency of the crane operation depends on the number of DC instructions, so forming as many DC instructions as possible can improve the efficiency of the cranes per time unit. The sub-feasible solutions are combined into a complete feasible solution. This is the process of initializing the population. The decode process is as follows:

In Eqs (35) and (36), the function ceil is the round toward positive infinity

4.1.3.   Copy operation: tournament selection method

The basic idea of the tournament selection method is to randomly select a determined number of individuals (feasible solutions) from the parent population, and put them into the fitness function for calculation. Then, the results are compared and individuals with good performance are selected. The operation is repeated until the number of selected individuals reaches the specified number, and these are sent to the next operation.

4.1.4.   Crossover operation: single-point self-crossover

A complete instruction consists of two bits. If the cross-breakpoint is in the middle of an instruction, the crossover process will produce a bad solution. Thus, only two consecutive gene fragments are allowed to cross during the crossover process. Considering that the crossing of two random chromosomes (solutions) in the population will cause duplicate location coordinates, this paper would splits the solution into two sub-solutions for crossing. Taking a solution containing twelve instructions as an example, the solution is divided into the sub-solution A and the sub-solution B of equal length. Here, the solution contains an even number of instructions. If there is an odd number of the solutions, then the sub-solution A would be rendered one more instruction than the sub-solution B. It is necessary to ensure that any instruction remains complete without being split. Then, perform the single point crossing. The specific process is shown in Figure 5. The rack number in the sub-solution of the crane is calculated while A = 8, M = 15 and C = 36. The instructions of different colors in the figure represent the instructions belonging to the same crane. Taking three cranes as an example, each crane only contains four instructions, whose rack numbers are within the condition range.

Under real working conditions, the number of instructions that a crane needs to process in an order will be much more than four, and the output instructions and input instructions are not necessarily equal. The simplification made here is to better explain the rules of the algorithm. Only if the sub-solutions A and B come from the same solution, it can be ensured that the same location genes do not appear during the crossover process.

4.1.5.   Mutation operation

The crossover operation restricts the point of the cross and packs the valid gene fragment. This operation does not create new genes in the crossover process, making it easier for the algorithm to reach the local optimum. Therefore, the mutation operation can continue to escape the touch of local optimum by appropriately increasing the mutation probability. The mutation process is valid only for the input instructions, while the output instructions are determined by the demand of the order. The rack numbers of the items in the output order are determined by the locations stored in the input process. Except for collating inventory, warehouse transfers, etc., the stored items cannot be moved by default. The basic idea of the mutation process is to change the rack numbers in the input instructions. It needs to be explained that the changed rack numbers cannot be used, including those that have already been put into use or been stored in the input instructions prior to the designated input instructions in the process.

Practically, if the number of remaining racks is relatively large, the efficiency of randomly exchanging the remaining rack numbers will be comparatively low. Even worse, it would be difficult to find a rack number with better performance. Therefore, in the mutation process, the algorithm randomly seeks for a new rack number nearby the initial one to see whether it shows a better performance. After the exchange process, if the new rack number is found to be once used before this entire order, then this process would be judged as invalid. If the new rack number is found to be once used in the input instruction, then the rack number with a better performance is preserved regardless of the initial sequences. The previous rack number is replaced if its performance is inferior, or retained if it is superior.

4.2.   Procedure of permutation-combination heuristics

The genetic algorithm is inevitably subject to slow convergence speed and would easily fall into the local optimum. Therefore, we propose a permutation-combination heuristic (PCH) based on genetic algorithms while at the same time similar to most heuristic algorithms. The basic processes of PCH are: encoding, initialization of feasible solution, arrange operation, offset operation, cross operation and change operation, as shown in Figure 6. Encoding and crossover operation are the same as those of the genetic algorithm described above.

4.2.1.   Arrange operation

This operation is basically the same as the process of initializing a feasible solution in the GA as mentioned above, but it is not a random combination. First, the instructions belonging to the same crane in the solution are collected. Second, the instructions are taken out one by one according to the serial number of the crane and combined into a feasible solution, until all instructions are taken, as shown in Part I of Figure 7.

4.2.2.   Offset operation

The purpose of the offset operation is to search for the best combination of output instructions and input instructions. The specific operation includes fixing the output instructions, collecting input instructions orderly, randomly generating the offset bit (the italic offset represents the number of digits), letting the input instructions offset to the left, then filling in the offset input instructions on the left to the right position in turn and finally filling the input instructions back in the original positions after the offset process is completed. As shown in Part II of Figure 7, the offset is 2 bits. The restriction of the offset is: 0 < offset < = N_S, where N_S represents the number of input instructions in a batch of orders.

4.2.3.   Change operation

The mutation operator here is also valid only for input instructions. The operation is divided into two parts: The first part is the exchange of instructions between the cranes; the second part is the change of instructions belonging to the same crane.

The specific process of the first part includes first counting the number of instructions processed by each crane in this batch of orders, and randomly selecting an input instruction from the crane with the largest number of instructions and changing this input instruction into a new instruction that belongs to the crane with the smallest number of instructions. This is to balance the intensity of work between cranes.

From the modeling point of view, it can be obtained that for multiple items of determined mass, the larger the mass stored in the lower tiers and the smaller the mass stored in the higher tiers, the less energy the crane would consume. The energy demand is also lower. The closer the storage location is to the I/O location, the shorter the vertical and horizontal movement distance of the crane and the shorter the order completion time.

Based on the above analysis, the second part of the process is obtained. From a three-dimensional perspective, the second part of the process can be concluded to be: Moving the input storage location downward while moving back and forth on that row of the rack to reach the optimum. From the perspective of digital information, moving downward means rack number(N) minus 2AC, and moving forward or backward means rack number(N) adding 1 or subtract 1. When the storage location is on the first tier (z = 1), the rack number cannot subtract 2AC; when the storage location is at the head of the column (x = 1) or the end of the column (x = C), the rack number cannot minus 1 or add 1. After the change, it is necessary to verify that the changed rack number whether has been used. If it has been used by the stored items, this change would be judged as invalid. If the new rack number is found to be once used in the input instruction, then the rack number with a better performance would be preserved. The previous rack number would be replaced if its performance is inferior, or retained if it is superior. The specific process is shown in Figure 8.

5.   Numerical study

In this section, numerical simulation experiments are carried out using MATLAB. The basic parameters of the racks are: A = 8, M = 15, C = 36, W_S = 2m, He = 1m. The basic parameters come from a reducer assembly warehouse in China. The basic parameters of the cranes are: Vxmax = 2m/s, a_x = 0.3m²/s, Vymax = 0.5m/s, a_y = 0.3 m²/s, m_c = 3400kg, m_v = 600kg, provided by crane manufacturers. The energy conversion efficiency is e = 0.8, the gravitational acceleration is g = 9.8N/kg and the storage utilization rate is 25%. The mass m_s of the items are generated by the algorithm. The range of the mass is between 800 and 1,100kg conforming a uniform distribution. Considering the uncertainty of the number of input instructions and output instructions in a batch of orders, the experiment adopts three matches of different N_S and N_R, specifically Group A: N_S = 50, N_R = 100; Group B: N_S = 100, N_R = 100; and Group C: N_S = 100, N_R = 50.

Finally, considering the unstableness of the results solved by the heuristic algorithm, the genetic algorithm and the PCH algorithm perform repetitive experiments for ten times. The mass of the input and output items and the location of the output items remain constant in the same group of repetitive experiments. Generally speaking, the settings of parameters in the heuristic algorithm have a great impact on its performance. Therefore, after literature review ^[17,18], multiple tests and discussions with relevant researchers, we present the settings of the parameters as follows.

For both the GA and the PCH, the population size is 100, the number of iterations is 100 and the crossover probability is 0.8. For the GA, the mutation probability is 0.3. The experiment results are shown in Table 1. T_F is measured in seconds (S) and EC is measured in kilowatt-hours (kWh). In this table, the standard values of EC and T_F are the minimum values that are solved by the genetic algorithm ten times. The EC uses energy consumption as the single-objective function, and the T_F uses order fulfillment time as the single-objective function. The T_F is rounded to one decimal place and the EC is rounded to two decimal places. For Diff, it is displayed as a percentage and rounded to three decimal places. These experiments were performed on a computer with Intel Core i5-12400F (64 bits and Performance-core Base Frequency 2.5 GHz) and 8*2 GB DRAMs.

Tables 1 and 2 are the experiment results of the GA and the PCH at 25% storage utilization rate of groups A, B and C, respectively. The fonts for the highest and lowest values of each group in the experiments in the table have been set in bold. The Diff is the percentage difference between each result and the standard value. The smaller the Diff, the better the results. According to the data in the table, the Diff of the order fulfillment time (T_F) obtained by the PCH of groups A, B and C is lower than that obtained by the GA, and the Diff value of the energy consumption (EC) obtained by the PCH is also lower than the result obtained by the GA.

It is observed that the order fulfillment time and the energy consumption obtained by the GA fluctuate greatly, while the PCH results are more stable. In order to visualize the difference in results, Figures 9 and 10 show the box-whisker plots of T_F and EC based on the data in Tables 1 and 2, where the light-blue box represents the PCH and the light- brown box represents the GA. Each box contains 50% of the data results of each experiment. It can be clearly observed that the PCH results of the groups A, B and C are lower than those of the GA. Except for the box obtained by the PCH for group C in Figure 10, which is not flatter than the GA, the boxes obtained by the PCH in all groups are flatter than those obtained by the GA under the same configuration, which further indicates that the results obtained by the PCH are more stable. In summary, the PCH not only outperforms the GA at 25% storage utilization rate, but is also more stable.

In order to explore the influence of different storage utilization rates on the algorithm solution, the storage racks with 50 and 75% storage utilization rates are simulated separately. Tables 3 and 4 are the experiment results of the GA and PCH with 50% storage utilization.

Tables 3 and 4 show that the T_F and EC results of group A and group B by the PCH are lower than those of the GA. In group C, the minimum Diff of T_F and EC obtained by the GA solution is smaller than the maximum Diff obtained by the PCH, indicating that it cannot be directly decided which algorithm is better by comparing the experiment results. However, the mean Diff obtained by the PCH solution in group C is lower than that of the GA: The mean Diff of T_F obtained by the PCH is –0.834%, which is less than 4.218% obtained by the GA; the mean Diff of EC is 1.438%, which is also less than 1.830% obtained by the GA.

Figures 11 and 12 are box-whisker plots of T_F and EC based on the data in Tables 3 and 4, respectively. It is observed that the performance of groups A and B is consistent with the previous analysis, and the boxes of T_F and EC obtained by the PCH are not only lower than those obtained by the GA, but also flatter. This indicates that the results obtained by the PCH are better and more stable than those obtained by the GA. In group C, the box of T_F obtained by the PCH is lower than the box obtained by the GA, indicating that the PCH is still better than the GA in T_F. The box of EC obtained by the PCH is slightly lower than the GA, but as a whole it is taller than the box obtained by the GA. The PCH is better than the GA based on the median of the two boxes (indicated by the horizontal line in the box) and the mean of the results. Through the previous analysis, it can be concluded that the optimization performance of the PCH in group C has weakened at 50% storage utilization rate.

Tables 5 and 6 show the experiment results of the GA and the PCH at 75% storage utilization rate, respectively. In group A, it is observed that the Diff of T_F and EC obtained by the PCH are lower than that obtained by the GA in each experiment. In group B, the Diff of T_F obtained by the PCH is lower than that of the GA in each experiment, but the Diff of EC obtained by the PCH is not entirely lower than that of the GA. This means that the optimal solution obtained by the GA solution is better than the worst solution solved by the PCH. In group C, the Diff of T_F and EC obtained by the PCH is also not entirely lower than those of the GA. At this point, the optimization ability of the PCH has been further weakened.

Figures 13 and 14 demonstrate box-whisker plots of T_F and EC based on the data in Tables 5 and 6. In group A, it is observed that the boxes of T_F and EC obtained by the PCH are not only lower than the boxes obtained by the GA, but also flatter. In group B, the box of T_F obtained by the PCH is lower than the box obtained by the GA, but the EC solved by the PCH has a numerical overlap with the box obtained by the GA. In group C, the boxes of T_F and EC obtained by the PCH both have a numerical overlap with the boxes obtained by the GA. With the exception of the box of EC obtained by the PCH in group C, all boxes obtained by the PCH are lower than those of the GA. The box of EC obtained by the PCH in group C is higher than the box obtained by the GA, but the median of the PCH box (indicated by the horizontal line in the box) is slightly lower than that of the GA, whereas the value remains very close. It can be concluded that the optimization ability of the PCH for EC in group C at 75% storage utilization rate is in the same level with or even inferior to the GA. Considering the two factors of T_F and EC, the PCH can still be considered slightly better than the GA in group C.

In Table 6, the T_F of group A has multiple duplicate results. Theoretically speaking, the results solved by the heuristic algorithm are generally different. Duplicate results can infer a local optimal solution or the global optimal solution in the value. It is generally very difficult to determine whether the solution is the local optimal solution or the global optimal solution. Therefore, we do not put a major emphasis on this topic.

Figures 15 and 16 are the optimization rate bar charts of order fulfillment time and energy consumption. The optimization rate calculation in the figures is done by taking the difference between the average values of the Diff obtained by the PCH and the GA. In group A, it can be observed that the optimization rate of T_F increases slightly with the increases of storage utilization, and the optimization rate of EC decreases first and then increases. In groups B and C, the optimization rates of T_F and EC decrease to varying degrees with the increase of storage utilization. Further analysis shows that the optimization ability of the PCH in groups B and C weakens with the increase of storage utilization.

The comprehensive optimization rate is drawn by weighting the average of the optimization rates of T_F and EC, as shown in Figure 17. It can be observed that the comprehensive optimization rates of groups B and C decrease, and the comprehensive optimization rate of group A first decreases and then increases. A possible reason might be that the increase in storage utilization has reduced the number of storage racks that can be selected for the input items and the reduction of optimization space would lead to a decrease in the optimization rate. The number of input instructions in groups B and C is 100. The number of input instructions in group A is 50, which has less optimization potential than groups B and C.

In Table 7, the solution time consumed by the PCH is less than that of the GA at 25% storage utilization rate. Especially in group A, the solution time of the PCH is saved by 48.81%, in which case the PCH has a huge advantage. At 50% storage utilization rate, the solution times of the PCH and GA is roughly the same. At 75% storage utilization rate, the PCH solution time is significantly higher than the GA solution time. Especially in group C, the solution time of the PCH is 51.14% more than that of the GA. However, the comprehensive optimization rate of the PCH is 1.55% higher than the GA. At this point, it is up to the AS/RS's manager to determine whether to spend additional time to solve for better system performance, or to use an upgraded computer with higher frequency and more cores to solve it.

Comparing the solution time of groups A and C, it takes more time to solve group C. The reason for this is that there are more N_S in the input instructions in group C. From the algorithmic level, the mutation operation of the GA and the change operation of the PCH are only valid for the input instructions. The increase of N_S in the input instructions has a more significant impact on the solution time of the algorithm than on N_R in the output instructions.

6.   Conclusions

In CB-AS/RS, the energy efficiency of the cranes is one of the most important performance metric. In this paper, we establish a location optimization model for the CB-AS/RS to minimize energy consumption and order fulfillment time. Here, we use spatial mapping coding to implement the GA to find the solution. A new PCH algorithm is proposed at the same time, whose basic process includes encoding, initializing a feasible solution, arrange operation, offset operation, cross operation and change operation. In real AS/RS operations, the items are constantly removed and replenished at the same time. It is difficult to predict the specific output and input quantities of the items in a batch of orders. Therefore, we conduct experiments with three sets of different output and input quantities, namely groups A, B and C. The essence of the PCH and the GA is not only permutation and combination of the existing operators, but also the exploration for the unknown operators. The PCH algorithm focuses more on the permutation and combination of input operators and takes into account the scheduling problem of the cranes (reflected in the change operation). Based on the experiment results, the following conclusions can be drawn.

1) At 25 and 50% storage utilization rates, the PCH performs better than the GA for groups A (N_S = 50, N_R = 100), B (N_S = 100, N_R = 100) and C (N_S = 100, N_R = 50) in solving T_F and EC, and is more stable. At 75% storage utilization rate, the PCH performs at an average level when solving EC in group C. After comprehensive consideration of T_F, the overall performance of the PCH is still considered better than that of the GA. Analyzing the T_F alone, the PCH solution results are better than those of the GA by 1.52% to 13.6%. From the EC analysis, the PCH is better than the GA by 0.02 to 5.58%. From the analysis of comprehensive optimization rate, the PCH is better than the GA by 1.55 to 9.59% in the specific value. All in all, the PCH is overall better than the GA without considering algorithm solution time.

2) In group A (N_S = 50, N_R = 100), with the increase of storage utilization rate at 25, 50 and 75%, the optimization ability of the PCH first decreases and then increases. In groups B (N_S = 100, N_R = 100) and C (N_S = 100, N_R = 50), the optimization ability of the PCH has been continuously weakening as the storage utilization rates increase.

3) Except for group C at 50% storage utilization rate, the solution time of the PCH is less than the GA at 25 and 50% storage utilization rates. At 75% storage utilization rate, the solution time of the PCH is higher than with the GA, and the largest increase in group C reaches 51.14%.

4) The solution time of the GA and the PCH is increasing with the storage utilization rate. The reason is that the climbing storage utilization increases the number of shelves that are used, which will make the mutation operation in the GA and the change operation in the PCH produce more repeated position judgment loops.

The major contributions of this paper are as follows:

1) Taking the existing CB-AS/RS as the research object, a location allocation model is established to minimize order completion time and energy consumption.

2) A Permutation-Combination Heuristics is proposed, which follows the coding method and crossover operation of the genetic algorithm, while at the same time it adds arrange operation and change operation. We offer an optimization idea for the researchers and administrators who study the CB-AS/RS.

Many experts have also pointed out the shortcomings of this study, such as the model may deviate from reality, and the hypothesis given may not be in line with the actual situation. We will continue to optimize the existing model to better fit the reality in future research, and we are very grateful to the anonymous reviewers for their contributions to this article.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the Zhejiang Provincial Key R & D Program (No.2020C01084), the Zhejiang Provincial Key R & D Program (No.2021C01071), and the Longgang Institute of Zhejiang Sci-Tech University Project (LGYJY2021004).

Conflict of interest

The authors declare that there are no conflicts of interest.

Appendix

Basic parameters

A: the number of aisles in the CBS /RS

M: the number of tiers in the CBS /RS

C: the number of columns in the CBS /RS

N_S: the number of cargo input in a batch of orders

N_R: the number of cargo output in a batch of orders

i, j: the locations in the racks

n_SC: the amount of time consumed in the single command cycle

n_DC: the amount of time consumed in the double command cycle

Ws: the width of a storage location

He: the height of a tier

Vx: horizontal velocity of the crane

Vxmax: max horizontal velocity of the crane

a_x: horizontal acceleration/deceleration of the crane

Vy: vertical velocity of the crane

Vymax: max vertical velocity of the crane

a_y: vertical acceleration/deceleration of the crane

G: acceleration of gravity

DH_ij: horizontal distance between i and j

DV_ij: vertical distance between i and j

k_r: coefficient of friction (0.01)

e: energy conversion efficiency

m_s: quality of the items and pallets

m_c: overall quality of a crane (including: upper crossbeam、lower crossbeam、left column、right column、carriage platform、forks、vertical travel motor、horizontal travel motor、electrical control cabinet)

m_v: quality of the lifting device (including: carriage platform, forks)

Auxiliary variables

t_Hij: horizontal running time between i and j

t_Vij: vertical running time between i and j

t_Haij: horizontal acceleration time between i and j

t_Hdij: horizontal deceleration time between i and j

t_Vcij: horizontal constant speed time between i and j

t_Vaij: vertical acceleration time between i and j

t_Vdij: vertical deceleration time between i and j

t_Vcij: vertical constant speed time between i and j

t_SC: single command cycle time

t_DC: double command cycle time

T_a: the running time of a crane

T_F: order fulfillment time

F_Hc: horizontal traction force of a crane moving horizontally at a constant speed

F_Ha: horizontal traction force of a crane moving horizontally at acceleration

F_Hd: horizontal traction force of a crane moving horizontally at deceleration

F_Vc: vertical traction force of a crane moving horizontally at a constant speed

F_Va: vertical traction force of a crane moving horizontally at acceleration

F_Vd: vertical traction force of a crane moving horizontally at deceleration

P_Hc: horizontal power of a crane moving horizontally at a constant speed

P_Ha: horizontal power of a crane moving horizontally at acceleration

P_Hd: horizontal power of a crane moving horizontally at deceleration

P_Vc: vertical power of a crane moving vertically at a constant speed

P_Va: vertical power of a crane moving vertically at acceleration

P_Vd: vertical power of a crane moving vertically at deceleration

EC_HFij: horizontal energy consumption of a crane with full load moving between i and j

EC_VFij: vertical energy consumption of a crane with full load moving between i and j

EC_HEij: horizontal energy consumption of a crane with empty load moving between i and j

EC_VEij: vertical energy consumption of a crane with empty load moving between i and j

EC_SC: energy consumption of the single command cycle pattern

EC_DC: energy consumption of the double command cycle pattern

ECa: energy consumption of a crane

EC: energy consumption of a batch of orders

Decision variables

N: Rack number

(identifier, N): Crane operation instructions (When the identifier is 1, it means that the item is stored in a location where the rack number is N; When the identifier is 0, it means that the item is retrieved from the rack numbered N)