# Multiobjective Parallel Algorithms for Solving Biobjective Open Shop Scheduling Problem

Authors: Seyed Hassan Shams Lahroudi; Farzaneh Mahalleh; Seyedsaeid Mirkamali
Journal: Complexity
 (2022)

---

## Abstract

Open Shop Scheduling Problem (OSSP) is one of the most important scheduling problems in the field of engineering and industry. This kind of problem includesm machines and n jobs, each job contains a certain number of operations, and each operation has a predetermined processing time on its corresponding machine. The order of processing of these operations affects the completion times of all jobs. Therefore, the purpose of OSSP is to achieve a proper order of processing of jobs using specified machines, so that the completion time of all jobs is minimized. In this paper, the strengths and limitations of three methods are evaluated by comparing the results of solving the OSSP in large-scale and small-scale benchmarks. In this case, the minimized completion time and total tardiness are considered the objective functions of the adapted methods. To solve small-scale problems, we adapt a mathematical model called Multiobjective Mixed Linear Programming (MOMILP). To solve large-scale problems, two metaheuristic algorithms including Multiobjective Parallel Genetic Algorithm (MOPGA) and Multiobjective Parallel Simulated Annealing (MOPSA) are adapted. In experimental results, we randomly generated small-scale problems to compare MOMILP with the Genetic Algorithm (GA) and Simulate Annealing (SA). To compare MOPSA and MOPGA with the state of the art and show their strengths and limitations, we use a standard large-scale benchmark. The simulation results of the proposed algorithms show that although the MOPSA algorithm is faster, the MOPGA algorithm is more efficient in achieving optimal solutions for large-scale problems compared with other methods.


---

## Body

## 1. Introduction

Optimal job scheduling problems are classified into single machine scheduling, job shop scheduling, open shop scheduling, flow shop scheduling, and hybrid scheduling problems; where the single machine scheduling is the simplest problem of scheduling. The problem is to schedule the processing ofn jobs with varying processing times on a single machine [1]. In the Job Shop Scheduling Problem (JSSP), the problem is more complex. Here, each job consists of a set of operations (O1, O2, …, On) that need to be processed in a specific order. Each operation has a unique machine on which it must be performed, and only one operation per task may be performed at the moment [2]. A more complex problem than JSSP is the Open Shop Scheduling Problem (OSSP). OSSP involves m machines (M1,M2,…,Mm) to perform n jobs (J1,J2,…,Jn), each job involving m operations. Jth operation of Jj is indicated by Oij which should be processed on machine Mi for pij≥0 units of time. Therefore, we will have a total of n×m operations (O11,O12,…,Onm) [3]. Like OSSP, each job in the Flow Shop Scheduling Problem (FSSP) consists of precisely m operations. However, the ith operation of the job must be executed on the ith machine. No machine is capable of doing many operations simultaneously. The execution time for each action of each task is provided earlier [4]. The most complex problem in the class of job scheduling problems is the hybrid scheduling problem (HSFP). In this problem, n jobs (J1,J2,…,Jn), must be processed in a multistage manufacturing facility in which each step has several concurrent and identical machines [5].This paper addresses the problem of Biobjective Open Shop Scheduling. In case of complexity of the problem, the OSSP is in between the job shop scheduling and the flow shop scheduling problems. However, all the problems in this class are considered as NP-Hard problems. The OSSP has a broad range of applications across different sectors that include timetabling, Satellite communication, health care management, transportation, tourism, computer vision, and many other applications. In the following subsection, we introduce some related works that tried to solve the problem of open shop scheduling.


### 1.1. Related Works

Production scheduling problems have been widely studied in the past years. A detailed review of this type of research work can be found in [6, 7]. Breit et al. [8] solved the OSSP of two machines with a predetermined stop time for one of the machines and proposed a heuristic algorithm. Next, they set some predetermined stop times for both machines and again used heuristic algorithms to approximate the problem. Hsu et al. [9] examined two different periodic maintenance approaches simultaneously in the single-machine scheduling problem. Periodic maintenance approaches considered for their problem are such that the machine stops after a certain period of time (T) or after processing a certain number of jobs (K) in order to perform the maintenance process. Obviously, if either of these conditions occurs earlier, the machine will stop.To quantitatively characterize the production and distribution optimization issue, Fu et al. [10] formulate a mixed-integer programming model to minimize the maximum completion time. To best handle the presented issue, an improved black widow optimization method is created to address the researched problem. In their suggested method, the solution representation, population initiation, procreation, cannibalism, and mutation, as well as a simulated annealing method, are uniquely constructed.Another condition that arises in the presence of machines’ unavailability constraint is the impermissibility of cutting work due to which in the absence of sufficient time to perform an operation and face the machine stop, the desired operation is not loaded on the machine and its processing is postponed until machine stop is finished. Transportation and movement times between machines are other criteria that have been considered in scheduling problems. Strusevich [11] considered a time interval between the completion of the processing of one activity of a job on one machine and the beginning of the next activity of the same job on another machine for two-machine OSSP and due to occurring transportation in real conditions named it transportation time. He also points out that in chemical and metallurgical applications, these transportation times are equivalent to the times of the heating or cooling processes.Recently, Fu et al. [12] offer a stochastic biobjective two-stage open shop scheduling problem that mimics a car repair process where duties are delegated to different third-party organizations with specialized equipment. The optimization issue was formulated by reducing the overall lateness and processing costs, subject to different resource restrictions. In order to tackle the open shop issue, a hybrid multiobjective optimization of migratory birds, along with a genetic operation and discrete event system, is developed by considering problem features.The OSSP is known as an NP-hard problem [13–15] and it is not possible to solve these problems in polynomial time except in small dimensions. Therefore, approximate solution achievements including heuristic and metaheuristic methods can be more efficient than exact methods. Panahi et al. [16] used multiobjective simulated annealing and Ant Colony Optimization (ACO) algorithm to solve Multiobjective Open Shop Scheduling (MOOSS) with the objectives of total delay and the longest completion time. Panahi and Tavakkoli-Moghaddam [17] combined general and multiobjective simulated annealing algorithms with ACO and used it to simultaneous optimization of two objective functions including the longest completion time and total delay in OSSP. Drezner [18] used the Genetic Algorithm (GA) and Simulated Annealing (SA) algorithm for batch scheduling in the case of multifunction parallel machines.Task scheduling has been studied in much different research works with an array of different applications from task scheduling in cloud computing [19] to the flow shops in factories [20, 21]. Most studies on scheduling problems have been conducted under the basic assumption that machines are available on all planning horizons, while this is not the case in the real world. Machines may not be available for reasons such as preventive maintenance, failure, rest periods, and residual work from the previous planning period that should be processed at the beginning of the new period, [22]. In many cases, unavailability times are known in advance, so the decision-maker can make a more efficient decision by considering them. Karimi et al. [23] considered a predetermined time for machine unavailability in the single machine problem and proposed two heuristic methods to solve it.In this paper, an OSSP is addressed in which stops occur at different times on machines and its duration is different but specific for different machines.The intended objective functions include minimizing the longest completion time (Cmax) and the sum of delays, which are presented as a weight combination in a single function. Also, due to the possibility of using different devices to transport different items in a shop, the movement time in a fixed route is different for each job. Also, the working path from one machine to another is considered different from its return path due to the minimization of material flow interference. Therefore, the work transferring matrix is asymmetric between different workstations. In order to solve problems with small dimensions and when the weight of each target is known, the exact approach of the mixed biobjective linear mathematical model is adapted to use as a problem-solving method, and later the results are compared with the results of GA and SA in solving small-scale OSSP. On the contrary, when the weight of each target is unknown, the achievement of Pareto optimal solutions using metaheuristic algorithms including Multiobjective Parallel Genetic Algorithm (MOPGA) and Multiobjective Parallel Simulated Annealing (MOPSA) algorithm are adapted to solve large-scale problems. The parameters of the metaheuristic algorithms are adjusted using the Taguchi experimental design method, as a result of which, algorithms show a state of stability in the face of various problems.


### 1.2. The Contributions of the Study to the Literature

This study has various contributions to the literature from various aspects. These are as follows:(i)
The exact approach of the mixed multiobjective linear mathematical model is adapted to use as a problem-solving method, and later the results are compared with the results of GA and SA in solving small-scale OSSP.(ii)
Multiobjective Parallel Genetic Algorithm (MOPGA) and Multiobjective Parallel Simulated Annealing (MOPSA) methods are adapted to solve large-scale OSSP.(iii)
To improve the exploitation and exploration capabilities of MOPGA, the new candidate generation stage of MOPGA was developed by changing crossover and mutation operators.(iv)
Variations in MOPGA and MOPSA were compared in detail with various heuristic algorithms selected from the literature according to the average RPD, execution time, and best solution criteria.In the following, Section2 presents a mixed linear mathematical programming model. The adapted metaheuristic algorithms and their details are described in Section 3. Section 4 describes the design of experiments and computational evaluation. Section 5 provides a conclusion and future studies.


## 1.1. Related Works

Production scheduling problems have been widely studied in the past years. A detailed review of this type of research work can be found in [6, 7]. Breit et al. [8] solved the OSSP of two machines with a predetermined stop time for one of the machines and proposed a heuristic algorithm. Next, they set some predetermined stop times for both machines and again used heuristic algorithms to approximate the problem. Hsu et al. [9] examined two different periodic maintenance approaches simultaneously in the single-machine scheduling problem. Periodic maintenance approaches considered for their problem are such that the machine stops after a certain period of time (T) or after processing a certain number of jobs (K) in order to perform the maintenance process. Obviously, if either of these conditions occurs earlier, the machine will stop.To quantitatively characterize the production and distribution optimization issue, Fu et al. [10] formulate a mixed-integer programming model to minimize the maximum completion time. To best handle the presented issue, an improved black widow optimization method is created to address the researched problem. In their suggested method, the solution representation, population initiation, procreation, cannibalism, and mutation, as well as a simulated annealing method, are uniquely constructed.Another condition that arises in the presence of machines’ unavailability constraint is the impermissibility of cutting work due to which in the absence of sufficient time to perform an operation and face the machine stop, the desired operation is not loaded on the machine and its processing is postponed until machine stop is finished. Transportation and movement times between machines are other criteria that have been considered in scheduling problems. Strusevich [11] considered a time interval between the completion of the processing of one activity of a job on one machine and the beginning of the next activity of the same job on another machine for two-machine OSSP and due to occurring transportation in real conditions named it transportation time. He also points out that in chemical and metallurgical applications, these transportation times are equivalent to the times of the heating or cooling processes.Recently, Fu et al. [12] offer a stochastic biobjective two-stage open shop scheduling problem that mimics a car repair process where duties are delegated to different third-party organizations with specialized equipment. The optimization issue was formulated by reducing the overall lateness and processing costs, subject to different resource restrictions. In order to tackle the open shop issue, a hybrid multiobjective optimization of migratory birds, along with a genetic operation and discrete event system, is developed by considering problem features.The OSSP is known as an NP-hard problem [13–15] and it is not possible to solve these problems in polynomial time except in small dimensions. Therefore, approximate solution achievements including heuristic and metaheuristic methods can be more efficient than exact methods. Panahi et al. [16] used multiobjective simulated annealing and Ant Colony Optimization (ACO) algorithm to solve Multiobjective Open Shop Scheduling (MOOSS) with the objectives of total delay and the longest completion time. Panahi and Tavakkoli-Moghaddam [17] combined general and multiobjective simulated annealing algorithms with ACO and used it to simultaneous optimization of two objective functions including the longest completion time and total delay in OSSP. Drezner [18] used the Genetic Algorithm (GA) and Simulated Annealing (SA) algorithm for batch scheduling in the case of multifunction parallel machines.Task scheduling has been studied in much different research works with an array of different applications from task scheduling in cloud computing [19] to the flow shops in factories [20, 21]. Most studies on scheduling problems have been conducted under the basic assumption that machines are available on all planning horizons, while this is not the case in the real world. Machines may not be available for reasons such as preventive maintenance, failure, rest periods, and residual work from the previous planning period that should be processed at the beginning of the new period, [22]. In many cases, unavailability times are known in advance, so the decision-maker can make a more efficient decision by considering them. Karimi et al. [23] considered a predetermined time for machine unavailability in the single machine problem and proposed two heuristic methods to solve it.In this paper, an OSSP is addressed in which stops occur at different times on machines and its duration is different but specific for different machines.The intended objective functions include minimizing the longest completion time (Cmax) and the sum of delays, which are presented as a weight combination in a single function. Also, due to the possibility of using different devices to transport different items in a shop, the movement time in a fixed route is different for each job. Also, the working path from one machine to another is considered different from its return path due to the minimization of material flow interference. Therefore, the work transferring matrix is asymmetric between different workstations. In order to solve problems with small dimensions and when the weight of each target is known, the exact approach of the mixed biobjective linear mathematical model is adapted to use as a problem-solving method, and later the results are compared with the results of GA and SA in solving small-scale OSSP. On the contrary, when the weight of each target is unknown, the achievement of Pareto optimal solutions using metaheuristic algorithms including Multiobjective Parallel Genetic Algorithm (MOPGA) and Multiobjective Parallel Simulated Annealing (MOPSA) algorithm are adapted to solve large-scale problems. The parameters of the metaheuristic algorithms are adjusted using the Taguchi experimental design method, as a result of which, algorithms show a state of stability in the face of various problems.


## 1.2. The Contributions of the Study to the Literature

This study has various contributions to the literature from various aspects. These are as follows:(i)
The exact approach of the mixed multiobjective linear mathematical model is adapted to use as a problem-solving method, and later the results are compared with the results of GA and SA in solving small-scale OSSP.(ii)
Multiobjective Parallel Genetic Algorithm (MOPGA) and Multiobjective Parallel Simulated Annealing (MOPSA) methods are adapted to solve large-scale OSSP.(iii)
To improve the exploitation and exploration capabilities of MOPGA, the new candidate generation stage of MOPGA was developed by changing crossover and mutation operators.(iv)
Variations in MOPGA and MOPSA were compared in detail with various heuristic algorithms selected from the literature according to the average RPD, execution time, and best solution criteria.In the following, Section2 presents a mixed linear mathematical programming model. The adapted metaheuristic algorithms and their details are described in Section 3. Section 4 describes the design of experiments and computational evaluation. Section 5 provides a conclusion and future studies.


## 2. Mathematical Model

In this section, the problem is introduced using simplifying assumptions as well as mathematical symbols and then a multiobjective mixed linear programming model is presented and analyzed to accurately solve small-scale problems.To achieve the mathematical model, the available time of each machine is considered as a packet of a certain lengthTi and jobs should be placed in such a way that their total processing time does not exceed the length of the packet. To clarify the issues, see Figure 1, in which Jj indicates jth job in sequence and Bil indicates lth packet of machine i.Figure 1
A sequence of jobs on the machine.


### 2.1. Problem Assumptions

Simplifying assumptions are considered as follows:(i)
Each machine may process up to one task at a time.(ii)
Each machine may process up to one task at a time.(iii)
Jobs can be processed in any order by all or some machines.(iv)
All jobs are available in the shop from the beginning of the planning horizon.(v)
Job cutting is not allowed during processing.(vi)
There is only one machine of each type in the shop (operations of each job are performed by only one specific machine).(vii)
Unavailable times occur at specific times on each machine.The length of available (and unavailable) times of machines is specific and depends on the machine and is fixed during the planning period.


### 2.2. Indexes and Parameters

j, k: indexes of job (j, k = 1,2, …, n)i, h: indexes of machines (i, h = 1,2, …, m)l, l′: indexes of packets (l, l′ = 1,2, …, b)pij: processing time of job j on machine iTrihj: job transferring time from machine i to machine hTi: the length of the available time of machine iti: the length of unavailable time of machine idj: delivery date of job j


### 2.3. Decision Variables

stij: start time of preparing job j on machine iCij: completion time of job j on machine iCj: completion time of job j on the last machine according to the sequence.Tardj: tardiness of completion time of job j from the delivery dateYihj=1if jobjon machineiis processed before machineh0otherwiseXijk=1,if jobjis processed on machineibefore jobk0,otherwiseAijl=1,if jobjis located in packetlof machineiBil0,otherwise


### 2.4. Multiobjective Mixed Linear Programming Model

In this section, using the above assumptions and symbols, a Multiobjective Mixed Linear Programming Model (MOMILP) is adapted to solve the open shop scheduling problem, taking into account the unavailability of machines, preparation and separation times of jobs and machines, and transmission times between machines to simultaneously minimize the criteria of the longest completion time and the sum of delays as a single function. The MOMILP model is as follows:(1)minz=θ1Cmax+θ2∑j=1nTardj,(2)stij+pij=cij∀j,jstij+pij=Trihj∀j,h,j;(3)M1−Yihj≤sthji≠h,(4)≤stik−M1−Aikl−M1−Xijkj≠kstij+pij−M1−Aijl∀i,j,k,l;(5)stij≥l−1Ti+tiAijl∀i,j,l,(6)Cij−M1−Aijl≤l−1ti+l.Ti∀i,j,l,(7)Cij+M1−Aijl≥l−1ti+Ti∀i,j,l,(8)M1−Aijl+M1−Aikl‵+Xijk≥1∀i,j,k,l,l‵j≠k,l<b,l‵(9)Xijk+Xikj=1∀i,k,j;j<n,k<j,(10)Yihj+Yhij=1∀i,j,h;im,hi,(11)∑l=1bAijl=1∀i,j,(12)∑j=1npij.Aijl≤Ti∀i,l,(13)Ci≥Cij∀i,j,(14)Cmax≥Cij∀i,(15)Tardj≥Cj−dj∀j,(16)stij,Cij,Cj,Cmax,Tardj≥0∀i,j,(17)Xijk,Yihj,Aijl∈0,1∀i,h,j,k,l≠h,j≠kEquation (1) is the objective function of the problem, which includes the convex linear combination of the two criteria: (1) the longest completion time and (2) the sum of the delays.θ1 and θ2 are positive coefficients that represent the weight of each target and θ1+θ2=1. The completion time of each operation (Cij) is obtained by equation (2). If job j on machine i is processed before machine h, we have Yihj=1, then the start time of preparing job j on machine h (sthj=1) will be larger than the completion time of that job on machine i, this is shown using equation (3). If job j is processed on machine i before job k (Xijk=1) and both jobs are located in the same packet l (Aijl=Aikl=1), then the start time of preparing job k on machine i (stik=1) will be larger than the completion time of job j on that machine, this criterion is shown in equation (4). Consequently, equation (5) ensures that no preparation takes place when machines are unavailable. Equations (6) and (7) together prevent the completion of operations during machines are unavailable. Equations (8) and (9) simultaneously specify the sequence of jobs on each machine and packets. Equation (10) determines the order of the machines for the operations of each job. Equation (11) ensures that each operation takes place in exactly one packet of each machine. Equation (12) controls the sum of the total processing times of the operations contained in a packet according to the maximum time of that packet, which is the access time of the relevant machine (Ti). The completion time of each job is calculated through equation (13) and according to the objective function. Equations (14) and (15) calculate the longest completion time and the delay time of each operation, respectively. Finally, equations (16) and (17) represent the negative variables and binary variables (0 or 1).


### 2.5. An Example

Here, to clarify the concept of the problem, we provide an example and solve it using the above MOMILP model. This example includes five jobs (n = 5), two machines (m = 2), and five packets for each machine (b = 5). The rest of the parameters are as follows:(18)pij=611121615141115915,Trih1=0020,Trih2=0000,Ti=2616,ti=24,dj,Trih3=0220,Trih4=0020,Trih5=0020Ti=2616,ti=24,djTi=2616,ti=24,dj=1948603728Considering the weight of the targetsθ1=θ2=0.5, the global optimal solution Z = 90.5 is obtained by solving the MOMILP model in Lingo 8 software; the Gantt chart of which is depicted in Figure 2. As can be seen, no operation is performed during the unavailable times of machines, which are shown in gray. Thus, our main goal of modeling this design has been met. In addition, the sum of the total processing times associated with the operations contained in each packet does not exceed the time of that packet. Also, the obtained sequence indicates that the delivery date has been observed so that the jobs with an earlier delivery time are at the beginning of the sequence and the jobs with a later delivery time are at the end of the sequence. Another noteworthy point is the observance of transfer times; as an example, we can point to the distance between the completion of job J1 on machine M2 and starting this job on the next machine (Tr211 = 2).Figure 2
Gantt chart of optimal solving of the example.


### 2.6. Model Analysis

In this section, the sensitivity of the MOMILP model to the number of jobs (n), machines (m), and packets (b) is investigated. To do so, all indexes related to variables and constraints are shown in Tables 1 and 2, respectively. The effect of increasing the number of jobs, machines, and packets on the number of variables and constraints can also be seen in Table 3.Table 1
Number of variables.
VariableNumberZ, Cmax1Cj, TardjnStij, CijmnXijkmn (n − 1)Yijkmn (m − 1)AijlmnbTable 2
Number of constraints.
ConstraintNumber(1)1(2), (11), (13)mn(3)mn (m − 1)(4)mbn (n − 1)(5), (6), (7)mnb(8)mn (n − 1) b (b-1)/2(9)mn (n − 1)/2(10)mn (m-1)/2(12)mb(14), (15)nTable 3
Number of variables and constraints according to the model MOMILP.
Problem size (m.n.b)Number of variablesNumber of constraints2, 4, 4613933, 4, 4956072, 5, 5918263, 5, 514112615, 30, 307156203992610, 30, 301503140821015, 40, 4012541642590110, 40, 402604112854801


## 2.1. Problem Assumptions

Simplifying assumptions are considered as follows:(i)
Each machine may process up to one task at a time.(ii)
Each machine may process up to one task at a time.(iii)
Jobs can be processed in any order by all or some machines.(iv)
All jobs are available in the shop from the beginning of the planning horizon.(v)
Job cutting is not allowed during processing.(vi)
There is only one machine of each type in the shop (operations of each job are performed by only one specific machine).(vii)
Unavailable times occur at specific times on each machine.The length of available (and unavailable) times of machines is specific and depends on the machine and is fixed during the planning period.


## 2.2. Indexes and Parameters

j, k: indexes of job (j, k = 1,2, …, n)i, h: indexes of machines (i, h = 1,2, …, m)l, l′: indexes of packets (l, l′ = 1,2, …, b)pij: processing time of job j on machine iTrihj: job transferring time from machine i to machine hTi: the length of the available time of machine iti: the length of unavailable time of machine idj: delivery date of job j


## 2.3. Decision Variables

stij: start time of preparing job j on machine iCij: completion time of job j on machine iCj: completion time of job j on the last machine according to the sequence.Tardj: tardiness of completion time of job j from the delivery dateYihj=1if jobjon machineiis processed before machineh0otherwiseXijk=1,if jobjis processed on machineibefore jobk0,otherwiseAijl=1,if jobjis located in packetlof machineiBil0,otherwise


## 2.4. Multiobjective Mixed Linear Programming Model

In this section, using the above assumptions and symbols, a Multiobjective Mixed Linear Programming Model (MOMILP) is adapted to solve the open shop scheduling problem, taking into account the unavailability of machines, preparation and separation times of jobs and machines, and transmission times between machines to simultaneously minimize the criteria of the longest completion time and the sum of delays as a single function. The MOMILP model is as follows:(1)minz=θ1Cmax+θ2∑j=1nTardj,(2)stij+pij=cij∀j,jstij+pij=Trihj∀j,h,j;(3)M1−Yihj≤sthji≠h,(4)≤stik−M1−Aikl−M1−Xijkj≠kstij+pij−M1−Aijl∀i,j,k,l;(5)stij≥l−1Ti+tiAijl∀i,j,l,(6)Cij−M1−Aijl≤l−1ti+l.Ti∀i,j,l,(7)Cij+M1−Aijl≥l−1ti+Ti∀i,j,l,(8)M1−Aijl+M1−Aikl‵+Xijk≥1∀i,j,k,l,l‵j≠k,l<b,l‵(9)Xijk+Xikj=1∀i,k,j;j<n,k<j,(10)Yihj+Yhij=1∀i,j,h;im,hi,(11)∑l=1bAijl=1∀i,j,(12)∑j=1npij.Aijl≤Ti∀i,l,(13)Ci≥Cij∀i,j,(14)Cmax≥Cij∀i,(15)Tardj≥Cj−dj∀j,(16)stij,Cij,Cj,Cmax,Tardj≥0∀i,j,(17)Xijk,Yihj,Aijl∈0,1∀i,h,j,k,l≠h,j≠kEquation (1) is the objective function of the problem, which includes the convex linear combination of the two criteria: (1) the longest completion time and (2) the sum of the delays.θ1 and θ2 are positive coefficients that represent the weight of each target and θ1+θ2=1. The completion time of each operation (Cij) is obtained by equation (2). If job j on machine i is processed before machine h, we have Yihj=1, then the start time of preparing job j on machine h (sthj=1) will be larger than the completion time of that job on machine i, this is shown using equation (3). If job j is processed on machine i before job k (Xijk=1) and both jobs are located in the same packet l (Aijl=Aikl=1), then the start time of preparing job k on machine i (stik=1) will be larger than the completion time of job j on that machine, this criterion is shown in equation (4). Consequently, equation (5) ensures that no preparation takes place when machines are unavailable. Equations (6) and (7) together prevent the completion of operations during machines are unavailable. Equations (8) and (9) simultaneously specify the sequence of jobs on each machine and packets. Equation (10) determines the order of the machines for the operations of each job. Equation (11) ensures that each operation takes place in exactly one packet of each machine. Equation (12) controls the sum of the total processing times of the operations contained in a packet according to the maximum time of that packet, which is the access time of the relevant machine (Ti). The completion time of each job is calculated through equation (13) and according to the objective function. Equations (14) and (15) calculate the longest completion time and the delay time of each operation, respectively. Finally, equations (16) and (17) represent the negative variables and binary variables (0 or 1).


## 2.5. An Example

Here, to clarify the concept of the problem, we provide an example and solve it using the above MOMILP model. This example includes five jobs (n = 5), two machines (m = 2), and five packets for each machine (b = 5). The rest of the parameters are as follows:(18)pij=611121615141115915,Trih1=0020,Trih2=0000,Ti=2616,ti=24,dj,Trih3=0220,Trih4=0020,Trih5=0020Ti=2616,ti=24,djTi=2616,ti=24,dj=1948603728Considering the weight of the targetsθ1=θ2=0.5, the global optimal solution Z = 90.5 is obtained by solving the MOMILP model in Lingo 8 software; the Gantt chart of which is depicted in Figure 2. As can be seen, no operation is performed during the unavailable times of machines, which are shown in gray. Thus, our main goal of modeling this design has been met. In addition, the sum of the total processing times associated with the operations contained in each packet does not exceed the time of that packet. Also, the obtained sequence indicates that the delivery date has been observed so that the jobs with an earlier delivery time are at the beginning of the sequence and the jobs with a later delivery time are at the end of the sequence. Another noteworthy point is the observance of transfer times; as an example, we can point to the distance between the completion of job J1 on machine M2 and starting this job on the next machine (Tr211 = 2).Figure 2
Gantt chart of optimal solving of the example.


## 2.6. Model Analysis

In this section, the sensitivity of the MOMILP model to the number of jobs (n), machines (m), and packets (b) is investigated. To do so, all indexes related to variables and constraints are shown in Tables 1 and 2, respectively. The effect of increasing the number of jobs, machines, and packets on the number of variables and constraints can also be seen in Table 3.Table 1
Number of variables.
VariableNumberZ, Cmax1Cj, TardjnStij, CijmnXijkmn (n − 1)Yijkmn (m − 1)AijlmnbTable 2
Number of constraints.
ConstraintNumber(1)1(2), (11), (13)mn(3)mn (m − 1)(4)mbn (n − 1)(5), (6), (7)mnb(8)mn (n − 1) b (b-1)/2(9)mn (n − 1)/2(10)mn (m-1)/2(12)mb(14), (15)nTable 3
Number of variables and constraints according to the model MOMILP.
Problem size (m.n.b)Number of variablesNumber of constraints2, 4, 4613933, 4, 4956072, 5, 5918263, 5, 514112615, 30, 307156203992610, 30, 301503140821015, 40, 4012541642590110, 40, 402604112854801


## 3. The Genetic Algorithm-Based Method

As mentioned in the introduction, the OSSP is known as an NP-hard problem and only small-scale problems can be solved accurately in a reasonable amount of time [24]. Although the proposed mathematical model achieves an accurate solution, its efficiency decreases as the dimensions of the problem increase. Hence, heuristic methods or approximate algorithms are more effective for solving medium- and large-scale problems that commonly occur in the real world.Genetic algorithms (GAs) are a class of general-purpose search strategies based on natural selection and genetics. GAs have been used effectively to a wide array of optimization challenges [19, 25, 26]. In contrast to local search algorithms such as SA, which are often focused on manipulating a single viable solution and are extremely quick, GAs retain and modify a population of possible solutions. Despite the fact that GAs have shown to be a diverse and successful search tool for addressing optimization issues, there are still several scenarios in which the basic GA performs poorly. Therefore, several hybridization methods have been reported in the literature. In this section, we also propose to use the Multiobjective Parallel Genetic Algorithm to solve the OSSP.


### 3.1. Chromosome Display

Chromosomes are in fact encoded solutions and points in solution space. In this paper, each operation is considered as a gene and each chromosome (solution) is represented as a permutation of genes (operations). For example, consider a problem with five jobs and 2 machines, and so 10 operations. A chromosome (or an encoded solution) can be looked at in Table4 for reference.where Oij is the operation of job j on machine i. According to this chromosome, operation O21 will be processed before operation O11 (note that these two operations are related to the same job) and operation O14 will be processed before operation O11 (note that these two operations are related to the same machine).Table 4
An example of a chromosome for a problem with 5 jobs and 2 machines, and 10 operations.
O11O12O13O14O15O21O22O23O24O2536825191047


### 3.2. Initial Population

Scattering in generations prevents rapid convergence and local optimums. Therefore, we randomly generate the initial population (initial generation) to maintain the scattering of solutions in the solution space as much as possible. That is, a random permutation from the set {1, 2, …,m.n} is considered as an individual where m.n indicates the number of operations.


### 3.3. Objective Function

In order to calculate the objective function for each chromosome (solution space points), these so-called coded solutions should be decoded. For this purpose, we place the operations on the machines in ascending order from the first operation to the last one; taking into account their unavailable times as well, we calculate their completion time. Then, taking into account the completion time of all operations, the longest completion time is obtainedCmax=maxi,jCij.Also, in order to calculate the second expression of the objective function, i.e., the sum of the delays, first the completion time of each job is calculated asCmax=Cj=maxiCij and then the delay of that job is obtained by Tardj=max0,Cj−dj. Finally, the sum of delays is calculated as ∑j=1nTardj. How to determine the values of the coefficients θ1 and θ2 is explained in the next sections. Therefore, the objective function of each chromosome is calculated as Z=θ1Cmax+θ2∑j=1nTardj.


### 3.4. Selection

Selection is the process by which individuals present in a generation are selected in pairs for mating. There are various approaches for selection in the literature [27]. In this paper, we use the roulette wheel selection approach, according to which individuals with higher fitness have a better chance of mating and this is in fact the basis of Darwin’s theory. In order to apply the roulette wheel selection, two parameters are required: individual selection probability and individual cumulative probability, which are obtained from equations (18) and (19), respectively:(19)PSi=fi∑j=1popsizefj,(20)CPi=∑j=1iPSj.According to the roulette wheel selection approach, a number is randomly generated in the range [0, 1], then individuali that meets the condition CPi−1<r<CPi is selected as one of the pairs. This process is repeated for the second one.


### 3.5. Crossover Operation

After selecting a pair of parents using one of the selection methods, a genetic crossover operator with the probabilitypc is used to combine the two parents and generate two children. There are several crossover techniques and combining parents: single-point crossover, two-point crossover, multipoint crossover, uniform crossover, and so on [28]. However, the main problem with the crossover operation is that the feasibility of new solutions (newborns) may not be guaranteed. In cases where new chromosomes or newborns produce infeasible solutions, corrective action is usually taken to turn them into feasible solutions, which will lengthen the solving time, by maintaining the feasibility of new chromosomes in each generation and without any modification process [29]. In the following, the crossover operation is described in four steps:Step 1: Select two operations randomly from chromosome 1 (parent).Step 2: Transfer the genes from the two selected operations along with all the genes between them from parent 1 to child 1 while maintaining their location.Step 3: Select the rest of the genes needed by child 1 from left to right until creating a complete chromosome and place them the same way from left to right in the empty spaces of child 1.Step 4: To create the second child, replace the two parents 1 and 2 and repeat steps 2 and 3 with the genes of the same operations selected in step 1.For example, consider the two parents in Tables5–8.Table 5
Chromosomes of parents 1 and 2.
Parent 136825191047Parent 249562171083Table 6
Transfer the genes from the two selected operations along with all the genes between them from parent 1 to child 1 while maintaining their location.Table 7
Select the rest of the genes needed by child 1 from parent 2 and place them in the empty spaces of child 1.Table 8
Exchange the two parents 1 and 2 and repeat steps 2 and 3 until creating child 2.Step 1. Suppose the two operations selected by the parent 1 areO14 and O22, in other words, genes 2 and 9, respectively.Step 2. Transfer genes 2, 5, 1, and 9 from parent 1 to child 1 while maintaining their locations.Step 3. Select the rest of the genes needed by child 1 (from left to right) from parent 2 and place them (from left to right) in the empty spaces of child 1.Step 4. Exchange the two parents 1 and 2 and repeat steps 2 and 3 until creating child 2 as follows:


### 3.6. Mutation

In order to prevent the generated generations from being directed toward the local optimums, a mutation operator with a probabilitypm is applied to each of the produced children. This operator exchanges the two randomly selected operations and their genes [29]. For example, suppose that operations O12 and O24 are selected from child 1. Then, genes 6 and 8 are exchanged and mutated child will be simillar to the result given in Table  9)Table 9
The genes 6 and 8 are exchanged and mutated is generated.


### 3.7. Termination Condition

The GA continues until a termination criterion is met, which ism.n.εs from the beginning and produces new generations [30]. Where ε is a coefficient whose value is adjusted through experiments. In the following sections, we explain how to set it up. Also, the reason for choosing this termination condition is that it gives the algorithm more time to solve larger problems.


### 3.8. Primary Genetic Algorithm

Before presenting the GA for solving multiobjective problems, here we describe the structure and performance of the original GA using the above explanations. The following is how the basic genetic algorithm works in 12 steps:Step 1. Generate popsize individuals (chromosome) randomly as the initial population.Step 2. Calculate the fitness of each individualStep 3. Select a pair of parents from individuals based on the selection strategy.Step 4. Apply crossover operator with the probabilitypc on the parents and generate two new children.Step 5. Apply mutation operator with the probabilitypm on each child.Step 6. Calculate the fitness of childrenStep 7. Repeat steps 3–6 until generating popsize new children.Step 8. Arrange all existing individuals, both old (i.e., the initial population) and new (i.e., generated children), that reach2×popsize according to their fitness.Step 9. Transfer popsize number of individuals with the best fitness to the new generation.Step 10. Check the termination condition, if it is not met, go to step 3; otherwise, step 11.Step 11. Select the best solution (the most fitted one) of the last generation as the answer to the algorithm.Step 12. End.Values related to the parameters of the GA includingpc, pm, popsize, and ε will be determined by designing statistical tests in future sections.


### 3.9. Multiobjective Parallel Genetic Algorithm

There is a lot of research in the literature in which GA has been used to solve multiobjective problems due to its high ability to deal with multiobjective optimization. Since the GAs deal with a population of points (solutions), Pareto optimal multiple solutions can be found in a population of GAs [27]. In the following, we explain how to deal with a multiobjective problem.In the mathematical model MOMILP,θ1 and θ2 are introduced according to the weight of each objective Cmax and ∑j=1nTardj, and the objective function is written as a single function (1). Now, our biobjective problem has become a single-objective problem. However, determining weights is a delicate and sensitive task. To tackle this problem, we can divide the population into several groups, with individuals in each group searching for one of the possible combinations of goal coefficients, while all groups search in parallel for Pareto optimal solutions. In this way, we can turn a multiobjective problem into a single-objective problem with different weights. Therefore, we consider a set of λ different weights with a little difference between the two consecutive weights. Equation (20) shows the set of possible combinations.(21)θ=θ11,θ21,θ12,θ22,…,θ1λ,θ2λ.Similarly, we divide the population intoλ subpopulations and assign weights of the above set to them. Thus, each subpopulation searches for Pareto optimal solutions in a separate path. Subpopulations act independently and unrelatedly as parallel GAs, but their other parameters and operators, including crossover and mutation, the size of the subpopulation, and termination conditions, are all the same.In this way, the scattering of solutions is effectively strengthened. Also, Pareto solutions are stored to keep dominant solutions. This prevents losing dominant solutions during the optimization process. This set is updated frequently to get closer to Pareto’s optimal solution. The process of MOPGA is as follows:Step 1. Generate the initial population randomly, includingλ×popsize membersStep 2. Divide the initial population intoλ subpopulationsStep 3. Create the weights related to the single objective function and assign to the subpopulations(22)θ=sup×ρ,1−sub×ρ:sub=0,1,…,λ−1.For instance, if we setλ=21 and ρ=0.05, we will have:(23)θ=0.0,1.0,0.05,0.95,…,1.0,0.0.Step 4. Perform the initial genetic algorithm process from step 2 for all subpopulations and save the dominant solutionsStep 5. Find the best solution among the dominant onesStep 6. Finish


### 3.10. The Simulated Annealing (SA)-Based Method

Another common metaheuristic algorithm for solving NP-hard problems is Simulated Annealing. This algorithm is inspired by physical systems in which the energy of atoms is reduced through the cooling process to a minimum level [30–32]. SA analyzes the complete surface of the goal function and attempts to optimize it while going both uphill and downhill. Consequently, it is essentially independent of the initial values, which are often a crucial input in traditional algorithms. In addition, it can escape from local optimums and locate the global optimum via uphill and downhill maneuvers. In addition, SA makes fewer rigorous assumptions about the objective function than traditional algorithms. Because of these lenient assumptions, it can cope with objective functions with ridges and plateaus more readily. Finally, it is capable of optimizing objective functions that are not even specified for certain parameter values [33].In this section, we first provide an overview of the SA and then describe the MOPSA.


#### 3.10.1. The Primary Simulated Annealing (SA) Algorithm

SA starts from a solution (energy level) as the basis of an initial temperature (temp 0) and searches in the solution space in the neighborhood of that solution. The search continues until reaching a termination condition, which we set asm.n.φ iterations. Parameter φ is a fixed coefficient in range (0, 1). In this algorithm, we also use the method of GA to encode the problem solutions and introduce chromosomes. We also considered a random permutation of operations as the base solution for the beginning. In order to move from one solution to another, the neighborhood search method described in Algorithm 1 is used. After completing the neighborhood search process, the current temperature is reduced by multiplying it by the parameter μϵ0,1 and the neighborhood search process is repeated.Algorithm 1: Neighborhood search simulated annealing algorithm.
(1)
Choose a gene at random.(2)
If the selected gene is located in the first location of the chromosome, then replace it with the next gene.Otherwise,(3)
If the selected gene is located at the last location of the chromosome replace it with the previous gene.Otherwise,(6)
Evaluate exchanging with the two lateral genes (previous and next).(7)
Replace the selected gene with the lateral gene that gets the most improvement.(8)
End of condition.The algorithm continues until reaching the final temperature (tempf). During iterations and searches, the new solution is accepted if it improves the value of the objective function, otherwise, the new solution is accepted with the probability exp−Δ/tempi, in which tempi is the current temperature and Δ=100×Znew−Zold/Znew. SA is summarized in Algorithm 2. We consider the final temperature as 1 and the value of other parameters of the algorithm including temp,μandφ set by the statistical experiment design method is discussed in Section 4.Algorithm 2: Steps of the primary simulated annealing algorithm.
(1)
Generate a random permutation of operations.(2)
Calculate the value of the objective function.(3)
Set the parameterstemp,tempf,μ,φ.(4)
Check the termination condition, then go to step 5 if not met and go to step 15 if met.(5)
Setrepeat=0.(6)
repeat<Intm.n.φ go to step 7, otherwise, go to step 14.(7)
Call the neighborhood search process.(8)
If the value of the objective function is improved, go to step 9; otherwise go to step 10.(9)
Accept the new solution, update the value of the objective function, and go to step 13.(10)
Select a random number in (0, 1) and put it inrand.(11)
Ifexp−Δ/tempi<rand go to step 12; otherwise, go to step 13.(12)
Accept the new solution, and update the value of the objective function.(13)
Setrepeat=repeat+1.(14)
Modify the new temperature astempnew=μ.tempold and go to step 4.(15)
End.


#### 3.10.2. Multiobjective Parallel Simulated Annealing (SA)

MOPSA process is very similar to the process of MOPGA and both use the same approach to deal with multiobjective optimization. Here, the objective function is considered as a single function (1) and several SA examines the different weights of each target in parallel and independently. The difference between starting MOPSA and MOPGA is that, in MOPGA, the initial population is divided into several subpopulations, but since SA starts with only one chromosome, such a division is not possible. The process of MOPSA is as follows:Step 1. Produceλ chromosomes randomlyStep 2. Produce the weights of the single objective function as equation (22) and assign them to the chromosomesStep 3. Perform the primary SA as shown in Algorithm 2 for all chromosomes from step 2 and store the dominant solutionsStep 4. Find the best solution among the dominant onesStep 5. End


## 3.1. Chromosome Display

Chromosomes are in fact encoded solutions and points in solution space. In this paper, each operation is considered as a gene and each chromosome (solution) is represented as a permutation of genes (operations). For example, consider a problem with five jobs and 2 machines, and so 10 operations. A chromosome (or an encoded solution) can be looked at in Table4 for reference.where Oij is the operation of job j on machine i. According to this chromosome, operation O21 will be processed before operation O11 (note that these two operations are related to the same job) and operation O14 will be processed before operation O11 (note that these two operations are related to the same machine).Table 4
An example of a chromosome for a problem with 5 jobs and 2 machines, and 10 operations.
O11O12O13O14O15O21O22O23O24O2536825191047


## 3.2. Initial Population

Scattering in generations prevents rapid convergence and local optimums. Therefore, we randomly generate the initial population (initial generation) to maintain the scattering of solutions in the solution space as much as possible. That is, a random permutation from the set {1, 2, …,m.n} is considered as an individual where m.n indicates the number of operations.


## 3.3. Objective Function

In order to calculate the objective function for each chromosome (solution space points), these so-called coded solutions should be decoded. For this purpose, we place the operations on the machines in ascending order from the first operation to the last one; taking into account their unavailable times as well, we calculate their completion time. Then, taking into account the completion time of all operations, the longest completion time is obtainedCmax=maxi,jCij.Also, in order to calculate the second expression of the objective function, i.e., the sum of the delays, first the completion time of each job is calculated asCmax=Cj=maxiCij and then the delay of that job is obtained by Tardj=max0,Cj−dj. Finally, the sum of delays is calculated as ∑j=1nTardj. How to determine the values of the coefficients θ1 and θ2 is explained in the next sections. Therefore, the objective function of each chromosome is calculated as Z=θ1Cmax+θ2∑j=1nTardj.


## 3.4. Selection

Selection is the process by which individuals present in a generation are selected in pairs for mating. There are various approaches for selection in the literature [27]. In this paper, we use the roulette wheel selection approach, according to which individuals with higher fitness have a better chance of mating and this is in fact the basis of Darwin’s theory. In order to apply the roulette wheel selection, two parameters are required: individual selection probability and individual cumulative probability, which are obtained from equations (18) and (19), respectively:(19)PSi=fi∑j=1popsizefj,(20)CPi=∑j=1iPSj.According to the roulette wheel selection approach, a number is randomly generated in the range [0, 1], then individuali that meets the condition CPi−1<r<CPi is selected as one of the pairs. This process is repeated for the second one.


## 3.5. Crossover Operation

After selecting a pair of parents using one of the selection methods, a genetic crossover operator with the probabilitypc is used to combine the two parents and generate two children. There are several crossover techniques and combining parents: single-point crossover, two-point crossover, multipoint crossover, uniform crossover, and so on [28]. However, the main problem with the crossover operation is that the feasibility of new solutions (newborns) may not be guaranteed. In cases where new chromosomes or newborns produce infeasible solutions, corrective action is usually taken to turn them into feasible solutions, which will lengthen the solving time, by maintaining the feasibility of new chromosomes in each generation and without any modification process [29]. In the following, the crossover operation is described in four steps:Step 1: Select two operations randomly from chromosome 1 (parent).Step 2: Transfer the genes from the two selected operations along with all the genes between them from parent 1 to child 1 while maintaining their location.Step 3: Select the rest of the genes needed by child 1 from left to right until creating a complete chromosome and place them the same way from left to right in the empty spaces of child 1.Step 4: To create the second child, replace the two parents 1 and 2 and repeat steps 2 and 3 with the genes of the same operations selected in step 1.For example, consider the two parents in Tables5–8.Table 5
Chromosomes of parents 1 and 2.
Parent 136825191047Parent 249562171083Table 6
Transfer the genes from the two selected operations along with all the genes between them from parent 1 to child 1 while maintaining their location.Table 7
Select the rest of the genes needed by child 1 from parent 2 and place them in the empty spaces of child 1.Table 8
Exchange the two parents 1 and 2 and repeat steps 2 and 3 until creating child 2.Step 1. Suppose the two operations selected by the parent 1 areO14 and O22, in other words, genes 2 and 9, respectively.Step 2. Transfer genes 2, 5, 1, and 9 from parent 1 to child 1 while maintaining their locations.Step 3. Select the rest of the genes needed by child 1 (from left to right) from parent 2 and place them (from left to right) in the empty spaces of child 1.Step 4. Exchange the two parents 1 and 2 and repeat steps 2 and 3 until creating child 2 as follows:


## 3.6. Mutation

In order to prevent the generated generations from being directed toward the local optimums, a mutation operator with a probabilitypm is applied to each of the produced children. This operator exchanges the two randomly selected operations and their genes [29]. For example, suppose that operations O12 and O24 are selected from child 1. Then, genes 6 and 8 are exchanged and mutated child will be simillar to the result given in Table  9)Table 9
The genes 6 and 8 are exchanged and mutated is generated.


## 3.7. Termination Condition

The GA continues until a termination criterion is met, which ism.n.εs from the beginning and produces new generations [30]. Where ε is a coefficient whose value is adjusted through experiments. In the following sections, we explain how to set it up. Also, the reason for choosing this termination condition is that it gives the algorithm more time to solve larger problems.


## 3.8. Primary Genetic Algorithm

Before presenting the GA for solving multiobjective problems, here we describe the structure and performance of the original GA using the above explanations. The following is how the basic genetic algorithm works in 12 steps:Step 1. Generate popsize individuals (chromosome) randomly as the initial population.Step 2. Calculate the fitness of each individualStep 3. Select a pair of parents from individuals based on the selection strategy.Step 4. Apply crossover operator with the probabilitypc on the parents and generate two new children.Step 5. Apply mutation operator with the probabilitypm on each child.Step 6. Calculate the fitness of childrenStep 7. Repeat steps 3–6 until generating popsize new children.Step 8. Arrange all existing individuals, both old (i.e., the initial population) and new (i.e., generated children), that reach2×popsize according to their fitness.Step 9. Transfer popsize number of individuals with the best fitness to the new generation.Step 10. Check the termination condition, if it is not met, go to step 3; otherwise, step 11.Step 11. Select the best solution (the most fitted one) of the last generation as the answer to the algorithm.Step 12. End.Values related to the parameters of the GA includingpc, pm, popsize, and ε will be determined by designing statistical tests in future sections.


## 3.9. Multiobjective Parallel Genetic Algorithm

There is a lot of research in the literature in which GA has been used to solve multiobjective problems due to its high ability to deal with multiobjective optimization. Since the GAs deal with a population of points (solutions), Pareto optimal multiple solutions can be found in a population of GAs [27]. In the following, we explain how to deal with a multiobjective problem.In the mathematical model MOMILP,θ1 and θ2 are introduced according to the weight of each objective Cmax and ∑j=1nTardj, and the objective function is written as a single function (1). Now, our biobjective problem has become a single-objective problem. However, determining weights is a delicate and sensitive task. To tackle this problem, we can divide the population into several groups, with individuals in each group searching for one of the possible combinations of goal coefficients, while all groups search in parallel for Pareto optimal solutions. In this way, we can turn a multiobjective problem into a single-objective problem with different weights. Therefore, we consider a set of λ different weights with a little difference between the two consecutive weights. Equation (20) shows the set of possible combinations.(21)θ=θ11,θ21,θ12,θ22,…,θ1λ,θ2λ.Similarly, we divide the population intoλ subpopulations and assign weights of the above set to them. Thus, each subpopulation searches for Pareto optimal solutions in a separate path. Subpopulations act independently and unrelatedly as parallel GAs, but their other parameters and operators, including crossover and mutation, the size of the subpopulation, and termination conditions, are all the same.In this way, the scattering of solutions is effectively strengthened. Also, Pareto solutions are stored to keep dominant solutions. This prevents losing dominant solutions during the optimization process. This set is updated frequently to get closer to Pareto’s optimal solution. The process of MOPGA is as follows:Step 1. Generate the initial population randomly, includingλ×popsize membersStep 2. Divide the initial population intoλ subpopulationsStep 3. Create the weights related to the single objective function and assign to the subpopulations(22)θ=sup×ρ,1−sub×ρ:sub=0,1,…,λ−1.For instance, if we setλ=21 and ρ=0.05, we will have:(23)θ=0.0,1.0,0.05,0.95,…,1.0,0.0.Step 4. Perform the initial genetic algorithm process from step 2 for all subpopulations and save the dominant solutionsStep 5. Find the best solution among the dominant onesStep 6. Finish


## 3.10. The Simulated Annealing (SA)-Based Method

Another common metaheuristic algorithm for solving NP-hard problems is Simulated Annealing. This algorithm is inspired by physical systems in which the energy of atoms is reduced through the cooling process to a minimum level [30–32]. SA analyzes the complete surface of the goal function and attempts to optimize it while going both uphill and downhill. Consequently, it is essentially independent of the initial values, which are often a crucial input in traditional algorithms. In addition, it can escape from local optimums and locate the global optimum via uphill and downhill maneuvers. In addition, SA makes fewer rigorous assumptions about the objective function than traditional algorithms. Because of these lenient assumptions, it can cope with objective functions with ridges and plateaus more readily. Finally, it is capable of optimizing objective functions that are not even specified for certain parameter values [33].In this section, we first provide an overview of the SA and then describe the MOPSA.


### 3.10.1. The Primary Simulated Annealing (SA) Algorithm

SA starts from a solution (energy level) as the basis of an initial temperature (temp 0) and searches in the solution space in the neighborhood of that solution. The search continues until reaching a termination condition, which we set asm.n.φ iterations. Parameter φ is a fixed coefficient in range (0, 1). In this algorithm, we also use the method of GA to encode the problem solutions and introduce chromosomes. We also considered a random permutation of operations as the base solution for the beginning. In order to move from one solution to another, the neighborhood search method described in Algorithm 1 is used. After completing the neighborhood search process, the current temperature is reduced by multiplying it by the parameter μϵ0,1 and the neighborhood search process is repeated.Algorithm 1: Neighborhood search simulated annealing algorithm.
(1)
Choose a gene at random.(2)
If the selected gene is located in the first location of the chromosome, then replace it with the next gene.Otherwise,(3)
If the selected gene is located at the last location of the chromosome replace it with the previous gene.Otherwise,(6)
Evaluate exchanging with the two lateral genes (previous and next).(7)
Replace the selected gene with the lateral gene that gets the most improvement.(8)
End of condition.The algorithm continues until reaching the final temperature (tempf). During iterations and searches, the new solution is accepted if it improves the value of the objective function, otherwise, the new solution is accepted with the probability exp−Δ/tempi, in which tempi is the current temperature and Δ=100×Znew−Zold/Znew. SA is summarized in Algorithm 2. We consider the final temperature as 1 and the value of other parameters of the algorithm including temp,μandφ set by the statistical experiment design method is discussed in Section 4.Algorithm 2: Steps of the primary simulated annealing algorithm.
(1)
Generate a random permutation of operations.(2)
Calculate the value of the objective function.(3)
Set the parameterstemp,tempf,μ,φ.(4)
Check the termination condition, then go to step 5 if not met and go to step 15 if met.(5)
Setrepeat=0.(6)
repeat<Intm.n.φ go to step 7, otherwise, go to step 14.(7)
Call the neighborhood search process.(8)
If the value of the objective function is improved, go to step 9; otherwise go to step 10.(9)
Accept the new solution, update the value of the objective function, and go to step 13.(10)
Select a random number in (0, 1) and put it inrand.(11)
Ifexp−Δ/tempi<rand go to step 12; otherwise, go to step 13.(12)
Accept the new solution, and update the value of the objective function.(13)
Setrepeat=repeat+1.(14)
Modify the new temperature astempnew=μ.tempold and go to step 4.(15)
End.


### 3.10.2. Multiobjective Parallel Simulated Annealing (SA)

MOPSA process is very similar to the process of MOPGA and both use the same approach to deal with multiobjective optimization. Here, the objective function is considered as a single function (1) and several SA examines the different weights of each target in parallel and independently. The difference between starting MOPSA and MOPGA is that, in MOPGA, the initial population is divided into several subpopulations, but since SA starts with only one chromosome, such a division is not possible. The process of MOPSA is as follows:Step 1. Produceλ chromosomes randomlyStep 2. Produce the weights of the single objective function as equation (22) and assign them to the chromosomesStep 3. Perform the primary SA as shown in Algorithm 2 for all chromosomes from step 2 and store the dominant solutionsStep 4. Find the best solution among the dominant onesStep 5. End


## 3.10.1. The Primary Simulated Annealing (SA) Algorithm

SA starts from a solution (energy level) as the basis of an initial temperature (temp 0) and searches in the solution space in the neighborhood of that solution. The search continues until reaching a termination condition, which we set asm.n.φ iterations. Parameter φ is a fixed coefficient in range (0, 1). In this algorithm, we also use the method of GA to encode the problem solutions and introduce chromosomes. We also considered a random permutation of operations as the base solution for the beginning. In order to move from one solution to another, the neighborhood search method described in Algorithm 1 is used. After completing the neighborhood search process, the current temperature is reduced by multiplying it by the parameter μϵ0,1 and the neighborhood search process is repeated.Algorithm 1: Neighborhood search simulated annealing algorithm.
(1)
Choose a gene at random.(2)
If the selected gene is located in the first location of the chromosome, then replace it with the next gene.Otherwise,(3)
If the selected gene is located at the last location of the chromosome replace it with the previous gene.Otherwise,(6)
Evaluate exchanging with the two lateral genes (previous and next).(7)
Replace the selected gene with the lateral gene that gets the most improvement.(8)
End of condition.The algorithm continues until reaching the final temperature (tempf). During iterations and searches, the new solution is accepted if it improves the value of the objective function, otherwise, the new solution is accepted with the probability exp−Δ/tempi, in which tempi is the current temperature and Δ=100×Znew−Zold/Znew. SA is summarized in Algorithm 2. We consider the final temperature as 1 and the value of other parameters of the algorithm including temp,μandφ set by the statistical experiment design method is discussed in Section 4.Algorithm 2: Steps of the primary simulated annealing algorithm.
(1)
Generate a random permutation of operations.(2)
Calculate the value of the objective function.(3)
Set the parameterstemp,tempf,μ,φ.(4)
Check the termination condition, then go to step 5 if not met and go to step 15 if met.(5)
Setrepeat=0.(6)
repeat<Intm.n.φ go to step 7, otherwise, go to step 14.(7)
Call the neighborhood search process.(8)
If the value of the objective function is improved, go to step 9; otherwise go to step 10.(9)
Accept the new solution, update the value of the objective function, and go to step 13.(10)
Select a random number in (0, 1) and put it inrand.(11)
Ifexp−Δ/tempi<rand go to step 12; otherwise, go to step 13.(12)
Accept the new solution, and update the value of the objective function.(13)
Setrepeat=repeat+1.(14)
Modify the new temperature astempnew=μ.tempold and go to step 4.(15)
End.


## 3.10.2. Multiobjective Parallel Simulated Annealing (SA)

MOPSA process is very similar to the process of MOPGA and both use the same approach to deal with multiobjective optimization. Here, the objective function is considered as a single function (1) and several SA examines the different weights of each target in parallel and independently. The difference between starting MOPSA and MOPGA is that, in MOPGA, the initial population is divided into several subpopulations, but since SA starts with only one chromosome, such a division is not possible. The process of MOPSA is as follows:Step 1. Produceλ chromosomes randomlyStep 2. Produce the weights of the single objective function as equation (22) and assign them to the chromosomesStep 3. Perform the primary SA as shown in Algorithm 2 for all chromosomes from step 2 and store the dominant solutionsStep 4. Find the best solution among the dominant onesStep 5. End


## 4. Experimental Design and Computation Evaluation

In this section, we seek to establish the adapted algorithms through experiments and statistical analysis. For this purpose, different levels of the relevant factors are considered and their different combinations are investigated. Therefore, the Taguchi approach is used as an experimental design method in this paper.In addition, after adjusting the parameters of the proposed algorithms using the Taguchi method and by generating random data, two problem groups, one with small dimensions and the other with large dimensions, are randomly generated and the results of the solution methods are compared.We also compare the results of our experiments on a large-scale benchmark with other methods that solved the OSSP.


### 4.1. Taguchi Experiment Design

Taguchi method is one of the statistical analysis methods in which using orthogonal arrays the number of experiments is significantly reduced compared to complete factorial designs [30]. For example, to design an experiment with four factors, each including three levels, a complete factorial design performs 34=81 experiments, but Taguchi reduces them to nine experiments.Taguchi believed that the key to improving quality was to reduce deviations, not just controlling placement within the specified range. Therefore, he tried to reduce the deviations to zero and paid less attention to being on a fixed slope. This is the main reason why the Taguchi method is mentioned as a reliable method [34]. According to the Taguchi method, factors are divided into two categories: controllable factors or signals and uncontrollable factors or noise. In order to reduce the deviations, a fraction called signal to noise (S/N) is defined which indicates the sensitivity of the response variable to uncontrollable factors and shows the deviation around a certain amount of signal. High values of S/N indicate low deviation, which means that controllable factors are more effective than uncontrollable ones. Using the optimal levels of factors obtained by this method, close to the mean value with the least deviation is expected [30].To convert the objective function to the fractionS/N, three types of criteria can be imagined in the Taguchi method: “the less the better,” “the more the better,” and “the closer to the nominal value the better” [30]. According to the measurement criterion of the problem discussed in this, which is expressed as a minimization expression, the “less is better” mode is appropriate. The relevant S/N formula is as follows:(24)SN=−10log10Z2.Also, a well-known performance criterion named RPD is used:(25)RPD=Algsol−MinsolMinsol.100,Where Algsol is the value of an objective function of an example which is obtained from an algorithm of experiment design and Minsol is the minimum value of the objective function of that example which is obtained from all the algorithms of experiments. So, it is obvious that low RPD values are more valuable.


#### 4.1.1. Data Generation

To perform Taguchi experiments, a problem with dimensionsm,n,b=10,5,20 is considered, and five iterations of the same are generated randomly. Also, other literature-inspired [9, 35–37] parameters of the problem are as follows:(i)
Processing times have a continuous uniform distribution in (1, 99).(ii)
Transfer times between machines have a continuous uniform distribution in (1) and (20).(iii)
Consecutive unavailable durations of each machine have a uniform distribution in (1, 50).(iv)
Consecutive machine access times (packet size of each machine) are calculated as follows:(26)Ti=maxa∑j=1npij,maxipij,wherea is selected from 1/5,1/4,1/3 randomly.Delivery time is calculated as follows:(27)dj=∑i=1mpij+∑i=1m∑h=1mTrihjm+U01.n−1.∑i=1mtim.In whichdj is the total time that a job is on all machines. ∑i=1mpij is the average movement between all machines. ∑i=1m∑h=1mTrihj/m is the average unavailable times of all machines.


#### 4.1.2. Adjusting the Parameters of the MOPGA

Factors to be set in the MOPGA algorithm are as follows: crossover probability (pc), mutation probability pm, population size (popsize), and termination criterion coefficient (ε). According to Table 10, three levels are considered for each factor [29, 30]. So, Orthogonal Array Design L9 is the most appropriate case because it has exactly the same number of factors and levels as the experiment and meets all objectives. Combinations of different levels of factors in each experiment of the design L9 for the MOPGA algorithm are shown in Table 11. MOPGA and all experiment designs are programmed in C#.Net visual language and run on a computer with GHz Intel® Core™ and 4 GB ram. According to the results, optimal levels of factors are as follows: pc=0.8, pm=0.2, popsize=30 and ε=0.4.Table 10
Factors of MOPGA algorithm and their levels.
FactorLevel(1)(2)(3)Pc0.40.60.8Pm0.00.10.2popsize203040ϵ0.20.30.4Table 11
Experiments related to the arrayL9 in the MOPGA algorithm.
ExperimentFactorPcPmpopsizeϵ10.40.0200.220.40.1400.330.40.2300.440.60.0400.450.60.1300.260.60.2200.370.80.0300.380.80.1200.490.80.2400.2In order to find the effect of each factor on the response variable, we perform an analysis of variance (ANOVA). It is important to note that in the designed experiment, the degree of error freedom is 0. To deal with this problem, the factors that have the lowest mean squares are considered errors. In the above experiment, factorspc and popsize with mean squares of 0.92 and 0.70, respectively, are of the least mean squares and so have the least effect on the response variable. Therefore, these two factors are considered errors. Analysis of the variance of the S/N fraction is given in Table 12. Therefore, the most effective factors are pm and ε with the percentage of effect of 86.92% and 7.81%, respectively.Table 12
Analysis of variance of theS/N fraction for the factors of the MOPGA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuepopsize2108.5654.2867.0786.9386.930.00ε211.235.616.947.8194.740.05Error (pc+popsize43.240.815.26100Total8123.02


#### 4.1.3. Adjusting the Parameters of the MOPSA

According to the description of the SA, three factors should be adjusted: initial temperature (temp0), temperature reduction coefficient (μ), and the coefficient of the number of iterations per temperature (φ). Considering two levels for each factor, orthogonal array L4 is suitable and its factors and levels are depicted in Table 13. Also, the combinations of different levels of factors in each experiment of design L4 of the MOPSA are shown in Table 14. MOPSA and all experiment designs are programmed in C#.Net and run on a computer with GHz Intel® Core™ and 4 GB ram. Optimal levels of factors are as follows: temp0=100, μ=0.7, and φ=0.3.Table 13
Factors of MOPSA algorithm and their levels.
FactorLevel(1)(2)temp050100μ0.70.8φ0.30.6Table 14
Array experimentsL4 in MOPSA algorithm.
ExperimentFactortemp0μϕ1500.70.32500.80.631000.70.641000.80.3An analysis of variance is performed to find the effect of each factor. Also, to prevent the error degree of freedom from zeroing, we consider the factor with the lowest mean squares as the error. The results indicate that the factorμ with the mean square 0.03 has the lowest effect on the response variable. Analysis of variance of S/N for the factors of the MOPSA algorithm is shown in Table 15 in which factor φ with the percentage of effect of 66.62% has the highest effect on the response variable. As well, temp0 is the second effective factor with a percentage of effect of 21.89.Table 15
Analysis of variance of theS/N fraction for the factors of the MOPSA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuetemp010.170.176.7221.8921.890.23ϕ10.470.4718.4066.6288.510.15Error (μ)10.030.0311.49100Total30.66


### 4.2. Computational Evaluation

In order to evaluate the performance of the MOMILP model and the proposed algorithms, two sets of problems including small-scale problems (n = 4, 5, 6; m = 2, 3) and large-scale problems (n = 10, 20; m = 5, 10) have been generated [26]. Considering one iteration for each of the small-scale problems and five iterations for each of the large-scale problems, we will have a total of 26 problems. The rest of the parameters are considered exactly as in the Data Generation Section. Considering the weights as θ1=θ2=0.5, small-scale problems are solved using the primary GA and SA and the results will be compared to the optimal solutions obtained by the MOMILP model with the mentioned weights. MOMILP model is coded in Lingo 8 and C#.Net is used to code the primary algorithms GA and SA.The performance of the MOMILP model and primary SA and GA are shown in Table16. As seen in the table, although the MOMILP model is a suitable method to solve small-scale problems, by increasing the computational time of the MOMILP model while magnifying the problem, the GA method becomes more efficient with a deviation of 4.85% from the optimum. As well, although the SA algorithm shows a relatively large deviation from the optimal solution, its time saving compared to GA is clear.Table 16
Performance of MOMILP model and primary GA and SA algorithms in dealing with small-scale problems (θ1=θ2=0.5.).
ProblemZDeviation from optimum (%)Solving time (s)m, n, bMOMILPGASAMOMILPGA (%)SA (%)MOMILPGASA12, 4, 8449.74449.74510.5300.0013.521184123, 4, 8347.31350.83668.9001.0292.6060105132, 5, 10492.95492.95726.7100.0047.42145106143, 5, 10382.65466.38737.00021.8892.611951147252, 6, 12848.09859.431415.2301.3466.871386105363, 6, 12—1026.972375.36———387441684Average04.8562.60To solve large-scale problems, MOPGA and MOPSA algorithms are used and each problem is solved five times using each of the algorithms. In addition, in order to compare algorithms, obtained results are converted to relative percentage difference (RPD) criteria and are shown in Table17. Therefore, MOPGA with an average RPD of 0.05% is more efficient and stable than MOPSA in solving large-scale problems.Table 17
Average RPD for algorithms MOPGA and MOPSA in solving large-scale problems.
Problemm.n.bMOPGAMOPSAProblemm.n.bMOPGAMOPSA75, 10, 200.030.69175, 20, 400.050.3985, 10, 200.030.62185, 20, 400.010.4795, 10, 200.071.10195, 20, 400.050.48105, 10, 200.020.76205, 20, 400.040.65115, 10, 200.040.80215, 20, 400.040.541210, 10, 200.080.732210, 20, 400.010.201310, 10, 200.110.802310, 20, 400.030.241410, 10, 200.070.732410, 20, 400.040.281510, 10, 200.060.802510, 20, 400.070.451610, 10, 200.080.612610, 20, 400.060.41Average RPD of MOPGA: 0.05Average RPD of MOPSA: 0.59In addition, MOPGA and MOPSA are programmed in C#.Net programming language, and all algorithms including models MOMILP, GA, SA, MOPGA, and MOPSA are performed on a computer with CPU 3 GHz Intel® Core™ I and 4 GB RAM.In order to show the efficiency of the proposed algorithm (MOPGA), we tested the above algorithm on benchmarks by Taillard [38] and compared it with other metaheuristic algorithms including SA [39], GA [40], Hybrid GA [40], GA [41], Hybrid GA [42], Ant Colony System (ACS) [43], Cuckoo Search (CS) [43], Cat Swarm Optimization (CSO) algorithm [44], Extended Genetic Algorithm Open Shop (EGA_OS) [14], and Bat Algorithm Open Shop (BA_OS) [22]. Table 18 shows the simulation results of the proposed algorithm and the comparable algorithms for different benchmarks. As mentioned in this table, the proposed algorithms were able to achieve the optimal solution in most benchmarks and solve the problem well.Table 18
Simulation results of the proposed algorithm and comparable algorithms for the Taillard [38] benchmarks.
BenchmarksOptimal solutionSAGAHybrid GAHybrid GAGAEGA_OSACSCSCSOBA_OSMOPSAMOPGA4 × 4 − 11931931932131931931931931931931931931934 × 4 − 22362362362402362392392362362362362362364 × 4 − 32712712712932712712712712712712712712714 × 4–42502502502532502502502502522502502502504 × 4 − 52952952953032952952952952952952952952954 × 4 − 61891891892091891891891891891891891891894 × 4 − 72012012012032012012012012012012012012014 × 4 − 82172172172242172172172172172172172172174 × 4 − 92612612612812612612612612612612612612614 × 4 − 102172172172302172172172172172172172172175 × 5 − 13003003013233003013003013013003003003005 × 5 − 22622622622692622632622622622622622622625 × 5 − 33233233313533233353233313353233233233235 × 5 − 4310310N/AN/A3103163103153143103103103105 × 5 − 5326326N/AN/A3263303263313293263263263265 × 5 − 63123123123273123123123173183123123123125 × 5 − 7303303N/AN/A3033083033083053033033033035 × 5 − 8300300N/AN/A3003043003043033003003003005 × 5 − 93533533533733533583533583583533533533535 × 5 − 103263263263413263283263293293263263263267 × 7 − 14354354384474354364354354364354354354357 × 7 − 24434434554544434474434454474434434434437 × 7 − 3468468N/AN/A4684724684794724684684684687 × 7 − 4463463N/AN/A4634634634674664634634634637 × 7 − 5416416N/AN/A4164174164194164164164164167 × 7 − 6451451N/AN/A4514554514604544524514524517 × 7 − 74224224434504224264224354254224224224227 × 7 − 8424424N/AN/A4244244244244244264244244247 × 7 − 94584584654674584584584584584584584594587 × 7 − 1039839840540639839839839839939839839939810 × 10 − 163763766765563763763763863964563764163710 × 10 − 2588588N/AN/A58858858858868858858858858810 × 10 − 3598598N/AN/A59859859859960059959859859810 × 10 − 457757758658157757757757757757757757757710 × 10 − 5640640N/AN/A64064064064064064064064064010 × 10 − 653853855554153853853853853853853853853810 × 10 − 7616616N/AN/A61661661661661661661661661610 × 10 − 8595595N/AN/A59559559559559559559559559510 × 10 − 959559562759859559559559559559559559559510 × 10 − 1059659662360559659659659659659659659659615 × 15 − 193793796793793793793793793793793793793715 × 15 − 2918918N/AN/A91891891891891891891891991815 × 15 − 387187190487187187187187187187187187187115 × 15 − 493493496993493493493493493493493493493415 × 15 − 5946946N/AN/A94694694694694694694694894615 × 15 − 6933933N/AN/A93393393393393393393393393315 × 15 − 7891891N/AN/A89189189189189189189189189115 × 15 − 889389392889389389389389389389389389389315 × 15 − 9899899N/AN/A89989989989990290289991389915 × 15 − 10902902N/AN/A90290290290290290290290290220 × 20 − 1115511551230116511551155115511551155115511551166115520 × 20 − 212411241N/AN/A12411241124212411242124312411260124120 × 20 − 3125712821292125712571257125712571257125712571257125920 × 20 − 412481274N/AN/A12481248124812481248124812481253124820 × 20 − 5125612891315125612561256125612561256125612561256125620 × 20 − 6120412041266120712041204120412041204120412041204120420 × 20 − 712941294N/AN/A12941294129412941295129412941310129420 × 20 − 811691169N/AN/A11731171117311701176117511701173117620 × 20 − 9128913071339128912891289128912891289128912891289128920 × 20 − 101241N/A13071241124112411241N/A12411241124112411241As you can see in Table19, the MOPSA algorithm has less execution time than the other two algorithms and was able to respond in less time. The MOPGA algorithm was able to achieve the optimal answer in less time than the CSO algorithm.Table 19
Execution time of the proposed algorithm and compared algorithms in terms of seconds.
BenchmarksCSOMOPSAMOPGABenchmarksCSOMOPSAMOPGA4 × 4 − 10.0010.0010.00310 × 10 − 1160.762158.514160.2164 × 4 − 20.0010.0010.00710 × 10 − 2474.782472.948474.3094 × 4 − 30.0010.0010.01810 × 10 − 368.56765.78166.1184 × 4 − 40.1560.1010.14810 × 10 − 441.0440.01440.7834 × 4 − 50.1250.0970.12910 × 10 − 5679.018678.091678.6154 × 4 − 60.0030.0010.00910 × 10 − 6126.936122.726123.5524 × 4 − 70.1090.0740.10110 × 10 − 722.87720.28121.2274 × 4 − 80.3590.2310.28710 × 10 − 8375.292371.557374.5914 × 4 − 90.2190.1420.18610 × 10 − 9263.801257.890262.1754 × 4 − 100.1560.0990.12110 × 10 − 10112.65110.03211.7585 × 5 − 10.6100.4510.57315 × 15 − 1250.323247.080248.9035 × 5 − 20.7350.6180.69115 × 15 − 2501.669494.118486.3375 × 5 − 32.2351.8762.11615 × 15 − 3498.766492.016493.7685 × 5 − 48.9066.5698.59215 × 15 − 4112.994109.85411.2435 × 5 − 51.3441.1231.33815 × 15 − 5834.63827.603831.7285 × 5 − 613.81311.09712.58915 × 15 − 6149.578147.932148.0055 × 5 − 714.68813.45813.97615 × 15 − 7867.539863.590564.3295 × 5 − 835.23631.23134.81115 × 15 − 8176.073175.327175.8245 × 5 − 94.0322.6423.62815 × 15 − 9658.615652.187655.4425 × 5 − 102.4532.3212.40415 × 15 − 10891.272887.849890.1197 × 7 − 124.37521.96823.21120 × 20 − 1457.839453.202455.6017 × 7 − 226.06524.37624.89720 × 20 − 2618.409607.114610.1057 × 7 − 351.65847.61250.62620 × 20 − 3340.266337.621338.2907 × 7 −227.313221.874224.52920 × 20 − 4816.467811.479814.4717 × 7 − 5169.751167.496168.00420 × 20 − 5432.391426.719431.7727 × 7 − 6640.466633.012137.83120 × 20 − 6806.064805.431805.7127 × 7 − 7349.208347.890348.21020 × 20 − 7859.499857.386859.1167 × 7 − 821.56720.00120.90120 × 20 − 8695.772694.329694.7077 × 7 − 983.95378.01681.37820 × 20 − 9207.641203.989206.7527 × 7 − 10112.39110.731112.00420 × 20 − 10489.458484.307488.262The accuracy of the proposed algorithm in programming things well on machines is shown in Figure3. This figure shows the Gantt chart for the 15×15−1 benchmark. In this case, a few free times left for the machine shows the accuracy of the proposed algorithm.Figure 3
Gantt chart of the MOPGA for the15×15−1 benchmark.Figure4 shows the fitness diagram of the MOPGA for the 15×15−1 benchmark. In this diagram, the best fitness is getting close to 900 while the walk number crosses 900 walks.Figure 4
Fitness diagram of the MOPGA for the15×15−1 benchmark.Figure5 shows the dispersion diagram for the 15×15−1 benchmark. Population dispersion is always maintained in the MOPGA. As you can see in the above algorithm, premature convergence of chromosomes is prevented.Figure 5
Dispersion diagram of the MOPGA for the15×15−1 benchmark.


## 4.1. Taguchi Experiment Design

Taguchi method is one of the statistical analysis methods in which using orthogonal arrays the number of experiments is significantly reduced compared to complete factorial designs [30]. For example, to design an experiment with four factors, each including three levels, a complete factorial design performs 34=81 experiments, but Taguchi reduces them to nine experiments.Taguchi believed that the key to improving quality was to reduce deviations, not just controlling placement within the specified range. Therefore, he tried to reduce the deviations to zero and paid less attention to being on a fixed slope. This is the main reason why the Taguchi method is mentioned as a reliable method [34]. According to the Taguchi method, factors are divided into two categories: controllable factors or signals and uncontrollable factors or noise. In order to reduce the deviations, a fraction called signal to noise (S/N) is defined which indicates the sensitivity of the response variable to uncontrollable factors and shows the deviation around a certain amount of signal. High values of S/N indicate low deviation, which means that controllable factors are more effective than uncontrollable ones. Using the optimal levels of factors obtained by this method, close to the mean value with the least deviation is expected [30].To convert the objective function to the fractionS/N, three types of criteria can be imagined in the Taguchi method: “the less the better,” “the more the better,” and “the closer to the nominal value the better” [30]. According to the measurement criterion of the problem discussed in this, which is expressed as a minimization expression, the “less is better” mode is appropriate. The relevant S/N formula is as follows:(24)SN=−10log10Z2.Also, a well-known performance criterion named RPD is used:(25)RPD=Algsol−MinsolMinsol.100,Where Algsol is the value of an objective function of an example which is obtained from an algorithm of experiment design and Minsol is the minimum value of the objective function of that example which is obtained from all the algorithms of experiments. So, it is obvious that low RPD values are more valuable.


### 4.1.1. Data Generation

To perform Taguchi experiments, a problem with dimensionsm,n,b=10,5,20 is considered, and five iterations of the same are generated randomly. Also, other literature-inspired [9, 35–37] parameters of the problem are as follows:(i)
Processing times have a continuous uniform distribution in (1, 99).(ii)
Transfer times between machines have a continuous uniform distribution in (1) and (20).(iii)
Consecutive unavailable durations of each machine have a uniform distribution in (1, 50).(iv)
Consecutive machine access times (packet size of each machine) are calculated as follows:(26)Ti=maxa∑j=1npij,maxipij,wherea is selected from 1/5,1/4,1/3 randomly.Delivery time is calculated as follows:(27)dj=∑i=1mpij+∑i=1m∑h=1mTrihjm+U01.n−1.∑i=1mtim.In whichdj is the total time that a job is on all machines. ∑i=1mpij is the average movement between all machines. ∑i=1m∑h=1mTrihj/m is the average unavailable times of all machines.


### 4.1.2. Adjusting the Parameters of the MOPGA

Factors to be set in the MOPGA algorithm are as follows: crossover probability (pc), mutation probability pm, population size (popsize), and termination criterion coefficient (ε). According to Table 10, three levels are considered for each factor [29, 30]. So, Orthogonal Array Design L9 is the most appropriate case because it has exactly the same number of factors and levels as the experiment and meets all objectives. Combinations of different levels of factors in each experiment of the design L9 for the MOPGA algorithm are shown in Table 11. MOPGA and all experiment designs are programmed in C#.Net visual language and run on a computer with GHz Intel® Core™ and 4 GB ram. According to the results, optimal levels of factors are as follows: pc=0.8, pm=0.2, popsize=30 and ε=0.4.Table 10
Factors of MOPGA algorithm and their levels.
FactorLevel(1)(2)(3)Pc0.40.60.8Pm0.00.10.2popsize203040ϵ0.20.30.4Table 11
Experiments related to the arrayL9 in the MOPGA algorithm.
ExperimentFactorPcPmpopsizeϵ10.40.0200.220.40.1400.330.40.2300.440.60.0400.450.60.1300.260.60.2200.370.80.0300.380.80.1200.490.80.2400.2In order to find the effect of each factor on the response variable, we perform an analysis of variance (ANOVA). It is important to note that in the designed experiment, the degree of error freedom is 0. To deal with this problem, the factors that have the lowest mean squares are considered errors. In the above experiment, factorspc and popsize with mean squares of 0.92 and 0.70, respectively, are of the least mean squares and so have the least effect on the response variable. Therefore, these two factors are considered errors. Analysis of the variance of the S/N fraction is given in Table 12. Therefore, the most effective factors are pm and ε with the percentage of effect of 86.92% and 7.81%, respectively.Table 12
Analysis of variance of theS/N fraction for the factors of the MOPGA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuepopsize2108.5654.2867.0786.9386.930.00ε211.235.616.947.8194.740.05Error (pc+popsize43.240.815.26100Total8123.02


### 4.1.3. Adjusting the Parameters of the MOPSA

According to the description of the SA, three factors should be adjusted: initial temperature (temp0), temperature reduction coefficient (μ), and the coefficient of the number of iterations per temperature (φ). Considering two levels for each factor, orthogonal array L4 is suitable and its factors and levels are depicted in Table 13. Also, the combinations of different levels of factors in each experiment of design L4 of the MOPSA are shown in Table 14. MOPSA and all experiment designs are programmed in C#.Net and run on a computer with GHz Intel® Core™ and 4 GB ram. Optimal levels of factors are as follows: temp0=100, μ=0.7, and φ=0.3.Table 13
Factors of MOPSA algorithm and their levels.
FactorLevel(1)(2)temp050100μ0.70.8φ0.30.6Table 14
Array experimentsL4 in MOPSA algorithm.
ExperimentFactortemp0μϕ1500.70.32500.80.631000.70.641000.80.3An analysis of variance is performed to find the effect of each factor. Also, to prevent the error degree of freedom from zeroing, we consider the factor with the lowest mean squares as the error. The results indicate that the factorμ with the mean square 0.03 has the lowest effect on the response variable. Analysis of variance of S/N for the factors of the MOPSA algorithm is shown in Table 15 in which factor φ with the percentage of effect of 66.62% has the highest effect on the response variable. As well, temp0 is the second effective factor with a percentage of effect of 21.89.Table 15
Analysis of variance of theS/N fraction for the factors of the MOPSA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuetemp010.170.176.7221.8921.890.23ϕ10.470.4718.4066.6288.510.15Error (μ)10.030.0311.49100Total30.66


## 4.1.1. Data Generation

To perform Taguchi experiments, a problem with dimensionsm,n,b=10,5,20 is considered, and five iterations of the same are generated randomly. Also, other literature-inspired [9, 35–37] parameters of the problem are as follows:(i)
Processing times have a continuous uniform distribution in (1, 99).(ii)
Transfer times between machines have a continuous uniform distribution in (1) and (20).(iii)
Consecutive unavailable durations of each machine have a uniform distribution in (1, 50).(iv)
Consecutive machine access times (packet size of each machine) are calculated as follows:(26)Ti=maxa∑j=1npij,maxipij,wherea is selected from 1/5,1/4,1/3 randomly.Delivery time is calculated as follows:(27)dj=∑i=1mpij+∑i=1m∑h=1mTrihjm+U01.n−1.∑i=1mtim.In whichdj is the total time that a job is on all machines. ∑i=1mpij is the average movement between all machines. ∑i=1m∑h=1mTrihj/m is the average unavailable times of all machines.


## 4.1.2. Adjusting the Parameters of the MOPGA

Factors to be set in the MOPGA algorithm are as follows: crossover probability (pc), mutation probability pm, population size (popsize), and termination criterion coefficient (ε). According to Table 10, three levels are considered for each factor [29, 30]. So, Orthogonal Array Design L9 is the most appropriate case because it has exactly the same number of factors and levels as the experiment and meets all objectives. Combinations of different levels of factors in each experiment of the design L9 for the MOPGA algorithm are shown in Table 11. MOPGA and all experiment designs are programmed in C#.Net visual language and run on a computer with GHz Intel® Core™ and 4 GB ram. According to the results, optimal levels of factors are as follows: pc=0.8, pm=0.2, popsize=30 and ε=0.4.Table 10
Factors of MOPGA algorithm and their levels.
FactorLevel(1)(2)(3)Pc0.40.60.8Pm0.00.10.2popsize203040ϵ0.20.30.4Table 11
Experiments related to the arrayL9 in the MOPGA algorithm.
ExperimentFactorPcPmpopsizeϵ10.40.0200.220.40.1400.330.40.2300.440.60.0400.450.60.1300.260.60.2200.370.80.0300.380.80.1200.490.80.2400.2In order to find the effect of each factor on the response variable, we perform an analysis of variance (ANOVA). It is important to note that in the designed experiment, the degree of error freedom is 0. To deal with this problem, the factors that have the lowest mean squares are considered errors. In the above experiment, factorspc and popsize with mean squares of 0.92 and 0.70, respectively, are of the least mean squares and so have the least effect on the response variable. Therefore, these two factors are considered errors. Analysis of the variance of the S/N fraction is given in Table 12. Therefore, the most effective factors are pm and ε with the percentage of effect of 86.92% and 7.81%, respectively.Table 12
Analysis of variance of theS/N fraction for the factors of the MOPGA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuepopsize2108.5654.2867.0786.9386.930.00ε211.235.616.947.8194.740.05Error (pc+popsize43.240.815.26100Total8123.02


## 4.1.3. Adjusting the Parameters of the MOPSA

According to the description of the SA, three factors should be adjusted: initial temperature (temp0), temperature reduction coefficient (μ), and the coefficient of the number of iterations per temperature (φ). Considering two levels for each factor, orthogonal array L4 is suitable and its factors and levels are depicted in Table 13. Also, the combinations of different levels of factors in each experiment of design L4 of the MOPSA are shown in Table 14. MOPSA and all experiment designs are programmed in C#.Net and run on a computer with GHz Intel® Core™ and 4 GB ram. Optimal levels of factors are as follows: temp0=100, μ=0.7, and φ=0.3.Table 13
Factors of MOPSA algorithm and their levels.
FactorLevel(1)(2)temp050100μ0.70.8φ0.30.6Table 14
Array experimentsL4 in MOPSA algorithm.
ExperimentFactortemp0μϕ1500.70.32500.80.631000.70.641000.80.3An analysis of variance is performed to find the effect of each factor. Also, to prevent the error degree of freedom from zeroing, we consider the factor with the lowest mean squares as the error. The results indicate that the factorμ with the mean square 0.03 has the lowest effect on the response variable. Analysis of variance of S/N for the factors of the MOPSA algorithm is shown in Table 15 in which factor φ with the percentage of effect of 66.62% has the highest effect on the response variable. As well, temp0 is the second effective factor with a percentage of effect of 21.89.Table 15
Analysis of variance of theS/N fraction for the factors of the MOPSA algorithm.
FactordfSSMSFPercentage of effect (%)Cumulative percentage (%)p-valuetemp010.170.176.7221.8921.890.23ϕ10.470.4718.4066.6288.510.15Error (μ)10.030.0311.49100Total30.66


## 4.2. Computational Evaluation

In order to evaluate the performance of the MOMILP model and the proposed algorithms, two sets of problems including small-scale problems (n = 4, 5, 6; m = 2, 3) and large-scale problems (n = 10, 20; m = 5, 10) have been generated [26]. Considering one iteration for each of the small-scale problems and five iterations for each of the large-scale problems, we will have a total of 26 problems. The rest of the parameters are considered exactly as in the Data Generation Section. Considering the weights as θ1=θ2=0.5, small-scale problems are solved using the primary GA and SA and the results will be compared to the optimal solutions obtained by the MOMILP model with the mentioned weights. MOMILP model is coded in Lingo 8 and C#.Net is used to code the primary algorithms GA and SA.The performance of the MOMILP model and primary SA and GA are shown in Table16. As seen in the table, although the MOMILP model is a suitable method to solve small-scale problems, by increasing the computational time of the MOMILP model while magnifying the problem, the GA method becomes more efficient with a deviation of 4.85% from the optimum. As well, although the SA algorithm shows a relatively large deviation from the optimal solution, its time saving compared to GA is clear.Table 16
Performance of MOMILP model and primary GA and SA algorithms in dealing with small-scale problems (θ1=θ2=0.5.).
ProblemZDeviation from optimum (%)Solving time (s)m, n, bMOMILPGASAMOMILPGA (%)SA (%)MOMILPGASA12, 4, 8449.74449.74510.5300.0013.521184123, 4, 8347.31350.83668.9001.0292.6060105132, 5, 10492.95492.95726.7100.0047.42145106143, 5, 10382.65466.38737.00021.8892.611951147252, 6, 12848.09859.431415.2301.3466.871386105363, 6, 12—1026.972375.36———387441684Average04.8562.60To solve large-scale problems, MOPGA and MOPSA algorithms are used and each problem is solved five times using each of the algorithms. In addition, in order to compare algorithms, obtained results are converted to relative percentage difference (RPD) criteria and are shown in Table17. Therefore, MOPGA with an average RPD of 0.05% is more efficient and stable than MOPSA in solving large-scale problems.Table 17
Average RPD for algorithms MOPGA and MOPSA in solving large-scale problems.
Problemm.n.bMOPGAMOPSAProblemm.n.bMOPGAMOPSA75, 10, 200.030.69175, 20, 400.050.3985, 10, 200.030.62185, 20, 400.010.4795, 10, 200.071.10195, 20, 400.050.48105, 10, 200.020.76205, 20, 400.040.65115, 10, 200.040.80215, 20, 400.040.541210, 10, 200.080.732210, 20, 400.010.201310, 10, 200.110.802310, 20, 400.030.241410, 10, 200.070.732410, 20, 400.040.281510, 10, 200.060.802510, 20, 400.070.451610, 10, 200.080.612610, 20, 400.060.41Average RPD of MOPGA: 0.05Average RPD of MOPSA: 0.59In addition, MOPGA and MOPSA are programmed in C#.Net programming language, and all algorithms including models MOMILP, GA, SA, MOPGA, and MOPSA are performed on a computer with CPU 3 GHz Intel® Core™ I and 4 GB RAM.In order to show the efficiency of the proposed algorithm (MOPGA), we tested the above algorithm on benchmarks by Taillard [38] and compared it with other metaheuristic algorithms including SA [39], GA [40], Hybrid GA [40], GA [41], Hybrid GA [42], Ant Colony System (ACS) [43], Cuckoo Search (CS) [43], Cat Swarm Optimization (CSO) algorithm [44], Extended Genetic Algorithm Open Shop (EGA_OS) [14], and Bat Algorithm Open Shop (BA_OS) [22]. Table 18 shows the simulation results of the proposed algorithm and the comparable algorithms for different benchmarks. As mentioned in this table, the proposed algorithms were able to achieve the optimal solution in most benchmarks and solve the problem well.Table 18
Simulation results of the proposed algorithm and comparable algorithms for the Taillard [38] benchmarks.
BenchmarksOptimal solutionSAGAHybrid GAHybrid GAGAEGA_OSACSCSCSOBA_OSMOPSAMOPGA4 × 4 − 11931931932131931931931931931931931931934 × 4 − 22362362362402362392392362362362362362364 × 4 − 32712712712932712712712712712712712712714 × 4–42502502502532502502502502522502502502504 × 4 − 52952952953032952952952952952952952952954 × 4 − 61891891892091891891891891891891891891894 × 4 − 72012012012032012012012012012012012012014 × 4 − 82172172172242172172172172172172172172174 × 4 − 92612612612812612612612612612612612612614 × 4 − 102172172172302172172172172172172172172175 × 5 − 13003003013233003013003013013003003003005 × 5 − 22622622622692622632622622622622622622625 × 5 − 33233233313533233353233313353233233233235 × 5 − 4310310N/AN/A3103163103153143103103103105 × 5 − 5326326N/AN/A3263303263313293263263263265 × 5 − 63123123123273123123123173183123123123125 × 5 − 7303303N/AN/A3033083033083053033033033035 × 5 − 8300300N/AN/A3003043003043033003003003005 × 5 − 93533533533733533583533583583533533533535 × 5 − 103263263263413263283263293293263263263267 × 7 − 14354354384474354364354354364354354354357 × 7 − 24434434554544434474434454474434434434437 × 7 − 3468468N/AN/A4684724684794724684684684687 × 7 − 4463463N/AN/A4634634634674664634634634637 × 7 − 5416416N/AN/A4164174164194164164164164167 × 7 − 6451451N/AN/A4514554514604544524514524517 × 7 − 74224224434504224264224354254224224224227 × 7 − 8424424N/AN/A4244244244244244264244244247 × 7 − 94584584654674584584584584584584584594587 × 7 − 1039839840540639839839839839939839839939810 × 10 − 163763766765563763763763863964563764163710 × 10 − 2588588N/AN/A58858858858868858858858858810 × 10 − 3598598N/AN/A59859859859960059959859859810 × 10 − 457757758658157757757757757757757757757710 × 10 − 5640640N/AN/A64064064064064064064064064010 × 10 − 653853855554153853853853853853853853853810 × 10 − 7616616N/AN/A61661661661661661661661661610 × 10 − 8595595N/AN/A59559559559559559559559559510 × 10 − 959559562759859559559559559559559559559510 × 10 − 1059659662360559659659659659659659659659615 × 15 − 193793796793793793793793793793793793793715 × 15 − 2918918N/AN/A91891891891891891891891991815 × 15 − 387187190487187187187187187187187187187115 × 15 − 493493496993493493493493493493493493493415 × 15 − 5946946N/AN/A94694694694694694694694894615 × 15 − 6933933N/AN/A93393393393393393393393393315 × 15 − 7891891N/AN/A89189189189189189189189189115 × 15 − 889389392889389389389389389389389389389315 × 15 − 9899899N/AN/A89989989989990290289991389915 × 15 − 10902902N/AN/A90290290290290290290290290220 × 20 − 1115511551230116511551155115511551155115511551166115520 × 20 − 212411241N/AN/A12411241124212411242124312411260124120 × 20 − 3125712821292125712571257125712571257125712571257125920 × 20 − 412481274N/AN/A12481248124812481248124812481253124820 × 20 − 5125612891315125612561256125612561256125612561256125620 × 20 − 6120412041266120712041204120412041204120412041204120420 × 20 − 712941294N/AN/A12941294129412941295129412941310129420 × 20 − 811691169N/AN/A11731171117311701176117511701173117620 × 20 − 9128913071339128912891289128912891289128912891289128920 × 20 − 101241N/A13071241124112411241N/A12411241124112411241As you can see in Table19, the MOPSA algorithm has less execution time than the other two algorithms and was able to respond in less time. The MOPGA algorithm was able to achieve the optimal answer in less time than the CSO algorithm.Table 19
Execution time of the proposed algorithm and compared algorithms in terms of seconds.
BenchmarksCSOMOPSAMOPGABenchmarksCSOMOPSAMOPGA4 × 4 − 10.0010.0010.00310 × 10 − 1160.762158.514160.2164 × 4 − 20.0010.0010.00710 × 10 − 2474.782472.948474.3094 × 4 − 30.0010.0010.01810 × 10 − 368.56765.78166.1184 × 4 − 40.1560.1010.14810 × 10 − 441.0440.01440.7834 × 4 − 50.1250.0970.12910 × 10 − 5679.018678.091678.6154 × 4 − 60.0030.0010.00910 × 10 − 6126.936122.726123.5524 × 4 − 70.1090.0740.10110 × 10 − 722.87720.28121.2274 × 4 − 80.3590.2310.28710 × 10 − 8375.292371.557374.5914 × 4 − 90.2190.1420.18610 × 10 − 9263.801257.890262.1754 × 4 − 100.1560.0990.12110 × 10 − 10112.65110.03211.7585 × 5 − 10.6100.4510.57315 × 15 − 1250.323247.080248.9035 × 5 − 20.7350.6180.69115 × 15 − 2501.669494.118486.3375 × 5 − 32.2351.8762.11615 × 15 − 3498.766492.016493.7685 × 5 − 48.9066.5698.59215 × 15 − 4112.994109.85411.2435 × 5 − 51.3441.1231.33815 × 15 − 5834.63827.603831.7285 × 5 − 613.81311.09712.58915 × 15 − 6149.578147.932148.0055 × 5 − 714.68813.45813.97615 × 15 − 7867.539863.590564.3295 × 5 − 835.23631.23134.81115 × 15 − 8176.073175.327175.8245 × 5 − 94.0322.6423.62815 × 15 − 9658.615652.187655.4425 × 5 − 102.4532.3212.40415 × 15 − 10891.272887.849890.1197 × 7 − 124.37521.96823.21120 × 20 − 1457.839453.202455.6017 × 7 − 226.06524.37624.89720 × 20 − 2618.409607.114610.1057 × 7 − 351.65847.61250.62620 × 20 − 3340.266337.621338.2907 × 7 −227.313221.874224.52920 × 20 − 4816.467811.479814.4717 × 7 − 5169.751167.496168.00420 × 20 − 5432.391426.719431.7727 × 7 − 6640.466633.012137.83120 × 20 − 6806.064805.431805.7127 × 7 − 7349.208347.890348.21020 × 20 − 7859.499857.386859.1167 × 7 − 821.56720.00120.90120 × 20 − 8695.772694.329694.7077 × 7 − 983.95378.01681.37820 × 20 − 9207.641203.989206.7527 × 7 − 10112.39110.731112.00420 × 20 − 10489.458484.307488.262The accuracy of the proposed algorithm in programming things well on machines is shown in Figure3. This figure shows the Gantt chart for the 15×15−1 benchmark. In this case, a few free times left for the machine shows the accuracy of the proposed algorithm.Figure 3
Gantt chart of the MOPGA for the15×15−1 benchmark.Figure4 shows the fitness diagram of the MOPGA for the 15×15−1 benchmark. In this diagram, the best fitness is getting close to 900 while the walk number crosses 900 walks.Figure 4
Fitness diagram of the MOPGA for the15×15−1 benchmark.Figure5 shows the dispersion diagram for the 15×15−1 benchmark. Population dispersion is always maintained in the MOPGA. As you can see in the above algorithm, premature convergence of chromosomes is prevented.Figure 5
Dispersion diagram of the MOPGA for the15×15−1 benchmark.


## 5. Conclusion

In this paper, a multiobjective mixed-integer linear programming model is adapted to solve biobjective OSSP considering unavailable times for machines. Also, “maximum completion time” and “total delays” are considered in optimization criteria simultaneously. The efficiency of the MOMILP model is investigated through a number of small-scale problems and compared with the efficiency and execution time of GA and SA. Although the results indicate the high ability of the MOMILP model to solve small-scale problems, the analyzes performed on the model highlight the need to develop other methods such as heuristic and metaheuristic methods to solve large-scale problems. Therefore, MOPGA and MOPSA algorithms were introduced for multiobjective problems with unknown weights for objectives. Also, the reliability of multiobjective algorithms is evaluated by designing Taguchi experiments. Results show that the MOPGA algorithm is more efficient than the MOPSA algorithm in dealing with large-scale problems. Finally, variations in MOPGA and MOPSA were compared in detail with various heuristic algorithms selected from the literature according to the Average RPD, execution time, and best solution criteria. Our experiments proved the success of the proposed solutions by running MOPGA and MOPSA not only on low- and medium-scaled open shop scheduling samples but also on large-scale open shop scheduling samples.In the following, some areas of development and expanding the present study are mentioned:(i)
Developing other solution methods, including new heuristic and metaheuristic methods to solve large-scale problems and comparison with the algorithms proposed in this paper can be a good idea for making the model more practical.(ii)
Considering more diverse optimization criteria according to the needs of today’s industry can also provide the basis for future studies.(iii)
In this paper, optimizing factors of the algorithms in discrete space are explored. Therefore, another area of development of the present paper can be the optimization of factors in the continuous space.(iv)
A sequence-dependent preparation and separation times can be considered and then a mathematical model can be developed.(v)
The number of vehicles can be used as a constraint on the problem.(vi)
Considering variable speeds for machines seems practical. That is, the processing speed of the machines decreases after performing a certain number of jobs or a certain period of time and then increases again after performing the maintenance process.



---

*Source: 5043058-2022-08-31.xml*