# A Hybrid Quantum Evolutionary Algorithm with Improved Decoding Scheme for a Robotic Flow Shop Scheduling Problem

Authors: Weidong Lei; Hervé Manier; Marié-Ange Manier; Xinping Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We aim at solving the cyclic scheduling problem with a single robot and flexible processing times in a robotic flow shop, which is a well-known optimization problem in advanced manufacturing systems. The objective of the problem is to find an optimal robot move sequence such that the throughput rate is maximized. We propose a hybrid algorithm based on the Quantum-Inspired Evolutionary Algorithm (QEA) and genetic operators for solving the problem. The algorithm integrates three different decoding strategies to convert quantum individuals into robot move sequences. The Q-gate is applied to update the states of Q-bits in each individual. Besides, crossover and mutation operators with adaptive probabilities are used to increase the population diversity. A repairing procedure is proposed to deal with infeasible individuals. Comparison results on both benchmark and randomly generated instances demonstrate that the proposed algorithm is more effective in solving the studied problem in terms of solution quality and computational time.


---

## Body

## 1. Introduction

Since the 1970s, with the development of the robotics and automation technologies, computer-controlled robots instead of workers have been gradually used in many manufacturing industries to perform high frequency or dangerous transportation jobs. The advantages of robotic or automated manufacturing systems include higher production rate, better product quality, and more efficient use of materials, improved safety, and reduced labor intensity. Besides, highly robotic or automated manufacturing systems can effectively meet the requirement of mass production and respond to global competition.In automated manufacturing systems, the production line with transporting device usually consists of several processing machines arranged in a line and one or more robots for transporting parts among machines, as shown in Figure1. Due to the industrial applications [1], the part processing time on each machine is usually flexible and imposed with a pair of minimum and maximum time bounds. Besides, the cyclic production mode is usually adopted in such systems because of easy implementation in a mass production environment. This leads to a repetitive schedule performed by the robot in every certain time. The duration of the robot to perform the repetitive schedule is called the cycle time or cycle length [2].Figure 1
A typical robotic production line.For such problem, scheduling of robot movements is very critical because the line productivity and the product quality extremely depend on it. Therefore, the decision concerns how to sequence the robot movements such that the cycle time is minimized, which is equivalent to maximizing the throughput rate of the production line in steady-state. Note that such scheduling problems in different industrial contexts may have different appellations, such as hoist scheduling problem [3], robotic cell, or robotic flow shop scheduling problem [4, 5].Lei and Wang [6] first proved that the cyclic flow shop scheduling problem with a single robot and flexible processing times is NP-complete. Many researchers have been constantly dedicated to this area and proposed various efficient methods for solving the relevant problems. The first fundamental work in this area was done by Phillips and Unger [7], who formulated a Mixed Integer Programming (MIP) model for solving the single robot cyclic scheduling problem with one degree, which means that each machine is emptied exactly once within a cycle. After that, various exact methods or heuristics were developed for the related problem in the literature [2, 8–10]. Moreover, Liu et al. [11] developed an MIP model for solving the problem with complex lines (i.e., including multifunction stations and parallel stations). Che and Chu [12] developed an efficient branch-and-bound (B&B) algorithm for the same problem as Liu et al. [11].Besides the above studies, researchers have also examined the problem with multiple robots. For instance, Lei and Wang [13] and Zhou and Liu [14] proposed various heuristic algorithms based on zone-partition approach for solving different kinds of multiple robots cyclic scheduling problems. Manier et al. [15] first formulated a mathematical model with disjunctive constraints and then proposed a branch-and-bound like algorithm for finding the optimal schedules for multiple robots. Manier and Lamrous [16] proposed an evolutionary approach to minimize the number of robots on parallel tracks. Moreover, Che and Chu [17] and Leung et al. [18] presented a specific B&B algorithm and an MIP approach for the multiple robots scheduling problem with unidirectional part flow, where the part processing route is the same as the machine arrange sequence. Recently, Elmi and Topaloglu [19] proposed an MIP model for the multiple robots scheduling problem with multidegree.Since a higher degree of cyclic schedule would generally improve the production rate, different MIP approaches were suggested for the single robot cyclic scheduling problem with 2 degrees, which means that each machine is emptied exactly 2 times during a cycle [20]. Feng et al. [21] proposed an MIP approach for the dynamic hoist scheduling problem with multiparts. Che et al. [22] proposed a B&B algorithm for the problem withk-degree cycles. More recently, Elmi and Topaloglu [23] developed a metaheuristic algorithm based on ant colony optimization for solving the robotic flow shop scheduling problem with multidegree. Furthermore, Li and Fung [24] proposed an MIP approach for the single robot multidegree cyclic scheduling with reentrance.In recent years, Quantum-Inspired Evolutionary Algorithm (QEA) has been receiving much attention from researchers. It can be seen as a probability optimization algorithm based on the concepts and principles of quantum computation [25, 26]. Since 1990s, QEA has been successfully applied to solve several well-known optimization problems, such as travelling salesman problem [26], knapsack problems [27, 28], and lane reservation problem [29]. Particularly, various variants of QEA have been proposed for solving different production scheduling problems. For instance, Li and Wang [30] employed QEA to solve the multiobjective flow shop scheduling problem. Zheng and Yamashiro [31] proposed different QEAs for the classic flow shop scheduling problem. Moreover, Gu et al. [32] proposed a parallel QEA for solving the stochastic job shop scheduling problem. Later on, Gu et al. [33] proposed a coevolutionary QEA for the same problem.Due to the NP-completeness of the studied problem, the computation time spent by exact methods usually increases exponentially with its size. Thus, it is a wise choice to use metaheuristics to find sufficiently good schedules within a reasonable time. Because of its several unique advantages, such as better population diversity and very well global search ability, QEA has gained great success in many different optimization problems. Up to now, there is no work reported on using QEA to solve the robotic flow shop scheduling problem. So in this article, we propose a new scheduling algorithm based on QEA and genetic operators for the single robot cyclic scheduling problem with flexible processing times in a robotic flow shop.A solution is usually represented by a permutation of job input sequence in classic flow shop or job shop scheduling problems. However, a chromosome is encoded by Q-bits in QEA, which is then converted into a binary chromosome. That is, QEA is generally based on a binary encoding. For this reason, a key issue in the development of QEA for production scheduling problems is to design an efficient decoding mechanism to convert a binary representation into a permutation-based representation. Typically, there are mainly two decoding schemes used in QEAs in the literature for solving various scheduling problems: binary-decimal decoding and shifting decoding [30]. For the binary-decimal decoding, it first uses a binary segment for each job and then converts it into a decimal number. After that, all jobs are sequenced based on their corresponding converted decimal numbers. It is understandable that the chromosome under such a scheme is usually very long, especially when the problem size is large. As a result, the search efficiency of the algorithm may be reduced. As for shifting decoding, it uses a permutation chromosome as a parent pattern and shifts its genes with the direction of a binary chromosome so as to generate a new permutation chromosome. Such a decoding usually has a better computational efficiency than binary-decimal decoding. But it cannot make full use of the advantage of QEA due to its permutation-based representation.To overcome the above drawbacks, we propose a new decoding scheme in this study. In our scheme, a binary chromosome is directly converted into permutation chromosome (i.e., a robot move sequence) using several different decoding procedures. Our decoding scheme can efficiently exploit the solution diversity due to Q-bits chromosome compared to shifting decoding and has a shorter chromosome than binary-decimal decoding [30].The main contribution of this article is summarized as follows. Firstly, we propose a new decoding scheme as mentioned above. Secondly, we propose a more effective repairing procedure than the one in [10]. Note that, in [10], the reparation is conducted by randomly swapping any two moves for each infeasible sequence. In this article, we first identify the move segment that causes infeasibility of the entire move sequence and then repair it. Finally, to increase the population diversity, crossover and mutation operators with adaptive probabilities are also implanted into our algorithm.The remainder of this article is organized as follows. In Section2, we introduce the problem description and show an illustrative example of the problem as well as the problem formulation. The proposed algorithm with a repairing procedure is the subject of Section 3. The experimental results as well as comparisons of the proposed algorithm with the existing approach are given in Section 4. Last section gives some conclusions and perspectives.


## 2. Problem Statement and Mathematical Model

### 2.1. Problem Statement

For description consistency, we briefly give a problem description and notation, which are similar to those existing in the literature [2, 10, 12]. Given N processing machines (i.e., M1,M2,…,MN) in a production line and a single robot for part transportation, machine 0 (i.e., M0) and machine N+1 (i.e., MN+1) are the input station and the output station, respectively. After a part is unloaded from M0, it is to be successively processed through M1 to MN. The robot moving a part from Mi to Mi+1 (0≤i≤N) is called (loaded) move i, which includes three suboperations: (1) unloading a part from a machine; (2) carrying the part to the next machine; (3) loading the part into the machine. The robot movement without carrying a part between a pair of machines is called an empty move.Moreover, the part processing time at each machine is flexible and should be confined within a pair of minimum and maximum time bounds. If the actual processing time violates the time limits, defective parts would be produced. Furthermore, at any time, each machine can process only one part. When a processing operation on a machine ends, the part must be moved by the robot to the next one without delay. The production lines run in a cyclic mode. In each cycle, each machine is emptied exactly once during a cycle, which involves cyclic schedules with one-degree. This article studies the one-degree cyclic scheduling problem with a single robot, and the decision concerns how to optimize the sequence of robot movements so as to minimize the cycle time.For consistency of the problem formulation, we use the following notations and variables in this article, which are similar to [2, 12]:[ai,bi]: the minimum and maximum bounds of the processing time on Mi, respectively, 1≤i≤N.θi: the time needed to perform move i, 0≤i≤N.δi,j: the travel time for empty robot from Mi to Mj; note that δi,i=0 and δi,j=δj,i, 0≤i,j≤N+1. Note that δi,j satisfies the triangle inequalities; that is, δi,j≤δi,k+δk,j, k≠{i,j}, i≠j, 0≤i,j≤N+1.The decision variables are the following ones:T: cycle time. It is the duration of a cycle.Si: the start time of (loaded) move i within a cycle, 0≤i≤N. Without loss of generality, move 0 is supposed to be the first move of a cycle; thus S0=0.To facilitate the formulation, we define the following intermediate variables:ci: if ci=0, then Mi is empty at the beginning of a cycle; else ci=1; then Mi is occupied by a part, 0≤i≤N. Define CN={c0,c1,…,cN}. It is called the initial part distribution at the beginning of a cycle. Without loss of generality, we assume c0=1 and c1=0, since M0 is always occupied by part and move 0 is the first move of a cycle.r[i]: the i+1th move performed by the robot within a cycle, 0≤i≤N. As mentioned above, we have r[0]=0. Define RN=r[0],r[1],r[2],…,r[N], which represents the sequence of moves during a cycle. An example of RN with N=3 is R3=0,2,3,1, where r[1]=2, r[2]=3, and r[3]=1, as shown in Figure 2. Here, r[1]=2 means that the second move transfers a part from M2 to M3.Figure 2
An example of cyclic schedule with a single robot.Figure2 shows an illustrative example of the studied problem with N=3. In this example, there are three processing machines (i.e., M1, M2, and M3) with a single robot for part transportation as well as the loading station (i.e., M0) and the unloading station (i.e., M4). In Figure 2, the inclined solid arrows and the broken arrows represent the loaded moves and the empty moves, respectively. The start point and end point of an inclined solid arrow (resp., a broken arrow) represent the start time and the end time of corresponding loaded (resp., empty) move, respectively. Furthermore, the horizontal solid line represents the duration of the part processing operation. The production line is supposed to be in steady-state. As can be seen from Figure 2, at time 0, M2 is the only machine to be occupied (and implicitly M0). So the initial part distribution is C3={1,0,1,0}. For this distribution, the optimal robot move sequence is R3=0,2,3,1 (i.e., S0<S2<S3<S1). When move 1 finishes, the robot comes back to M0 and performs move 0 of the next cycle. We can also see that the robot performs the same loaded (or empty) move sequence in time interval [T,2T] as those ones in time interval [0,T]. This is called cyclic production mode. The duration of the repetitive sequence (i.e., R3) is the cycle time T.


### 2.2. Mathematical Model

As mentioned above, move 0 is supposed to start at time 0; then the start times of other moves are all greater than 0. Thus, we have(1)S0=0,Si>0,1≤i≤N.In Figure2, we notice that the start time of processing operation i is the same as the end time of loaded move i-1 (i.e., Si-1+θi-1); the end time of processing operation i is the same as the start time of loaded move i (i.e., Si). Moreover, there are in total two possible states (empty or occupied) for each machine at the beginning of a cycle. Based on the above observations, the actual processing time on Mi can be represented as Si-(Si-1+θi-1) and T+Si-(Si-1+θi-1), for ci=0 and ci=1, respectively. Considering the processing time requirements, we have(2)ai≤ciT+Si-Si-1+θi-1≤bi,1≤i≤N.Furthermore, the robot must have enough time to perform any two successive moves (i.e.,r[i] and r[i+1]); thus, the following relation holds:(3)Sri+θri+δri+1,ri+1≤Sri+1,0≤i≤N-1.When the last move (i.e.,r[N]) finishes, the robot must come back to M0 for executing move 0 of the next cycle. Hence, we have(4)SrN+θrN+δrN+1,0≤T,1≤rN≤N.Based on the above works, the mathematical model for the single robot one-degree cyclic scheduling problem with flexible processing times can be formulated as the following [2]:(5)min.Ts.t.(1)–(4).


## 2.1. Problem Statement

For description consistency, we briefly give a problem description and notation, which are similar to those existing in the literature [2, 10, 12]. Given N processing machines (i.e., M1,M2,…,MN) in a production line and a single robot for part transportation, machine 0 (i.e., M0) and machine N+1 (i.e., MN+1) are the input station and the output station, respectively. After a part is unloaded from M0, it is to be successively processed through M1 to MN. The robot moving a part from Mi to Mi+1 (0≤i≤N) is called (loaded) move i, which includes three suboperations: (1) unloading a part from a machine; (2) carrying the part to the next machine; (3) loading the part into the machine. The robot movement without carrying a part between a pair of machines is called an empty move.Moreover, the part processing time at each machine is flexible and should be confined within a pair of minimum and maximum time bounds. If the actual processing time violates the time limits, defective parts would be produced. Furthermore, at any time, each machine can process only one part. When a processing operation on a machine ends, the part must be moved by the robot to the next one without delay. The production lines run in a cyclic mode. In each cycle, each machine is emptied exactly once during a cycle, which involves cyclic schedules with one-degree. This article studies the one-degree cyclic scheduling problem with a single robot, and the decision concerns how to optimize the sequence of robot movements so as to minimize the cycle time.For consistency of the problem formulation, we use the following notations and variables in this article, which are similar to [2, 12]:[ai,bi]: the minimum and maximum bounds of the processing time on Mi, respectively, 1≤i≤N.θi: the time needed to perform move i, 0≤i≤N.δi,j: the travel time for empty robot from Mi to Mj; note that δi,i=0 and δi,j=δj,i, 0≤i,j≤N+1. Note that δi,j satisfies the triangle inequalities; that is, δi,j≤δi,k+δk,j, k≠{i,j}, i≠j, 0≤i,j≤N+1.The decision variables are the following ones:T: cycle time. It is the duration of a cycle.Si: the start time of (loaded) move i within a cycle, 0≤i≤N. Without loss of generality, move 0 is supposed to be the first move of a cycle; thus S0=0.To facilitate the formulation, we define the following intermediate variables:ci: if ci=0, then Mi is empty at the beginning of a cycle; else ci=1; then Mi is occupied by a part, 0≤i≤N. Define CN={c0,c1,…,cN}. It is called the initial part distribution at the beginning of a cycle. Without loss of generality, we assume c0=1 and c1=0, since M0 is always occupied by part and move 0 is the first move of a cycle.r[i]: the i+1th move performed by the robot within a cycle, 0≤i≤N. As mentioned above, we have r[0]=0. Define RN=r[0],r[1],r[2],…,r[N], which represents the sequence of moves during a cycle. An example of RN with N=3 is R3=0,2,3,1, where r[1]=2, r[2]=3, and r[3]=1, as shown in Figure 2. Here, r[1]=2 means that the second move transfers a part from M2 to M3.Figure 2
An example of cyclic schedule with a single robot.Figure2 shows an illustrative example of the studied problem with N=3. In this example, there are three processing machines (i.e., M1, M2, and M3) with a single robot for part transportation as well as the loading station (i.e., M0) and the unloading station (i.e., M4). In Figure 2, the inclined solid arrows and the broken arrows represent the loaded moves and the empty moves, respectively. The start point and end point of an inclined solid arrow (resp., a broken arrow) represent the start time and the end time of corresponding loaded (resp., empty) move, respectively. Furthermore, the horizontal solid line represents the duration of the part processing operation. The production line is supposed to be in steady-state. As can be seen from Figure 2, at time 0, M2 is the only machine to be occupied (and implicitly M0). So the initial part distribution is C3={1,0,1,0}. For this distribution, the optimal robot move sequence is R3=0,2,3,1 (i.e., S0<S2<S3<S1). When move 1 finishes, the robot comes back to M0 and performs move 0 of the next cycle. We can also see that the robot performs the same loaded (or empty) move sequence in time interval [T,2T] as those ones in time interval [0,T]. This is called cyclic production mode. The duration of the repetitive sequence (i.e., R3) is the cycle time T.


## 2.2. Mathematical Model

As mentioned above, move 0 is supposed to start at time 0; then the start times of other moves are all greater than 0. Thus, we have(1)S0=0,Si>0,1≤i≤N.In Figure2, we notice that the start time of processing operation i is the same as the end time of loaded move i-1 (i.e., Si-1+θi-1); the end time of processing operation i is the same as the start time of loaded move i (i.e., Si). Moreover, there are in total two possible states (empty or occupied) for each machine at the beginning of a cycle. Based on the above observations, the actual processing time on Mi can be represented as Si-(Si-1+θi-1) and T+Si-(Si-1+θi-1), for ci=0 and ci=1, respectively. Considering the processing time requirements, we have(2)ai≤ciT+Si-Si-1+θi-1≤bi,1≤i≤N.Furthermore, the robot must have enough time to perform any two successive moves (i.e.,r[i] and r[i+1]); thus, the following relation holds:(3)Sri+θri+δri+1,ri+1≤Sri+1,0≤i≤N-1.When the last move (i.e.,r[N]) finishes, the robot must come back to M0 for executing move 0 of the next cycle. Hence, we have(4)SrN+θrN+δrN+1,0≤T,1≤rN≤N.Based on the above works, the mathematical model for the single robot one-degree cyclic scheduling problem with flexible processing times can be formulated as the following [2]:(5)min.Ts.t.(1)–(4).


## 3. Hybrid Method

### 3.1. QEA Concepts

Quantum Evolutionary Algorithm (QEA) is formed based on the principles of quantum computation [25]. In a quantum computer, Q-bit is the smallest unit of information. Each Q-bit may be in “0” state, “1” state, or any superposition of the two. The following equation is usually used to define a Q-bit [27, 30]:(6)ψ=α0+β1,whereα2+β2=1.In (6), α and β are two complex numbers, which represent the probability amplitudes of states 0 and 1, respectively. As a result, α2 and β2 represent the probabilities that the Q-bit would be found in the “0” and “1” states, respectively. Suppose that a quantum individual is a string of N Q-bits [as shown below in (7)]; this individual can represent 2N states at the same time. Such a representation commonly has a better population diversity compared with other representations [30]. However, each Q-bit finally collapses to a single state by using the way presented in Section 3.3, so QEA can be seen as a probabilistic algorithm. Besides, QEA uses Q-gate to change the values of α and β so as to influence the Q-bit state. Until now, several Q-gates have been proposed in the literature [30], such as the NOT gate, controlled NOT gate, and rotation gate.QEA is like other evolutionary algorithms, such as genetic algorithm and annealing evolution algorithm, and it has a similar evolution paradigm. It begins with an initial population, in which each individual corresponds to a possible solution of a specific problem. After evaluating the population fitness, it applies Q-gate to update individuals for generating new offspring and then evaluates the new population. When the stop condition is satisfied, it ends and outputs the best solution. In what follows, we present in Section3.2 the representation; in Section 3.3, we determine the states of Q-bits in each individual; in Section 3.4, we present the decoding scheme; in Sections 3.5 and 3.6, we describe the fitness evaluation function and the repairing procedure, respectively; in Section 3.7, we introduce the rotation gate and the genetic operators to update individuals; finally, in Section 3.8, we present the pseudocode of the proposed algorithm.


### 3.2. Representation

Indeed, we notice that machine state and Q-bit state have the same characteristics. That is, they both are either 0 or 1. Since precedence relations need to be determined betweenN moves in this article, we let Q-bit i corresponding to machine i, for 1≤i≤N, and use Rule 1 and Rule 2 introduced in the following section to determine each Q-bit state. If Q-bit i is in state “0” (i.e., ci=0), move i-1 is performed before move i during a cycle; otherwise (i.e., ci=1), move i is performed before move i-1 during a cycle. Hence, an individual Ψ containing N Q-bits is used to represent N machine states and is defined as follows:(7)Ψ=α1α2⋯αNβ1β2⋯βN,where |αi|2+βi2=1, 1≤i≤N. Note that, in the initialization step, all Q-bits in Ψ are initialized as the equal probability (i.e., 1/2) of being 0 or 1. From above, we can know that each quantum individual corresponds to a complete part distribution CN, more precisely the state of each Mi (i.e., empty or occupied).In more classical and direct representations for the studied problem, each individual represents a moving sequence, so the value of genej gives the index of the machine from which the jth move starts during one cycle. In such representations, the solution space contains N! individuals. With our representation, we handle in a first step only 2N-1 individuals. This number may be further reduced for some instances with Rule 1, as explained in the following. Moreover, each Q-bits individual generally corresponds to several moving sequences, which we consider in a second step. Each time Rule 1 enables us to determine that an individual is not good, then all the associated moving sequences are unfeasible ones and there is no use evaluating them.


### 3.3. Initialization

For each specific instance, some machine states (or Q-bits states) may be directly determined by the following method. Specifically, we first suppose thatci=1; therefore move i occurs before move i-1 within a cycle. Moreover, let us suppose that move i-1 and move i are the last move and the second move of a cycle, respectively. Correspondingly, the minimum processing time of a part on machine Mi with ci=1 is δi,0+θ0+δ1,i. Indeed, move i would be the first move to be performed after move 0, and move i-1 would be the last move of the cycle. Else, the processing time on Mi would be greater than δi,0+θ0+δ1,i, which would make the following assertion even more true. Then we can compare this processing time with bi which is the maximum authorized time on Mi.(a1)
Ifbi<δi,0+θ0+δ1,i (hereafter called Rule 1) happens, then we can know that the processing time requirement on Mi is violated. Consequently, all sequences relevant with ci=1 are infeasible ones. So ci must be 0.(a2)
Else,ci may be 0 or 1.Take move 1 and move 2 in Figure2, for example, the processing time of a part on M2 is equal to δ2,0+θ0+δ1,2. Rule 1 can be used to reduce the enumerating space of CN and thus improve the search efficiency. Indeed, if Rule 1 enables us to fix 0 to the values of p variables ci, then the search space of CN can be reduced to 2N-p-1 individuals.For the state of Q-biti in a quantum individual Ψ that is not determined by Rule 1, a random number rdi is generated from the uniform distribution [0, 1). If rdi>|αi|2, then Q-bit i is in state “1” (i.e., ci=1); else, Q-bit i is in state “0” (i.e., ci=0). This method is called Rule 2. Based on the above, the states of all Q-bits in one individual can be easily determined by Rule 1 and Rule 2, that is to say, the initial part distribution CN.


### 3.4. Decoding Scheme

In what follows, we present how we derive the robot move sequence from a quantum individual. For a better diversification, three different decoding procedures described in the following are used to convert a quantum individual into possible robot move sequenceRN, providing that the states of all Q-bits (i.e., CN) in a quantum individual are already determined.


#### 3.4.1. Decoding Procedure 1

For ease of description, we first defineλi as a copy of ci and λi=ci. Let Φ be a set that records the performed moves. It should be note that λi can be seen as an indicator that indicates the state (i.e., empty or occupied) of Mi in the process. Thus, the value of λi is dynamically modified in the process. That is, when move i finishes, both the states of Mi and Mi+1 are changed; that is, Mi becomes empty and Mi+1 is occupied by a part. Thus, we set λi=0, λi+1=1 and put move i into set Φ.Procedure1 mainly depends on the probability sizes of Q-bits in Ψ to derive the move sequence. In particular, for each given CN, when move r[k] finishes, for 0≤k≤N, we first determine the possible candidates for next move r[k+1]; if move i satisfies the condition of λi=1 and λi+1=0 (note that if i=N, the output station can be seen as always empty) and i∉Φ, then we set Ω=Ω∪{i}. Notice that we use set Ω to store all possible candidates for the next mover[k+1]. Finally, we choose move j with the highest probability (i.e., αj2) in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters and use a similar way to derive the following move (i.e., r[k+2]). Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 1 illustrated the process of Decoding Procedure 1.Pseudocode 1:Pseudocode of Decoding Procedure 1.
Input:  CN=c0,c1,…,cN
Begin
fori=0toN
λi←ci;
r[0]←0; Φ←⌀; Ω←⌀;
λ1←1; λN+1←0; Φ←Φ∪{r[0]};
fori=1toN
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
forj=1toN
find  max⁡αk2,k∈Ω;
rj←k, Φ←Φ∪k,λk←0;
if(k!=N)   then  λk+1←1;
Ω←⌀;
form=1toN
if  (λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For example, a complete part distribution (corresponding to a quantum individual)CN with N=5 is C5={1,0,1,0,1,0}. When the first move (i.e., r[0]) finishes, by using the above decoding procedure, we have λ1=1, λ2=1, λ3=0, λ4=1, λ5=0, and Φ={0}, from which we can know that M1, M2, and M4 are currently occupied by a part. As the robot cannot unload a part from an empty machine and also cannot load a part into an occupied machine, we have Ω={2,4}. Finally, according to the selection rule, if |α2|2≥|α4|2, we have r[1]=2;else, r[1]=4. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


#### 3.4.2. Decoding Procedure 2

For ease of description, we keep the intermediate parametersλi, Φ, and Ω defined in Procedure 1. Furthermore, we define ti as the start time of move i in the process of deriving the whole sequence, for 0≤i≤N, and let t0=0. To derive a move sequence from each given CN, Procedure 2 mainly depends on the rule of minimal time unit increment between tr[k] and tr[k+1], for 1≤k≤N, while respecting the processing time requirement, since the objective of the problem is to minimize the cycle time T. That is to say, in each step, we have a set of several moves and choose one move with the earliest starting time as move r[k+1] from the set.In particular, on one hand, when mover[k] finishes, as similarly done in Procedure 1, we derive the value of Ω from each given CN. On the other hand, we design a different strategy to determine move r[k+1] compared with the Procedure 1. At first, we calculate each tj in set Ω; that is, tj=tr[i]+θr[i]+δr[i]+1,j. Then, for each move j in set Ω, we check whether move j-1 exists in the partial sequence r[0],r[1],…,r[i]. If it exists and tj-tj-1-θj-1<aj happens, then we update tj=tj-1+θj-1+aj so as to meet the processing requirement. Finally, we choose move j with the smallest value of tj in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters so as to derive move r[k+2]. Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 2 illustrated the process of Decoding Procedure 2.Pseudocode 2:Pseudocode of Decoding Procedure 2.
Input: CN={c0,c1,…,cN}
Begin
r[0]←0; Φ←⌀; Ω←⌀; t0←0;
for  i=0  to  N
λi←ci;
λ1←1; λN+1←0; Φ←Φ∪r0;
for  i=1  to  N
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
for  i=1  to  N
compute  tj=tr[i]+θr[i]+δr[i]+1,j, j∈Ω;
form=0toi-1,j∈Ω
if((j-r[m])==1&&tj-tj-1-θj-1<aj)
then  tj=tj-1+θj-1+aj;
find  min⁡tj, j∈Ω;
r[i]←j,Φ←Φ∪{j},λj←0;
if(j!=N) then  λj+1←1;
Ω←⌀;
for  m=1  to  N
if(λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For instance, an example ofCN with N=5 is C5={1,0,0,1,0,1}. When the first move (i.e., r[0]) finishes, by definitions, we have λ1=1, λ2=0, λ3=1, λ4=0, λ5=1, and Φ={0} as well as Ω={1,3,5}. We first calculate t1 (note that if t1-θ0<a1, then we set t1=θ0+a1) and then t3 and t5 by t0+θ0 plus δ1,1, δ1,3, δ1,5, respectively, and then choose the move with the smallest starting time among the three candidates as r[1]. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


#### 3.4.3. Decoding Procedure 3

Procedure3 mainly depends on the precedence relationship between move i-1 and move i (i.e., the value of ci) to derive the move sequence. For each given CN (i.e., a quantum individual), if ci=1, then move i is set before move i-1 in RN; else, move i is set after move i-1 in RN. Pseudocode 3 illustrated the process of Decoding Procedure 3.Pseudocode 3:Pseudocode of Decoding Procedure 3.
Input:  CN={c0,c1,…,cN}, RN=r0,r1,…,rN
Begin
for  i=1  to  N-1
p←0, q←0, k←0, l←0, m←0;
for  j=1  to  N
if (r[j]==i)  then  p←j, l←l+1;
if (r[j]==i+1)  then  q←j, m←m+1;
else if (l!=0&&m!=0) break;
if (ci+1==0) then
if (p>q) then
k←r[q], and  for  j=q  to  N-1
r[j]←r[j+1],r[N]←k;
else
if (p<q) then
k←r[q], and  for  j=q  to  2
r[j]←r[j-1],r[1]←k;
End
Output: RN=r[0],r[1],r[2],…,r[N]For instance, examples ofCN and RN with N=5 are, respectively, C5={1,0,1,1,0,1} and R5=0,2,1,4,3,5, from which we can easily derive a new sequence; that is, R5=0,3,2,1,5,4. Note that, at the initial step, we set r[i]=i, 0≤i≤N.Based on the above descriptions, we first apply the three proposed decoding procedures to each quantum individual and then select the best sequence (i.e., the best fitness) from the three generated sequences to represent this individual.


### 3.5. Fitness Evaluation

To facilitate the description,fitX is defined to represent the fitness value of each individual X. Basically the higher the fitness of an individual (i.e., fit(X)), the better the quality of individual X. Therefore, fit(X) is computed as follows: fit(X)=F/T, where F is a positive parameter and set as 1 (note that in fact F can also be set as other positive integer values, and its value does not impact the result or the procedure running) in this article. From this definition, we see that the smaller the cycle time T (T>0), the greater the fitness value. For each individual relevant with a robot move sequence, it is evaluated by using the graph-based polynomial algorithm developed by Chen et al. [2]. In particular, if the sequence is proved to be feasible, then the algorithm returns a positive value for the cycle time T and the individual fitness can be calculated; otherwise, the individual fitness is set to be 0.


### 3.6. Repairing Procedure

It should be noted that constraints (2)~(4) formulated in Section 2.2 can be regarded as two classes. One is flexible processing time constraints and the other is robot transportation capacity constraints, which are (2), (3), and (4), respectively. Generally, if a sequence RN is infeasible, the following cases happen:(C1)
the flexible processing time constraint is violated.(C2)
the robot transportation capacity constraint is violated.In what follows, we present the repairing procedure based on the above cases. For an individual with an associated robot move sequenceRN, we identify each partial sequence in a whole robot move sequence RN which is either in sequence of i-1→•→•→i (which means move i-1 is performed before move i within a cycle) or of i→•→•→i-1 (which means move i is performed before move i-1 within a cycle). That is to say, a complete robot sequence RN consists of N pieces of such a partial sequence. For ease of description, we define the following parameters:zi-1,i: the duration between the finish time of move i-1 and the start time of move i for a partial sequence i-1→•→•→i, for 1≤i≤N. Note that zi-1,i generally is equal to the sum of all loaded move (denoted by •) times and relevant empty move times. If there exists a pair of moves j-1 and j in the sequence i-1→•→j-1→•→j→i and zj-1,j<aj, then we let zi-1,i=zi-1,i+aj-zj-1,j. Note that zi-1,i may span the cycle or be within a cycle. For example, in Figure 2, the two consecutive sequences are 0→2→3→1→0→2→3→1. From it, we can see that z0,1 and z2,3 are within a cycle, but z1,2 spans the cycle. Therefore, zi-1,i can be used to check the satisfaction of flexible processing requirements no matter ci=0 or ci=1.d: the mean time of all loaded move times, d=∑i=0Nθi/(N+1).e: the mean time of all empty move times, e=2(∑i=0N+1∑j=iN+1δi,j-δ0,1)/(N2+5N). Note that the possible number of empty moves is (N2+5N+6)/2. Since the empty moves between M0 and M0, M0 and M1, and MN+1 and MN+1 do not actually happen, the number is reduced to (N2+5N)/2.For an infeasible sequenceRN, we first use parameters zi-1,i to check the sequence RN.(1)
Ifzi-1,i is verified to be greater than its upper bound bi, then we remove one or more moves from the corresponding partial sequence, so as to make the partial sequence feasible; else if zi-1,i is verified to be smaller than its lower bound ai and the time gap between ai and zi-1,i is greater than the sum of d and 2e, then we insert possible moves into the partial sequence.(2)
Then, we identify the violated robot capacity constraints by the start times of all moves (i.e.,Si, 1≤i≤N) given by the evaluation process. For ease of description, let moves i and j be the identified two moves violating the robot capacity constraints; that is, Si+θi+δi+1,j>Sj, with Si<Sj. If these two moves are two consecutive moves, we set move j before move i in sequence RN so as to make the sequence be feasible; else, we remove one or more moves between moves i and j so as to make the two moves satisfy the robot capacity constraints.


### 3.7. Update Individuals

#### 3.7.1. Rotation Gate

In this paper, the rotation gateU(Δωi) is adopted as the variation operator to update the Q-bits in (7). ω0 is set to be as the initial rotation angle. For individual X, the Q-bit i in it can be updated as the following way [27, 30]:(8)αi′βi′=cos⁡Δωi-sin⁡Δωisin⁡Δωicos⁡Δωiαiβi.We define fit_b as the fitness of the best individual found in population. The rotation angle Δωi is defined according to the respective values of the corresponding parameter ci in the individual X (labeled ci-X) and in the best one (labeled ci-best). If the condition fitX<fit_b holds, then consider the following conditions [27].Case 1.
If biti is in the first or the third quadrant, then consider the following.Case
1.1. If ci-best=1 and ci-X=0, then Δωi=-ω0; here the rotation angle Δωi is set negative so as to increase the probability that Q-bit i is in state “1.”Case
1.2. If ci-best=0 and ci-X=1, then Δωi=ω0, and the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “0.”Case
1.3. Else, Δωi=0.Case 2.
If biti is in the second or the fourth quadrant, then consider the following.Case
2.1. If ci-best=1 and ci-X=0, then Δωi=ω0; here the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “1.”Case
2.2. If ci-best=0 and ci-X=1, then Δωi=-ω0, and the rotation angle Δωi is set as negative so as to increase the probability that Q-bit i is in state “0.”Case
2.3. Else, Δωi=0;Besides, since the probability of a Q-biti in state “0” may be equal to 1 or 0, the updated Q-bit i may be trapped in state “0” or “1,” which may lead to the premature convergence of population. Thus, a small constant μ is applied to ensure that the probabilities of the two states both belonged to the range [μ,1-μ]. As a result, the following equation must be considered [34]:(9)αi′βi′=μ,1-μT,ifαi′<μ,1-μ,μT,ifαi′>μ,α′β′T,else.By applying the decoding procedures given in Section3.4 to each updated quantum individual, robot move sequences can be generated from it.


#### 3.7.2. Genetic Operators

In this subsection, selection, crossover, and mutation operators are applied to further evolve the population [35]. To facilitate the description, the following notations are given:Pc, Pm: crossover and mutation probabilities, respectively.fit_a: the average fitness of the entire population.fit_0: the maximum fitness of a specific instance, which is computed as follows:fit_0=1/TL. TL is the lower bound on cycle time T for the instance. It can be obtained by the following way, which is taken from Chen et al. [2]:(10)TL≥max⁡ai+θi+θi-1+δi+1,i-1,1≤i≤N.According to Srinivas and Patnaik [36], Pc and Pm are defined, respectively, in a similar way:(11)Pc=0.7×fit_0-fit_bfit_0-fit_a,Pm=0.5×fit_0-fitXfit_0-fit_a.Adaptively adjustingPc and Pm (i.e., (11)) can prevent divergence and escape from the local optimal, since (11) can dynamically reduce Pc and Pm for individuals with high fitness or increase Pc and Pm for individuals with low fitness.In this research, two-point crossover operator is applied to generate the offspring. First, two individuals are chosen by the binary tournament method as parents 1 and 2; then, for parent 1, two different positionsp and q are randomly chosen (p<q), p,q∈[1,N]. For i∈[1,p) and (q,N], the values of r[i] for the new offspring 1 are inherited from parent 1. For i∈[p,q], the new r[i] is sequentially chosen from parent 2, on condition that its value was not already chosen from parent 1. The same operations are done, starting with parent 2 and then parent 1, to generate offspring 2. This operation is depicted in Figure 3(a), in which ∣ is the chosen position.Figure 3
Crossover and mutation operators.
(a)
Two-point crossover operation
(b)
Mutation operationBesides, a mutation operator is adopted to prevent a solution falling into a local optimum, which is designed as follows. For a chosen individualRN=r[0],r[1],r[2],…,r[N], first, we randomly choose a positionp, p∈1,N and then reorder the partial move sequence r[p+1],r[p+2],…,r[N]. This operation is depicted in Figure 3(b), in which ∣ is the chosen position.


### 3.8. Pseudocode and Complexity Analysis of the Proposed HQEA

#### 3.8.1. Pseudocode of the Proposed HQEA

Based on the above works presented in Sections3.1–3.7, the pseudocode of HQEA for solving the considered problem can be written as follows.Input:
Popsize (population size).MaxIter (maximum number of generations).ω0 (initial rotation angle).Rt (maximum repairing times).N (problem size)Output: the best cycle time T and its related robot move sequence RN.Begin
(a)
t←0;(b)
Initialize Q-bits chromosomes;(c)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(d)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(e)
Evaluate all obtained sequences and repair the infeasible ones;(f)
Assign fitness to individuals;(g)
Sort all individuals based on their fitness from best to worse;(h)
Store the best solution at generationt;While (t<MaxIter)doBegin
(i)
t←t+1;(j)
Update each Q-bits chromosome by using quantum rotation gate;(k)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(l)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(m)
Evaluate all obtained sequences and repair the infeasible ones;(n)
Assign fitness to individuals;(o)
Sort all individuals based on their fitness from best to worse;(p)
Store the best solution at generationt;(q)
Apply the proposed adaptive genetic operators to the obtained sequences;(r)
Evaluate the obtained sequences and repair the infeasible ones;(s)
Assign fitness to individuals;(t)
Sort individuals and update the best solution;EndIf the stopping condition is reached, then the algorithm output the best solution, including cycle timeT and its related robot move sequence RN.End


#### 3.8.2. Complexity Analysis of the Proposed HQEA

The computational complexity of the proposed HQEA is analyzed as follows. Step (a) and Step (h), respectively, takeO(1) and O(N) times. In step (b), initializing the entire population takes O(2PopsizeN) time. In step (c), determining the Q-bit state for all quantum individuals takes O(Popsize(N-1)) time. In step (d), decoding all the quantum individuals takes O(3PopsizeN(N-1)) time. According to [2], evaluating all the individuals in step (e) takes O(27N6Popsize) time. Besides, the repairing procedure in step (e) takes O(3N2(N-1)RtPopsize)time. The fitness assignment operation in step (f) takes O(3Popsize)time. Sequencing all the individuals in step (g) takes O(Popsizelog⁡Popsize) time. Consequently, the computational complexity for step (a)–step (h) is approximately O(27N6Popsize).In thewhile loop, updating all the quantum individuals in step (j) takes O(PopsizeN) time. Step (q) takes O(2PopsizeN) time. Similar to step (a)–step (h), the computational complexity for other steps in thewhile loop is approximately O(36N6Popsize). Since thewhile loop executesMaxIter times, the computational complexity for thewhile loop is about O(36N6PopsizeMaxIter). Then, the computational complexity is approximately O(36N6PopsizeMaxIter)+O(27N6Popsize), that is, O(9N6Popsize(4MaxIter+3)).


## 3.1. QEA Concepts

Quantum Evolutionary Algorithm (QEA) is formed based on the principles of quantum computation [25]. In a quantum computer, Q-bit is the smallest unit of information. Each Q-bit may be in “0” state, “1” state, or any superposition of the two. The following equation is usually used to define a Q-bit [27, 30]:(6)ψ=α0+β1,whereα2+β2=1.In (6), α and β are two complex numbers, which represent the probability amplitudes of states 0 and 1, respectively. As a result, α2 and β2 represent the probabilities that the Q-bit would be found in the “0” and “1” states, respectively. Suppose that a quantum individual is a string of N Q-bits [as shown below in (7)]; this individual can represent 2N states at the same time. Such a representation commonly has a better population diversity compared with other representations [30]. However, each Q-bit finally collapses to a single state by using the way presented in Section 3.3, so QEA can be seen as a probabilistic algorithm. Besides, QEA uses Q-gate to change the values of α and β so as to influence the Q-bit state. Until now, several Q-gates have been proposed in the literature [30], such as the NOT gate, controlled NOT gate, and rotation gate.QEA is like other evolutionary algorithms, such as genetic algorithm and annealing evolution algorithm, and it has a similar evolution paradigm. It begins with an initial population, in which each individual corresponds to a possible solution of a specific problem. After evaluating the population fitness, it applies Q-gate to update individuals for generating new offspring and then evaluates the new population. When the stop condition is satisfied, it ends and outputs the best solution. In what follows, we present in Section3.2 the representation; in Section 3.3, we determine the states of Q-bits in each individual; in Section 3.4, we present the decoding scheme; in Sections 3.5 and 3.6, we describe the fitness evaluation function and the repairing procedure, respectively; in Section 3.7, we introduce the rotation gate and the genetic operators to update individuals; finally, in Section 3.8, we present the pseudocode of the proposed algorithm.


## 3.2. Representation

Indeed, we notice that machine state and Q-bit state have the same characteristics. That is, they both are either 0 or 1. Since precedence relations need to be determined betweenN moves in this article, we let Q-bit i corresponding to machine i, for 1≤i≤N, and use Rule 1 and Rule 2 introduced in the following section to determine each Q-bit state. If Q-bit i is in state “0” (i.e., ci=0), move i-1 is performed before move i during a cycle; otherwise (i.e., ci=1), move i is performed before move i-1 during a cycle. Hence, an individual Ψ containing N Q-bits is used to represent N machine states and is defined as follows:(7)Ψ=α1α2⋯αNβ1β2⋯βN,where |αi|2+βi2=1, 1≤i≤N. Note that, in the initialization step, all Q-bits in Ψ are initialized as the equal probability (i.e., 1/2) of being 0 or 1. From above, we can know that each quantum individual corresponds to a complete part distribution CN, more precisely the state of each Mi (i.e., empty or occupied).In more classical and direct representations for the studied problem, each individual represents a moving sequence, so the value of genej gives the index of the machine from which the jth move starts during one cycle. In such representations, the solution space contains N! individuals. With our representation, we handle in a first step only 2N-1 individuals. This number may be further reduced for some instances with Rule 1, as explained in the following. Moreover, each Q-bits individual generally corresponds to several moving sequences, which we consider in a second step. Each time Rule 1 enables us to determine that an individual is not good, then all the associated moving sequences are unfeasible ones and there is no use evaluating them.


## 3.3. Initialization

For each specific instance, some machine states (or Q-bits states) may be directly determined by the following method. Specifically, we first suppose thatci=1; therefore move i occurs before move i-1 within a cycle. Moreover, let us suppose that move i-1 and move i are the last move and the second move of a cycle, respectively. Correspondingly, the minimum processing time of a part on machine Mi with ci=1 is δi,0+θ0+δ1,i. Indeed, move i would be the first move to be performed after move 0, and move i-1 would be the last move of the cycle. Else, the processing time on Mi would be greater than δi,0+θ0+δ1,i, which would make the following assertion even more true. Then we can compare this processing time with bi which is the maximum authorized time on Mi.(a1)
Ifbi<δi,0+θ0+δ1,i (hereafter called Rule 1) happens, then we can know that the processing time requirement on Mi is violated. Consequently, all sequences relevant with ci=1 are infeasible ones. So ci must be 0.(a2)
Else,ci may be 0 or 1.Take move 1 and move 2 in Figure2, for example, the processing time of a part on M2 is equal to δ2,0+θ0+δ1,2. Rule 1 can be used to reduce the enumerating space of CN and thus improve the search efficiency. Indeed, if Rule 1 enables us to fix 0 to the values of p variables ci, then the search space of CN can be reduced to 2N-p-1 individuals.For the state of Q-biti in a quantum individual Ψ that is not determined by Rule 1, a random number rdi is generated from the uniform distribution [0, 1). If rdi>|αi|2, then Q-bit i is in state “1” (i.e., ci=1); else, Q-bit i is in state “0” (i.e., ci=0). This method is called Rule 2. Based on the above, the states of all Q-bits in one individual can be easily determined by Rule 1 and Rule 2, that is to say, the initial part distribution CN.


## 3.4. Decoding Scheme

In what follows, we present how we derive the robot move sequence from a quantum individual. For a better diversification, three different decoding procedures described in the following are used to convert a quantum individual into possible robot move sequenceRN, providing that the states of all Q-bits (i.e., CN) in a quantum individual are already determined.


### 3.4.1. Decoding Procedure 1

For ease of description, we first defineλi as a copy of ci and λi=ci. Let Φ be a set that records the performed moves. It should be note that λi can be seen as an indicator that indicates the state (i.e., empty or occupied) of Mi in the process. Thus, the value of λi is dynamically modified in the process. That is, when move i finishes, both the states of Mi and Mi+1 are changed; that is, Mi becomes empty and Mi+1 is occupied by a part. Thus, we set λi=0, λi+1=1 and put move i into set Φ.Procedure1 mainly depends on the probability sizes of Q-bits in Ψ to derive the move sequence. In particular, for each given CN, when move r[k] finishes, for 0≤k≤N, we first determine the possible candidates for next move r[k+1]; if move i satisfies the condition of λi=1 and λi+1=0 (note that if i=N, the output station can be seen as always empty) and i∉Φ, then we set Ω=Ω∪{i}. Notice that we use set Ω to store all possible candidates for the next mover[k+1]. Finally, we choose move j with the highest probability (i.e., αj2) in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters and use a similar way to derive the following move (i.e., r[k+2]). Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 1 illustrated the process of Decoding Procedure 1.Pseudocode 1:Pseudocode of Decoding Procedure 1.
Input:  CN=c0,c1,…,cN
Begin
fori=0toN
λi←ci;
r[0]←0; Φ←⌀; Ω←⌀;
λ1←1; λN+1←0; Φ←Φ∪{r[0]};
fori=1toN
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
forj=1toN
find  max⁡αk2,k∈Ω;
rj←k, Φ←Φ∪k,λk←0;
if(k!=N)   then  λk+1←1;
Ω←⌀;
form=1toN
if  (λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For example, a complete part distribution (corresponding to a quantum individual)CN with N=5 is C5={1,0,1,0,1,0}. When the first move (i.e., r[0]) finishes, by using the above decoding procedure, we have λ1=1, λ2=1, λ3=0, λ4=1, λ5=0, and Φ={0}, from which we can know that M1, M2, and M4 are currently occupied by a part. As the robot cannot unload a part from an empty machine and also cannot load a part into an occupied machine, we have Ω={2,4}. Finally, according to the selection rule, if |α2|2≥|α4|2, we have r[1]=2;else, r[1]=4. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


### 3.4.2. Decoding Procedure 2

For ease of description, we keep the intermediate parametersλi, Φ, and Ω defined in Procedure 1. Furthermore, we define ti as the start time of move i in the process of deriving the whole sequence, for 0≤i≤N, and let t0=0. To derive a move sequence from each given CN, Procedure 2 mainly depends on the rule of minimal time unit increment between tr[k] and tr[k+1], for 1≤k≤N, while respecting the processing time requirement, since the objective of the problem is to minimize the cycle time T. That is to say, in each step, we have a set of several moves and choose one move with the earliest starting time as move r[k+1] from the set.In particular, on one hand, when mover[k] finishes, as similarly done in Procedure 1, we derive the value of Ω from each given CN. On the other hand, we design a different strategy to determine move r[k+1] compared with the Procedure 1. At first, we calculate each tj in set Ω; that is, tj=tr[i]+θr[i]+δr[i]+1,j. Then, for each move j in set Ω, we check whether move j-1 exists in the partial sequence r[0],r[1],…,r[i]. If it exists and tj-tj-1-θj-1<aj happens, then we update tj=tj-1+θj-1+aj so as to meet the processing requirement. Finally, we choose move j with the smallest value of tj in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters so as to derive move r[k+2]. Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 2 illustrated the process of Decoding Procedure 2.Pseudocode 2:Pseudocode of Decoding Procedure 2.
Input: CN={c0,c1,…,cN}
Begin
r[0]←0; Φ←⌀; Ω←⌀; t0←0;
for  i=0  to  N
λi←ci;
λ1←1; λN+1←0; Φ←Φ∪r0;
for  i=1  to  N
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
for  i=1  to  N
compute  tj=tr[i]+θr[i]+δr[i]+1,j, j∈Ω;
form=0toi-1,j∈Ω
if((j-r[m])==1&&tj-tj-1-θj-1<aj)
then  tj=tj-1+θj-1+aj;
find  min⁡tj, j∈Ω;
r[i]←j,Φ←Φ∪{j},λj←0;
if(j!=N) then  λj+1←1;
Ω←⌀;
for  m=1  to  N
if(λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For instance, an example ofCN with N=5 is C5={1,0,0,1,0,1}. When the first move (i.e., r[0]) finishes, by definitions, we have λ1=1, λ2=0, λ3=1, λ4=0, λ5=1, and Φ={0} as well as Ω={1,3,5}. We first calculate t1 (note that if t1-θ0<a1, then we set t1=θ0+a1) and then t3 and t5 by t0+θ0 plus δ1,1, δ1,3, δ1,5, respectively, and then choose the move with the smallest starting time among the three candidates as r[1]. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


### 3.4.3. Decoding Procedure 3

Procedure3 mainly depends on the precedence relationship between move i-1 and move i (i.e., the value of ci) to derive the move sequence. For each given CN (i.e., a quantum individual), if ci=1, then move i is set before move i-1 in RN; else, move i is set after move i-1 in RN. Pseudocode 3 illustrated the process of Decoding Procedure 3.Pseudocode 3:Pseudocode of Decoding Procedure 3.
Input:  CN={c0,c1,…,cN}, RN=r0,r1,…,rN
Begin
for  i=1  to  N-1
p←0, q←0, k←0, l←0, m←0;
for  j=1  to  N
if (r[j]==i)  then  p←j, l←l+1;
if (r[j]==i+1)  then  q←j, m←m+1;
else if (l!=0&&m!=0) break;
if (ci+1==0) then
if (p>q) then
k←r[q], and  for  j=q  to  N-1
r[j]←r[j+1],r[N]←k;
else
if (p<q) then
k←r[q], and  for  j=q  to  2
r[j]←r[j-1],r[1]←k;
End
Output: RN=r[0],r[1],r[2],…,r[N]For instance, examples ofCN and RN with N=5 are, respectively, C5={1,0,1,1,0,1} and R5=0,2,1,4,3,5, from which we can easily derive a new sequence; that is, R5=0,3,2,1,5,4. Note that, at the initial step, we set r[i]=i, 0≤i≤N.Based on the above descriptions, we first apply the three proposed decoding procedures to each quantum individual and then select the best sequence (i.e., the best fitness) from the three generated sequences to represent this individual.


## 3.4.1. Decoding Procedure 1

For ease of description, we first defineλi as a copy of ci and λi=ci. Let Φ be a set that records the performed moves. It should be note that λi can be seen as an indicator that indicates the state (i.e., empty or occupied) of Mi in the process. Thus, the value of λi is dynamically modified in the process. That is, when move i finishes, both the states of Mi and Mi+1 are changed; that is, Mi becomes empty and Mi+1 is occupied by a part. Thus, we set λi=0, λi+1=1 and put move i into set Φ.Procedure1 mainly depends on the probability sizes of Q-bits in Ψ to derive the move sequence. In particular, for each given CN, when move r[k] finishes, for 0≤k≤N, we first determine the possible candidates for next move r[k+1]; if move i satisfies the condition of λi=1 and λi+1=0 (note that if i=N, the output station can be seen as always empty) and i∉Φ, then we set Ω=Ω∪{i}. Notice that we use set Ω to store all possible candidates for the next mover[k+1]. Finally, we choose move j with the highest probability (i.e., αj2) in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters and use a similar way to derive the following move (i.e., r[k+2]). Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 1 illustrated the process of Decoding Procedure 1.Pseudocode 1:Pseudocode of Decoding Procedure 1.
Input:  CN=c0,c1,…,cN
Begin
fori=0toN
λi←ci;
r[0]←0; Φ←⌀; Ω←⌀;
λ1←1; λN+1←0; Φ←Φ∪{r[0]};
fori=1toN
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
forj=1toN
find  max⁡αk2,k∈Ω;
rj←k, Φ←Φ∪k,λk←0;
if(k!=N)   then  λk+1←1;
Ω←⌀;
form=1toN
if  (λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For example, a complete part distribution (corresponding to a quantum individual)CN with N=5 is C5={1,0,1,0,1,0}. When the first move (i.e., r[0]) finishes, by using the above decoding procedure, we have λ1=1, λ2=1, λ3=0, λ4=1, λ5=0, and Φ={0}, from which we can know that M1, M2, and M4 are currently occupied by a part. As the robot cannot unload a part from an empty machine and also cannot load a part into an occupied machine, we have Ω={2,4}. Finally, according to the selection rule, if |α2|2≥|α4|2, we have r[1]=2;else, r[1]=4. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


## 3.4.2. Decoding Procedure 2

For ease of description, we keep the intermediate parametersλi, Φ, and Ω defined in Procedure 1. Furthermore, we define ti as the start time of move i in the process of deriving the whole sequence, for 0≤i≤N, and let t0=0. To derive a move sequence from each given CN, Procedure 2 mainly depends on the rule of minimal time unit increment between tr[k] and tr[k+1], for 1≤k≤N, while respecting the processing time requirement, since the objective of the problem is to minimize the cycle time T. That is to say, in each step, we have a set of several moves and choose one move with the earliest starting time as move r[k+1] from the set.In particular, on one hand, when mover[k] finishes, as similarly done in Procedure 1, we derive the value of Ω from each given CN. On the other hand, we design a different strategy to determine move r[k+1] compared with the Procedure 1. At first, we calculate each tj in set Ω; that is, tj=tr[i]+θr[i]+δr[i]+1,j. Then, for each move j in set Ω, we check whether move j-1 exists in the partial sequence r[0],r[1],…,r[i]. If it exists and tj-tj-1-θj-1<aj happens, then we update tj=tj-1+θj-1+aj so as to meet the processing requirement. Finally, we choose move j with the smallest value of tj in set Ω as move r[k+1], and let λj=0, λj+1=1 (for j≠N), Φ=Φ∪{j}, and Ω=⌀. In the next step, we update all related parameters so as to derive move r[k+2]. Once the whole sequence (i.e., RN) is determined, this procedure stops. Pseudocode 2 illustrated the process of Decoding Procedure 2.Pseudocode 2:Pseudocode of Decoding Procedure 2.
Input: CN={c0,c1,…,cN}
Begin
r[0]←0; Φ←⌀; Ω←⌀; t0←0;
for  i=0  to  N
λi←ci;
λ1←1; λN+1←0; Φ←Φ∪r0;
for  i=1  to  N
if (λi==1&&λi+1==0) then
Ω←Ω∪{i};
for  i=1  to  N
compute  tj=tr[i]+θr[i]+δr[i]+1,j, j∈Ω;
form=0toi-1,j∈Ω
if((j-r[m])==1&&tj-tj-1-θj-1<aj)
then  tj=tj-1+θj-1+aj;
find  min⁡tj, j∈Ω;
r[i]←j,Φ←Φ∪{j},λj←0;
if(j!=N) then  λj+1←1;
Ω←⌀;
for  m=1  to  N
if(λm==1&&λm+1==0&&m∉Φ)
then  Ω←Ω∪{m};
End
Output:  RN=r[0],r[1],r[2],…,r[N]For instance, an example ofCN with N=5 is C5={1,0,0,1,0,1}. When the first move (i.e., r[0]) finishes, by definitions, we have λ1=1, λ2=0, λ3=1, λ4=0, λ5=1, and Φ={0} as well as Ω={1,3,5}. We first calculate t1 (note that if t1-θ0<a1, then we set t1=θ0+a1) and then t3 and t5 by t0+θ0 plus δ1,1, δ1,3, δ1,5, respectively, and then choose the move with the smallest starting time among the three candidates as r[1]. The similar ways are used to update λi, Φ, Ω and then determine r[k], 2≤k≤5.


## 3.4.3. Decoding Procedure 3

Procedure3 mainly depends on the precedence relationship between move i-1 and move i (i.e., the value of ci) to derive the move sequence. For each given CN (i.e., a quantum individual), if ci=1, then move i is set before move i-1 in RN; else, move i is set after move i-1 in RN. Pseudocode 3 illustrated the process of Decoding Procedure 3.Pseudocode 3:Pseudocode of Decoding Procedure 3.
Input:  CN={c0,c1,…,cN}, RN=r0,r1,…,rN
Begin
for  i=1  to  N-1
p←0, q←0, k←0, l←0, m←0;
for  j=1  to  N
if (r[j]==i)  then  p←j, l←l+1;
if (r[j]==i+1)  then  q←j, m←m+1;
else if (l!=0&&m!=0) break;
if (ci+1==0) then
if (p>q) then
k←r[q], and  for  j=q  to  N-1
r[j]←r[j+1],r[N]←k;
else
if (p<q) then
k←r[q], and  for  j=q  to  2
r[j]←r[j-1],r[1]←k;
End
Output: RN=r[0],r[1],r[2],…,r[N]For instance, examples ofCN and RN with N=5 are, respectively, C5={1,0,1,1,0,1} and R5=0,2,1,4,3,5, from which we can easily derive a new sequence; that is, R5=0,3,2,1,5,4. Note that, at the initial step, we set r[i]=i, 0≤i≤N.Based on the above descriptions, we first apply the three proposed decoding procedures to each quantum individual and then select the best sequence (i.e., the best fitness) from the three generated sequences to represent this individual.


## 3.5. Fitness Evaluation

To facilitate the description,fitX is defined to represent the fitness value of each individual X. Basically the higher the fitness of an individual (i.e., fit(X)), the better the quality of individual X. Therefore, fit(X) is computed as follows: fit(X)=F/T, where F is a positive parameter and set as 1 (note that in fact F can also be set as other positive integer values, and its value does not impact the result or the procedure running) in this article. From this definition, we see that the smaller the cycle time T (T>0), the greater the fitness value. For each individual relevant with a robot move sequence, it is evaluated by using the graph-based polynomial algorithm developed by Chen et al. [2]. In particular, if the sequence is proved to be feasible, then the algorithm returns a positive value for the cycle time T and the individual fitness can be calculated; otherwise, the individual fitness is set to be 0.


## 3.6. Repairing Procedure

It should be noted that constraints (2)~(4) formulated in Section 2.2 can be regarded as two classes. One is flexible processing time constraints and the other is robot transportation capacity constraints, which are (2), (3), and (4), respectively. Generally, if a sequence RN is infeasible, the following cases happen:(C1)
the flexible processing time constraint is violated.(C2)
the robot transportation capacity constraint is violated.In what follows, we present the repairing procedure based on the above cases. For an individual with an associated robot move sequenceRN, we identify each partial sequence in a whole robot move sequence RN which is either in sequence of i-1→•→•→i (which means move i-1 is performed before move i within a cycle) or of i→•→•→i-1 (which means move i is performed before move i-1 within a cycle). That is to say, a complete robot sequence RN consists of N pieces of such a partial sequence. For ease of description, we define the following parameters:zi-1,i: the duration between the finish time of move i-1 and the start time of move i for a partial sequence i-1→•→•→i, for 1≤i≤N. Note that zi-1,i generally is equal to the sum of all loaded move (denoted by •) times and relevant empty move times. If there exists a pair of moves j-1 and j in the sequence i-1→•→j-1→•→j→i and zj-1,j<aj, then we let zi-1,i=zi-1,i+aj-zj-1,j. Note that zi-1,i may span the cycle or be within a cycle. For example, in Figure 2, the two consecutive sequences are 0→2→3→1→0→2→3→1. From it, we can see that z0,1 and z2,3 are within a cycle, but z1,2 spans the cycle. Therefore, zi-1,i can be used to check the satisfaction of flexible processing requirements no matter ci=0 or ci=1.d: the mean time of all loaded move times, d=∑i=0Nθi/(N+1).e: the mean time of all empty move times, e=2(∑i=0N+1∑j=iN+1δi,j-δ0,1)/(N2+5N). Note that the possible number of empty moves is (N2+5N+6)/2. Since the empty moves between M0 and M0, M0 and M1, and MN+1 and MN+1 do not actually happen, the number is reduced to (N2+5N)/2.For an infeasible sequenceRN, we first use parameters zi-1,i to check the sequence RN.(1)
Ifzi-1,i is verified to be greater than its upper bound bi, then we remove one or more moves from the corresponding partial sequence, so as to make the partial sequence feasible; else if zi-1,i is verified to be smaller than its lower bound ai and the time gap between ai and zi-1,i is greater than the sum of d and 2e, then we insert possible moves into the partial sequence.(2)
Then, we identify the violated robot capacity constraints by the start times of all moves (i.e.,Si, 1≤i≤N) given by the evaluation process. For ease of description, let moves i and j be the identified two moves violating the robot capacity constraints; that is, Si+θi+δi+1,j>Sj, with Si<Sj. If these two moves are two consecutive moves, we set move j before move i in sequence RN so as to make the sequence be feasible; else, we remove one or more moves between moves i and j so as to make the two moves satisfy the robot capacity constraints.


## 3.7. Update Individuals

### 3.7.1. Rotation Gate

In this paper, the rotation gateU(Δωi) is adopted as the variation operator to update the Q-bits in (7). ω0 is set to be as the initial rotation angle. For individual X, the Q-bit i in it can be updated as the following way [27, 30]:(8)αi′βi′=cos⁡Δωi-sin⁡Δωisin⁡Δωicos⁡Δωiαiβi.We define fit_b as the fitness of the best individual found in population. The rotation angle Δωi is defined according to the respective values of the corresponding parameter ci in the individual X (labeled ci-X) and in the best one (labeled ci-best). If the condition fitX<fit_b holds, then consider the following conditions [27].Case 1.
If biti is in the first or the third quadrant, then consider the following.Case
1.1. If ci-best=1 and ci-X=0, then Δωi=-ω0; here the rotation angle Δωi is set negative so as to increase the probability that Q-bit i is in state “1.”Case
1.2. If ci-best=0 and ci-X=1, then Δωi=ω0, and the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “0.”Case
1.3. Else, Δωi=0.Case 2.
If biti is in the second or the fourth quadrant, then consider the following.Case
2.1. If ci-best=1 and ci-X=0, then Δωi=ω0; here the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “1.”Case
2.2. If ci-best=0 and ci-X=1, then Δωi=-ω0, and the rotation angle Δωi is set as negative so as to increase the probability that Q-bit i is in state “0.”Case
2.3. Else, Δωi=0;Besides, since the probability of a Q-biti in state “0” may be equal to 1 or 0, the updated Q-bit i may be trapped in state “0” or “1,” which may lead to the premature convergence of population. Thus, a small constant μ is applied to ensure that the probabilities of the two states both belonged to the range [μ,1-μ]. As a result, the following equation must be considered [34]:(9)αi′βi′=μ,1-μT,ifαi′<μ,1-μ,μT,ifαi′>μ,α′β′T,else.By applying the decoding procedures given in Section3.4 to each updated quantum individual, robot move sequences can be generated from it.


### 3.7.2. Genetic Operators

In this subsection, selection, crossover, and mutation operators are applied to further evolve the population [35]. To facilitate the description, the following notations are given:Pc, Pm: crossover and mutation probabilities, respectively.fit_a: the average fitness of the entire population.fit_0: the maximum fitness of a specific instance, which is computed as follows:fit_0=1/TL. TL is the lower bound on cycle time T for the instance. It can be obtained by the following way, which is taken from Chen et al. [2]:(10)TL≥max⁡ai+θi+θi-1+δi+1,i-1,1≤i≤N.According to Srinivas and Patnaik [36], Pc and Pm are defined, respectively, in a similar way:(11)Pc=0.7×fit_0-fit_bfit_0-fit_a,Pm=0.5×fit_0-fitXfit_0-fit_a.Adaptively adjustingPc and Pm (i.e., (11)) can prevent divergence and escape from the local optimal, since (11) can dynamically reduce Pc and Pm for individuals with high fitness or increase Pc and Pm for individuals with low fitness.In this research, two-point crossover operator is applied to generate the offspring. First, two individuals are chosen by the binary tournament method as parents 1 and 2; then, for parent 1, two different positionsp and q are randomly chosen (p<q), p,q∈[1,N]. For i∈[1,p) and (q,N], the values of r[i] for the new offspring 1 are inherited from parent 1. For i∈[p,q], the new r[i] is sequentially chosen from parent 2, on condition that its value was not already chosen from parent 1. The same operations are done, starting with parent 2 and then parent 1, to generate offspring 2. This operation is depicted in Figure 3(a), in which ∣ is the chosen position.Figure 3
Crossover and mutation operators.
(a)
Two-point crossover operation
(b)
Mutation operationBesides, a mutation operator is adopted to prevent a solution falling into a local optimum, which is designed as follows. For a chosen individualRN=r[0],r[1],r[2],…,r[N], first, we randomly choose a positionp, p∈1,N and then reorder the partial move sequence r[p+1],r[p+2],…,r[N]. This operation is depicted in Figure 3(b), in which ∣ is the chosen position.


## 3.7.1. Rotation Gate

In this paper, the rotation gateU(Δωi) is adopted as the variation operator to update the Q-bits in (7). ω0 is set to be as the initial rotation angle. For individual X, the Q-bit i in it can be updated as the following way [27, 30]:(8)αi′βi′=cos⁡Δωi-sin⁡Δωisin⁡Δωicos⁡Δωiαiβi.We define fit_b as the fitness of the best individual found in population. The rotation angle Δωi is defined according to the respective values of the corresponding parameter ci in the individual X (labeled ci-X) and in the best one (labeled ci-best). If the condition fitX<fit_b holds, then consider the following conditions [27].Case 1.
If biti is in the first or the third quadrant, then consider the following.Case
1.1. If ci-best=1 and ci-X=0, then Δωi=-ω0; here the rotation angle Δωi is set negative so as to increase the probability that Q-bit i is in state “1.”Case
1.2. If ci-best=0 and ci-X=1, then Δωi=ω0, and the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “0.”Case
1.3. Else, Δωi=0.Case 2.
If biti is in the second or the fourth quadrant, then consider the following.Case
2.1. If ci-best=1 and ci-X=0, then Δωi=ω0; here the rotation angle Δωi is set as positive so as to increase the probability that Q-bit i is in state “1.”Case
2.2. If ci-best=0 and ci-X=1, then Δωi=-ω0, and the rotation angle Δωi is set as negative so as to increase the probability that Q-bit i is in state “0.”Case
2.3. Else, Δωi=0;Besides, since the probability of a Q-biti in state “0” may be equal to 1 or 0, the updated Q-bit i may be trapped in state “0” or “1,” which may lead to the premature convergence of population. Thus, a small constant μ is applied to ensure that the probabilities of the two states both belonged to the range [μ,1-μ]. As a result, the following equation must be considered [34]:(9)αi′βi′=μ,1-μT,ifαi′<μ,1-μ,μT,ifαi′>μ,α′β′T,else.By applying the decoding procedures given in Section3.4 to each updated quantum individual, robot move sequences can be generated from it.


## 3.7.2. Genetic Operators

In this subsection, selection, crossover, and mutation operators are applied to further evolve the population [35]. To facilitate the description, the following notations are given:Pc, Pm: crossover and mutation probabilities, respectively.fit_a: the average fitness of the entire population.fit_0: the maximum fitness of a specific instance, which is computed as follows:fit_0=1/TL. TL is the lower bound on cycle time T for the instance. It can be obtained by the following way, which is taken from Chen et al. [2]:(10)TL≥max⁡ai+θi+θi-1+δi+1,i-1,1≤i≤N.According to Srinivas and Patnaik [36], Pc and Pm are defined, respectively, in a similar way:(11)Pc=0.7×fit_0-fit_bfit_0-fit_a,Pm=0.5×fit_0-fitXfit_0-fit_a.Adaptively adjustingPc and Pm (i.e., (11)) can prevent divergence and escape from the local optimal, since (11) can dynamically reduce Pc and Pm for individuals with high fitness or increase Pc and Pm for individuals with low fitness.In this research, two-point crossover operator is applied to generate the offspring. First, two individuals are chosen by the binary tournament method as parents 1 and 2; then, for parent 1, two different positionsp and q are randomly chosen (p<q), p,q∈[1,N]. For i∈[1,p) and (q,N], the values of r[i] for the new offspring 1 are inherited from parent 1. For i∈[p,q], the new r[i] is sequentially chosen from parent 2, on condition that its value was not already chosen from parent 1. The same operations are done, starting with parent 2 and then parent 1, to generate offspring 2. This operation is depicted in Figure 3(a), in which ∣ is the chosen position.Figure 3
Crossover and mutation operators.
(a)
Two-point crossover operation
(b)
Mutation operationBesides, a mutation operator is adopted to prevent a solution falling into a local optimum, which is designed as follows. For a chosen individualRN=r[0],r[1],r[2],…,r[N], first, we randomly choose a positionp, p∈1,N and then reorder the partial move sequence r[p+1],r[p+2],…,r[N]. This operation is depicted in Figure 3(b), in which ∣ is the chosen position.


## 3.8. Pseudocode and Complexity Analysis of the Proposed HQEA

### 3.8.1. Pseudocode of the Proposed HQEA

Based on the above works presented in Sections3.1–3.7, the pseudocode of HQEA for solving the considered problem can be written as follows.Input:
Popsize (population size).MaxIter (maximum number of generations).ω0 (initial rotation angle).Rt (maximum repairing times).N (problem size)Output: the best cycle time T and its related robot move sequence RN.Begin
(a)
t←0;(b)
Initialize Q-bits chromosomes;(c)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(d)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(e)
Evaluate all obtained sequences and repair the infeasible ones;(f)
Assign fitness to individuals;(g)
Sort all individuals based on their fitness from best to worse;(h)
Store the best solution at generationt;While (t<MaxIter)doBegin
(i)
t←t+1;(j)
Update each Q-bits chromosome by using quantum rotation gate;(k)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(l)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(m)
Evaluate all obtained sequences and repair the infeasible ones;(n)
Assign fitness to individuals;(o)
Sort all individuals based on their fitness from best to worse;(p)
Store the best solution at generationt;(q)
Apply the proposed adaptive genetic operators to the obtained sequences;(r)
Evaluate the obtained sequences and repair the infeasible ones;(s)
Assign fitness to individuals;(t)
Sort individuals and update the best solution;EndIf the stopping condition is reached, then the algorithm output the best solution, including cycle timeT and its related robot move sequence RN.End


### 3.8.2. Complexity Analysis of the Proposed HQEA

The computational complexity of the proposed HQEA is analyzed as follows. Step (a) and Step (h), respectively, takeO(1) and O(N) times. In step (b), initializing the entire population takes O(2PopsizeN) time. In step (c), determining the Q-bit state for all quantum individuals takes O(Popsize(N-1)) time. In step (d), decoding all the quantum individuals takes O(3PopsizeN(N-1)) time. According to [2], evaluating all the individuals in step (e) takes O(27N6Popsize) time. Besides, the repairing procedure in step (e) takes O(3N2(N-1)RtPopsize)time. The fitness assignment operation in step (f) takes O(3Popsize)time. Sequencing all the individuals in step (g) takes O(Popsizelog⁡Popsize) time. Consequently, the computational complexity for step (a)–step (h) is approximately O(27N6Popsize).In thewhile loop, updating all the quantum individuals in step (j) takes O(PopsizeN) time. Step (q) takes O(2PopsizeN) time. Similar to step (a)–step (h), the computational complexity for other steps in thewhile loop is approximately O(36N6Popsize). Since thewhile loop executesMaxIter times, the computational complexity for thewhile loop is about O(36N6PopsizeMaxIter). Then, the computational complexity is approximately O(36N6PopsizeMaxIter)+O(27N6Popsize), that is, O(9N6Popsize(4MaxIter+3)).


## 3.8.1. Pseudocode of the Proposed HQEA

Based on the above works presented in Sections3.1–3.7, the pseudocode of HQEA for solving the considered problem can be written as follows.Input:
Popsize (population size).MaxIter (maximum number of generations).ω0 (initial rotation angle).Rt (maximum repairing times).N (problem size)Output: the best cycle time T and its related robot move sequence RN.Begin
(a)
t←0;(b)
Initialize Q-bits chromosomes;(c)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(d)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(e)
Evaluate all obtained sequences and repair the infeasible ones;(f)
Assign fitness to individuals;(g)
Sort all individuals based on their fitness from best to worse;(h)
Store the best solution at generationt;While (t<MaxIter)doBegin
(i)
t←t+1;(j)
Update each Q-bits chromosome by using quantum rotation gate;(k)
Determine Q-bit state (i.e.,CN) by Rule 1 and Rule 2;(l)
Apply decoding scheme to generate possible move sequences (i.e.,RN);(m)
Evaluate all obtained sequences and repair the infeasible ones;(n)
Assign fitness to individuals;(o)
Sort all individuals based on their fitness from best to worse;(p)
Store the best solution at generationt;(q)
Apply the proposed adaptive genetic operators to the obtained sequences;(r)
Evaluate the obtained sequences and repair the infeasible ones;(s)
Assign fitness to individuals;(t)
Sort individuals and update the best solution;EndIf the stopping condition is reached, then the algorithm output the best solution, including cycle timeT and its related robot move sequence RN.End


## 3.8.2. Complexity Analysis of the Proposed HQEA

The computational complexity of the proposed HQEA is analyzed as follows. Step (a) and Step (h), respectively, takeO(1) and O(N) times. In step (b), initializing the entire population takes O(2PopsizeN) time. In step (c), determining the Q-bit state for all quantum individuals takes O(Popsize(N-1)) time. In step (d), decoding all the quantum individuals takes O(3PopsizeN(N-1)) time. According to [2], evaluating all the individuals in step (e) takes O(27N6Popsize) time. Besides, the repairing procedure in step (e) takes O(3N2(N-1)RtPopsize)time. The fitness assignment operation in step (f) takes O(3Popsize)time. Sequencing all the individuals in step (g) takes O(Popsizelog⁡Popsize) time. Consequently, the computational complexity for step (a)–step (h) is approximately O(27N6Popsize).In thewhile loop, updating all the quantum individuals in step (j) takes O(PopsizeN) time. Step (q) takes O(2PopsizeN) time. Similar to step (a)–step (h), the computational complexity for other steps in thewhile loop is approximately O(36N6Popsize). Since thewhile loop executesMaxIter times, the computational complexity for thewhile loop is about O(36N6PopsizeMaxIter). Then, the computational complexity is approximately O(36N6PopsizeMaxIter)+O(27N6Popsize), that is, O(9N6Popsize(4MaxIter+3)).


## 4. Experimental Results

To verify the applicability of the proposed algorithm, both benchmark and randomly generated instances were used in the experimental studies. All computational experiments were conducted on an ASUS Laptop with an Intel Core i5-3210M Processor 2.50 GHZ and on a Windows 8 environment. The parameters were set as follows: population size:Popsize=50; the maximum number of generations: MaxIter=200; initial rotation angle ω0=0.05π; μ=0.008. The maximum repairing times: Rt=6. For evaluating the quality of the solution obtained with our HQEA, the problem was also formulated by the MIP approach and solved by the ILOG CPLEX (Version 12.4).


### 4.1. Benchmark Instances

The proposed algorithm was verified by using five well-known benchmark instances in the literature: Mini Phillips (Mini,N=8), Black and Oxide 2 (BO2, N=11), Phillips and Unger (P&U, N=12), and Ligne 1 (N=12) and Ligne 2 (N=14), which can be found in [7, 18, 37–39].Table1 gives the experimental results for five benchmark instances obtained with our algorithm and CPLEX software, in terms of the number of remaining possible CN after applying Rule 1 (Nb. CN after Rule 1 for short), the convergence generation (Con.gen. for short), the best cycle times, and the CPU times (measured in seconds). The “Con.gen.” refers to how many generations are needed for our algorithm to find the best solution and no improvement on the solution in the later iterations. Consequently, the subcolumn “Con. time” represents the time needed by the “Con.gen.” and is computed as: Con.time = Con.gen. × (Our’sCPUtime/MaxIter).Table 1
Results for the benchmark instances.
Instances
Nb.CN after Rule 1
Con.gen.
Best cycle times
CPU times (In seconds)
Our
CPLEX
SD
Our
Con. time
CPLEX
Mini
2
6
2
287
287
0
4.75
0.048
0.16
BO2
210
13
279.3
279.3
0
5.26
0.342
0.25
P&U
2
10
29
521
521
0
5.65
0.819
0.47
Ligne1
211
24
411
392
4.84%
7.35
0.882
0.72
Ligne2
213
26
712
712
0
6.71
0.872
0.48In Table1, we can see that Rule 1 works well on two benchmark instances (i.e., Mini and P&U) as shown in column “Nb. CN after Rule 1,” as the enumerating space of CN is, respectively, reduced 50% for the two instances (note that there are in total 2N-1 individuals for each instance with given value of N.). In column “best cycle times,” our algorithm finds the same solutions as the optimal ones obtained with the CPLEX (see “Our” and “CPLEX”), except for Ligne 1. The standard deviation (SD) of the best cycle time obtained with our algorithm from the optimal cycle time obtained with CPLEX for Ligne 1 is less than 5%. Although the CPU times spent by our algorithm are generally longer than those spent by CPLEX, we can also see in column “Con.gen.” that our algorithm finds the optimal solutions for most cases in very early generations. Note that the time gaps between Con. time and CPLEX are very narrow, less than 1 s. Due to this very small amount of gaps, the difference in CPU times between CPPEX and our algorithm is meaningless and can be negligible. In summary, our algorithm is an effective method for solving the benchmark instances in terms of solution quality and CPU times.


### 4.2. Random Instances

In this subsection, random instances are generated to further test the performance of the proposed algorithm. We compare our algorithm with the QEA with shifting decoding scheme to demonstrate the effectiveness of our decoding scheme. We also compare it with commercial software CPLEX and Tabu search (TS) algorithm [10]. The random instances are generated as follows. We set N which belongs to {10,15,18,20,22}, and let U(c1,c2) be a uniform distribution between parametersc1 andc2. The random tests were set as two different groups. One (called Group 1) was defined as the following way: the time windows were ai=U(30,120) and bi=ai+U(10,750), 1≤i≤N; the empty moves and the loaded moves were, respectively, defined as δi,i+1=U(3,6), δi,j=∑k=ij-1δk,k+1, 0≤i,j≤N+1, and θi=20+δi,i+1,0≤i≤N. The other (called Group 2) was defined as ai=U(40,120), bi=30+U(1,8)×ai, 1≤i≤N, δi,i+1=U(2,5), δi,j=∑k=ij-1δk,k+1, 0≤i,j≤N+1, θi=15+δi,i+1, and 0≤i≤N. These defined parameters were based on the magnitude of the data from real production lines [7, 38]. For each given N, 5 instances were randomly generated.Table2 reports the remaining number of CN for each randomly generated instance after applying Rule 1. As mentioned before, there are in total 2N-1 individuals for each instance with a given value of N. As presented in Table 2, Rule 1 is efficient on 22 random instances (i.e., the numbers in bold font). We can also see in Table 2 that the enumerating space of CN for each instance among the 22 instances is reduced at least 50% and at most 87.5% after applying Rule 1. Based on these results, Rule 1 seems efficient for the studied problem.Table 2
Results for the remaining number ofCN for each instance after applying Rule 1.
N
Group 1
Group 2
1
2
3
4
5
1
2
3
4
5
10
29
2
8
29
29
29
29
29
29
29
29
15
214
214
214
2
13
214
214
214
214
214
214
18
2
16
217
217
2
16
217
2
16
217
217
2
16
2
15
20
2
18
219
219
2
17
2
18
2
18
2
18
2
18
2
18
219
22
2
18
221
2
20
2
20
2
18
2
19
2
20
221
2
19
2
19Firstly, we compare our algorithm with the QEA with shifting decoding scheme (i.e., SQEA). Table3 presents the comparison results between our decoding scheme and shifting decoding scheme on Group 1 and Group 2. We can see that our decoding scheme generally outperforms the shifting decoding scheme in terms of solution quality and CPU times for all random instances. In particular, the deviations (i.e., AD) of our algorithm from that with shifting decoding generally decrease with the problem size. Besides, our algorithm spent less time than that with shifting decoding for all random instances.Table 3
Comparison results between our decoding scheme and shifting decoding scheme on Group 1 and Group 2.
N
Group 1
Group 2
Average cycle times
Average CPU times
Average cycle times
Average CPU times
Our
SQEA
AD
Our
SQEA
Our
SQEA
AD
Our
SQEA
10
400.4
401.2
−0.20%
6.74
10.12
318.4
318.4
0
6.83
17.8
15
607.2
628
−3.31%
24.56
51.79
470.6
470.8
−0.04%
37.44
146.68
18
808.8
817.4
−1.05%
54.88
286.55
627.4
638.2
−1.69%
49.46
267.57
20
897.2
927.8
−3.30%
117.53
360.14
678.6
690.2
−1.68%
141.02
275.59
22
1058.6
1351.2
−21.65%
274.43
315.62
802.6
878.2
−8.61%
190.16
373.43Secondly, Tables4 and 5 report the comparison results for randomly generated instances using our algorithm, Yan’s algorithm, and commercial optimization software CPLEX. Columns AD1 and AD2, respectively, represent the standard deviation of our solution from those obtained with CPLEX and Yan’s algorithm; that is, AD1=(Our-CPLEX)/CPLEX×100% and AD2=(Our-Yan)/Yan×100%. As presented in Tables 4 and 5, our algorithm and Yan’s algorithm find the same solutions as the optimal ones obtained with CPLEX for random instances with N=10. For the remaining random instances, the average cycle times obtained with our algorithm are smaller than those obtained with Yan’s algorithm. As a result, the deviations (i.e., AD2) of our algorithm from Yan’s algorithm are all negative, which range from −5.89% to −1.9% in Table 4 and from −3.93% to −0.42% in Table 5. Note that the smaller the AD2, the better the solution quality our algorithm obtained over Yan’s algorithm. Therefore, our algorithm has a better solution quality than Yan’s algorithm. We also notice that CPLEX has a better solution quality than our algorithm and Yan’s algorithm but it spent much longer CPU times, which will be discussed later. Moreover, the values of AD1 in Tables 4 and 5 both increase with the problem size but are less than 4% and 3%, respectively, which are generally small and acceptable.Table 4
Comparison results for the randomly generated instances Group 1.
N
Average cycle times
Average CPU times (In seconds)
Our
Yan
CPLEX
AD1
AD2
Our
Yan
CPLEX
10
400.4
400.4
400.4
0
0
6.74
2.7
1.44
15
607.2
624.6
602.4
0.8%
−2.79%
24.56
19.95
42.95
18
808.8
859.4
797.6
1.4%
−5.89%
54.88
32.16
1351.53
20
897.2
914.6
865.8
3.63%
−1.90%
117.53
114.51
1692.12
22
1058.6
1122.4
1025
3.28%
−5.68%
274.43
211.34
2712.38Table 5
Comparison results for the randomly generated instances Group 2.
N
Average cycle times
Average CPU times (In seconds)
Our
Yan
CPLEX
AD1
AD2
Our
Yan
CPLEX
10
318.4
318.4
318.4
0
0
6.83
4.58
1.38
15
470.6
472.6
466.4
0.9%
−0.42%
37.44
62.35
51.50
18
627.4
636.4
612.6
2.42%
−1.41%
49.46
92.84
324.24
20
678.6
684
661.8
2.54%
−0.79%
141.02
53.52
1077.9
22
802.6
835.4
779.8
2.92%
−3.93%
190.16
102.62
1897.76For the average CPU times, we can see from Tables4 and 5 that both our algorithm and Yan’s algorithm perform much better than CPLEX for each value of N, except for N=10. We also notice that Yan’s algorithm has a better performance than our algorithm in terms of CPU times except for N=15 and N=18 in Table 5. But their gaps are not so large. Moreover, although the CPU times spent by the three approaches generally increase with the instance size N, the CPU times spent by CPLEX generally have a very sharper growth than those spent by our algorithm and Yan’s algorithm, especially for large-size instances. From these results, we can see that our proposed algorithm has a better computational performance than CPLEX.Furthermore, Table6 reports the comparison results of our proposed algorithm with and without a repairing procedure. From it, we can see that our algorithm without a repairing procedure (denoted by Our#) generally performs worse on instances Group 1 and Group 2 than that with a repairing procedure (denoted by Our), in terms of solution quality. From column “AD,” we can see that when problem size (i.e., N) increases, our algorithm has a much better solution than that without a repairing procedure, although it costs relatively much more CPU time. In other words, our designed repairing procedure for infeasible solutions is effective for the studied problem. Based on the above results and analysis, we can conclude that our proposed algorithm works well and is effective to solve the studied scheduling problem.Table 6
Comparison results of our proposed algorithm with and without repairing procedure on randomly generated instances Group 1 and Group 2.
N
Group 1
Group 2
Average cycle times
Average CPU times
Average cycle times
Average CPU times
Our
Our#
AD
Our
Our#
Our
Our#
AD
Our
Our#
10
400.4
406.8
−1.57%
6.74
3.94
318.4
320
−0.50%
6.83
4.09
15
607.2
642.6
−5.51%
24.56
9.86
470.6
518.4
−9.22%
37.44
10.77
18
808.8
966.8
−16.34%
54.88
17.33
627.4
722.2
−13.07%
49.46
17.85
20
897.2
1058.2
−15.21%
117.53
36.31
678.6
1095.8
−38.07%
141.02
33.27
22
1058.6
1400
−24.39%
274.43
87.25
802.6
1199.2
−33.07%
190.16
63.25
Note. Our# represents our algorithm without repairing procedure; AD=(Our-Our#)/Our#×100%.


## 4.1. Benchmark Instances

The proposed algorithm was verified by using five well-known benchmark instances in the literature: Mini Phillips (Mini,N=8), Black and Oxide 2 (BO2, N=11), Phillips and Unger (P&U, N=12), and Ligne 1 (N=12) and Ligne 2 (N=14), which can be found in [7, 18, 37–39].Table1 gives the experimental results for five benchmark instances obtained with our algorithm and CPLEX software, in terms of the number of remaining possible CN after applying Rule 1 (Nb. CN after Rule 1 for short), the convergence generation (Con.gen. for short), the best cycle times, and the CPU times (measured in seconds). The “Con.gen.” refers to how many generations are needed for our algorithm to find the best solution and no improvement on the solution in the later iterations. Consequently, the subcolumn “Con. time” represents the time needed by the “Con.gen.” and is computed as: Con.time = Con.gen. × (Our’sCPUtime/MaxIter).Table 1
Results for the benchmark instances.
Instances
Nb.CN after Rule 1
Con.gen.
Best cycle times
CPU times (In seconds)
Our
CPLEX
SD
Our
Con. time
CPLEX
Mini
2
6
2
287
287
0
4.75
0.048
0.16
BO2
210
13
279.3
279.3
0
5.26
0.342
0.25
P&U
2
10
29
521
521
0
5.65
0.819
0.47
Ligne1
211
24
411
392
4.84%
7.35
0.882
0.72
Ligne2
213
26
712
712
0
6.71
0.872
0.48In Table1, we can see that Rule 1 works well on two benchmark instances (i.e., Mini and P&U) as shown in column “Nb. CN after Rule 1,” as the enumerating space of CN is, respectively, reduced 50% for the two instances (note that there are in total 2N-1 individuals for each instance with given value of N.). In column “best cycle times,” our algorithm finds the same solutions as the optimal ones obtained with the CPLEX (see “Our” and “CPLEX”), except for Ligne 1. The standard deviation (SD) of the best cycle time obtained with our algorithm from the optimal cycle time obtained with CPLEX for Ligne 1 is less than 5%. Although the CPU times spent by our algorithm are generally longer than those spent by CPLEX, we can also see in column “Con.gen.” that our algorithm finds the optimal solutions for most cases in very early generations. Note that the time gaps between Con. time and CPLEX are very narrow, less than 1 s. Due to this very small amount of gaps, the difference in CPU times between CPPEX and our algorithm is meaningless and can be negligible. In summary, our algorithm is an effective method for solving the benchmark instances in terms of solution quality and CPU times.


## 4.2. Random Instances

In this subsection, random instances are generated to further test the performance of the proposed algorithm. We compare our algorithm with the QEA with shifting decoding scheme to demonstrate the effectiveness of our decoding scheme. We also compare it with commercial software CPLEX and Tabu search (TS) algorithm [10]. The random instances are generated as follows. We set N which belongs to {10,15,18,20,22}, and let U(c1,c2) be a uniform distribution between parametersc1 andc2. The random tests were set as two different groups. One (called Group 1) was defined as the following way: the time windows were ai=U(30,120) and bi=ai+U(10,750), 1≤i≤N; the empty moves and the loaded moves were, respectively, defined as δi,i+1=U(3,6), δi,j=∑k=ij-1δk,k+1, 0≤i,j≤N+1, and θi=20+δi,i+1,0≤i≤N. The other (called Group 2) was defined as ai=U(40,120), bi=30+U(1,8)×ai, 1≤i≤N, δi,i+1=U(2,5), δi,j=∑k=ij-1δk,k+1, 0≤i,j≤N+1, θi=15+δi,i+1, and 0≤i≤N. These defined parameters were based on the magnitude of the data from real production lines [7, 38]. For each given N, 5 instances were randomly generated.Table2 reports the remaining number of CN for each randomly generated instance after applying Rule 1. As mentioned before, there are in total 2N-1 individuals for each instance with a given value of N. As presented in Table 2, Rule 1 is efficient on 22 random instances (i.e., the numbers in bold font). We can also see in Table 2 that the enumerating space of CN for each instance among the 22 instances is reduced at least 50% and at most 87.5% after applying Rule 1. Based on these results, Rule 1 seems efficient for the studied problem.Table 2
Results for the remaining number ofCN for each instance after applying Rule 1.
N
Group 1
Group 2
1
2
3
4
5
1
2
3
4
5
10
29
2
8
29
29
29
29
29
29
29
29
15
214
214
214
2
13
214
214
214
214
214
214
18
2
16
217
217
2
16
217
2
16
217
217
2
16
2
15
20
2
18
219
219
2
17
2
18
2
18
2
18
2
18
2
18
219
22
2
18
221
2
20
2
20
2
18
2
19
2
20
221
2
19
2
19Firstly, we compare our algorithm with the QEA with shifting decoding scheme (i.e., SQEA). Table3 presents the comparison results between our decoding scheme and shifting decoding scheme on Group 1 and Group 2. We can see that our decoding scheme generally outperforms the shifting decoding scheme in terms of solution quality and CPU times for all random instances. In particular, the deviations (i.e., AD) of our algorithm from that with shifting decoding generally decrease with the problem size. Besides, our algorithm spent less time than that with shifting decoding for all random instances.Table 3
Comparison results between our decoding scheme and shifting decoding scheme on Group 1 and Group 2.
N
Group 1
Group 2
Average cycle times
Average CPU times
Average cycle times
Average CPU times
Our
SQEA
AD
Our
SQEA
Our
SQEA
AD
Our
SQEA
10
400.4
401.2
−0.20%
6.74
10.12
318.4
318.4
0
6.83
17.8
15
607.2
628
−3.31%
24.56
51.79
470.6
470.8
−0.04%
37.44
146.68
18
808.8
817.4
−1.05%
54.88
286.55
627.4
638.2
−1.69%
49.46
267.57
20
897.2
927.8
−3.30%
117.53
360.14
678.6
690.2
−1.68%
141.02
275.59
22
1058.6
1351.2
−21.65%
274.43
315.62
802.6
878.2
−8.61%
190.16
373.43Secondly, Tables4 and 5 report the comparison results for randomly generated instances using our algorithm, Yan’s algorithm, and commercial optimization software CPLEX. Columns AD1 and AD2, respectively, represent the standard deviation of our solution from those obtained with CPLEX and Yan’s algorithm; that is, AD1=(Our-CPLEX)/CPLEX×100% and AD2=(Our-Yan)/Yan×100%. As presented in Tables 4 and 5, our algorithm and Yan’s algorithm find the same solutions as the optimal ones obtained with CPLEX for random instances with N=10. For the remaining random instances, the average cycle times obtained with our algorithm are smaller than those obtained with Yan’s algorithm. As a result, the deviations (i.e., AD2) of our algorithm from Yan’s algorithm are all negative, which range from −5.89% to −1.9% in Table 4 and from −3.93% to −0.42% in Table 5. Note that the smaller the AD2, the better the solution quality our algorithm obtained over Yan’s algorithm. Therefore, our algorithm has a better solution quality than Yan’s algorithm. We also notice that CPLEX has a better solution quality than our algorithm and Yan’s algorithm but it spent much longer CPU times, which will be discussed later. Moreover, the values of AD1 in Tables 4 and 5 both increase with the problem size but are less than 4% and 3%, respectively, which are generally small and acceptable.Table 4
Comparison results for the randomly generated instances Group 1.
N
Average cycle times
Average CPU times (In seconds)
Our
Yan
CPLEX
AD1
AD2
Our
Yan
CPLEX
10
400.4
400.4
400.4
0
0
6.74
2.7
1.44
15
607.2
624.6
602.4
0.8%
−2.79%
24.56
19.95
42.95
18
808.8
859.4
797.6
1.4%
−5.89%
54.88
32.16
1351.53
20
897.2
914.6
865.8
3.63%
−1.90%
117.53
114.51
1692.12
22
1058.6
1122.4
1025
3.28%
−5.68%
274.43
211.34
2712.38Table 5
Comparison results for the randomly generated instances Group 2.
N
Average cycle times
Average CPU times (In seconds)
Our
Yan
CPLEX
AD1
AD2
Our
Yan
CPLEX
10
318.4
318.4
318.4
0
0
6.83
4.58
1.38
15
470.6
472.6
466.4
0.9%
−0.42%
37.44
62.35
51.50
18
627.4
636.4
612.6
2.42%
−1.41%
49.46
92.84
324.24
20
678.6
684
661.8
2.54%
−0.79%
141.02
53.52
1077.9
22
802.6
835.4
779.8
2.92%
−3.93%
190.16
102.62
1897.76For the average CPU times, we can see from Tables4 and 5 that both our algorithm and Yan’s algorithm perform much better than CPLEX for each value of N, except for N=10. We also notice that Yan’s algorithm has a better performance than our algorithm in terms of CPU times except for N=15 and N=18 in Table 5. But their gaps are not so large. Moreover, although the CPU times spent by the three approaches generally increase with the instance size N, the CPU times spent by CPLEX generally have a very sharper growth than those spent by our algorithm and Yan’s algorithm, especially for large-size instances. From these results, we can see that our proposed algorithm has a better computational performance than CPLEX.Furthermore, Table6 reports the comparison results of our proposed algorithm with and without a repairing procedure. From it, we can see that our algorithm without a repairing procedure (denoted by Our#) generally performs worse on instances Group 1 and Group 2 than that with a repairing procedure (denoted by Our), in terms of solution quality. From column “AD,” we can see that when problem size (i.e., N) increases, our algorithm has a much better solution than that without a repairing procedure, although it costs relatively much more CPU time. In other words, our designed repairing procedure for infeasible solutions is effective for the studied problem. Based on the above results and analysis, we can conclude that our proposed algorithm works well and is effective to solve the studied scheduling problem.Table 6
Comparison results of our proposed algorithm with and without repairing procedure on randomly generated instances Group 1 and Group 2.
N
Group 1
Group 2
Average cycle times
Average CPU times
Average cycle times
Average CPU times
Our
Our#
AD
Our
Our#
Our
Our#
AD
Our
Our#
10
400.4
406.8
−1.57%
6.74
3.94
318.4
320
−0.50%
6.83
4.09
15
607.2
642.6
−5.51%
24.56
9.86
470.6
518.4
−9.22%
37.44
10.77
18
808.8
966.8
−16.34%
54.88
17.33
627.4
722.2
−13.07%
49.46
17.85
20
897.2
1058.2
−15.21%
117.53
36.31
678.6
1095.8
−38.07%
141.02
33.27
22
1058.6
1400
−24.39%
274.43
87.25
802.6
1199.2
−33.07%
190.16
63.25
Note. Our# represents our algorithm without repairing procedure; AD=(Our-Our#)/Our#×100%.


## 5. Conclusions

This paper proposed a hybrid QEA with improved decoding scheme to solve a single robot one-degree cyclic scheduling problem with flexible processing times in a robotic flow shop. In particular, three different decoding procedures were proposed to convert Q-bit individual into robot move sequences. A repairing procedure was designed to repair the infeasible sequence. Both Q-gate and adaptive genetic operators as variant operators were applied to evolve the population. The effectiveness of the proposed algorithm was demonstrated by solving benchmark instances and randomly generated instances compared with commercial software CPLEX and Yan’s algorithm. Experimental results indicate that our decoding scheme outperforms the shifting decoding scheme, and the proposed algorithm can provide high-quality solutions within a reasonable time. The results also imply that the proposed algorithm generally has a shorter computation time than CPLEX, especially for large-size instances and has a better solution quality than Yan’s algorithm.For further research direction, the hybrid QEA algorithm developed in this study can be extended to the robotic scheduling problem with multiobjective optimization (such as simultaneously optimizing of cycle time and production cost). This extension can be accomplished by applying the well-known Pareto-based ranking approach to evaluate each individual. In further research, we may extend the proposed algorithm to solve the robotic scheduling problem in a multidegree context.



---

*Source: 3064724-2017-04-16.xml*