Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 3064724, 13 pages
https://doi.org/10.1155/2017/3064724

Research Article

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

Weidong Lei,¹ Hervé Manier,² Marié-Ange Manier,² and Xinping Wang¹

¹School of Management, Xi’an University of Science and Technology, Xi’an 710054, China
²Université Bourgogne Franche-Comté, UTBM, OPERA, 90100 Belfort, France

Correspondence should be addressed to Weidong Lei; leiweidong@xust.edu.cn

Received 8 October 2016; Accepted 28 February 2017; Published 16 April 2017

Academic Editor: Calogero Orlando

Copyright © 2017 Weidong Lei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 Figure 1. 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 with k-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 Section 2, 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 processing machines (i.e., ) in a production line and a single robot for part transportation, machine 0 (i.e., ) and machine (i.e., ) are the input station and the output station, respectively. After a part is unloaded from , it is to be successively processed through to . The robot moving a part from to () is called (loaded) move , which includes three suboperations: () unloading a part from a machine; () carrying the part to the next machine; () 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]:

: the minimum and maximum bounds of the processing time on , respectively, .: the time needed to perform move , .: the travel time for empty robot from to ; note that and , . Note that satisfies the triangle inequalities; that is, , , , .

The decision variables are the following ones:

: cycle time. It is the duration of a cycle.: the start time of (loaded) move within a cycle, . Without loss of generality, move 0 is supposed to be the first move of a cycle; thus .

To facilitate the formulation, we define the following intermediate variables:

: if , then is empty at the beginning of a cycle; else ; then is occupied by a part, . Define . It is called the initial part distribution at the beginning of a cycle. Without loss of generality, we assume and , since is always occupied by part and move 0 is the first move of a cycle.: the th move performed by the robot within a cycle, . As mentioned above, we have . Define , which represents the sequence of moves during a cycle. An example of with is , where , , and , as shown in Figure 2. Here, means that the second move transfers a part from to .

Figure 2: An example of cyclic schedule with a single robot.

Figure 2 shows an illustrative example of the studied problem with . In this example, there are three processing machines (i.e., , , and ) with a single robot for part transportation as well as the loading station (i.e., ) and the unloading station (i.e., ). 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, is the only machine to be occupied (and implicitly ). So the initial part distribution is . For this distribution, the optimal robot move sequence is (i.e., ). When move 1 finishes, the robot comes back to 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 as those ones in time interval . This is called cyclic production mode. The duration of the repetitive sequence (i.e., ) is the cycle time .

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

In Figure 2, we notice that the start time of processing operation is the same as the end time of loaded move (i.e., ); the end time of processing operation is the same as the start time of loaded move (i.e., ). 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 can be represented as and , for and , respectively. Considering the processing time requirements, we have

Furthermore, the robot must have enough time to perform any two successive moves (i.e., and ); thus, the following relation holds:

When the last move (i.e., ) finishes, the robot must come back to for executing move 0 of the next cycle. Hence, we have

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]:

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]:

In (6), and are two complex numbers, which represent the probability amplitudes of states 0 and 1, respectively. As a result, and 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 Q-bits [as shown below in (7)]; this individual can represent 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 Section 3.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 between moves in this article, we let Q-bit corresponding to machine , for , and use Rule and Rule introduced in the following section to determine each Q-bit state. If Q-bit is in state “0” (i.e., ), move is performed before move during a cycle; otherwise (i.e., ), move is performed before move during a cycle. Hence, an individual containing Q-bits is used to represent machine states and is defined as follows:

where , . 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 , more precisely the state of each (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 gene gives the index of the machine from which the th move starts during one cycle. In such representations, the solution space contains individuals. With our representation, we handle in a first step only individuals. This number may be further reduced for some instances with Rule , 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 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 that ; therefore move occurs before move within a cycle. Moreover, let us suppose that move and move are the last move and the second move of a cycle, respectively. Correspondingly, the minimum processing time of a part on machine with is Indeed, move would be the first move to be performed after move 0, and move would be the last move of the cycle. Else, the processing time on would be greater than , which would make the following assertion even more true. Then we can compare this processing time with which is the maximum authorized time on .

(a1)If (hereafter called Rule ) happens, then we can know that the processing time requirement on is violated. Consequently, all sequences relevant with are infeasible ones. So must be 0.(a2)Else, may be 0 or 1.

Take move 1 and move 2 in Figure 2, for example, the processing time of a part on is equal to Rule can be used to reduce the enumerating space of and thus improve the search efficiency. Indeed, if Rule enables us to fix 0 to the values of variables , then the search space of can be reduced to individuals.

For the state of Q-bit in a quantum individual that is not determined by Rule , a random number is generated from the uniform distribution [0, 1). If , then Q-bit is in state “1” (i.e., ); else, Q-bit is in state “0” (i.e., ). This method is called Rule . Based on the above, the states of all Q-bits in one individual can be easily determined by Rule and Rule , that is to say, the initial part distribution .

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 sequence , providing that the states of all Q-bits (i.e., ) in a quantum individual are already determined.

3.4.1. Decoding Procedure 1

For ease of description, we first define as a copy of and . Let be a set that records the performed moves. It should be note that can be seen as an indicator that indicates the state (i.e., empty or occupied) of in the process. Thus, the value of is dynamically modified in the process. That is, when move finishes, both the states of and are changed; that is, becomes empty and is occupied by a part. Thus, we set , and put move into set Φ.

Procedure mainly depends on the probability sizes of Q-bits in to derive the move sequence. In particular, for each given , when move finishes, for , we first determine the possible candidates for next move ; if move satisfies the condition of and (note that if , the output station can be seen as always empty) and , then we set Notice that we use set Ω to store all possible candidates for the next move r[k+1]. Finally, we choose move with the highest probability (i.e., ) in set Ω as move , and let , (for ), , and In the next step, we update all related parameters and use a similar way to derive the following move (i.e., ). Once the whole sequence (i.e., ) is determined, this procedure stops. Pseudocode 1 illustrated the process of Decoding Procedure .

Input:  

Begin

                         ;

            ; ; ;

            ; ; ;

                    if () then

                                ;

                    find  ;

                    , ;

                    if ()   then  ;

                    ;

                            if  

                                     then  ;

End

Output:  

Pseudocode 1: Pseudocode of Decoding Procedure .

For example, a complete part distribution (corresponding to a quantum individual) with is When the first move (i.e., ) finishes, by using the above decoding procedure, we have , , , , , and , from which we can know that , , and 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 Finally, according to the selection rule, if , we have else, The similar ways are used to update , Φ, Ω and then determine ,

3.4.2. Decoding Procedure 2

For ease of description, we keep the intermediate parameters , , and Ω defined in Procedure . Furthermore, we define as the start time of move in the process of deriving the whole sequence, for , and let To derive a move sequence from each given , Procedure mainly depends on the rule of minimal time unit increment between and , for , while respecting the processing time requirement, since the objective of the problem is to minimize the cycle time . 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 from the set.

In particular, on one hand, when move finishes, as similarly done in Procedure , we derive the value of Ω from each given . On the other hand, we design a different strategy to determine move compared with the Procedure . At first, we calculate each in set Ω; that is, Then, for each move in set Ω, we check whether move exists in the partial sequence . If it exists and happens, then we update so as to meet the processing requirement. Finally, we choose move with the smallest value of in set Ω as move , and let , (for ), , and In the next step, we update all related parameters so as to derive move . Once the whole sequence (i.e., ) is determined, this procedure stops. Pseudocode 2 illustrated the process of Decoding Procedure .

Input:

Begin

            ; ; ; ;

            for    to  

                        ;

            ; ; ;

            for    to  

                 if () then

                                    ;

            for    to  

                 compute  , ;

                    for     to   ,

                          then  ;

                 find  , ;

                    ;

                 if () then  ;

                  ;

                 for    to  

                      if

                             then  ;

End

Output:  

Pseudocode 2: Pseudocode of Decoding Procedure .

For instance, an example of with is When the first move (i.e., ) finishes, by definitions, we have , , , , , and as well as We first calculate (note that if , then we set ) and then and by plus , , , respectively, and then choose the move with the smallest starting time among the three candidates as . The similar ways are used to update , Φ, Ω and then determine ,

3.4.3. Decoding Procedure 3

Procedure mainly depends on the precedence relationship between move and move (i.e., the value of ) to derive the move sequence. For each given (i.e., a quantum individual), if , then move is set before move in ; else, move is set after move in . Pseudocode 3 illustrated the process of Decoding Procedure .

Input:  ,

Begin

            for    to  

                    , , , , ;

                    for    to  

                            if ()  then  , ;

                        if ()  then  , ;

                            else if () break;

                    if () then

                                    if () then

                                        , and  for    to  

                                            ;

                    else

                                    if () then

                                        , and  for    to  

                                        ;

End

Output:

Pseudocode 3: Pseudocode of Decoding Procedure .

For instance, examples of and with are, respectively, and , from which we can easily derive a new sequence; that is, Note that, at the initial step, we set ,

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, is defined to represent the fitness value of each individual . Basically the higher the fitness of an individual (i.e., ), the better the quality of individual . Therefore, is computed as follows: , where is a positive parameter and set as 1 (note that in fact 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 (), 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 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 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 sequence , we identify each partial sequence in a whole robot move sequence which is either in sequence of (which means move is performed before move within a cycle) or of (which means move is performed before move within a cycle). That is to say, a complete robot sequence consists of pieces of such a partial sequence. For ease of description, we define the following parameters:

: the duration between the finish time of move and the start time of move for a partial sequence , for Note that 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 and in the sequence and , then we let . Note that may span the cycle or be within a cycle. For example, in Figure 2, the two consecutive sequences are From it, we can see that and are within a cycle, but spans the cycle. Therefore, can be used to check the satisfaction of flexible processing requirements no matter or .: the mean time of all loaded move times, .: the mean time of all empty move times, . Note that the possible number of empty moves is . Since the empty moves between and , and , and and do not actually happen, the number is reduced to .

For an infeasible sequence , we first use parameters to check the sequence .

(1)If is verified to be greater than its upper bound , then we remove one or more moves from the corresponding partial sequence, so as to make the partial sequence feasible; else if is verified to be smaller than its lower bound and the time gap between and is greater than the sum of and , 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., , ) given by the evaluation process. For ease of description, let moves and be the identified two moves violating the robot capacity constraints; that is, , with If these two moves are two consecutive moves, we set move before move in sequence so as to make the sequence be feasible; else, we remove one or more moves between moves and 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 gate ) is adopted as the variation operator to update the Q-bits in (7). is set to be as the initial rotation angle. For individual , the Q-bit in it can be updated as the following way [27, 30]:

We define fit_ as the fitness of the best individual found in population. The rotation angle is defined according to the respective values of the corresponding parameter in the individual (labeled ) and in the best one (labeled ). If the condition holds, then consider the following conditions [27].

Case 1. If bit is in the first or the third quadrant, then consider the following.

Case 1.1. If and , then ; here the rotation angle is set negative so as to increase the probability that Q-bit is in state “1.”

Case 1.2. If and , then , and the rotation angle is set as positive so as to increase the probability that Q-bit is in state “0.”

Case 1.3. Else, .

Case 2. If bit is in the second or the fourth quadrant, then consider the following.

Case 2.1. If and , then ; here the rotation angle is set as positive so as to increase the probability that Q-bit is in state “1.”

Case 2.2. If and , then , and the rotation angle is set as negative so as to increase the probability that Q-bit is in state “0.”

Case 2.3. Else, ;

Besides, since the probability of a Q-bit in state “0” may be equal to 1 or 0, the updated Q-bit 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 . As a result, the following equation must be considered [34]:

By applying the decoding procedures given in Section 3.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:

, : 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: . is the lower bound on cycle time for the instance. It can be obtained by the following way, which is taken from Chen et al. [2]:

According to Srinivas and Patnaik [36], and are defined, respectively, in a similar way:

Adaptively adjusting and (i.e., (11)) can prevent divergence and escape from the local optimal, since (11) can dynamically reduce and for individuals with high fitness or increase and 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 positions and are randomly chosen (), For and , the values of for the new offspring 1 are inherited from parent 1. For , the new 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.

(a) Two-point crossover operation

(b) Mutation operation

Figure 3: Crossover and mutation operators.

Besides, a mutation operator is adopted to prevent a solution falling into a local optimum, which is designed as follows. For a chosen individual , first, we randomly choose a position p, and then reorder the partial move sequence . 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 Sections 3.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). (initial rotation angle). (maximum repairing times). (problem size)

Output: the best cycle time and its related robot move sequence .Begin

(a);(b)Initialize Q-bits chromosomes;(c)Determine Q-bit state (i.e., ) by Rule and Rule ;(d)Apply decoding scheme to generate possible move sequences (i.e., );(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 generation ;While () doBegin

(i);(j)Update each Q-bits chromosome by using quantum rotation gate;(k)Determine Q-bit state (i.e., ) by Rule and Rule ;(l)Apply decoding scheme to generate possible move sequences (i.e., );(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 generation ;(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 time and its related robot move sequence .

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, take and times. In step (b), initializing the entire population takes time. In step (c), determining the Q-bit state for all quantum individuals takes time. In step (d), decoding all the quantum individuals takes time. According to [2], evaluating all the individuals in step (e) takes time. Besides, the repairing procedure in step (e) takes time. The fitness assignment operation in step (f) takes time. Sequencing all the individuals in step (g) takes time. Consequently, the computational complexity for step (a)–step (h) is approximately .

In the while loop, updating all the quantum individuals in step (j) takes time. Step (q) takes time. Similar to step (a)–step (h), the computational complexity for other steps in the while loop is approximately . Since the while loop executes MaxIter times, the computational complexity for the while loop is about . Then, the computational complexity is approximately , that is, .

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: ; the maximum number of generations: ; initial rotation angle ; The maximum repairing times: 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, ), Black and Oxide 2 (BO2, ), Phillips and Unger (P&U, ), and Ligne 1 () and Ligne 2 (), which can be found in [7, 18, 37–39].

Table 1 gives the experimental results for five benchmark instances obtained with our algorithm and CPLEX software, in terms of the number of remaining possible after applying Rule (Nb. after Rule 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. after Rule
Con.gen.
Best cycle times
CPU times (In seconds)

Our
CPLEX
SD
Our
Con. time
CPLEX

Mini
2 ⁶
2
287
287
0
4.75
0.048
0.16

BO2
2¹⁰
13
279.3
279.3
0
5.26
0.342
0.25

P&U
2 ¹⁰
29
521
521
0
5.65
0.819
0.47

Ligne
2¹¹
24
411
392
4.84%
7.35
0.882
0.72

Ligne
2¹³
26
712
712
0
6.71
0.872
0.48

In Table 1, we can see that Rule works well on two benchmark instances (i.e., Mini and P&U) as shown in column “Nb. after Rule ,” as the enumerating space of is, respectively, reduced 50% for the two instances (note that there are in total individuals for each instance with given value of .). 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 which belongs to , and let be a uniform distribution between parameters c1 and c2. The random tests were set as two different groups. One (called Group 1) was defined as the following way: the time windows were and , ; the empty moves and the loaded moves were, respectively, defined as , , , and , The other (called Group 2) was defined as , , , , , , , and These defined parameters were based on the magnitude of the data from real production lines [7, 38]. For each given , 5 instances were randomly generated.

Table 2 reports the remaining number of for each randomly generated instance after applying Rule . As mentioned before, there are in total individuals for each instance with a given value of . As presented in Table 2, Rule 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 for each instance among the 22 instances is reduced at least 50% and at most 87.5% after applying Rule . Based on these results, Rule seems efficient for the studied problem.

Table 2: Results for the remaining number of for each instance after applying Rule .

Group 1
Group 2

1
2
3
4
5
1
2
3
4
5

10
2⁹
2 ⁸
2⁹
2⁹
2⁹
2⁹
2⁹
2⁹
2⁹
2⁹

15
2¹⁴
2¹⁴
2¹⁴
2 ¹³
2¹⁴
2¹⁴
2¹⁴
2¹⁴
2¹⁴
2¹⁴

18
2 ¹⁶
2¹⁷
2¹⁷
2 ¹⁶
2¹⁷
2 ¹⁶
2¹⁷
2¹⁷
2 ¹⁶
2 ¹⁵

20
2 ¹⁸
2¹⁹
2¹⁹
2 ¹⁷
2 ¹⁸
2 ¹⁸
2 ¹⁸
2 ¹⁸
2 ¹⁸
2¹⁹

22
2 ¹⁸
2²¹
2 ²⁰
2 ²⁰
2 ¹⁸
2 ¹⁹
2 ²⁰
2²¹
2 ¹⁹
2 ¹⁹

Firstly, we compare our algorithm with the QEA with shifting decoding scheme (i.e., SQEA). Table 3 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.

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.43

Secondly, Tables 4 and 5 report the comparison results for randomly generated instances using our algorithm, Yan’s algorithm, and commercial optimization software CPLEX. Columns AD¹ and AD², respectively, represent the standard deviation of our solution from those obtained with CPLEX and Yan’s algorithm; that is, and 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 . 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., AD²) 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 AD², 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 AD¹ 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.

Average cycle times
Average CPU times (In seconds)

Our
Yan
CPLEX
AD¹
AD²
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.38

Table 5: Comparison results for the randomly generated instances Group 2.

Average cycle times
Average CPU times (In seconds)

Our
Yan
CPLEX
AD¹
AD²
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.76

For the average CPU times, we can see from Tables 4 and 5 that both our algorithm and Yan’s algorithm perform much better than CPLEX for each value of , except for . We also notice that Yan’s algorithm has a better performance than our algorithm in terms of CPU times except for and 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 , 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, Table 6 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., ) 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.

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;

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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was partially supported by the scientific research program funded by Xi’an University of Science and Technology (2015QDJ052, 201749), the National Natural Science Foundation of China (71401133, 71601154), and China Scholarship Council.

References

Q. Zhang, H. Manier, M.-A. Manier, and C. Bloch, “Heuristics for predictive hoist scheduling problems,” European Journal of Industrial Engineering, vol. 8, no. 5, pp. 695–715, 2014.

H. Chen, C. Chu, and J.-M. Proth, “Cyclic scheduling of a hoist with time window constraints,” IEEE Transactions on Robotics and Automation, vol. 14, no. 1, pp. 144–152, 1998.

M.-A. Manier and C. Bloch, “A classification for hoist scheduling problems,” International Journal of Flexible Manufacturing Systems, vol. 15, no. 1, pp. 37–55, 2003.

P. Yan, A. Che, X. Cai, and X. Tang, “Two-phase branch and bound algorithm for robotic cells rescheduling considering limited disturbance,” Computers and Operations Research, vol. 50, pp. 128–140, 2014.

J. Feng, A. Che, and N. Wang, “Bi-objective cyclic scheduling in a robotic cell with processing time windows and non-Euclidean travel times,” International Journal of Production Research, vol. 52, no. 9, pp. 2505–2518, 2014.

L. Lei and T. Wang, “A proof: the cyclic hoist scheduling problem is NP-complete,” Working Paper, Rutgers University, 1989.

L. W. Phillips and P. S. Unger, “Mathematical programming solution of a hoist scheduling program,” AIIE Transactions, vol. 8, no. 2, pp. 219–225, 1976.

A. Che, J. Feng, H. Chen, and C. Chu, “Robust optimization for the cyclic hoist scheduling problem,” European Journal of Operational Research, vol. 240, no. 3, pp. 627–636, 2015.

L. Lei and T. J. Wang, “Determining optimal cyclic hoist schedules in a single-hoist electroplating line,” IIE Transactions, vol. 26, no. 2, pp. 25–33, 1994.

P. Yan, A. Che, N. Yang, and C. Chu, “A tabu search algorithm with solution space partition and repairing procedure for cyclic robotic cell scheduling problem,” International Journal of Production Research, vol. 50, no. 22, pp. 6403–6418, 2012.

J. Liu, Y. Jiang, and Z. Zhou, “Cyclic sduling of a single hoist in extended electroplating lines: a comprehensive integer programming solution,” IIE Transactions, vol. 34, no. 10, pp. 905–914, 2002.

A. Che and C. Chu, “Cyclic hoist scheduling in large real-life electroplating lines,” OR Spectrum, vol. 29, no. 3, pp. 445–470, 2007.

L. Lei and T. Wang, “The minimum common-cycle algorithm for cyclic scheduling of two material handling hoists with time window constraints,” Management Science, vol. 37, no. 12, pp. 1629–1638, 1991.

Z. Zhou and J. Liu, “A heuristic algorithm for the two-hoist cyclic scheduling problem with overlapping hoist coverage ranges,” IIE Transactions, vol. 40, no. 8, pp. 782–794, 2008.

M.-A. Manier, C. Varnier, and P. Baptiste, “Constraint-based model for the cyclic multi-hoists scheduling problem,” Production Planning and Control, vol. 11, no. 3, pp. 244–257, 2000.

M.-A. Manier and S. Lamrous, “An evolutionary approach for the design and scheduling of electroplating facilities,” Journal of Mathematical Modelling and Algorithms, vol. 7, no. 2, pp. 197–215, 2008.

A. Che and C. Chu, “Single-track multi-hoist scheduling problem: a collision-free resolution based on a branch-and-bound approach,” International Journal of Production Research, vol. 42, no. 12, pp. 2435–2456, 2004.

J. M. Y. Leung, G. Zhang, X. Yang, R. Mak, and K. Lam, “Optimal cyclic multi-hoist scheduling: a mixed integer programming approach,” Operations Research, vol. 52, no. 6, pp. 965–976, 2004.

A. Elmi and S. Topaloglu, “Multi-degree cyclic flow shop robotic cell scheduling problem with multiple robots,” International Journal of Computer Integrated Manufacturing, pp. 1–17, 2016.

A. El Amraoui, M.-A. Manier, A. El Moudni, and M. Benrejeb, “A mixed linear program for a multi-part cyclic hoist scheduling problem,” International Journal of Science and Techniques of Automatic Control & Computer Engineering, vol. 11, pp. 612–623, 2008.

J. Feng, A. Che, and C. Chu, “Dynamic hoist scheduling problem with multi-capacity reentrant machines: a mixed integer programming approach,” Computers and Industrial Engineering, vol. 87, article 4074, pp. 611–620, 2015.

A. Che, Z. Zhou, C. Chu, and H. Chen, “Multi-degree cyclic hoist scheduling with time window constraints,” International Journal of Production Research, vol. 49, no. 19, pp. 5679–5693, 2011.

A. Elmi and S. Topaloglu, “Multi-degree cyclic flow shop robotic cell scheduling problem: ant colony optimization,” Computers and Operations Research, vol. 73, pp. 67–83, 2016.

X. Li and R. Y. K. Fung, “A Mixed integer linear programming solution for single hoist Multi-degree cyclic scheduling with reentrance,” Engineering Optimization, vol. 46, no. 5, pp. 704–723, 2014.

D. Deutsch, “Quantum theory, the Church-Turing principle and the universal quantum computer,” Proceedings of the Royal Society. London. Series A, vol. 400, no. 1818, pp. 97–117, 1985.

A. Narayanan and M. Moore, “Quantum-inspired genetic algorithms,” in Proceedings of the IEEE International Conference on Evolutionary Computation, p. 6166, Nagoya, Japan, May 1996.

K.-H. Han and J.-H. Kim, “Quantum-inspired evolutionary algorithm for a class of combinatorial optimization,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 6, pp. 580–593, 2002.

T.-C. Lu and G.-R. Yu, “An adaptive population multi-objective quantum-inspired evolutionary algorithm for multi-objective 0/1 knapsack problems,” Information Sciences, vol. 243, pp. 39–56, 2013.

A. Che, P. Wu, F. Chu, and M. Zhou, “Improved quantum-inspired evolutionary algorithm for large-size lane reservation,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 45, no. 12, pp. 1535–1548, 2015.

B.-B. Li and L. Wang, “A hybrid quantum-inspired genetic algorithm for multiobjective flow shop scheduling,” IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, vol. 37, no. 3, pp. 576–591, 2007.

T. Zheng and M. Yamashiro, “Solving flow shop scheduling problems by quantum differential evolutionary algorithm,” International Journal of Advanced Manufacturing Technology, vol. 49, no. 5, pp. 643–662, 2010.

J. Gu, X. Gu, and M. Gu, “A novel parallel quantum genetic algorithm for stochastic job shop scheduling,” Journal of Mathematical Analysis and Applications, vol. 355, no. 1, pp. 63–81, 2009.

J. Gu, M. Gu, C. Cao, and X. Gu, “A novel competitive co-evolutionary quantum genetic algorithm for stochastic job shop scheduling problem,” Computers and Operations Research, vol. 37, no. 5, pp. 927–937, 2010.

K.-H. Han and J.-H. Kim, “Quantum-inspired evolutionary algorithms with a new termination criterion, H ε gate, and two-phase scheme,” IEEE Transactions on Evolutionary Computation, vol. 8, no. 2, pp. 156–169, 2004.

S. Akpınar and G. Mirac Bayhan, “A hybrid genetic algorithm for mixed model assembly line balancing problem with parallel workstations and zoning constraints,” Engineering Applications of Artificial Intelligence, vol. 24, no. 3, pp. 449–457, 2011.

M. Srinivas and L. M. Patnaik, “Adaptive probabilities of crossover and mutation in genetic algorithms,” IEEE Transactions on Systems, Man, and Cybernetics, vol. 24, no. 4, pp. 656–667, 1994.

P. Yan, G. Wang, A. Che, and Y. Li, “Hybrid discrete differential evolution algorithm for biobjective cyclic hoist scheduling with reentrance,” Computers and Operations Research, vol. 76, pp. 155–166, 2016.

M.-A. Manier, Contribution à l'ordonnancement cyclique du système de manutention d'une ligne de galvanoplastie [PhD thesis in Automatics and Computer Science], Université de Franche-Comté, Besançon, France, 1994.

Z. Zhou, A. Che, and P. Yan, “A mixed integer programming approach for multi-cyclic robotic flowshop scheduling with time window constraints,” Applied Mathematical Modelling, vol. 36, no. 8, pp. 3621–3629, 2012.