Multiobjective particle swarm optimization with direction search and differential evolution for distributed flow-shop scheduling problem

1.   Introduction

Production scheduling plays a vital role in many manufacturing systems, and effective production planning is an important factor that enables the industry to improve production efficiency and resource utilization ^[1]. Therefore, it is of great significance to study production scheduling problems, especially to develop effective solving algorithms. Production scheduling is the allocation of limited production resources to tasks within a specific time period to optimize one or more goals determined by decision makers ^[2]. With the rapid development of manufacturing technology and equipment, scheduling in production systems has become more and more complex. Under the trend of globalization, since the operation of a single factory has been unable to fully and flexibly respond to the high changes in market demand ^[3], the structure of the factory has changed from centralized to decentralized, and the multi-shop co-production pattern has become increasingly popular. Compared with a single factory production center, distributed manufacturing is conducive to resource sharing, which can reduce costs and risks, improve product quality ^[4], improve economic benefits, and respond quickly to market changes. Among all kinds of distributed manufacturing problems, the distributed flow-shop scheduling problem (DFSP) is one of the most widely studied scheduling problems and has strong application value ^[5].

To the best of our knowledge, the DFSP with the makespan criterion is currently the most studied because it is challenging ^[6]. The DFSP can be regarded as a generalized version of flow-shop scheduling problem (FSP), which involves two sub-problems ^[7]. The first is to assign jobs to factories, and the second is to sort jobs in each factory, and these two sub-problems affect each other, so they must be solved together ^[8]. However, the two sub-problems of job sequence and factory allocation are coupled with each other, and the optimization of each sub-problem does not necessarily mean that the coupled problem is optimized. Therefore, although there have been many research results on the single-shop scheduling problem, it cannot be directly applied to the DFSP. In addition, from a mathematical point of view, the DFSP is an NP-hard problem, and the increase of sub-problems increases the solution space, which in turn makes the problem more difficult to solve. To sum up, DFSP is much more difficult to solve than single-shop scheduling, and has strong research value, requiring researchers to develop new methods.

Today, there are more and more ways to solve DFSP. Traditional mathematical methods, such as enumeration, branch and bound, are suitable for small, simple problems due to their high computational complexity. Therefore, heuristics and meta-heuristics are used to solve complex combinatorial optimization problems. In recent years, existing optimization algorithms include genetic algorithm (GA) ^[9], particle swarm optimization algorithm (PSO) ^[10,11], tabu search ^[12], distribution estimation algorithm ^[13], hybrid immune algorithm ^[14], scatter search ^[15], chemical reaction optimization ^[16] and iterative greedy algorithm ^[17], artificial bee colony algorithm (ABC) ^[18], water wave optimization algorithm ^[19], and cooperative optimization algorithm ^[20]. These algorithms work differently than traditional mathematical methods. The unique algorithm theory and operation mechanism guarantee that the exploration and utilization of the algorithm reach a certain balance in the iterative process. In addition, addressing other objectives of DFSP such as the total processing time criterion ^[21] is also considered.

To better address two closely related sub-problems of DFSP, two industrial indicators are proposed ^[21]. One indicator is the most studied makespan ^[22] and the other is total processing time. Minimizing the total processing time can reduce the job in progress and lead to stable utilization of resources. The advantage is even higher in distributed environments where the imbalances among resources and jobs in progress may increase due to the existence of several factories. Total processing time has been pointed out as a more relevant and significant objective for today's dynamic manufacturing environment. Therefore, this study addresses the DFSP in order to minimize makespan and total processing time criteria. When solving the problem, one of these two indicators can be biased, so the problem model can be clarified, and it is very convenient to solve the DFSP. But this method can only take care of one indicator, it is difficult to optimize another indicator. Another method that can be used to solve these two closely related problems is to arrange the process sequence of each job completely randomly and assign it to the factory randomly ^[16]. This method can take into account both the makespan indicator and the total processing time indicator. Therefore, we use this method to solve DFSP.

The evolutionary algorithm may be not as accurate as traditional mathematics methods, but its calculation speed has been greatly improved. Therefore, evolutionary algorithm is more popular than traditional mathematical method. PSO and GA are more classic algorithms in evolutionary algorithm. Despite that PSO has some similar properties with GA, PSO does not apply evolution operators like mutation and crossover. PSO has a simple structure, is easy to calculate, and has obvious advantages in calculation speed. The scale of DFSP may be very large, and the calculation speed of PSO is still very fast in dealing with complex problems, so PSO is a good choice to solve DFSP.

PSO is a swarm intelligence optimization algorithm developed from bird foraging. PSO will first initialize a particle swarm, each particle will search the solution space according to its own historical optimal and global optimal, and finally find an excellent solution ^[11,23]. In this process, the algorithm does not need much calculation, so its calculation speed has an advantage in many evolutionary algorithms. PSO was first proposed to solve the problem of continuous solution space. The solution space of DFSP is discontinuous, which needs to be solved with a discrete PSO ^[24,25,26]. The structure of discrete PSO is too simple, and the disadvantage of discrete PSO is that it is very easy to fall into local optimum. Therefore, it is necessary to add a mechanism to prevent the local optimum ^[27,28].

Differential evolution (DE) algorithm is also a kind of evolutionary algorithm. The calculation quantity of the algorithm is small and the calculation speed is fast. DE is based on the difference between particles to evolve bad individual to excellent individual ^[29]. Therefore, DE can strengthen the connection between particles in a particle swarm and is suitable for improving the PSO.

DFSP usually have multiple constraints ^[30], and single-objective PSO is difficult to solve such problems. Therefore, a multiobjective particle swarm optimization algorithm (MoPSO) is studied. Since too many constraints on the DFSP will affect the performance of the algorithm, it is necessary to improve MoPSO ^[31].

In this study, a multiobjective particle swarm optimization (MoPSODS-DE) based on enhanced direction search and differential evolution is proposed. The enhanced direction search strategy explores the particle swarm in multiple directions of the Pareto front (PF), and speeds up the convergence of the upper, central and lower areas of the PF. The differential evolution strategy performs differential operations on the updated particle swarm to improve the quality of the particles and prevent falling into local optimum. Therefore, combining these two strategies can further improve the convergence performance and distribution performance of the algorithm, as well as improve the diversity of the particle swarm. In addition, a method of two-vector representation of particles is designed for DFSP for encoding and decoding, and the particles are updated using exchange sequences, crossover, and mutation. The proposed MoPSODS-DE are compared with NSGA-Ⅱ, SPEA2, MoPSO, MOEA/D, a fast multiobjective hybrid evolutionary algorithm (MOHEA) and hybrid multiobjective evolutionary algorithm based on DE (MOHEA/DE). The results show that MoPSODS-DE not only has excellent convergence performance, but also can guarantee that the obtained solutions have good distribution performance and diversity.

The remainder of this study is organized as follows. In Section 2, related works about solving the DFSP is shown. In Section 3, the DFSP formulation is introduced. Section 4 describes the MoPSODS-DE algorithm. In Section 5, experimental data of MoPSODS-DE are presented and discussed. Finally, Section 6 concludes this study with some conclusions and future work.

2.   Related work

The DFSP is a variant of FSP, and the characteristic of this type of problem is that if the scale of the problem becomes larger and larger, the difficulty of solving it will increase exponentially. Naderi and Ruiz as pioneers in introducing DFSP, established 6 mathematical models, developed 2 factory assignment rules, and proposed 2 variable neighborhood algorithms ^[32]. Furthermore, effective constructive heuristics and meta-heuristics are proposed to solve it. Since then, other authors also proposed many high-performance methods to deal with this problem.

In many realistic FSP, multiple interrelated objectives need to be considered, such as makespan, total processing time, machine idle time, and total delay ^[33]. Therefore, it becomes a trend for DFSP to consider multiple objective functions. For example, Rifai et al. ^[34] dealt with a distributed reentrant PFSP considering three objective functions: minimization of makespan, minimization of the total cost, and minimization of expected delays. These authors developed a multiobjective adaptive large neighborhood search algorithm to estimate the PF. In addition, Deng and Wang ^[35] proposed a competitive memetic algorithm for solving multiobjective DFSP considering the minimization of makespan and total delay. Lu et al. ^[36] established a new mathematical model for distributed energy-saving FSP with the objective of maximum lifetime and minimum total energy consumption, and designed a hybrid multiobjective optimization algorithm that combined iterative greedy and efficient local search, which provides a set of compromise solutions for this problem.

Furthermore, besides the classic DFSP, other types of distributed production scheduling problems have been studied in recent years. Most of them are extensions of the classic DFSP, such as distributed no-wait FSP ^[37], distributed no-idle permutation FSP ^[38], distributed hybrid no-idle FSP ^[39] and distributed finite buffer FSP ^[40] consider no-wait constraints, no-idle constraints, hybrid no-wait constraints, and finite buffer constraints in DFSP, respectively. Others come from extensions of other classic production scheduling problems in multi-factory production environments, such as distributed two-phase scheduling problems, distributed reentrant FSP ^[34], and distributed assembly permutation FSP ^[7]. For example, Zhao et al. proposed a two-stage cooperative evolutionary algorithm with problem-specific knowledge to solve the energy-efficient scheduling problem of the no-waiting FSP, whose criterion is to minimize the makespan and total energy consumption ^[41]. To solve the scheduling problem of distributed assembly without idle flow shop, Zhao et al. proposed a collaborative water wave optimization algorithm to solve the distributed assembly no-idle FSP, with the objective of minimizing the maximum assembly completion time ^[42]. In addition, Zhao et al. proposed a self-learning discrete Jaya algorithm to solve the energy-saving distributed no-idle FSP in heterogeneous factory systems with the criterion of maximizing minimize total tardiness, total energy consumption and factory load balancing ^[43].

GA is the most classic evolutionary algorithm, which is used in various complex problems. Jia et al. studied the DFSP under different criteria and employed a standard GA to solve the problem ^[44]. Zhang et al. proposed an enhanced GA for the distributed assembly permutation FSP, and a crossover strategy is designed based on local search to speed up convergence ^[45]. Wang and Wang proposed a mixed integer linear programming model and a novel iterative greedy algorithm to solve DFSP, and designed a local intensive method based on variable neighborhood search ^[46]. In this paper, the two operations of crossover and mutation in differential evolution algorithm are used to update the factory sequence of particles in the particle swarm.

PSO is a kind of swarm intelligence algorithm, which has a simple structure and fast calculation speed, and is often used in various NP-hard problems. For PSO, Tasgetiren et al. studied improved variable neighborhood search PSO, finding many optimal solutions for the top 90 Taillard benchmark instances ^[47]. On the other hand, Pan et al. ^[48] applied discrete PSO to solve the no-wait flow-shop scheduling problem. DFSP usually have multiple constraints ^[30], and single-objective PSO is difficult to solve such problems. To this end, a MoPSO is studied, which mainly solves multiobjective problems and is more practical than PSO. Since too many constraints on DFSP will affect the performance of the algorithm, it is necessary to improve MoPSO. Rahimi-Vahed and Mirghorbani ^[49] studied a PSO method targeting weighted average completion time and weighted average delay. The proposed MoPSO is compared with a multiobjective GA. The algorithm outperforms multiobjective GA on some specific performance metrics. Su and Chi ^[50] proposed a MoPSO that combines PSO with DE and utilizes a simulated annealing strategy to generate a scale factor to dynamically tune PSO and DE. Cui et al. ^[51] proposed a hybrid MoPSO, which uses different operators to generate new offspring, and then uses the environmental selection mechanism to screen out the next generation of individuals. To improve the efficiency of solving a bi-objective mixed no-idle FSP that minimizes the makespan and total processing time, Zhang et al. proposed a multiobjective particle swarm optimization based on multi-direction update ^[52].

ABC is a swarm intelligence algorithm like MoPSO, which is also widely used in solving DFSP. Li et al. proposed a hybrid ABC to solve DFSP^[53]. In order to enable decision makers to conveniently choose the final solution according to actual needs, the purpose is to obtain a set of optimal solutions on PF, rather than finding a single solution ^[35]. Pan et al. used total processing time as the optimization goal of DFSP, designed different heuristic algorithms and improved ABC ^[54]. Huang et al. proposed three constructive heuristic methods and effective discrete ABC to solve the DFSP ^[55]. To solve the multiobjective distributed fuzzy permutation FSP, Emin et al. developed an artificial bee colony (ABC) algorithm ^[56].

However, whether it is MoPSO, ABC or other swarm intelligence algorithms, there is a common disadvantage, that is, it is difficult to find the global optimal solution. Because the DE has a single structure, simple operation, and can reflect the relationship between two individuals, so many researchers also use DE to solve scheduling problems. To address the electric vehicle charging scheduling (EVCS) problem and consider the dependencies between site selection, charging options at each site, and charging capacity settings. Liu et al. designed a mixed variable differential evolution as a scheduling algorithm for EVCS systems ^[57]. Zhou et al. developed an adaptive differential evolution algorithm for a single batch machine scheduling problem with varying release times and job sizes ^[58]. Zhang et al. proposed a HMOEA/DE to solve the FSP with the criteria of minimizing the maskspan and tardiness simultaneously. The DE strategy is used as a local search mechanism, which further improves the convergence and distribution performance of the algorithm on elite population ^[59]. For the FSP with blocking, Han et al. proposed a new discrete ABC, which was combined with DE and solved using the makespan criterion ^[60]. This paper adopts the combination of improved MoPSO (MoPSODS) and DE search to further improve the performance of the algorithm in all aspects.

The DFSP can be solved with the above algorithms. By comparison, the search speed of PSO algorithm is faster and the PSO algorithm based on enhanced direction search is better than other algorithms. Since the PSO uses pbest and gbest to tune search direction, while there is no operation between individuals. Therefore, combining the DE search with the particle swarm algorithm based on enhanced direction search can further improve the convergence performance of the algorithm and the distribution and diversity of solutions.

3.   Problem formulation

The DFSP is described as follows: suppose the problem has $n$ jobs, $m$ machines, and $F$ factories. The machines in the $F$ factories are the same, and each job can be processed independently in each factory. When a job is assigned to a factory, all the processes of the job are processed within that factory without being transferred to other factories. Each job needs to be processed on $m$ machines. Each processing of each job can only be completed on one machine, and each machine can only process one job at a time. It does not stop after each operation is started. No need to consider the preparation time before processing, and finally find a job order to arrange the production of the factory ^[54].

The mathematical model of DFSP is expressed in the following notations.

Indices

$i$: index of job, ($i = 1, 2, \ldots, n$);

$j$: index of machine, ($j = 1, 2, \ldots, m$);

$f$: index of factory, ($f = 1, 2, \ldots, F$).

Parameters

$n$: the number of jobs;

$m$: the number of machines;

$F$: the number of factories;

$n_{f}$: the number of jobs in the $f$-th factory;

$t_{(i, j, f)}$: the processing time of the job located on the $i$-th position in the processing sequence in the $f$-th factory, on the machine $j$;

$C_{(i, j, f)}$: the completion time of the job located on the $i$-th position in the processing sequence in the $f$-th factory, on the machine $j$.

$X_{(i, k, f)}$: the binary variable, if job $i$ occupies the position $k$ in the $f$-th factory $X_{(i, k, f)}$ = 1. Otherwise, $X_{(i, k, f)}$ = 0.

Mathematical model

The mathematical model for minimization of makespan and total processing time can be expressed as follows:

In this paper, two common objectives are used, one is the makespan and the other is the total processing time, and then minimize the two objectives. Constraint (3.1) guarantees that each job is only assigned to one factory and only appeared once in the factory. Constraint (3.2) explains that the possibility for n positions where all jobs should be allocated to each factory. Constraint (3.3) controls and guarantees that a job is both a predecessor and successor of another job at the same time. Constraints (3.4)–(3.6) illustrate that the completed time of the first job processed on the first machine in each factory, which is guaranteed that the start time of the first machine in each factory is 0. Constraints (3.7) and (3.8) are that each job is only start to process after the previous job processed on the machine in an identical factory, which guarantees that multiple jobs are not processed on one machine. Constraint (3.9) means that the makespan must be greater than or equal to the makespan of each job. Therefore, Eq (3.10) can be used to calculate the makespan of each job. The makespan and total processing time of DFSP can be obtained by Eqs (3.11) and (3.12), respectively. Figure 1 shows the Gantt chart of DFSP. Taking two factories and five jobs as an example, the values of the two objectives of DFSP can be obtained through the calculation of the above formula. The makespan of factory 1 is 20 and the makespan of factory 2 is 22, so the makespan of DFSP is 22. In addition, the total processing time of DFSP is determined by the makespan of the jobs in the factory. The total processing time of factory 1 is $12 + 20 = 32$, and the total processing time of factory 2 is $12 + 20 + 22 = 54$, so the total processing time of DFSP can be obtained as 54.

4.   Solution algorithm

4.1.   Overview of the MoPSODS-DE

The MoPSODS-DE algorithm is improved on the basis of the MoPSO. The first step is to consider the encoding and decoding strategies and adapt the MoPSO algorithm for the DFSP. Secondly, the enhanced direction search strategy is used to adapt the MoPSO to greatly improve the convergence performance without reducing the distribution performance. Third, the above improved algorithm, namely MoPSODS, is mixed with differential evolution search to improve the local convergence ability of the algorithm and prevent the algorithm from premature convergence. Then use the method of two-vector representation of particles to perform reasonable and effective update operations on the particles, which further enhances the adaptability of the algorithm. The algorithmic framework of MoPSODS-DE is shown in Figure 2.

MoPSO is the same as the ordinary PSO algorithm, without grouping the particle swarm. In order to strengthen the multi-direction search ability of MoPSO in PF, the PDDR-FF fitness function and VEGA strategy are introduced to group the particle swarm, thereby further improving the strong convergence ability of PF.

Using PDDR-FF fitness function and VEGA strategy, the particles are divided into three sub-particle swarms ^[61]. The VEGA strategy is used to select two single objectives, and the PDDR-FF is used to select particles in the central area of the PF. The three sub-particle swarms performed well on objective 1, objective 2 and central area of the PF respectively. The sum of the number of particles in these three sub-particle swarms is equal to the size of the particle swarm.

Then, use the PSO formula to update the velocity and position of the particles. Update particles use a combination of exchange sequence, crossover, and mutation to operate on two vectors respectively. The sub-particle swarm that perform well on objective 1 is updated towards objective 1, the sub-particle swarm that perform well on objective 2 is updated towards objective 2, and the sub-particle swarm near the central area of PF is updated towards the PF. In this way, each sub-particle swarm will be updated in a different direction. The updated particles are mixed into a new particle swarm, which is a temporary particle swarm, and its size is the same as that of the particle swarm.

Apply the differential search strategy to the temporary particle swarm. Randomly select two particles to compare the advantages and disadvantages, and then perform the difference operation to generate a new particle. If the size of the particle swarm is N, then the difference operation is performed N times, and there are a total of N difference particles. Then select N excellent particles from the difference particles and the temporary particle swarm as the new generation particle swarm.

The MoPSODS-DE algorithm can be divided into the above steps. The pseudo-code of MoPSODS-DE is shown in the algorithm 1.

4.2.   Encoding and decoding

Encoding and decoding is the first step in the algorithm, which is used to solve compatibility problems ^[62]. The DFSP requires two-vector to represent a solution. One is that the job sequence is represented by job sequence (JS) ${\textbf{v}}_1$, and the other is that the factory vector is represented by factory allocation (FA) ${\textbf{v}}_2$. As shown in Figure 3. The first vector contains a simple arrangement of 10 jobs. Each number in the vector is the job number. In the case of only 3 factories, there can only be 3 numbers in the second vector, each number is the number of the factory. These two vectors constitute a particle. This encoding method can effectively solve the two key problems of DFSP. At the same time, transforming these two vectors can not only achieve the assignment of jobs, but also achieve the sorting of jobs within the factory ^[6].

In the case of using the above encoding scheme, the decoding scheme is very simple. The first vector is the total sequence of jobs. When decoding, first select the jobs of each factory according to the second vector, and then sort the jobs according to the first vector. After the sequence is determined, the entire processing flow is determined, and the processing schedule of the machine can be calculated by using the Eqs (3.6)–(3.12). As shown in Figure 4, three JS sequences will be obtained after decoding, each of which represents the JS of a factory, so the two objective values for the problem can finally be calculated.

4.3.   Two strategies of the MoPSODS-DE

MoPSODS-DE is based on the MoPSO, which adopts an enhanced direction search strategy for mixed sampling of particle swarm, and a local search for the updated particle swarm through a differential evolution search strategy. The enhanced direction search strategy can speed up the problem solving, and the differential evolution search strategy can further improve the quality of the solution. This section describes two strategies in detail.

4.3.1.   Direction search strategy

Enhanced direction search plays a very important role in this algorithm, which is an ingenious search strategy. It sets the search direction according to the PF, making the algorithm search in the specified direction, which can speed up the algorithm search. Enhanced direction search strategy requires two selection strategies, one selection strategy is PDDR-FF and the other is VEGA. PDDR-FF is a new particle evaluation method based on dominant and non-dominant relationships. PDDR-FF can analyze whether each particle in the current particle swarm is excellent, and its principle is also very simple and easy to implement ^[61].

In this algorithm, PDDR-FF is used to filter out the particles in the center area of the PF. These particles are excellent in each objective, and can be used to guide other particles to evolve toward the center area of the PF. The mathematical formula of PDDR-FF is as follows,

where s is the evaluated particle, q(s) is the number of particles that dominate s, p(s) is the number of the particles which are dominated by s. When the particle is not dominated, the q(s) value is 0, and the $eval$ value is always less than or equal to 1. When the particle is not dominated by other particles, the p(s) value is 0, and the $eval$ value is always greater than or equal to 1. By comparing the $eval$ value, the particles on the first PF can be selected.

In most cases, the multiobjective problem will be decomposed into multiple single-objective problems to be solved. VEGA search technology is used to decompose multiobjective problems into multiple single-objective problems, and can find excellent solutions under one objective. An excellent solution is a sub-particle swarm under each objective. The role of VEGA is to find particles that satisfy a single objective, and form sub-particle swarm. When there are three objectives, there are three sub-particle swarms.

As shown in Figure 5, the particle swarms found by the VEGA search technology are in the edge area of the PF. The particles of these sub-particle swarms can guide the particle swarm towards the edge area of the PF, which guarantees that the searched solution meets each objective as much as possible. In general, the particle swarm is divided into different sub-particle swarms by the enhanced direction search strategy, and each sub-particle swarm is updated in specific directions, making the algorithm converge faster. The specific directions and corresponding sub-swarms are to divide sub-swarm 1 to the upper edge of the PF, sub-swarm 2 to the lower edge of the PF, and sub-swarm 3 to the central area of the PF.

4.3.2.   Differential evolution search strategy

This section will describe the differential evolution search strategy in detail, and explain the combination of the differential evolution search strategy and the MoPSO.

The differential evolution search strategy is an inspiration from the DE. The DE is an evolution algorithm based on population, and the core of the algorithm is the differential idea. The difference idea is to find the difference between two particles, and then use this difference to evolve the less excellent particle to make it better. The point of differential evolution search strategy is also the idea of difference. The differential evolution search strategy applies to each generation of particles. First, two particles in the current particle swarm are randomly selected, then the difference is calculated, and finally the difference is used to evolve the non-excellent particle. Allow the current particle not to do the differential strategy, and allow the particle to be selected multiple times. Random selection is to not interfere excessively with the operation of the algorithm, making the algorithm more robust. When two identical particles are selected, in order to save time, no differential operation is performed, because there is no difference, and even if the differential operation is performed, it is the same. Allowing the particles to be selected multiple times is for a more random differential strategy to guarantee that the differential strategy can run efficiently.

The difference formula is as follows. ${\textbf{x}}_1$ and ${\textbf{x}}_2$ are the selected pair of particles, where ${\textbf{x}}_1$ is better than ${\textbf{x}}_2$, and the generated difference particle is ${\textbf{x}}$. $r$ is a random number between 0 and 1, this study makes $r$ = 1. $P$ is the DE probability, $rP$ is the difference probability, and its range is between 0–1, which is used to control the degree of difference. The larger the $rP$, the closer ${\textbf{x}}_2$ is to ${\textbf{x}}_1$.

The diagram of the differential evolution search strategy is shown in Figure 6. The characteristics of the PSO determine that the particles do not use the relationship between particles within the particle swarm when they are updated, so a differential evolution search strategy is needed to strengthen the connection between the particles. The differential evolution search strategy makes the particles not too scattered when searching for solutions, thus speeding up the search speed of the algorithm and improving the convergence performance of the algorithm.

The pseudo-code of the differential evolution search strategy is as follows.

4.4.   Update operation of particles based on two-vector representation

Since the solution space of DFSP is not continuous, the sorting method is the simplest method to solve DFSP. The encoding and decoding methods are based on the sorting method. The proposed method of representation particles based on two-vector is to perform better encode and decode of the DFSP, and to update the individual more appropriately. In addition, the details of the algorithm updating the particles are described below.

Marco et al. proposed an algebraic PSO for the permutations search space ^[63]. This algorithm is improved from the classical PSO and is used to solve sequence problems. This algorithm rewrites the formula of the PSO algorithm to be suitable for solving the permutations search space, and gives some explanations. However, the algorithm formula remains basically unchanged, and the performance of the algorithm needs to be improved for the permutations search space.

Exchange sequence is a very important part here. The solution space of the permutation search space problem is a sequence, and the algorithm is always calculating the sequence. The exchange sequence plays an important role in the sequence calculation process, it can represent the difference between particles, and it is very convenient to use, and does not require other operations.

There are also some problems when using exchange sequences. The exchange sequence is the difference between the two sequences, and it represents the conversion step from one sequence to another sequence, so an exchange sequence is only related to the current two particles, and it cannot be applied to the third particle. Therefore, the use of the exchange sequence has the following restrictions. The order of the exchange elements in the exchange sequence cannot be changed, and it cannot be read at random. It can only be read in order.

In order to better solve the DFSP, the algebraic PSO algorithm and exchange sequence will be used in combination, and the algorithm formula will be improved. The improved PSO formula is as follows:

The initial velocity of traditional PSO is determined by the inertia of the particles, but it is not applicable to the exchange sequence. Because of the existence of the exchange sequence, the previous speed is only applicable to the previous update process. In this study, the particle velocity update is divided into three parts, and rs is a randomly generated exchange sequence. rs is used to replace the initial speed in the original PSO formula, which can be expressed as Eq (4.3). Its role is to change the search direction of the particles and strengthen the search range of the particles. $pbest - {{\textbf{x}}_1}(t)$ can get the exchange sequence of the current particle and $pbest$. $gbest - {{\textbf{x}}_2}(t)$ can get the exchange sequence of the current particle and $gbest$. ${{\textbf{v}}'}(t + 1)$, ${{\textbf{v}}''}(t + 1)$, and ${{\textbf{v}}'''}(t + 1)$ are three exchange sequences corresponding to different reference particles for particle update, which can be represented by Eqs (4.3), (4.5) and (4.7) respectively. ${{\textbf{x}}_1}(t + 1)$ and ${{\textbf{x}}_2}(t + 1)$ are intermediate particles, which are used to update the current position of the particles, which can be represented by Eqs (4.4) and (4.6). Equation (4.8) shows that ${{\textbf{x}}}(t + 1)$ is the new generation particles.

In this study, aiming at the encoding of the problem, the particles are updated by crossover and mutation of the FA with a certain probability. There are three kinds of factory selection: 1, 2, and 3, which respectively represent three different factories, and then randomly generate a random decimal less than 1. If the random number is less than the crossover rate, a random single-point crossover will be performed in the current FA, otherwise, no crossover will be performed. At the same time, it is also necessary to randomly generate a decimal number less than 1. If this number is less than the mutation rate, the particles will be updated randomly by mutation, otherwise, no mutation will be performed. However, the update range of the serial number is still between 1, 2, and 3.

Therefore, the specific operation for the PSO update formula can be shown in Figures 7–9 as examples. Figure 7 corresponds to Eqs (4.3) and (4.4), which is also the first step of PSO. For the JS of the particle, suppose the number of exchange pairs in the exchange sequence is 1, and rs = (3, 5). That is, the actual executing sequence is (3, 5), the new sequence ${{\textbf{x}}_1}(t + 1)$ ($v_1$) is $\{ J_3$, $J_1$, $J_9$, $J_5$, $J_4$, $J_2$, $J_6$, $J_8$, $J_7$, $J_{10} \}$. And then for the FA of the particle, perform a random exchange sequence of (5, 10) on the ($v_2$) sequence, and then perform random mutation, and perform mutation on the sequence number 8. Finally, the job and FA have been changed, ${{\textbf{x}}_1}(t + 1)$ ($v_2$) is $\{1, 2, 1, 3, 3, 1, 3, 1, 2, 2\}$.

Figure 8 corresponds to Eqs (4.5) and (4.6), which is also the second step of PSO. For the JS of the particle, suppose $c_1$ = 1.0, $r_1$ = 0.4, the difference between the particle in this state and the $pbest$ particle is that the number of exchange pairs in the exchange sequence is 2, and 1.0 $\times$ 0.4 $\times$ 2 = 0.8, so the number of exchange pairs of executing the exchange sequence is 0. That is, the actual executing sequence is empty, and the new sequence ${{\textbf{x}}_2 }(t + 1)$ ($v_1$) is the same as the old individual ${{\textbf{x}}_1}(t + 1)$ ($v_1$). For the FA of the particle, suppose the starting sequence number of the single point crossover starts from 8, and the random mutation is mutation at the sequence number 6, then the ${{\textbf{x}}_2 }(t + 1)$ ($v_2$) sequence after the operation is completed is $\{1, 2, 1, 3, 3, 2, 3, 3, 2, 1\}$.

Figure 9 corresponds to Eqs (4.7) and (4.8), which is also the last step of PSO. Suppose $c_2$ = 1.0, $r_2$ = 0.4, for the difference between the JS in this state and the $gbest$ sequence, the number of exchange pairs in the exchange sequence is 2, and 1.0 $\times$ 0.4 $\times$ 3 = 1.2, so the number of exchange pairs of executing the exchange sequence is 1. That is, the actual executing sequence is (1, 4). Finally, for the FA crossover and mutation, the sequence number for random single-point crossover is 6, and the sequence number for random mutation is 4, so that a new individual particle can be obtained. The two sequences ${{\textbf{x}}}(t + 1)$ ($v_1$) and ${{\textbf{x}}}(t + 1)$ ($v_2$) of the particle are $\{ J_5$, $J_1$, $J_9$, $J_3$, $J_4$, $J_2$, $J_6$, $J_8$, $J_7$, $J_{10} \}$ and $\{1, 2, 1, 2, 3, 1, 2, 1, 2, 3\}$, respectively.

Since MoPSODS has updated two vectors of the particle, in order to guarantee the quality of the particle and the convergence performance of the algorithm, the difference operation does not operate on the factory vector of the particle. Therefore, the difference operation is only updated for the workpiece vector, and the factory sequence of the factory vector remains unchanged, so that the bad particles can move to the good particles more accurately.

The Figure 10 shows an example of an update operation on particles using the difference operation. On the one hand, the difference between particles $x_2$ and $x_1$ in the JA vector is ($J_2$, $J_6$), ($J_5$, $J_8$). Therefore, the JA vector of $x_2$ is improved in the direction of $x_1$ by way of exchange sequence, that is, the positions of $J_2$ and $J_6$, and the positions of $J_5$ and $J_8$ are exchanged. Suppose r = 1.0, P = 0.8, and the length of exchange sequence is 2, and 1.0 $\times$ 0.8 $\times$ 2 = 1.6, so the length of executing the exchange sequence is 1. That is, the actual executing sequence is ($J_2$, $J_6$), the new individual $x$ (JS) is $\{ J_3$, $J_2$, $J_4$, $J_8$, $J_9$, $J_1$, $J_6$, $J_5$, $J_{10}$, $J_7 \}$

5.   Experimental results

This section analyzes the performance of the algorithm from an experimental perspective.

Naderi and Ruiz presented 720 large instances, increasing the number of factories $F$ to the known 120 Taillard instances ^[64]. Data sets range from 20 jobs and 5 machines to 500 jobs and 20 machines. Small instances are easily solved by most existing DFSP meta-heuristics and are therefore not used in this work. This study will use large instances and choose data sets where the number of factories is 3. The name of the problem consists of the number of jobs, the number of machines and the number of factories. 20_5_3 is a problem that contains 20 jobs, 5 machines, and 3 factories, and $500_20_3$ is a problem that contains 500 jobs, 20 machines, and 3 factories.

The proposed algorithm MoPSODS-DE will be compared with NSGA-Ⅱ, SPEA2, MoPSO, MOEA/D, MOHEA, and MOHEA/DE. In the experimental part of this study, indicators hypervolume ($HV$), generational distance ($GD$), inverted generational distance ($IGD$), Spread, and coverage ($C$) are used to evaluate various performances of the algorithm. The results of significance analysis for algorithm performance are stored in the table, and excellent solutions will be highlighted in bold font style. In the table, "*" indicates that the comparison result is not significant, "+" indicates that our algorithm is better than other algorithms, and "-" indicates that our algorithm is worse than other algorithms.

The $HV$ index evaluates the dominated space of the obtained solution set ^[65]. The larger the $HV$ value, the better the convergence performance of the algorithm. The reference point of $HV$ in this experiment is set as combining the worst fitness values in each objective, which found by all algorithms run multiple times. The GD index and the IGD index are very similar. They both evaluate the distance between the solution and the PF ^[66]. The smaller values of GD and IGD might mean that the algorithm has better convergence performance and distribution performance. GD is an index to evaluate the convergence performance of the algorithm, and the IGD index can evaluate the convergence and distribution performance of the algorithm at the same time. A smaller value of GD means that the algorithm has better convergence performance, while a smaller value of IGD means that the algorithm not only has better convergence performance but also better distribution performance. The Spread index is an index to evaluate the diversity of the solution set generated by the algorithm, and to judge whether the distribution of the solution is widespread ^[59]. The smaller the Spread value, the better the distribution performance and the more diverse the solution set.

5.1.   Parameter settings

In this section, the main purpose is to find the optimal parameter settings for the velocity and position update formulas and the differential operation formulas in MoPSODS-DE. Because the algorithm is an improvement for MoPSO, so the parameters of MoPSO is first designed to provide a reliable parameter range for MoPSODS-DE parameters, and then the parameters of MoPSODS-DE are designed. The exchange sequence is used to represent the difference between particles, so the values of $c_1$ and $c_2$ are both 1. Three parameters of MoPSO is first considered: the range of $r_1$ and $r_2$, the crossover rate ($v_2$), and the mutation rate ($v_2$). The levels of these three parameters are set as follows: the range of $r_1$ and $r_2$ at eight levels (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, and 1.0), the crossover rate ($v_2$) at five levels (0.1, 0.2, 0.3, 0.4, 0.5), and Mutation rate ($v_2$) at eight levels (0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08 and 0.1). There are 320 (8 x 5 x 8) different combinations of the three parameters. The design of experiments (DOE) is performed on instances with 100_10 (100 jobs and 10 machines), and the instance ran 30 independent replicates ^[67]. Therefore, a full factorial experimental design with 9600 (320 x 30) treatments is verified to decide parameter settings. Particularly, the size of the three subparticle swarms is set to one-third of the entire particle swarm size. Figure 11 shows the main effect of each parameter, where the best parameter settings are the ranges of $r_1$ and $r_2$ = 0.5, crossover rate ($v_2$) = 0.1, and mutation rate ($v_2$) = 0.02.

Then, design for the parameters of MoPSODS-DE. The levels of the four parameters related to the algorithm are set as follows: the range of $r_1$ and $r_2$ at three levels (0.4, 0.5, 0.6), the crossover rate ($v_2$) at three levels (0.1, 0.2, 0.3), the mutation rate ($v_2$) at eight levels (0.01, 0.02, 0.03, and 0.04), and five levels (0.2, 0.4, 0.6, 0.8 and 0.1) for DE probability $P$. There are 180 (3 x 3 x 4 x 5) different combinations of the four parameters. The full factorial experimental design parameters were still used. Figure 12 shows the main effect of each parameter, where the best parameter settings are the ranges of $r_1$ and $r_2$ = 0.4, crossover rate ($v_2$) = 0.1, mutation rate ($v_2$) = 0.01, and DE probability $P$ = 0.8.

Because the parameters not only include the ranges of $r_1$ and $r_2$, crossover rate ($v_2$), mutation rate ($v_2$), and DE probability $P$, but also population size, elite population size, sub-particle swarm1 size, sub-particle swarm2 size, sub-particle swarm3 size, crossover rate ($v_1$) and mutation rate ($v_1$) required by other algorithms, and maximum generation. Therefore, this study uses the full factorial experiments to design the parameters for all comparison algorithms of NSGA-Ⅱ, SPEA2, MoPSO, MOEA/D, MOHEA, and MOHEA/DE. After the full factorial experiments for each algorithm, the optimal parameter combination of each algorithm after optimization is obtained. The design of all parameters corresponding to all algorithms is shown in the Table 1.

For all algorithms in this study, the population size and the maximum generation are fixed, and the performance of each algorithm can be compared when the algorithm runs to the maximum generation. MOHEA and MOHEA/DE have sub-populations, and MoPSODS-DE have sub-particle swarms. The MOHEA and MOHEA/DE has two sub-populations and an elite population, while the MoPSODS-DE has three sub-particle swarms, and the sub-population size is set according to the population size, and the sum of the three sub-particle swarm sizes is equal to the particle swarm size.

5.2.   Time complexity analysis of algorithms

In the MoPSODS-DE algorithm, assuming the size of the particle swarm is $n$. Since the encoding method of two vectors is used, jobs and factories need to be one-to-one correspondence to calculate the makespan and the total processing time during decoding. Therefore, the computational complexity of calculating makespan and total processing time is $O(n^2)$ respectively. In the initial stage, the particle swarm is initialized, and the computational complexity of setting the historical optimal individual and the global optimal individual is $O(n^2)$ and $O(n)$, respectively. When updating the particle swarm, the time complexity of swarming the particle swarm is $O(n)$ and the time complexity of updating the particle position is $O(n^2)$. In the local search stage, two particles are randomly selected, and the time complexity of the differential evolution update of the particles is $O$($n^2$). Finally, when replacing old and new particles, the time complexity of quickly sorting the particle swarm according to the fitness function value of the particles and select the first $n$ particles to enter the next iteration is $O$($n^2$+$2nlog4n$). In summary, in $M$ iterations, the time complexity of MoPSODS-DE is shown as follows:

Also for DFSP, it can be known that the time complexities of NSGA-Ⅱ, SPEA2, MoPSO, MOEA/D, MOHEA, and MOHEA/DE are $O$($n^3$), $O$($n^2$logn), $O$($n^2$), $O$($n^2$), $O$($n^2$), and $O$($n^2$), respectively.

5.3.   Effectiveness of exchange sequence

The encoding method based on two vectors to represent particles uses integer encoding. When using this encoding method to update particles, exchange sequence, crossover and mutation are used. These update operations can not only improve the convergence performance of the algorithm but also increase the diversity of the particle swarm. To test the effectiveness of exchange sequence under integer encoding, this study encodes the following two discrete algorithms for the single-objective (makespan) DFSP: 1) PSO-noES based on position-order real number encoding and without the use of exchange sequence. 2) PSO-ES for updating particles based on integer encoding exchange sequence. The two algorithms are compared on 12 standard problems, and the comparison of the objective value (makespan) is used to illustrate the effectiveness of the exchange sequence. Figure 13 presents the comparison results, which show that it is more efficient to use the PSO of the exchange sequence to solve the DFSP.

In addition, this study codes the following two discrete algorithms for multiobjective DFSP to test the effectiveness of exchange sequence: (i) MoPSO-noES based on real-order encoding and without the use of exchange sequence. (ii) MoPSO-ES for updating particles based on integer encoding exchange sequence. The two algorithms are compared on 12 standard problems, and the performance of the algorithm is shown by two comprehensive evaluation indexes, HV and IGD. Table 2 presents the comparison results, which show that: 1) in the given 12 problems, MoPSO-ES can get 9 better results; and 2) the significance analysis results show that MoPSODS-ES has significant convergence and distribution performance compared with MoPSO-noES. It can be seen that after using exchange sequence, particles can be updated more efficiently, thereby improving the convergence performance and distribution performance of the algorithm.

It can be seen from the above that the position-order real number encoding method is suitable for the original particle swarm to solve the continuous optimization problem, and the individual obtained by encoding represents the direction vector, and this encoding method cannot well record the sequence relationship of the two-vector encoding method. Therefore, it is not suitable for solving discrete optimization problems. The integer encoding method is suitable for discrete PSO to solve discrete optimization problems. This encoding method can use exchange sequence to update individuals. The flexible exchange and update of the encoded individual sequence can not only preserve the context of the sequence, so that the good job sorting and factory allocation remain unchanged, but also can move to the good individual more reasonably.

5.4.   Effectiveness of direction search strategy

The direction search strategy explores the particle swarm in multiple directions of the PF, which can enhance the strong convergence ability of the algorithm in different areas of the PF. Therefore, it is necessary to test the effect of the direction search strategy on the performance of the algorithm. In this study, MoPSO and MoPSODS are compared on 12 standard problems, and the results of HV and IGD are used as the criterion for evaluating the performance of the algorithm, so that the effectiveness of the direction search strategy can be verified. Table 3 shows the HV and IGD values of MoPSO and MoPSODS on 12 problems. It can be seen that the index values of MoPSODS on all problems are significantly better than MoPSO, which fully proves the effectiveness of the direction search strategy. The direction search strategy can not only improve the convergence performance of the algorithm, but also improve the distribution performance.

5.5.   Effectiveness of DE search strategy

The DE search strategy can strengthen the connection between particles in the particle swarm, and perform differential operations on the updated particle swarm to improve the quality of the particles and prevent falling into local optimum. Therefore, the effectiveness of DE search strategy needs to be tested. This study also compares MoPSODS and MoPSODS-DE on 12 standard problems, and uses the results of two frequently used comprehensive evaluation indexes HV and IGD as the criterion for evaluating whether the DE search strategy is effective. The mean HV and mean IGD values of MoPSODS and MoPSODS-DE after 30 runs on 12 problems are shown in Figure 4. The two indexes results in Figure 4 illustrate that MoPSODS-DE outperforms MoPSODS on 10 of the 12 problems, and the performance is significant on large-scale problems. Such results fully demonstrate the effectiveness of the DE search strategy. Therefore, it can be shown that the addition of DE search strategy further improves the convergence performance and distribution performance of the algorithm.

5.6.   Results and discussion

This section analyzes algorithm performance through various indexes. HV, GD, IGD and Spread mainly reflect the convergence performance and the distribution performance of the algorithm.

The HV index results are shown in Table 5. It can be seen that in 10 of the 12 problems, the HV value of the MoPSODS-DE is the largest, which shows that the HV value of other algorithms is significantly smaller than that of the MoPSODS-DE algorithm. In addition, in the two problems that the main algorithm does not perform best, the significant difference between the algorithm is not obvious. Therefore, in terms of HV value results, MoPSODS-DE is better than all comparison algorithms, so its convergence performance is the most outstanding.

The results of the GD index are shown in Table 6, and MoPSODS-DE performs the best on 10 problems. On the other 2 problems, although the GD value of MOHEA is the largest, it is not obvious compared with the significant difference analysis of MoPSODS-DE. Table 7 shows the results of the IGD index, where MoPSODS-DE performs the best on all problems. Therefore, it can be seen from the results of GD index and IGD index that MoPSODS-DE has the best convergence performance and distribution performance among all the compared algorithms.

From the results of Spread index in Table 8, MoPSODS-DE performs best on 6 of the 12 questions. The only problems where MoPSODS-DE performs significantly worse are small-scale problems such as 20_5_3 and 50_5_3. For other problems, the distribution performance of MoPSODS-DE is still very outstanding. In conclusion, MoPSODS-DE has good convergence performance, and MoPSODS-DE can also make the distribution of the obtained solutions more diverse and extensive.

Besides, in order to clearly analyze the significant difference between the MoPSODS-DE and the comparison algorithms on all 12 problems, the experiment conducted the multi-problem Wilcoxon signed rank test ^[68], and the experimental results are listed in Table 9. Further, we set the significance level $\alpha$ equal to 0.05. Among them, $"w/t/l"$ means that MoPSODS-DE is significantly better than, similar to and worse than comparison algorithms. $R+$ is the number of times the main algorithm was better than the comparison algorithms in each test, and $R-$ is the opposite. The p-value denotes the importance of the result of a test. The smaller the p-value, the greater the significance of the result.

From the experimental results in Table 9, it is clear in terms of the four indexes of HV, GD, IGD, and Spread, that all of MoPSODS-DE were significantly better than the compared algorithms, and they also had an absolute advantage in $R+$. The same result is also reflected in the significant difference $\alpha$. It can be seen from these results that MoPSODS-DE has achieved excellent overall performance, and in terms of distribution performance, MoPSODS-DE also performed very well.

The $C$ index is an index that measures the convergence of algorithm performance ^[20]. $C$($S_1$, $S_2$) is the percent of the individuals in $S_2$ who are weakly dominated by $S_1$. The larger $C$($S_1$, $S_2$) is, the better $S_1$ outperforms $S_2$ in $C$. In Table 10, the $C$ index results show that MoPSODS-DE is definitely better than other comparison algorithms in 11 problems. The $C$ index shows that MoPSODS-DE is well realized in terms of the convergence and distribution of the solution to the problem.

Finally, due to space limitations, only the convergence speed comparison of all algorithms for large-scale problems with more than 200 jobs is listed. Figure 14 shows the values of the comprehensive evaluation indexes HV and IGD obtained by the algorithm under different iterations. It can be seen from the Figure 14 that the MoPSODS-DE algorithm has the fastest convergence speed, and the convergence speed of MoPSO and MOEA/D is also a little faster than other algorithms. At the same time, it can be observed that the index values obtained by MoPSODS-DE are better. It fully illustrates the convergence performance and distribution performance of MoPSODS-DE due to all other compared algorithms.

In conclusion, from the indicators of HV, GD, and IGD, MoPSODS-DE has obvious advantages, indicating that MoPSODS-DE has the best convergence performance and distribution performance. From the Spread index, it can be seen that MoPSODS-DE outperforms other algorithms in large-scale problems, and the solutions generated by MoPSODS-DE are more diverse and widely distributed. Besides, MoPSODS-DE generates solutions that outperform other algorithms on all problems in terms of $C$ index. MoPSODS-DE uses the grouping strategy of enhanced direction search to make the particle swarm rapidly move to the PF in multiple directions, thereby accelerating the convergence of the particle swarm; at the same time, MoPSODS-DE uses the combination of exchange sequence, crossover and mutation to update the particles, which further improves the search ability of the algorithm and the diversity of the particle swarm. In addition, the algorithm uses the DE search strategy to perform local search on the particle swarm, which improves the quality of the solution and further enhances the convergence performance and distribution performance of the algorithm. Therefore, it can be proved that the MoPSODS algorithm combining the grouping operation, the exchange sequence operation, the crossover operation and the mutation operation performs well in the convergence performance and distribution performance, and with the addition of DE search, the MoPSODS-DE algorithm performs the best, which is superior to other compare algorithms.

6.   Conclusions

This study proposes the MoPSODS-DE to solve the DFSP. The objective is to minimize the makespan and total processing time. The MoPSODS-DE uses a combination of two strategies. One is the enhanced direction search strategy, which divides the particle swarm into different sub-particle swarms according to the upper area of the PF, the central area of the PF, and the lower area of the PF, so that individuals in the particle swarm can be in multiple directions on the PF at the same time. This strategy can specifically enhance the strong convergence ability of particle swarm optimization in different areas of PF, and improve the convergence of the algorithm. The second is the differential evolution search strategy, which can prevent MoPSODS from converging prematurely and avoid getting stuck in local optimum to find better solutions. Furthermore, an encoding-decoding method based on two-vector representation is designed and uses exchange sequences, crossover, and mutation to operate on particles reasonably efficiently. On 12 benchmark problems, MoPSODS-DE is compared with other evolutionary algorithms, and these comparison algorithms are optimized using DOE, and the best parameter combination of each algorithm is obtained. The evaluation indicators of the evaluation algorithm in this study are HV, GD, IGD, Spread, and $C$. The proposed MoPSODS-DE has good distribution performance, convergence performance, and can maintain the diversity of the particle swarm.

The MoPSODS-DE also has shortcoming, that is, the CPU running time is not the shortest. The main computational time overhead is due to the additional computational time required for DE search, but experimental analysis shows that DE search can improve the convergence performance and diversity performance of MoPSODS-DE. Therefore, how to reduce the computation time of the optimization algorithm is a problem worth studying in the future. In addition, other extended problems related to the DFSP, such as more practical problems such as green energy and job assembly, are also worth investigating.

Acknowledgments

This research work is supported by the National Natural Science Foundation of China (62172141, 61772173, U1904107), Innovative Funds Plan of Henan University of Technology (2020ZKCJ02), Science and Technology Research Project of Henan Province (212102210357), Innovative Research Team (in Science and Technology) in University of Henan Province (21IRTSTHN018), Zhengzhou Science and Technology Collaborative Innovation Project (21ZZXTCX19), and the Grant-in-Aid for Scientific Research (C) of Japan Society of Promotion of Science (19K12148).

Conflict of interest

The authors declare there is no conflict of interest.