Memetic Algorithm With Meta-Lamarckian Learning and Simplex Search for Distributed Flexible Assembly Permutation Flowshop Scheduling Problem

This paper studies a novel and practical distributed flexible assembly permutation flowshop scheduling problem with makespan criterion, which has attracted wide attention due to important applications in modern manufacturing. The problem integrates two machine environments of distributed production and flexible assembly, which can process and assemble the jobs into customized products. We first present a mixed integer linear programming model to characterize the problem essence and to solve small-size problems. Due to the NP-hard, we further propose an efficient memetic algorithm, which consists of a global exploration optimizer designed based on improved social spider optimization and two local exploitation optimizers designed based on meta-Lamarckian learning and simplex search, respectively. To implement the algorithm, a problem-specific encoding scheme is presented. Algorithmic parameters are calibrated by a design of experiments, and a comprehensive computational campaign is conducted to evaluate the performance of the mathematical model and algorithms. Statistical results show that their problem-solving abilities are effective, and especially the proposed memetic algorithm outperforms the existing algorithms significantly.


I. INTRODUCTION
Flowshop scheduling problem (FSP) has been extensively studied in operational research and computer science due to its theoretical significance and practical applications [1]- [5]. In flowshop layout, the jobs are processed on a set of machines disposed in a fixed order. The typical FSP assumes that there is only a manufacturing unit or factory to be used to complete all tasks. With the deepening of economic globalization, however, modern enterprises are more and more dragged into the market with cutthroat competitiveness. Therefore, it becomes popular to duplicate multiple manufacturing factories to add flexibility into the production campaign [6]. This leads to the manufacturing reform from The associate editor coordinating the review of this manuscript and approving it for publication was Chih-Yu Hsu .
typical single FSP to the cooperation of multiple FSP, i.e., distributed FSP (DFSP). DFSP can decompose total production tasks into distinct manufacturing units, which helps enterprises to harvest many advantages such as better corporate reputation, higher product quality, and lower production cost and period. Compared to FSP, the scheduling in distributed environment is more complicated considering two relevant decisions must be addressed at the same time: assign the jobs to candidate factories and solve a FSP at each factory. Naderi and Ruiz [7] handled such a distributed permutation flowshop scheduling problem (DPFSP) for the first time.
The economic trend towards customized manufacturing leads to increasing prevalence for assembly manufacturing system in recent decades. As a generalization of FSP, the assembly flowshop scheduling problem (AFSP) was proposed by Lee et al. [8], in which the layout is a flowshop for VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ job-processing plus a single machine for product-assembling. Due to the limitation of single machine assembly, it is detrimental to the concurrent and mass production. Therefore, more flexible or multiple machine assembly comes into view. Mozdgir et al. [9] and Navaei et al. [10] are the pioneers in this field. They considered a two stage assembly flowshop scheduling with unrelated second stage assembly machines. Recently, the integration scheduling problem for distributed flowshops plus an extra single-machine assembly stage was first handled by Hatami et al. [11], which is referred to as the distributed assembly permutation flowshop scheduling problem (DAPFSP). In addition to two decisions of factory assignment and job sequencing at each factory, DAPFSP must address the third decision: product sequencing on assembly machine. The disadvantage of single-machine assembly, however, still exists and even becomes more severe due to the higher productivity from the distributed production stage. For this reason, most enterprises adopt the multi-machine flexible assembly mode in actual production. This paper adopts the flexible assembly layout to replace the single machine assembly in DAPFSP, which generates a more practical distributed scheduling problem with flexible assembly, termed distributed flexible assembly permutation flowshop scheduling problem (DFAPFSP). An example of DFAPFSP can be found in the manufacturing of automobile engines. Different parts such as bent axles, cylinder blocks and cylinder heads are first processed using one of multiple flowshops, and then assembled into the complete engine on one of some unrelated assembly machines. The DFAPFSP is NP-hard due to the NP-hard of its special cases such as FSP, AFSP, DFSP, and DAPFSP. Hence, the evolutionary algorithm (EA) is one of the best candidates for solving the DFAPFSP. A high-quality solution in optimization problem can be gained by EA mainly depending on the good balance of exploration and exploitation abilities. In this sense, EA is not competent due to its insufficient local search ability. So, researchers devised new optimization model, referred to as memetic algorithm (MA). Usually, MA consists of a global exploration optimizer and a group of local exploitation optimizers, which are constructed based on the problem-specific knowledge and implemented in a complementary way [12]. According to the unique design preferences and evolutionary framework of MA [13], this paper presents an efficient MA for the DFAPFSP, and its main contributions are summarized as follows: a) Propose a global exploration optimizer by refining an efficient social spider optimization, which possesses the abilities of bi-population co-evolution and cross-population interaction; b) Propose a meta-Lamarckian learning strategy-based local exploitation optimizer to refine the best solution in population; c) Propose a simplex search-based local exploitation operator to refine the diversity and the overall quality of the population; d) Formulate a novel MA to cope with the DFAPFSP based on the above global and local optimizers. The rest of this paper proceeds as follows. Section II reviews the closely related literature. Section III formulates the problem under study. Section IV details the proposed memetic algorithm. Section V conducts a comprehensive computational experiment. Finally, Section VI concludes this paper with some discussions and future research topics.

II. LITERATURE REVIEW
The problem studied in the paper is a new and promising research topic. To the best of our knowledge, the work on this problem published at present is very limited although it has been focused on widely. We review the closely related work as follows.

A. DISTRIBUTED PRODUCTION SCHEDUULING PROBLEM
In distributed production scheduling problems, DFSP has obtained the most attention and gained the most extensive research achievements. In the pioneer work of Naderi and Ruiz [7], they proposed multiple mathematical models and heuristics to solve the problem. Since then, a variety of the heuristic and metaheuristic methods have been devised for it, such as estimation of distributed algorithm (EDA) [14], [15], chemical reaction optimization [16], iterated greedy algorithm [17]- [20], scatter search algorithm [21], differential evolution (DE) algorithm [22], [23].
The first attempt that addressed the AFSP in distributed manufacturing setting, i.e., DAPFSP, was made by Hatami et al. [11], and a mathematical model and some heuristics and metaheuristics were proposed to solve it. For DAPFSP, Wang and Wang [24] constructed an EDA-based memetic algorithm. Lin and Zhang [25] established a biogeography-based optimization algorithm. Lin et al. [26] presented a backtracking search hyperheuristic method. Sang et al. [27] proposed an invasive weed optimization algorithm. Ferone et al. [28] devised a biased-randomized iterated local search algorithm. Ochi and Driss [29] presented a bounded-search iterated greedy algorithm. To add more flexibility, Zhang et al. [30] first addressed the distributed flexible assembly permutation flowshop scheduling problem, i.e., DFAPFSP under study. In this pioneer work, a constructive heuristic and two hybrid algorithms were designed based on variable neighborhood search and particle swarm optimization. In addition to DFSP, the distributed scheduling problems with non-flowshop layout are also extensively concerned [6], [31] in recent decades.

B. MEMETIC ALGORITHM AND RELATED STRATEGIES
Due to the remarkable optimization ability, MA has gained more and more attention along with successful applications in distinct fields [12]. Recently, a variety of MAs have been built to solve the production scheduling problems. Ishibuchi et al. [32] proposed a memetic genetic algorithm for solving the multiobjective permutation FSP. Wang and Wang [24] devised an EDA-based MA for DAPFSP. Wu and Che [33] constructed a memetic DE for the energy-efficient parallel machine scheduling problem. Liu et al. [34] established a PSO-based MA for FSP. Pan et al. [35] proposed a high-performing MA for FSP with blocking constraints. Qian et al. [36] devised a DE-based MA for multiobjective jobshop scheduling problem. Gonzalez and Vela [37] established an efficient MA for the single machine scheduling problem. Shao et al. [38] developed a MA based on node and edge histogram for no-idle flowshop scheduling problem. Deng and Wang [39] gave a competitive MA for multiobjective DPFSP.
The social spider optimization (SSO) is a population-based EA, which is inspired by the cooperative behavior among social spider colony [40]. At present, it has been successfully applied to a lot of optimization problems [31]. The key feature of SSO is that its population is divided into two subpopulations, which first evolves based on different mechanisms, and then interacts between them to keep the population diversity. With such bi-population co-evolution and cross-population interaction, it shows excellent global exploration ability. For this reason, we improve the basic SSO to perform global search in the proposed MA. One of the features in MA is the application of multiple strategies to be local search. But how to arrange them in the evolution significantly affects the efficiency of MA. Ong and Keane [41] pointed out that meta-Lamarckian learning strategy can realize such arrangement in an adaptive way. The simplex search is an effective local search technique, which is easy to implement and has been successfully applied to various domains [42]. Because the original simplex search is weak in global exploration and its performance depends on the initial value, it is usually used by the hybridization with EA. Inspired by these facts, we establish two local exploitation optimizers for the proposed MA based on the meta-Lamarckian learning strategy and simplex search method, respectively.
From the above literature review, it can be seen that DFSP has been a hot research topic, and a large number of achievements have been reported in recent years. However, the investigation on the DFAPFSP is still very preliminary although it has shown important and extensive applications in modern manufacturing. Therefore, we urgently need to present more effective and efficient algorithms to address this problem.

III. DISTRIBUTED FLEXIBLE ASSEMBLY PERMUTATION FLOWSHOP SCHEDULING PROBLEM
A. PROBLEM DESCRIPTION Figure 1 provides the structure of DFAPFSP, and Table 1 lists all necessary notations throughout the paper. In DFAPFSP, the layout consists of two consecutive stages: distributed production and flexible assembly. The production stage has f identical factories working in parallel, each of which is set m production machines in the flowshop manner. The assembly stage can be viewed as a single factory with g unrelated assembly machines. In the manufacturing process, there are t different products to be manufactured by assembling n jobs  according to a predefined assembly plan. To be specific, n jobs are first assigned to f distributed factories and at the same time specified a processing order in each factory. After all jobs in the assembly plan of a product are finished at the production stage, they are assembled into the final product by any a machine at the assembly stage. The objective is to find a schedule for jobs and products to minimize makespan. In addition to the constraints in traditional FSP, this paper also assumes: each product has an assembly plan and each job belongs to a special product. In distributed factories, each job must exactly select one factory and each factory can operates all jobs. If a job is allotted to one factory, the transfer to other factories is not allowed. In flexible assembly, there are a few unrelated machines, where each product must exactly select VOLUME 8, 2020 one machine and each machine can assemble all products. A product can start assembly only when its assembly plan is finished in distributed factories. There are infinite buffers between machines and between manufacturing stages.

B. MAKESPAN CALCULATION
Let λ k = [λ k (1), λ k (2), . . ., λ k (n k )] be the job processing sequence in factory k at the production stage, where n k is the number of jobs assigned to factory k. Then, the job schedule in distributed factories can be represented as λ = {λ 1 , λ 2 , . . . , λ f }, i.e., a set of job sequences. Denote P i,j and CP i,j as the processing time and the completion time of job i on processing machine j, respectively. It follows that: Denote C l as the latest completion time of jobs in assembly plan AP l of product l at the production stage. It is calculated as follows: Let ξ s = [ξ s (1), ξ s (2), . . ., ξ s (t s )] be the assembling sequence of t s products assigned to assembly machine s. The product schedule on all assembly machines can be represented as ξ = {ξ 1 , ξ 2 , . . . , ξ q }, i.e., a set of product sequences. Denote T l,s and CA l,s as the assembly time and completion time of product l on assembly machine s, respectively. Then, the following calculation can be available: Finally, we combine the job schedule λ and product schedule ξ to gain a complete DFAPFSP schedule π = (λ; ξ ). This schedule consists of four consecutive decisions: assign jobs to candidate factories, sequence jobs in each of distributed factory, assign products to assembly machines and arrange products on each assembly machine. Denote the makespan of π as C max (π ). The objective is to find a schedule π with minimum C max (π ). It can be calculated by:

C. ILLUSTRATIVE EXAMPLE
To understand the problem in depth, an example is presented to illustrate it. It consists of n = 9, t = 3, f = 3, m = 3 and g = 2. The assembly plan is that AP 1 = {J 1 , J 6 , J 7 }, AP 2 = {J 2 , J 5 , J 9 } and AP 3 = {J 3 , J 4 , J 8 } corresponding to products P 1 , P 2 and P 3 , respectively. Table 2 presents the operating times (P ij and T ls ) for jobs and products on different machines. Consider problem schedule π = (λ; ξ ), where With the recursive calculation (1)- (7), it is easy to calculate the starting and completion time of each job and product in this schedule on corresponding machines. They are also provided in Table 2 by columns st, CP ij and CA ls , respectively. Based on these data, the makespan is calculated as 24 by Eq. (8), and the Gantt chart is presented in Figure 2 to show the scheduling details.

D. MATHEMATICAL MODEL
To formulate the mathematical description, a mixed integer linear programming (MILP) model is built in this section. It is inspired by the research in [11], which is a special case of the work in this paper. The indices, parameters and variables used in this model see Table 1, and the MILP model can be described as follows: The objective function (9) shows the makespan to be minimized. In the distributed production stage, the model should subject to the constraints: Constraint (10) shows that each job has only one predecessor. Constraint (11) controls that each job has at most one successor. Constraint (12) assures that dummy job 0 must be a predecessor of another job f times, while constraint (13) shows that dummy job 0 should be a successor f − 1 times. The existence of dummy job 0 in such way can divide the job sequence into f parts, each of which represents a factory's job schedule. Thus, there is no factory index to be used in the constraints. Constraint (14) points out that a job cannot be another job's predecessor and successor at the same time. Constraints (15) and (16) indicate that at a time a job is processed by at most one machine and a machine can process at most one job.
In the flexible assembly stage, the model has to subject to the constraints: CA l,s ≥ C l + T l,s ∀l, s Constraints (17) and (18) indicate that each product has only one predecessor and at most one successor. Constraints (19) and (20) enforce that every product must be exactly assigned to one assembly machine. Constraint (21) specifies that a product cannot be at the same time another product's predecessor and successor. Constraint (22) states that each product cannot start assembly before its assembly plan is finished in distributed factories. Constraint (23) assures that if product l is assembled immediately after product v, its assembly on machine s cannot start until machine s finishes product v. Constraint (24) defines the makespan.
Finally, the domain of the decision variables is defined by the following constraints:

IV. MEMETIC ALGORITHM WITH META-LAMARCKIAN LEARNING AND SIMPLEX SEARCH
Due to the NP-hardness, the above mathematical model is only able to solve some small-size problems effectively. Whereas for the large-size problems, it is hard to obtain the optimal solutions within a reasonable computational time. For solving the DFAPFSP with different sizes, we propose a memetic algorithm that mainly consists of a global exploration optimizer devised by improving social spider optimization and two local exploitation optimizers devised separately based on meta-Lamarckian learning and simplex search. It is denoted as MAGL in short and the details can be depicted as follows.

A. ENCODING AND DECODING SCHEMES
Encoding scheme is a key factor in the algorithm design [43]. First, it affects the number and distribution of solutions in the search space. Second, it affects the neighborhood structure and size of a solution such that affecting the time cost. Third, along with the decoding scheme, it affects the solution quality. Thus, the encoding and decoding schemes should be well designed based on the problem knowledge.
In this paper, a problem solution is represented as a vector x = [x 1 , x 2 , . . . , x n ] of n real numbers, each of which are randomly generated in interval [1, f +1). To decode x as a complete DFAPFSP schedule, the following two steps are performed: Step 1: Construct the job schedule λ in distributed factories. Let x i correspond to job J i and denote |x i | as the integer part of x i . Assign job J i to factory |x i |, and the jobs in the same factory are processed in ascending order of corresponding numbers in x. In this way, the job schedules in distributed factories are obtained. Consider f = 2 and an example x = [1.26, 2.33, 1.67, 2.75, 2.14, 1.91, 1.41] to explain this step. Since x 1 = 1.26, x 3 = 1.67, x 6 = 1.91 and x 7 = 1.41 have the same integer part 1, the jobs J 1 , J 3 , J 6 and J 7 are assigned VOLUME 8, 2020 to factory 1. According to the order of x 1 < x 7 < x 3 < x 6 , the job processing sequence in factory 1 can be determined to be λ 1 = [J 1 , J 7 , J 3 , J 6 ]. Similarly, there is λ 2 = [J 5 , J 2 , J 4 ], and therefore the job schedule in distributed factories is λ = {λ 1 , λ 2 }.
Step 2: Construct the product schedule ξ on flexible assembly machines. For any two products, the one whose assembly plan is completed earlier in the distributed factories is assigned first. That is to say, if C l is less than C h , product l is assigned prior to product h. For each product to be assigned first, it is assembled by any one assembly machine that completes this product at the earliest possible time. If there are more than one choice, we select one from the candidates randomly. With this scheme, the product schedule ξ and then a complete problem schedule π = (λ; ξ ) can be determined.
As stated before, a complete DFAPFSP schedule should contain four interrelated decisions. If all these decisions are considered in encoding scheme, the solution representation may be very complex, resulting in the search inefficiency due to large search space. Thus, we propose the incomplete encoding scheme above. The decisions missed in encoding scheme will be dynamically retrieved during the decoding process.

B. INITIALIZATION AND ACCEPTANCE CRITERION
The initial population is generated according to the idea of avoiding any idle factories, and therefore the following steps are adopted: Step 1: Randomly generate vector x = [x 1 , x 2 ,. . . , x n ] from the interval [1, f + 1) according to Eq. (28), i.e., set x min = 1 and x max = f + 1.
where rand is a random number in [0, 1].
Step 3: Replace y r by the value randomly generated from the interval [r, r + 1), r = 1, 2, . . ., f , which ensures that each factory is assigned at least one job.
Step 4: Repeat the above steps until meeting the specified population size.
The acceptance criteria determine which of the old and new solutions can survive in next population. So, its choice affects an algorithmic performance to a certain extent. In this work, we give two acceptance criteria and denote them as AC 1 and AC 2 , respectively. AC 1 is a biased acceptance criterion and is applied only in the global search module. For the minimization problem under study, it is formulated as follows: where y is the survivor, τ ∈ (0, 1) is the bias acceptance constant, and RE is the relative fitness error of x new regarding x old and is defined by RE = (C max (x new )− C max (x old )) / C max (x old ). As can be seen, x new surely survives, i.e., y = x new , if it is absolutely better than x old , i.e., RE < 0. When x new is only slightly worse than x old , i.e., 0 < RE < τ , it may also survive by a small probability τ − RE.
In this way, τ can control several solutions with relatively poor quality to survive in next population. The population diversity is maintained to some extent. AC 2 is a simple greedy or descent acceptance criterion, i.e., accepting new solution only if it is absolutely better than the old one. This criterion is used only in the local search module.

C. GLOBAL EXPLORATION OPTIMIZER 1) IDEA OF BASIC SSO
Basic SSO models the search space as a spider web, where spiders interact with each other by vibrating the web. Each position the spiders go through corresponds to a problem solution. All spider positions make up the SSO population. Each spider is given a weight that is calculated by its related fitness function, and the greater the weight, the better the quality of the solution. During the search process, the SSO uses two search agents: male and female spiders. Accordingly, the population is divided into two subsets based on the spiders' gender. The spider positions in these two subsets are updated by different operators, and then a mating operator is carried out between two subsets. The detailed steps of the basic SSO can refer to the work in [40]. Compared to traditional EAs, the SSO has good global search ability, which derived from its ingenious design: bi-population co-evolution and cross-population interaction.

2) IMPROVED SSO FOR MAGL GLOBAL OPTIMIZER
According to the template of the basic SSO, we propose an improved version of it, which has five modules: population division, weight calculation, updating female positions, updating male positions, and mating operation. The improved SSO is used for the MAGL global optimizer. Next, we detail each module as follows.

a: POPULATION DIVISION AND WEIGHT CALCULATION
In social spider colony, the number of females N f accounts for 65%-90% of total population number N p , and of course the other is the number of males N m . They are calculated by: where floor(·) represents a function that maps real number to the integer. The first N f members of total population S is divided into female subset F, and the other members are the male subset M . That is, S = F ∪ M .
Each spider x receives a weight w which reflects the solution quality that corresponds to this spider. To model the weight, the following formula is adopted: where c 1 and c 2 are two constants and used to avoid zero denominator and zero weight. In this paper, we set c 1 = 0.01 and c 2 = 0.001.

b: UPDATING FEMALE POSITIONS
The social spiders update their positions in the web by interacting with each other and learning the external stimulus. The spiders with different gender update the position in different way. The research in Zhang and Xing [31] showed that a female spider migrates its position depending on: the cooperation with the closest and more weighted spider, the cooperation with the spider with the highest weight in the current population, and the uncertain external stimulus. Thus, a female spider moves its position from x to y to be modeled as follows: where r 1 , r 2 , r 3 and r 4 are four random numbers in [0, 1], x c and w c are the position and weight of the closest and more weighted spider, x b and w b are the position and weight of the spider with the highest weight in the population, notation || · || calculates the Euclidean distance, v c,x is the vibration perceived by the spider at x from the spider at x c , v b,x is the vibration perceived by the spider at x from the spider at x b , the last random item in Eq. (33) produces an unbiased perturbation along any directions, and κ controls the move near or far away, which works in this paper according to the scheme: if rand < 0.7, then κ = 1, else κ = −1.

c: UPDATING MALE POSITIONS
According to the role played in the colony, male spiders are divided into two categories: dominant and non-dominant male spiders. Dominant males have larger weight than the non-dominant. Dominant male spiders usually move towards the nearest female spider for reproduction, whereas nondominant male spiders prefer to concentrate in the center of the male colony for making use of the resources wasted by dominant males. In this paper, half of the male spiders with greater weight are referred to as dominant and the other half are called non-dominant. For male spiders, the position migration from x to y is modeled by: x mean where x a and w a represent the position and weight of the nearest female spider to male spider at x, v a,x is the vibration perceived by the male spider at x from the female spider at x a , M is the subset of all male spiders, x mean is the weighted mean of all male positions.

d: IMPROVED ADAPTIVE MATING OPERATOR
After both female and male positions are updated, a mating operator will be executed on dominant male spiders in current colony. In the original mating operator, each dominant male spider has a mating radius and mates with the females within this radius to produce new baby spider. The position value of the baby spider at every dimension is randomly chosen from the corresponding dimensional values of the mating candidates by roulette-wheel technique based on their weights. If there are no female spiders within this radius, the mating operator is not performed. However, such mating operation has two limitations: a) there may be no female spiders in the specified mating radius for a long time; b) it is difficult to generate better baby spiders in the later search stage due to lack of population diversity. Therefore, we improve the original mating operator to avoid such two shortcomings, and the new version can be depicted as follows.
where i = 1, 2, . . . , n, t is the current iteration number, t max is the allowable maximum iteration number, Pi ∈ [P min , P max ] can be viewed as the inheritance rate for the mating parents. As can be seen from Figure 3, Pi dynamically drops as sigmoid curve, which make the baby spider inherit more genes of mating parents in the beginning search stage, and gradually increase variation with iteration. In this way, we can maintain the population diversity as much as possible during the evolutionary process. New baby spider is used to replace the worst individual in population, and inherit the same gender.

D. LOCAL EXPLOITATION OPTIMIZERS
The global optimizer proposed above focuses more on the global exploration, but its local exploitation capacity is relatively insufficient. In order to balance them effectively, two local exploitation optimizers are further proposed as follows. Based on the four job moves, two compound job moves are further defined by an integration of them: CM 1 is JM 1 plus JM 3 , i.e., perform JM 1 and JM 3 simultaneously, and CM 2 is JM 2 plus JM 4 .
In the search of MLS, we apply JM k , k = 1, 2, 3, 4, to perform a slight perturbation on the search candidates for preventing falling into a local trap, and then utilize CM k , k = 1, 2, to conduct the local exploitation on the perturbation results. Inspired by the work in [34], we not only consider two distinct compound job moves in local exploitation, but also adopt an adaptive meta-Lamarckian learning technique [41] to decide which compound job move is more suitable for current local exploitation.
This meta-Lamarckian learning method consists of a training phase and a working phase. In the training phase, each of CM k , k = 1, 2, is executed Tim = n× (n−1) times on the best member within the initial or first-generation population of our proposed MAGL algorithm. Then for each move, the reward value and selection probability are calculated using Eq. (41) and Eq. (42), respectively.
where C bef is the makespan value before performing the move, and C aft is the best makespan value generated during Tim times executions by the AC 2 acceptance criterion. The working phase is performed on the best member in the non-initial MAGL population. For the search candidate, we first apply one of JM k , k = 1, 2, 3, 4, to perturb it, and then use the roulette-wheel strategy to decide which one of CM 1 and CM 2 to be selected for improving the perturbing solution based on selection probability. If CM k is chosen, it will be carried out Tim times to update its reward value via η k = η k + η new , where η new is the newly obtained reward also calculated according to Eq. (41). Accordingly, the selection probability of CM 1 and CM 2 is updated by Eq. (42).
It should be noted that the search result generated from MLS is job permutations, which cannot be evolved by the proposed MAGL global optimizer. So, a conversion from λ to x is needed after implementing MLS. The conversion is stated by Algorithm 1, and the process of MLS is presented by Algorithm 2. find the best member λ in current population 7: k = 1 8: while k ≤ 4 do{ 9: λ ← perturb λ by JM k 10: l = 0 11: while l < Tim do{ 12: select CM 1 or CM 2 based on ζ 1 and ζ 2 by roulette wheel rule 13: λ ← perform the selected move Tim times on λ 14: λ ← select the better of λ and λ 15: calculate η new by Eq. (41) 16: η k ← η k + η new 17: update ζ 1 and ζ 2 by Eq. (42) 18: l = l + 1} 19: λ ← select the better of λ and λ 20: if λ is improved then k = 1 21: else k = k + 1}

2) SLS
Simplex-based local search. During the process of iterating the above global and local optimizers, the optimal individual of each iteration may fall into local optimum. So, it is necessary to find new solutions to jump out of the local optimum. Simplex method is a direct optimization method, which does not need to know gradient information and has strong local search ability [42]. Thus, a simplex-based local search (SLS) method is proposed herein. It first constructs a simplex using current population, and rescales it to generate a solution (spider position) with better fitness value based on four operators: reflection, expansion, contraction and shrink. Then, the worst solution in population is replaced by this solution generated newly. In this way, the population's overall quality can be gradually improved to approximate the global optimum. One iteration of SLS can be described as follows: Step 1: Reflection. Find the spiders with the global best and the worst position in current population. Denote them as x gb and x w , respectively. Calculate x cent as the average of the swarm positions excluding x w . Generate a new position x ref by reflecting x w according to Eq. (43).
, go to step 4.
Step 2: Expansion. Expand the reflection position x ref along the same direction as the reflection. This operation is performed by Eq. (44). If C max (x exp ) < C max (x gb ), the best one of x exp and x ref replaces Step 3: Contraction. Perform the contraction operation using Eq. (45). If C max (x con ) < C max (x w ), x con replaces x w .
Step 4: Shrink. Perform the shrink operation by Eq. (46). If C max (x shr ) < C max (x w ), the best one of x shr and x ref replaces x w , otherwise x ref replaces x w .
where α = 1, β = 2, and γ = δ = 0.5. For a population with only three individuals x gb , x w and x, these four operators is shown in Figure 4.

E. OVERALL MAGL ALGORITHM
With the steps designed previously, the overall MAGL algorithm is shown in Figure 5. It should be noted that due to the position moves, partial spiders may move out of the specified range resulting in infeasible solutions. Once that happens, each illegal variable in solution will be reset by a random value in the encoding interval stipulated in Section IV-A.

V. NUMERICAL EXPERIMENTS A. EXPERIMENTAL PREPARATIONS
To demonstrate the algorithmic performance, a variety of numerical experiments are carried out in this section. Some necessary preparations for carrying out experiments are first introduced as follows. Benchmark problem. The studied DFAPFSP is a novel scheduling problem, and no benchmark is presented for it in    processing time of jobs and the assembly time of products are deterministic integers, which are randomly generated in interval [1,99] with uniform distribution.
Algorithmic calibration. For the proposed MAGL, there are four parameters that need to be tuned: Np, τ , Tim, and P min -P max . We calibrate them by Taguchi method of design of experiment (DOE) [44]. Each parameter is set four levels that are shown in Table 3. According to the number of parameters and the levels of each parameter, 16 different combinations of parameter levels are tested in DOE, see Table 4. For the instances to be used for calibration, the one with combination n/f/m/t/g = 25/2/4/8/2 is selected for small-size scenario. For each of 16 combinations of parameter levels, we carry out the calibrated algorithm 5 times independently for the chosen instance and the current combination. Then, the average makespan value among 5 runs is used for the first response variable (RV 1 ). For large-size scenario, we use the same scheme to gain the second response variable (RV 2 ) by the instance with combination n/f/m/t/g =100/6/6/20/4. The statistic results are listed in Table 5, where the ''Rank'' line shows the influence priority of the parameters on MAGL performance and the ''Value'' line provides the suggested parameter values. From Table 5, it can be concluded: 1) For small-size instances, the influence priority on MAGL is that τ ranks first, followed by Tim, P min -P max , and Np, and the suggested parameter values are Np = 30, τ = 0.04, Tim = 200, and P min -P max = 0.2-0.9; 2) For large-size instances, the influence priority on MAGL is that P min -P max ranks first, followed by Tim, τ and Np, and the suggested parameter values are Np = 40, τ = 0.02, Tim = 200, and P min -P max = 0.1-1.0. Other preparations. To prepare more convergence time for those testing instances with larger size, the termination criterion used in algorithmic calibration and comparison is the maximum elapsed runtime ρ× n × m × t milliseconds, where ρ is set at the value 10 for small-size instances and the value 3 for large-size instance. All algorithms in the experiments are coded in Matlab 2012a and implemented on the PC: Intel(R) Core(TM) i3-4150 CPU @ 3.50 GHz with 4.00 GB RAM in the 32-bit Windows 7 operating system. In order to measure the algorithmic performance, the relative percentage deviation (RPD) is applied and calculated by: where C max is the makespan of a schedule obtained by the given algorithm for an instance, and C best is the best known makespan found throughout the experiments. Due to inherent randomness, the computational results of metaheuristic methods may vary with different executions. Thus, instead of directly applying RPD as the performance metric, the metaheuristics are run 20 times independently for each instance. Then, the best RPD (bRPD) and the average RPD (aRPD) among 20 runs are used for the final performance metrics.

B. VALIDATION OF MILP MODEL AND MAGL ALGORITHM
To validate our MILP model and MAGL algorithm, we compare them against three algorithms [30]: TPHS, HVNS and HPSO. The three algorithms are originally proposed for the DFAPFSP but with setup constraints. So, we adapt them to the problem in this paper by the way of our makespan calculation. The MILP model and the TPHS method are deterministic, and so they are only implemented once for each instance to generate the metric RPD. Whereas for the algorithms HVNS, HPSO and MAGL, two metrics of bRPD and aRPD are calculated. The MILP model is solved by the CPLEX 12.2 with the time limits of 5 ×n×m×t seconds for small size instances and 3600 seconds for large size instances. If CPLEX cannot gain an optimal solution within the stipulated time, it will continue to be executed until finding a feasible solution (which may be not optimal). The computational results have been grouped in Tables 6 and 7 by problem parameters n, f , m, t and g, and the best metrics are marked in bold.
For small size instances, Table 6 shows that the proposed MILP model and the MAGL algorithm are significantly better than the compared three competitors TPHS, HVNS and HPSO for each instance group. It can be concluded that the proposed MILP model and MAGL algorithm are very effective when solving the DFAPFSP with small size. From Table 6, we can also see that the performance of the MAGL algorithm is very close to that of the MILP model   (i.e., exact method), which further demonstrates the effectiveness of MAGL algorithm. In addition, Table 6 shows that the RPD values of MILP model are slightly larger than zero, which implies that the optimal solution has not been found in some small instances. The fundamental reason is that the scheduling process in DFAPFSP is complex, which needs to involve four inter-related decisions. So, even using MILP model to solve small-size problems is still challenging. The direct reason is that a time limit is set for CPLEX tool when solving MILP model. If an optimal solution in the allowed time is not obtained, it will continue until finding a feasible solution (which may not be optimal). For large size instances, it can be seen from Table 7 that the MAGL algorithm significantly surpasses all competitors compared, which demonstrates its excellent ability to solve large-size problems. Whereas for the MILP model,  Table 7 shows that the mean value of RPD reaches 19.685 together with a time cost of 3600 seconds. This RPD value is much larger than that generated from the algorithms compared. It suggests that the MILP model is not effective when solving the DFAPFSP with large size. Furthermore, the above comparison results are analyzed by using the non-parametric Wilcoxon signed-rank test [45] with confidence level 0.05. It can be seen from Table 8 that these comparisons are statistically significant.

C. VALIDATION OF MAGL COMPONENTS
The MAGL proposed in this paper is a hybrid method consisting of different components. The problem-solving ability of the whole algorithm has been proved. Next, we analyze the effectiveness of its main components. For the purpose, the other four versions of MAGL algorithm are tested, i.e., MAGL G0 , MAGL G1 , MAGL SS and MAGL ML . MAGL G0 is MAGL without any local optimizers and with the original mating operation presented in [40]. MAGL G1 is MAGL G0 with the improved mating operator presented in this paper. MAGL SS is MAGL only with local search SLS, whereas MAGL ML is MAGL only with local search MLS. We compare them with TPHS and MAGL by 12 differentsize instances from the proposed benchmark problems. The comparison results are shown in Table 9.
From Table 9, it can be seen that MAGL G1 outperforms TPHS significantly for each instance and each performance metric, which proves that the global exploration optimizer we proposed is very effective. Also, MAGL G1 gets smaller bRPD and aRPD values than MAGL G0 , demonstrating the effectiveness of our improved mating operator. With the application of local search strategies, the performance of MAGL SS and MAGL ML surpasses those of TPHS, MAGL G0 and MAGL G1 at a considerable margin. It suggests that the devised two local search strategies are of great significance to the MAGL algorithm. Furthermore, from the comparison of MAGL SS and MAGL ML , we can draw a conclusion that the local search MLS is more critical than the local search SLS. As expected, MAGL algorithm performs best among all the candidates being compared, the success of which is mainly due to the effectiveness of the MAGL components and their organic hybridization in MAGL. To be specific, the global exploration optimizer is efficient due to the well-designed mechanism: bi-population co-evolution and cross-population interaction. The local search optimizer based on meta-Lamarckian learning strategy can effectively enhance the exploitation in the promising region, while the simplex search-based local search optimizer is helpful to refine the population's overall quality. The hybridization of the global and local search optimizers within MA framework can well balance the exploration and exploitation abilities. Besides, Table 10 shows that the above comparison results are statistically significant.

VI. CONCLUSIONS
This paper studies a novel distributed flexible assembly permutation flowshop scheduling problem (DFAPFSP) to minimize makespan criterion. It is a practical and promising research topic due to the successful applications in modern manufacturing. The configuration consists of the distributed flowshop factories for job-processing at the first stage and the flexible assembly machines for product-assembling at the second stage. To solve the DFAPFSP, we first construct a MILP model, and then put forward a memetic algorithm (MAGL) that has been proved to be very effective based on a comprehensive computational campaign. For the MAGL, the contributions includes four aspects: a) Present a global exploration optimizer with bi-population co-evolution and cross-population interaction; b) Present a meta-Lamarckian learning-based local exploitation optimizer to improve the best solution in population; c) Present a simplex search-based local exploitation optimizer to refine the population's overall quality; and d) Formulate a novel and efficient MA for the DFAPFSP based on the encoding scheme designed by introducing problem-specific knowledge.
For future research, some aspects deserve our attention. The proposed MA algorithm is efficient to cope with the DFAPFSP, which inspires us to apply it to more complex scenarios such as restricted buffers, batch processing, multiobjective, heterogeneous factories, and so on. It is also very meaning to design more sophisticated global and local optimizers based on problem related knowledge and other excellent optimization technologies.