An Improved Version of Discrete Particle Swarm Optimization for Flexible Job Shop Scheduling Problem with Fuzzy Processing Time

The fuzzy processing time occasionally exists in job shop scheduling problem of flexible manufacturing system. To deal with fuzzy processing time, fuzzy flexible job shop model was established in several papers and has attracted numerous researchers’ attention recently. In our research, an improved version of discrete particle swarm optimization (IDPSO) is designed to solve flexible job shop scheduling problem with fuzzy processing time (FJSPF). In IDPSO, heuristic initial methods based on triangular fuzzy number are developed, and a combination of six initial methods is applied to initialize machine assignment and random method is used to initialize operation sequence. Then, some simple and effective discrete operators are employed to update particle’s position and generate new particles. In order to guide the particles effectively, we extend global best position to a set with several global best positions. Finally, experiments are designed to investigate the impact of four parameters in IDPSO by Taguchi method, and IDPSO is tested on five instances and compared with some state-of-the-art algorithms. The experimental results show that the proposed algorithm can obtain better solutions for FJSPF and is more competitive than the compared algorithms.


Introduction
Flexible job shop scheduling problem (FJSP) is a general form of job shop scheduling problem (JSP), which is completely NP-hard problem [1].Up to date, many different mathematical methods and techniques, such as mixed integer programming [2], disjunctive graph [3], priority dispatch rules [4,5], and neural networks [6], are developed to optimize FJSP.Nature-inspired algorithms have made great progress in solving optimization problems in the past two decades [7,8] and are extended to solve a wide range of applications such as scheduling optimization problems [9][10][11][12], transportation problems [13,14], traveling salesman problems [15][16][17], and engineering optimization design [18].Due to their validity and robustness, various nature-inspired algorithms (particle swarm optimization (PSO) [19], differential evolution (DE) [20], firefly algorithm (FA) [21], artificial bee colony (ABC) [22,23], harmony search (HS) [24,25], evolutionary algorithm (EA) [26], and biogeography-based optimization (BBO) [27]) also attract a large number of researches' attention on the application of FJSP.By defining the fitness of particles using Pareto ranking and crowding distance, Shao et al. presented hybrid discrete algorithms for addressing FJSP, which had a global search procedure of particle swarm optimization and a local search procedure of simulated annealing [19].Utilizing differential evolution and local search, Yuan and Xu presented a combined differential evolution (HDE) and speed-up local search based on critical path to improve the efficiency [20].Recently, firefly algorithm attracted many researches' attention and Karthikeyan et al. proposed a hybrid discrete firefly algorithm (HDFA) to optimize the FJSP [21].In HDFA, some discrete transformations of attractiveness, distance, and movement were developed and local search was applied to improve the accuracy of the solutions of FJSP.Focusing on scheduling problems with maintenance activities, Li et al. proposed discrete artificial bee colony 2 Mathematical Problems in Engineering (DABC) [22].In DABC, an effective decoding method was studied to decode the representation of FJSP with maintenance activities and food sources of employed bees, onlooker bees, and scout bees were produced by Tabu search.Wang et al. presented an artificial bee colony to minimize the makespan [23].A population updating mechanism for generating the scout bees and local search strategy based on critical path for onlooker bees were employed in ABC to achieve good performance.To minimize the makespan and the mean of earliness and tardiness, a discrete harmony search (DHS) was developed by Gao et al. [24].In this algorithm, a discrete function to generate new harmonies was designed and some local search approaches were used to strengthen exploitation ability.To minimize makespan, Yuan et al. developed a hybrid harmony search (HHS), which employed local search process to perform exploitation [25].In HHS, a new decoding method was applied to reduce the search space and neighborhood structure of local search was based on common critical operations.
In most job shop models, processing time is assumed as a deterministic value.However, in uncertain environment, processing time is frequently acquired in a certain range.Therefore, flexible job shop scheduling model with fuzzy processing time (FJSPF) is a closer approximation to the real manufacturing system with uncertain situation and the model of FJSPF was built in some researches to address this situation.As a FJSP model with fuzzy processing time, FJSPF is also highly difficult to solve and several intelligent algorithms based on fuzzy number have been developed to optimize it.Lei presented a decomposition-integration genetic algorithm (DIGA), which provided a two-string representation and decomposed the chromosome into two parts (job sequencing part and machine assigning part), and the two parts of the population evolved independently in DIGA [28].To minimize the fuzzy makespan, Lei also proposed a coevolutionary genetic algorithm (CGA), which designed a novel representation, crossover operator, and a modified tournament selection.In CGA, job sequencing and machine assignment evolved independently and then cooperated to determine the scheduling scheme [29].To solve the FJSPF, an estimation of distribution algorithm (EDA) based on probability model was proposed by Wang et al. [30].In EDA, probability matrix was firstly generated and then updated with the elites.Moreover, roulette selection was used to produce new individuals.
In this paper, an improved version of PSO (called IDPSO), which consists of discrete operators, is introduced to solve FJSPF.First, triangular fuzzy number (TFN) is used to express fuzzy processing time and several initial methods based on TFN are developed in initialization phase.Then six initial methods are applied to initialize machine assignment and random method is used to initialize operation sequence.Then simple discrete operators are embedded in IDPSO to update the particles' positions.To guide particles by the personal best positions and the global best position in IDPSO, discrete operator  2 is applied to update the current particles' positions using the personal best positions or the global best position.Moreover, the global best position is extended to contain several global best positions in the predefined global best set and that can make the guidance more effective.
The remainder of this paper is organized as follows: Section 2 introduces the problem statement of FJSPF.Section 3 describes the implementation of IDPSO for FJSPF.The simulation results and analysis of the proposed algorithm are provided in Section 4. A conclusion is finally given in Section 5. (

FJSPF Model
For FJSPF model, processing time on any machine is fuzzy and it can be called fuzzy processing time.Generally, fuzzy processing time can be treated as TFN and is expressed as follows: where  1 ,  2 , and  3 are the minimal probable processing time, the most probable processing time, and the maximal probable processing time.The optimization task is to determine an assignment on machines and a sequence of operations to optimize one or more criteria.The maximal completion time is a significantly important factor to impact due date.Therefore, the criterion is to minimize the maximal fuzzy completion time in this paper and its expression is as follows: where   is a TFN and means the fuzzy completion time on machine .

Operations on Triangular Fuzzy
Number.Some operations for real number are hardly applied to TFN since TFN includes three values, and then special operations should be designed to operate TFNs.Therefore, researchers define three essential operations (addition operation, max operation, and rank operation) to dispose TFNs [31].Specifically, two TFNs are added by addition operation; max operation is required to obtain the maximal TFN and rank operation is applied to sort TFNs.Max operation and rank operation need to compare two TFNs.In order to compare TFNs, three criteria should be defined through a TFN  = ( 1 ,  2 ,  3 ) as follows: Criterion 1 (1): (1) Input: , Criterion 3 (3): According to the above three criteria, we can use two TFNs,  = ( 1 ,  2 ,  3 ) and  = ( 1 ,  2 ,  3 ), to demonstrate the comparing procedure as shown in Pseudocode 1.

IDPSO for FJSPF
where  1 and  2 are cognitive coefficient and social coefficient.
is inertia weight. 1 and  2 are random numbers in (0, 1). is current iteration number.The particles will obtain new promising positions until a satisfactory solution is found or the maximal iteration number is met.

Detail of IDPSO. PSO was firstly introduced by Kennedy
and Eberhart in 1995 and many variants have strong abilities to solve nonlinear and continuous optimization problems [32].FJSPF is a class of discrete optimization problem, which has many discrete variables that need to operate.Therefore, an improved PSO with several discrete operators is introduced in discrete PSO version.According to the mechanism of PSO, personal best positions and global best position are the guiders of the population and they make all particles move to global optimum.Therefore, the position x +1  of IDPSO can be manipulated according to the following: where ⊗ is probability operation, indicating that the following operator carries on with the corresponding probability.+ is an operation, which indicates that the current operator is completed and the next operator is carrying on.,  1 , and  2 are the probabilities and  1 and  2 are two discrete operators to operate discrete variables.In detail,  2 equals 1 −  1 in our algorithm. 2 is firstly applied to the current position x   and the personal best position p   .Then,  2 is applied to the current position x   and the global best position g  best .This makes the current particle guided by personal best position and global best position in IDPSO.In addition, the details of  1 and  2 are presented in Section 3.5.Perturbation operator  3 is then adopted and will also be introduced in Section 3.5.
Two positions, personal best position p  and global best position g best , are the key factors to affect particles'   convergence curve.According to PSO, the personal best position of particle  is also updated as follows: In

3.3.
Encoding and Decoding.Two-vector representation (operation sequence vector and machine assignment vector) regularly represents operation sequence and machine assignment.Numerous researches adopted operation-based representation to represent the operation sequence vector since it is simple and no illegal representation needs to be repaired.For the operation-based representation, the length of operation sequence vector equals the total number of operations.Each operation is denoted by job number and the th appearance of job number denotes the th operation of the corresponding job.For machine assignment vector, assign the machines to all operations with ascending order of job number.One example of two-vector representation is illustrated in Figure 1.For instance, the sequence [1 2 1 2 3 3] in Figure 1 1 shows the processing times of 3 jobs and 3 machines' instance.Figure 2 shows the Gantt chart of this representation scheme.

Initialization.
In IDPSO, operation sequence vector is initialized randomly and several initial methods based on TFNs are developed to initialize machine assignment.It is obvious that initial methods from FJSP cannot be directly applied to FJSPF due to the particularity of TFNs.For machine assignment, we deal with the minimal probable processing time, the most probable processing time, and the maximal probable processing time of TFN separately and  then new initial methods (Rules 1, 2, and 3) are formed.These rules can decompose a TFN into three processing times and find a promising machine assignment for each processing time separately.
Rule 1. Select all minimal probable processing times of TFNs in fuzzy processing timetable to form a new processing timetable.In the new table, find the minimal processing time for each operation, and assign the corresponding machine to the operation.Repeat the procedure until all operations are assigned.
For instance, the first element "1" of  11 in Table 2 is selected from the minimal probable processing time of the first element (1,4,5) in Table 3.In Table 3, the number "1" represents the machine assigned on the corresponding operation.
Rule 2. The procedures of Rules 2 and 3 are the same as that of Rule 1.The difference is that Rule 2 selects the most probable processing times of TFNs to form a new processing timetable and Rule 3 selects the maximal probable processing times of TFNs to form a new processing timetable.Three existing rules, global minimal processing time rule (Rule 4) [33], local minimal processing time rule (Rule 5) [34], and earliest completion machine (Rule 6) [35], were proved to be effective initial methods in FJSP and are also modified to initialize the population in this paper.Finally, empirical combination of the above six heuristic methods (10% by Rule 1, 10% by Rule 2, 10% by Rule 3, 30% by Rule 4, 30% by Rule 5, and 10% by Rule 6) is used to initialize machine assignment vector of the population.The details of the three existing rules are as follows.
Rule 4. Firstly, find the machine with the minimal TFN in the processing timetable, and assign the operation on this machine.Then, add the assigned processing time to other processing times in the same column and delete the row in which the operation is assigned.Repeat the procedure until all operations are assigned.
Rule 5.For each job, find the machine with the minimal TFN for every operation   of the job, and assign the corresponding operation on this machine.Then, add the assigned time to other processing times in the same column and delete the row in which the operation is assigned.Repeat the procedure until all operations are assigned.Rule 6.For each operation in order of operation sequence vector, assign the operation to the machine that can complete the operation with the earliest fuzzy completion time.

Designed Operators.
According to formula ( 9), the function of  1 is to make the particles unchanged and  is the unchanged probability.Otherwise, the particles are the changed particles.In IDPSO, these unchanged particles are only operated by  2 and the other particles are only operated by perturbation operator  3 .The function of  2 is to exchange the information from personal best positions and global best position with probabilities  1 and  2 , respectively.Since many promising particles are produced in the initial phase, exchanging multipoints information can significantly change the scheduling of FJSPF and may destroy the useful information of promising particles.Consequently, exchanging singlepoint information from the best positions contributes to the reservation of particles' self-information and acquiring helpful information from the best positions.Therefore,  2 adopts simple and effective single-point information exchanging method.In addition,  2 is used to the current particle's position and its personal best position with the probability  1 .Then,  2 is used to the current particle's position and the global best position with the probability  2 .It is worth pointing out that  2 is only implemented on machine assignment vector. 2 works as follows.For the particle PA 1 , keep the operation sequence unchanged.Then randomly select an element in machine assignment vector of PA 1 , and replace this element with the element of the personal best position p  or the global best position g best in the same place.How  2 works is illustrated in Figure 3.
2 with single-point information interchange has a weak exploration ability.To avoid this weakness, perturbation operator  3 is used to explore other space and can make the algorithm have a good ability to balance exploitation and exploration.Perturbation operation  3 can be executed as follows.
Method 1 and Method 2 should be defined in perturbation operator  3 .For the particle PA 1 , keep machine assignment unchanged.Randomly select an operation   in operation sequence vector and record its place pl in operation sequence vector of PA 1 , and then find the place pl  where the operation   locates in operation sequence vector of PA 2 .Finally, Execute perturbation operator f 3 on x t i f(x i ) < f(p i ) Update global best set and g best Produce the new position as follows: Figure 5: Flowchart of IDPSO algorithm.
Method 1 is exchanging the operation   with the operation in the place pl  of operation sequence vector of particle PA Figure 6: Trend of factor levels for each parameter.

Parameter Setting.
In this section, the impacts of four parameters in our algorithm are investigated and some experiments are conducted using Taguchi method.The experiment design involves four parameters: , , ( 1 ,  2 ), and Na, and each parameter is divided into four levels.Value ranges of four parameters are listed in Table 4.According to the experiment design of four factors and four levels, experiments of an orthogonal array (4 4 ) are obtained.The factor levels of each parameter are listed in Table 5, and the orthogonal experiments are listed in Table 6.According to the orthogonal table, the designed experiments are conducted on Instance 1 [28] and each designed experiment runs 20 times independently.
The maximal generation number is 1000.After the experiments, the results of 1 for each designed experiment are also listed in the last column of Table 6.According to the results of 1 in Table 6, we can obtain the average value of 1 (AVC1) for each factor level of each parameter in Table 7.For each parameter, the effect on performance is analyzed by "Rank" in Table 7.In Table 7, "Delta" denotes the maximal AVC1 minus the minimal AVC1 for each parameter, and "Rank" denotes the rank of "Delta.""Rank" reflects the significance of each parameter.For comparisons, ( 1 ,  2 ) rank first, and  rank second.Therefore, the parameters ( 1 ,  2 ) are the most significant factor on the performance of our algorithm.According to Table 7, the trends of each factor level are illustrated in Figures 6(a)-6(d).In Figure 6, it is obvious that parameter combinations of the proposed

Numerical Computation and Comparison.
Five instances ranging from 10 jobs-40 operations to 15 jobs-80 operations (Instance 1-Instance 5) [28,29] are considered to evaluate the proposed algorithm and IDPSO is compared with some state-of-the-art algorithms.The test is implemented in Matlab 7.1 on Lenovo computer with an Intel 2.4 G processor and 4 G RAM.The proposed algorithm will run 30 times on each instance.The parameters are selected as suggested in Section 4.1.The suitable parameter settings of the proposed algorithm are summarized as follows: the population size is 100; the size of global best set is 10; the maximal generation number is 1000; the values of ,  1 , and  2 are 0.94, 0.4, and 0.6, respectively; the compared algorithms are DIGA [28], CGA [29], EDA [30], PEGA [36], PSO+SA [37], and HABC [38] and their results on the five instances are from the corresponding literature.
As shown in Table 8, it lists the average value (Average value), the best solution (Best solution), and the worst solution (Worst solution) obtained by seven algorithms.For Instances 1-4, the proposed algorithm outperforms all other six algorithms in terms of "Average value," "Best solution," and "Worst solution."For Instance 5, the results of "Average value," "Best solution," and "Worst solution" obtained by the proposed algorithm rank second among these compared algorithms.Besides, the average CPU time of 20 runs is listed in Table 9.For Instances 1-4, the proposed algorithm spends the second shortest average CPU time among all seven algorithms.For Instance 5, the proposed algorithm costs the shortest average CPU time.In summary, for all five instances, the proposed algorithm is the most effective algorithm among these compared algorithms.In addition, the Gantt charts of best solution obtained by the proposed algorithm are illustrated in Figures 7(a)-7(e).
Wilcoxon rank sum test can confirm whether one data set is superior to another data set in statistics and two values ( value, ℎ value) can be acquired by the test to assess the  (5, 7, 9) (0, 0, 0) (5, 6, 9) (0, 0, 0) (7, 9, 11) (0, 0, 0) (4, 6, 9) (0, 0, 0) (7,10,16   AV(1) and STD(1) in Table 10, it is clearly seen that IDPSO obtains the best value among three algorithms for all five instances.For comparing AV(1) and AV(1) noise , STD(1) and STD(1) noise obtained by IDPSO, the results under noise conditions are slightly worse than the results with no noise.Therefore, the noise hardly damages the performance of IDPSO to some degree.Also, since all ℎ value 1 and ℎ value 2 are equal to 1, this indicates that IDPSO statistically outperforms EDA and DIGA with the probability of 95%.From the above analysis, IDPSO has a better comprehensive performance than other algorithms to optimize FJSPF and also can solve FJSPF under noise conditions.Some superiorities of the proposed algorithm are summarized as follows: Several initial methods based on TFNs and three initial methods from other papers are applied to produce abundant initial particles and more promising solutions.High-quality and high-diversity initial population can make the particles reach more promising area in the initial phase.Simple and effective discrete operators are introduced to IDPSO for exchanging the information of particles.Discrete operator  2 is to exchange the machine information from personal best positions and global best positions, and then perturbation operator  3 can find better solution from other space and increase the exploration of the population.The appropriate combination of discrete operators  2 and  3 is a valid method to balance the exploration and exploitation capabilities.Moreover, global best position is a significant factor to guide the population to the global optimum.Several global best positions are used to update the particles' positions and that contributes to correcting some global best positions' errors in search process and avoiding being trapped into local optima.Consequently, this comprehensive strategy makes IDPSO have a better performance on FJSPF.

Conclusion
In this work, an improved version of PSO with discrete operators is provided to solve FJSPF.Firstly, several heuristic initial methods based on TFNs are developed to initialize machine assignment and random method is used to initialize operation sequence.Secondly, simple and effective discrete operators are employed in IDPSO to update particle's position and perturbation operator  3 is applied to explore more area.Thirdly, several global best positions retained in global best set can make the population obtain more information about global best area and correct some global best positions' errors.Finally, the impacts of four parameters are investigated by Taguchi method and parameter selection is suggested by designed experiments on Instance 1. Furthermore, five instances of FJSPF are used to evaluate IDPSO and the comparisons with several previous published algorithms are also performed.The experimental results show that IDPSO is able to obtain better solutions on Instance 1 to Instance 5 and is a competitive method to solve FJSPF.How to address dynamic scheduling problem and establish more valid models for job shop manufacturing system needs our future efforts.

Figure 2 :
Figure 2: Gantt chart of one solution of 3 jobs and 3 machines.
several particles to search the optimal solution.Each particle has three attributes: velocity, position, and personal best position, and the population has a global best position.For a specific issue, suppose  particles are in  dimensional space.Position represents a potential solution of the issue and velocity represents the move of particle.Velocity and position of particle  can be denoted by k  = (V 1 , V 2 , . . ., V  ) and x  = ( 1 ,  2 , . . .,   ).Personal best position (denoted by p  = ( ,1 ,  ,2 , . . .,  , )) and global best position (denoted by g best = ( best,1 ,  best,2 , . . .,  best, )) represent the best position obtained by particle  and the best position obtained by all particles.The velocity k +1 3.1.Theory ofParticle Swarm Optimization.Particle swarm optimization (PSO) is an intelligent algorithm and employs a population with order to avoid the weakness of PSO being trapped into local optima, the global best position is extended to a predefined global best set which contains several global best positions in IDPSO.Suppose global best set has Na global best positions.Then the global best set updates itself as follows: add current population and global best set to form a new set, and then find Na best positions from new set to the next global best set.Finally, global best position g best is randomly selected from global best set.

Table 1 :
Processing times for 3 jobs and 3 machines instance.
of  11 in Table 1, and the second element "2" of  11 in Table 2 is selected from the minimal probable processing time of the second element (2, 3, 4) of  11 in Table 1.Therefore, Table 2 is formed by selecting all the minimal probable processing times of Table 1.And then the machine assignment is obtained by Rule 1 and is shown

Table 2 :
Minimal probable processing times for 3 jobs and 3 machines instance.

Table 3 :
Machine assignment for 3 jobs and 3 machines instance.

Table 4 :
Value ranges of four factor parameters.

Table 5 :
Factor levels of parameters.

Table 6 :
Orthogonal table of designed experiments.

Table 7 :
Rank of each parameter by 1.

Table 9 :
CPU times of the five instances.

Table 10 :
Results of statistical analysis and noise conditions of five instances.