A new AGV scheduling algorithm based on harmony search for material transfer in a real-world manufacturing system

Abstract

With the development of advanced manufacturing technology, more and more enterprises utilize automated guided vehicles to transfer materials for computer numerical control machines in the workshop. Because automated guided vehicle scheduling problem is non-deterministic polynomial–hard, highly efficient automated guided vehicle scheduling algorithms are necessary to improve productivity of the workshop. In order to solve this problem, a mathematical model with a new objective function is established firstly, which incorporates two indicators, that is, the standard deviation of the waiting time of computer numerical control material buffers and the total travel distance of automated guided vehicle, reflecting the overall capacity of computer numerical control machines and the energy efficiency of automated guided vehicles during the actual production, respectively. Then, an improved harmony search algorithm is proposed, which includes an effective discrete encoding scheme of harmony, a new initialization method for harmony memory based on opposition-based learning strategy, a dynamic harmony memory considering rate parameter and a local search strategy. The property of the proposed algorithm is evaluated through three cases from the real-world manufacturing system for producing back cover of smart phone. Compared with the competing algorithms, the results demonstrate the superiority of proposed algorithm for solving automated guided vehicle scheduling problem obviously.

Introduction

Automated guided vehicles (AGVs) are widely used to transport commodities and materials between workstations and warehouse. ¹ AGVs are among the most efficient methods for material handling due to their supreme routing flexibility, space utilization, product quality, and safety. The delivery demands of computer numerical control (CNC) machines are unpredictable and highly dependent on the product. In general, an optimized transfer scheduling can effectively decrease the waiting time of CNC machines and improve the production efficiency. Thus, optimally addressing the problem of AGV scheduling for transferring materials to CNC machines (AGVSP) is instrumental to developing an enterprise into a manufacturing giant.

At present, for almost all the enterprises, the transfer AGV scheduling for CNC machines is based on “First Come First Served (FCFS),” that is, transfer materials to CNC machines according to their sequence of calling materials. However, this scheduling method may not be the optimal solution; sometimes, it tends to seriously compromise the production efficiency of enterprises. The count of feasible scheduling solutions increases exponentially with the number of CNC machines, so the AGVSP is a non-deterministic polynomial(NP)-hard problem. Therefore, it is practically impossible to use a brute-force approach for identifying the optimal solution to this problem. However, it is feasible to solve this problem with meta-heuristic algorithms.

To achieve a significant improvement concerning the overall performance of AGVSP, the scheduling process must consider the interdependencies between CNC machines and AGVs. In this article, we research on a part of the real complex workshop where there has been single AGV feed for cells and an optimal methodology for solving AGVSP is proposed. Due to the discreteness of AGVSP, harmony search (HS) algorithm is one of the most appropriate meta-heuristics for settling the AGVSP. By a judicious AGV scheduling from three cases in a real-world manufacturing system for producing back cover of smart phone, the standard deviation of the waiting time of CNC material buffers can be reduced to a great extent, that is, the overall capacity of CNC machines can be improved significantly. Moreover, this effort can also potentially shorten an AGV’s travel routes, that is, improve the energy efficiency of AGV.

Related works

Dispatching AGVs to transport the assigned jobs to destinations is a kind of AGV scheduling. An appropriate route should be searched when the AGV scheduling is planned. The AGV scheduling problems belong to a kind of vehicle routing problem (VRP) or sequencing problem.^2,3 The research of the scheduling problem by French ⁴ is classic. They not only introduced the job-shop scheduling but also used material to illustrate the ideas behind combinatorial optimization.

A variety of literatures researched on AGV scheduling problems for different applications. Soylu et al. ⁵ developed an artificial neural network algorithm to seek the shortest routing for the single, multi-capacitated AGV carrying out pick and deliver requests. Tsirimpas et al. ⁶ proposed dynamic programming algorithms and illustrated three practical variations to find the optimal routing of the single, limited capacity vehicle transporting the assigned product to the appointed client. Minis and Tatarakis ⁷ proposed a dynamic programming algorithm to minimize the routing cost of the single vehicle delivering random products and picking up items. Hu et al. ⁸ combined the static problem model and solution approach with a rolling-horizon approach to research the dynamic routing problem. In the matter of workshop layout, the multi-capacitated Rail Guided Vehicle (RGV) working on a linear track and the two-sided loading/unloading operations are similar to ours. There is also a resemblance between the objective of minimizing the energy consumption of RGV and the total travel distance of AGV. However, they regarded the routing problem as a mixed-integer linear program. Qudeiri et al. ⁹ proposed an artificial neural network approach to predict the products’ routes’ group and then combined with genetic algorithms (GAs) to optimize proper routes for all product types in flexible manufacturing systems. Gultekin et al. ¹⁰ developed a mathematical model for the scheduling problem of flexible robotic cell considering fixed process time and sequence for the CNC machines. They regarded the problem as a special traveling salesman problem. Nejad et al. ¹¹ developed a simulated annealing algorithm to study the problem of transferring materials to machines by robot, which is similar to this study. They considered the scheduling problem of a flexible robotic cell in which machines were identical and parallel, meanwhile the sequence of the actions was also determined, which was regarded as the traveling salesman problem. The objective is to minimize the cycle time of the robot.

The above literatures research on AGV scheduling problem, but they are single objective optimization problems; however, the single-objective scheduling model cannot fully address the requirements of the real-world production system. To the best of our knowledge, few researches have been performed to solve the AGVSP using HS algorithm. In this article, we investigate a multi-objective AGVSP problem derived from a realistic manufacturing case and then propose an Improved Harmony Search–based Material Transfer Scheduling Algorithm (IHSMTSA) to solve AGVSP.

The remainder of this article is organized as follows. The related works are elucidated in section “Related works.” Section “Mathematical model of AGVSP” describes the mathematical model of AGVSP. The proposed material transfer algorithm with sufficient detail is presented in section “Proposed algorithm of AGVSP.” In section “Experimental results and analysis,” an empirical evaluation through computational experiments is conducted. Section “Conclusion” presents the conclusions and outlines the potential future direction of this research.

Mathematical model of AGVSP

Problem description

The production line in the workshop illustrated in Figure 1 is popular at present. The workshop is composed of multiple CNC machine cells, and each cell contains multiple CNC machines. The production and processing operation is carried out in the form of CNC machine cells and each cell sets a buffer. Materials needed for CNC machine cells will be delivered by AGV. Then, when CNC machines are brought into production operations, materials will be obtained from the buffer. In this article, we research on one of the production lines in the workshop shown in Figure 1, which is briefly depicted in Figure 2. Here, we only study the optimization problem of AGV transferring materials to the buffer of each CNC machine cell, which aims to reduce the standard deviation of the waiting time of CNC material buffers and the total travel distance of AGV. With the purpose of avoiding wide variance between cells, the standard deviation of waiting time of CNC material buffers is considered. The risk of starved feeding will increase by the wide variance, which will influence capacity equalization of cells and then reduce overall productivity.

Figure 1.

Schematic diagram of a CNC production line with cells and buffers.

Figure 2.

Schematic diagram of CNC material buffers calling material.

Assuming a CNC production line is shown in Figure 2, and the CNC buffers’ sequence of calling material is (1, 10, 11, 19, 2, 20), if an AGV feeds CNC material buffers in this order, it would need to travel repeatedly from one end of the track to the other. Consequently, it is clear that this transfer scheduling approach will make the total waiting time of the buffers greatly extended because the AGV spends too much time on traveling. More seriously, if a CNC material buffer waits too long, some CNC machines in this cell may run but lack of materials will reduce the productivity greatly. According to the specific situation of material calling depicted in Figure 2, a reasonable solution should make AGV feed the CNC material buffers at one end of the track followed by moving to the other end of the track to feed other CNC material buffers, a process being able to greatly reduce the total waiting time of the CNC material buffers so as to improve the productivity. However, during the actual machining, the CNC buffers’ sequence of calling material is dynamic and unpredictable. Therefore, based on the above introduction and marked signs, the problem that we will research in this article can be described as how to arrange scientifically the optimal sequence of tasks of having sent out demand signals for materials within record time slice (RTS) applied an effective modern heuristic algorithm and then produce the delivery sequence of AGV. RTS stands for a period of time. In other words, tasks of sending out demand signals for materials by buffers during this time period are researched in this article, and for the sake of research, RTS is assigned as 1000 in which there has been enough cells that sent out demand signals for materials. However, it must be specified that the method presented and validated in this article can also be applied to other situations when RTS is assigned other values.

IHSMTSA optimization model

As mentioned above, the standard deviation of the waiting time of CNC material buffers and the total travel distance of the AGV are the two indicators should be considered when solving AGVSP. Therefore, we propose the mathematical model as follows

$$\mathrm{Min}f(X)=w_1{f}_1+w_2{f}_2$$

(1)

$$f_1=\sqrt{\frac{1}{n}\sum _{i=1}^n{(\mathrm{wT}(x_i)-w{T}_{\mathrm{aver}})}^2}$$

(2)

$$f_2=\sum _{i=1}^n\mathrm{abs}(\mathrm{pos}(x_i)-\mathrm{pos}(x_{i-1}))$$

(3)

$$\mathrm{wT}\left(x_i\right)=\{\begin{array}{c}\mathrm{RTS}-\mathrm{calltime}\left(x_1\right)+\mathrm{rT}+\frac{\mathrm{abs}(\mathrm{pos}\left(x_1\right)-\mathrm{pos}\left(x_0\right))}{V_{\mathrm{agv}}}+T_{\mathrm{sl}}\mathrm{if}i=1\hfill \\ \mathrm{calltime}\left(x_{i-1}\right)-\mathrm{calltime}\left(x_i\right)+\mathrm{wT}\left(x_{i-1}\right)+\frac{\mathrm{abs}(\mathrm{pos}\left(x_i\right)-\mathrm{pos}\left(x_{i-1}\right))}{V_{\mathrm{agv}}}+T_{\mathrm{sl}}\mathrm{otherwise}\hfill \end{array}$$

(4)

$$w{T}_{\mathrm{aver}}=\frac{1}{n}\sum _{i=1}^n\mathrm{wT}(x_i)$$

(5)

$$\mathrm{pos}(x_0)=\mathrm{pos}(\mathrm{AGV})$$

(6)

The objective function (1) is to minimize the standard deviation of the waiting time of CNC material buffers and the total travel distance of the AGV. $X=(x_1,x_2,\dots ,x_i,\dots ,x_n)$ denotes a solution, which represents a sequence of CNC material buffers fed by AGV; $f_1$ denotes the standard deviation of the waiting time of CNC material buffers corresponding to the solution $X$ ; $f_2$ represents the total travel distance of the AGV; and $w_1$ and $w_2$ are the weights of $f_1$ and $f_2$ , respectively. The values of $f_1$ and $f_2$ are calculated in equations (2) and (3), respectively. Equation (4) presents the waiting time of CNC material buffers calling material. Where, $\mathrm{calltime}(x_i)$ denotes the time of the buffer $x_i$ calling material, $\mathrm{pos}(x_i)$ indicates the position of the buffer $x_i$ , and $\mathrm{pos}(\mathrm{AGV})$ indicates the current location of the AGV. It should be noted that the location of the AGV needs to be updated in accordance with the transfer operation. It can be seen from equation (4) that the waiting time of the current buffer consists of five parts: The first part denotes the duration from the calling time to the end of RTS, the second part corresponds to the algorithm running time, the third part is the waiting time of the former buffer calling material, the fourth part is the time it takes for AGV to travel to the current buffer, and the fifth part is the loading time. Equation (5) is the average waiting time of CNC material buffers calling material. $w{T}_{\mathrm{aver}}$ is the average waiting time of CNC material buffers calling material. Equation (6) presents that $x_0$ is the initial location of AGV.

The conditions associated with the AGVs-CNC system are specified in the following:

1. The distance between the buffer centers of two neighboring CNC cells on the same track (m): 5.5.

2. The initial position of AGV (m): $\mathrm{iniPos}$ .

3. The average moving speed of AGV (m/s): V_agv = 0.45.

4. Time needed for conducting a single loading of AGV to a buffer (s): T_sl = 30.

5. The sum of the algorithm running time (s): rT.

6. RTS = 1000 s.

Proposed algorithm of AGVSP

Brief introduction of HS algorithm

As a new meta-heuristics search algorithm, the HS algorithm ¹² is impressed by the searching process of a perfect harmony for music. HS employs a suite of newness stochastic derivatives employed to the discrete variables that leverages a musician’s experience to identify the searching direction so as to ensure the convergence of the algorithm. It does not need differential gradients and initial guesses to be applied to the discrete and continuous variables. The HS algorithm can also extract the optimal solution when it is in the local optima. It is able to overcome the shortcoming associated with the GA’s theory in building module and take the relationship into consideration definitely in an integration operation.

HS algorithm is increasingly used for optimization of problem-solving schemes in many fields, such as a vehicle dispatching problem in the NC workshop which is equipped with the just-in-time (JIT) distribution, ¹³ a flexible job-shop scheduling problem, ¹⁴ the comprehensive production and transport scheduling problem associated with the make-to-order (MTO) manufacturing, ¹⁵ the routing problem of wireless sensor network (WSN), ¹⁶ the optimization of continuous problems, ¹⁷ realization of two-dimensional (2D) nearest neighbor quantum circuits, ¹⁸ discrete sizing optimization of truss structures, ¹⁹ fine-tuning deep belief networks, ²⁰ the clustering problem of WSN, ²¹ image reconstruction from projections, ²² and single area unit commitment problem. ²³

HS first initializes the harmony memory, followed by populating the memory with randomly generated new harmony. Only when the newly generated harmony is better than the worst harmony in the original harmony memory, the new harmony is put in the memory and replaced with the worst harmony of the original memory. This replacement process repeats itself until the stop criterion is met. The detailed workflow of HS is outlined in the following: ²⁴

• Step 1. Initialize the algorithm parameters and specify the optimization problem

The size of the harmony memory size (HMS) refers to the number of solution vectors in HM and is usually set to 5. The algorithm parameters like the harmony memory considering rate (HMCR), pitch adjusting rate (PAR), bandwidth (bw), and the stopping criterion (maximum iteration number) are applied to explore better solution vectors and are all defined in Step 3.

The optimization problem is defined as

$$\mathrm{Minimize}f\left(X\right)\mathrm{subject}\mathrm{to}{x}_i\in X_i,i=1,2,\dots ,n$$

(7)

where $X$ denotes the set of designed variables; $f(X)$ denotes the objective function; $X^l$ denotes the set of possible values of $x_i$ , which, in the case of continuous optimization, follows $X^i=\left\{x_i|x_i^L\le x_i\le x_i^U\right\}$ with $x_i^L$ being the lower bounds and $x_i^U$ being the upper bounds of $x_i$ respectively; and $n$ denotes the number of designed variables.

• Step 2. Initialize the harmony memory HM with random values drawn from vectors $[X^L,X^U]$

The HM, as shown in equation (8), is a memory space storing solution vectors, that is, harmony, which is generated randomly following the uniform distribution

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{HM}=\left[\begin{array}{c}{X}_1\hfill \\ \multicolumn{1}{c}{X_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{X_i}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{X_{\mathrm{HMS}}}\end{array}\right]}\\ =\left[\begin{array}{cccccc}{x}_{1,1}\hfill & x_{1,2}\hfill & \cdots \hfill & x_{1,j}\hfill & \cdots \hfill & x_{1,n}\hfill \\ \multicolumn{1}{c}{x_{2,1}}& x_{2,2}\hfill & \cdots \hfill & x_{2,j}\hfill & \cdots \hfill & x_{2,n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{x_{i,1}}& x_{i,2}\hfill & \cdots \hfill & x_{i,j}\hfill & \cdots \hfill & x_{i,n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{x_{\mathrm{HMS},1}}& x_{\mathrm{HMS},2}\hfill & \cdots \hfill & x_{\mathrm{HMS},j}\hfill & \cdots \hfill & x_{\mathrm{HMS},n}\hfill \end{array}\right]\hfill \end{array}$$

(8)

• Step 3. Improvise a new harmony vector

A new harmony $X\prime =({x\prime }_1,{x\prime }_2,\dots ,{x\prime }_i,\dots ,{x\prime }_n)$ is produced following two steps, that is, harmony memory considerations and pitch adjustments. The first step is shown as follows

$$x{\prime }_i\leftarrow \{\begin{array}{cc}{x\prime }_i\in \{x_{1,i},x_{2,i},\dots ,x_{\mathrm{HMS},i}\}\hfill & \mathrm{if}{P}_1<\mathrm{HMCR}\hfill \\ \multicolumn{1}{c}{{x\prime }_i\in X_i}& \mathrm{else}\hfill \end{array}$$

(9)

where HMCR is a constant between 0 and 1, and $P_1$ is a random number generated from the linear space spanning $[0,1]$ . Equation (9) indicates that $x{\prime }_i$ should be selected from the ith column of $\mathrm{HM}$ , that is, $\{x_{1,i},x_{2,i},\dots ,x_{\mathrm{HMS},i}\}$ , if $P_1$ is smaller than HMCR; otherwise, $x{\prime }_i$ can be a feasible value belonging to $X^l$ .

We then need to determine whether or not to conduct pitch adjustment on these new harmonies selected from the HM. The criterion is shown as follows

$$\mathrm{Pitch}\mathrm{adjusting}\mathrm{decision}\mathrm{for}x{\prime }_i\leftarrow \{\begin{array}{c}\mathrm{Yes}\hfill \\ \multicolumn{1}{c}{\mathrm{No}}\end{array}\begin{array}{c}\mathrm{if}(P_2<\mathrm{PAR})\hfill \\ \multicolumn{1}{c}{\mathrm{else}}\end{array}$$

(10)

where PAR is a constant between 0 and 1 and $P_2$ is a random number between 0 and 1. The formula above indicates that if PAR is greater than the random number $P_2$ , the new harmony $x\prime$ must be pitch adjusted; otherwise, the original stays. In the case where a pitch adjustment is needed, $x\prime$ would be subject to the following process

$$x{\prime }_i\leftarrow x{\prime }_i\pm \mathrm{rand}\ast \mathrm{bw}$$

(11)

where bw denotes the arbitrary distance bandwidth and $\mathrm{rand}$ is a random number generated between 0 and 1 following uniform distribution.

• Step 4. If the newly improvised harmony ${\overline{x}}_i$ is superior to the worst harmony in HM, the former would replace the latter.

• Step 5. Terminate if the current iteration number exceeds $\mathrm{MAXIMP}$ . Otherwise, repeat Steps 3 and 4. At last, the optimal solution is derived from the best harmony in HM.

The proposed IHSMTSA for AGVSP

The IHSMTSA for solving AGVSP is proposed in this section. Figure 3 describes the specific flow of IHSMTSA, which shows that the key steps of IHSMTSA include the initialization of HM, the process of the improvisation for a new harmony, and the local search. This section describes the encoding of harmony first and subsequently illustrates the initialization of HM and the improvisation of a new harmony. Finally, the local search scheme is presented.

Figure 3.

Flowchart of IHSMTSA.

Encoding of harmony

Considering that AGVSP is a discrete optimization problem, it is more appropriate to use discrete coding method, which can greatly narrow down the search range in the solution space. We adopt the following vector for representing the ith harmony in the HM

$$X_i=(x_{i,1},x_{i,2},\dots ,x_{i,j},\dots ,x_{i,n})$$

(12)

where $X_i$ denotes the ith harmony of HM, n denotes the length of the harmony, that is, the total number of CNC machines which need to be fed, and x_i,j represents the jth pitch of harmony $X_i$ .

Since the data associated with the material calling differ from time slice to time slice, the following data structure is defined to facilitate the encoding of algorithm (Figure 4).

Figure 4.

Structure type.

A callMaterial type array named call is defined to record the information pertaining to the material called by CNC material buffers. The serial number of each CNC material buffer calling material and the corresponding time at which the material is called are saved in the array call following the order in which the materials are called, as shown in Figure 5.

Figure 5.

Structure of the array call.

Supposing the CNC buffers’ sequence of calling material and the corresponding time at which the material is called associated within the RTS are stored in the array call as shown in Figure 5. The tone in a harmony corresponds to the subscript number of the elements in the array call. For example, the ith harmony of HM is $X_i=(2,3,5,7,8,4,0,9,6,1)$ , indicating that the order of AGV loading for CNC machines is $(17,11,19,13,10,6,8,14,2,4)$ .

Initialization of harmony

As a new method, opposition-based learning (OBL) can improve the performance of machine learning at a large extent, which is proposed by Tizhoosh. ²⁵ OBL increases the probability of seeking the global optimal solution by extending the domain of the initial solution in the searching area. Building the initial HM by OBL can therefore assist in the identification of a more appropriate initial solution, which is illustrated in a previous study.

Supposing $P=(x_1,x_2,\dots ,x_i,\dots x_D)$ is a D-dimensional solution vector, and $x_i=[x_i^L,x_i^U]$ , $i=1,2,\cdots ,D$ , the opposing solution vector $\mathrm{OP}=(x{\prime }_1,x{\prime }_2,\dots ,x{\prime }_i,\dots ,x{\prime }_D)$ is defined as follows

$$x{\prime }_i=x_i^L+x_i^U-x_i$$

(13)

Furthermore, the harmony i is defined as $X_i(x_{i,1},x_{i,2},\dots ,x_{i,j},\dots ,x_{i,n})$ , whose opposing solution vector is $O{X}_i=(x{\prime }_{i,1},x{\prime }_{i,2},\dots ,x{\prime }_{i,j},\dots ,x{\prime }_{i,n})$ . The element $x{\prime }_{i,j}$ in the opposing solution vector can be computed from equation (14)

$$x{\prime }_{i,j}=\mathrm{call}.\mathrm{size}-1-x_{i,j}$$

(14)

where $\mathrm{call}.\mathrm{size}$ denotes the number of elements in the array call.

Improvisation of a new harmony

To apply the original classic HS optimization algorithm for solving a continuous optimization problem, one can use the Tone Tuning Bandwidth $\mathrm{bw}$ as a tuning parameter to adjust the harmonies selected from the HM. The AGVSP, however, is a discrete optimization problem, which may lead to tremendous amount of operational cost for adjusting harmonics based on the traditional HS algorithm. As a result, we choose to skip the step of pitch adjustment while a new harmony is improvised. Consequently, HMCR becomes the sole algorithm-specific parameter in the IHSMTSA. Specifically, the dynamic adjustment of HMCR is conducted to prevent IHSMTSA from falling into the local optima in the early stage, giving rise to an enhanced global search ability. Meanwhile, the ability of local search in the late stage is also reinforced.

An excessively small value of HMCR may compromise the local search ability of the HS algorithm, making it difficult to converge to the global optimal solution; an excessively large value of HMCR, however, may lead to an overly strong local search ability, rendering the algorithm prone to land on the local optimal solution instead of converging to the global optimal solution. Therefore, in this article, the HS algorithm conducts a dynamic update to the HMCR parameter. In the early stage, HMCR takes on a relatively small value to ensure a good global search ability. In the end, HMCR is assigned a relatively large value to improve the local search ability, which facilitates the convergence to the global optimal solution.

The relationship between HMCR and the number of iterations is shown in equation (15)

$$\mathrm{HMCR}\left(\mathrm{iter}\right)=\{\begin{array}{cc}\frac{3\ast ({\mathrm{HMCR}}_{max}-{\mathrm{HMCR}}_{\min })}{\mathrm{ite}{r}_{max}}\ast \mathrm{iter}+{\mathrm{HMCR}}_{\min },\hfill & \mathrm{if}(\mathrm{iter}<\frac{\mathrm{ite}{r}_{max}}{3})\hfill \\ \multicolumn{1}{c}{{\mathrm{HMCR}}_{\max }}& \mathrm{otherwise}\hfill \end{array}$$

(15)

where $\mathrm{ite}{r}_{max}$ denotes the maximum number of iterations, $\mathrm{iter}$ represents the current iteration number, and ${\mathrm{HMCR}}_{min}$ and ${\mathrm{HMCR}}_{max}$ denote the minimum and maximum values of $\mathrm{HMCR}$ , respectively. In general, the values of ${\mathrm{HMCR}}_{min}$ and ${\mathrm{HMCR}}_{max}$ are set to be 0.2 and 0.99, respectively. The relationship between $\mathrm{HMCR}$ and the iteration number is shown in Figure 6.

Figure 6.

Strategy for dynamically adjusting HMCR.

Supposing $X\prime =(x{\prime }_1,x{\prime }_2,\dots ,x{\prime }_i,\dots ,x{\prime }_n)$ is a new harmony, the element $x{\prime }_i$ of $X\prime$ can be generated using equation (16)

$$x{\prime }_i\leftarrow \{\begin{array}{cc}{x\prime }_i{\to }_{\mathrm{rd}}\{x_{1,i},x_{2,i},\dots ,x_{\mathrm{HMS},i}\}\hfill & \mathrm{if}(P_1<\mathrm{HMCR})\hfill \\ \multicolumn{1}{c}{{x\prime }_i{\to }_{\mathrm{rd}}{X}^i}& \mathrm{Otherwise}\hfill \end{array}$$

(16)

where ${\to }_{\mathrm{rd}}$ denotes the process of random selection. Equation (16) dictates that $x{\prime }_i$ is randomly selected from $\{x_{1,i},x_{2,i},\dots ,x_{\mathrm{HMS},i}\}$ in the case where the random number $P_1$ is smaller than HMCR; otherwise, $x{\prime }_i$ is chosen from $X^i$ randomly.

Local search

As a heuristic algorithm for solving optimization problems, local search is performed within the vicinity of the current solution. Specifically, it searches new solution within the neighborhood of the current solution. When the new solution is superior to the current solution, the former will substitute the latter. In view of above features of the AGVSP, the local search in IHSMTSA is implemented by moving a randomly selected tone to a new position in the harmony, as shown in Figure 7. If a new harmony with an improved fitness value is identified, it will replace the existing one. With this local search strategy, the quality of the optimal solution will be improved greatly.

Figure 7.

Schematic diagram of local search.

Experimental results and analysis

Experimental settings

The experimental results are presented in this section. The algorithms compared here are all implemented using C++ programming language by Microsoft visual studio 2015 and executed on a personal computer with 2.5 GHz Intel Core i7 3610QM processor, 8 GB RAM and Microsoft Windows 10 operating system. In order to verify the superiority of the proposed IHSMTSA, the comparisons are conducted by incorporating the particle swarm optimization (PSO) with inertia weight; GA with elite selection; single-point crossover and two bit-swap mutation strategies; memetic algorithm (MA) with elite selection, single-point crossover, two bit-swap mutation, and local search of inserting a gene into another position; and HS algorithms. The CNC production line used in our experiments in shown in Figure 8, which includes 30 CNC machine cells, that is, 30 material buffers and 1 AGV. To make the comparison comprehensive, three cases from the real-world manufacturing system for producing back cover of smart phone are tested here involving the aforementioned five algorithms.

Figure 8.

Schematic diagram of the CNC production line used in the present experiments.

Case 1

There are 15 CNC material buffers calling material within the RTS, whose information is shown in Table 1.

Table 1.

Case 1.

No.	1	16	28	21	2	8	3	5	15	30	18	19	12	9	6
Time	10	150	160	190	250	280	450	560	650	760	820	880	950	980	990

Case 2

There are 20 CNC material buffers calling material within the RTS, whose information is shown in Table 2.

Table 2.

Case 2.

No.	2	14	30	16	7	19	29	5	20	11
Time	6	50	80	100	150	200	260	280	320	380
No.	1	21	8	26	3	6	17	21	22	15
Time	420	440	480	520	580	660	750	820	860	910

Case 3

There are 25 CNC material buffers calling material within the RTS, whose information is shown in Table 3.

Table 3.

Case 3.

No.	8	19	28	6	12	20	24	2	30	25	13	3	1
Time	10	48	82	110	146	180	212	250	278	325	356	388	421
No.	29	22	14	7	26	10	23	16	5	15	9	18
Time	451	482	522	546	587	622	678	755	787	825	866	970

Parameter settings

The parameter value of AGVs-CNC system and the experiments are listed in Table 4. It is worth noting that the experimental results were obtained by taking an average of 25 independent runs for fairness. The parameters of IHSMTSA and the comparison of algorithms are set as shown in Table 5.

Table 4.

Parameter settings of experiments.

Items	Value
Distance between two buffer centers	5.5 (m)
iniPos	0 (m)
Vagv	0.45 (m/s)
Tsl	30 (s)
rT	10 (s)
RTS	1000 (s)
Independent runs	25
FEs	10,000
w 1	0.7
w 2	0.3

RTS: record time slice.

Table 5.

Parameter settings of algorithms.

Algorithms	Parameters
IHSMTSA	HMS = 6, HMCRmin = 0.2, HMCRmax = 0.99
HS	HMS = 6, HMCR = 0.95, PAR = 0.05
PSO	PopSize = 100, ωmax = 0.9, ωmin = 0.1, c1 = 0.9, c2 = 0.1
GA	PopSize = 100, Pc = 0.9, Pm = 0.01
MA	PopSize = 100, Pc = 0.9, Pm = 0.01

IHSMTSA: Improved Harmony Search based Material Transfer Scheduling Algorithm; HMS: harmony memory size; HMCR: harmony memory considering rate; PAR: pitch adjusting rate; HS: harmony search; PSO: particle swarm optimization; GA: genetic algorithm; MA: memetic algorithm; PopSize: population size; ω_max: upper limit of inertia weight; ω_min: lower limit of inertia weight; c₁, c₂: acceleration factors; P_c: crossover probability; P_m: mutation probability.

Results and analysis

To test the efficiency of the present algorithm, the IHSMTSA is compared with other meta-heuristic-based transfer scheduling algorithms in terms of the convergence rate for the cases 1–3, respectively. Figures 9–11 show the convergence behaviors of all the algorithms, plotting the average fitness value and the optimal fitness value as a function of evaluation number for each algorithm. Tables 6–8 show the standard deviation of the waiting time of CNC material buffers, the average waiting time of CNC material buffers, and the total travel distance of AGV.

Figure 9.

Convergence behaviors of various algorithms in case 1: (a) average fitness value versus FEs and (b) best fitness value versus FEs.

Figure 10.

Convergence behaviors of various algorithms in case 2: (a) average fitness value versus FEs and (b) best fitness value versus FEs.

Figure 11.

Convergence behaviors of various algorithms in case 3: (a) average fitness value versus FEs and (b) best fitness value versus FEs.

Table 6.

Comparison of algorithms in case 1.

Methods	Standard deviation of the waiting time of CNC material buffers	Average waiting time of CNC material buffers (s)	Total travel distance of AGV (m)
IHSMTSA	149.712	1146.89	313.06
HS	156.735	1155.43	316.58
PSO	228.282	1155.56	346.28
MA	154.56	1206.63	317.46
GA	150.959	1174.07	326.26
FCFS	150.773	1404.37	555.5

AGV: automated guided vehicle; CNC: computer numerical control; IHSMTSA: Improved Harmony Search based Material Transfer Scheduling Algorithm; HS: harmony search; PSO: particle swarm optimization; GA: genetic algorithm; MA: memetic algorithm; FCFS: First Come First Served.

Note: The significance of the values given in boldface is reflected in the experimental analysis paragraph of case 1.

Table 7.

Comparison of algorithms in case 2.

Methods	Standard deviation of the waiting time of CNC material buffers	Average waiting time of CNC material buffers (s)	Total travel distance of AGV (m)
IHSMTSA	161.959	1061.03	138.38
HS	186.565	1143.38	180.18
PSO	326.832	1292.59	314.38
MA	215.986	1179.49	208.34
GA	194.398	1158.83	216.7
FCFS	256.223	1596.76	594

AGV: automated guided vehicle; CNC: computer numerical control; IHSMTSA: Improved Harmony Search based Material Transfer Scheduling Algorithm; HS: harmony search; PSO: particle swarm optimization; GA: genetic algorithm; MA: memetic algorithm; FCFS: First Come First Served.

Note: The significance of the values given in boldface is reflected in the experimental analysis paragraph of case 2.

Table 8.

Comparison of algorithms in case 3.

Methods	Standard deviation of the waiting time of CNC material buffers	Average waiting time of CNC material buffers (s)	Total travel distance of AGV (m)
IHSMTSA	160.854	1277.79	223.08
HS	228.171	1323.59	253
PSO	439.593	1491.77	436.7
MA	250.121	1325.58	284.24
GA	236.749	1308.04	272.58
FCFS	388.008	1780.56	687.5

AGV: automated guided vehicle; CNC: computer numerical control; IHSMTSA: Improved Harmony Search based Material Transfer Scheduling Algorithm; HS: harmony search; PSO: particle swarm optimization; GA: genetic algorithm; MA: memetic algorithm; FCFS: First Come First Served.

Note: The significance of the values given in boldface is reflected in the experimental analysis paragraph of case 3.

Case 1: the experimental results of 15 buffers calling material

Figure 9(a) shows the average fitness value versus FEs associated with five algorithms in case 1. The solution found by IHSMTSA is much better than those found by other algorithms including PSO, MA, and GA. Although the HS converges faster than IHSMTSA in early iterations as HS tends to get trapped in local optima, the HS converges slower than IHSMTSA in late iterations. Compared to the rate of convergence, it is of more importance to examine the quality of the solution. We can get the summary that IHSMTSA yields the highest probability of searching the best solution among the algorithms assessed here. Figure 9(b) shows the convergence curve of the best fitness value. The first population fitness value of IHSMTSA is worse than those of other algorithms due to its limited population capacity which makes it yield a lower probability of obtaining good solutions. However, the IHSMTSA can eventually converge to a more optimal solution than other algorithms after an increased number of iterations due to the fact that the dynamic changed HMCR can avoid falling into the local optima in early iterations, giving rise to an enhanced global search ability, and improve the ability of local search in later iterations. Besides, the proposed local search strategy has a great contribution on improving the quality of the optimal solution. Based on case 1 results, it can be concluded that IHSMTSA is superior to other four algorithms in terms of efficiency and accuracy.

In this case of AGVSP, the comparison involving IHSMTSA, HS, PSO, MA, GA algorithms and FCFS has also been conducted. Table 6 illustrates the statistical results associated with the five algorithms and FCFS in case 1. The data used in FCFS are collected from engineering project in order to examine the science and validity of the proposed algorithm. Meanwhile, the data used in other cases are also credible. As we can see from Table 6, the proposed IHSMTSA performs best in terms of all the indicators, that is, the standard deviation of the waiting time of CNC material buffers, the average waiting time of CNC material buffers, and the total travel distance of AGV. Especially, the average waiting time of CNC material buffers obtained by IHSMTSA has been decreased by 18.33% compared with those obtained by FCFS, which means that IHSMTSA performs far better than FCFS in terms of the overall capacity. Moreover, the total travel distance of AGV obtained by IHSMTSA has been decreased by 43.64% compared with those obtained by FCFS, which means that IHSMTSA can improve the energy efficiency of AGV significantly.

Case 2: the experimental results of 20 buffers calling material

As we can see from Figure 10, IHSMTSA performs far better than other algorithms compared in terms of the rate of convergence and the accuracy of the best solution, similar to case 1. It can be seen from Table 7 that the standard deviation of the waiting time of CNC material buffers, the average waiting time of CNC material buffers, and the total travel distance of AGV associated with IHSMTSA are all shorter than those of the other four algorithms and FCFS. Especially, the values of the three indicators obtained by IHSMTSA have been decreased by 36.79%, 33.55%, and 76.7% compared with those obtained by FCFS, indicating that the advantages of IHSMTSA can stand out when solving larger scale AGVSP.

Case 3: the experimental results of 25 buffers calling material

As we can see from Figure 11, IHSMTSA also performs far better than other algorithms compared in terms of the rate of convergence and the accuracy of the best solution, similar to case 1 and case 2, which further confirms the advantages of IHSMTSA on the rate of convergence and accuracy. It can be seen from Table 8 that the values of the three indicators associated with IHSMTSA are far smaller than those of the other four algorithms and FCFS. Especially, the values of the three indicators obtained by IHSMTSA have been decreased by 58.54%, 28.24%, and 67.55% compared with those obtained by FCFS, which further the advantages of IHSMTSA on solving larger scale AGVSP.

By running the program 25 times in three cases, separately, the optimal sequences of CNC material buffers fed by AGV can be obtained by applying the method of IHSMTSA, as shown in Table 9.

Table 9.

The optimization results.

Item	Best sequences
Case 1	8, 14, 30, 1, 6, 17, 12, 2, 19, 15, 16, 13, 28, 29, 18
Case 2	2, 16, 19, 7, 20, 30, 5, 8, 14, 1, 26, 29, 21, 11, 6, 17, 3, 15, 22, 21
Case 3	28, 20, 6, 19, 2, 24, 8, 1, 13, 12, 22, 30, 25, 26, 3, 29, 14, 16, 23, 7, 18, 10, 5, 9, 15

From the three cases mentioned above, it is obvious that the result obtained by IHSMTSA is much better than “FCFS” in real-world production. And the proposed IHSMTSA has better performance than other algorithms compared in terms of both the quality of solutions and the rate of convergence, when solving the AGVSP described in this article. Especially, the advantages of IHSMTSA will stand out when dealing with larger scale AGVSP such as cases 2 and 3. Moreover, the performance of IHSMTSA is much better than other algorithms and FCFS in terms of the three indicators, that is, the standard deviation of the waiting time of CNC material buffers, the average waiting time of CNC material buffers, and the total travel distance of AGV, which means that IHSMTSA has a better performance than other algorithms and FCFS compared to increasing the overall capacity of CNC machines, as well as improving the energy efficiency of AGVs. Therefore, the IHSMTSA algorithm can greatly improve the profitability of enterprises.

Conclusion

In this article, we first proposed an effective objective function model for AGVSP, which incorporates two indicators, that is, the standard deviation of the waiting time of CNC material buffers and the total travel distance of AGV, reflecting the overall capacity of CNC machines and the energy efficiency of AGVs during the actual production, respectively. Then, a novel meta-heuristic method based on HS algorithm, that is, IHSMTSA, is applied to solve the AGVSP.

The current design of IHSMTSA incorporates four major improvements to the traditional HS algorithm: (1) Considering the characteristics of the problem of interest here, the discrete coding method is adopted which can greatly reduce the search range in the solution space. (2) The OBL approach is adopted here to initialize the HM, which can expand the initial solution space in the searching area, thereby enhancing the probability of landing on the global optimum. (3) The parameter HMCR is dynamically changed, which can prevent IHSMTSA from falling into the local optima in early iterations and can enhance the local search ability in later iterations. (4) A local search strategy is recommended to improve the quality of the optimal solution.

As demonstrated in the simulation study, our IHSMTSA is notably superior to other algorithms in terms of the performance in solving AGVSP.

The contributions of this article are summarized as follows:

1. An improved HS algorithm, IHSMTSA, is developed to address the AGVSP, which incorporates a multitude of major improvements to the traditional HS algorithm.

2. The standard deviation of the waiting time of CNC material buffers and the total travel distance of AGV are reduced greatly. The workshop capacity is also increased. As a result, the production efficiency is increased and the energy consumption is reduced to a great extent.

On a going-forward basis, attempts will be made to implement the IHSMTSA in a real-world intelligent production workshop so as to further assess its effectiveness. Then, this algorithm can be applied to more production workshop. On this basis, we can continue to study and improve the algorithm so that the algorithm can be applied to more areas.