Solving multi-objective Modified Distributed Parallel Machine and Assembly Scheduling Problem (MDPMASP) with eligibility constraints using metaheuristics

ABSTRACT

We present a new generalization and metaheuristics for the distributed parallel machine and assembly scheduling problem (DPMASP), namely MDPMASP with eligibility constraints. There is a set of unrelated factories or production lines. Each has a series of non-identical parallel machines with varying processing speeds. Then, it was disposed of as a single assembly machine in a series. Each product is assembled from a set of components (jobs). Each item necessitates multiple unidentical jobs. The objectives are to minimize mean flow time and the number of tardy jobs. We suggest four basic heuristics and three metaheuristics to tackle the problem. The Taguchi approach is used to discuss and calibrate various metaheuristic parameters. Algorithms are compared using the four performance measures. The computational results show that the proposed algorithms can solve moderate-sized problems efficiently and near-optimal solutions in most cases. Moreover, based on the four performance measures, MGA is the best method, followed by MSA, MPSO, SH2, SH4, SH1 and SH3.

KEYWORDS:

• MDPMASP
• eligibility constraints
• multi-objective problem
• mean flow time
• the number of tardy
• heuristic
• metaheuristic

1. Introduction

The global economy has undergone radical changes in the last two decades. This is influenced by the growth of global trade and explosive international competition. Globalization has encouraged natural selection, and the strongest will survive. Market success will be obtained by companies that are able to respond quickly to customer needs and desires. One way to respond to these conditions is to increase the variety of products offered. A diversified product portfolio and diversified offering are the company’s main assets.

Providing multiple attributes in a product is an imt5 portant way to attract customers, but it often increases the complexity of planning and scheduling processes. In the automotive industry sector, where there is variation (variation offered to customers) and internal variation (variation involved in making the product). Forecasting and dealing with internal variations is a challenge for companies to build products to order. The effectiveness of strategies to mitigate negative effects, such as modularity, mutability, late configuration, and option bundling, depends on the order fulfillment strategy the company follows. Several kinds of manufacturing production strategies such as: make to order or make to assembly can reduce these negative effects.

Suppliers play a very important role at every stage of the product life cycle. From sourcing raw materials to helping increase production, and to finding better raw material options when the market gets saturated, companies need to work closely with their suppliers to get the best out of their products. Automotive suppliers are defined as companies that produce goods that are used in the production process of cars or become parts of cars, so they supply these goods directly or indirectly to car manufacturers, which is identical to Original Equipment Manufacturer (OEM). These items can be individual components, such as screws, or whole assemblies, such as pre-assembled door modules. Therefore, automotive suppliers are part of the automotive supply chain.

Production scheduling has become a must for manufacturing operations looking to take their production facilities to the next level. Production scheduling refers to the allocation of resources, operations, and processes needed to make goods and services. The purpose of a production schedule is to balance customer requirements with available resources while operating in a cost-effective manner. If the production schedule is inaccurate and not feasible with the available resources, we will have problems getting the goods produced and delivered to the customer on time. In the automotive industry sector, production scheduling is crucial and difficult to do due to the high variety of products and components, each component must be produced by a different supplier. However, it is possible that a supplier is able to produce more than one component. Therefore, a scheduling problem arises called the Distributed Parallel Machine and Assembly Scheduling Problem (DPMASP) with eligibility constraints (Sara Hatami, Ruiz, & Romano, 2015).

The DPMASP with eligibility constraints, as first described by Hatami et al. (2015), is to optimize the processing of n jobs in a series of two manufacturing stages, namely, the production stage and assembly stage. The production stage consists of a set F of f identical factories with a set M of m unrelated parallel machines at one factory. Besides, the eligibility constraints are considered in this problem. The eligibility constraints indicate some conditions where job j can only be assigned to an eligible factory. In the second stage, a single assembly machine is assigned to assemble all the final products. All factories are identical and have the same number of machines.

This article discusses an expansion of the DPMASP to an unrelated multi-factory manufacturing with an unidentical parallel machine scenario. A set of products were manufactured through an assembly of several jobs, each fabricated in the PFSP factories of the production stage. A product could be assembled at the single machine assembly stage only after all the jobs in the production stage are completed. Moreover, each factory has several unidentical parallel machines with different specifications and capacities that influence processing speed. Therefore, the assumption of an identical parallel machine in DPMASP is no longer valid. In conclusion, to make the DPMASP more realistic to implement on a more complex system, some assumptions about identical parallel machines must be changed.

The previous study on the distributed assembly scheduling problem only considers the optimization of a single objective, i.e. minimizing the makespan (Duan et al., 2016; S Hatami et al., 2015; Pan, Gao, Xin-Yu et al., 2019; Shao et al., 2019; K. Wang et al., 2017). Thus, the concept of optimality for a single objective is meaningless since the best solution for one objective may not be the best for another (Glover & Sörensen, 2015). Besides, most studies present MILP and local search metaheuristics to solve the distributed assembly flow shop scheduling problem and its variants (Framinan et al., 2019). Local search metaheuristics iteratively make small changes to a single solution. Constructive metaheuristics construct solutions from their constituting parts. Population-based metaheuristics iteratively combine solutions into new ones. However, these classes are not mutually exclusive, and many metaheuristic algorithms combine ideas from different classes. Such methods are called hybrid metaheuristics (Glover & Sörensen, 2015).

As a result, we analyze the Modified Distributed Parallel Machine and Assembly Flowshop Scheduling Problem, abbreviated as MDPMASP, with eligibility constraints. Two manufacturing stages are considered: production and assembly. There is a set of F distributed unrelated factories with a set M of m unidentical parallel machines in each factory as a shop configuration. In contrast, the assembly stage consists of a single assembly machine referred to as M_A. There are four interrelated decisions to be made to solve this problem. These decisions generate an assembly product sequence as the initial solution, complete the factory selection for each job, schedule within each factory, and assign jobs to one of the machines set at the specific factory. Additional details will be described further. In addition, we consider the mean flow time and the number of tardy products as optimal criteria. To our best knowledge, this is the first attempt at resolving such a critical scheduling issue. The main contributions are constructing a new multi-objective scheduling problem and developing several algorithms that were often implemented in simple cases previously. As regards the computational complexity, we may deduce that the Multi-Objective MDPMASP with eligibility constraints is an NP-hard problem because the standard parallel machine problem (even when two identical machines are used, i.e. the P2/Cmax problem) is already NP-Hard according to Lenstra et al. (1977).

In this paper, we proposed four simple heuristic algorithms and three high-performing population-based algorithms such as genetic algorithm (GA), simulated annealing (SA), and particle swarm optimization (PSO) to solve the multi-objective MDPMASP with eligibility constraints. One of the population-based approaches’ advantages was being able to explore multiple regions of the fitness landscape concurrently, thus mitigating search bias in the initial position or sampling process. They also provide a good balance between exploration and exploitation in the search landscape (Prugel-Bennett, 2010). In addition, population-based algorithms allow obtaining better (feasible) solutions that do not dominate each other in a shorter time than local search metaheuristics. Although it has a slightly more complicated procedure. Significant problem-specific knowledge and accelerations are investigated in order to improve the algorithms’ effectiveness. The small size problems have been solved using MILP and total enumeration sampling to validate the proposed algorithms. Four performance metrics, i.e. overall non-dominated vector generation (ONVG), overall number of non-dominated solutions (ONSN), average quality (AQ), and distance metric (DI_R) proposed by Qian et al. (2009) are used to measure the algorithms’ performance. It was done because the multi-objective problems result in multiple non-dominated solutions.

The rest of our paper is organized as follows: Section 2 is a literature review, and the problem definition is explained in Section 3. In Section 4, we present a Mixed Integer Linear Programming (MILP) model to solve the considered problem. In Section 5, four simple heuristics and three metaheuristic algorithms are proposed to solve the problem. Section 6 presents a comprehensive computational evaluation of the proposed algorithms. Finally, some concluding remarks and future research directions are provided in Section 9, which shows the computational results.

2. Literature review

Due to the present economic trend toward customized manufacturing, assembly production is becoming more prevalent. Typically, an assembly floor is divided into two stages: production and assembly (Hatami et al., 2015). Assembly Scheduling Problems (ASPs) are gaining increasing attention from academics and practitioners as a subset of scheduling decisions on which parts/components/subsets of products or services must be manufactured in parallel and then assembled in a final stage (Framinan & Perez-Gonzalez, 2017). Among ASPs, the Assembly Flowshop Scheduling Problem (AFSP) is a prevalent type that occurs in a wide variety of production systems (S.-Y. Wang & Wang, 2015). The layout for this problem includes a flowshop for part processing and a single assembly machine for product processing. Lee et al. (1993) were the first to solve the AFSP problem, followed by Yokoyama (2001), Yokoyama and Santos (2005), Sung and Juhn (2009), and Sheikh et al. (2018), and Framinan et al. (2019), and numerous others. AFSP is used in a variety of industries, including engine assembly (Lee et al., 1993) and personal computer manufacture (Hariri & Potts, 1997), to mention a few.

Numerous difficulties confront the industrial sector in reconciling manufacturing complexity and customer desires. In an uncertain and volatile global market, several organizations must quickly grow their product portfolio while maintaining a shorter product life cycle (Elmaraghy et al., 2013). To address these difficulties, diverse firms are required to collaborate to assist in manufacturing and supply chain activities. Unlike the AFSP assumption that all goods and components are processed in a single manufacturing facility. Due to the need for distributed manufacturing in today’s fractured and globalized economy, large corporations frequently operate many production centers or factories (Moon et al., 2002). Managers must decide how to distribute jobs to manufacturers and how to properly schedule jobs inside each facility to remain competitive in continually changing markets. Thus, the so-called Distributed Scheduling Problem (DSP), which concerns the assignment and scheduling of jobs across multiple factories, DSP has become a popular study area (Fernandez-Viagas & Framinan, 2015; Bahman Bahman Naderi & Ruiz, 2014). Collaboration with other factories has resulted in numerous benefits, including increased adaptability and resilience, increased product variety and quality, decreased manufacturing costs, and less risk.

For instance, B. Naderi and Ruiz (2010) study scheduling issues in a distributed permutation flowshop (DFSP) composed of a collection of similar factories, each of which is a permutation flowshop. Ruiz et al. (2019) describe an iterated greedy algorithm for the DFSP that greatly improves results. (Pan, Gao, Xin-Yu, et al., 2019) discuss the DFSP with the total flowtime requirement and offer three successful constructive heuristics as well as four meta-heuristics. Li et al. (2019) offer an improved artificial bee colony algorithm in a prefabricated system for the DFSP with distance coefficient. Hatami et al. (2013) propose a hybrid DFSP-assembly line with a distributed flowshop as the first stage and a single assembly machine as the second stage. Xiong et al. (2014) use distributed manufacturing to solve a two-stage assembly flowshop problem.

To the extent of our knowledge, the DAPFSP was introduced by Hatami et al. (2013) for the first time, taking into the eligibility constraints. In DAPFSP, there are two different stages. The first stage consists of a set of identical factories (all factories are capable of processing all components and each consists of a permutation flowshop). The second stage entails the establishment of an assembly facility. The assembly factory assembles all of the components into finished goods. The idea is to keep the makespan as short as possible. To solve this problem, Hatami et al. (2013) proposed an MILP model, three constructive heuristics, and a VND algorithm. Zhang et al. (2020) construct a genetic algorithm (GA) with an enhanced crossover strategy and three local search schemas that outperform Hatami’s constructive heuristics. Later, an estimation of the distribution algorithm-based memetic algorithm (EDA-MA) with a bi-vector-based representation was investigated by Wang and Wang (2015) to improve the potential solution. Lin and Zhang (2016) present a hybrid biogeography-based optimization (HBBO) algorithm that uses a path-relinking heuristic as an approach for product local search strategy and an insertion-based heuristic to optimize part permutations. Lin et al. (2017) propose a Backtracking Search Hyper-Heuristic (BSHH) consisting of 10 low-level heuristics to perform operations on the solution space and a backtracking search to modify the heuristics at the lowest level. In addition, Hatami et al. (2015) present the DPMASP with eligibility constraints. DPMASP consists of the distributed unrelated factories in the first stage and the single assembly machine in the second stage. In the first stage, there is a set of unrelated factories with identical parallel machines in the production stage. In this case, parallel machines within a factory have the same specifications and production capacity. It is shown by the same job processing times on parallel machines at the specific factory. Jobs have to be assigned to factories and to machines. Finished jobs in the first stage are assembled into final products in the second stage. Due to eligible constraints, machines cannot be idle, and some jobs can be processed only in certain factories. To solve this problem, Hatami et al. (2015) proposed a mathematical model and two high-performing heuristics to minimize the makespan in the assembly stage. There are three assumptions made in DPMASP with eligibility constraints. First, the assembly stage is separated from the assembly operation at each factory. Second, allow the different jobs that compose a product to be produced in different factories. Third, each product might have a number of jobs (components) different from m. The problem studied by Hatami et al. (2015) is different from our proposed problem, i.e. MDPMASP with eligibility constraints. Here, we consider the use of non-identical parallel machines with varying specifications and production capacity. This increases the problem’s complexity and difficulty in solving it.

To the best of our knowledge, no further literature exists on DPMASP with eligibility constraints since 2019, so this is the first effort that considers the non-identical parallel machine in a distributed manufacturing setting. Moreover, we present four simple heuristics and three metaheuristic algorithms to minimize the mean flow time and number of tardy products in MDPMASP with eligibility constraints.

3. Problem definition

The modified distributed parallel machine and assembly scheduling problem (MDPMASP) with eligibility constraints is now described in detail. This problem consists of two manufacturing stages: production and assembly. Both stages are assigned to produce a set T of t different products manufactured through an assembly process of a set N of n jobs in series according to the requested demand during the planning horizon. There is a set F of f distributed non-identical factories with a set M of m unrelated – parallel machines in each one. In the assembly stage, there is a single assembly facility denoted as MA. Furthermore, eligibility constraints are considered and represented by LF_j, LF_j ⊆ F, which is the subset of factories where job j can be assigned, where f ≥ | LF_j | ≥ 1, j = 1, … , n., and job j can only be processed in the eligible factory.

The related job is assembled into a product by the assembly program of each product N_s, N_s ⊆ N. Therefore, product s consists of |N_s| jobs. Each job is associated with a single assembly program for a specific product and ∑|N_s| = n, where s = 1, 2, … , t. A product can be assembled at a single assembly machine after all jobs in N_s have been completed in the first stage. For the first stage, PT_jqi denotes the production time of the job j at a machine i of factory q. The second stage is to set AT_s to be the assembly time of product s. Note that all factories are non-identical and have the same number of machines. All processing times are positive, known integer quantities, and deterministic or random within a certain range.

Each product requires two or more jobs. Jobs are different in kinds and batch sizes. At first, all jobs are processed in the production stage, including several unrelated parallel machines with different specifications, capacity, and speed for processing. All machines in the specific factory can process all jobs according to LFj, and each machine can also process jobs in the specific factory. Still, they would be processed by exactly one machine at a time. Subsequently, all jobs are disposed to the second stage and are processed by the MA machine.

A simple problem is used to illustrate the proposed algorithms. It consists of fifteen products (t = 15), six jobs (n = 6), and three different factories (f = 3) with two unrelated parallel machines (m = 2) in each factory. Each machine in the factory runs at a different speed, resulting in differences in the completion time of jobs on a machine. Every product s is composed of several jobs that exist in the related assembly program Ns. The assembly program of the six products are N₁ = {9,10}, N₂ = {7,12,13,14}, N₃ = {8,15}, N₄ = {6}, N₅ = {5}, and N₆ = {1,2,3,4,11}, meaning job 9 and job 10 must be finished in order to assemble product 1.

There are three kinds of sequences: a product sequence (π), a job to factory sequence (π_Fq, ∀q), and a job to factory machine sequence (π_FMqi, ∀q, i). In this example, there are several data sets of each sequence, π = {2,6,1,4,5,3}, π_F1 = {14,12,13,1,4,8}; π_F2 = {11,2,6,5,15}; π_F3 = {7,3,10,9}, and π_FM11 = {14,1,8}; π_FM12 = {12,13,4}; π_FM21 = {11,5,15}; π_FM22 = {2,6}; π_FM31 = {7,9}; π_FM32 = {3,10}. While all jobs in the assembly program of a product s are completed in the production stage, it can be assembled in the assembly stage. The earliest assembling time of product s is denoted as E_s. Figure 1 shows the problem with three factories, two machines, 15 jobs, and a single assembly machine.

Figure 1. Representative Solution of MDPMASP with eligibility constraints.

The purpose of MDPMASP with eligibility constraints is to assign jobs to one machine at a specific factory, schedule all assigned jobs in each machine at each factory, and schedule all products in the single assembly machine while minimizing the mean flow time and the number of tardy products simultaneously. Due to the computational complexity of the problem, it is an NP-Hard problem. The MDPMASP with eligibility constraints is an essential generalization of existing problems and has not been studied before to the best of our knowledge,

(1) $F=\frac{1}{t}\sum _{s=1}^t{F}_s=\frac{1}{t}\sum _{s=1}^t{E}_s+A{T}_s$(1)

(2) $E_s=\underset{j^{\ast }\in N_s}{max}\{C_{j^{\ast }}\}$(2)

where $\bar{F}$ is the mean flow time; $E_s$ is the maximum completion time of job j* (C_j*, j*∈N_s), and it is equal to the starting time of product s in an assembly stage; $A{T}_s$ is the time of product $s$ at the assembly stage. The number of tardy products ($n_T$) can be calculated by

$n_T={\sum }_{s=1}^t{n}_T^{s}where{n}_T^s=\left\{\begin{array}{c}1,if(c_s-D{D}_s)>0\\ 0,otherwiese\end{array}\right.$(3)

$n_T$ denotes the number of tardy products; $C_s$ denotes the maximum completion time of a product $s$; $D{D}_s$ denotes the product’s due date. Four simple heuristics and three metaheuristic algorithms are proposed to find efficient schedules that minimize the mean flow time and the number of tardy products on the MDPMASP with eligibility constraints system.

The mathematical model for this problem was developed based on (Sara Hatami et al., 2015). The distinction was that, in this model, the factories are not identical, and the machines in each factory are unrelated. The machines in a factory have different processing times. As a result, we need to add an index q to the production time variable $P{T}_{jk}$, so that it becomes $P{T}_{jkq}$. A mixed-integer linear programming model for the Multi-Objective Modified Distributed Parallel Machine and Assembly Scheduling Problem (MDPMASP) with eligibility constraints is explained as follows. This model was developed based on

Table

Index	Description
$i,j$	Jobs, $i,j=0,1,\dots ,n$, where $0$ is a dummy job
$k$	Machines, $k=1,\dots ,m$
$q$	Factories, $q=1,\dots ,f$
$l,s$	Products, $l,s=0,1,\dots ,t$, where $0$ is a dummy product
$M$	A sufficiently large positive number

Table

Parameter	Description
$n$	Number of jobs
$m$	Number of machines
$f$	Number of factories
$t$	Number of products
$P{T}_{jqk}$	Production time of job $j$ on machine $k$ in factory $q$
$A{T}_s$	Assembly time of product $s$ at the assembly stage
$D{D}_s$	Due date of product $s$ at the assembly stage
$G_{js}$	Binary parameter equal to $1$ if job $j$ belongs to product $s$, and 0 otherwise
$L{F}_j$	A subset of factories where job $j$ can be assigned

Table

Variable	Description
$X_{ijkq}$	Binary variable equal to 1 if job $i$ is an immediate predecessor of job $j$ on machine $k$ in factory $q$, and 0 otherwise, where $j\ne k$
$Y_{ls}$	Binary variable equal to 1 if product $l$ is an immediate predecessor of product $s$ at the assembly stage, and 0 otherwise, where $s\ne l$
$C{P}_j$	Completion time of job $j$ at the production stage
$C{A}_s$	Completion time of product $s$ at the assembly stage
$\bar{F}$	Mean flow time
$n_T$	Number of tardy product

The objective function of the model is to minimize the mean flow time and the number of tardy products, respectively,

(4) $Min\bar{F}\hspace{0.17em}$(4)

(5) $\bar{F}=\frac{1}{t}\sum _{s=1}^t{F}_s$(5)

(6) $F_s=C{A}_s=\underset{1\le j\le n}{max}\{C{P}_j.G_{js}\}+A{T}_s$(6)

(7) $C{P}_j=C{P}_i+P{T}_{jqk}.X_{ijqk},\mathrm{\forall }i,j\ne i,k,q\in L{F}_i{\cap }^{\hspace{0.17em}}L{F}_j$(7)

(8) $Min{n}_T$(8)

(9) $n_T=\sum _{s=1}^n\delta (T_s),\delta (x)=1ifx>0,and\delta (x)=10otherwise$(9)

(10) $T_s=\underset{1\le s\le t}{max}\{0,L_s\}$(10)

(11) $L_s=C{A}_s-D{D}_s$(11)

(12) $C{A}_s=\underset{1\le j\le n}{max}\{C{P}_j.G_{js}\}+A{T}_s$(12)

(13) $C{P}_j=C{P}_i+P{T}_{jqk}.X_{ijqk},\mathrm{\forall }i,j\ne i,q,k\in L{F}_i{\cap }^{\hspace{0.17em}}L{F}_j.$(13)

The constraints of the model are as follows:

1. Control and ensure that every job must be processed on exactly on one machine at a certain factory and must have exactly one preceding and succeeding job,

(14) $\sum _{i=0,i\ne j}^n\sum _{i=1}^m\sum _{q=1}^f{X}_{ijqk}=1\mathrm{\forall }j,$(14)

(15) $\sum _{j=0,j\ne i}^n\sum _{i=1}^m\sum _{q=1}^f{X}_{ijqk}=1\mathrm{\forall }i.$(15)

• (2) State that each machine at each factory has to have a dummy job 0 as a predecessor and successor, respectively.

(16) ${\sum }_{\begin{array}{c}j=1\\ q\in L{F}_j\end{array}}^n{X}_{0jqk}\mathrm{\forall }q,k,$(16)

(17) ${\sum }_{\begin{array}{c}i=1\\ q\in L{F}_i\end{array}}^n{X}_{i0qk}\mathrm{\forall }q,k.$(17)

• (3) Determine that if a job is sequenced on a machine, then its predecessor and successor must be processed on the same machine.

(18) ${\sum }_{\begin{array}{c}j=1,j\ne i\\ q\in L{F}_i\end{array}}^n(X_{ijqk}-X_{jiqk})=0\mathrm{\forall }i,q,k.$(18)

• (4) Avoid the occurrence of cross-precedences, which means that a job cannot be both a predecessor and a successor of another job at the same time,

(19) $\sum _{k=1}^m\sum _{\begin{array}{c}q=1\\ q\in L{F}_i\\ q\in L{F}_j\end{array}}^f{X}_{ijqk}+\sum _{k=1}^m\sum _{\begin{array}{c}q=1\\ q\in L{F}_i\\ q\in L{F}_j\end{array}}^f{X}_{jiqk}\le 1\mathrm{\forall }i\in \{1,\dots ,n-1\},j>1.$(19)

• (5) Specify that if a job $j$ is placed immediately after job $i$, its processing on the machine $k$ cannot start before the processing of the job $i$ on the machine $k$ finishes,

(20) $C{P}_j\ge C{P}_i+P{T}_{jqk}+M(1-X_{ijqk})\mathrm{\forall }i,j,q,k\in L{F}_i{\cap }^{\hspace{0.17em}}L{F}_j.$(20)

• (6) In the assembly stage, each product must have one predecessor and not more than one succeeding product,

(21) $\sum _{s=0,s\ne l}^t{Y}_{ls}=1\mathrm{\forall }s$(21)

(22) $\sum _{s=0,s\ne l}^t{Y}_{ls}=1\mathrm{\forall }l.$(22)

• (7) Prevent a product from being both a predecessor and a successor to another product at the same time,

(23) $Y_{ls}+Y_{sl}\le 1\mathrm{\forall }l\in \{1,\dots ,t-1\},s>1.$(23)

• (8) Imply that each product cannot begin to assemble before all its jobs in its assembly program are completed in the first stage,

(24) $C{A}_s\ge (C_j.G_{js})+A{T}_s\mathrm{\forall }j,s.$(24)

• (9) Determine that if product $s$ is placed immediately after product $l$, it cannot start to be assembled on the assembly machine before the assembling of product $l$ in the assembly machine has finished,

(25) $C{A}_s\ge C{A}_l+A{T}_s+M(1-Y_{ls})\mathrm{\forall }l,s\ne l.$(25)

• (10) Define the makespan and the domain of decision variables, respectively,

(26) $C_{max}\ge C{A}_s\mathrm{\forall }s$(26)

(27) $X_{ijqk}\in \{0,1\}\mathrm{\forall }i,j\ne i,q,k\in L{F}_i{\cap }^{\hspace{0.17em}}L{F}_j$(27)

(28) $Y_{ls}\in \{0,1\}\mathrm{\forall }l,s\ne l$(28)

(29) $C{P}_j\ge 0\mathrm{\forall }j$(29)

(30) $C{A}_s\ge 0\mathrm{\forall }s.$(30)

Note that $C{P}_0=C{A}_0=0$.

4. Proposed algorithms

MDPMASP with eligibility constraints is an NP-hard problem If $n>\left(f\times m\right)$. To solve large problems, a heuristic and metaheuristic-based solution search technique is required.

4.1. Simple heuristic algorithms

Heuristic methods can be used to find solutions to optimization problems, but they do not guarantee that they always produce optimal solutions. In principle, the heuristic method is used to obtain the best possible solution in a required time span. Heuristic methods can be used to find solutions to optimization problems, but they do not guarantee that they always produce optimal solutions. In principle, the heuristic method is used to obtain the best possible solution in the required time span. We present four simple heuristic algorithms, namely SH-1, SH-2, SH-3, and SH-4. The four algorithms will produce three solution structures, i.e. product sequence (π), job to factory sequences (π_Fq, ∀q), and job to factory-machine sequence (π_FMqk, ∀q, k). The goal to be achieved is the simultaneous minimization of the mean flow time and the number of tardy products. In general, a simple heuristic algorithm was developed based on two algorithms developed by Erenay et al. (2010), namely, independent beam search (BS-I) and dependent beam search (BS-D), and the two algorithms presented by B. Naderi and Ruiz (2010), namely, NR₁ and NR₂. The BS-I and BS-D algorithms are proven to be efficient in minimizing the mean flow time and number of tardy jobs simultaneously on a single machine. Meanwhile, NR₁ and NR₂ proved efficient in minimizing the makespan in the distributed permutation flow shop scheduling problem (DPFSP). Table 1 shows the structure of the proposed simple heuristic algorithms.

Table 1. The proposed simple heuristics

Proposed Simple Heuristic (SH)	Product Sequence (π)	Job to factory sequences (πFq,∀q)	Job to factory-machine sequences (πFMqk,∀q,k)
SH-1	BS-I Algorithm	Proposed job to the factory assignment procedure	NR1
SH-2	BS-I Algorithm		NR2
SH-3	BS-D Algorithm		NR1
SH-4	BS-D Algorithm		NR2

The following is the general procedure of the proposed simple heuristic algorithm: First, the BS-I or BS-D algorithm presented by Erenay et al. (2010) is used to determine the best product sequence (π) based on AT_s and DD_s. This can reduce the average flow time and the number of late products in a hierarchical machine. After identifying an efficient product sequence, focus should be given to minimizing E_s. The minimum E_s can be achieved by forming a job to factory sequence and a job to factory-machine sequence with the lowest makespan. Second, weights representing the priority level of each job are required to determine the job to factory sequences (π_Fq, ∀q) for each product sequence (π). Third, the job to factory-machine sequence (π_FMqk, ∀q, k) is then built using two job to factory assignment rules described in B. Naderi & Ruiz (2010). The first one, referred to as NR₁, assigns job j to the machine with the lowest current Cmax, without considering job j. The second rule, NR₂, assigns job j to the factory with the lowest Cmax after scheduling job j. Fourth, each product’s earliest assembly time (E_s), completion time (C_s), and flow time (F_s) are known. Fifth, the mean flow time and the number of tardy products can be computed. The process is applied until all jobs in the assembly program of product s are scheduled. The pseudocode of proposed simple heuristic algorithms is shown in Figure 2.The processing time and weight calculation are listed in the Table 2.

Figure 2. Pseudocode of the simple heuristic algorithms.

Table 2. Weight calculation for the simple example

j*	PTj*qk	wj*			min{wj*}
7 12 13 14	0 0 0 0 7 6 9 7 0 0 8 9 5 5 8 7 0 0 10 9 0 0 0 0	0 8.0000 5.0000 9.5000	0 0 7.5000 0	6.5000 8.5000 0 0	6.5000 8.0000 5.0000 9.5000

4.2. Metaheuristic algorithms

4.2.1. The Modified Genetic Algorithm (MGA)

Genetic algorithm (GA) belongs to the class of evolutionary algorithms based on a population of individuals. GA are numerical optimization algorithms inspired by both natural selection and natural genetics. This method is a general one and capable of being applied to an extremely wide range of problems. GA generally consists of the following three operators: selection, crossover, and mutation (chosen in part by analogy with the natural world) to direct population (over a series of time steps or generations) towards convergence at the global optimum (Santosa & Willy, 2011).

The modified genetic algorithm (MGA) tries to find efficient schedules by minimizing the mean flow time and the number of tardy products simultaneously using Pareto’s dominant concept. After the objective values are found, an efficient schedule is determined using the non-dominated sorting and crowding distance procedure (Deb et al., 2002). The following is the pseudocode of the proposed MGA (see Figure 3).

Figure 3. Pseudocode of the Modified Genetic Algorithm (MGA).

4.2.2. The Modified Simulated Annealing (MSA)

Simulated annealing is a generalization of a Monte Carlo method for statistically finding the global optimum between different local optima. Simulated Annealing (SA) is more powerful than simple local search by accepting worse solutions with a certain probability. There are several steps that the simulated annealing algorithm go through while a specified rate decreases the temperature. The annealing process runs through several iterations at each temperature to sample the search space. It generates a new solution at each cycle by applying a random change to the current solution (Santosa & Willy, 2011). In Simulated Annealing, a system is initialized at a temperature T with some configuration whose energy is evaluated to be E.

A new configuration is constructed by applying a random chance, and the change in energy ∆E is computed. The new configuration is unconditionally accepted if it lowers the energy of the system. If the change increases the system’s energy, the new configuration is accepted with some random probability. In the original Metropolis scheme (Santosa & Willy, 2011), the probability is given by the Boltzmann factor exp(-∆E/kT). This process is repeated sufficient times at the current temperature to sample the search space, and then the temperature is decreased. The process is repeated at successively lower temperatures until a frozen state is achieved. This procedure allows the system to move to lower energy states while still jumping out of local minima (especially at higher temperatures) due to the probabilistic acceptance of some upward moves. The following is the pseudocode of the proposed MSA (see Figure 4).

Figure 4. Pseudocode of the Modified Simulated Annealing (MSA).

4.2.3. The Modified Particle Swarm Optimization (MPSO)

Particle swarm optimization (PSO) was developed based on the behavior of a swarm of insects, such as a flock of birds or a school of fish. The word particle denotes, for example, a bee in a colony or a bird in a flock. Each individual or particle in a swarm behaves in a distributed way using its own intelligence and the collective or group intelligence of the swarm. As such, if one particle discovers a good path, the rest of the swarm will also be able to follow the good path instantly even if their location is far away in the swarm. In the context of multivariable optimization, the swarm is assumed to be of specified or fixed size with each particle located initially at random locations and velocity in the multidimensional design space. Each particle wanders around in the design space and remembers the best position (in terms of the food source or objective function value) it has discovered. Each particle keeps track of its coordinates in the problem space, which are associated with the best solution (fitness) it has achieved so far. This value is called Lbest. When a particle takes all the population as its topological neighbors, the best value is a global best and is called Gbest. At each time step, changing the velocity of (accelerating) each particle toward its pbest and lbest locations (local version of PSO). Acceleration is weighted by a random term, with separate random numbers being generated for acceleration toward Lbest and Gbest locations. In the past several years, PSO has been successfully applied in many research and application areas. It is demonstrated that PSO gets better results in a faster and cheaper manner compared to other methods. The following is the pseudocode of the proposed MPSO (see Figure 5).

Figure 5. Pseudocode of the Modified Particle Swarm Optimization (MPSO).

4.3. Proposed the job to factory procedure

The job to factory procedure is used to assign each job to one of the eligible factories and to define the job priority scale by considering the job sequence (π) and assembly program of the related product. This procedure worked in the following manner. As the example, for π = {2,6,1,4,5,3}, start from the first product in π (s= 2). First, find j*, where j*∈ N₂. Second, for all j*, determine the weight of each factory (w_j*qi, ∀q, i) by an average value of PT_j*qk, where k = 1 to m. Third, find a factory that has the least weight (except: w_j* = 0), then set it as the chosen factory. The last step is to define a job to factory sequence (π_Fq, ∀q), so j* has to be sorted by descending order of min{w_j*). Repeat the same steps until all jobs for all products are scheduled. The results are: π_F1 = {14,12,13,1,4,8}; π_F2 = {11,2,6,5,15}; π_F3 = {7,3,10,9}.

5. Result and discussion

In general, this research is conducted in three stages: 1) verifying and validating the proposed algorithm; 2) identifying the optimal parameter configuration for the proposed algorithm; and 3) measuring, comparing, and testing the proposed algorithm’s performance.

All proposed algorithms are coded in MATLAB R2014a and run on a computer with Intel Core i3-6006 U CPU at 2.0 GHz, 8Gb RAM, and 1TB HDD in Windows 7 system. Meanwhile, the solution of the mathematical model using an exact algorithm will be obtained by using LINGO solver. MDPMASP problem is an NP-hard optimization problem. LINGO software is not able to solve large-scale problems in an acceptable time. To get a fair performance comparison between the exact and proposed algorithms, the elapsed time to solve each instance is equal to 60 minutes.

5.1. Verification and validation

A model that is built must be credible and representative. It is used to ensure the reproducibility of measurement results in real conditions. In software testing, there are two terms: verification and validation. Verification is the process of determining whether algorithms and programs run as intended by repeatedly running the program. The verification ensures that the conceptual model (pseudocode and assumptions) is translated correctly into the programming language (Law & Kelton, 2000). Validation is the process of determining whether the stated model accurately represents the system (Hoover & Perry, 1989). Validation is the process of determining whether a conceptual model (as opposed to a computer program) accurately represents the modeled system (Law & Kelton, 2000).

In this study, verification and validation were carried out by performing manual calculations for small cases. Verification is carried out to ensure whether the MDPMASP with eligibility constraints has been represented in the real system properly, whether all mathematical operations and relations are correct, whether there are errors in coding, and whether the structure of the solution search for each algorithm is systematic and logical. While the validation is done by reviewing the structure of the solution and the proximity of the solution generated by the proposed algorithms to the optimal solution. The optimal solution is obtained by running the MILP model and the total enumeration technique.

For example, MDMASP with eligibility constraints with performance criteria $\bar{F}$ and $n_T$, consists of three products ($t=3$), eight jobs ($n=8$) and two unrelated distributed factories ($f=2$) with two unidentical parallel machines in each factory ($m=2$). Then by performing the procedure for generating AT_s and DD_s randomly at the intervals in Table 3, the assembly time and lead time for each product are obtained (see, Table 4). Every product s is composed of several jobs that exist in the related assembly program N_s. The assembly program of the three products are N₁ = {1,4}, N₂ = {2,6,7}, and N₃ = {3,5,8} see in Table 5. This means that job 1 and job 4 must be finished in order to assemble product 1 Table 2.

Table 3. Problem setting for verification and validation process: a small case

	t	n	f	m	ATINT	DDINT
Problem 1	3	8	2	2	INT[20,35]	INT[5,20]
Problem 2	3	8	2	2	INT[5,20]	INT[20,35]

Table 4 The production time for each job

Product (s)	Problem 1		Problem 2
	ATs	DDs	ATs	DDs
1	24	13	13	59
2	28	77	20	44
3	31	77	6	11

Table 5 Eligibility constraints and production time for each job

Job (j)	Eligibility Constraints (LFj)	Product Index (IPj)	Production Time (PTjqk)
			Problem 1				Problem 2
			PTj1k		PTj2k		PTj1k		PTj2k
			PTj11	PTj12	PTj21	PTj22	PTj11	PTj12	PTj21	PTj22
1	1	1	11	14	0	0	28	26	0	0
2	1,2	2	6	6	7	8	21	22	21	25
3	1	3	13	17	0	0	23	26	0	0
4	1	1	14	12	0	0	20	34	0	0
5	1	3	0	0	5	10	0	0	27	25
6	1,2	2	7	17	9	13	34	25	21	32
7	1	2	7	14	0	0	26	32	0	0
8	2	3	0	0	9	15	0	0	26	21

There are 73,728 feasible solutions based on the computational results of the MILP model and total enumeration technique. The best non-dominated solution can be selected from all feasible solutions, which is expressed as the Pareto Front (PF). Non-dominated solution of each problem is given in Table 6.

Table 6. Pareto front of each problem

	Mean Flow Time ($\bar{F}$)	Number of Tardy Product ($n_T$)
Problem 1	65,6667	1
	63,3333	2
Problem 2	55,3333	2
	55,6667	3

Two significant findings were drawn as a result of the verification process: 1) the procedure for finding solutions was in accordance with the logic of the flow chart and conceptual model that had been defined previously. There are four interrelated decisions that must be made, namely: determining the feasible and effective factory for each job; determining job priorities for selecting machines in a factory; assigning jobs to each machine; and determining the sequence of the job assembling process into the final product at the final stage. Figure 6 and Table 7 illustrate the problem, best feasible solution, and optimal solution. 2) The computer program contains no coding errors, and the solution structures are mathematically correct and logically consistent. Validation is done to ensure that the algorithm and programming code can represent and solve the proposed problem properly.

Figure 6. Illustration of the SH-1 solution structure for problem 1.

Table 7. Solution representation of MDPMASP with eligibility constraints

	Alg	π	πF1	πF2	πFM11	πFM12	πFM21	πFM22	Objective Value	Optimal Solution
Problem 1	SH-1	1 2 3	4 1 7 2 3	6 8 5	4 7	1 2 3	6 5	8	67.0000 2.0000	65.6667 1.0000
		2 3 1	7 2 3 4 1	6 8 5	7 4 1	2 3	6 5	8	65.6667 1.0000	63.3333 2.0000
	SH-2	1 2 3	4 1 7 2 3	6 8 5	1 7 3	4 2	6 5	8	65.0000 2.0000
		2 3 1	7 2 3 4 1	6 8 5	7 3 1	2 4	6 5	8	65.6667 1.0000
	SH-3	1 2 3	4 1 7 2 3	6 8 5	4 7	1 2 3	6 5	8	67.0000 2.0000
		2 3 1	7 2 3 4 1	6 8 5	7 4 1	2 3	6 5	8	65.6667 1.0000
	SH-4	1 2 3	4 1 7 2 3	6 8 5	1 7 3	4 2	6 5	8	65.0000 2.0000
		2 3 1	7 2 3 4 1	6 8 5	7 3 1	2 4	6 5	8	65.6667 1.0000
	MGA	2 1 3	2 7 4 1 3	6 8 5	2 7 1	4 3	6 5	8	63.3333 2.0000
		2 3 1	7 2 3 4 1	6 8 5	7 3 1	2 4	6 5	8	65.6667 1.0000
		2 1 3	7 2 4 1 3	6 8 5	7 1 3	2 4	6 5	8	63.3333 2.0000
	MSA	2 3 1	7 2 3 4 1	6 8 5	7 3 1	2 4	6 5	8	65.6667 1.0000
		2 1 3	2 7 4 1 3	6 8 5	2 7 1	4 3	6 5	8	63.3333 2.0000
		2 1 3	7 2 1 4 3	6 8 5	7 1 3	2 4	6 5	8	63.3333 2.0000
		2 3 1	2 7 3 4 1	6 8 5	2 7 4	3 1	6 5	8	65.6667 1.0000
		2 3 1	2 7 3 1 4	6 8 5	2 7 1	3 4	6 5	8	65.6667 1.0000
		2 1 3	2 7 1 4 3	6 8 5	2 7 3	1 4	6 5	8	63.3333 2.0000
		2 1 3	7 2 4 1 3	6 8 5	7 1 3	2 4	6 5	8	63.3333 2.0000
		2 3 1	7 2 3 1 4	6 8 5	7 3	2 1 4	6 5	8	65.6667 1.0000
	MPSO	2 3 1	7 2 3 4 1	6 8 5	7 3 1	2 4	6 5	8	65.6667 1.0000
		2 1 3	7 2 1 4 3	6 8 5	7 1 3	2 4	6 5	8	63.3333 2.0000
		2 3 1	2 7 3 4 1	6 8 5	2 7 4	3 1	6 5	8	65.6667 1.0000
		2 1 3	2 7 1 4 3	6 8 5	2 7 3	1 4	6 5	8	63.3333 2.0000
		2 3 1	7 3 1 4	6 2 8 5	7 4	3 1	6 8	2 5	65.6667 1.0000
		2 3 1	7 3 1 4	2 6 8 5	7 3 1 4		6 5	2 8	65.6667 1.0000
		2 1 3	2 7 4 1 3	6 8 5	2 7 1	4 3	6 5	8	63.3333 2.0000
		2 3 1	2 7 3 1 4	6 8 5	2 7 1	3 4	6 5	8	65.6667 1.0000
		2 1 3	7 2 4 1 3	6 8 5	7 1 3	2 4	6 5	8	63.3333 2.0000
		2 3 1	7 2 3 1 4	6 8 5	7 3	2 1 4	6 5	8	65.6667 1.0000
		2 1 3	7 2 1 4 3	6 5 8	7 1 3	2 4	6 8	5	63.3333 2.0000
Problem 2	SH-1	3 1 2	3 1 4 7 2	5 8 6	3 4 2	1 7	5	8 6	57.6667 2.0000	55.3333 2.0000
	SH-2	3 1 2	3 1 4 7 2	5 8 6	3 4 2	1 7	8 6	5	57.3333 2.0000
	SH-3	3 1 2	3 1 4 7 2	5 8 6	3 4 2	1 7	5	8 6	57.6667 2.0000
	SH-4	3 1 2	3 1 4 7 2	5 8 6	3 4 2	1 7	8 6	5	57.3333 2.0000
	MGA	3 1 2	3 4 1 7	8 5 6 2	3 4	1 7	8 6	5 2	55.3333 2.0000
	MSA	3 1 2	3 1 4 7	5 8 6 2	3 4	1 7	8 2	5 6	55.3333 2.0000
	MPSO	3 1 2	3 1 4 7	5 8 2 6	3 4	1 7	8 2	5 6	55.3333 2.0000
		1 3 2	1 4 3 7	8 5 6 2	4 3	1 7	5 6	8 2	55.3333 2.0000

Table 8 shows the convergence and divergence values of each proposed algorithm. The convergence value represents the proximity of the resulting solution to the optimal solution. The smaller the convergence value, the better the algorithm’s performance. While the value of divergence is used to assess the diversity of the resulting solutions. The smaller the divergence value (closer to zero), the more diverse solutions the algorithm can produce. The computational results show that all the proposed algorithms perform very well in solving multi-objective MDPMSP problems with eligibility constraints, although there are differences in the level of accuracy. So it can be concluded that the proposed model (MILP, algorithm, and program) has been verified, is valid and can represent and solve problems very well. As a result, we can perform further experiments to find out how the algorithm performs on larger problem sizes.

Table 8. Convergency and divergency matrix

Alg	Problem 1		Problem 2
	Convergency	Divergency	Convergency	Divergency
SH-1	0.8333	1	2.3334	1
SH-2	0.6009	1	2.0000	1
SH-3	0.8333	1	2.3334	1
SH-4	0.6009	1	2.0000	1
MGA	0	1	0	1
MSA	0	1	0	1
MPSO	0	1	0	1

5.2. Parameter tuning

The quality and effectiveness of the algorithm are affected by the parameter configuration. Therefore, to obtain optimal performance, these parameters must be adjusted by using an efficient design of experiments. The most commonly used and intuitive strategy of experimental design is the full factorial design William G and Gertrude M (1992) which loses its efficiency by increasing the number of parameters. Here, we applied the Taguchi approach Ross (1988), where a large number of decision variables would be tuned through a small number of experiments (Karimi & Davoudpour, 2014).

Taguchi’s experimental design method was used to analyze the impacts of the whole set of parameters using a small number of experiments based on an extraordinary orthogonal array design. Numerous control parameters are examined during factor analysis in order to optimize the performance of the proposed metaheuristic algorithm, i.e. population size for MGA (N_MGA), mutation probability (P_m), maximum number of generation (Gmax), initial temperature (T₀), maximum number of iteration for MSA (Itmax_MSA), temperature reduction constant (Fr), minimum temperature (ε), population size for MPSO (N_MPSO), maximum number of iteration for MPSO (Itmax_MPSO), and initial velocity (v₀). Each parameter is categorized into three levels, and DoF = 6. An orthogonal array $L_9\left(3^4\right)$ with one dummy factors in MGA and MPSO was chosen as the experimental design and summarized in Table 9. The Taguchi approach results in nine different treatments, which have to be run five times.

Table 9. Parameter configurations

Trial	Algorithm	Configuration of Parameter Value
1	Modified Genetic Algorithm	NMGA (50)	Pm (0.3)	Gmax (5)	-
2		NMGA (50)	Pm (0.5)	Gmax (10)	-
3		NMGA (50)	Pm (0.7)	Gmax (20)	-
4		NMGA (100)	Pm (0.3)	Gmax (10)	-
5		NMGA (100)	Pm (0.5)	Gmax (20)	-
6		NMGA (100)	Pm (0.7)	Gmax (5)	-
7		NMGA (200)	Pm (0.3)	Gmax (20)	-
8		NMGA (200)	Pm (0.5)	Gmax (5)	-
9		NMGA (200)	Pm (0.7)	Gmax (10)	-
1	Modified Simulated Annealing	T0 (10^2)	ItmaxMSA (10)	Fr (0.3)	ε (10^−3)
2		T0 (10^2)	ItmaxMSA (20)	Fr (0.5)	ε (10^−3)
3		T0 (10^2)	ItmaxMSA (50)	Fr (0.7)	ε (10^−3)
4		T0 (10^4)	ItmaxMSA (10)	Fr (0.5)	ε (10^−3)
5		T0 (10^4)	ItmaxMSA (20)	Fr (0.7)	ε (10^−3)
6		T0 (10^4)	ItmaxMSA (50)	Fr (0.3)	ε (10^−3)
7		T0 (10^6)	ItmaxMSA (10)	Fr (0.7)	ε (10^−3)
8		T0 (10^6)	ItmaxMSA (20)	Fr (0.3)	ε (10^−3)
9		T0 (10^6)	ItmaxMSA (50)	Fr (0.5)	ε (10^−3)
1	Modified Particle Swarm Optimization	NMPSO (50)	ItmaxMPSO (5)	v0 (0)	-
2		NMPSO (50)	ItmaxMPSO (10)	v0 (0)	-
3		NMPSO (50)	ItmaxMPSO (20)	v0 (0)	-
4		NMPSO (100)	ItmaxMPSO (5)	v0 (0)	-
5		NMPSO (100)	ItmaxMPSO (10)	v0 (0)	-
6		NMPSO (100)	ItmaxMPSO (20)	v0 (0)	-
7		NMPSO (200)	ItmaxMPSO (5)	v0 (0)	-
8		NMPSO (200)	ItmaxMPSO (10)	v0 (0)	-
9		NMPSO (200)	ItmaxMPSO (20)	v0 (0)	-

The signal-to-noise (S/N) ratio is used to analyze variations of proposed algorithms performance in response to each set of parameters. The signal describes the desired impact of the yield, while the noise represents the undesired impact. Because the S/N ratio comes from the quadratic loss function, three examples are widely used which are (1) the nominal is best, (2) the smaller is better and (3) the bigger is better. Based on the scheduling problem’s features, the S/N ratio is applied for the smaller-the-better of the convergence and divergence metric,

(31) $S/Nratio=-10log\frac{1}{R}\sum _{r=1}^R{(con{v}_r\ast di{v}_r)}^2$(31)

where $n$ is the number of replications and $y$ is the output of the $i$ th replication.

While conducting each experiment, the obtained objective function, i.e. convergence and divergence metrics, must be converted according to the Taguchi method and proportional to the signal-to-noise ratio, and analysis is provided according to its variations. In the Taguchi method, the S/Ns ratio indicates the variable ratio and the objective function in each experiment which is converted to this ratio in order to make a decision. Moreover, the convergence value will be measured based on the Pareto front (PF) value obtained by combining the entire solution set from each algorithm, then from that set the non-dominated solution set (Deb et al., 2002).

At this stage, the experiment was carried out with as many as 30 replications for each trial with t = 3, n = 8, f = 2, m = 2, PT_INT = INT[10 30] and AT_INT = INT[30 50]. All processing time and other information are generated randomly. All of the obtained solutions from the Taguchi design of experiments for adjustment of parameters have been shown as S/Ns ratio in Figure 7 and the best factor levels for the proposed metaheuristic algorithms are provided in Table 10.

Figure 7. Mean S/Ns ratio.

Table 10. The optimum level of each parameter

Algorithm.	Parameter	Optimum level
MGA	Population size for MGA, NMGA	200
	Mutation probability, Pm	0.3
	The maximum number of generations, Gmax	20
MSA	Initial temperature, T0	10^4
	The maximum number of iterations or cycles for MSA, ItmaxMSA	50
	Temperature reduction constant, Fr	0.7
	The minimum temperature, ε	10^−3
MPSO	Population size for MPSO, NMPSO	200
	The maximum number of iterations or cycles for MPSO, ItmaxMPSO	20
	Initial velocity, v0	0

5.3. Comparison algorithm and evaluation

Performance measurement is carried out to determine the quality of the proposed algorithm. Because the multi-objective problems result in multiple non-dominated solutions, the performance of several algorithms compared is a complex problem and requires some performance metrics. Through this research, all of the proposed algorithms are measured by four performance metrics, i.e. overall non-dominated vector generation (ONVG), overall number of non-dominated solution (ONSN), average quality (AQ), and distance metric (DI_R) proposed by (Qian et al., 2009). The experimental parameters are informed in Table 11. It was generated based on the L₈(2⁷) orthogonal array of the Taguchi experimental design, with five two-level factors and DoF = 5. Thirty runs were generated for each set of problems using Matlab software.

Table 11. The experimental parameters

Set of Problem Parameters	t	n	f	M	ATintv	PTintv
SPP 1	5	25	3	2	radint (30 50)	radint(10 30)
SPP 2	5	25	3	10	radint(10 30)	radint (30 50)
SPP 3	5	70	6	2	radint (30 50)	radint(10 30)
SPP 4	5	70	6	10	radint(10 30)	radint (30 50)
SPP 5	7	25	6	2	radint(10 30)	radint (30 50)
SPP 6	7	25	6	10	radint (30 50)	radint(10 30)
SPP 7	7	70	3	2	radint (30 50)	radint(10 30)
SPP 8	7	70	3	10	radint(10 30)	radint (30 50)

The comparison of proposed algorithms is presented in Figure 7. The results show that all proposed algorithms can obtain satisfactory solutions for given MDPMASP with eligibility constraints. But the statistical performance metrics show that there exist some differences between them.

The ONVG matrix is used to measure the number of non-dominated solution points generated by the proposed algorithm. The non-dominated solution set generated by the MSA algorithm produces a larger average ONVG metric than that generated by the other algorithms. This means that the MSA algorithm is able to produce more non-dominated solutions than the other proposed algorithms even though MSA has a smaller number of feasible solutions than MGA and MPSO. From the results of the Kruskal–Wallis statistical test (see, Table 12), each proposed algorithm has a significant effect on alpha = 0.05 on the ONVG matrix.

Table 12. Kruskal–Wallis test

Instance	Mean Rank							Kruskal–Wallis (Asym.Sig < 0.05)
	SH-1	SH-2	SH-3	SH-4	MGA	MSA	MPSO
SPP1	20.950	20.950	20.950	20.950	43.700	55.500	45.800	Y(0.000)
SPP2	20.500	20.500	20.500	20.500	45.500	55.500	49.350	Y(0.000)
SPP3	22.000	22.000	22.000	22.000	40.800	54.200	42.600	Y(0.000)
SPP4	22.900	22.900	22.900	22.900	35.900	55.500	41.200	Y(0.000)
SPP5	26.450	26.450	26.450	26.450	21.700	55.500	31.520	Y(0.000)
SPP6	24.200	24.200	24.200	24.200	30.700	55.500	39.950	Y(0.000)
SPP7	26.800	26.800	26.800	26.800	20.300	55.500	28.450	Y(0.000)
SPP8	22.000	22.000	22.000	22.000	39.500	55.500	43.540	Y(0.000)

*Y if there is a significant effect on α = 0.05 and N if not.

Therefore, the Post-Hoc test (Mann-Whitney U test) needs to be done to find out which algorithm gives the best performance. Based on the results of the Mann–Whitney test (Table 13), it can be concluded that among the presented algorithms, the MSA and MPSO algorithms have a statistically significant number of dominant solutions at all instances.

Table 13. Mann-Whitney U test

No	Proposed Algorithms		Asymp. Sig (2-tailed)
			SPP1	SPP2	SPP3	SPP4	SPP5	SPP6	SPP7	SPP8
1	SH-1	SH-2	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
2		SH-3	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
3		SH-4	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
4		MGA	Y(0.000)	Y(0.000)	Y(0.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
5		MSA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
6		MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
7	SH-2	SH-3	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
8		SH-4	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
9		MGA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
10		MSA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
11		MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
12	SH-3	SH-4	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
13		MGA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
14		MSA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
15		MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
16	SH-4	MGA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	N(1.000)	N(1.000)	N(1.000)	N(1.000)
17		MSA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
18		MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
19	MGA	MSA	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
20		MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)
21	MSA	MPSO	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)	Y(0.000)

*Y if there is a significant effect on α = 0.05 and N if not.

Meanwhile, based on the overall non-dominated solution number (ONSN) metrics, it can be seen that the most competitive non-dominated solution can be obtained from the modified genetic algorithm (MGA). From DI_R metric, it can be concluded that the modified GA can obtain a solution closer to the Pareto Front than the other. In addition, DI_R can point out proximity and diversity, which are two important features for multi-objective optimization, and they may influence each other's modified GA to some extent. Thus, the AQ metric is used for comprehensive comparison. From Figure 8, it can be seen that the AQ values of the modified MGA are less than those of the other. This means that the obtained solutions of the MGA are closer to the Pareto front than those of the other and the obtained solutions of the MGA distribute more uniformly and cover more area of the Pareto front. So, it is concluded that the MGA is a most effective and efficient multi-objective optimization algorithm compared to the other.

Figure 8. Comparison of proposed algorithms.

5.4. Practical insights

MDMASP problems with eligibility constraints are often encountered in manufacturing companies with the following characteristics: 1) intermittent manufacturing process; 2) production process with low volume and high product variety; 3) consists of several dedicated fabrication lines (factories) and a single assembly line. Each component can only be produced at a specific factory and assembled when all components of a product are finished processing. Each factory has several non-identical parallel machines that have different capacities and production speeds. From another perspective, the factory can also be interpreted as a supplier that produces semi-finished goods to be assembled later.

In addition, to minimize the mean flow time, we can reduce the inventory of semi-finished products, improve machine utility, and speed up order completion time. Moreover, by minimizing the number of tardy jobs, we can prevent delays in order completion, increase customer satisfaction, and avoid the risk of penalty fees.

6. Conclusion

We proposed a new model of Multi-Objective Modified Distributed Parallel Machine and Assembly Scheduling Problem (MDPMASP) with Eligibility Constraints and four simple heuristic algorithms and two iterative algorithms to solve it. The objective functions are minimizing mean flow time and the number of tardy products simultaneously. The final results of the proposed algorithms were dependent on its initial parameters. Therefore, the Taguchi design of experiment was used to find the best combination of parameters. Then to measure the performance of the proposed algorithms, computational experiments were performed on a set of experimental parameters.

The optimal combination of parameters is N_MGA(200), P_m(0.3), G_max(20), T₀(10⁴), Itmax_MSA(50), F_r(0.7), ε(10⁻³), N_MPSO(200), Itmax_MPSO(20), and v₀(0). Based on the evaluation of small cases, it can be concluded that the proposed algorithms were able to solve all problems with their complexity, although there were varying levels of performance. The comparison of each algorithm showed that the MGA is the best algorithm, followed by MSA, MPSO, SH-2, SH-4, SH-1 and SH-3. For the future research, there should be a focus on the development of algorithms based on other metaheuristic approaches, such as ant colony optimization (ACO) or other evolutionary methods.

This problem has a lot of potential to be applied and developed, it is possible to add some practical constraints. For example, sequence-dependent setup times, parallel assembly lines or flowshop assembly lines, objective functions related to due date, energy consumption or machine utility levels, as well as multiple objective functions. Moreover, various optimization models and heuristics also have the potential to be developed. This can significantly improve the quality and speed of the solution search process. We expect this can be a trigger for ideas for more realistic future work.

Disclosure statement

No potential conflict of interest was reported by the author(s).