Multi-objective AGV scheduling in an FMS using a hybrid of genetic algorithm and particle swarm optimization

Flexible manufacturing system (FMS) enhances the firm’s flexibility and responsiveness to the ever-changing customer demand by providing a fast product diversification capability. Performance of an FMS is highly dependent upon the accuracy of scheduling policy for the components of the system, such as automated guided vehicles (AGVs). An AGV as a mobile robot provides remarkable industrial capabilities for material and goods transportation within a manufacturing facility or a warehouse. Allocating AGVs to tasks, while considering the cost and time of operations, defines the AGV scheduling process. Multi-objective scheduling of AGVs, unlike single objective practices, is a complex and combinatorial process. In the main draw of the research, a mathematical model was developed and integrated with evolutionary algorithms (genetic algorithm (GA), particle swarm optimization (PSO), and hybrid GA-PSO) to optimize the task scheduling of AGVs with the objectives of minimizing makespan and number of AGVs while considering the AGVs’ battery charge. Assessment of the numerical examples’ scheduling before and after the optimization proved the applicability of all the three algorithms in decreasing the makespan and AGV numbers. The hybrid GA-PSO produced the optimum result and outperformed the other two algorithms, in which the mean of AGVs operation efficiency was found to be 69.4, 74, and 79.8 percent in PSO, GA, and hybrid GA-PSO, respectively. Evaluation and validation of the model was performed by simulation via Flexsim software.


General
In today's competitive market, customer satisfaction plays a key role in companies' market share. Therefore, organizations have shifted their strategy from manufacturing large quantities of a single product to a range of products, and improving the quality and on-time delivery. To meet these challenges, organizations should have a flexible manufacturing platform. In automated manufacturing environments, FMS provides wide flexibility and concurrent production of a wide variety of parts in small capacities. It comprises material handling a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Model derivation
This section explains the mathematical model development for AGV scheduling using the three criteria selected based on the literature review. The three criteria are categorized into two main objectives: (1) minimizing the makespan and (2) minimizing the number of AGVs while considering the AGV's battery charge. Each sub-section explains the mathematical definitions used to develop the model. However, prior to the model development it is necessary to define the conditions and limitations considered in the model framework. Thus, the following conditions were defined for the model development: • All AGVs have unit-load capacity.
• AGVs and machines operate continuously without breakdown.
• There are no traffic problems, collision, deadlock or conflict.
• AGV loading and unloading times are fixed and considered in travel times.
• AGVs can always park at their pick-up/drop-off (P/D) locations.
• The velocity of the vehicles is constant and vehicles move forward only.
• Loading/unloading (L/U) equipment such as pallets are sufficiently allocated as well as input and output buffer for machines.
• The machine-to-machine distance and L/U point-to-machine distances are known.
• Each machine operates only one product at a time.
• The setup times are included in the time of production.
• The AGVs are stored in the home until dispatching commands are allocated.

Notations n
Total number of jobs j, j 0 Indexes of jobs, genes', and dimensions' code, j, j 0 = 1, 2, . . ., n m j Total number of operations for each job j i, i 0 Indexes of operations, i, i 0 = 1, 2, . . ., m j, j 0 θ Total number of operations for all of jobs The velocity of α th particle on d th dimension at iteration (t) v tþ1 ad The velocity of α th particle on d th dimension at iteration (t+1) V t a The velocity of α th particle in the swarm at iteration (t) q t ad The position of α th particle on d th dimension at iteration (t) q tþ1 ad The position of α th particle on d th dimension at iteration (t+1) Q t a The position of α th particle in the swarm at iteration (t) B t ad The best position of α th particle on d th dimension found so far G t d The global best position of the swarm on d th dimension found so far φ 1 and φ 2 Uniformly distributed random numbers in the interval [0, 1] C 1 Self-confidence C 2 Swarm The β th weight of the β th objective function ψ A ratio to make balance among objectives with different ranges of value f(x) Fitness function Minimizing the makespan. This step involves calculating makespan (MS) which is the time required for all operations to be completed. Makespan can be expressed by Subject to: Constraint number 3 ensures the feasibility of completion time of the first operation of each job. Constraints number 4 and 5 ensure the feasibility of pick up time of operations. Inequality number 6 describes the operations precedency constraint. Inequalities number 7 and 8 represent the operation and the AGV un-overlapping constraints respectively.
Minimizing the number of AGVs. This step involves calculating the number of AGVs, which is denoted by NA by considering the AGVs battery charge sufficiency. Number of AGV can be expressed by Subject to where Eq 10 makes sure the assigned AGV has enough battery charge to do the job and return home, while it chooses the AGV which takes less time to reach the point. Eq 11 calculates the charge that AGV needs to do the job and return home. As battery-run-time of an AGV and battery-charging-time can be defined depending on the type of batteries used, charging methods, charge rate, application, manufacturer, and assignments the vehicles perform, γ has been defined to adopt to any kind of battery, charging method, etc. The automatic and opportunity battery charging considered here, which on average, an AGV charges for 10-12 minutes every hour [30,31]. Eq 12 determines the suitable time for adding a new AGV. Multi-objective evaluation. Decision makers refer to choosing a solution out of all the efficient solutions as a posteriori approach. Pareto is one well-known pioneer in multi-objective optimization problems. In this method, Pareto-optimal set is a group of best trade-off schedules, and Pareto-front refers to a set of Pareto solutions [32]. Overall fitness function formulation for two objectives is described by Where δ 1 is the weight of first objective function and ψ is a ratio to make balance among objectives with different ranges of value [33][34][35], which is defined by AGV specifications This step involves calculating specifications of AGV number a including its total running time denoted by RtA a (loaded (LtA a ) + unloaded time (UtA a )), its waiting time (WtA a ), its idle time (ItA a ), its consumed battery charge (ChA a ), its battery usage percentage (BU), and its efficiency (AE) by Eqs 15 to 27.
WT a ji ¼ ChA a ¼ ChT

Proposed algorithms
Three different evolutionary algorithms (EA) have been developed to optimize the mathematical AGV scheduling model. The three algorithms (GA, PSO, and hybrid GA-PSO) are later evaluated in terms of their strength and suitability for the scheduling problem.

Genetic algorithm
GA is a search algorithm based on the mechanics of the natural selection process. The major steps of GA algorithm development according to the study objective are described in this section. However, readers for a thorough review of the GA are referred to publications of [36][37][38][39][40].
Step 1. Initializing parameters. It involves setting the parameters of the GA and creating the first generation of chromosomes based on the notations section. The general schematic for reading data for the problem is presented in Table 1. The first column shows a chromosome (C r ) and the second one shows the genes (G e ) of the chromosome. The encoding of each gene is presented in the third column, which will be discussed later.
Step 2. Initializing population. A set of chromosomes is needed to create a population. For constructing a chromosome, it is necessary to define a proper genetic representation (encoding) due to its significant effects on all the subsequent steps of the GA.
Chromosome representation and encoding. As it is shown in Table 1, each chromosome is formed by genes. The order of genes represents the priority of operations, which decreases from left to right; and the genes' code defines operations related to each job. Gene's code are the same as their job number so that all the genes related to J 1 operations have the code '1' and subsequently the code '2' is given to all the genes related to the operations of J 2 , and so on. As the operations of each job are expected to be performed sequentially, the repetition of genes' code represents the corresponding operation number of the job as clearly described in the following example.
The number of genes in each chromosome equals the number of total operations in a jobset, which is expressed by: Chromosome generating. A chromosome (C r ) is a random construct of operations, which is expressed by where j, j' are indexes of jobs, j, j' = 1, 2, . . ., n and m j , j' = number of operations for each job. O ji is the operation i of job j. Chromosome coding and generating is explained below via an example of 3 jobs (J 1 , J 2 , and J 3 ). Each job has 4, 3, and 5 operations respectively. Overall, there is θ = 4+3+5 = 12 operations. Therefore, the chromosome is a random construct of ½1111 |ffl{zffl} . A sample could be [221132313133]. Here, code '1', '2', and '3' imply operations of J 1 , J 2 , and J 3 respectively. From the left, the first '2' represents the first operation of J 2 , the second '2' represents the second operation of J 2 , the first '1' represents the first operation of J 1 , and so on.
Step 3. Multi-objective evaluation. After initializing the population size, each chromosome is evaluated by minimizing the makespan and the AGV numbers, while considering the battery charge of AGV that are defined by Eqs 1 to 12. Then, the total fitness values of the efficient frontiers will be calculated based on Eq 13.
Step 4. New population. New population will be produced based on the below sub-steps: selection, crossover, elitism, and mutation operation.
Selection. To constantly enhance the population overall fitness, selection helps to discard the bad/weak designs and only keep the best ones in the population. It increases the likelihood of selection of fitter individuals for the next generation. There are a few different selection methods but their basis is the same. The tournament candidate selection, which is a proportionate random selection method, is used in this study. In this method, every individual in the population is paired at random with another. The fitness values of each pair is compared. The fitter individual of the pair moves on to the next round, while the other is disqualified. This continues until there are a number of winners equal to the desired number of parents. Then, this last group of winners is paired as the parents for new individuals [41].
Crossover. Crossover operator generates two new chromosomes for the next generation out of two selected chromosomes by exchanging some of their genes. This study employs two crossover operators based on partial strings exchange; a one-point crossover and a two-point crossover [42], where the one-point crossover is illustrated in Fig 1 based on the example in step 2.
The offspring of crossover between the strings may not produce a legal encoding, for example, uncorrected number of operations per job may be seen. Therefore, they should be repaired and legalized. For repairing mechanism, counting from the left, the redundant genes will be deleted and compensate the missing ones, in order for each offspring to comprise all the operations of all the jobs. Repair mechanism is shown in Fig 2. Legal chromosomes for the example in step 2 should include four code '1', three code '2', and five code '3'. In Fig 2a, counting from the left, in offspring 1, code '2' is repeated four times, but there are three operations for J 2 , so it should repeat three times. There is one code '2' that is redundant and should be replaced by the missing code. Code '1' is repeated four times, which is correct, but code '3' is repeated only four times, which should be 5 times. So in Fig 2b, the 4 th code '2' will be replaced by number '3'. In offspring 2, number '2' is repeated two times and number '3' is repeated six times, so the last code '3' will be changed to code '2'.
The number of crossovers is calculated based on the crossover rate (CR) and population size (PS) using Mutation. Mutation is another important operator of GA that initiates extra variability in a population to create and maintain the diversity. Mutation is not applied on chromosomes that are immune. The number of mutations in each generation is calculated using Eq (31) based on the mutation rate (Pm), population size (PS), and maximum gene code (G max ).

Number of mutations ffi PS
Shift mutation is used in this study [43] and it is shown in Fig 3. Based on the coding used in this study, chromosomes produced out of shift mutation are legal and no need to be repaired.
Elitism. The first three best chromosomes from each generation are transferred directly to the next generation in the elitism step to avoid annihilation. It is possible to maintain a fixed fitness value in some generations, but elitism makes sure they will never deteriorate.
Step 5. Termination. The loop of chromosome generation is terminated when the number of generation reaches its maximum, then the elite chromosome returns as the best solution.

Particle swarm optimization
PSO is a population based stochastic technique inspired by social behaviour of bird flocking or fish schooling. Extensive reviews on PSO algorithm development can be found in [44][45][46]. PSO with limited information has been studied to avoid the waste of information [47]. Whereby population topologies and their performance had be studied [48]. Gao & collaborators proposed a method called Selectively-informed PSO (SIPSO) to allow the particles to learn at difference strategies based on their connections [49]. The PSO configuration for the mathematical model is described in details in the following steps: Step 1. Initializing parameters. Initialization involves setting the parameters of the PSO and creating a group of particles to make the initial swarm. The general scheme for reading the  data in the problem is presented in Table 1. The third column shows a particle (PR α ) and the forth one shows dimensions of the particle (d). The dimensions' codes are presented in the third column, which will be discussed later.
Step 2. Initializing population (swarm). A group of particles are needed to create a swarm. Each particle has position (Q) and velocity (V) in the search space at iteration (t), where they are described briefly in the following sub-steps: Particle position. First position of particle is filled by two digit numbers for 'd' dimensions of the particle using Eq 32. The number of dimensions is equal to the total number of operations, which is calculated by Eq 28.
where q min = 0, q max = 10 and φ 1 is a uniform random number between 0 and 1. Particle velocity. Initial velocities for the PSO particles are generated by the formula below: where v min = 0, v max = 10 and φ 2 is a uniform random number between 0 and 1.
Step 3. Particle representation and encoding. Every possible sequence of operations is considered as a particle, where the dimension of the particle represents each operation. Three sub-steps for encoding a particle are as follows: applying smallest position value (SPV) rule, assigning the dimensions' codes to the particles, and identifying sequence of operations in each job.

Applying smallest position value (SPV) rule. SPV is a rule that facilitates transformation
of the continuous PSO algorithm to discrete cases applicable to all types of the scheduling problem [50]. As an example for better understanding of SPV rule, the corresponding sequence of a given continuous position like [0.3, 1.2, 0.9, 2.4] would be [4,2,3,1]. In a descending order, '0.3' is the smallest value and its sequence will be '4'; '2.4' is the largest so its order in the group will be '1'.
Assigning the dimensions' codes to the particles. In this stage, the dimension's codes as it is shown in the 6 th column of Table 1 are assigned to the particles. Dimension's codes are based on the job number.
Identifying sequence of operations in each job. From the left side, the first appearance of a job number is assumed the first operation of that job (i.e., O j1 ). Similarly, the second time repetition of the same job number stands for the second operation of the same job (i.e., O j2 ) and so on. Once the first encountered generated number is assigned to the first operation of a job, this technique automatically handles the precedence constraints.
The stages of encoding an example with 3 jobs are shown in Table 2. Each job has 4, 3, and 5 operations respectively. The total operations are 12, which means the particle sample will have 12 dimensions being randomly generated using Eq 32 and shown in the first row of Table 2. In the second row of the Table, based on SPV rule, the numbers of 1 to 12 are assigned to the particles in an ascending order. In the third row, the dimensions' codes based on the job numbers are given to the particles as follows: first four numbers are assigned to the first job, so their code is '1', followed by the second three numbers assigned to the second job, so their code is '2' and the remaining five numbers are assigned to the third job and their code is '3'. The sequence of operations in each job is shown in the fourth row of the Table 3. From the left, the first particle has the code '1', so it belongs to job 1 and it is the first code '1', which makes it the first operation of job 1 denoted by O 11 ; the next code is '2', so it belongs to job 2, but as it is the first code '2', it is the first operation of job 2. The same structure is followed for the remained 10 operations.
Step 4. Multi-objective evaluation. Once the swarm is generated, each particle is evaluated by minimizing the makespan and AGV numbers, while considering the AGV battery charge, which are defined by Eqs 1 to 12. Then, the total fitness values of the efficient frontiers will be calculated based on Eq 13.
Personal best. B t a represents the best position associated with the best permutation and fitness value of the particle α obtained so far and is called the personal best. For each particle, B t a can be determined and updated at each iteration.
Global best. G t denotes the best position of the globally best particle achieved so far in the whole swarm.
Step 5. New swarm. To produce a new swarm, the position and velocity of the particles should be updated. Updated particles will be evaluated again according to the step four and the best local and global particle will be determined. This procedure will be repeated up to a point where the termination criterion is satisfied. The updating procedure is explained as follows: 1. Updating the velocity of each particle. The velocity of each particle is updated using where C 1 is self-confidence while C 2 is swarm confidence (common values of C 1 and C 2 varies between 0.1 and 0.5 but the values between 0.1 and 1 has been tested as well. In some literature, the value of 2 have also been observed). Inertia weight (ω) is a parameter to control the impact of the previous velocity on the current velocity [51,52]. Let ω be varying with time by the following linear decreasing function. 2. Updating the position of each particle. The position of particle is updated using the updated velocity as below: Step 6. Termination. The loop of swarm groups is terminated when it reaches the maximum number of iteration, then the particle with global best returns as the best solution.

Hybrid GA and PSO
The PSO algorithm is a more robust optimization algorithms compared to many other algorithms as it can work almost independent from the problem. It does not require extensive prior-knowledge regarding the problem except the fitness evaluation of each particle [53]. On the other hand, GA has the capability of simultaneous evaluation of many points in the search area, which increases the probability of finding the global solution of the problem. Hybridization of EAs has been studied in many researches [54][55][56]. Generally, PSO functions based on the social interaction knowledge and all the individual particles will be considered in each generation. Unlike PSO, fitter chromosomes will be chosen in GA and the weaker ones will fade away from generation to generation [57]. Hence, by integrating the advantages of the compensatory properties of PSO and GA, their hybrid is used to obtain better result [58][59][60]. In the proposed GA-PSO algorithm for this study, after generating and evaluating the initial swarm and after position and velocity updating, the crossover operation has been used in the GA segment to avoid premature convergence; and a mutation operation was applied to maintain the diversity of the swarms. Elitism step was also performed to improve immune particle filter. Fig  4 illustrates the steps of hybrid GA-PSO and some parts of programming codes are listed in S1 Appendix.

Computational results and discussion
To validate the model, two numerical examples have been used. The first example had 6 jobs (J 1 , . . ., J 6 ) processing on 6 machines (M 1 , . . ., M 6 ), and each job with 2 to 5 operations. The second one with 15 jobs (J 1 , . . ., J 15 ) processing on 10 machines (M 1 , . . ., M 10 ), and each job with 1 to 5 operations [23,61]. Tables 3 and 4  The makespan of scheduling before optimization by a random sequence and assigning one AGV to each of the six jobs for example 1 is shown in Fig 5. However, illustration of before optimization for example 2 was not possible due to its big figure size and detail.
Based on the experimental approach, the best setting of hybrid GA-PSO parameters was found to be the crossover and mutation rates of 0.2 and 0.08 respectively, C 1 = 0.01, C 2 = 0.9, ω min = 0.01, and ω max = 0.5. The algorithms were run 30 times, each run with a population size of 100 in 100 iterations, and their first two best results based on different AGV numbers are shown in Table 7.
The third column of Table 7 shows the best result of each algorithm, and the forth column shows the best result of each algorithm using a different number of AGVs compared with the first column. The fitness value in Table 7 has been calculated based on Eq 13, c ¼ max ðMSÞ max ðNAÞ , and d 1 ¼ 2 3 . Max (MS) was presumed to be equal to the sum of travel times and operation times. Max (NA) was presumed equal to the whole number of operations. All the steps were repeated To investigate the effect of optimization methods on AGV scheduling, AGVs' specification both before and after the optimization for example 1 were explored. The studied specifications are AGVs' total running time (loaded and unloaded), idle time, battery usage, and operation efficiency computed using Eqs 15 to 27. In Fig 9, prior to the optimization, the AGVs total running time is low because a higher number of AGVs are employed with no intention to use their highest potential, compared to the optimized schedule, thus the idle time of AGVs has increased dramatically. The AGV number four (AGV 4 ) had the highest operation efficiency (37.3%) before the optimization; although the scheduling model was designed to sequentially appoint tasks to AGVs based on their numbers' order, so that AGV number one (AGV 1 ) would have the highest operation efficiency level.
In GA-PSO, the makespan, number of AGVs and their idle time have been reduced, and consequently efficiency of AGVs' operation has enhanced (Fig 9). Potency of hybrid GA-PSO in solving scheduling problems and its superiority against its constituting algorithms have also  been largely mentioned in other published studies [62][63][64][65]. Overall, application of the hybrid GA-PSO in scheduling studies is concluded to be more effective than its constituting EAs.

Simulation by Flexsim
In order to prove the feasibility of the proposed model, a simulation practice based on the above example has been performed using the Flexsim software. Fig 10 shows a scene from the simulation space. The simulation outcome confirmed the optimization results by obtaining equal makespan magnitude to all the three algorithms. The experimental results proved the validity and feasibility of the model, which provide useful reference for further research on the

Conclusion
This research focused on the multi-objective AGV scheduling in an FMS using GA, PSO, and hybrid GA-PSO algorithms. A model was developed for the task scheduling of AGVs considering multiple objectives of minimizing makespan and number of AGVs, while considering the battery charge of AGVs. Using the numerical examples, near-optimum schedules for the combined objective functions were obtained. The inter-comparison of the three algorithms results showed that the hybrid GA-PSO yields the least makespan and AGV numbers. Literature has also largely exhibited the excellence of hybrid GA-PSO over its constituents in solving the scheduling problems [62][63][64][65], however the scheduling model proposed in this study distinguishes it from previous studies. The scheduling problem was further scrutinized by comparing the AGVs characteristics such as the total running time (loaded and unloaded), idle time, battery usage, and operation efficiency-before and after the optimization. It was found that after the optimization, despite a small rise in AGVs' total running time (loaded and unloaded), the AGVs' idle time reduction enhanced the AGVs' operation efficiency. This leads to effective utilization of AGVs and hence the overall efficiency of the system will be enhanced. In line with the experimental results, the AGV system simulation using the Flexsim software has also proved the feasibility of the developed model and suitability of the optimization algorithms for  the scheduling problem. The developed model can be adopted to any FMS with different configuration and environment, and it can be applied for optimizing the objectives separately or in a combinatorial fashion. This research framework can be stretched out to examine newer algorithms and hybrids, and also employ more criteria in the scheduling model for further studies in this context.
Supporting information S1 Appendix. Programming codes for hybrid GA-PSO. (PDF)