A mathematical model and simulated annealing algorithm for solving the cyclic scheduling problem of a flexible robotic cell

Flexible robotic cells are used to produce standardized items at a high production speed. In this study, the scheduling problem of a flexible robotic cell is considered. Machines are identical and parallel. In the cell, there is an input and an output buffer, wherein the unprocessed and the finished items are kept, respectively. There is a robot performing the loading/unloading operations of the machines and transporting the items. The system repeats a cycle in its long run. It is assumed that each machine processes one part in each cycle. The cycle time depends on the order of the actions. Therefore, determining the order of the actions to minimize the cycle time is an optimization problem. A new mathematical model is presented to solve the problem, and as an alternative, a simulated annealing algorithm is developed for large-size problems. In the simulated annealing algorithm, the objective function value of a given solution is computed by solving a linear programming model which is the first case in the literature to the best of our knowledge. Several numerical examples are solved using the proposed methods, and their performances are evaluated.


Introduction
Cell manufacturing indicates a connective system among product-oriented and process-oriented systems. Using a robot in such cells helps to produce standardized items at a high production speed. 1 A cell with a number of computer numerical control (CNC) machines and a robot is called flexible robotic cell (FRC). 2 In FRCs, the CNC machines perform manufacturing processes, and the robot transports the items from the input buffer to the machines, loads/unloads the CNC machines, and transports the items to the output buffer. 3 The same group of processes is performed on all the CNC machines. Hence, each item is processed only on one machine. The considered system repeats a cycle in its long run. If the system is in a specific state at the beginning of a cycle, it reaches the same state at the end of the cycle and then repeats the same actions in the same order in the subsequent cycles. The duration of a cycle is called cycle time. Each machine processes one part in each cycle. Decreasing the cycle time in such a system means increasing the production rate. The cycle time depends on the order of the actions. Thus, determining the order of the actions to minimize the cycle time is an optimization problem.
A thorough review of inflexiable robotic cell scheduling problem with single and multiple robots including a single and dual gripper can be found in the survey by Dawande et al. 4 Gultekin et al. 5 suggested a new cycle for FRC that performs better in comparison to the classical robot move cycles for two-machine cells. Moreover, they showed that a robot-centered layout reduces the cycle time compared to an in-line layout and found an optimal number of machines to minimize the cycle time of m-machine cells. In another study, Gultekin et al. 6 presented a mathematical formulation to determine the minimum cycle time for a parallel machine cell. Jolai et al. 7 studied an FRC scheduling problem with identical part types, machines are flexible and able to swap. They determined all one-unit cycle times and proposed a new sequence of robot movements that dominates all robot move cycles. Yildiz et al. 8 proposed two pure cycles and showed that these two cycles jointly dominate all other pure cycles for a broad range of the process times. They also presented the worst case for minimizing the cycle time. Foumani and Jenab 9 developed one-unit cycles for line layout robotic cells and presented a robot move sequence that minimizes the cycle time. They also introduced the optimality regions when all parts met the first machine twice and determined the optimality conditions for different cycles when each part meets both machines twice. They carried out the sensitivity analysis for both cases and suggested the best and the worst cycles mathematically. Foumani and Jenab 10 extended their research to m-unit pure cycles when the robot is able to swap. They presented a lower bound and introduced a pure cycle that always dominates the others. Jiang et al. 11 applied two heuristics to minimize the makespan of a job scheduling problem. They considered a two-machine system where the machines are parallel and identical, and the machines are loaded/unloaded by a server. Gultekin et al. 12 studied on an FRC in which a dual-gripper robot serves the machines. They considered a twomachine FRC and found five feasible pure cycles to maximize the throughput rate. Foumani et al. 13 focused on maximizing the throughput rate of FRC problems including multi-function robotic cells, and in another study, they considered the scheduling problem of n-unit production in the FRC and found that one-unit cycles dominate the rest. Furthermore, they considered an FRC including two machines with three different scenarios of inspections including in-process and postprocess inspection, the cell with a multi-function robot, and the linear layout FRC. 14 They extended their studies on two-machine FRCs considering different pick-up scenarios. They converted a multiple sensor system to a single sensor and found the cycle times based on a geometric distribution. 15 Recently, simulated annealing algorithm (SAA) is used to solve a wide range of optimization problems. SAA is prominence from high solution performance, fine results in short times among meta-heuristics approach. In the literature, the SAA has been used for solving the traveling salesman problem (TSP), 16 the location-routing problem, 17 the emergency logistics problem, 18 the assembly line balancing problem, 19,20 disassembly scheduling problem, 21 the production and preventive maintenance problem, 22 the flow shop scheduling problem, 23 the clustering problem, 24 the facility layout problem, 25 the cell formation problem, 26 for distributed job shop problem, 27 and so on.
In the literature, some researchers employed TSP approaches for modeling scheduling problem of the FRC with m-machine and a robot. As examples, Foumani et al. 28 formulated the FRC problem including a multi-function robot, as a TSP aimed to minimize the cycle time or to maximize the production rate. Additionally, to calculate the productivity of the cell, they found the lower bound for the objective function considering the both uphill and downhill permutations. Gultekin et al. 6 developed a mathematical model for the scheduling problem of FRC considering fixed process time for the machines using TSP's Miller-Tucker-Zemlin (MTZ) method and expressed that the problem is non-deterministic polynomial-time (NP)-hard. They used CPLEX 9.0 and reported some evidence for solving the problem considering four, five and six machines in the cell. Similar to any other type of NPhard problems, computational time for solving the problem exponentially rises in case the number of the machines in the cell is increased. 29 Using twenty random inputs, they concluded that the computational time for the cell with four-machine is 7.72, for the fivemachine case is 1866.7 and with only a single run for a six-machine case, it needs 805,184.4 s to solve the problem. They did not use any meta-heuristic algorithm for solving large-sized problems. In their model, they faced with a non-linear constraint (i.e. constraint 4) and in order to propose a linear model, they defined a new auxiliary variable, namely, Y lUj and six extra liner constraints (i.e. constraints [9][10][11][12][13][14] to convert the non-linear model to a linear one. Y lUj and its related constraints have a notable adverse effect on computational time (performance) of the model especially for the FRC with more machines. These issues motivated the authors to model the problem using other exist TSP based modeling approaches along with considering variable process time for the machines. Since, the scheduling is more vital and complicated for the cells with short process time in the machines or for the cells with a busy robot (the machines have idle time, and robot activity order determine the cycle time), and also aiming to employ fewer number of variables and constraint in the model, ''Network Flow'' modeling approach is used. 30 For solving the scheduling problem of large-size cells a more feasible model from computational time points of view is needed. In addition, the use of a powerful metaheuristic algorithm from time complexity, space complexity, the advantages and disadvantages of the calculation results is essential.
In this article, TSP's ''Network Flow'' modeling approach is used to model the scheduling problem of FRC with m-machine. Moreover, SAA is examined to solve the large-size problem in the model. For doing this, first, a new mathematical model for the scheduling problem is proposed; next, for justification of the proposed model, the developed model and an existing model in the literature are solved by CPLEX under normalized conditions, and the results are compared. Finally, SAA is used to solve the large-sized problem and performances of the proposed approach using several numerical instances.
The considered problem is defined and formulated in section ''Problem definition and formulation.'' Section ''The developed SAA'' describes the proposed SAA. Section ''Experimental results'' includes experimental results about the performance of the proposed methods. The study is concluded in section ''Conclusion.''

Problem definition and formulation
A schematic of an FRC can be seen in Figure 1. Let the number of the machines in the considered FRCexplained in the introduction be m. The process time of a part on a machine is p. (In the rest of the article, the machines are assumed to be identical. Thus, they have the same process time p. In the case of non-identical machines, instead of using p for all machines, using p i for machine i as the process time of machine i will be enough to modify all the presented methods.) The robot performs all loading/unloading activities and transports all parts inside the cell. The loading activity of machine i (L i ) consists of picking, transferring, and loading a part from the input buffer to the machine i. Similarly, the unloading activity of machine i (U i ) includes getting the processed part from machine i, and transferring and putting it into the output buffer. Note that the robot stays beside of machine i at the end of L i and beside of the output buffer at the end of any unloading activity. A cycle time is a duration spanning from the starting of the system from a specific state and returning to the same state. In order to start such a cyclic production, the system needs a setup. Each machine may be loaded or emptied at the beginning of the cycle. During a cycle, each machine must be loaded and unloaded once. Let L is the set of loading activities, U is the set of unloading activities, and A is the set of all loading and unloading activities.
Let e be the loading/unloading time for each machine and each buffer and d be the robot travel time between the input buffer and the first machine, between two consecutive machines, and between the last machine and the output buffer. d ab in the following formula gives the certain time needed between the completion times of activities a and b for the robot's operations such as taking part, moving, and putting the part when activity b follows activity a, that is, d ab is not related with the process times on the machines When the process times on the machines are considered, some uncertain amount of waiting time for the robot can be needed. Let's consider machine i and its loading activity L i and unloading activity U i . At the completion time of L i , the machine starts its operation and finishes it after p time unit. Then, the robot may start unloading this part. During the unloading operation, the robot takes part from machine i (it takes e time units), moves to the output buffer (it takes ((m + 1 À i)d) time units) and puts the part into the output buffer (it takes e time units). Thus, the time between the completion of L i and U i must be at least (2e + (m + 1 À i)d + p). There may be several other activities between L i and U i , and the total time for performing those activities may not be big enough to complete the process on machine i. In such a case the robot must wait for the end of the process on machine i. The waiting times depend on the order of the activities. Note that d ab does not contain this uncertain amount of waiting time.
Since the robot performs the same order of activities in a cycle, to prevent permutation and have a fixed cycle, we consider L 1 as the first activity. Thus, the time from L 1 to the next L 1 is the cycle time (T), and the problem is to determine the order of all loading and unloading activities to minimize T. Decision variables x ab = 1 if activity b is performed after activity a by the robot 0 otherwise & t ab : the completion time of activity b when it is performed just after activity a, it is zero if activity b is not performed just after activity a; w ab : the time that the robot waits before starting activity b when it is performed just after activity a, it is zero if activity b is not performed just after activity a The objective is to minimize the cycle time, which is the time that the robot performs L 1 after the last activity. Note that the cycle starts at the time that L 1 is completed. That time is considered as time zero. During the cycle, all the activities, including L 1 , must be completed. So, the end of a cycle is the completion time of L 1 , which is also the beginning of the next cycle. By constraints (2) and (3), it is guaranteed that the robot performs all the activities. It passes from one activity to another activity. Constraints (4) and (5) fix t ab and w ab variables to zero if the robot does not perform activity b just after activity a. If activity b is performed just after activity a, then x ab is 1 and the corresponding t ab and w ab variables are allowed to be positive by constraints (4) and (5). Note that the waiting time for unloading a part cannot be more than the process time. Because of this in constraint (5), p is used as the coefficient of x ab instead of a big number M. Constraint (6) is the balance constraint. The completion time of an activity equals to the completion time of the previous activity plus the robot operation times between these two successive activities and plus the waiting time before performing the later activity. Constraint (7) is the balance constraint for L 1 . Constraints (8), (9), and (10), together, guaranteed to have enough time between a load of a part and unload of it for finishing its process. Constraint (11) does the same thing for the part processed on the first machine.

The developed SAA
In the attempts of solving the problem using the mathematical models, it is seen that the solution time increases very rapidly when the number of the machines increases, and mathematical model based exact solution methods fails to solve the problems. The SAA is a wellknown and efficient meta-heuristic approach. It runs using a single solution at a time, hence does not cause memory shortage problems for even very large-size problem instances. Moreover, it produces a feasible neighboring solution and does not need a repair algorithm, which may lead to highly diversified solutions and deteriorates intensification. Therefore, especially in order to solve large-size problems, an SAA is developed. The method starts with an initial solution which is constructed by any constructive algorithm. At any iteration, the algorithm generates a neighboring candidate solution by making a randomly chosen small change on the current solution. If the candidate solution is better than the current one, the candidate solution is adopted as the new current solution. However, if the candidate solution turns out to be worse than the current solution, the algorithm may either adopt the candidate solution as the next current solution with some acceptance probability or reject it. By giving a chance to move to inferior solutions, the algorithm obtains some capability for escaping from the local minimums. The function that gives an acceptance probability of a bad solution is where F is the evaluation function, and T is the control parameter of the algorithm called temperature. The probability of accepting an inferior solution decreases if the difference between the current solution and an inferior candidate solution increases, or the temperature drops. At the beginning of the algorithm, the value of temperature is higher, and it falls during the search according to a function known as the cooling schedule that provides intensification over time. Because of the higher value of temperature, initially the algorithm searches the space roughly; however, because of the cooling effect over time, it focuses on some good solution regions. The algorithm stops when a termination criterion is satisfied. Either the number of iterations or the running time or the final value of the control parameter T can be used as the termination criterion. The details of the developed algorithm are given in the following sections.

Representation
A solution is presented by an array having 2m elements in which the numbers of 1 to m correspond to the loading of the first machine to the mth machine, and the numbers of m + 1 to 2m correspond to the unloading of the first machine to the mth machine, respectively. To prevent permutations, the first element of each array is always 1. For example, L 1 L 3 L 4 U 2 U 3 U 1 U 4 L 2 order for a four-machine is presented in Figure 2.

Initial solution
After some preliminary experiments and testing different policies, it is decided to generate the initial solutions, randomly. Number 1 is fixed to the first order and then each of the remaining numbers up to 2m assigned to an empty order, randomly. Generating random solutions helps to start from different initial solutions in each run that avoid entrapment in local optima.

Computing cycle time for a given solution
All of the parameters, except waiting times, are easy to compute for finding the objective value of a given solution (an order of the activities). Let's consider a fourmachine case where e = 1, d = 2, and p = 80, and the L 1 L 3 L 4 U 2 U 3 U 1 L 2 U 4 order for this case. The times needed for performing the robot operations, such as loading, unloading, and transportation, between the successive activities in this order can easily be computed by the d ab formula. For example, d L1L3 is 2e + 4d = 10. Similarly, times between L 3 L 4 , L 4 U 2 , U 2 U 3 , U 3 U 1 , U 1 L 2 , L 2 U 4 , and U 4 L 1 are calculated as 16, 12, 10, 18, 16, 8, and 14, respectively.
However, as it is explained in section ''The developed SAA,'' if the robot turns to a machine for unloading the part earlier than its completion, it should wait. The robot may wait before each of the unloading activities, and these waiting times are not so trivial to compute. Let w i be the waiting time before performing U i . Note that these waiting times may be zero. Then, the updated times between the completions of the successive activities in the given order are 10 for L 1 L 3 , 16 for L 3 L 4 , (12 + w 2 ) for L 4 U 2 , (10 + w 3 ) for U 2 U 3 , (18 + w 1 ) for U 3 U 1 , 16 for U 1 L 2 , (8 + w 4 ) for L 2 U 4 and finally 14 for U 4 L 1 . Thus, the cycle time for this order is CT = 10 + 16 + 12 + 10 + 18 + 16 + 8 + 14 + w 1 + w 2 + w 3 + w 4 = 104 + w 1 + w 2 + w 3 + w 4 .
As it is explained in section ''The developed SAA,'' the time between the completions of L i and U i must be at least (2e + (m + 1 À i)d + p) for machine i. When we consider the above-given order, the robot performs the first L 1 and then L 3 , L 4 , U 2 , U 3 and then U 1 . Hence, the total time between the completions of L 1 and U 1 is 10 + 16 + (12 + w 2 ) + (10 + w 3 ) + (18 + w 1 ) = 66 + w 1 + w 2 + w 3 . This time should be at least (2e + (m + 1 À 1)d + p) = 90. Therefore, on w 1 , w 2 , and w 3 , we have the following condition When we consider machine 2, U 2 is earlier than L 2 in the given order. So, a part is loaded to machine 2 in a cycle, and it is unloaded in the following cycle. It is shown in Figure 3.
The time between L 2 and U 2 is then (8 + w 4 ) + 14 + 10 + 16 + (12 + w 2 ) = 60 + w 2 + w 4 . This time should be at least (2e + (m + 1 À 2)d + p) = 88. Then, we have the condition Considering machines 3 and 4, the following conditions are obtained, respectively 64 + w 1 + w 2 + w 3 + w 4 ! (2e + (m + 1 À 4)d + p) Consequently, the cycle time for the given order is CT = 104 + w 1 + w 2 + w 3 + w 4 , and in order to find the minimum CT, the values of w 1 , w 2 , w 3 , and w 4 must be determined considering the above conditions on them. This can be done by solving the following linear programming (LP) model s:t: The objective function values of the created orders in the developed SAA are computed by solving such LP models.

Generating candidate solutions
Shift, swap, and reverse operators are used to generate neighboring solutions of the current solution. By the shift operator, the order of a randomly selected activity is changed randomly when the order of the other activities remained the same. Using the swap operator, two activities are chosen randomly, and only their orders are replaced with each other. In the reverse operator similar to the swap operator, two activities are randomly selected. Then, these two activities and all of the activities between these two activities are reversed. Figure 4 gives examples for each of these mechanisms.
In each iteration of the developed SAA, three new solutions are generated from the current solution using each of the shift, swap, and reverse operators. The objective function of each is computed, and the best of these three new solutions is selected as the generated candidate neighboring solution. This solution is adopted as the next current solution if it is better than the current solution or it is accepted using the acceptance probability function.

Cooling
The geometric cooling, which is the most common cooling method, is applied. According to this method, the developed SAA starts its search at an initial temperature. Then, after each iteration, a certain percent of T is counted as the value of T for the next iteration, that is, T = a*T where 0 \ a \ 1.

Stopping criterion
A limit on the solution time is used as the stopping criterion. The starting time is kept, and the total time from the starting time to the current time is computed after each iteration. When it exceeds the limit, the search is stopped. Figure 3. The time from L 2 to U 2 for the order

The pseudocode of the developed SAA
The pseudocode of the developed SAA is given below.

Experimental results
The performances of the proposed methods are evaluated on several problem instances. In these experiments, an Intel Ò Coreä i5-3320 CPU at 2.60 GHz with a RAM of 4.0 GB computer is used for the runs.

Justification of the proposed mathematical model by an exist reference model
The proposed mathematical model is compared with the model presented in Jiang et al. 11 Both models have been coded in CPLEX 12.6 software. Table 1 shows the related results, and Figures 5-7 display the solution times of the models. In the test instances, e and d are 1 and 2 time units, respectively.
According to the results in Table 1 and Figures 5-7, the model proposed in this study solves the problem in a shorter time than the model reported in reference. The solution time of the reference model increases rapidly when the number of the machines increases. The proposed model solves the problem in a concise time when the process time is small. However, much longer times are needed to solve these cases using the other (c) Generate a neighboring solution of the current solution using the reverse mechanism.
(d) Compute the cycle time of these neighboring solutions and select the best of them as the candidate solution.
Step Record the current time: Tnow. If the duration from Tstart to Tnow is more than Time Limit, then STOP and present the Best Solution, otherwise T = a*T and go to Step 3. model. Considering higher process times, solution times of the models converge to each other. In these cases, the solution time of the proposed model increases rapidly when the number of the machines increases too. It should be noted that any of the models could not solve problem instances having more than six machines in a cell with high process times.
In Table 1, when the process time is large, the optimal cycle time is larger than the process time by 24; it is larger by 28 in five-machine cell, and by 32 in sixmachine cell. By analyzing the robot task sequences derived in experiments, it seems that an optimal cycle time can be derived from a large number of machines. In order to validate or deny this claim 10 instances of FRC problem with four and five machines are solved.
In these examples, the process times of the machines are different and randomly selected in [1,100] and [1,300] ranges, where the range of the e and d are randomly selected from the [1,4] and [1,5] intervals, respectively. As shown in Table 2, the resulted cycle time indicates that the cyclic times are not harmonic, and it is not predictable from the results of the cell with the fewer machine for a small and large range of process time.

Performance of the developed SAA
The performance of the proposed SAA depends on the values of its parameters, which are the initial value of T, Time Limit as the stopping criterion, and a. To determine these values, some experiments are needed. Since Taguchi method proposes fractional factorial experiments, it is very effective in parameter setting. 31 Based on noise minimization, this method selects the best level of the parameters. Using the following equation, the deviation of the response is examined, wherein Y designates the value of reply, and n characterizes the number of orthogonal ranges 32 For each of initial value of T, Time Limit, and a, ranges are determined. These ranges are  for the initial value of T, [1,5] for Time Limit, and [0.993,0.997] for a. For each of these parameters, three different values are used: (1) the lower bound, (2) the average, and (3) the upper bound of the corresponding range. Thus, nine different combinations of these values are tested. In these tests, the developed SAA is used for solving five different instances which are (4, 75), (6,150), (8,250), (10,500), and (12,750), where the first entry shows the number of the machines and the second one shows the process time (i. e. (m, p)). For each of the nine combinations of the parameter values, each of these five test instances is solved by the proposed SAA 10 times. The average of the cycle times of the   best solutions found by these 10 runs recorded and presented in Table 3.
Then, the S/N ratios are computed using the results in the last column of Table 3. Figure 8 shows the results of S/N ratios. According to these ratios, when the initial value of T is set to its lower bound (which is 90), Time Limit is set to its upper bound (which is 5 min), and a is set to its upper bound (which is 0.997), the best results are obtained. Thus, these values are used in the following tests.
First, the performance of the proposed SAA is tested on the above test instances whose optimal solutions are found by the mathematical models. Table 4 contains the cycle times and solution times of all the examined test problems for four-to six-machine cells found by the proposed SAA.
According to the results in Table 4 and Figures 9-11, the proposed SAA found almost all optimal solutions. Only 6 of the 33 instances could not be solved optimally. In these instances, the gap between the optimal cycle times and the best cycle times found by the proposed SAA is less than 10% of the optimal cycle times. Thus, it may be concluded that the proposed SAA has a very good performance and it may be used to find good solutions to larger instances.

Conclusion
Using TSP's network flow modeling approach, a novel mathematical model is proposed for solving the cyclic scheduling problem of FRCs with identical and parallel      machines which are located on a line. Since the computation time for solving large-sized problems in the model was high, an SAA is examined. The numerical experiments show that the proposed mathematical model performs better in comparison with the one existing in the literature. Solution times by the mathematical models increases very rapidly when the size of the problem increases. The proposed SAA finds almost the optimal solutions of the considered problem instances in affordable short times. This approach might be seen in future studies as a real-time scheduling approach in flexible manufacturing systems. Considering these problems under uncertainty of each parameter can be regarded as another future subject. Finally, it is interesting to consider circular layouts. The cases where machines have buffers with limited capacity may also be considered for future research.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.