Two meta-heuristic algorithms for optimizing a multi-objective supply chain scheduling problem in an identical parallel machines environment

Article history: Received January 25 2021 Received in Revised Format March 17 2021 Accepted March 21 2021 Available online March, 21 2021 Optimizing the trade-off between crucial decisions has been a prominent issue to help decisionmakers for synchronizing the production scheduling and distribution planning in supply chain management. In this article, a bi-objective integrated scheduling problem of production and distribution is addressed in a production environment with identical parallel machines. Besides, two objective functions are considered as measures for customer satisfaction and reduction of the manufacturer’s costs. The first objective is considered aiming at minimizing the total weighted tardiness and total operation time. The second objective is considered aiming at minimizing the total cost of the company’s reputational damage due to the number of tardy orders, total earliness penalty, and total batch delivery costs. First, a mathematical programming model is developed for the problem. Then, two well-known meta-heuristic algorithms are designed to spot near-optimal solutions since the problem is strongly NP-hard. A multi-objective particle swarm optimization (MOPSO) is designed using a mutation function, followed by a nondominated sorting genetic algorithm (NSGA-II) with a one-point crossover operator and a heuristic mutation operator. The experiments on MOPSO and NSGA-II are carried out on small, medium, and large scale problems. Moreover, the performance of the two algorithms is compared according to some comparing criteria. The computational results reveal that the NSGA-II performs highly better than the MOPSO algorithm in small scale problems. In the case of medium and large scale problems, the efficiency of the MOPSO algorithm was significantly improved. Nevertheless, the NSGA-II performs robustly in the most important criteria. © 2021 by the authors; licensee Growing Science, Canada


Introduction
It is essential for supply chains that customers and manufacturers work together in a coordinated approach in order to achieve information sharing, communication, and high efficiency. For the sake of these objectives, there should be a quick flow of information between manufacturers, retailers, distribution centers, delivery systems, and customers (Simchi-Levi, 2004). The challenge of the integrated supply chain scheduling at different levels has been one of the most significant issues in Supply Chain Management (SCM) and Operations Research (OR) (Thomas & Griffin, 1996). In classical production environments, schedules and delivery dates are fixed and certain, and the transportation unit is not taken into account. Moreover, decisions on production scheduling and delivery schedules are considered separately (Ho & Chang, 1995). Thus, the Integrated Production and Distribution Scheduling (IPDS) model refers to issues that consider production and distribution scheduling simultaneous and integrated (Chen & Vairaktarakis, 2005;Wang et al., 2014). In IPDS problems concerning the batch delivery system, dispatching jobs should wait until all of them are processed in the same batch. Accordingly, the shipping date of each batch has to be equal to the last completed job, which belongs to the batch, and the number of batches that are sent for the customer affects delivery costs. On the one hand, the waiting time of the jobs for processing the batches may raise the number of tardy jobs, resulting in an increase in tardiness costs. On the other hand, sending several jobs in batches brings about a decline in delivery costs and earliness penalties (Chen, 2010). Therefore, there is a trade-off between the tardiness penalties, the operation time of the entire process (makespan), and the earliness penalties and delivery costs. In this regard, simultaneous optimization of this trade-off can make the problem more complicated, but makes it closer to reality. One of the conventional production systems in the supply chain scheduling problems is the "parallel machines" environment. A parallel production planning problem consists of a set of jobs which should be processed and scheduled on a set of parallel machines. Hence, a sequence of processing jobs must be found in order to optimize the total operation time and other supply chain costs using one or more objective functions (Cheng & Sin, 1990;Sivrikaya-Şerifoǧlu & Ulusoy, 1999). In this article, the sequential processing of jobs is considered as the production ability, in which each job has only one action that can be processed by each machine. Then, the other jobs from the same batch are sequentially processed; thereby, the batch processing time equals the sum of processing times of all jobs belonging to it. Wang and Cheng (2000) performed the first research on parallel machines scheduling with batch delivery aiming at minimizing the total operation time and delivery costs. Lenstra et al. (1977) have shown that parallel production scheduling is a strongly NP-hard problem for ≥ 3. Furthermore, if the problem is considered with multiple objectives, it will become more complicated (Wang & Liu, 2015). Hence, researchers have begun to develop approximated solution methods, and as a result, several meta-heuristic algorithms have been introduced for solving multi-objective scheduling problems (Coello, 2002;Tyagi & Gupta, 2018). As claimed by Deb (2014), three solution methods, including a priori, a posteriori, and interactive approaches are suggested to solve multi-objective optimization problems (MOPs). According to the second approach, the Pareto curve is created by a set of Pareto solutions. So, decision-makers are able to choose some solutions from the Pareto curve as the optimal solution. In the last decade, researchers have taken into account the second approach in an attempt to resolve multiple objectives problems (Ojstersek et al., 2020). In this study, a bi-objective scheduling for production and distribution will be integrated in a production environment with identical parallel machines. Moreover, minimizing total weighted tardiness and total operation time is taken into account as the first objective function. In addition, minimizing the total weighted number of tardy orders, total earliness penalty, and total batch delivery costs is considered as the second objective function. Then, a bi-objective problem is modeled using a mixed integer linear programming (MILP) approach. Subsequently, two meta-heuristics algorithms are proposed to solve the problem. The first one is the multi-objective particle swarm optimization (MOPSO), and uses a mutation function. The second one is the non-dominated sorting genetic algorithm (NSGA-II), and uses a one-point crossover operator and a heuristic mutation operator. The Taguchi parameter tuning method is applied to adjust the algorithms' parameters. After setting the parameters, the performance of the two algorithms is examined in small, medium, and large scale problems, using five comparing criteria. Finally, the interpretations of the experiments are presented.
The rest of the article is organized as follows: In Section 2, the literature on the IPDS problem is reviewed to identify the research gap and explain the main contributions of this article. In Section 3, the problem definition is presented, and the notations used throughout the paper are given. Section 4 provides the mathematical model for this problem, including the objective functions and constraints, along with their explanations. In Section 5, the details of the MOPSO and the NSGA-II algorithms are explained. In Section 6, experimentation and computational results are provided. The last section (Sect.7) draws conclusions and proposes future research directions.

Literature review
A review of the studies on the IPDS problem in the parallel machines environment is provided. Furthermore, the researchers' achievements in the meta-heuristic algorithms in this field, as well as the main contributions of this article, are also presented in this section.

Integrated production and distribution scheduling
Most studies have focused on the IPDS problem with respect to maximizing customer service levels and minimizing customer delivery time, where delivery costs were not considered. Chen (2010) conducted a proper review of scheduling problems with regard to dispatching. He introduced the symbol / / / / to illustrate related problems. Here, represents the structure of machines, represents the constraints and specific conditions. Besides, denotes the characteristics of the delivery process shown as , , which and represent the number of vehicles and their capacity, respectively. Furthermore, indicates the number of customers, and finally, denotes the objective function(s). To cite an example, Potts (1980) was among the first to investigate an IPDS problem. This problem involves a type of delivery time regardless of transportation costs, and is in the category of single and immediate delivery. As one of the first studies on batch delivery and transportation costs, Cheng and Kahlbacher (1993) examined the 1// (∞, ∞), /1/ ∑ + problem aiming at minimizing the total weighted earliness and delivery costs, and revealed that it is an ordinarily NP-hard problem. Chang et al. (2014) investigated an IPDS problem in order to minimize total job delivery time and batch delivery costs. Besides, an ant colony optimization (ACO) was designed to solve the problem. Gao et al. (2015) studied the IPDS problem concerning batch delivery with limited vehicle capacity to multiple customers aiming at minimizing the makespan. They also proved that the general version of the IPDS problem is strongly NP-hard. Aminzadegan et al. (2019) investigated the 1/ ℎ/ (∞, ∞), problem with a single machine and two types of customers. To this end, they integrated production, distribution, and resource allocation to minimize tardiness, lost sale and batch delivery costs, and resource allocation. Moreover, an ant lion optimization (ALO) algorithm and adaptive genetic algorithm (AGA) were customized to solve the problem. In a make-to-order (MTO) strategy, Liu et al. (2020) considered the IPDS problem with vehicle routing problem as transportation stage aiming at minimizing the total order delivery time. To solve the problem, the variable neighborhood search (VNS) algorithm was applied using ten neighborhood structures and a local search algorithm. As a result, the effectiveness of the algorithm was superior to memetic and large neighborhood search (LNS) algorithms.
The problem of supply chain scheduling in parallel machines environment has been studied by researchers in recent years (Munoz-Villamizar et al., 2019;Nikabadi & Naderi, 2016). Lin and Jeng (2004) examined the / / (∞, ∞), /1/ and / / (∞, ∞), /1/ ∑ problems to minimize tardiness and the total number of late orders. Their computational results indicated that the dynamic programming algorithm was able to solve the small size problems, and problems up to 500 jobs and five machines were solved using heuristic algorithms. Considering the different setup times for unrelated parallel machines, Liu and Lu (2016) investigated the → , ℎ/ − , = 1, ≥ 1/ and → , ℎ/ − , = 1, ≥ 1/ problems. In their model, orders are processed on parallel machines and then they are delivered to customers. They also considered the availability and unavailability of the machines and interruption of jobs constraints for these problems. Joo and Kim (2017) considered the IPDS problem with sequence-dependent setup time and minimized the makespan to ascertain machine scheduling, batching, and distribution schedule. The problem was solved by several types of the GA framework and the analysis indicated the outperformance of the GA_CRC type algorithm. The work of Ekici et al. (2019) proposed an MTO strategy aiming at minimizing total tardiness and earliness penalties in the attendance of sequencedependent setup time and the other features of unrelated parallel machines. To address the problem, three heuristic methods and a Tabu search-based meta-heuristics were designed.
In studies carried out on the IPDS problem with identical parallel machines, Shim and Kim (2008) examined an identical parallel machines scheduling problem to minimize total tardiness of orders with the notation of // (∞, ∞), , / 1/ ∑ . Ullrich (2013) considered the IPDS problem in an identical parallel machines environment with machine ready times, and delivery time windows for the vehicle routing problem. Afterward, a GA was adopted to handle the problem. Schaller (2014) solved the scheduling problem with an objective function in an attempt to entirely eliminate total tardiness in identical parallel machines by taking into account setup times for the machines. By using Tabu-search and GA, and the optimal branch and bound method, he showed that GA performed much better in problem-solving. Chen et al. (2015) studied the / / problem, which was a two-stage scheduling problem. That is, the orders are processed on identical parallel machines, and afterward are delivered to customers by a single vehicle. They solved the problem by using an approximate algorithm with a constrained boundary to minimize the maximum delivery time. Moreover, Chen et al. (2016) developed the work of Chen et al. (2015) by considering a set of due dates assigned to jobs. Wang et al. (2019) solved an IPDS problem with some characteristics such as machine-dependent ready times, and multi-trip vehicle routing with time windows and uncertain travel times. They proposed a robustness approach to handle uncertain travel times, and a memetic algorithm for minimizing the tardiness and delivery costs. Wang et al. (2019) considered an IPDS problem in the hybrid flow-shop environment including identical parallel machines in three stages. To address the problem, a four-layered constructive ( ), a hybrid heuristic method ( ), and a VNS algorithm were proposed to minimize the delivery completion time.

Use of the meta-heuristic algorithms
Synchronization of production and distribution scheduling is an NP-hard problem. Accordingly, several meta-heuristics algorithms (i.e., PSO, GA, ACO) have been designed and customized to solve these supply chain scheduling problems. Subsequently, multi-objective meta-heuristics algorithms (i.e., NSGA-II, MOPSO, SPEA-II, MOACO) have also been proposed by researchers to optimize the trade-off between various objectives in supply chain scheduling problems. Cakici et al. (2012) proposed a bi-objective GA with different heuristics to minimize total weighted tardiness and delivery costs for an MTO scheduling problem. Besides, the work of Cakici et al. (2014) solved an IPDS problem in a parallel machines production environment with a single vehicle and multiple customers. Hamidinia et al. (2012) examined the 1// ∑ ∑ + + ℎ + problem with a single machine for processing jobs and multiple customers. Also, mathematical programming was designed to minimize total earliness, tardiness, holding penalties, and batch delivery cost. Furthermore, they used a GA and integer programming to solve the problem. The findings demonstrated that the proposed GA had high performance. Rajkanth et al. (2017) examined the single machine and parallel machines scheduling problems, in which the parallel machines problem was solved using two models including the jobs sequence and jobs assignment. As a result, an appropriate upper limit was achieved from a GA to assess the efficiency of the meta-heuristic algorithm provided for parallel machine scheduling. Raghavan et al. (2018) examined a scheduling problem for challenges faced by factories that produced certain electronic components, including reworking and re-processing. They showed that, compared to the CPLEX solver and modified shortest total estimated processing time (MSTEPT) algorithm, the modified GA had the best quality solutions in small to medium size problems. However, the MSTEPT outperformed the modified GA in large-scale problems. Assarzadegan and Rasti-Barzoki (2016) investigated the 1/ (∞, ∞), scheduling problem based on the due date assignment for a single machine and several customers. Then, they solved the problem by using the parallel simulated annealing (PSA) and AGA algorithms. Jiang et al. (2017) examined the / ℎ, , , ≤ / problem, which is a scheduling coordination problem on identical parallel machines along with batch delivery. They applied some constraints such as different processing times and sizes, as well as a smaller number of jobs than the capacity of the machine. To solve the problem, they used a combination of the DPSO and GA algorithms. Saeidi (2016) solved a multi-objective scheduling problem on parallel machines with the intention to minimize both the makespan and machine costs by using MOPSO and NSGA-II. By comparing the obtained solutions, the NSGA-II demonstrated a higher performance than the MOPSO algorithm. Hassanzadeh et al. (2016) addressed a bi-objective IPDS problem in the flow-shop environment using the HCMOSPO and HBNSGA-II algorithms. Then they compared the results with the MOPSO algorithm and the NSGA-II. They reported that the developed algorithms had a superior performance than the other algorithms. Sheikh et al. (2018) developed multiple-objectives scheduling problem in a flow-shop environment. Additionally, three objective functions were considered, including minimizing total tardiness, total completion time, and the makespan. Also, the problem was solved using heuristic and meta-heuristics algorithms, i.e., the NSGA-III and the MOPSO algorithm. The results indicated the superiority of the meta-heuristics algorithms, with the NSGA-III exceeding the MOPSO algorithm in one of the three performance criteria. Ganji et al. (2020) applied three meta-heuristic algorithms, including MOACO, NSGA-II, and MOPSO, to address a green multi-objective IPDS problem to minimize customer dissatisfaction, total tardiness penalty, fuel costs, emitted carbon from vehicles, and distribution costs. The outcomes of experiments revealed the outperformance of NSGA-II compared to the other two algorithms. Some studies were performed in the context of multiobjective optimization parallel machines scheduling problems with meta-heuristic algorithms. Guo et al. (2016) studied a biobjective IPDS problem in an unrelated parallel machines environment. The study aims to minimize the total number of delayed orders and the sum of holding costs, labor costs, tardiness and earliness costs, and batch delivery costs. A bi-level evolutionary optimization (BLEO) algorithm was used to solve the problem. Shahidi-Zadeh et al. (2017) developed a multiobjective problem in the unrelated parallel machines environment, and considered characteristics such as ready jobs, release time, and the batch's capacity constraint. Two objective functions were considered, including minimization of total tardiness, total earliness, and costs of purchasing the machinery, and minimization of the makespan. The results showed a better performance of the bi-objective harmonic search (MOHS) algorithm compared to the other algorithms. Zhou et al. (2018) considered the / , / , problem. They defined two objective functions, including minimization of total electricity consumption costs, which is an environmental sustainability indicator, and minimization of the makespan. They solve the large instances by using a multiple-objectives discrete differential evolution algorithm. Afterward, the performance was compared to the NSGA-II and AMGA (the archive-based micro-genetic algorithm). It was revealed that the suggested algorithm had a preferable performance regarding the quality of solutions. Wu and Che (2019) considered a multi-objective energy-efficient scheduling problem for unrelated parallel machines. Afterward, they employed two objective functions in an attempt to minimize the total energy consumption and the makespan. They developed a memetic differential evaluation (MDE) algorithm and showed that the MDE algorithm was considerably superior than the multi-objective SPEA-II and NSGA-II.

Main contributions
There is a significant research gap according to the literature review. As can be seen from Table 1, optimizing the trade-off between crucial criteria of decision-making in the multi-objective IPDS problem has not been considered, to the best of our knowledge. So, a new bi-objective IPDS problem in the identical parallel machine environment will be optimized to fill the research gap. The first objective is considered as a measure to increase customer satisfaction (minimizing tardiness penalty and makespan). The second one is considered as a measure to reduce some manufacturer's costs such as warehouse cost, increasing the company's reputation, and transportation costs. Accordingly, the main contributions are as follows:  Concept: The idea of integrated production and distribution scheduling is conceptualized so as to help decision-makers for evaluating the trade-off between customer service costs and manufacturer's costs. After that, the idea is operationalized in a problem of supply chain scheduling in an identical parallel machines environment.
 Model: A new multi-objective MILP model is developed to assess the trade-off between two crucial objective functions. The first objective minimizes the total weighted lateness and total operation time. The second one minimizes the total earliness penalty, total weighted number of tardy orders, and total batch delivery costs.
 Solution approach: Two meta-heuristics algorithms are suggested to address the small, medium, and large scale problems. First, the MOPSO algorithm with a mutation function is proposed, followed by the NSGA-II with a onepoint crossover operator and a heuristic mutation operator. Afterward, the algorithms are compared using five comparing criteria. It should be noted that a trade-off between batch delivery costs and the customer's contractual due date is examined in detail as using these algorithms. Table 1 The literature review of associated researches and this study

Problem description
In this section, we explain the notations applied in the whole article, followed by describing the details and assumptions of the problem.

Notations
Indices

Details and assumptions of the problem
This study presents a bi-objective scheduling problem of production and distribution in the production environment with identical parallel machines, as defined by Chen (2010) It means that there is a processing system for jobs that includes parallel machines and the sequence-independent setup time to form each sub-batch of batch ; here, orders are delivered directly in batches to the customer. Finally, two objective functions are addressed. The following assumption are considered for the problem:  The number of independent jobs must be processed. Each job is processed at a time on one of the machines.  At time zero, all jobs are available.  The processing time for each job ( ) is definite and above zero.  Each machine can only process one job at a time and all machines are accessible at all time.
 The preemption of the job is not allowed.  Setup times are independent of the sequence for forming the sub-batches of each batch on the parallel machines.
Consider a supply chain scheduling problem in an identical parallel machines environment, which includes a production system with machines and jobs. A customer orders jobs to a manufacturer which each job has a priority weight of , a processing time of , and a contractual due date of . If the job's completion time ( ) is longer than its contractual due date ( ), the job is tardy, and a job is early if the job's completion time ( ) is lower than its contracted due date ( ). Subsequently, an earliness penalty of for early jobs and a tardiness cost of for tardy jobs are considered. Besides, the number of tardy jobs delivered to the customer has a notable impact on the manufacturer's reputation. Hence, there is another type of tardiness cost called "the cost of the company's reputational damage" that only depends on the number of tardy jobs based on their priority weight ( ) and not depends on the amount of tardiness.
In the manufacturing unit, jobs are sequentially processed in batches, and each batch is delivered to the customer by a vehicle with the capacity of . In this problem, when a job is completed in a batch, it waits for all jobs to be completed. Each batch can include one or several sub-batches that are processed on different machines. Afterward, the sequence-independent setup time ( ) is imposed when a sub-batch on any machine is created from each batch and it is applied before the beginning of the sub-batch processing. After this setup time, jobs in the same sub-batch are sequentially processed. So, the variable , is a time when all the sub-batch (SB) jobs of the batch on the machine are processed and equals to the last completed job which belongs to the sub-batch of batch . Moreover, variable is a time that all jobs in batch are processed and it equals to the latest completion time of the sub-batch ( , ). Here, is a binary variable and if the job is in the batch and on machine , its value is 1, and 0 otherwise. Similarly, in the case of , if the job is tardy, the value is 1 and otherwise 0. Moreover, in the case of , if batch is formed, i.e., it contains the job, its value is 1 and 0 otherwise. The batch delivery date to the customer is the maximum value of the jobs' completion times ( ) in that batch. So, the maximum completion time of all orders is equal to the total operation time (makespan). Later on, delivery costs can be reduced by batch delivery, in which the delivery system has one or several batches. In this problem, there are sufficient vehicles with the capacity of , and the delivery cost for each batch is equal to . For instance, if the batch delivery cost ( ) is raised by a third-party logistics company, the manufacturer will decide to decrease the dispatching number of batches; thereby, the tardiness penalty can increase, and vice versa. As a result, finding an optimal sequence of jobs processing on the parallel machines and assigning them to the batch are the primary purposes of this IPDS problem in order to optimize the trade-off between the decision variables. Fig. 1 indicates the schematic view of the problem.
Several jobs with different priority weights are ordered by the customer to the manufacturer The manufacturing unit with identical parallel machines with sequence-independent setup times to form each batch Sufficient capacityconstrained vehicles with different batch delivery costs for the customer One customer with different contractual due dates

Mathematical modeling
In this study, we consider the production environment with identical parallel machines and present a multi-objective model with mixed integer linear programming. Moreover, we attempt to minimize the most important costs based on the completion time of the orders and batch delivery. Hence, we attempt to minimize the two objective functions simultaneously. The objectives are (1) minimizing the total tardiness penalty and makespan, and (2) minimizing the total cost of the company's Production stage Delivery stage reputational damage due to the number of tardy orders, total earliness penalty and total batch delivery costs. The multiobjective problem is as follows: , , , , , , , , , , , The first objective function, presented in Relation (1), minimizes the total weighted tardiness and total operating time. Tardy jobs are crucial for producer in terms of the level of customer service and customer satisfaction with the timely delivery of orders. Therefore, one of our goals is to minimize total weighted tardiness ( ) and total operation time ( ). The second objective function, presented in Relation (2), minimizes the total cost of the company's reputational damage ( ); if the job is tardy, it imposes a reputation cost on the company based on its weight or priority, total earliness penalty ( ) and total batch delivery cost to the customer. Relations (3) and (4) calculate the amount of tardiness and earliness for job , respectively. Relation (5) determines whether or not job is tardy. Equation (6) assures that each job is only processed in one batch and on one machine because no preemption is allowed. Relation (7) specifies whether a sub-batch of batch is located on machine . Constraint (8) calculates the completion time for a sub-batch of batch processed on machine . Equation (9) considers the starting schedule at time zero. The completion time of batch is calculated in Constraint (10), and the delivery time of job is calculated in equation (11). Relation (12) calculates the makespan and constraints (13) and (14) guarantee that if the number of members in the batch is zero, it does not form and otherwise, the batch is formed and the cost of delivery is included in the objective function (2). Relations (15) and (16) are used to demonstrate the natures of the decision variables.
As observed above, the model is nonlinear since some variables are multiplied in relations (10) and (11). Thus, auxiliary variables are defined, including , = , .
in equation (11), and some relations are used to linearize constraints (10) and (11). The relations are defined as follows: The nonlinear model converts to a MILP model by replacing the above relations instead of constraints (10) and (11).

Solution approach
Solving real-life problems is impossible or takes a very long time, due to their large scales. Moreover, the mentioned problem is a bi-objective optimization problem and is strongly NP-hard (Lenstra et al., 1977). Hence, two multi-objective metaheuristics, including MOPSO and NSGA-II, are designed for solving the problem as approximated solution methods since the exact solution methods cannot efficiently solve this problem.
This section includes defining and explaining the details of the MOPSO and the NSGA-II algorithms presented in this article as efficient meta-heuristics algorithms for optimizing the mentioned bi-objective mathematical model.

The MOPSO approach
Eberhart and Kennedy (1995) first introduced the particle swarm optimization (PSO) algorithm. This algorithm is based on an uncertain search method to optimize the function inspired by the collective movement of birds or fishes seeking food. In the PSO algorithm, each particle has a fitness value calculated by a fitness function. Whatever a particle is closer to the objective (optimization) in the solution space, it has more competency. Coello et al. (2004) introduced a multi-objective type of the PSO algorithm called the MOPSO in 2004 and developed it in 2006 and 2007 (Coello et al., 2007;Coello, 2006). Since a unique optimal solution cannot be defined in the multi-objective optimization structure, the MOPSO uses a non-dominated solution procedure where particles randomly choose their leaders from an approximated Pareto curve. In this algorithm, a grid structure is defined to search for the solution space. Moreover, the best non-dominated solution is stored in an external memory (repository). After the particle updates its velocity and location, the best solution should be optimized. The pseudo-code of the MOPSO in detail is demonstrated in Appendix A.

The MOPSO algorithm characteristics
In this sub-section, the particles' structure is defined, and then the fitness function, and the mutation operator applied on the MOPSO are explained in detail.

Particles' structure
The MOPSO algorithm begins its work with a set of particles, which are the initial solution. Due to the need for a feasible solution, at first, the feasible solution must be saved according to a specific structure called particles' structure. Since the PSO algorithm has a continuous solution space, an initial solution is generated with random numbers between 0 and 1. Fig. 2 shows three strings generated by random numbers between 0 and 1 for an initial solution. The first string indicates the job number; the second string indicates the machine number allocated to the job. The batch number allocated to the job is also displayed in the third string; the number of jobs from 1 to is equal to the length of the strings.

Fitness function
Particles are created for the initial solution in a continuous space. Thus, the numbers between 0 and 1 must be transformed to integer numbers to turn the continuous solution space into a discrete solution space. Finally, the values of the objective functions can be computed.The first decision is to assign the machine number to the job. To do so, the random numbers in interval [0,1] are generated and equation (17) is used to transform random numbers between 0 and 1 to integer numbers for the string of machine numbers: According to Eq. (17), is the number of machines and is the random number between 0 and 1, which assigned to the job in the second string. Thus, illustrates the integer machine number, which is assigned to the job in the second string. The second decision is to assign the batch number to the job. To do so, the random numbers in interval [0,1] are generated and Eq. (18) is used to transform numbers between 0 and 1 to integer numbers for the string of batch number: According to Eq. (18), is the number of batches and is the random number between 0 and 1, which assigned to the job in the third string. So, indicates the integer batch number, which is assigned to the job in the third string. According to Fig. 2, suppose that the customer orders six jobs that must be processed on three machines and delivered to the customer in three batches. Figure 3 shows the method is used to transform the random numbers to integer numbers in the strings of Fig. 2 and also, both the machine number and batch number are assigned to each job.

Use of a mutation function
A mutation operator is utilized to intercept premature convergence, and better Pareto solutions can be obtained. In the MOPSO algorithm, after a mutation operation is performed on any of the genes in the structure, a new particle is created, in which the mutation rate ( ) is applied to each gene and the particle mutates with the probability of mutation ( ). In the mutation, a gene may be removed from a set of particle genes or a gene is added to the collection that has not so far existed in the population. Figure 4 indicates the mutation of genes in the strings of Fig. 2. Moreover, Fig. 5 shows a new mutated particle after the conversion from continuous to discrete operation.
The mutation function is performed in accordance with the following steps: 1. The genes are selected randomly in the strings based on a function of the number of jobs.

A random number with a uniform distribution is
3. The batch and machine numbers associated with the job are determined, and then allocated to their positions in the strings.  Goldberg and Holland (1988) were the first scholars to generalize the concept of multi-objective fitness problems. They used the non-dominated solution sorting concept and approached the algorithm for the Pareto solution set in consecutive iterations. Goldberg's idea in the concept of the Pareto set was directly used by Deb et al. (2002). Actually, by using the idea of preference for non-dominated solutions and assigning them more fitness, a GA was created with non-dominated sorting. The NSGA-II is among the most well-known algorithms for multiple objectives optimization. Following the introduction of the first version of the algorithm to provide a variety of optimal solutions for the Pareto curve in 2002, the developers of the algorithm generated an elitist mechanism based on the importance of non-dominated queues under the name NSGA-II. The pseudo-code of the NSGA-II in detail is shown in Appendix A.

The NSGA-II characteristics
In this sub-section, the structure of the chromosome is defined, and then a one-point crossover operator and a mutation operator are explained in detail.

Chromosomes' structure
Each repetitive meta-heuristic algorithm requires a structure (encoding) for representing solutions. A GA begins with the initial solutions of a set of chromosomes. Each of these chromosomes has a value for the objective function, called the fitness value. Subsequently, the chance of survival and reproduction is higher for a chromosome if it has a better fitness function value. These chromosomes are referred to as initial populations and can be produced in various ways, such as random reproduction, heuristic methods, and so on. In this paper, each chromosome has a length equal to the number of jobs from 1 to and has three rows: the first row indicates the job number, the second and third rows illustrate the machine number and batch number, respectively, which are randomly determined and assigned to the job. Therefore, the initial population is randomly generated. Figure 6 indicates the chromosome structure.

Fig. 6 A sample chromosome
Suppose there is a problem with six jobs, three batches, and two machines. Figure 6 represents the job is processed on the machine and delivered to the customer in the batch. For example, the second job is processed on the second machine and delivered to the customer in the first batch.

The crossover operator
The most important operator in genetic algorithms is the crossover operator. Crossover is a process in which the old chromosome generation is mixed and combined to form a new chromosome generation. The couples considered as the parent in the selection section will share their genes and create new members. The crossover in a GA eliminates the genetic dispersion or genetic diversity of the population, as it provides a condition so that good genes can be found by other genes. In this study, we use a one-point crossover operator used in the literature on the scheduling problems (please see (Bose et al., 2019;Shen,  2019)). The one-point crossover is applied to both second and third rows for machine assignment and batching the jobs. The one-point crossover operator on chromosomes is shown in Fig. 7.   Fig. 7 The crossover operator

The mutation operator
A mutation is another operator that creates other feasible solutions. In genetic algorithms, after generating a member in a new population, each gene mutates with the probability of mutation and produces a mutated population. In this paper, the mutation operator is applied according to the following steps: 1. The genes are selected randomly on the chromosome and the number of selected genes is determined by a function of the number of jobs. 2. To mutate the genes in the second row, random numbers are selected from the interval [0, ], where is the number of machines and assigned to the machine number for the job. 3. To mutate the genes in the third row, random numbers are selected from the interval [0, ], where is the number of batches and assigned to the batch number for the job. 4. Previous steps are repeated to generate new offsprings until a new mutated population is completed. Fig. 8 indicates the mechanism of mutation operator for machine assignment to the job and batching the jobs.

Fig. 8
The mutation operator

Experimentation and computational results
In this section, computational results are used to compare the efficiency and effectiveness of the NSGA-II and the MOPSO algorithm to solve small, medium, and large scale problems for nine classes presented in the following section (please see Sect. 6.2.). Firstly, the algorithm's parameters are tuned using the Taguchi design method. Secondly, the effectiveness of both algorithms is compared using the comparison criteria, and the results of the meta-heuristic algorithms are finally analyzed.

Parameter tuning
The efficiency of meta-heuristic algorithms is directly related to their parameters so that the correct choice of parameters increases the efficiency of the algorithm and also increases the quality of solutions. There are several methods for setting the parameters of the algorithms. In this research, the Taguchi method is used to ascertain the optimum value of parameters with their levels in algorithms. The controlled parameters of the Taguchi method are shown in Tables 2 and 3 for tuning parameters with their levels for the NSGA-II and the MOPSO algorithm, respectively. The problem is divided into three scales and the experiments are implemented for nine classes. Moreover, Minitab 16 was used to carry out all experiments related to the Taguchi method. For this purpose, a problem with 15 jobs and five machines is considered as a small scale problem, and a problem with 100 jobs and 15 machines is considered as a medium and large scale problem to adjust the parameters. Each Job number 1 2 3 4 5 6 7 Machine number 1 2 2 3 3 1 2 Job Batch 2 3 2 1 4 3 1

Parent 2
Job number 1 2 3 4 5 6 7 Machine number 1 2 1 1 2 3 3 Job Batch 2 3 3 2 3 1 2 Offspring 1 Job number 1 2 3 4 5 6 7 Machine number 2 3 2 3 3 1 2 Job Batch 1 4 2 1 4 3 1 Offspring 2 Job number 1 2 3 4 5 6 7 Machine number 2 3 1 1 2 3 3 Job Batch 1 4 3 2 3 1 2 Job number 1 2 3 4 5 6 7 Machine number 2 2 1 3 1 3 2 Job Batch 1 2 3 3 4 1 4 experiment number included nine classes and five random instances were generated for each class; each sample was solved five times with algorithms, and the average is presented. Therefore, 405 small scale and 405 medium and large scale sample problems are created to tune the parameters of the NSGA-II. In addition, the number of 1215 small scale and 1215 medium and large scale sample problems are created to adjust the parameters of the MOPSO. Finally, a total of 3240 sample problems were created. To use the Taguchi method, several comparison criteria are used so that the efficiency and effectiveness of each algorithm can be assessed (please see Sect. 6.3.). Moreover, a weight is allocated to each criterion based on its importance. Finally, each criterion is unscaled by the related deviation index (RDI), and then a response variable for the Taguchi method is achieved. Fig. 9 shows the signals to noise charts from the implementation of these tests. Based on the charts in Fig. 9, the selected NSGA-II and MOPSO algorithm factors levels are demonstrated in Tables 4 and 5, respectively.

Data Means for MOPSO Medium-and large-scale
Signal-to-noise: Smaller is better

Data generation
For the aim of evaluating the effectiveness of the NSGA-II and the MOPSO algorithm, problems are designed into two scales, then the performance of the algorithms is compared and analyzed using the five criteria. The meta-heuristics algorithms MOPSO and NSGA-II were encoded in MATLAB 2015b, and a Core i7 CPU and 8 GB RAM computer with Windows 10 (64 bit) was used as a test system. In small scale sample problems, the number of jobs is assumed to be 8, 12, 16, 20 and 24, and the number of machines is assumed to be 2, 4 and 6. In the medium and large scale sample problems, the number of jobs is considered to be 40, 80, 120, 180 and 250, and the number of machines is considered to be 8, 16, and 27. Other parameters were produced randomly based on a uniform distribution in the following intervals.

/ ]
A summary of the information above can be found in Table 6. Additional problem information, such as the number of machines and jobs, is also shown in Table 7.  The batch delivery costs can be associated with the priority of jobs ( ) and therefore, the delivery costs are introduced according to the weight of jobs. In addition, since the complexity of the problem could be affected by the delivery cost, three classes of the problem with low, medium, and high costs of delivery are generated: Classes A, B, and C represent low, medium, and high costs of delivery ( ), respectively. Furthermore, due to the impact of contractual due dates of jobs on the complexity of the problem, three sub-classes of problems of classes A, B and C are considered. Each sub-class 1, 2, and 3 represents the loose, moderate, and tight contractual due date ( ), respectively. Table 8 shows all nine classes. For any combination of jobs and machines (five jobs×three machines) in each class; from A1 to C3, there are 15 problems that five random instances generated for each problem. In the case of small scale, the number of sample problems is 675 (9 × 15 × 5), and 675 sample problems are also generated for medium and large scale. Each sample is run five times by the algorithms. The ultimate result equals the average of five runs for each class.

Comparing criteria
In single objective problems, an optimal solution is chosen for the purpose of the problem, while in multiple objective problems, we encounter a collection of Pareto solutions at the end of solving problems and the performance of the algorithms cannot be easily evaluated. Evaluation criteria are important since they can be used to examine the performance of algorithms.
In this research, we apply five comparison criteria from the literature to assess the performance of algorithms (Attar et al., 2014;Piroozfard et al., 2018;Zarei & Rasti-Barzoki, 2018). The comparing criteria are including quality metric criterion (QM), spacing metric criterion (SM), diversification metric (DM), CPU time (CPUT), and mean ideal distance (MID). A lower value of SM, MID, and CPUT criteria, and a higher percentage of QM criterion, and a higher value of DM criterion indicate better efficiency of the algorithm. The characteristics of comparing criteria are described in detail in Appendix B.

Comparison results and analysis for small scale samples
In this sub-section, the performance of NSGA-II and MOPSO algorithms is investigated for small scale problems according to the preceding criteria. Here, the charts of the three classes A, B, and C are examined, each being derived from the mean of classes A1 to C3 (The resulting data of each class is provided in Appendix C).

Comparison of algorithms by Standard Error and average value
The Standard Error (SE) and the average value (Mean) for each of the criteria are calculated for better analysis and assessment of the quality of the presented solutions. The quality of statistics is usually measured by using the. Thus, the performances of the two algorithms are compared in each of the nine classes A1 to C3 using the SE and the average value. The lower value of the SE is preferable. According to Fig. 15 and Table 9, the mean value and SE for the SM criterion are lower in all the classes in the NSGA-II than in the MOPSO algorithm. The NSGA-II has a higher average value in the QM criterion than in the MOPSO algorithm, and the SE is the same for both algorithms. In the MID criteria, it can be deduced by comparing the two algorithms that the MOPSO algorithm has a lower average value except for classes A2, B3, and C3. Moreover, the MOPSO algorithm has a lower SE in all the classes compared to the NSGA-II, meaning that the MOPSO algorithm is preferable in overall performance in this criterion. In the DM criterion, although the NSGA-II is more diverse except for classes A2, B3, C1, and C3, the SE in this criterion is lower for the MOPSO algorithm than for the NSGA-II in the rest of the classes. This means that the MOPSO algorithm has a more robust solution with a lower mean value for the DM criterion. Eventually, the MOPSO algorithm is shown to perform well on both factors in the CPUT criterion except for class C3.     SM  C3  B2  QM  B3  A3  MID  B3  A3  DM  C3  A3  CPUT C3 A1

Comparison results and analysis for medium and large scales samples
With regard to the mentioned criteria, this sub-section examines the performance of the NSGA-II and the MOPSO algorithm for problems with the medium and large scales. The charts of the three classes A, B, and C, each respectively, derived from the mean of classes A1 to C3, are investigated (The resulting data of each class is provided in Appendix C).     According to Figs. 16, 17, and 18, it could be noticed that the performance of the NSGA-II in the SM criterion is better than that of the MOPSO algorithm except in the tests of 8|180 and 8|250 for all classes. In addition, considering the MOPSO diagram, it can be inferred that the performance of this algorithm is enhanced with an increasing problem scale to medium and large. With increasing the number of jobs and machines in medium and large scale in the QM criterion, the difference between the two algorithms is less than that in the small scale for covering the non-dominated solutions, meaning the performance of the MOPSO algorithm improved greatly in high dimensions. In other words, with moving from medium to large in the number of jobs and machines, the efficiency of the two algorithms became closer to each other. However, at 60% problems of class A, 77% problems of class B, and 73% problems of class C, the NSGA-II had better performance than the MOPSO algorithm. In the MID criterion, the NSGA-II included 86%, 80%, and 86% of the solutions in classes A, B, and C, respectively, and had better performance than the MOPSO algorithm. In experiments such as 8|180 and 8|250, the MOPSO algorithm was able to perform better in all the classes than the NSGA-II, demonstrating that the performance of the NSGA-II was reduced with an increase in the number of jobs. However, when the number of machines increases, the NSGA-II again showed better performance. Compared to small scale problems, it can also be concluded that the performance of the NSGA-II improved in the MID criterion with increasing problem dimensions. Regarding Figs. 19 and 20, in the DM criterion, the MOPSO algorithm had a greater diversity of solutions than the NSGA-II in 66% problems of class A, 73% problems of class B and 66% problems of class C in all the tests. On the other hand, the NSGA-II had more diversification than the MOPSO algorithm in all the three classes in the tests such as 8|120, 8|180, 8|250, and 16|250, meaning that the NSGA-II had better efficiency in the fewer number of machines and more jobs, but had decreased efficiency with increasing the number of machines. At last, it can be realized that the MOPSO algorithm solves problems six times faster than the NSGA-II in all the classes.

. Comparison of algorithms by Standard Error and average value
The performance of the two algorithms in each of the nine classes A1 to C3 is compared using the SE and the average value for medium and large scale sample problems.
According to Fig. 21 and Table 11, the NSGA-II has a lower average value and lower SE in the SM criterion than the MOPSO algorithm and is more efficient. The NSGA-II allocates more non-dominated solutions than the MOPSO algorithm in the QM criterion, but both algorithms have the same SE. In the MID criterion, the comparison of the two algorithms shows that the average value in the NSGA-II is the lowest except for classes A3, B3, and C3. However, the MOPSO algorithm has more robust solutions with higher mean value since the SE value for all MOPSO classes is lower than the NSGA-II. In the DM criterion, it could be noticed that the MOPSO algorithm provides solutions with more robust features than the NSGA-II because of having a lower SE with a lower average. By comparing the two algorithms in the CPUT criterion, it can be  perceived that the NSGA-II has robust Pareto solutions with a higher average due to having a lower SE in all the classes except for classes A2 and A3, as compared to the MOPSO algorithm. Table 12 shows the best class in each of the criteria based on the SE and average value of the NSGA-II and the MOPSO algorithm for medium and large scale problems.

Fig. 21
Mean value for the two algorithms in medium and large scale sample problems for each class

Conclusion and suggestions
In this research, a bi-objective scheduling problem of production and distribution was integrated in an identical parallel machines environment. First, a mathematical model was presented for the problem. Then, two meta-heuristics algorithms, including the MOPSO algorithm and the NSGA-II were customized to solve the problem. The NSGA-II was used with a onepoint crossover operator and a mutation operator, in which the number of mutations is a function of the number of jobs to generate more mutated offsprings to spot the best optimum solution in the solution space. Also, the MOPSO algorithm was used with a mutation function for generating new particles to find better quality solutions. In the second step, the parameters tuning of the two algorithms was performed for small, medium, and large scale problems using the Taguchi method so that the algorithms could improve problem-solving. After tuning the parameters, the sample problems were solved in small, medium, and large scale problems. Besides, five comparing criteria including SM, QM, MID, DM, and CPUT were utilized to compare the performance and efficiency of the meta-heuristic algorithms in nine classes of problems by multiplying three levels of delivery costs and three levels of the contractual due dates. Finally, the Standard Error (SE) and the average value (Mean) of each criterion were used for the precise analysis of the classes.

CPUT(s)
Mean MOPSO Mean NSGA-II showed that the NSGA-II was completely superior to the MOPSO algorithm in the five criteria except for the CPUT and diversity of Pareto solutions in most of the generated samples. Moreover, in the medium and large scale sample problems, the results revealed the significantly improved performance and efficiency of the MOPSO algorithm in some performance criteria such as QM, and DM. However, generally, the NSGA-II performed better than the MOPSO algorithm in solving problems in all the tests except for the CPUT and DM. For future research, the corresponding model can be extended with constraints better describing real-life conditions. Additionally, a vehicle routing problem for transportation can be investigated as a more detailed examination of distribution costs. Moreover, other multi-objective algorithms can be used to address the problem and their performance can be evaluated by comparing with the algorithms presented in the current study.