# Multipopulation Genetic Algorithms with Different Interaction Structures to Solve Flexible Job-Shop Scheduling Problems: A Network Science Perspective

Authors: Ding-Shan Deng; Wei Long; Yan-Yan Li; Xiao-Qiu Shi
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Populations of multipopulation genetic algorithms (MPGAs) parallely evolve with some interaction mechanisms. Previous studies have shown that the interaction structures can impact on the performance of MPGAs to some extent. This paper introduces the concept of complex networks such as ring-shaped networks and small-world networks to study how interaction structures and their parameters influence the MPGAs, where subpopulations are regarded as nodes and their interaction or migration of elites between subpopulations as edges. After solving the flexible job-shop scheduling problem (FJSP) by MPGAs with different parameters of interaction structures, simulation results were measured by criteria, such as success rate and average optimal value. The analysis reveals that (1) the smaller the average path length (APL) of the network is, the higher the propagation rate will be; (2) the performance of MPGAs increased first and then decreased along with the decrease of APL, indicating that, for better performance, the networks should have a proper APL, which can be adjusted by changing the structural parameters of networks; and (3) because the edge number of small-world networks remains unchanged with different rewiring possibilities of edges, the change in performance indicates that the MPGA can be improved by a more proper interaction structure of subpopulations as other conditions remain unchanged.


---

## Body

## 1. Introduction

Genetic algorithm (GA) [1], an original metaheuristic, is easy to fall into local optima when employed to solve resource-constrained project scheduling problems [2] such as flexible job-shop schedule problems (FJSPs) due to the complexity of searched space and high dimensions [3].To maintain population diversity (enhancing the search diversity) and avoid premature convergence, a multipopulation genetic algorithm (MPGA) [4] is one feasible method where subpopulations are generated and individuals migrate periodically among them. In the evolutionary process, when intra-subpopulation evolution pushes individuals towards different local optima, migration can introduce new genes into the subpopulations [5]. As in [5, 6], the migration can be implemented in different topologies that define how subpopulations are connected, which can influence the outcomes of MPGAs to some extent.To more effectively research on the performance with different interaction structures between subpopulations, we can consider subpopulations as nodes and their interaction or migration of elites between subpopulations as edges so that MPGA can be regarded as complex networks. Similarly, but not identically, Du et al. [7] proposed the networked evolutionary algorithm where nodes represent information process units, i.e., individuals, and connections denote information transmission links. Payne et al. [8] shed light on dynamic population structures, wherein edges are dynamically rewired according to the genotypic or phenotypic properties of individuals or according to the success of prior interindividual interaction. Chao et al. [9, 10] proposed a multiobjective cellular grey wolf optimizer with a topological structure for hybrid flow-shop scheduling, in which each wolf or individual can be regarded as a grid of a lattice structure, and the interaction among them is restricted to the neighborhood. In the aforementioned three papers, the nodes of networks were regarded as individuals; however, this paper mainly discusses the interactions of subpopulations. For example, Leitao et al. [5] conducted simulations to study how the topology of island models (subpopulations) impacts the effectiveness and diversity of evolutionary algorithms with three types of networks, i.e., the ring and torus, as well as the fully connected networks. Besides, Zhang and Li [11] proposed a multiobjective evolutionary algorithm based on decomposition (MOEA/D). It can also be regarded as a multipopulation algorithm, and each subpopulation has only one individual. Fu et al. [12] extended and used the idea of MOEA/D to realistic problems, e.g., flow-shop scheduling problem, and each subpopulation contains multiple individuals. MOEA/D decomposes a multiobjective optimization problem into a number of scalar optimization subproblems, and the offspring individual to one subproblem can be generated upon the parent individuals from its neighboring subproblems. This neighborhood relation can also be treated as a ring-shaped network. Unlike MPGA discussed in this paper, the subpopulations of MOEA/D are heterogeneous, representing different subproblems, and cannot evolve independently to get an optimal solution. Additionally, other networks such as 4D hypercubes [13] and a 4×4 toroidal mesh are also utilized [4].When it comes to scheduling problems [14, 15], MPGAs are widely used to get the optimal solution. Kimms et al. [16] introduced an MPGA as a procedure to solve the synchronized and integrated two-level lot sizing and scheduling problem; in the migration process, a copy of each best individual is inserted into the next population replacing a random individual. Therefore, the interaction of subpopulations can be abstracted as a ring-shaped network. Zandieh and Karimi [17] proposed an MPGA to search the Pareto optimal solution for a multiobjective group scheduling problem in hybrid flexible flow-shop with sequence-dependent setup times, where the interaction of subpopulations can be abstracted as fully connected networks.We specifically choose the FJSP [18–21] as a measurement to test the performance of the MPGA as interaction structures of subpopulations vary. Research on the FJSP is essential for enterprises, especially for small-medium enterprises, to carry out production planning and scheduling in order to meet the delivery dates under a complex market [22]. For instance, Zhang et al. [23] put forward the MPGA based on the multiobjective scheduling of flexible job-shop, in which the abstracted topology of migration is a centralized network.As summarized in Table1, complex networks, such as scale-free networks and small-world networks, have been employed in cellular evolutionary algorithms to control interaction behaviors between individuals. However, for MPGAs, the focus of interaction structures among subpopulations is mainly rings, tori, hypercubes, and so on, and complex networks have rarely been utilized to represent subpopulations and their interaction purposefully.Table 1
Summary of background studies.
ReferencesNetworks or topologiesDenotation of nodes[4]Toroidal mesh networks; hypercubesSubpopulations[5]Fully connected networks; rings; toriSubpopulations[7]BA networksIndividuals[8]Scale-free networks; small-world networksIndividuals[9, 10]LatticesIndividuals[11, 12]RingsSubpopulations[13]4D hypercubesSubpopulations[16]Ring-shaped networksSubpopulations[17]Fully connected networksSubpopulations[23]Centralized networksSubpopulationsTherefore, in our previous studies [24, 25], we addressed how seven different network structures, including the ring-shaped network and the small-world network, influence the propagation rate of advantageous genes and thus affect the performance of MPGA for solving the FJSP. However, only the scale of networks is discussed, and the influence of other structural parameters is not researched. Thus, in this paper, we mainly take advantage of ring-shaped networks and small-world networks [26] to study the structural influence on the performance of MPGAs to solve the FJSP because ring-shaped networks are elementary topology widely used in MPGAs intentionally or unintentionally, and small-world networks with a definite node number will possess the same edge number as the rewiring possibility of edges varies, which would make them convenient and easy to study the influence of interaction structures as other conditions remain unchanged.Concretely, we change the inherent structural parameter of the corresponding network constantly and record the simulation results of MPGAs with the network. The success rate (SR) and average optimal value (AOV) are utilized to measure the results of MPGAs. Besides, the Hamming distance index (HDI) is introduced to evaluate the difference between elite individuals and characterize the propagation rate of advantageous genes, which provide an insight into how different interaction structures of subpopulations affect the performance of MPGAs.The sections below are organized as follows. In Section2, basic knowledge of complex networks and MPGA with networks for solving the FJSP, as well as the evaluation index, are introduced. Then, experiments and results analysis are reported in Section 3. Conclusions follow in Section 4.


## 2. Preliminaries

### 2.1. Network Models

A network consists of some nodes connected by some edges with a certain topology (structure) [27]. Some basic concepts utilized to characterize the network are introduced in the following. On top of that, typic models of the ring-shaped network and small-world network are elaborated.


#### 2.1.1. Degree and Average Path Length

The concept of degree is the most fundamental character and measure of a node in a network (in this paper, networks are undirected networks), and the degree of a node is defined to be the numberki of its existing edges. The average degree of a network is the average value of all such node degrees ki over the entire network, denoted by‾k.The average distance or average path length (APL) of a network is defined to be the average value of all distances over the network:(1)APL=2N×N−1∑i<jdij,where N is the total number of nodes in the network and dij is the distance between node i and node j, i.e., the total number of edges connecting them through the shortest route [27].


#### 2.1.2. Ring-Shaped Networks

A network in which every node has the same degree, i.e., the same number of connecting edges, is called a regular network. A typical sparse regular network is a nearest-neighbor coupled network, where every node is connected to 2K-nearest neighbors, K nodes on each side. Particularly, as shown in Figure 1, such a network with a periodic boundary connectivity condition is a ring-shaped network [27].Figure 1
Ring-shaped networks.


#### 2.1.3. Small-World Networks

The so-called WS small-world networks [28] were first proposed by Watts and Strogatz in 1998, which have both features of large clustering coefficients and short average path lengths. This model is generated by rewiring the edges of a ring-shaped network one by one with probability p, in which the case of p=0 corresponds to a regular network, and p=1 corresponds to a kind of ER random-graph networks.A WS small-world network can be generated by the following algorithm:(1)
Start from a ring-shaped network withN nodes, where each node is connected to its 2KK>0 neighbors, K nodes on each side(2)
For every pair of connected nodes in the ring-shaped network, reconnect the edge with possibilityp in such a way that the beginning end of the edge is kept, but the other end is disconnected, and then rewire it to a node randomlyThis operation is conducted edge by edge on the original ring-shaped network, once and once only, either clockwise or counterclockwise. Additionally, self-loop and multiple edges are avoided. Figure2 shows an example of WS small-world networks. Since the WS small-world network is the original and most widely used topology, the term “small-world network” mentioned in this paper is always the WS small-world network.Figure 2
WS small-world networks.Additionally, compared to WS small-world networks, the algorithm can be modified by replacing “random rewiring edges” with “random adding edges,” resulting in the NW small-world network model [29]. However, the edge number of NW small-world networks is changed with the rewiring process, so we discard the use of NW small-world networks.


### 2.2. Flexible Job-Shop Scheduling Problem

FJSP, which is also introduced in [24, 30], can be divided into two subproblems: machine subproblem and operation subproblem, namely, selecting a specific machine for each operation and arranging a proper processing order of all operations. Usually, there exist a job set within n jobs, labeled as J=Jii=1n, and a machine set within m machines, labeled as M=Mkk=1m. For each job Ji, there is an operation set within li operations, labeled as Oi=Oijj=1li. Each operation Oij can only be processed on a specific group of machines, labeled as Sij⊆M. The goal of the scheduling problem is to assign operations with proper orders to proper machines at the proper time, to pursue some given objectives such as minimizing maximum machine workload and minimizing makespan. For simplicity without losing generality, we choose to minimize makespan as the objective. The mathematical model is as follows [25]:(2)minFmax=maxFij,∀i,j,s.t.(3)Fijk≥0,∀i,j,k,(4)Pijk≥0,∀i,j,k,(5)Bijk≥0,∀i,j,k,(6)Xijk∈0,1,∀i,j,k,(7)i,j,l,i′,j′∈1,2,3,…,,(8)∑k∈SijXijk=1∧Fijk−Bijk=Pijk,∀i,j,(9)Fij−Fij−1≥Pijk⋅Xijk,∀i,j,k,(10)Fi′j′k≤Bijk∨Fijk≤Bi′j′k,∀i′,j′≠i,j.Table2 lists the notations used in the FJSP. For equation (6), as Oij is arranged to be processed on Mk, Xijk=1; otherwise, 0. Equation (7) denotes the domains of variables. Equation (8) guarantees that each operation can only be processed at one machine without disruption, and ‘∧’ represents logical AND. Therefore, Fij=Fijk,∀i,j,∃k,k∈Sij. Equation (9) ensures each job will be processed in the correct order. Equation (10) ensures that one machine can only process one operation at the same time, and ‘∨’ represents logical OR.Table 2
Notations related to the FJSP.
SymbolMeaningnThe number of jobsJThe job setJiTheith jobmThe number of machinesMThe machine setMkThekth machineliThe number of operations of theith jobOiThe operation set of theith jobOijThejth operation of the ith jobSijThe candidate machine setOij can be processed onFmaxMakespanFijCompletion time ofOijFijkFinish time ofOij on MkPijkProcessing time ofOij on MkBijkStart time ofOij on MkXijkIndicate ifOij is processed on MkBesides, Table3 lists a four-step process of the FJSP for 2 jobs with a total of 4 operations 4×4, where ‘—’ means the corresponding operation cannot be processed in that machine. The fourth number “19” in the first line, for example, represents that O11 can be processed in M4 within time horizon 19.Table 3
One instance of the FJSP.
JobsOperationsMachinesM1M2M3M4J1O1124—1619O12—21—13J2O219—68O22——86


### 2.3. MPGA with Networks for Solving the FJSP

For using an MPGA with networks to solve the FJSP, multiple operations are included, such as encoding, decoding, crossing, mutation, and migration.


#### 2.3.1. Encoding

In this algorithm, we employ an integer encoding method [31] to generate individuals, each of which can represent a feasible solution of the FJSP. The process of encoding can be divided into two stages, to wit, encoding for machine subproblem and operation subproblem. For the former, an encoded solution is denoted by a string of integers whose length is equal to the total operation number of all jobs, in which an integer in each position represents an operation, and the value represents the ordinal number of machines in the candidate machine set Sij of this operation. For the FJSP mentioned in Table 3, suppose there is an encoded solution (1 2 3 1) within four integers which reflect that there are a total of four operations in the two jobs. The first integer 1 means that O11 is arranged with M1. For the second integer, because the candidate machine set of O12 is S12=M2,M4, integer 2 denotes that O12 will be processed on the fourth machine, namely, M4, rather than M2. Similarly, O21 is arranged to be processed on M4 and O22 on M3.For the latter one, a string of integer represents an encoded solution, and the number of integers equals the total operation number of all jobs too. The value of an integer denotes the serial number of jobs, and if a job hasli operations, the job number will appear li times. The sequence of integers represents the processing order of the corresponding jobs. For example, there is an encoded solution (2 1 1 2) for the FJSP as shown in Table 3, in which integers “1” and “2” appear twice, respectively, due to both J1 and J2 including two operations. In detail, the fourth integer is “2,” and it is the second occurrence of “2” in this string, indicating the processing of the second operation of the second job, namely, O22. Therefore, the processing order of operations indicated by this solution is O21, O11, O12, and O22.Therefore, we combine (1 2 3 1) and (2 1 1 2) together, getting an individual (1 2 3 1 2 1 1 2), which meansO21 will be processed on M4, O11 on M1, O12 on M4, and O22 on M3 sequentially.


#### 2.3.2. Decoding

An encoded individual must be decoded into an original solution of the FJSP to calculate the fitness. The decoding algorithm proposed by Shi et al. [30] was adopted in this paper as follows:At first, we convert an individual into a matrix, labeled asHd. For instance, [3 1 2 1 2 1 1 2], a given individual for the FJSP in Table 3, can be converted into a matrix as shown in Figure 3.In matrixHd, the first column stores the sequence numbers of jobs; the second column stores the sequence numbers of operations of the corresponding job; the third column stores the sequence numbers of machines in Sij; and the fourth column stores the time horizon as the corresponding operation processed in the specific machine. Besides, the last two columns are reserved to store the start time and completion time, respectively. For instance, the first row of Hd represents that O21 is processed on M3 within time horizon 19.Next, we update the start time and completion time row by row. IfOij represented by the corresponding row is the first operation of Ji and Mk selected to process Oij has not processed any other operation, then Bijk=0,andFijk=Bijk+Pijk. Assign Bijk and Fijk into columns 5 and 6 of Hd.IfOij is the first operation of Ji and Mk has processed other operations, then find all idle-time intervals of Mk denoted by SqEqq=1,2,…,, check all idle intervals one by one to find the first one whose duration is longer than Pijk, and set Bijk=SqandFijk=Bijk+Pijk.IfOij is not the first operation of Ji, the former operation Oij−1 is finished on Mk′, and Mk has not been assigned to any operation, then Bijk=Fij−1k′andFijk=Bijk+Pijk.IfOij is not the first operation of Ji and Mk has been assigned to other operations, then search for and check all idle-time intervals of Mk from left to right on the timeline. If Eq−Sq≥PijkandFij−1k′≤Sq, then Bijk=SqandFijk=Bijk+Pijk; if Eq−Sq≥Pijk,Fij−1k′≥Sq,andEq−Fij−1k≥Pijk, then Bijk=Fij−1k′andFijk=Bijk+Pijk. Since Eq of the last interval is positive infinity, an idle-time interval can always be found to satisfy one of these conditions. In this way, there is no earlier idle-time interval on Mk for Ji to insert without delaying or interrupting other operations, indicating that an active schedule is achieved [32]. The optimal schedule for any regular measure of performance, including makespan, will be a member of the active schedule set [33]. Therefore, the optimal value of makespan can be found by MPGAs with this decoding method.At last, after updating columns 5 and 6 of all rows, the maximum value of column 6 can be found as the makespan of the FJSP, and our goal is to minimize makespan through MPGA. As shown in Figure3, the makespan of our example is 35.Figure 3
Illustration of decoding.


#### 2.3.3. Crossing

The processing of crossing is also divided into two stages. For the machine subproblem, we implement the standard two-point crossing process [30]. In this process, we choose two individuals randomly, called parent 1 and parent 2, both of which will be divided into three parts with the same size at random common breakpoints.We partially select parts from parent 1 and parent 2 to shape offspring 1, and remaining parts compose offspring 2, which are indicated in Figure4.Figure 4
Crossing of the machine subproblem.For the code of the operation subproblem [34], all jobs are randomly divided into two groups: group 1 and group 2. Offspring 1 inherits the integers of parent 1 belonging to group 1 and those of parent 2 belonging to group 2; in the same method, offspring 2 inherits the integers of parent 1 belonging to group 2 and those of parent 2 belonging to group 1, respectively, meanwhile preserving the sequence of these integers. For example, as indicated in Figure 5, jobs 1 and 2 belong to group 1, and jobs 3 and 4 belong to group 2. To get offspring 1, we first preserve the integers belonging to group 1, i.e., job 1 and job 3. Next, we insert integers of group 2, i.e., job 2 and job 4, from parent 2 into parent 1. Meanwhile, the original sequence of integers representing jobs 2 and 4 is also preserved. In the same way, we get offspring 2.Figure 5
Crossing of the operation subproblem.


#### 2.3.4. Mutation

For MPGAs, mutation can be a supplementary strategy to maintain diversity [35]. Based on the characteristic of integer encoding, it is divided into two stages. In the process of the machine subproblem, we randomly choose several individuals based on the mutation probability. Then, in several random positions of these selected individuals, the original integers can be replaced by alternative integers, which should be smaller than the total number of the corresponding candidate machines. For the operation subproblem, in the same way, some individuals are randomly selected. Within each selected individual, several pairs of integers are randomly selected, and the positions of two integers, within each pair, are swapped.


#### 2.3.5. Migration

The migration strategy is often utilized to mitigate the premature convergence of evolutionary algorithms [36]. To obtain an MPGA with networks, we divide the population of standard genetic algorithms into subpopulations, denoted by nodes. Communications between nodes, denoted by edges, occur when certain individuals in one node migrate to another periodically based on different networks. Figure 6 displays an MPGA with a certain network.Figure 6
Schematic diagram of the MPGA with networks.In the process of migration, we select a subpopulation randomly in an MPGA, denoted by one node of the network, and find all neighbor nodes (subpopulations). Next, we find the best elite of these subpopulations which will replace a random individual in each of these selected subpopulations.


#### 2.3.6. Implementation of the MPGA with Networks to Solve the FJSP

As shown in Algorithm1, we initialize Ns subpopulations, and each of them can be run parallel [37]. Within each subpopulation, there are S individuals, namely, Ns×S individuals, in total, and we first decode the individuals to get fitness, i.e., makespan. Then, tournament selection strategy is employed to generate the next subpopulation and maintain the subpopulation size unchanged. After that, crossing and mutation operators are implemented to generate new individuals. After all of the subpopulations are updated, we calculate the fitness again and select the elites of each subpopulation into the elite set. At last, our proposed migration strategy is utilized to propagate advantageous genes. As iterating to the maximum generation, the best fitness and individual are output.Algorithm 1: MPGA with networks.
Input:FJSP: the FJSP to solveNs: the number of subpopulations;S: the size of each subpopulation;Pm: mutation probability;itermax: the maximum generation to run;Output:R∗: the best fitness or makespan;g∗: the best individual;G=G1,G2,…,GM ← encodingFJSP,Ns,S;Whilegenration≤itermax doFor i_sub = 1:Ns doR ← decodingGi_sub;Gsi_sub ← selectionGi_sub based on R;Gci_sub ← crossing Gsi_sub;Gi_sub ← mutation Gci_sub,Pm;End forR ← decodingG;Select the elites of each subpopulation to compose elite setE;G ← migrationG,E;End whileFind the best individualg∗ from E and the corresponding fitness R∗.


### 2.4. Evaluation Indexes

An index to evaluate the propagation rate of advantageous genes among subpopulations under different network structures is needed. All the best individuals of each subpopulation are selected into an elite set. With the accumulation of advantageous genes through the operators of MPGA, the difference between individuals becomes small. Therefore, if the propagation rate of advantageous genes is larger, the difference between these elites will be smaller. Thus, we introduce the HDI to evaluate the difference between elite individuals. The HDI is calculated as follows [24]:(11)HDI=∑i=1100∑j=12×JT1−δhij1,hij2100×2×JT,where we only sample 100 pairs of individuals to limit time consumption and δ… is the Kronecker function. As the two independent variables are unequal, the value of δ… is 0; otherwise, 1. JT represents the total operation number of all jobs, and hij1andhij2 denote the ith pair’s two jth integers in this sample. Therefore, smaller HDI indicates less difference between elites, suggesting the propagation among subpopulations is faster.Besides, we use the SR to assess the performance of MPGA, better performance with higher SR. This measure is defined as follows:(12)SR=NsNt,where Ns denotes the number of times the algorithm finds the optimal value and Nt denotes the total number of times it runs. The optimal value, mentioned in this paper, is defined as the best value found by all current algorithms by now, rather than the best value found by the currently running algorithm.Also, due to the stochastic nature of MPGA, we have to run the proposed algorithm repeatedly, and statistical features are often utilized to analyze simulation results [38–40]. Therefore, the AOV is introduced and calculated as follows:(13)AOV=1Nt∑iNtOVi,in which OVi is the optimal value obtained by the ith run of the algorithm. The smaller the AOV is, the better the performance will be.


## 2.1. Network Models

A network consists of some nodes connected by some edges with a certain topology (structure) [27]. Some basic concepts utilized to characterize the network are introduced in the following. On top of that, typic models of the ring-shaped network and small-world network are elaborated.


### 2.1.1. Degree and Average Path Length

The concept of degree is the most fundamental character and measure of a node in a network (in this paper, networks are undirected networks), and the degree of a node is defined to be the numberki of its existing edges. The average degree of a network is the average value of all such node degrees ki over the entire network, denoted by‾k.The average distance or average path length (APL) of a network is defined to be the average value of all distances over the network:(1)APL=2N×N−1∑i<jdij,where N is the total number of nodes in the network and dij is the distance between node i and node j, i.e., the total number of edges connecting them through the shortest route [27].


### 2.1.2. Ring-Shaped Networks

A network in which every node has the same degree, i.e., the same number of connecting edges, is called a regular network. A typical sparse regular network is a nearest-neighbor coupled network, where every node is connected to 2K-nearest neighbors, K nodes on each side. Particularly, as shown in Figure 1, such a network with a periodic boundary connectivity condition is a ring-shaped network [27].Figure 1
Ring-shaped networks.


### 2.1.3. Small-World Networks

The so-called WS small-world networks [28] were first proposed by Watts and Strogatz in 1998, which have both features of large clustering coefficients and short average path lengths. This model is generated by rewiring the edges of a ring-shaped network one by one with probability p, in which the case of p=0 corresponds to a regular network, and p=1 corresponds to a kind of ER random-graph networks.A WS small-world network can be generated by the following algorithm:(1)
Start from a ring-shaped network withN nodes, where each node is connected to its 2KK>0 neighbors, K nodes on each side(2)
For every pair of connected nodes in the ring-shaped network, reconnect the edge with possibilityp in such a way that the beginning end of the edge is kept, but the other end is disconnected, and then rewire it to a node randomlyThis operation is conducted edge by edge on the original ring-shaped network, once and once only, either clockwise or counterclockwise. Additionally, self-loop and multiple edges are avoided. Figure2 shows an example of WS small-world networks. Since the WS small-world network is the original and most widely used topology, the term “small-world network” mentioned in this paper is always the WS small-world network.Figure 2
WS small-world networks.Additionally, compared to WS small-world networks, the algorithm can be modified by replacing “random rewiring edges” with “random adding edges,” resulting in the NW small-world network model [29]. However, the edge number of NW small-world networks is changed with the rewiring process, so we discard the use of NW small-world networks.


## 2.1.1. Degree and Average Path Length

The concept of degree is the most fundamental character and measure of a node in a network (in this paper, networks are undirected networks), and the degree of a node is defined to be the numberki of its existing edges. The average degree of a network is the average value of all such node degrees ki over the entire network, denoted by‾k.The average distance or average path length (APL) of a network is defined to be the average value of all distances over the network:(1)APL=2N×N−1∑i<jdij,where N is the total number of nodes in the network and dij is the distance between node i and node j, i.e., the total number of edges connecting them through the shortest route [27].


## 2.1.2. Ring-Shaped Networks

A network in which every node has the same degree, i.e., the same number of connecting edges, is called a regular network. A typical sparse regular network is a nearest-neighbor coupled network, where every node is connected to 2K-nearest neighbors, K nodes on each side. Particularly, as shown in Figure 1, such a network with a periodic boundary connectivity condition is a ring-shaped network [27].Figure 1
Ring-shaped networks.


## 2.1.3. Small-World Networks

The so-called WS small-world networks [28] were first proposed by Watts and Strogatz in 1998, which have both features of large clustering coefficients and short average path lengths. This model is generated by rewiring the edges of a ring-shaped network one by one with probability p, in which the case of p=0 corresponds to a regular network, and p=1 corresponds to a kind of ER random-graph networks.A WS small-world network can be generated by the following algorithm:(1)
Start from a ring-shaped network withN nodes, where each node is connected to its 2KK>0 neighbors, K nodes on each side(2)
For every pair of connected nodes in the ring-shaped network, reconnect the edge with possibilityp in such a way that the beginning end of the edge is kept, but the other end is disconnected, and then rewire it to a node randomlyThis operation is conducted edge by edge on the original ring-shaped network, once and once only, either clockwise or counterclockwise. Additionally, self-loop and multiple edges are avoided. Figure2 shows an example of WS small-world networks. Since the WS small-world network is the original and most widely used topology, the term “small-world network” mentioned in this paper is always the WS small-world network.Figure 2
WS small-world networks.Additionally, compared to WS small-world networks, the algorithm can be modified by replacing “random rewiring edges” with “random adding edges,” resulting in the NW small-world network model [29]. However, the edge number of NW small-world networks is changed with the rewiring process, so we discard the use of NW small-world networks.


## 2.2. Flexible Job-Shop Scheduling Problem

FJSP, which is also introduced in [24, 30], can be divided into two subproblems: machine subproblem and operation subproblem, namely, selecting a specific machine for each operation and arranging a proper processing order of all operations. Usually, there exist a job set within n jobs, labeled as J=Jii=1n, and a machine set within m machines, labeled as M=Mkk=1m. For each job Ji, there is an operation set within li operations, labeled as Oi=Oijj=1li. Each operation Oij can only be processed on a specific group of machines, labeled as Sij⊆M. The goal of the scheduling problem is to assign operations with proper orders to proper machines at the proper time, to pursue some given objectives such as minimizing maximum machine workload and minimizing makespan. For simplicity without losing generality, we choose to minimize makespan as the objective. The mathematical model is as follows [25]:(2)minFmax=maxFij,∀i,j,s.t.(3)Fijk≥0,∀i,j,k,(4)Pijk≥0,∀i,j,k,(5)Bijk≥0,∀i,j,k,(6)Xijk∈0,1,∀i,j,k,(7)i,j,l,i′,j′∈1,2,3,…,,(8)∑k∈SijXijk=1∧Fijk−Bijk=Pijk,∀i,j,(9)Fij−Fij−1≥Pijk⋅Xijk,∀i,j,k,(10)Fi′j′k≤Bijk∨Fijk≤Bi′j′k,∀i′,j′≠i,j.Table2 lists the notations used in the FJSP. For equation (6), as Oij is arranged to be processed on Mk, Xijk=1; otherwise, 0. Equation (7) denotes the domains of variables. Equation (8) guarantees that each operation can only be processed at one machine without disruption, and ‘∧’ represents logical AND. Therefore, Fij=Fijk,∀i,j,∃k,k∈Sij. Equation (9) ensures each job will be processed in the correct order. Equation (10) ensures that one machine can only process one operation at the same time, and ‘∨’ represents logical OR.Table 2
Notations related to the FJSP.
SymbolMeaningnThe number of jobsJThe job setJiTheith jobmThe number of machinesMThe machine setMkThekth machineliThe number of operations of theith jobOiThe operation set of theith jobOijThejth operation of the ith jobSijThe candidate machine setOij can be processed onFmaxMakespanFijCompletion time ofOijFijkFinish time ofOij on MkPijkProcessing time ofOij on MkBijkStart time ofOij on MkXijkIndicate ifOij is processed on MkBesides, Table3 lists a four-step process of the FJSP for 2 jobs with a total of 4 operations 4×4, where ‘—’ means the corresponding operation cannot be processed in that machine. The fourth number “19” in the first line, for example, represents that O11 can be processed in M4 within time horizon 19.Table 3
One instance of the FJSP.
JobsOperationsMachinesM1M2M3M4J1O1124—1619O12—21—13J2O219—68O22——86


## 2.3. MPGA with Networks for Solving the FJSP

For using an MPGA with networks to solve the FJSP, multiple operations are included, such as encoding, decoding, crossing, mutation, and migration.


### 2.3.1. Encoding

In this algorithm, we employ an integer encoding method [31] to generate individuals, each of which can represent a feasible solution of the FJSP. The process of encoding can be divided into two stages, to wit, encoding for machine subproblem and operation subproblem. For the former, an encoded solution is denoted by a string of integers whose length is equal to the total operation number of all jobs, in which an integer in each position represents an operation, and the value represents the ordinal number of machines in the candidate machine set Sij of this operation. For the FJSP mentioned in Table 3, suppose there is an encoded solution (1 2 3 1) within four integers which reflect that there are a total of four operations in the two jobs. The first integer 1 means that O11 is arranged with M1. For the second integer, because the candidate machine set of O12 is S12=M2,M4, integer 2 denotes that O12 will be processed on the fourth machine, namely, M4, rather than M2. Similarly, O21 is arranged to be processed on M4 and O22 on M3.For the latter one, a string of integer represents an encoded solution, and the number of integers equals the total operation number of all jobs too. The value of an integer denotes the serial number of jobs, and if a job hasli operations, the job number will appear li times. The sequence of integers represents the processing order of the corresponding jobs. For example, there is an encoded solution (2 1 1 2) for the FJSP as shown in Table 3, in which integers “1” and “2” appear twice, respectively, due to both J1 and J2 including two operations. In detail, the fourth integer is “2,” and it is the second occurrence of “2” in this string, indicating the processing of the second operation of the second job, namely, O22. Therefore, the processing order of operations indicated by this solution is O21, O11, O12, and O22.Therefore, we combine (1 2 3 1) and (2 1 1 2) together, getting an individual (1 2 3 1 2 1 1 2), which meansO21 will be processed on M4, O11 on M1, O12 on M4, and O22 on M3 sequentially.


### 2.3.2. Decoding

An encoded individual must be decoded into an original solution of the FJSP to calculate the fitness. The decoding algorithm proposed by Shi et al. [30] was adopted in this paper as follows:At first, we convert an individual into a matrix, labeled asHd. For instance, [3 1 2 1 2 1 1 2], a given individual for the FJSP in Table 3, can be converted into a matrix as shown in Figure 3.In matrixHd, the first column stores the sequence numbers of jobs; the second column stores the sequence numbers of operations of the corresponding job; the third column stores the sequence numbers of machines in Sij; and the fourth column stores the time horizon as the corresponding operation processed in the specific machine. Besides, the last two columns are reserved to store the start time and completion time, respectively. For instance, the first row of Hd represents that O21 is processed on M3 within time horizon 19.Next, we update the start time and completion time row by row. IfOij represented by the corresponding row is the first operation of Ji and Mk selected to process Oij has not processed any other operation, then Bijk=0,andFijk=Bijk+Pijk. Assign Bijk and Fijk into columns 5 and 6 of Hd.IfOij is the first operation of Ji and Mk has processed other operations, then find all idle-time intervals of Mk denoted by SqEqq=1,2,…,, check all idle intervals one by one to find the first one whose duration is longer than Pijk, and set Bijk=SqandFijk=Bijk+Pijk.IfOij is not the first operation of Ji, the former operation Oij−1 is finished on Mk′, and Mk has not been assigned to any operation, then Bijk=Fij−1k′andFijk=Bijk+Pijk.IfOij is not the first operation of Ji and Mk has been assigned to other operations, then search for and check all idle-time intervals of Mk from left to right on the timeline. If Eq−Sq≥PijkandFij−1k′≤Sq, then Bijk=SqandFijk=Bijk+Pijk; if Eq−Sq≥Pijk,Fij−1k′≥Sq,andEq−Fij−1k≥Pijk, then Bijk=Fij−1k′andFijk=Bijk+Pijk. Since Eq of the last interval is positive infinity, an idle-time interval can always be found to satisfy one of these conditions. In this way, there is no earlier idle-time interval on Mk for Ji to insert without delaying or interrupting other operations, indicating that an active schedule is achieved [32]. The optimal schedule for any regular measure of performance, including makespan, will be a member of the active schedule set [33]. Therefore, the optimal value of makespan can be found by MPGAs with this decoding method.At last, after updating columns 5 and 6 of all rows, the maximum value of column 6 can be found as the makespan of the FJSP, and our goal is to minimize makespan through MPGA. As shown in Figure3, the makespan of our example is 35.Figure 3
Illustration of decoding.


### 2.3.3. Crossing

The processing of crossing is also divided into two stages. For the machine subproblem, we implement the standard two-point crossing process [30]. In this process, we choose two individuals randomly, called parent 1 and parent 2, both of which will be divided into three parts with the same size at random common breakpoints.We partially select parts from parent 1 and parent 2 to shape offspring 1, and remaining parts compose offspring 2, which are indicated in Figure4.Figure 4
Crossing of the machine subproblem.For the code of the operation subproblem [34], all jobs are randomly divided into two groups: group 1 and group 2. Offspring 1 inherits the integers of parent 1 belonging to group 1 and those of parent 2 belonging to group 2; in the same method, offspring 2 inherits the integers of parent 1 belonging to group 2 and those of parent 2 belonging to group 1, respectively, meanwhile preserving the sequence of these integers. For example, as indicated in Figure 5, jobs 1 and 2 belong to group 1, and jobs 3 and 4 belong to group 2. To get offspring 1, we first preserve the integers belonging to group 1, i.e., job 1 and job 3. Next, we insert integers of group 2, i.e., job 2 and job 4, from parent 2 into parent 1. Meanwhile, the original sequence of integers representing jobs 2 and 4 is also preserved. In the same way, we get offspring 2.Figure 5
Crossing of the operation subproblem.


### 2.3.4. Mutation

For MPGAs, mutation can be a supplementary strategy to maintain diversity [35]. Based on the characteristic of integer encoding, it is divided into two stages. In the process of the machine subproblem, we randomly choose several individuals based on the mutation probability. Then, in several random positions of these selected individuals, the original integers can be replaced by alternative integers, which should be smaller than the total number of the corresponding candidate machines. For the operation subproblem, in the same way, some individuals are randomly selected. Within each selected individual, several pairs of integers are randomly selected, and the positions of two integers, within each pair, are swapped.


### 2.3.5. Migration

The migration strategy is often utilized to mitigate the premature convergence of evolutionary algorithms [36]. To obtain an MPGA with networks, we divide the population of standard genetic algorithms into subpopulations, denoted by nodes. Communications between nodes, denoted by edges, occur when certain individuals in one node migrate to another periodically based on different networks. Figure 6 displays an MPGA with a certain network.Figure 6
Schematic diagram of the MPGA with networks.In the process of migration, we select a subpopulation randomly in an MPGA, denoted by one node of the network, and find all neighbor nodes (subpopulations). Next, we find the best elite of these subpopulations which will replace a random individual in each of these selected subpopulations.


### 2.3.6. Implementation of the MPGA with Networks to Solve the FJSP

As shown in Algorithm1, we initialize Ns subpopulations, and each of them can be run parallel [37]. Within each subpopulation, there are S individuals, namely, Ns×S individuals, in total, and we first decode the individuals to get fitness, i.e., makespan. Then, tournament selection strategy is employed to generate the next subpopulation and maintain the subpopulation size unchanged. After that, crossing and mutation operators are implemented to generate new individuals. After all of the subpopulations are updated, we calculate the fitness again and select the elites of each subpopulation into the elite set. At last, our proposed migration strategy is utilized to propagate advantageous genes. As iterating to the maximum generation, the best fitness and individual are output.Algorithm 1: MPGA with networks.
Input:FJSP: the FJSP to solveNs: the number of subpopulations;S: the size of each subpopulation;Pm: mutation probability;itermax: the maximum generation to run;Output:R∗: the best fitness or makespan;g∗: the best individual;G=G1,G2,…,GM ← encodingFJSP,Ns,S;Whilegenration≤itermax doFor i_sub = 1:Ns doR ← decodingGi_sub;Gsi_sub ← selectionGi_sub based on R;Gci_sub ← crossing Gsi_sub;Gi_sub ← mutation Gci_sub,Pm;End forR ← decodingG;Select the elites of each subpopulation to compose elite setE;G ← migrationG,E;End whileFind the best individualg∗ from E and the corresponding fitness R∗.


## 2.3.1. Encoding

In this algorithm, we employ an integer encoding method [31] to generate individuals, each of which can represent a feasible solution of the FJSP. The process of encoding can be divided into two stages, to wit, encoding for machine subproblem and operation subproblem. For the former, an encoded solution is denoted by a string of integers whose length is equal to the total operation number of all jobs, in which an integer in each position represents an operation, and the value represents the ordinal number of machines in the candidate machine set Sij of this operation. For the FJSP mentioned in Table 3, suppose there is an encoded solution (1 2 3 1) within four integers which reflect that there are a total of four operations in the two jobs. The first integer 1 means that O11 is arranged with M1. For the second integer, because the candidate machine set of O12 is S12=M2,M4, integer 2 denotes that O12 will be processed on the fourth machine, namely, M4, rather than M2. Similarly, O21 is arranged to be processed on M4 and O22 on M3.For the latter one, a string of integer represents an encoded solution, and the number of integers equals the total operation number of all jobs too. The value of an integer denotes the serial number of jobs, and if a job hasli operations, the job number will appear li times. The sequence of integers represents the processing order of the corresponding jobs. For example, there is an encoded solution (2 1 1 2) for the FJSP as shown in Table 3, in which integers “1” and “2” appear twice, respectively, due to both J1 and J2 including two operations. In detail, the fourth integer is “2,” and it is the second occurrence of “2” in this string, indicating the processing of the second operation of the second job, namely, O22. Therefore, the processing order of operations indicated by this solution is O21, O11, O12, and O22.Therefore, we combine (1 2 3 1) and (2 1 1 2) together, getting an individual (1 2 3 1 2 1 1 2), which meansO21 will be processed on M4, O11 on M1, O12 on M4, and O22 on M3 sequentially.


## 2.3.2. Decoding

An encoded individual must be decoded into an original solution of the FJSP to calculate the fitness. The decoding algorithm proposed by Shi et al. [30] was adopted in this paper as follows:At first, we convert an individual into a matrix, labeled asHd. For instance, [3 1 2 1 2 1 1 2], a given individual for the FJSP in Table 3, can be converted into a matrix as shown in Figure 3.In matrixHd, the first column stores the sequence numbers of jobs; the second column stores the sequence numbers of operations of the corresponding job; the third column stores the sequence numbers of machines in Sij; and the fourth column stores the time horizon as the corresponding operation processed in the specific machine. Besides, the last two columns are reserved to store the start time and completion time, respectively. For instance, the first row of Hd represents that O21 is processed on M3 within time horizon 19.Next, we update the start time and completion time row by row. IfOij represented by the corresponding row is the first operation of Ji and Mk selected to process Oij has not processed any other operation, then Bijk=0,andFijk=Bijk+Pijk. Assign Bijk and Fijk into columns 5 and 6 of Hd.IfOij is the first operation of Ji and Mk has processed other operations, then find all idle-time intervals of Mk denoted by SqEqq=1,2,…,, check all idle intervals one by one to find the first one whose duration is longer than Pijk, and set Bijk=SqandFijk=Bijk+Pijk.IfOij is not the first operation of Ji, the former operation Oij−1 is finished on Mk′, and Mk has not been assigned to any operation, then Bijk=Fij−1k′andFijk=Bijk+Pijk.IfOij is not the first operation of Ji and Mk has been assigned to other operations, then search for and check all idle-time intervals of Mk from left to right on the timeline. If Eq−Sq≥PijkandFij−1k′≤Sq, then Bijk=SqandFijk=Bijk+Pijk; if Eq−Sq≥Pijk,Fij−1k′≥Sq,andEq−Fij−1k≥Pijk, then Bijk=Fij−1k′andFijk=Bijk+Pijk. Since Eq of the last interval is positive infinity, an idle-time interval can always be found to satisfy one of these conditions. In this way, there is no earlier idle-time interval on Mk for Ji to insert without delaying or interrupting other operations, indicating that an active schedule is achieved [32]. The optimal schedule for any regular measure of performance, including makespan, will be a member of the active schedule set [33]. Therefore, the optimal value of makespan can be found by MPGAs with this decoding method.At last, after updating columns 5 and 6 of all rows, the maximum value of column 6 can be found as the makespan of the FJSP, and our goal is to minimize makespan through MPGA. As shown in Figure3, the makespan of our example is 35.Figure 3
Illustration of decoding.


## 2.3.3. Crossing

The processing of crossing is also divided into two stages. For the machine subproblem, we implement the standard two-point crossing process [30]. In this process, we choose two individuals randomly, called parent 1 and parent 2, both of which will be divided into three parts with the same size at random common breakpoints.We partially select parts from parent 1 and parent 2 to shape offspring 1, and remaining parts compose offspring 2, which are indicated in Figure4.Figure 4
Crossing of the machine subproblem.For the code of the operation subproblem [34], all jobs are randomly divided into two groups: group 1 and group 2. Offspring 1 inherits the integers of parent 1 belonging to group 1 and those of parent 2 belonging to group 2; in the same method, offspring 2 inherits the integers of parent 1 belonging to group 2 and those of parent 2 belonging to group 1, respectively, meanwhile preserving the sequence of these integers. For example, as indicated in Figure 5, jobs 1 and 2 belong to group 1, and jobs 3 and 4 belong to group 2. To get offspring 1, we first preserve the integers belonging to group 1, i.e., job 1 and job 3. Next, we insert integers of group 2, i.e., job 2 and job 4, from parent 2 into parent 1. Meanwhile, the original sequence of integers representing jobs 2 and 4 is also preserved. In the same way, we get offspring 2.Figure 5
Crossing of the operation subproblem.


## 2.3.4. Mutation

For MPGAs, mutation can be a supplementary strategy to maintain diversity [35]. Based on the characteristic of integer encoding, it is divided into two stages. In the process of the machine subproblem, we randomly choose several individuals based on the mutation probability. Then, in several random positions of these selected individuals, the original integers can be replaced by alternative integers, which should be smaller than the total number of the corresponding candidate machines. For the operation subproblem, in the same way, some individuals are randomly selected. Within each selected individual, several pairs of integers are randomly selected, and the positions of two integers, within each pair, are swapped.


## 2.3.5. Migration

The migration strategy is often utilized to mitigate the premature convergence of evolutionary algorithms [36]. To obtain an MPGA with networks, we divide the population of standard genetic algorithms into subpopulations, denoted by nodes. Communications between nodes, denoted by edges, occur when certain individuals in one node migrate to another periodically based on different networks. Figure 6 displays an MPGA with a certain network.Figure 6
Schematic diagram of the MPGA with networks.In the process of migration, we select a subpopulation randomly in an MPGA, denoted by one node of the network, and find all neighbor nodes (subpopulations). Next, we find the best elite of these subpopulations which will replace a random individual in each of these selected subpopulations.


## 2.3.6. Implementation of the MPGA with Networks to Solve the FJSP

As shown in Algorithm1, we initialize Ns subpopulations, and each of them can be run parallel [37]. Within each subpopulation, there are S individuals, namely, Ns×S individuals, in total, and we first decode the individuals to get fitness, i.e., makespan. Then, tournament selection strategy is employed to generate the next subpopulation and maintain the subpopulation size unchanged. After that, crossing and mutation operators are implemented to generate new individuals. After all of the subpopulations are updated, we calculate the fitness again and select the elites of each subpopulation into the elite set. At last, our proposed migration strategy is utilized to propagate advantageous genes. As iterating to the maximum generation, the best fitness and individual are output.Algorithm 1: MPGA with networks.
Input:FJSP: the FJSP to solveNs: the number of subpopulations;S: the size of each subpopulation;Pm: mutation probability;itermax: the maximum generation to run;Output:R∗: the best fitness or makespan;g∗: the best individual;G=G1,G2,…,GM ← encodingFJSP,Ns,S;Whilegenration≤itermax doFor i_sub = 1:Ns doR ← decodingGi_sub;Gsi_sub ← selectionGi_sub based on R;Gci_sub ← crossing Gsi_sub;Gi_sub ← mutation Gci_sub,Pm;End forR ← decodingG;Select the elites of each subpopulation to compose elite setE;G ← migrationG,E;End whileFind the best individualg∗ from E and the corresponding fitness R∗.


## 2.4. Evaluation Indexes

An index to evaluate the propagation rate of advantageous genes among subpopulations under different network structures is needed. All the best individuals of each subpopulation are selected into an elite set. With the accumulation of advantageous genes through the operators of MPGA, the difference between individuals becomes small. Therefore, if the propagation rate of advantageous genes is larger, the difference between these elites will be smaller. Thus, we introduce the HDI to evaluate the difference between elite individuals. The HDI is calculated as follows [24]:(11)HDI=∑i=1100∑j=12×JT1−δhij1,hij2100×2×JT,where we only sample 100 pairs of individuals to limit time consumption and δ… is the Kronecker function. As the two independent variables are unequal, the value of δ… is 0; otherwise, 1. JT represents the total operation number of all jobs, and hij1andhij2 denote the ith pair’s two jth integers in this sample. Therefore, smaller HDI indicates less difference between elites, suggesting the propagation among subpopulations is faster.Besides, we use the SR to assess the performance of MPGA, better performance with higher SR. This measure is defined as follows:(12)SR=NsNt,where Ns denotes the number of times the algorithm finds the optimal value and Nt denotes the total number of times it runs. The optimal value, mentioned in this paper, is defined as the best value found by all current algorithms by now, rather than the best value found by the currently running algorithm.Also, due to the stochastic nature of MPGA, we have to run the proposed algorithm repeatedly, and statistical features are often utilized to analyze simulation results [38–40]. Therefore, the AOV is introduced and calculated as follows:(13)AOV=1Nt∑iNtOVi,in which OVi is the optimal value obtained by the ith run of the algorithm. The smaller the AOV is, the better the performance will be.


## 3. Experiments and Results’ Analysis

In this section, we design experiments to evaluate how the parameters of interaction structures of subpopulations, i.e., ring-shaped networks and small-world networks, affect the performance of MPGAs.We take two benchmark problems proposed by Kacem et al. [41] as examples. One possesses 8 jobs and 8 machines, labeled as P1; the other possesses 10 jobs and 10 machines, labeled as P2. For the basic parameters of MPGA, we set the number of subpopulations Ns as 100, the size of a subpopulation S as 40, the mutation possibility Pm as 0.08, and the total number of iterations itermax as 400.


### 3.1. MPGAs with Ring-Shaped Networks

For ring-shaped networks, the parameterK is changed constantly. The value of K slides from 1 to 24, and the interval is 1. Due to the stochastic nature of evolutionary algorithms, for each value of K, the MPGA with a certain network is run 50 times.Figure7 illustrates the curves of the HDI over iterations for different K. The X-axis represents the values of the HDI, and the Y-axis represents the times of iteration. As shown in the following, the HDI decreases faster as K becomes larger, which indicates the propagation rate of advantageous genes becomes larger as K increases.Figure 7
HDI over iterations with ring-shaped networks. (a) P1. (b) P2.
(a)(b)To further reveal the relationship between the propagation rate and the structural parameter of networks, Figure8 shows the curves of the HDI at 100th iteration, 200th iteration, and 300th iteration along with APL for different K, in which the X-axis represents the value of K, and the Y-axis (left) and Y-axis (right) represent the value of HDI and APL, respectively. It can be seen that, as K becomes larger, the HDI at the same iteration decreases in accordance with the trend of APL shown by cyan lines, namely, the smaller the APL is, the faster the propagation of advantageous genes will be. Also, the variation of parameter K can change the APL to affect the propagation rate of advantageous genes. In this way, the performance of MPGAs can be influenced by different values of K.Figure 8
HDI overK with ring-shaped networks. (a) P1. (b) P2.
(a)(b)To directly show the influence of variation ofK on the performance of MPGAs, curves of SR and AOV over K are shown in Figure 9. For both benchmark problems, when K increases from 1 to 25, the performance of MPGAs increases rapidly at first and then decreases with fluctuation. Concretely, in Figure 9(a), the SR of P1 begins with 24% and increases to the peak, i.e., 90%, when K=4. Then, SR decreases slowly with fluctuation. In contrast, the AOV of P1 begins with 14.9 and decreases to 14.12. Then, AOV rises slowly with fluctuation. For P2, shown in Figure 9(b), the trends of SR and AOV are similar to those in P1. The peak of the SR or the lowest value of AOV can also be found at K=4. For the worst performance of P2, the lowest SR is 4% at K=24, and the highest AOV is 8.42 at K=23.Figure 9
SR and AOV overK with ring-shaped networks. (a) P1. (b) P2.
(a)(b)Therefore, the variation ofK can influence the propagation rate of advantageous genes, consequently affecting the performance of algorithms. In detail, as K is very small, the APL is large, and the advantageous genes can hardly propagate to distant subpopulation in the network. Each subpopulation evolves with little communication and cannot benefit from elites of other subpopulations. However, as K is too large, the APL is very small, and the advantageous genes can be propagated rapidly. Therefore, MPGA is more likely to fall into local optimal, leading to premature convergence. For our benchmark problems, the moderate value of K is 4, so the corresponding propagation rate is neither too large nor too small, and the MPGA with networks can get the best performance.Additionally, because of the stochastic nature of metaheuristics, to show significant differences between MPGAs with different structural parametersK of the network, the nonparametric Wilcoxon signed-rank test of the 50 independent runs is conducted [42]. The null hypothesis, termed as H0, is set as “there exists no difference between MPGAs with different structural parameters of networks.” Accordingly, the alternative hypothesis H1 is “the MPGAs with two different values of the structural parameter are statistically different.” A significance level 0.05 is employed, and a p value less than 0.05 indicates that there exist significant differences between the two samples [43]. As shown in Table 4, we take K=4 vs. other values of K as examples, and the row titled with “total” indicates how many times h=1 within our two benchmark problems. Basically, the larger the difference between two values of K is, more likely h=1, namely, there exist differences between MPGAs with different K, especially, as the difference of K is large. The results further indicate parameter K can influence the performance of MPGAs.Table 4
Wilcoxon signed-rank test results for ring-shaped networks.
Parameter pairsK=4K=4K=4K=4K=4K=4K=4K=4K=4K=4vs.vs.vs.vs.vs.vs.vs.vs.vs.vs.K=1K=2K=6K=9K=11K=13K=15K=18K=21K=25P1p value1.2e − 070.01040.48500.24070.09207.2e − 040.05400.00170.00115.4e − 04h1100010111P2p value1.2e − 050.53160.07090.03390.02820.00510.00763.3e − 043.9e − 051.5e − 05h1001111111Total2101121222


### 3.2. MPGAs with Small-World Networks

For small-world networks, the structural parameters areK and p as mentioned in Section 2.1.3. Since the number of subpopulations (100) is very small, causing the scale of the corresponding network to be small, the rewiring possibility p can bring significant uncertainty, namely, the number of rewired edges fluctuates with the same p. Thus, we directly change the number of rewired edges P as the structural parameter to implement simulation. We suppose the value of P slides from 0 to 120 and the interval becomes larger as P increases. When P=0, the corresponding network is a ring-shaped network. It is noted that, with the same PP≠0, there still exists uncertainty because different edges can be selected to rewire to different locations. To limit this, we randomly generate three different small-world networks with the same P, and MPGA is run 20 times with each network, namely, a total of 60 times for each value of P. At last, for P1, the value of K is set as 2. For P2, the value of Kis set as 4.Figure10 illustrates curves of the HDI over iterations for different P. It can be noted that the HDI decreases faster as P increases, indicating the propagation rate of advantageous genes becomes larger with the increase of P.Figure 10
HDI over iterations with small-world networks. (a) P1. (b) P2.
(a)(b)To intuitively show the influence of variation ofP on the performance of MPGAs, curves of the HDI at 100th, 200th, and 300th iterations over P are shown in Figure 11, along with APL. The X-axis represents the value of P, and the Y-axis (left) and Y-axis (right) represent the value of HDI and APL, respectively. As shown by red, green, and blue lines, with the increase of P, HDI of the same iteration decreases in accordance with the trend of APL shown by cyan lines. It can be concluded that the variation of P can change the APL, thus affecting the propagation rate of advantageous genes, by which the performance of MPGAs is influenced. Compared with the ring-shaped network in Figure 8, the decrease of APL caused by the increase of P is gentler than that caused by the increase of K, resulting in a more gradual change of HDI in Figure 11.Figure 11
HDI overP with small-world networks. (a) P1. (b) P2.
(a)(b)Figure12 directly shows how the performance, measured by SR and AOV, is affected by P, in which the X-axis represents the value of P, and the Y-axis (left) and Y-axis (right) represent the value of SR and AOV, respectively. For P1, shown in Figure 12(a), the SR of P1 starting from 64.44% increases to the peak of 86.67% as P=4 and then decreases with fluctuation. Besides, the AOV starting from 14.41 decreases to the lowest value 14.12 as P=5 and then increases with fluctuation.Figure 12
SR and AOV overP with small-world networks. (a) P1. (b) P2.
(a)(b)In Figure12(b), because P2 is more sophisticated than P1, the success rate is already low. We choose to set K as 4, at which, as shown in Figure 9, MPGA with the ring-shaped network gets the best performance, i.e., SR=38%and AOV=7.64. Therefore, the improvement of performance is limited and subject to fluctuation due to the stochastic nature of evolutionary algorithms. However, the best performance can be found at P=6 at which the peak of the SR is 48.33%, and the lowest value of AOV is 7.5. Besides, for the worst performance, the SR is 18.33%, and AOV is 7.86.Therefore, a similar conclusion can be dawn that the variation of structural parameterP can influence the propagation rate of advantageous genes, thus affecting the performance of MPGAs with small-world networks. Besides, it is noteworthy that the edge number and node number are not changed as P varies, so the change in performance indicates that the MPGA can be improved by choosing a more proper interaction structure of subpopulations as other conditions remain unchanged.Besides, the Wilcoxon signed-rank test of the 60 independent runs is implemented again. As shown in Table5, we take P=6 vs. other values of P as examples. Although there are fewer times h=1, the results still indicate that parameter P can influence the performance of MPGAs with small-world networks.Table 5
Wilcoxon signed-rank test results for small-world networks.
Parameter pairsP=6P=6P=6P=6P=6P=6P=6P=6P=6P=6vs.vs.vs.vs.vs.vs.vs.vs.vs.vs.P=1P=2P=3P=4P=5P=8P=12P=20P=55P=90P1p value0.83480.14440.54850.03890.01640.25130.14440.04950.33590.0762h0001100100P2p value0.0480.07710.41910.01330.09460.01060.01330.02850.0410.0077h1001011111Total1002111211


### 3.3. Comparison between Ring-Shaped Networks and Small-World Networks

We can conclude that both parameterK of ring-shaped networks and unique parameter P of small-world networks can influence the performance of MPGAs, but the degree of influence is different. This is because due to the variation of K and P, the APL is changed to a different extent. For small-world networks, as K=2 and P slides from 0 to 120, APL varies from 12.88 to 3.53; as K=4 and P slides from 0 to 120, APL varies from 7.00 to 2.57. However, for ring-shaped networks, as K slides from 1 to 25, APL varies from 25.25 to 1.50. It can be seen that the decrease of APL is dramatic as K increases, bringing significant influence on the performance, while for P of small-world networks, the variation of P is gentle, causing the influence to be slighter. Therefore, to adjust the APL and the propagation rate of advantageous genes, we can first change parameter K, and the variation of P can be a supplementary method to adjust APL in a smaller step. In this way, we can get a more proper interaction structure to improve the performance of MPGA.


### 3.4. Efficacy of the MPGA with Networks to Solve the FJSP

To further demonstrate the efficacy of MPGA with networks to solve the FJSP, our algorithm is employed to solve small-size FJSPs (SFJ01–10), and larger-size FJSPs (MFJS01–08) [44]. Based on the above analysis, we set K=4 and P=3. Other parameters for MPGA are the same as above. Previous studies, such as AIA [45], HHS [46], M2 and M4 [47], MILP [48], and HA [49], are also experimented on the aforementioned instances, and their results are cited to compare with our algorithm.As shown in Table6, regarding MFJS02 and MFJS03, the optimal value of our algorithm is 446 and 466 (better than 448 and 468 found by AIA). Regarding MFJS04, the optimal value is 554 (better than 564 found by M2 and 590 found by M4). Regarding MFJS05, the optimal value is 514 (better than 527 and 546 found by AIA and M4, respectively). Regarding MFJS06, the optimal value is 634 (better than 635 and 666 found by AIA and M4, respectively). Regarding MFJS07, the optimal value is 879 (better than M2 and M4). At last, for MFJS08, the optimal value is 884, which is not given by M2, M4, and MILP.Table 6
Comparison with other algorithms.
ProblemsOptimal valueAIAHHSM2M4MILPHAThis articleSFJS0166666666666666SFJS02107107107107107107107SFJS03221221221221221221221SFJS04355355355355355355355SFJS05119119119119119119119SFJS06320320320320320320320SFJS07397397397397397397397SFJS08253253253253253253253SFJS09210210210210210210210SFJS10516516516516516516516MFJS01468468468468468468468MFJS02448446446466446446446MFJS03468466466466446466466MFJS04554554564590554554554MFJS05527514514546514514514MFJS06635634634666634634634MFJS078798799281990879879879MFJS08884884———884884


## 3.1. MPGAs with Ring-Shaped Networks

For ring-shaped networks, the parameterK is changed constantly. The value of K slides from 1 to 24, and the interval is 1. Due to the stochastic nature of evolutionary algorithms, for each value of K, the MPGA with a certain network is run 50 times.Figure7 illustrates the curves of the HDI over iterations for different K. The X-axis represents the values of the HDI, and the Y-axis represents the times of iteration. As shown in the following, the HDI decreases faster as K becomes larger, which indicates the propagation rate of advantageous genes becomes larger as K increases.Figure 7
HDI over iterations with ring-shaped networks. (a) P1. (b) P2.
(a)(b)To further reveal the relationship between the propagation rate and the structural parameter of networks, Figure8 shows the curves of the HDI at 100th iteration, 200th iteration, and 300th iteration along with APL for different K, in which the X-axis represents the value of K, and the Y-axis (left) and Y-axis (right) represent the value of HDI and APL, respectively. It can be seen that, as K becomes larger, the HDI at the same iteration decreases in accordance with the trend of APL shown by cyan lines, namely, the smaller the APL is, the faster the propagation of advantageous genes will be. Also, the variation of parameter K can change the APL to affect the propagation rate of advantageous genes. In this way, the performance of MPGAs can be influenced by different values of K.Figure 8
HDI overK with ring-shaped networks. (a) P1. (b) P2.
(a)(b)To directly show the influence of variation ofK on the performance of MPGAs, curves of SR and AOV over K are shown in Figure 9. For both benchmark problems, when K increases from 1 to 25, the performance of MPGAs increases rapidly at first and then decreases with fluctuation. Concretely, in Figure 9(a), the SR of P1 begins with 24% and increases to the peak, i.e., 90%, when K=4. Then, SR decreases slowly with fluctuation. In contrast, the AOV of P1 begins with 14.9 and decreases to 14.12. Then, AOV rises slowly with fluctuation. For P2, shown in Figure 9(b), the trends of SR and AOV are similar to those in P1. The peak of the SR or the lowest value of AOV can also be found at K=4. For the worst performance of P2, the lowest SR is 4% at K=24, and the highest AOV is 8.42 at K=23.Figure 9
SR and AOV overK with ring-shaped networks. (a) P1. (b) P2.
(a)(b)Therefore, the variation ofK can influence the propagation rate of advantageous genes, consequently affecting the performance of algorithms. In detail, as K is very small, the APL is large, and the advantageous genes can hardly propagate to distant subpopulation in the network. Each subpopulation evolves with little communication and cannot benefit from elites of other subpopulations. However, as K is too large, the APL is very small, and the advantageous genes can be propagated rapidly. Therefore, MPGA is more likely to fall into local optimal, leading to premature convergence. For our benchmark problems, the moderate value of K is 4, so the corresponding propagation rate is neither too large nor too small, and the MPGA with networks can get the best performance.Additionally, because of the stochastic nature of metaheuristics, to show significant differences between MPGAs with different structural parametersK of the network, the nonparametric Wilcoxon signed-rank test of the 50 independent runs is conducted [42]. The null hypothesis, termed as H0, is set as “there exists no difference between MPGAs with different structural parameters of networks.” Accordingly, the alternative hypothesis H1 is “the MPGAs with two different values of the structural parameter are statistically different.” A significance level 0.05 is employed, and a p value less than 0.05 indicates that there exist significant differences between the two samples [43]. As shown in Table 4, we take K=4 vs. other values of K as examples, and the row titled with “total” indicates how many times h=1 within our two benchmark problems. Basically, the larger the difference between two values of K is, more likely h=1, namely, there exist differences between MPGAs with different K, especially, as the difference of K is large. The results further indicate parameter K can influence the performance of MPGAs.Table 4
Wilcoxon signed-rank test results for ring-shaped networks.
Parameter pairsK=4K=4K=4K=4K=4K=4K=4K=4K=4K=4vs.vs.vs.vs.vs.vs.vs.vs.vs.vs.K=1K=2K=6K=9K=11K=13K=15K=18K=21K=25P1p value1.2e − 070.01040.48500.24070.09207.2e − 040.05400.00170.00115.4e − 04h1100010111P2p value1.2e − 050.53160.07090.03390.02820.00510.00763.3e − 043.9e − 051.5e − 05h1001111111Total2101121222


## 3.2. MPGAs with Small-World Networks

For small-world networks, the structural parameters areK and p as mentioned in Section 2.1.3. Since the number of subpopulations (100) is very small, causing the scale of the corresponding network to be small, the rewiring possibility p can bring significant uncertainty, namely, the number of rewired edges fluctuates with the same p. Thus, we directly change the number of rewired edges P as the structural parameter to implement simulation. We suppose the value of P slides from 0 to 120 and the interval becomes larger as P increases. When P=0, the corresponding network is a ring-shaped network. It is noted that, with the same PP≠0, there still exists uncertainty because different edges can be selected to rewire to different locations. To limit this, we randomly generate three different small-world networks with the same P, and MPGA is run 20 times with each network, namely, a total of 60 times for each value of P. At last, for P1, the value of K is set as 2. For P2, the value of Kis set as 4.Figure10 illustrates curves of the HDI over iterations for different P. It can be noted that the HDI decreases faster as P increases, indicating the propagation rate of advantageous genes becomes larger with the increase of P.Figure 10
HDI over iterations with small-world networks. (a) P1. (b) P2.
(a)(b)To intuitively show the influence of variation ofP on the performance of MPGAs, curves of the HDI at 100th, 200th, and 300th iterations over P are shown in Figure 11, along with APL. The X-axis represents the value of P, and the Y-axis (left) and Y-axis (right) represent the value of HDI and APL, respectively. As shown by red, green, and blue lines, with the increase of P, HDI of the same iteration decreases in accordance with the trend of APL shown by cyan lines. It can be concluded that the variation of P can change the APL, thus affecting the propagation rate of advantageous genes, by which the performance of MPGAs is influenced. Compared with the ring-shaped network in Figure 8, the decrease of APL caused by the increase of P is gentler than that caused by the increase of K, resulting in a more gradual change of HDI in Figure 11.Figure 11
HDI overP with small-world networks. (a) P1. (b) P2.
(a)(b)Figure12 directly shows how the performance, measured by SR and AOV, is affected by P, in which the X-axis represents the value of P, and the Y-axis (left) and Y-axis (right) represent the value of SR and AOV, respectively. For P1, shown in Figure 12(a), the SR of P1 starting from 64.44% increases to the peak of 86.67% as P=4 and then decreases with fluctuation. Besides, the AOV starting from 14.41 decreases to the lowest value 14.12 as P=5 and then increases with fluctuation.Figure 12
SR and AOV overP with small-world networks. (a) P1. (b) P2.
(a)(b)In Figure12(b), because P2 is more sophisticated than P1, the success rate is already low. We choose to set K as 4, at which, as shown in Figure 9, MPGA with the ring-shaped network gets the best performance, i.e., SR=38%and AOV=7.64. Therefore, the improvement of performance is limited and subject to fluctuation due to the stochastic nature of evolutionary algorithms. However, the best performance can be found at P=6 at which the peak of the SR is 48.33%, and the lowest value of AOV is 7.5. Besides, for the worst performance, the SR is 18.33%, and AOV is 7.86.Therefore, a similar conclusion can be dawn that the variation of structural parameterP can influence the propagation rate of advantageous genes, thus affecting the performance of MPGAs with small-world networks. Besides, it is noteworthy that the edge number and node number are not changed as P varies, so the change in performance indicates that the MPGA can be improved by choosing a more proper interaction structure of subpopulations as other conditions remain unchanged.Besides, the Wilcoxon signed-rank test of the 60 independent runs is implemented again. As shown in Table5, we take P=6 vs. other values of P as examples. Although there are fewer times h=1, the results still indicate that parameter P can influence the performance of MPGAs with small-world networks.Table 5
Wilcoxon signed-rank test results for small-world networks.
Parameter pairsP=6P=6P=6P=6P=6P=6P=6P=6P=6P=6vs.vs.vs.vs.vs.vs.vs.vs.vs.vs.P=1P=2P=3P=4P=5P=8P=12P=20P=55P=90P1p value0.83480.14440.54850.03890.01640.25130.14440.04950.33590.0762h0001100100P2p value0.0480.07710.41910.01330.09460.01060.01330.02850.0410.0077h1001011111Total1002111211


## 3.3. Comparison between Ring-Shaped Networks and Small-World Networks

We can conclude that both parameterK of ring-shaped networks and unique parameter P of small-world networks can influence the performance of MPGAs, but the degree of influence is different. This is because due to the variation of K and P, the APL is changed to a different extent. For small-world networks, as K=2 and P slides from 0 to 120, APL varies from 12.88 to 3.53; as K=4 and P slides from 0 to 120, APL varies from 7.00 to 2.57. However, for ring-shaped networks, as K slides from 1 to 25, APL varies from 25.25 to 1.50. It can be seen that the decrease of APL is dramatic as K increases, bringing significant influence on the performance, while for P of small-world networks, the variation of P is gentle, causing the influence to be slighter. Therefore, to adjust the APL and the propagation rate of advantageous genes, we can first change parameter K, and the variation of P can be a supplementary method to adjust APL in a smaller step. In this way, we can get a more proper interaction structure to improve the performance of MPGA.


## 3.4. Efficacy of the MPGA with Networks to Solve the FJSP

To further demonstrate the efficacy of MPGA with networks to solve the FJSP, our algorithm is employed to solve small-size FJSPs (SFJ01–10), and larger-size FJSPs (MFJS01–08) [44]. Based on the above analysis, we set K=4 and P=3. Other parameters for MPGA are the same as above. Previous studies, such as AIA [45], HHS [46], M2 and M4 [47], MILP [48], and HA [49], are also experimented on the aforementioned instances, and their results are cited to compare with our algorithm.As shown in Table6, regarding MFJS02 and MFJS03, the optimal value of our algorithm is 446 and 466 (better than 448 and 468 found by AIA). Regarding MFJS04, the optimal value is 554 (better than 564 found by M2 and 590 found by M4). Regarding MFJS05, the optimal value is 514 (better than 527 and 546 found by AIA and M4, respectively). Regarding MFJS06, the optimal value is 634 (better than 635 and 666 found by AIA and M4, respectively). Regarding MFJS07, the optimal value is 879 (better than M2 and M4). At last, for MFJS08, the optimal value is 884, which is not given by M2, M4, and MILP.Table 6
Comparison with other algorithms.
ProblemsOptimal valueAIAHHSM2M4MILPHAThis articleSFJS0166666666666666SFJS02107107107107107107107SFJS03221221221221221221221SFJS04355355355355355355355SFJS05119119119119119119119SFJS06320320320320320320320SFJS07397397397397397397397SFJS08253253253253253253253SFJS09210210210210210210210SFJS10516516516516516516516MFJS01468468468468468468468MFJS02448446446466446446446MFJS03468466466466446466466MFJS04554554564590554554554MFJS05527514514546514514514MFJS06635634634666634634634MFJS078798799281990879879879MFJS08884884———884884


## 4. Conclusion and Future Scope

For further analyzing the influence of interaction structures on the performance of MPGAs, MPGA based on complex networks, such as ring-shaped networks and small-world networks, has been presented in this work. Then, how structural parametersK and Pof these two types of networks affect the performance of MPGAs is discussed by solving the FJSP.First, the HDI is utilized to quantitatively measure the differences of elites between subpopulations. The smaller the HDI is, the higher the propagation rate of advantageous genes will be. The simulation results indicate that the HDI has a positive correlation with the APL. This means the smaller the APL is, the higher the propagation rate is.Next, as the structural parametersK and P increase, the curves of the HDI over iterations decrease faster, namely, the propagation rate becomes faster with the increases of K and P, resulting in the variation of the performance of MPGAs captured by SR and AOV. For instance, the highest and lowest SR of P2 are 48.33% and 4.00%. Since the propagation rate of advantageous genes can neither be too high nor too low, we should choose proper K and P for MPGAs with networks. Based on experiments, the proper intervals for K can be [2, 5], and the choice of P is dependent on the value of K because the variation of K can bring a significant change of APL, and the variation of P can be a supplementary method to adjust APL in a smaller step.At last, the edge number and node number of small-world networks are not changed asP varies, so the change in performance indicates that MPGA can be improved by choosing a more proper interaction structure of subpopulations as other conditions remain unchanged.Because of the limitation of computing capability, the number of nodes is not big enough compared with classic network science, but this study can still illustrate the importance of interaction structures and the meaning of introducing complex networks as a study tool. In the future, more complex or compound networks, such as multilayer networks, can be introduced into the research.



---

*Source: 8503454-2020-11-27.xml*