# Schedule-Allocate and Robust Sequencing in Three-Machine Robotic Cell under Breakdowns

Authors: Bahareh Vaisi; Hiwa Farughi; Sadigh Raissi
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The purpose of this paper is to model two problems comprising schedule-allocate (in case of producing identical parts) and sequencing of parts (in case of producing different parts). The first model is used for minimizing the cycle time and operational cost, and the second one for minimizing both the mean and standard deviation of the total production cost as well the cycle time, in an unreliable three-machine robotic cell which confronted with many uncertainty factors. In the current article, mathematical modelling and simulation-based optimization method have been presented to schedule-allocate similar parts and trace the optimal sequence of different parts. Several solution procedures, including epsilon-constraint method and multiobjective particle swarm optimization algorithm, for identical parts case and response surface methodology for different parts case are applied. The results derived from solving numerical examples revealed some advantages in terms of time to attain the optimal solution.


---

## Body

## 1. Introduction

Applying robots upturn the degree of automation and promote system productivity. In a general definition, a robotic cell consists of robot system, computer numerical control (CNC) machines, and other relevant systems. Robots are responsible for loading, unloading, and displacement of the parts within the cell.Because of the arising complexity on mathematical modelling, machine breakdowns so far have relaxed on scheduling and sequencing optimization process while this may result in changing the optimal solutions. Thus, there is noticeably a need for further study in the areas of optimizing scheduling-allocation of similar parts and determining a sequence of parts in multiple-part-type case that would address the real management objectives.Multicriteria scheduling was surveyed in a study and the most important results appeared in [1]. In a study, Gultekin et al. [2] relaxed the flow shop assumption in a three-machine robotic cell, and by considering the processing times of the parts on the machines as decision variables, they proposed new lower bounds for the one-unit and two-unit robot move cycles under the new problem domain for the flow shop-type robot move cycles. In a flow shop system, for producing the parts, each part meets all machines in the same sequence. In a flow shop robotic cell, one-unit cycles were proved to be optimal for three-machine cells. A one-unit robot move cycle is a robot move cycle that by starting the system from an initial state (input buffer), the robot performs each activity exactly once, and when the robot returns to the initial state, the cycle ends. In the one-unit robot move cycle, exactly one part is produced [2]. A bicriteria scheduling model that jointly minimizes the cycle time and total manufacturing cost in a robotic cell of two and three identical CNC machines producing identical parts was presented by Gultekin et al. [3]. They declared that manufacturing costs for the first time is being studied as an objective function in robotic cell scheduling problems. Alagöz and Azizoğlu [4] presented an optimization algorithm to minimize the stability measure considering the efficiency measure in its minimum level for rescheduling in a parallel machine environment under machine eligibility constraints. Minimization of total manufacturing cost and total weighted completion time were two objectives considered in [5, 6]. In some studies, such as [7], a list of qualitative and quantitative objectives was targeted. Again in [8], manufacturing cost and total weighted tardiness as objectives were minimized for a flexible manufacturing system. Restrepo and Balakrishnan [9] presented a fuzzy logic-based methodology. Shortest processing time and weighted earliest due date were the objectives to schedule a robot-centred cell. In [10], cycle time and total robot travel distance in a robotic cell with processing time windows and non-Euclidean travel times were minimized to schedule the cell. In studies [11, 12], the authors addressed a dynamic scheduling with multiobjective optimization in robotic cells. While Mansouri et al. [13] analysed the trade-off between minimizing make span and total energy consumption in a two-machine flow shop, Ghadiri Nejad et al. [14] considered the trade-off between scheduling of cell processes and the production cost. In another study, a fuzzy multiobjective flexible cell scheduling with two objectives including total cost of production system and total delivery tardiness of jobs under green and energy-efficient strategy was addressed by Far et al. [15]. Simultaneous analysis of throughput and cost of quality for inspection scenarios in a two-machine robotic cell was done by Foumani et al. [16]. The authors claimed that studying robotic cells with controllable inspection times is for the first time.Scheduling in most of the researches in the literature was done based on a single criterion. Basically the criterion is cycle time minimization or output rate maximization. In papers [17–19], the criterion was maximizing throughput rate. For robotic cell scheduling problems, a classification was done and discussed based on machine environment, processing restrictions, and objective function and provided in [20]. In [21], robotic cells with stochastic processing times which are common in semiconductor manufacturing were considered. Online scheduling for improving output was offered in this study. Multidegree cyclic robotic flow shop cell scheduling problem in a cell with multiple single gripper robots to maximize the throughput rate was considered in [22]. An ant colony algorithm was proposed to optimize the problem. In [23], the criterion was maximizing total profit in a fully automated robot spot-welding line which has prescheduled breaks (lunch and tea breaks). In this study, the authors designed the robotic line so that tool replacements can be performed during the abovementioned breaks. Makespan minimization for the robotic cell with controllable processing times was investigated in [24], through a novel genetic algorithm. Foumani et al. [25] minimized the partial cycle time for the robotic cell with an inspection process under different inspection scenarios. Cyclic scheduling in an m-machine robotic cell was studied in [26, 27]. A common approach was minimizing cycle time objective which is equal to maximizing the throughput rate. In [28], an improved discrete bacterial foraging optimization algorithm was proposed to minimize the cycle time in a two-machine robotic cell with sequence-dependent setup times. Cyclic scheduling is defined as the production of parts by repeating a specific sequence of robot movements. The purpose of scheduling in [29] was to determine a cyclic scheduling of jobs performed by machines and an operator, in an automated cell which minimizes the cycle time.From the part-type standpoint, there are three different groups of studies on the robotic cell’s part sequencing reported in the literature. In the first group, the problem is solved based on the well-known Gilmore and Gomory algorithm, such as in [30, 31]. In the second group, the problem is formulated based on TSP problem, such as in [17, 32–34]. Studies [35–40] belong to the third group. In this group, the authors employ heuristics or some especial algorithms to solve the problem. In a new study of “Batur and Erol [41],” a metaheuristic approach with simulated annealing tool was applied to schedule a two-machine robotic cell with inconstant processing time of part.Simulation-based optimization has not been applied to solve the problem of determining part input sequence to the machines heretofore.One of the basic assumptions in the previous studies in the field of a robotic cell is that robots and machines never experience breakdown and never require maintenance. But due to the impossibility of this issue in the real world, in order to make more practical the robotic cell scheduling topic, the assumption of no failure in the machines was relaxed. On the contrary, breakdown and maintenance make the machines’ availability subject to random processes of lifetime. Some previous studies focused on maintenance and reliability appraisement in flexible manufacturing systems, such as [42–44]. Moreover, in none of the previous studies, availability of the robotic cell is advised as a constraint because in all previous studies, machine/robot is supposed to be available. Additionally, the review of the related literature reveals that applying design of experiment (DOE), response surface methodology (RSM), and multicriteria decision making (MCDM) methods together as tools to establish a robust sequence for parts’ entrance to robotic manufacturing cell has not been done before. Therefore, to bridge this gap, scenario concept is being developed. Another novelty of the present research partly lies in consideration of scenario concept for parts input sequencing to the robotic cell where S6 cycle is planned as robot’s move cycle in the 3-machine robotic manufacturing cell. S6 represents a specific robot move cycle in a three-machine robotic cell. The complete process of this move cycle will be described in the Problem Statement section.In this study, there is uncertainty in robot’s waiting and cycle time considering machine’s stochastic breakdown. Furthermore, in case of producing dissimilar parts, different sequencing patterns for the parts entrance to the robotic cell make the robotic cell uncertain. The rest of the paper is structured as follows: first, the general information about the problem, definitions, and notations are stated. Then, a mathematical model for solving the problem in cases of producing identical parts, assumptions, and solution procedures for the problem are developed. With the purpose of overcoming different uncertainty, a computer simulation (CS) model in the case of producing different parts is developed. Afterwards, numerical examples for testing the applicability of the mathematical model (identical parts) and simulation-based optimization approach (different parts) are introduced. The results for the two cases of producing similar and dissimilar parts through solution procedures are demonstrated, and in the case of producing different parts, the simulation model is validated, and also RSM and MCDM methods are detailed to optimize robustly the interactions between robotic cell inputs and outputs. Besides, statistical analysis for comparing the solution methods through optimal schedules and discussion based on optimal part input sequences are revealed. Finally, the conclusion and suggestions for upcoming studies are prepared.


## 2. Problem Statement

Assume a common robotic cell towards organizing rapidly fluctuating customer appeals [42]. A group of machines, such as CNCs, which use robot or robots (to transport materials or parts) and do useful work, form a robotic manufacturing cell. Applying a robotic cell could reduce production cost and promote quality features. Such cells, as described before, consist of one or more machines and are supported by robots. A typical in-line three-machine robotic cell is shown in Figure 1.Figure 1
A typical layout for a 3-machine in-line robotic cell.In the production system, there are three identical CNC machines that each of which has no priority in operation. This means that all machines are capable of performing the required operations at the same time. A single gripper robot is responsible for loading and unloading the parts and also doing the movement of fresh material basket to the selected machine and from the machine to the products basket. Single gripper robots can hold one part at a time. Generally, in a 3-machine cell, six cycles, nominated asS1, S2, S3, S4, S5, and S6 cycles, may be applied for part movements [2, 17]. The present study concentrated on the S6 cycle since it is a one-unit robot move cycle, widely known and commonly used [2, 19], but less attention has been paid in the modelling due to more complexity than the other cycles.The overall process ofS6 cycle describes as follows: in the S6 cycle, primarily, the robot takes its place in front of the input buffer. Then, the following operations are followed by the robot successively: (1) picking up part i, (2) moving to the first machine, (3) loading the first machine, (4) moving to the third machine, (5) waiting for completing the part (i − 2) process (if necessary), (6) unloading the part from the third machine, (7) moving the product to the output buffer, (8) loading the output buffer, (9) moving to the second machine, (10) if necessary, waiting until the completion of the process on part (i − 1), (11) unloading from the second machine, (12) moving the part to the third machine, (13) loading the part on the third machine, (14) moving to the first machine, (15) if necessary, waiting until the completion of the process on part i, (16) then unloading the part from the first machine, (17) moving the part to the second machine, (18) loading part i on the second machine, and (19) finally, turning back to the input buffer. A01A34A23A12 implies the activity sequencing of S6 cycle as abovementioned; Apq is the robot’s activity sequence from station p to station q for p = 0, 1, 2, 3 and q = 1, 2, 3, 4 [2, 17].This robotic cell is studied in two situations including the production of identical parts or different parts.


### 2.1. Definitions

In real manufacturing process, scheduling is one of the key factors at the operational level that affects production efficiency, quality, and cost. Optimal scheduling in robotic cells and performing maintenance activities to balance various objectives, including cost and time, taking into account system availability, increase productivity in these manufacturing systems.Basic definitions for scheduling-allocation in the three-machine robotic cell under breakdown are very close to those in two-machine cell. By introducing a control limit to trigger PM actions for each machine (C), maintenance points and maintenance actions such as PM or CM are determined as explained in [45]. In the case of producing different parts as well, the S6 move cycle is measured.


### 2.2. Notations

Notations for scheduling-allocation and sequencing in the three-machine robotic cell are introduced. First, the following parameters are presented besides/instead of their equivalent in [45]:c¨: the processing time for a part on the 3rd machineP: the total processing time for a part (in 3-machine problem P = a + b + c¨)TS6: cycle time based on the S6 robot move cycleti: the processing times of part i based on the percentage of operations done by each machine j;i=1,…,nTo define the sequence of parts in the proposed unreliable robotic cell by minimizing robustly the cycle time (S6 cycle) and total production cost, the total production cost is formulated based on “Selim Akturk and Gurel [46]” description. In many cases such as the current one, tooling cost considered to be a constant cost and has no effects on optimization process. Cycle time is the time a manufacturing system spends on producing one part. There are succeeding parameters and variables in the mathematical model for computing the total production cost:Co: machining cost ($/min)CPM: cost of each preventive maintenance visit ($/each visit) ignoring setup costCTOOL: cost of tool ($/tool)—tool replacement is forbidden in an operational cyclepj: expected down rate of machine jH: observation periodn: the number of produced parts by type based on the percentage of operations done by each machine during the simulation periodN: the number of throughout products in the simulation periodF: total cost ($/times unit)TTF: time to failure (times unit)TTR: time to repair (times unit)λ: failure rateµ: repair rate


## 2.1. Definitions

In real manufacturing process, scheduling is one of the key factors at the operational level that affects production efficiency, quality, and cost. Optimal scheduling in robotic cells and performing maintenance activities to balance various objectives, including cost and time, taking into account system availability, increase productivity in these manufacturing systems.Basic definitions for scheduling-allocation in the three-machine robotic cell under breakdown are very close to those in two-machine cell. By introducing a control limit to trigger PM actions for each machine (C), maintenance points and maintenance actions such as PM or CM are determined as explained in [45]. In the case of producing different parts as well, the S6 move cycle is measured.


## 2.2. Notations

Notations for scheduling-allocation and sequencing in the three-machine robotic cell are introduced. First, the following parameters are presented besides/instead of their equivalent in [45]:c¨: the processing time for a part on the 3rd machineP: the total processing time for a part (in 3-machine problem P = a + b + c¨)TS6: cycle time based on the S6 robot move cycleti: the processing times of part i based on the percentage of operations done by each machine j;i=1,…,nTo define the sequence of parts in the proposed unreliable robotic cell by minimizing robustly the cycle time (S6 cycle) and total production cost, the total production cost is formulated based on “Selim Akturk and Gurel [46]” description. In many cases such as the current one, tooling cost considered to be a constant cost and has no effects on optimization process. Cycle time is the time a manufacturing system spends on producing one part. There are succeeding parameters and variables in the mathematical model for computing the total production cost:Co: machining cost ($/min)CPM: cost of each preventive maintenance visit ($/each visit) ignoring setup costCTOOL: cost of tool ($/tool)—tool replacement is forbidden in an operational cyclepj: expected down rate of machine jH: observation periodn: the number of produced parts by type based on the percentage of operations done by each machine during the simulation periodN: the number of throughout products in the simulation periodF: total cost ($/times unit)TTF: time to failure (times unit)TTR: time to repair (times unit)λ: failure rateµ: repair rate


## 3. Schedule-Allocate of Identical Parts in Three-Machine Robotic Cell

### 3.1. Problem Formulation

According to the notations in Problem Statement, the cycle time formula based on robot’sS6 move pattern comes here. It should be noted that ε is the loading/unloading time and δ is the travel time of the robot between two successive machines:(1)TS6=8ε+12δ+Max0,W1,W2,W3,(2)W1=a−4ε+8δ,W2=b−4ε+8δ,W3=c¨−4ε+8δ.The following model is proposed for scheduling-allocation of parts:(3)MinF=∑j=13Co∑i=1ntiOij+CPMPij+CTOOL,MinTS6=8ε+12δ+W,s.t.,(4)a=∑i=1ntioi1+WR1,(5)b=∑i=1ntioi2+WR2,(6)c¨=∑i=1ntioi3+WR3,(7)W=Max0,a−4ε+8δ+WR1,b−4ε+8δ+WR2,c¨−4ε+8δ+WR3,(8)WRj=tj,PMRj,PM+tj,CMRj,CM−tj,PMRj,PMtj,CMRj,CM,(9)Pij=AtiUij+AWRj+BUijtik−1+BWRjtikOijTS6CPM,(10)Uij=ti∗OijTS6AV∞,(11)AV∞=4µ3+12λµ2+12µλ24µ3+12λµ2+12µλ2+6λ3,(12)Oi1+Oi2+Oi3=1,(13)0<C<H,(14)W≥0,(15)Rj,PM,Rj,CM,Oij∈0,1.The reported literature reveals that many approaches such as [46] focused on minimizing total production cost consisted of machining, tooling, and preventive maintenance costs. In many cases as this study, tooling cost is considered to be constant value and has no effects on optimization process. Minimizing the total production cost is placed as the 1st objective function, and simultaneously minimizing cycle time is known as the 2nd objective function. Here, the preferred objective functions are expressed by equations in (3) respectively.In the proposed model, the robotic cell scheduling constraints are equations (4)–(7) and (14). These equations are provided based on the definitions of the S6 robot move cycle and redefined consistently with this problem.The constraints related to the maintenance are equations (8)–(11), (13), and (15). The explanations of equations (8)–(10), (13), and (15) are based on the descriptions of equations (7), (8), (11)–(13) one-to-one in [45], and redefined in accordance with S6 robot move cycle.As we may know, availability (AVt) is the probability that a component is functioning at time t by ignoring failure and repair history during [0, t]. Accordingly, steady-state availability (AV∞) is calculated based on equation (16):(16)AV∞=limt⟶∞AVt.Equation (11) expressed the steady-state availability for a system entailed three similar components; each one has exponential distribution lifetime, and two maintenance crews have responsibility to repair the failed machine based on exponential distribution repair time without any preemptions. Proof is accessible in [47]. The decision for allocating the operations to the machines is made through equation (12).


### 3.2. Assumptions 1

In addition to the assumptions in [45], by considering three machines in an in-line robotic cell, the basic assumptions for the current study are summed up as follows:(1)
The machines in the robotic cell are CNC machines. They are capable of performing all the operations.(2)
Every part to be processed has a number of operations to be completed on the machines.(3)
There is not any buffer storage between the machines, and each part is either on a machine or being handled by the robot.(4)
Produced parts are identical.(5)
Production time is continuous.(6)
No preemption is allowed in the processing of any operation.(7)
The focus is on theS6 cycle in the 3-machine robotic cell.(8)
Two maintenance crews have responsibility for repairing the failed machines.


### 3.3. Solution Procedures

Since in the process of solving multiobjective problems, rather than optimal points, nondominant set points are reasoned, so as to attain the set, particular procedures should be done. Until now, many methods have been presented to deal with them [1]. One of the prevalent methods is epsilon-constraint, which changes the problem into a one dimensional problem [48]. In this section, epsilon-constraint approach and a multiobjective particle swarm optimization (MOPSO) algorithm are applied to make different sets of nondominated solutions for the proposed model in Problem Formulation.


#### 3.3.1. Epsilon-Constraint Approach

The epsilon-constraint approach represented byεF/TS6 is applied to generate a representation of the Pareto optimal solutions. Accordingly, minimizing F as an objective gave an upper bound on TS6. The cycle time is acquired as a constraint, due to the independency of S6 cycle time to the allocation of operations on the machines (equation (3)). εF/TS6 is considered to determine the appropriate settings for the processing time on the machines by minimizing total cost for S6 cycle time (TS6). According to these settings, it is possible to mention the total cost as a function of the S6 cycle time, which means that F-compose-TS6 of processing time is made. Consequently, for any given total cost/S6 cycle time, determining the corresponding S6 cycle time/total cost and the processing times of the parts on the machines is practicable.In this study, to produce the Pareto optimal solutions, commercial software (GAMS) is utilized and the exact Pareto front of the bicriteria scheduling-allocation problem is addressed. The formulation of epsilon-constraint problem for theS6 cycle is a mixed-integer nonlinear programming (MINLP) model which simultaneously allocates the operations to the machines and establishes processing time values for them by considering breakdowns caused by machine failures. The formulation is presented as follows:(17)∈F\TS6:Min∑j=13∑i=1nCotiOij+CPMPij+CTOOL,s.t.TS6≤EPSM,Equation4,…,Equation15.Lower bound’s quality for a minimization model is a crucial part of finding a best solution. Mostly, relaxing some constraints and solving an easier problem can compute the lower bound [17]. Here, the lower bound commences finding the best attainable value of f1 in the absence of f2 (f1min) and the best attainable value of f2 in the absence of f1 (f2min). These values are the lower bounds for the objective functions (cycle time and total cost) of the problem.Due to the inefficiency of exact methods in terms of computational time, applying metaheuristic algorithms especially to solve big-size problems will be desirable. In the next step, we put MOPSO heuristics based entirely on Pareto dominance concepts to solve the model [49].


#### 3.3.2. MOPSO

Before describing the MOPSO algorithm, the particle swarm optimization (PSO) needs to be presented, shortly. PSO is a computational method to optimize a problem by iteratively improving a candidate solution, concerning a certain measure of quality. The similarity of PSO with evolutionary algorithms makes clear the notion of using a Pareto ranking scheme as a way to extend the approach to handle multiobjective optimization problems, [50].The fitness function value for algorithm is given identically as the determined score for the problem. After generating some predefined particles in each iteration, each particle should be moved based on other particles situations. Two types of solutions must be updated in each iteration of this algorithm. The best permutation experienced by each particle is denoted byPBest, while the best experienced arrangement by all particles is denoted by GBest. Particles movement direction and their final position in each iteration are calculated according to the following equations [51]:(18)Vit+1=wVit+c1r1PBestit−Xit+c2r2GBestt−Xit,(19)Xit+1=Xit+Vit+1,where w is the inertia weight and is employed to control the impact of the previous history of velocities on the current velocity of a given particle. c1 and c2 are the learning factors and are usually constants. To deal with MOPSO, few changes have to be made to PSO. The purpose of PSO is to find one best solution, while objective of MOPSO is finding a set of different solutions (Pareto optimal set). To achieve this, nondominated solutions obtained in each iteration are stored in an archive. The evaluation of archival nondominated solutions is done to guarantee diversity and spacing measures on the Pareto front [52]. Diversity and spacing metrics as comparison metrics are defined in the Results section. The pseudocode of the MOPSO algorithm is divulged in Figure 2.Figure 2
The pseudocode of the MOPSO.


## 3.1. Problem Formulation

According to the notations in Problem Statement, the cycle time formula based on robot’sS6 move pattern comes here. It should be noted that ε is the loading/unloading time and δ is the travel time of the robot between two successive machines:(1)TS6=8ε+12δ+Max0,W1,W2,W3,(2)W1=a−4ε+8δ,W2=b−4ε+8δ,W3=c¨−4ε+8δ.The following model is proposed for scheduling-allocation of parts:(3)MinF=∑j=13Co∑i=1ntiOij+CPMPij+CTOOL,MinTS6=8ε+12δ+W,s.t.,(4)a=∑i=1ntioi1+WR1,(5)b=∑i=1ntioi2+WR2,(6)c¨=∑i=1ntioi3+WR3,(7)W=Max0,a−4ε+8δ+WR1,b−4ε+8δ+WR2,c¨−4ε+8δ+WR3,(8)WRj=tj,PMRj,PM+tj,CMRj,CM−tj,PMRj,PMtj,CMRj,CM,(9)Pij=AtiUij+AWRj+BUijtik−1+BWRjtikOijTS6CPM,(10)Uij=ti∗OijTS6AV∞,(11)AV∞=4µ3+12λµ2+12µλ24µ3+12λµ2+12µλ2+6λ3,(12)Oi1+Oi2+Oi3=1,(13)0<C<H,(14)W≥0,(15)Rj,PM,Rj,CM,Oij∈0,1.The reported literature reveals that many approaches such as [46] focused on minimizing total production cost consisted of machining, tooling, and preventive maintenance costs. In many cases as this study, tooling cost is considered to be constant value and has no effects on optimization process. Minimizing the total production cost is placed as the 1st objective function, and simultaneously minimizing cycle time is known as the 2nd objective function. Here, the preferred objective functions are expressed by equations in (3) respectively.In the proposed model, the robotic cell scheduling constraints are equations (4)–(7) and (14). These equations are provided based on the definitions of the S6 robot move cycle and redefined consistently with this problem.The constraints related to the maintenance are equations (8)–(11), (13), and (15). The explanations of equations (8)–(10), (13), and (15) are based on the descriptions of equations (7), (8), (11)–(13) one-to-one in [45], and redefined in accordance with S6 robot move cycle.As we may know, availability (AVt) is the probability that a component is functioning at time t by ignoring failure and repair history during [0, t]. Accordingly, steady-state availability (AV∞) is calculated based on equation (16):(16)AV∞=limt⟶∞AVt.Equation (11) expressed the steady-state availability for a system entailed three similar components; each one has exponential distribution lifetime, and two maintenance crews have responsibility to repair the failed machine based on exponential distribution repair time without any preemptions. Proof is accessible in [47]. The decision for allocating the operations to the machines is made through equation (12).


## 3.2. Assumptions 1

In addition to the assumptions in [45], by considering three machines in an in-line robotic cell, the basic assumptions for the current study are summed up as follows:(1)
The machines in the robotic cell are CNC machines. They are capable of performing all the operations.(2)
Every part to be processed has a number of operations to be completed on the machines.(3)
There is not any buffer storage between the machines, and each part is either on a machine or being handled by the robot.(4)
Produced parts are identical.(5)
Production time is continuous.(6)
No preemption is allowed in the processing of any operation.(7)
The focus is on theS6 cycle in the 3-machine robotic cell.(8)
Two maintenance crews have responsibility for repairing the failed machines.


## 3.3. Solution Procedures

Since in the process of solving multiobjective problems, rather than optimal points, nondominant set points are reasoned, so as to attain the set, particular procedures should be done. Until now, many methods have been presented to deal with them [1]. One of the prevalent methods is epsilon-constraint, which changes the problem into a one dimensional problem [48]. In this section, epsilon-constraint approach and a multiobjective particle swarm optimization (MOPSO) algorithm are applied to make different sets of nondominated solutions for the proposed model in Problem Formulation.


### 3.3.1. Epsilon-Constraint Approach

The epsilon-constraint approach represented byεF/TS6 is applied to generate a representation of the Pareto optimal solutions. Accordingly, minimizing F as an objective gave an upper bound on TS6. The cycle time is acquired as a constraint, due to the independency of S6 cycle time to the allocation of operations on the machines (equation (3)). εF/TS6 is considered to determine the appropriate settings for the processing time on the machines by minimizing total cost for S6 cycle time (TS6). According to these settings, it is possible to mention the total cost as a function of the S6 cycle time, which means that F-compose-TS6 of processing time is made. Consequently, for any given total cost/S6 cycle time, determining the corresponding S6 cycle time/total cost and the processing times of the parts on the machines is practicable.In this study, to produce the Pareto optimal solutions, commercial software (GAMS) is utilized and the exact Pareto front of the bicriteria scheduling-allocation problem is addressed. The formulation of epsilon-constraint problem for theS6 cycle is a mixed-integer nonlinear programming (MINLP) model which simultaneously allocates the operations to the machines and establishes processing time values for them by considering breakdowns caused by machine failures. The formulation is presented as follows:(17)∈F\TS6:Min∑j=13∑i=1nCotiOij+CPMPij+CTOOL,s.t.TS6≤EPSM,Equation4,…,Equation15.Lower bound’s quality for a minimization model is a crucial part of finding a best solution. Mostly, relaxing some constraints and solving an easier problem can compute the lower bound [17]. Here, the lower bound commences finding the best attainable value of f1 in the absence of f2 (f1min) and the best attainable value of f2 in the absence of f1 (f2min). These values are the lower bounds for the objective functions (cycle time and total cost) of the problem.Due to the inefficiency of exact methods in terms of computational time, applying metaheuristic algorithms especially to solve big-size problems will be desirable. In the next step, we put MOPSO heuristics based entirely on Pareto dominance concepts to solve the model [49].


### 3.3.2. MOPSO

Before describing the MOPSO algorithm, the particle swarm optimization (PSO) needs to be presented, shortly. PSO is a computational method to optimize a problem by iteratively improving a candidate solution, concerning a certain measure of quality. The similarity of PSO with evolutionary algorithms makes clear the notion of using a Pareto ranking scheme as a way to extend the approach to handle multiobjective optimization problems, [50].The fitness function value for algorithm is given identically as the determined score for the problem. After generating some predefined particles in each iteration, each particle should be moved based on other particles situations. Two types of solutions must be updated in each iteration of this algorithm. The best permutation experienced by each particle is denoted byPBest, while the best experienced arrangement by all particles is denoted by GBest. Particles movement direction and their final position in each iteration are calculated according to the following equations [51]:(18)Vit+1=wVit+c1r1PBestit−Xit+c2r2GBestt−Xit,(19)Xit+1=Xit+Vit+1,where w is the inertia weight and is employed to control the impact of the previous history of velocities on the current velocity of a given particle. c1 and c2 are the learning factors and are usually constants. To deal with MOPSO, few changes have to be made to PSO. The purpose of PSO is to find one best solution, while objective of MOPSO is finding a set of different solutions (Pareto optimal set). To achieve this, nondominated solutions obtained in each iteration are stored in an archive. The evaluation of archival nondominated solutions is done to guarantee diversity and spacing measures on the Pareto front [52]. Diversity and spacing metrics as comparison metrics are defined in the Results section. The pseudocode of the MOPSO algorithm is divulged in Figure 2.Figure 2
The pseudocode of the MOPSO.


## 3.3.1. Epsilon-Constraint Approach

The epsilon-constraint approach represented byεF/TS6 is applied to generate a representation of the Pareto optimal solutions. Accordingly, minimizing F as an objective gave an upper bound on TS6. The cycle time is acquired as a constraint, due to the independency of S6 cycle time to the allocation of operations on the machines (equation (3)). εF/TS6 is considered to determine the appropriate settings for the processing time on the machines by minimizing total cost for S6 cycle time (TS6). According to these settings, it is possible to mention the total cost as a function of the S6 cycle time, which means that F-compose-TS6 of processing time is made. Consequently, for any given total cost/S6 cycle time, determining the corresponding S6 cycle time/total cost and the processing times of the parts on the machines is practicable.In this study, to produce the Pareto optimal solutions, commercial software (GAMS) is utilized and the exact Pareto front of the bicriteria scheduling-allocation problem is addressed. The formulation of epsilon-constraint problem for theS6 cycle is a mixed-integer nonlinear programming (MINLP) model which simultaneously allocates the operations to the machines and establishes processing time values for them by considering breakdowns caused by machine failures. The formulation is presented as follows:(17)∈F\TS6:Min∑j=13∑i=1nCotiOij+CPMPij+CTOOL,s.t.TS6≤EPSM,Equation4,…,Equation15.Lower bound’s quality for a minimization model is a crucial part of finding a best solution. Mostly, relaxing some constraints and solving an easier problem can compute the lower bound [17]. Here, the lower bound commences finding the best attainable value of f1 in the absence of f2 (f1min) and the best attainable value of f2 in the absence of f1 (f2min). These values are the lower bounds for the objective functions (cycle time and total cost) of the problem.Due to the inefficiency of exact methods in terms of computational time, applying metaheuristic algorithms especially to solve big-size problems will be desirable. In the next step, we put MOPSO heuristics based entirely on Pareto dominance concepts to solve the model [49].


## 3.3.2. MOPSO

Before describing the MOPSO algorithm, the particle swarm optimization (PSO) needs to be presented, shortly. PSO is a computational method to optimize a problem by iteratively improving a candidate solution, concerning a certain measure of quality. The similarity of PSO with evolutionary algorithms makes clear the notion of using a Pareto ranking scheme as a way to extend the approach to handle multiobjective optimization problems, [50].The fitness function value for algorithm is given identically as the determined score for the problem. After generating some predefined particles in each iteration, each particle should be moved based on other particles situations. Two types of solutions must be updated in each iteration of this algorithm. The best permutation experienced by each particle is denoted byPBest, while the best experienced arrangement by all particles is denoted by GBest. Particles movement direction and their final position in each iteration are calculated according to the following equations [51]:(18)Vit+1=wVit+c1r1PBestit−Xit+c2r2GBestt−Xit,(19)Xit+1=Xit+Vit+1,where w is the inertia weight and is employed to control the impact of the previous history of velocities on the current velocity of a given particle. c1 and c2 are the learning factors and are usually constants. To deal with MOPSO, few changes have to be made to PSO. The purpose of PSO is to find one best solution, while objective of MOPSO is finding a set of different solutions (Pareto optimal set). To achieve this, nondominated solutions obtained in each iteration are stored in an archive. The evaluation of archival nondominated solutions is done to guarantee diversity and spacing measures on the Pareto front [52]. Diversity and spacing metrics as comparison metrics are defined in the Results section. The pseudocode of the MOPSO algorithm is divulged in Figure 2.Figure 2
The pseudocode of the MOPSO.


## 4. Sequence of Different Parts in Three-Machine Robotic Cell

### 4.1. Proposed Simulation Model

By considering the features of main parameters in Problem Statement, the cycle time formula in equation (1) could also be presented in the following equation precisely:(20)TS6=8ε+12δ+W.Parameters and variables were introduced in the Problem Statement section for computing the cost; hence, the total cost per operational cycle could be formulated by the following equation:(21)F=Co∑i=1nnti+HCPM∑j=13Pj+NCTOOL.Currently, the robotic cells require complex choices regarding sequencing topics caused by uncertainty which rises in different parameters such as time to failure, time to repair, and cycle time that can be supported efficiently by CS models. CS has been established since the early 1960s, and it may be one of the most commonly used analytical tools of management science [53]. In order to investigate the robotic cell, we integrate CS with linear programming (LP) optimization approach to control different sources of uncertainties. Uncertainty in sequencing of arriving the parts in addition to randomly distributed cycle time because of considering stochastic machine breakdowns face this robotic cell with a lot of complexity. It is too difficult to optimize such systems only by using the mathematical programming methods. In these cases, there is definitely no approach to help us efficiently without applying CS models. “As discussed by Fu [54],” CS has been proved to be a successful interface between operations research and computer science.For simulation purposes, the enterprise dynamics (ED) tool is utilized to model the sequence of parts in this robotic manufacturing cell. Although using simulation tool in scheduling for the robotic cell has not been done so far, a comprehensive survey of scheduling on flexible manufacturing systems and analyzing future trend that employed simulation techniques as the analysing tool has been done in [55]. Literature review reveals some latest studies on the application of simulation in manufacturing systems such as simulating mass and small batch production [56], incremental sheet forming [57], and product processing in a mixed flow assembly line [58]. Kikolski [59] investigated some production scenarios using simulation models, built in a Tecnomatix Plant Simulation system. Simulation-based metaheuristic optimization was done by Diaz et al. [60]. A press shop was simulated using Arena [61]. To describe the operation of a job shop, a simulation model was developed [62]. Simulation-based optimization framework was the contribution for supporting strategic maintenance in the production system [63]. A combination of simulation and bow-tie modelling was proposed to study the disruption risks in oil and gas supply chains [64]. Additionally, a comprehensive review of papers between 2002 and 2013 concerning simulation in manufacturing has done in [65]. “Mourtzis et al. [66]” through a significant paper surveyed the major milestones in the development of simulation tools. They studied recent industrial practices of simulation, evolution, advances, and its future trends in the context of manufacturing system.Figure3 displays the layout of the CS model for the robotic cell. Here, three different products flow inside the system by the single griper robot based on the S6 cycle in turn, by means of ED elements.Figure 3
Simulation ED layout model for the three-machine robotic cell.In order to reach an unbiased estimation, the simulation model is run for a long time (more than 10000 hours) with 25 replications, and data are cleared after a warm-up period taking 50 hours. In a simulation software such as ED, the system has to warm-up for some time before it reaches a balanced situation. This is described as the warm-up period. Warm-up period is usually equalled to 10% of observation period in Experiment wizard settings. The results are gathered during the observation period. An observation period starts at moment zero (the default value) or when the warm-up period has ended.Hence, due to attaining more consistency in decision-making affected from different sources of uncertainties, four performance measures (PFMs) are defined as the CS model outputs. These set of objective functions look for solutions which minimize fluctuations of the two simulation outputs in terms of their standard deviations.Y1: expected S6 cycle time in terms of mean valueY2: expected operating cost per part in terms of mean valueY1SD: standard deviation of S6 cycle timeY2SD: standard deviation of part operating costNumerous researches advised on robust optimization especially when significant uncertainties arise from different sources. For the first time, “it was introduced by Soyster [67],” the category of robustness to the operations research. Actually, robustness leads to optimal solutions in uncertain circumstances to estimate a moderately consistent solution.


### 4.2. Assumptions 2

Since in robotic cells, scheduling and sequencing were usually done in deterministic atmosphere, few studies have been done on this issue under stochastic environments. Accepted that if machine’s failure/repair rates are constant over time, one may consider that time to failure and repair time deploy from an exponential distribution, otherwise, an alternate distribution such as Weibull, Gamma, or other density functions might be addressed for their modelling. In addition to the first seven assumptions in [68], the basic assumptions for the current study are summarized as follows:(1)
An in-line robotic cell is available with three machines and a single gripper robot(2)
Every part to be processed has a number of operations to be completed on the machines(3)
Sequence of operations followed throughS6 cycle for a 3-machine robotic manufacturing cell.(4)
Machines meet random failure and require maintenance.(5)
Time to failure and Time to repair follow any specified statistical density function with recognized reliable parameters.(6)
The processing time of each part per machine, without being influenced by stopovers, has been determined ultimately.(7)
Manufacturing each part needs a diversity of operations to be done on the machines, by percentage, some of which are done on first machine, the rest are done on the subsequent machines, correspondingly.


## 4.1. Proposed Simulation Model

By considering the features of main parameters in Problem Statement, the cycle time formula in equation (1) could also be presented in the following equation precisely:(20)TS6=8ε+12δ+W.Parameters and variables were introduced in the Problem Statement section for computing the cost; hence, the total cost per operational cycle could be formulated by the following equation:(21)F=Co∑i=1nnti+HCPM∑j=13Pj+NCTOOL.Currently, the robotic cells require complex choices regarding sequencing topics caused by uncertainty which rises in different parameters such as time to failure, time to repair, and cycle time that can be supported efficiently by CS models. CS has been established since the early 1960s, and it may be one of the most commonly used analytical tools of management science [53]. In order to investigate the robotic cell, we integrate CS with linear programming (LP) optimization approach to control different sources of uncertainties. Uncertainty in sequencing of arriving the parts in addition to randomly distributed cycle time because of considering stochastic machine breakdowns face this robotic cell with a lot of complexity. It is too difficult to optimize such systems only by using the mathematical programming methods. In these cases, there is definitely no approach to help us efficiently without applying CS models. “As discussed by Fu [54],” CS has been proved to be a successful interface between operations research and computer science.For simulation purposes, the enterprise dynamics (ED) tool is utilized to model the sequence of parts in this robotic manufacturing cell. Although using simulation tool in scheduling for the robotic cell has not been done so far, a comprehensive survey of scheduling on flexible manufacturing systems and analyzing future trend that employed simulation techniques as the analysing tool has been done in [55]. Literature review reveals some latest studies on the application of simulation in manufacturing systems such as simulating mass and small batch production [56], incremental sheet forming [57], and product processing in a mixed flow assembly line [58]. Kikolski [59] investigated some production scenarios using simulation models, built in a Tecnomatix Plant Simulation system. Simulation-based metaheuristic optimization was done by Diaz et al. [60]. A press shop was simulated using Arena [61]. To describe the operation of a job shop, a simulation model was developed [62]. Simulation-based optimization framework was the contribution for supporting strategic maintenance in the production system [63]. A combination of simulation and bow-tie modelling was proposed to study the disruption risks in oil and gas supply chains [64]. Additionally, a comprehensive review of papers between 2002 and 2013 concerning simulation in manufacturing has done in [65]. “Mourtzis et al. [66]” through a significant paper surveyed the major milestones in the development of simulation tools. They studied recent industrial practices of simulation, evolution, advances, and its future trends in the context of manufacturing system.Figure3 displays the layout of the CS model for the robotic cell. Here, three different products flow inside the system by the single griper robot based on the S6 cycle in turn, by means of ED elements.Figure 3
Simulation ED layout model for the three-machine robotic cell.In order to reach an unbiased estimation, the simulation model is run for a long time (more than 10000 hours) with 25 replications, and data are cleared after a warm-up period taking 50 hours. In a simulation software such as ED, the system has to warm-up for some time before it reaches a balanced situation. This is described as the warm-up period. Warm-up period is usually equalled to 10% of observation period in Experiment wizard settings. The results are gathered during the observation period. An observation period starts at moment zero (the default value) or when the warm-up period has ended.Hence, due to attaining more consistency in decision-making affected from different sources of uncertainties, four performance measures (PFMs) are defined as the CS model outputs. These set of objective functions look for solutions which minimize fluctuations of the two simulation outputs in terms of their standard deviations.Y1: expected S6 cycle time in terms of mean valueY2: expected operating cost per part in terms of mean valueY1SD: standard deviation of S6 cycle timeY2SD: standard deviation of part operating costNumerous researches advised on robust optimization especially when significant uncertainties arise from different sources. For the first time, “it was introduced by Soyster [67],” the category of robustness to the operations research. Actually, robustness leads to optimal solutions in uncertain circumstances to estimate a moderately consistent solution.


## 4.2. Assumptions 2

Since in robotic cells, scheduling and sequencing were usually done in deterministic atmosphere, few studies have been done on this issue under stochastic environments. Accepted that if machine’s failure/repair rates are constant over time, one may consider that time to failure and repair time deploy from an exponential distribution, otherwise, an alternate distribution such as Weibull, Gamma, or other density functions might be addressed for their modelling. In addition to the first seven assumptions in [68], the basic assumptions for the current study are summarized as follows:(1)
An in-line robotic cell is available with three machines and a single gripper robot(2)
Every part to be processed has a number of operations to be completed on the machines(3)
Sequence of operations followed throughS6 cycle for a 3-machine robotic manufacturing cell.(4)
Machines meet random failure and require maintenance.(5)
Time to failure and Time to repair follow any specified statistical density function with recognized reliable parameters.(6)
The processing time of each part per machine, without being influenced by stopovers, has been determined ultimately.(7)
Manufacturing each part needs a diversity of operations to be done on the machines, by percentage, some of which are done on first machine, the rest are done on the subsequent machines, correspondingly.


## 5. Numerical Examples

### 5.1. Identical Parts Case

In this section, three Test Problems comprising different groups of operations with their processing times in a 3-machine robotic cell for producing identical parts are included. The epsilon-constraint and MOPSO methods are applied for the designated Test Problems. The Test Problems are randomly generated and given in Table1. The parameters and user-defined values for the deliberated robotic cell are presented in Table 2. As described in [45], and in addition to the Problem Statement section’s notations, A, B, and K are the coefficients related to PM index function; and tj,PM and tj,CM are sequentially the duration of a PM and CM visit for machine j. The parameters were set after a number of trials, based on the CNC machines’ physical characteristics and their average failure/repair rates, as discussed by respective authors in [46]. It should be noted that there is no tool replacement for the operations while performing S6 cycle.Table 1
Test problems.
Test problemGroup sizeExample no.Processing times1(small size)(5–10)operations110, 8, 7, 4, 3.210, 7, 13, 8, 5, 2, 5.37, 4, 3, 7, 4, 8, 10, 3, 7.2(large size)(10–18)operations410, 8, 7, 4, 3, 7, 4, 8, 10, 3.510, 8, 7, 4, 3, 7, 4, 8, 10, 3, 7, 4, 8.610, 7, 13, 8, 5, 2, 5, 4, 10, 7, 4, 8, 10, 3.710, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7.810, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 5, 2.3(big size)(18–50)operations910, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 5, 2, 3.1010, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 31110, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 3, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 2.1210, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 3, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 2, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 1, 2, 3, 7, 9.Table 2
Features of necessary parameters.
CPM = 35CTOOL = 45B = 80C = 7µ = 2ε = 1tj,PM = 7C0 = 50A = 40K = 4H = 15λ = 3δ = 2tj,CM = 10


### 5.2. Different Parts Case

Two well-defined examples derived from previous works [33, 68] are applied to examine the relevancy to simulation-based optimization approach.Example 1.
Three dissimilar products are manufactured with consistenttis. For the blue part, tB equals 57, the red part’s tRis 84, and the Purple part has 87 minutes processing time (tP). ε and δ, respectively, are 1 and 2 minutes.Example 2.
A robotic cell is contemplated producing a single product with 57 minutes processing time for a blue coloured product. Here,ε and δ are, respectively, 1 and 2 minutes.
The parameters and user-defined values for the measured robotic cell are accessible in Table3. It should be renowned that there is no tool replacement for the operations during cycle.Table 3
Tabulated main parameters estimation.
Preventive cost;CPM in $ per visitMachining cost;Co in $ per minutesTools cost;CTOOL in $Time to failure density functionTime to repair density function355045Negative exponential with mean of 10 hoursNegative exponential with mean of 2 hours


## 5.1. Identical Parts Case

In this section, three Test Problems comprising different groups of operations with their processing times in a 3-machine robotic cell for producing identical parts are included. The epsilon-constraint and MOPSO methods are applied for the designated Test Problems. The Test Problems are randomly generated and given in Table1. The parameters and user-defined values for the deliberated robotic cell are presented in Table 2. As described in [45], and in addition to the Problem Statement section’s notations, A, B, and K are the coefficients related to PM index function; and tj,PM and tj,CM are sequentially the duration of a PM and CM visit for machine j. The parameters were set after a number of trials, based on the CNC machines’ physical characteristics and their average failure/repair rates, as discussed by respective authors in [46]. It should be noted that there is no tool replacement for the operations while performing S6 cycle.Table 1
Test problems.
Test problemGroup sizeExample no.Processing times1(small size)(5–10)operations110, 8, 7, 4, 3.210, 7, 13, 8, 5, 2, 5.37, 4, 3, 7, 4, 8, 10, 3, 7.2(large size)(10–18)operations410, 8, 7, 4, 3, 7, 4, 8, 10, 3.510, 8, 7, 4, 3, 7, 4, 8, 10, 3, 7, 4, 8.610, 7, 13, 8, 5, 2, 5, 4, 10, 7, 4, 8, 10, 3.710, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7.810, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 5, 2.3(big size)(18–50)operations910, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 5, 2, 3.1010, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 31110, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 3, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 2.1210, 7, 13, 8, 5, 2, 5, 4, 10, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 10, 8, 7, 4, 3, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 2, 10, 8, 7, 4, 3, 7, 4, 8, 10, 3, 1, 2, 3, 7, 9.Table 2
Features of necessary parameters.
CPM = 35CTOOL = 45B = 80C = 7µ = 2ε = 1tj,PM = 7C0 = 50A = 40K = 4H = 15λ = 3δ = 2tj,CM = 10


## 5.2. Different Parts Case

Two well-defined examples derived from previous works [33, 68] are applied to examine the relevancy to simulation-based optimization approach.Example 1.
Three dissimilar products are manufactured with consistenttis. For the blue part, tB equals 57, the red part’s tRis 84, and the Purple part has 87 minutes processing time (tP). ε and δ, respectively, are 1 and 2 minutes.Example 2.
A robotic cell is contemplated producing a single product with 57 minutes processing time for a blue coloured product. Here,ε and δ are, respectively, 1 and 2 minutes.
The parameters and user-defined values for the measured robotic cell are accessible in Table3. It should be renowned that there is no tool replacement for the operations during cycle.Table 3
Tabulated main parameters estimation.
Preventive cost;CPM in $ per visitMachining cost;Co in $ per minutesTools cost;CTOOL in $Time to failure density functionTime to repair density function355045Negative exponential with mean of 10 hoursNegative exponential with mean of 2 hours


## 6. Results

### 6.1. Identical Parts Case

#### 6.1.1. Epsilon-Constraint Results

(MinF-Min TS6) is acquired by solving the coded epsilon-constraint with solver GAMS-BARON. The codes were run on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU. Figure 4 distinctively plots twenty nondominated solutions of feasible region for objective functions in Test Problems (Table 1). Then, in Table 4, the upper (U) and lower (L) bounds of the generated solutions are presented. It should be noted that the solution method was run 10 times per example, and the small-size and large-size results were achieved by 1800 seconds average computational time. Because of the effect of the group size on the model’s computation time, producing the parts with more than 18 operations will be more than 3600 seconds computation’s time. So, the computational time was adjusted to 3600 seconds in the third test problem, and then the model was run 10 times per example.Figure 4
Epsilon-constraint results for the first twenty iterations in Test Problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 4
Upper and lower bounds of the generated Pareto front in the Test Problems.
Example no.BoundsCostS6 cycle time1L1676.51732U1721.089702L2550.03932U2599.119703L2735.44832U2813.112704L3255.75734U3327.174725L4200.63140U4237.42786L4866.2944U4940.935827L5201.57547U5242.042858L5631.05649U5706.741879L5705.43950U5746.8896910L8147.49266U8186.1358511L11450.8688U11481.4410712L15781.14116U15809135The allocation of operations to the machines based on the proposed model is achievable as well. As an instance, regarding the machines’ generated degradation level, the allocation of operations to the machines for Example 11 (Table1) is demonstrated in Table 5.Table 5
Allocation of operations to the machines, for Example 11, based on the epsilon-constraint method.
Machine no.123Generated degradation level13.880947.08100433.89859Allocated operations no.3, 6, 7, 16, 19, 20, 21, 22, 23, 25, 26, 28, 29, 30, 31, 32, 351, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, 349, 10, 18Example 11 includes 35 operations to manufacture a product. The best allocation of operations to the machines according to Table5 was revealed on the basis of the machines’ generated degradation level. Consequently, the operations numbers 9, 10, and 18 with processing times 10, 10, and 10, respectively, are allocated to the third machine; the operations numbers 1, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, and 34 with their processing times 10, 7, 8, 5, 4, 8, 7, 4, 3, 7, 8, 3, 7, 10, and 3, respectively, are allocated to the second machine; and the rest of the operations are allocated to the first machine. This allocation leads to an instantaneous (MinF–MinTS6). This result determines the maximum output rate subsequently.The reliability condition in comparison with the proposed model was also tested. Considering no breakdown for machines in the cell (i.e., an available robotic cell), the model was run again using commercial software (GAMS). The results for the first five Examples of Test Problems (Table1) are shown in Table 6.Table 6
Coordinates of the optimal point in the designated Examples for the robotic cell in the case of testing the availability condition.
Example no.12345Coordinates of optimal point:F, TS63245, 344845, 445195, 476445, 558045, 66It is clear from the results that in case of considering a reliable cell (i.e., a robotic cell with reliable machines where no failures occur on them), an optimal point for the bicriteria model comes as a specific value, whereas in case of considering machine failure and repair in the robotic cell, a Pareto optimal (nondominated) solution is attainable.


#### 6.1.2. MOPSO Results

To express the performance of the proposed model for solving the problem, the proposed MOPSO is examined on the designated Test problems (Table1) by means of MATLAB software on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU.(1) Solution Representation. The MOPSO solution chosen for the problem is a string of decimal numbers. The string as a particle consists of processing times of operations to produce part i (ti). The structure of the string is shown in Figure 5.Figure 5
Forming a solution to solve a problem.All solutions are randomly made in harmony with the above construct. This solution representation is divided into three sections. The first section is related to the allocation of operations to the machines which is illustrated in Figure6. Dimensions of this particle in the form of a matrix are equal to (number of machines∗number of operations). In each row, the maximum value of that row’s processing time is determined. For example, in the first row of the matrix, 0.4 is selected, which means that the cell corresponding to that value takes 1, and the rest of the cells in that row takes the value zero. This means that machine 3 is assigned to process operation 1.Figure 6
Solution representation for allocation of operations to the machines.The second and third sections of solution representation are associated with maintenance decisions for the machines. In these particles, to convert the decimal solution to the zero-one solution, if the value is above 0.5, then the solution is equal to 1, otherwise, it is equal to zero. Figures7 and 8 show maintenance decisions. Based on the example in Figure 7, second and third machines need PM actions. The example of Figure 8 shows that first and third machines require CM operations.Figure 7
Solution representation for PM actions.Figure 8
Solution representation for CM actions.(2) Results. Parameter setting for MOPSO is done by hand in a trial-and-error procedure considering convergence analysis. In this model, w is set as 0.5; r1 and r2 are the random variables between 0 and 1; c1 and c2 are 1 and 2, respectively. The repository size is usually set equal to the swarm (population) size. Here, the algorithm is run for a population size of 50 and a repository size of 50 particles. Number of grids per dimension is 5. The values of other parameters are tabulated in Table 7. It should be noted that grid inflation rate (α) and leader selection pressure (β) are the two predefined parameters. If β>α, there is a cumulative tendency to converge towards PBest for the particle. Higher amount of deletion selection pressure (γ) leads to an increase in probability of member deletion. Deletion selection pressure is a predefined parameter as well.Table 7
MOPSO parameter setting.
Population size = repository sizeαβγNumber of gridsMutation rateMaximum number of iterations500.12250.1100Accordingly, the model is analysed using the MOPSO algorithm, and its results are provided in Table8. In Figure 9, the best results of feasible region for objective functions in Test Problems (Table 1) are represented. Close to the epsilon-constraint, the allocation of operations to the machines for Example 11 (Table 1) is established in Table 9.Table 8
Summary of the objective functions results and elapsed time for solution methods.
Example no.CriteriaMinimum costMinimumS6 cycle timeElapsed time (in seconds)Epsilon-constraintMOPSOEpsilon-constraintMOPSOEpsilon-constraintMOPSO11676.5172281.962323210.8718.93522550.0393451.423323212.3418.58332735.4483871.081323213.9718.46943255.7574556.381343414.0418.95254200.6315909.7424041598.9516.70264866.296695.4024444620.6718.61475201.5757181.277474790017.16885631.0567802.1394949180019.14695705.4398085.5075050360028.957108147.49211245.8196667360026.8491111450.8616030.4278888360028.9431215781.1422353.457116116360032.526Figure 9
MOPSO results of feasible region for objective functions in Test problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 9
Allocation of operations to the machines, for Example 11, based on the MOPSO algorithm.
Machine no.123Allocated operations no.6, 10, 13, 16, 17, 18, 20, 21, 24, 25, 28, 29, 353, 7, 8, 9, 11, 22, 23, 26, 30, 321, 2, 4, 5, 12, 14, 15, 19, 27, 31, 33, 34In Example 11, the best allocation of operations to the machines according to Table9 and the reason of the machines’ generated degradation level was disclosed. With the declared allocation, minimum cost and cycle time simultaneously will be practicable, and also maximum output rate is reasonable.A comparison between the solution methods for Example 11, based on Table8 shows that, however, minimum cycle time is 88 based on both solution methods, and minimum cost was achieved regarding epsilon-constraint method with the amount of 11450.86 $/min. In other Examples of this Test Problem also, minimum cost is realisable considering epsilon-constraint approach. However, MOPSO’s elapsed time is less than that of epsilon-constraint according to this Example.(3) Comparison Metrics. Diversity and spacing metrics are provided for assessing the relative performance of MOPSO. The metrics definition is summarized in Table 10, where xi–yi is the Euclidean distance between two nondominated solutions xi and yi [69]. Comparison metrics were calculated for Test problems (Table 1) over 10 runs of the MOPSO algorithm. Table 11 presents the metrics obtained for small-, large-, and big-sized instances.Table 10
Comparison metrics.
Spacing metric (SM)d⇀=∑i=1ndi/nDiversity metric (DM)DM=∑i=1nmaxxi−yiSM=∑i=1nd⇀−di2/n−1Table 11
Computational results for comparison metrics in Test Problems.
Test ProblemExample no.DMSMSmall size1177.21502157.99903213.0850Large size4273.80505283.44906320.88107286.04308320.690Big size9326.375010386.960011475.863012568.5890


#### 6.1.3. Statistical Analysis

Objective functions and elapsed time required by each method are defined as criteria in Table8 so that the outcome of solution methods will be compared. To evaluate the performance of solution methods, epsilon-constraint method’s minimum cost was measured against that of MOPSO algorithm. F (i.e., cost), in exact solution, is consistently less than metaheuristic solution, regarding Figure 10. In the considered Test Problems (Table 1), difference between the objective function values is at least 605 and at most 6572 currencies.Figure 10
Comparison ofF (cost) based on the solution procedures in Test Problems.A statistical significant difference between the solution methods is not observed forTS6, and precisely describing, the test of equality between these solution procedures based on the result is significant at 99.7% significance level with 149 as nonparametric Mann–Whitney statistics and p value of 0.977 by applying Minitab™. As a result, there is no difference between the solution methods’ behaviour in terms of cycle time.From the standpoint of elapsed time, although the MOPSO algorithm has larger amounts on the first four Examples of Table1 thoroughly, on other Examples, this behaviour has changed significantly and the epsilon-constraint method has bigger values.


### 6.2. Different Parts Case

#### 6.2.1. Simulation Results

In order to compare the results from a common circumstance, the simulation model is run for the 2nd example, after relaxing machine failure. This accomplishment enables to equate simulation the simulated results with the analytic method. After running the model more than 10000 hours and replicating the model 25 times, the average S6 cycle time and operating cost damped on 86.719 and 187.5, respectively, at 95% confidence level.


#### 6.2.2. Model Validation Testing

In order to examine validation of the proposed simulation model, the Mann–Whitney hypothesis testing was applied to compare statistically significant difference between the mean of daily production on the simulation model after replicating the model 25 times and the analytic results obtained from equation (20). Hypothesis testing via applying the Minitab statistical package reported high amount of p value (greater than 0.84) which revealed that there is no significant difference to reject the simulation model; it will be concluded that there is no significant difference between results derived from the simulation model and the analytic method.Accordingly, some scenarios through the DOE method are designed to find the best sequence amongst experiments. This method could be conducted to quantitatively evaluate the desired objective functions through estimating fitness function using RSM on the presence of different sources of uncertainty arising from failure/repair process, stochastic setup or process times, and so on.The proposed CS model could easily be customized to any real cases, and there is no restriction on applying a special density function. Experts could straightforwardly extend the model to any other sources of uncertainty stand up from tools, operators, environment/climate, preemption, etc.


#### 6.2.3. RSM

RSM is an assemblage of statistical design and numerical optimization techniques and is used for optimizing processes and product designs [70]. RSM is applied in situations where some input variables influence performance measures of the system. This methodology could make a relationship between the variables. For the proposed simulated robotic cell, first, scenarios are designed and variables are defined.To design the scenarios, DOE is used, and the methodology of fractional factorial DOE is applied. DOE is a method to recognize the main factors in a process and, then, ascertain and fix the problem in a process. In real engineering situations, there are typically multiple factors involved, and it is usually significant to consider them together in the case they interact. A fractional factorial experiment is a disparity of the basic fractional design in a way that only a division of the runs is used [71]. This method reduces the number of experiments one has to perform, by studying multiple factors instantaneously. Factorial design works well when interactions between the variables are strong and important and where every variable contributes significantly.Therefore, input variables are defined includingx1, x2, x3, and x4 correspondingly as the percentage of part’s performed operations on the second machine, the percentage of part’s performed operations on the third machine, the colour of the first allocated Part to the second/third machines, and the sequence type. Table 12 presents input settings on Minitab and their three levels; the low, centre, and high; for designing random scenarios of Example 1 (Different Parts Case).Table 12
Process parameters with their levels.
InputsCode levelsInput variables(−1)0(+1)x1203040x2203040x3BluePurpleRedx4B-P-RR-B-PB-R-PIn Minitab, full factorial design conducted24=16 experiments. By replicating the experiments two times and considering one experiment based on the centre points, total runs are 33. In order to estimate the four mentioned objective functions, the following first-order regression equations are developed considering only the main effects of various process parameters:(22)Y1^=61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+e1,Y2^=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4+e2,Y1SD^=3352.909+655.84375x1+249.09375x2+153.96875−577.59375+e3,Y2SD^=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625+e4.Table13 exhibits the detailed results of the simulated responses as obtained from the 33 experimental trials. Moreover, the regression residuals are calculated via difference between the real and estimated output amounts for the objective functions and reported in Table 13.Table 13
Experimental details, predictions, and residuals of fitness functions for Example1.
RunFactorsPFMsResidualsS. No.M2 (x1)M3 (x2)Col. (x3)Seq. (x4)Y1Y2Y1SDY2SDe1e2e3e41−1−1−1169.989525094.724276310.5−3.3780.2501046.591−0.852−11−1150.458318121.958358210.57−17.188−4695.591367.404−0.533−1−11−170.046925115.143319311.980.132634.34313.4660.6784−11−1−168.270324480.632156710.643.9242291.933−1802.78−0.5815−111175.090624478.783110810.547.5981646.058−1414.53−0.3916111164.0820905.797550111.9110.7822682.4471666.7790.76671−1−1166.235523754.198154511.217.06223269.099−1483.1−0.352811−1−149.15517657.079533211.65−0.99677.755650.5280.2169−111−155.507318123.471467010.37−8.685−4080.4992.278−0.68210111−149.154417657.076827210.8−0.84462.5763282.591−0.46511−1−1−1177.084125125.851213611.73.71631.377419.5910.35012111155.935118262.415274010.462.63739.065−1094.22−0.683131−11−149.164317659.954475611.79−6.556−2211.47264.7780.27514−11−1175.146524497.28111211.687.5011679.731−1102.60.58015−111175.093324479.709191910.627.6011646.984−603.534−0.311161−11156.123618323.137486710.5−2.897−2177.141530.966−0.89317−111−150.462818123.472248711.75−13.730−4080.4−1190.720.6981811−1149.22117680.594583111.87−4.230−527.582304.7160.5581911−1155.962318271.136274011.622.51162.962−786.2840.30720−1−11−170.148925151.202225010.630.234670.402−929.534−0.67221−1−11170.040225113.195496111.62−3.1743.5452936.6540.43922−11−1−168.305724493.579157211.123.9602304.88−1797.78−0.101231−1−1−154.054817651.809372311.58−1.818−2204.44−460.284−0.10324111−149.189717670.027481311.62−0.80875.527−176.4090.35525−1−1−1−176.917125071.426644411.826.849605.8023572.4040.34926−1−11169.974125088.895181611.97−3.241−20.755−208.3470.78927−1−1−1−177.064525119.035252511.736.997653.411−346.5970.2592811−1−149.233317684.92533811.72−0.918105.596656.5280.57629000063.453422761.417144911.261.7701416.93−1903.91−0.047301−11−173.789426450.876249911.8118.0696579.451−1992.220.295311−1−1150.998718314.877152611.98−8.175−2170.22−1502.10.418321−1−1−149.172517662.833439910.9−6.701−2193.42215.716−0.783331−11151.017318321.587121011.22−8.003−2178.69−2126.03−0.173Statistical hypothesis testing for the residuals based on applying normal goodness of fit test revealed that there is no trace to reject the null hypothesis (allp values are greater than 0.76). Consequently, the validation of the estimated fitness functions could be accepted statistically.In the present study, an effort has been made to determine the optimal values of objective functions to get the best possible surface quality within the specific test range. Effective optimization methods including goal programming (GP), multichoice goal programming (MCGP), and augmented epsilon-constraint method (AUGMECON) are developed to solve the optimization problem; for this study, it is called as model (23):(23)MinY1^,MinY2^,MinY1SD,^MinY2SD^,s.t.,−1≤xi≤1,i=1,2,3,4.


#### 6.2.4. Solution Methods

(1) GP Method. The GP method as a MCDM technique is used for solving the optimization problem; model (23). Thus, modelling the GP technique comes in equations (24)–(29) which is the minimization of deviations from the set of goals. It should be noted that di+ and di− are over achievement and under achievement of the goal i one-to-one, in accordance with GP technique.(24)Min∑i=14di++di−,s.t.,(25)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=49.1544,(26)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=17651.809,(27)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=1108,(28)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=10.37,(29)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−≥0.The present work is an effort to effectively perform the Goal programming method for producing the optimal solution in a multicriteria sequencing problem.(2) MCGP Method. Following the idea of MCGP with utility functions “proposed by Chang [72],” first the utility function is depicted in Figure 11; then, model (23) is formulated as comes in the following equations:Figure 11
The linear and left triangular utility functions for the objectives.
(a)(b)(c)(d)(30)Min∑i=14widi++di−+βifi−,s.t.,(31)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=y1,(32)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=y2.(33)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=y3,(34)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=y4,(35)49.1544≤y1≤77.0841,(36)17651.809≤y2≤26450.876,(37)1108≤y3≤8272,(38)10.37≤y4≤11.98,(39)ƛ1≤77.0841−y127.9297,(40)ƛ2≤26450.876−y28799.067,(41)ƛ3≤8272−y37164,(42)ƛ4≤11.98−y41.61,(43)ƛi+fi−=1,i=1,2,3,4,(44)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−,fi−,ƛi≥0.(3) AUGMECON Method. In this study, to produce a representation of the Pareto optimal solutions for the problem, the epsilon-constraint approach is used. Following the idea of AUGMECON “proposed by Mavrotas [73],” firstly, payoff table is created as shown in Table 14. Then, the ranges are calculated and shown in Table 15. Finally, model (23) is formulated in equations (45)–(49).Table 14
Payoff table for the model (23).
Y1Y2Y1SDY2SDY151.64851.72563.33360.393Y217908.90017901.30021658.90020527.400Y1SD4411.8204257.8502775.3203178.320Y2SD11.20411.28811.24611.037Table 15
Calculation of the ranges in payoff table for the model (23).
GoalsRangesY111.685Y23757.600Y1SD1636.500Y2SD0.251(45)Min61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+epsilons23757.600+s31636.500+s40.251,s.t.,(46)Y21000=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x41000,(47)Y1SD1000=3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x41000,(48)Y2SD=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4,(49)−1≤xi≤1,i=1,2,x3,x4=−1or0or1.


## 6.1. Identical Parts Case

### 6.1.1. Epsilon-Constraint Results

(MinF-Min TS6) is acquired by solving the coded epsilon-constraint with solver GAMS-BARON. The codes were run on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU. Figure 4 distinctively plots twenty nondominated solutions of feasible region for objective functions in Test Problems (Table 1). Then, in Table 4, the upper (U) and lower (L) bounds of the generated solutions are presented. It should be noted that the solution method was run 10 times per example, and the small-size and large-size results were achieved by 1800 seconds average computational time. Because of the effect of the group size on the model’s computation time, producing the parts with more than 18 operations will be more than 3600 seconds computation’s time. So, the computational time was adjusted to 3600 seconds in the third test problem, and then the model was run 10 times per example.Figure 4
Epsilon-constraint results for the first twenty iterations in Test Problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 4
Upper and lower bounds of the generated Pareto front in the Test Problems.
Example no.BoundsCostS6 cycle time1L1676.51732U1721.089702L2550.03932U2599.119703L2735.44832U2813.112704L3255.75734U3327.174725L4200.63140U4237.42786L4866.2944U4940.935827L5201.57547U5242.042858L5631.05649U5706.741879L5705.43950U5746.8896910L8147.49266U8186.1358511L11450.8688U11481.4410712L15781.14116U15809135The allocation of operations to the machines based on the proposed model is achievable as well. As an instance, regarding the machines’ generated degradation level, the allocation of operations to the machines for Example 11 (Table1) is demonstrated in Table 5.Table 5
Allocation of operations to the machines, for Example 11, based on the epsilon-constraint method.
Machine no.123Generated degradation level13.880947.08100433.89859Allocated operations no.3, 6, 7, 16, 19, 20, 21, 22, 23, 25, 26, 28, 29, 30, 31, 32, 351, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, 349, 10, 18Example 11 includes 35 operations to manufacture a product. The best allocation of operations to the machines according to Table5 was revealed on the basis of the machines’ generated degradation level. Consequently, the operations numbers 9, 10, and 18 with processing times 10, 10, and 10, respectively, are allocated to the third machine; the operations numbers 1, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, and 34 with their processing times 10, 7, 8, 5, 4, 8, 7, 4, 3, 7, 8, 3, 7, 10, and 3, respectively, are allocated to the second machine; and the rest of the operations are allocated to the first machine. This allocation leads to an instantaneous (MinF–MinTS6). This result determines the maximum output rate subsequently.The reliability condition in comparison with the proposed model was also tested. Considering no breakdown for machines in the cell (i.e., an available robotic cell), the model was run again using commercial software (GAMS). The results for the first five Examples of Test Problems (Table1) are shown in Table 6.Table 6
Coordinates of the optimal point in the designated Examples for the robotic cell in the case of testing the availability condition.
Example no.12345Coordinates of optimal point:F, TS63245, 344845, 445195, 476445, 558045, 66It is clear from the results that in case of considering a reliable cell (i.e., a robotic cell with reliable machines where no failures occur on them), an optimal point for the bicriteria model comes as a specific value, whereas in case of considering machine failure and repair in the robotic cell, a Pareto optimal (nondominated) solution is attainable.


### 6.1.2. MOPSO Results

To express the performance of the proposed model for solving the problem, the proposed MOPSO is examined on the designated Test problems (Table1) by means of MATLAB software on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU.(1) Solution Representation. The MOPSO solution chosen for the problem is a string of decimal numbers. The string as a particle consists of processing times of operations to produce part i (ti). The structure of the string is shown in Figure 5.Figure 5
Forming a solution to solve a problem.All solutions are randomly made in harmony with the above construct. This solution representation is divided into three sections. The first section is related to the allocation of operations to the machines which is illustrated in Figure6. Dimensions of this particle in the form of a matrix are equal to (number of machines∗number of operations). In each row, the maximum value of that row’s processing time is determined. For example, in the first row of the matrix, 0.4 is selected, which means that the cell corresponding to that value takes 1, and the rest of the cells in that row takes the value zero. This means that machine 3 is assigned to process operation 1.Figure 6
Solution representation for allocation of operations to the machines.The second and third sections of solution representation are associated with maintenance decisions for the machines. In these particles, to convert the decimal solution to the zero-one solution, if the value is above 0.5, then the solution is equal to 1, otherwise, it is equal to zero. Figures7 and 8 show maintenance decisions. Based on the example in Figure 7, second and third machines need PM actions. The example of Figure 8 shows that first and third machines require CM operations.Figure 7
Solution representation for PM actions.Figure 8
Solution representation for CM actions.(2) Results. Parameter setting for MOPSO is done by hand in a trial-and-error procedure considering convergence analysis. In this model, w is set as 0.5; r1 and r2 are the random variables between 0 and 1; c1 and c2 are 1 and 2, respectively. The repository size is usually set equal to the swarm (population) size. Here, the algorithm is run for a population size of 50 and a repository size of 50 particles. Number of grids per dimension is 5. The values of other parameters are tabulated in Table 7. It should be noted that grid inflation rate (α) and leader selection pressure (β) are the two predefined parameters. If β>α, there is a cumulative tendency to converge towards PBest for the particle. Higher amount of deletion selection pressure (γ) leads to an increase in probability of member deletion. Deletion selection pressure is a predefined parameter as well.Table 7
MOPSO parameter setting.
Population size = repository sizeαβγNumber of gridsMutation rateMaximum number of iterations500.12250.1100Accordingly, the model is analysed using the MOPSO algorithm, and its results are provided in Table8. In Figure 9, the best results of feasible region for objective functions in Test Problems (Table 1) are represented. Close to the epsilon-constraint, the allocation of operations to the machines for Example 11 (Table 1) is established in Table 9.Table 8
Summary of the objective functions results and elapsed time for solution methods.
Example no.CriteriaMinimum costMinimumS6 cycle timeElapsed time (in seconds)Epsilon-constraintMOPSOEpsilon-constraintMOPSOEpsilon-constraintMOPSO11676.5172281.962323210.8718.93522550.0393451.423323212.3418.58332735.4483871.081323213.9718.46943255.7574556.381343414.0418.95254200.6315909.7424041598.9516.70264866.296695.4024444620.6718.61475201.5757181.277474790017.16885631.0567802.1394949180019.14695705.4398085.5075050360028.957108147.49211245.8196667360026.8491111450.8616030.4278888360028.9431215781.1422353.457116116360032.526Figure 9
MOPSO results of feasible region for objective functions in Test problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 9
Allocation of operations to the machines, for Example 11, based on the MOPSO algorithm.
Machine no.123Allocated operations no.6, 10, 13, 16, 17, 18, 20, 21, 24, 25, 28, 29, 353, 7, 8, 9, 11, 22, 23, 26, 30, 321, 2, 4, 5, 12, 14, 15, 19, 27, 31, 33, 34In Example 11, the best allocation of operations to the machines according to Table9 and the reason of the machines’ generated degradation level was disclosed. With the declared allocation, minimum cost and cycle time simultaneously will be practicable, and also maximum output rate is reasonable.A comparison between the solution methods for Example 11, based on Table8 shows that, however, minimum cycle time is 88 based on both solution methods, and minimum cost was achieved regarding epsilon-constraint method with the amount of 11450.86 $/min. In other Examples of this Test Problem also, minimum cost is realisable considering epsilon-constraint approach. However, MOPSO’s elapsed time is less than that of epsilon-constraint according to this Example.(3) Comparison Metrics. Diversity and spacing metrics are provided for assessing the relative performance of MOPSO. The metrics definition is summarized in Table 10, where xi–yi is the Euclidean distance between two nondominated solutions xi and yi [69]. Comparison metrics were calculated for Test problems (Table 1) over 10 runs of the MOPSO algorithm. Table 11 presents the metrics obtained for small-, large-, and big-sized instances.Table 10
Comparison metrics.
Spacing metric (SM)d⇀=∑i=1ndi/nDiversity metric (DM)DM=∑i=1nmaxxi−yiSM=∑i=1nd⇀−di2/n−1Table 11
Computational results for comparison metrics in Test Problems.
Test ProblemExample no.DMSMSmall size1177.21502157.99903213.0850Large size4273.80505283.44906320.88107286.04308320.690Big size9326.375010386.960011475.863012568.5890


### 6.1.3. Statistical Analysis

Objective functions and elapsed time required by each method are defined as criteria in Table8 so that the outcome of solution methods will be compared. To evaluate the performance of solution methods, epsilon-constraint method’s minimum cost was measured against that of MOPSO algorithm. F (i.e., cost), in exact solution, is consistently less than metaheuristic solution, regarding Figure 10. In the considered Test Problems (Table 1), difference between the objective function values is at least 605 and at most 6572 currencies.Figure 10
Comparison ofF (cost) based on the solution procedures in Test Problems.A statistical significant difference between the solution methods is not observed forTS6, and precisely describing, the test of equality between these solution procedures based on the result is significant at 99.7% significance level with 149 as nonparametric Mann–Whitney statistics and p value of 0.977 by applying Minitab™. As a result, there is no difference between the solution methods’ behaviour in terms of cycle time.From the standpoint of elapsed time, although the MOPSO algorithm has larger amounts on the first four Examples of Table1 thoroughly, on other Examples, this behaviour has changed significantly and the epsilon-constraint method has bigger values.


## 6.1.1. Epsilon-Constraint Results

(MinF-Min TS6) is acquired by solving the coded epsilon-constraint with solver GAMS-BARON. The codes were run on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU. Figure 4 distinctively plots twenty nondominated solutions of feasible region for objective functions in Test Problems (Table 1). Then, in Table 4, the upper (U) and lower (L) bounds of the generated solutions are presented. It should be noted that the solution method was run 10 times per example, and the small-size and large-size results were achieved by 1800 seconds average computational time. Because of the effect of the group size on the model’s computation time, producing the parts with more than 18 operations will be more than 3600 seconds computation’s time. So, the computational time was adjusted to 3600 seconds in the third test problem, and then the model was run 10 times per example.Figure 4
Epsilon-constraint results for the first twenty iterations in Test Problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 4
Upper and lower bounds of the generated Pareto front in the Test Problems.
Example no.BoundsCostS6 cycle time1L1676.51732U1721.089702L2550.03932U2599.119703L2735.44832U2813.112704L3255.75734U3327.174725L4200.63140U4237.42786L4866.2944U4940.935827L5201.57547U5242.042858L5631.05649U5706.741879L5705.43950U5746.8896910L8147.49266U8186.1358511L11450.8688U11481.4410712L15781.14116U15809135The allocation of operations to the machines based on the proposed model is achievable as well. As an instance, regarding the machines’ generated degradation level, the allocation of operations to the machines for Example 11 (Table1) is demonstrated in Table 5.Table 5
Allocation of operations to the machines, for Example 11, based on the epsilon-constraint method.
Machine no.123Generated degradation level13.880947.08100433.89859Allocated operations no.3, 6, 7, 16, 19, 20, 21, 22, 23, 25, 26, 28, 29, 30, 31, 32, 351, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, 349, 10, 18Example 11 includes 35 operations to manufacture a product. The best allocation of operations to the machines according to Table5 was revealed on the basis of the machines’ generated degradation level. Consequently, the operations numbers 9, 10, and 18 with processing times 10, 10, and 10, respectively, are allocated to the third machine; the operations numbers 1, 2, 4, 5, 8, 11, 12, 13, 14, 15, 17, 24, 27, 33, and 34 with their processing times 10, 7, 8, 5, 4, 8, 7, 4, 3, 7, 8, 3, 7, 10, and 3, respectively, are allocated to the second machine; and the rest of the operations are allocated to the first machine. This allocation leads to an instantaneous (MinF–MinTS6). This result determines the maximum output rate subsequently.The reliability condition in comparison with the proposed model was also tested. Considering no breakdown for machines in the cell (i.e., an available robotic cell), the model was run again using commercial software (GAMS). The results for the first five Examples of Test Problems (Table1) are shown in Table 6.Table 6
Coordinates of the optimal point in the designated Examples for the robotic cell in the case of testing the availability condition.
Example no.12345Coordinates of optimal point:F, TS63245, 344845, 445195, 476445, 558045, 66It is clear from the results that in case of considering a reliable cell (i.e., a robotic cell with reliable machines where no failures occur on them), an optimal point for the bicriteria model comes as a specific value, whereas in case of considering machine failure and repair in the robotic cell, a Pareto optimal (nondominated) solution is attainable.


## 6.1.2. MOPSO Results

To express the performance of the proposed model for solving the problem, the proposed MOPSO is examined on the designated Test problems (Table1) by means of MATLAB software on a portable PC with MS-Windows Vista, 3.0 GB of RAM, and 2.0 GHz Core 2 Duo CPU.(1) Solution Representation. The MOPSO solution chosen for the problem is a string of decimal numbers. The string as a particle consists of processing times of operations to produce part i (ti). The structure of the string is shown in Figure 5.Figure 5
Forming a solution to solve a problem.All solutions are randomly made in harmony with the above construct. This solution representation is divided into three sections. The first section is related to the allocation of operations to the machines which is illustrated in Figure6. Dimensions of this particle in the form of a matrix are equal to (number of machines∗number of operations). In each row, the maximum value of that row’s processing time is determined. For example, in the first row of the matrix, 0.4 is selected, which means that the cell corresponding to that value takes 1, and the rest of the cells in that row takes the value zero. This means that machine 3 is assigned to process operation 1.Figure 6
Solution representation for allocation of operations to the machines.The second and third sections of solution representation are associated with maintenance decisions for the machines. In these particles, to convert the decimal solution to the zero-one solution, if the value is above 0.5, then the solution is equal to 1, otherwise, it is equal to zero. Figures7 and 8 show maintenance decisions. Based on the example in Figure 7, second and third machines need PM actions. The example of Figure 8 shows that first and third machines require CM operations.Figure 7
Solution representation for PM actions.Figure 8
Solution representation for CM actions.(2) Results. Parameter setting for MOPSO is done by hand in a trial-and-error procedure considering convergence analysis. In this model, w is set as 0.5; r1 and r2 are the random variables between 0 and 1; c1 and c2 are 1 and 2, respectively. The repository size is usually set equal to the swarm (population) size. Here, the algorithm is run for a population size of 50 and a repository size of 50 particles. Number of grids per dimension is 5. The values of other parameters are tabulated in Table 7. It should be noted that grid inflation rate (α) and leader selection pressure (β) are the two predefined parameters. If β>α, there is a cumulative tendency to converge towards PBest for the particle. Higher amount of deletion selection pressure (γ) leads to an increase in probability of member deletion. Deletion selection pressure is a predefined parameter as well.Table 7
MOPSO parameter setting.
Population size = repository sizeαβγNumber of gridsMutation rateMaximum number of iterations500.12250.1100Accordingly, the model is analysed using the MOPSO algorithm, and its results are provided in Table8. In Figure 9, the best results of feasible region for objective functions in Test Problems (Table 1) are represented. Close to the epsilon-constraint, the allocation of operations to the machines for Example 11 (Table 1) is established in Table 9.Table 8
Summary of the objective functions results and elapsed time for solution methods.
Example no.CriteriaMinimum costMinimumS6 cycle timeElapsed time (in seconds)Epsilon-constraintMOPSOEpsilon-constraintMOPSOEpsilon-constraintMOPSO11676.5172281.962323210.8718.93522550.0393451.423323212.3418.58332735.4483871.081323213.9718.46943255.7574556.381343414.0418.95254200.6315909.7424041598.9516.70264866.296695.4024444620.6718.61475201.5757181.277474790017.16885631.0567802.1394949180019.14695705.4398085.5075050360028.957108147.49211245.8196667360026.8491111450.8616030.4278888360028.9431215781.1422353.457116116360032.526Figure 9
MOPSO results of feasible region for objective functions in Test problems. (a) Example 1. (b) Example 2. (c) Example 3. (d) Example 4. (e) Example 5. (f) Example 6. (g) Example 7. (h) Example 8. (i) Example 9. (j) Example 10. (k) Example 11. (l) Example 12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Table 9
Allocation of operations to the machines, for Example 11, based on the MOPSO algorithm.
Machine no.123Allocated operations no.6, 10, 13, 16, 17, 18, 20, 21, 24, 25, 28, 29, 353, 7, 8, 9, 11, 22, 23, 26, 30, 321, 2, 4, 5, 12, 14, 15, 19, 27, 31, 33, 34In Example 11, the best allocation of operations to the machines according to Table9 and the reason of the machines’ generated degradation level was disclosed. With the declared allocation, minimum cost and cycle time simultaneously will be practicable, and also maximum output rate is reasonable.A comparison between the solution methods for Example 11, based on Table8 shows that, however, minimum cycle time is 88 based on both solution methods, and minimum cost was achieved regarding epsilon-constraint method with the amount of 11450.86 $/min. In other Examples of this Test Problem also, minimum cost is realisable considering epsilon-constraint approach. However, MOPSO’s elapsed time is less than that of epsilon-constraint according to this Example.(3) Comparison Metrics. Diversity and spacing metrics are provided for assessing the relative performance of MOPSO. The metrics definition is summarized in Table 10, where xi–yi is the Euclidean distance between two nondominated solutions xi and yi [69]. Comparison metrics were calculated for Test problems (Table 1) over 10 runs of the MOPSO algorithm. Table 11 presents the metrics obtained for small-, large-, and big-sized instances.Table 10
Comparison metrics.
Spacing metric (SM)d⇀=∑i=1ndi/nDiversity metric (DM)DM=∑i=1nmaxxi−yiSM=∑i=1nd⇀−di2/n−1Table 11
Computational results for comparison metrics in Test Problems.
Test ProblemExample no.DMSMSmall size1177.21502157.99903213.0850Large size4273.80505283.44906320.88107286.04308320.690Big size9326.375010386.960011475.863012568.5890


## 6.1.3. Statistical Analysis

Objective functions and elapsed time required by each method are defined as criteria in Table8 so that the outcome of solution methods will be compared. To evaluate the performance of solution methods, epsilon-constraint method’s minimum cost was measured against that of MOPSO algorithm. F (i.e., cost), in exact solution, is consistently less than metaheuristic solution, regarding Figure 10. In the considered Test Problems (Table 1), difference between the objective function values is at least 605 and at most 6572 currencies.Figure 10
Comparison ofF (cost) based on the solution procedures in Test Problems.A statistical significant difference between the solution methods is not observed forTS6, and precisely describing, the test of equality between these solution procedures based on the result is significant at 99.7% significance level with 149 as nonparametric Mann–Whitney statistics and p value of 0.977 by applying Minitab™. As a result, there is no difference between the solution methods’ behaviour in terms of cycle time.From the standpoint of elapsed time, although the MOPSO algorithm has larger amounts on the first four Examples of Table1 thoroughly, on other Examples, this behaviour has changed significantly and the epsilon-constraint method has bigger values.


## 6.2. Different Parts Case

### 6.2.1. Simulation Results

In order to compare the results from a common circumstance, the simulation model is run for the 2nd example, after relaxing machine failure. This accomplishment enables to equate simulation the simulated results with the analytic method. After running the model more than 10000 hours and replicating the model 25 times, the average S6 cycle time and operating cost damped on 86.719 and 187.5, respectively, at 95% confidence level.


### 6.2.2. Model Validation Testing

In order to examine validation of the proposed simulation model, the Mann–Whitney hypothesis testing was applied to compare statistically significant difference between the mean of daily production on the simulation model after replicating the model 25 times and the analytic results obtained from equation (20). Hypothesis testing via applying the Minitab statistical package reported high amount of p value (greater than 0.84) which revealed that there is no significant difference to reject the simulation model; it will be concluded that there is no significant difference between results derived from the simulation model and the analytic method.Accordingly, some scenarios through the DOE method are designed to find the best sequence amongst experiments. This method could be conducted to quantitatively evaluate the desired objective functions through estimating fitness function using RSM on the presence of different sources of uncertainty arising from failure/repair process, stochastic setup or process times, and so on.The proposed CS model could easily be customized to any real cases, and there is no restriction on applying a special density function. Experts could straightforwardly extend the model to any other sources of uncertainty stand up from tools, operators, environment/climate, preemption, etc.


### 6.2.3. RSM

RSM is an assemblage of statistical design and numerical optimization techniques and is used for optimizing processes and product designs [70]. RSM is applied in situations where some input variables influence performance measures of the system. This methodology could make a relationship between the variables. For the proposed simulated robotic cell, first, scenarios are designed and variables are defined.To design the scenarios, DOE is used, and the methodology of fractional factorial DOE is applied. DOE is a method to recognize the main factors in a process and, then, ascertain and fix the problem in a process. In real engineering situations, there are typically multiple factors involved, and it is usually significant to consider them together in the case they interact. A fractional factorial experiment is a disparity of the basic fractional design in a way that only a division of the runs is used [71]. This method reduces the number of experiments one has to perform, by studying multiple factors instantaneously. Factorial design works well when interactions between the variables are strong and important and where every variable contributes significantly.Therefore, input variables are defined includingx1, x2, x3, and x4 correspondingly as the percentage of part’s performed operations on the second machine, the percentage of part’s performed operations on the third machine, the colour of the first allocated Part to the second/third machines, and the sequence type. Table 12 presents input settings on Minitab and their three levels; the low, centre, and high; for designing random scenarios of Example 1 (Different Parts Case).Table 12
Process parameters with their levels.
InputsCode levelsInput variables(−1)0(+1)x1203040x2203040x3BluePurpleRedx4B-P-RR-B-PB-R-PIn Minitab, full factorial design conducted24=16 experiments. By replicating the experiments two times and considering one experiment based on the centre points, total runs are 33. In order to estimate the four mentioned objective functions, the following first-order regression equations are developed considering only the main effects of various process parameters:(22)Y1^=61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+e1,Y2^=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4+e2,Y1SD^=3352.909+655.84375x1+249.09375x2+153.96875−577.59375+e3,Y2SD^=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625+e4.Table13 exhibits the detailed results of the simulated responses as obtained from the 33 experimental trials. Moreover, the regression residuals are calculated via difference between the real and estimated output amounts for the objective functions and reported in Table 13.Table 13
Experimental details, predictions, and residuals of fitness functions for Example1.
RunFactorsPFMsResidualsS. No.M2 (x1)M3 (x2)Col. (x3)Seq. (x4)Y1Y2Y1SDY2SDe1e2e3e41−1−1−1169.989525094.724276310.5−3.3780.2501046.591−0.852−11−1150.458318121.958358210.57−17.188−4695.591367.404−0.533−1−11−170.046925115.143319311.980.132634.34313.4660.6784−11−1−168.270324480.632156710.643.9242291.933−1802.78−0.5815−111175.090624478.783110810.547.5981646.058−1414.53−0.3916111164.0820905.797550111.9110.7822682.4471666.7790.76671−1−1166.235523754.198154511.217.06223269.099−1483.1−0.352811−1−149.15517657.079533211.65−0.99677.755650.5280.2169−111−155.507318123.471467010.37−8.685−4080.4992.278−0.68210111−149.154417657.076827210.8−0.84462.5763282.591−0.46511−1−1−1177.084125125.851213611.73.71631.377419.5910.35012111155.935118262.415274010.462.63739.065−1094.22−0.683131−11−149.164317659.954475611.79−6.556−2211.47264.7780.27514−11−1175.146524497.28111211.687.5011679.731−1102.60.58015−111175.093324479.709191910.627.6011646.984−603.534−0.311161−11156.123618323.137486710.5−2.897−2177.141530.966−0.89317−111−150.462818123.472248711.75−13.730−4080.4−1190.720.6981811−1149.22117680.594583111.87−4.230−527.582304.7160.5581911−1155.962318271.136274011.622.51162.962−786.2840.30720−1−11−170.148925151.202225010.630.234670.402−929.534−0.67221−1−11170.040225113.195496111.62−3.1743.5452936.6540.43922−11−1−168.305724493.579157211.123.9602304.88−1797.78−0.101231−1−1−154.054817651.809372311.58−1.818−2204.44−460.284−0.10324111−149.189717670.027481311.62−0.80875.527−176.4090.35525−1−1−1−176.917125071.426644411.826.849605.8023572.4040.34926−1−11169.974125088.895181611.97−3.241−20.755−208.3470.78927−1−1−1−177.064525119.035252511.736.997653.411−346.5970.2592811−1−149.233317684.92533811.72−0.918105.596656.5280.57629000063.453422761.417144911.261.7701416.93−1903.91−0.047301−11−173.789426450.876249911.8118.0696579.451−1992.220.295311−1−1150.998718314.877152611.98−8.175−2170.22−1502.10.418321−1−1−149.172517662.833439910.9−6.701−2193.42215.716−0.783331−11151.017318321.587121011.22−8.003−2178.69−2126.03−0.173Statistical hypothesis testing for the residuals based on applying normal goodness of fit test revealed that there is no trace to reject the null hypothesis (allp values are greater than 0.76). Consequently, the validation of the estimated fitness functions could be accepted statistically.In the present study, an effort has been made to determine the optimal values of objective functions to get the best possible surface quality within the specific test range. Effective optimization methods including goal programming (GP), multichoice goal programming (MCGP), and augmented epsilon-constraint method (AUGMECON) are developed to solve the optimization problem; for this study, it is called as model (23):(23)MinY1^,MinY2^,MinY1SD,^MinY2SD^,s.t.,−1≤xi≤1,i=1,2,3,4.


### 6.2.4. Solution Methods

(1) GP Method. The GP method as a MCDM technique is used for solving the optimization problem; model (23). Thus, modelling the GP technique comes in equations (24)–(29) which is the minimization of deviations from the set of goals. It should be noted that di+ and di− are over achievement and under achievement of the goal i one-to-one, in accordance with GP technique.(24)Min∑i=14di++di−,s.t.,(25)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=49.1544,(26)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=17651.809,(27)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=1108,(28)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=10.37,(29)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−≥0.The present work is an effort to effectively perform the Goal programming method for producing the optimal solution in a multicriteria sequencing problem.(2) MCGP Method. Following the idea of MCGP with utility functions “proposed by Chang [72],” first the utility function is depicted in Figure 11; then, model (23) is formulated as comes in the following equations:Figure 11
The linear and left triangular utility functions for the objectives.
(a)(b)(c)(d)(30)Min∑i=14widi++di−+βifi−,s.t.,(31)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=y1,(32)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=y2.(33)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=y3,(34)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=y4,(35)49.1544≤y1≤77.0841,(36)17651.809≤y2≤26450.876,(37)1108≤y3≤8272,(38)10.37≤y4≤11.98,(39)ƛ1≤77.0841−y127.9297,(40)ƛ2≤26450.876−y28799.067,(41)ƛ3≤8272−y37164,(42)ƛ4≤11.98−y41.61,(43)ƛi+fi−=1,i=1,2,3,4,(44)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−,fi−,ƛi≥0.(3) AUGMECON Method. In this study, to produce a representation of the Pareto optimal solutions for the problem, the epsilon-constraint approach is used. Following the idea of AUGMECON “proposed by Mavrotas [73],” firstly, payoff table is created as shown in Table 14. Then, the ranges are calculated and shown in Table 15. Finally, model (23) is formulated in equations (45)–(49).Table 14
Payoff table for the model (23).
Y1Y2Y1SDY2SDY151.64851.72563.33360.393Y217908.90017901.30021658.90020527.400Y1SD4411.8204257.8502775.3203178.320Y2SD11.20411.28811.24611.037Table 15
Calculation of the ranges in payoff table for the model (23).
GoalsRangesY111.685Y23757.600Y1SD1636.500Y2SD0.251(45)Min61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+epsilons23757.600+s31636.500+s40.251,s.t.,(46)Y21000=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x41000,(47)Y1SD1000=3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x41000,(48)Y2SD=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4,(49)−1≤xi≤1,i=1,2,x3,x4=−1or0or1.


## 6.2.1. Simulation Results

In order to compare the results from a common circumstance, the simulation model is run for the 2nd example, after relaxing machine failure. This accomplishment enables to equate simulation the simulated results with the analytic method. After running the model more than 10000 hours and replicating the model 25 times, the average S6 cycle time and operating cost damped on 86.719 and 187.5, respectively, at 95% confidence level.


## 6.2.2. Model Validation Testing

In order to examine validation of the proposed simulation model, the Mann–Whitney hypothesis testing was applied to compare statistically significant difference between the mean of daily production on the simulation model after replicating the model 25 times and the analytic results obtained from equation (20). Hypothesis testing via applying the Minitab statistical package reported high amount of p value (greater than 0.84) which revealed that there is no significant difference to reject the simulation model; it will be concluded that there is no significant difference between results derived from the simulation model and the analytic method.Accordingly, some scenarios through the DOE method are designed to find the best sequence amongst experiments. This method could be conducted to quantitatively evaluate the desired objective functions through estimating fitness function using RSM on the presence of different sources of uncertainty arising from failure/repair process, stochastic setup or process times, and so on.The proposed CS model could easily be customized to any real cases, and there is no restriction on applying a special density function. Experts could straightforwardly extend the model to any other sources of uncertainty stand up from tools, operators, environment/climate, preemption, etc.


## 6.2.3. RSM

RSM is an assemblage of statistical design and numerical optimization techniques and is used for optimizing processes and product designs [70]. RSM is applied in situations where some input variables influence performance measures of the system. This methodology could make a relationship between the variables. For the proposed simulated robotic cell, first, scenarios are designed and variables are defined.To design the scenarios, DOE is used, and the methodology of fractional factorial DOE is applied. DOE is a method to recognize the main factors in a process and, then, ascertain and fix the problem in a process. In real engineering situations, there are typically multiple factors involved, and it is usually significant to consider them together in the case they interact. A fractional factorial experiment is a disparity of the basic fractional design in a way that only a division of the runs is used [71]. This method reduces the number of experiments one has to perform, by studying multiple factors instantaneously. Factorial design works well when interactions between the variables are strong and important and where every variable contributes significantly.Therefore, input variables are defined includingx1, x2, x3, and x4 correspondingly as the percentage of part’s performed operations on the second machine, the percentage of part’s performed operations on the third machine, the colour of the first allocated Part to the second/third machines, and the sequence type. Table 12 presents input settings on Minitab and their three levels; the low, centre, and high; for designing random scenarios of Example 1 (Different Parts Case).Table 12
Process parameters with their levels.
InputsCode levelsInput variables(−1)0(+1)x1203040x2203040x3BluePurpleRedx4B-P-RR-B-PB-R-PIn Minitab, full factorial design conducted24=16 experiments. By replicating the experiments two times and considering one experiment based on the centre points, total runs are 33. In order to estimate the four mentioned objective functions, the following first-order regression equations are developed considering only the main effects of various process parameters:(22)Y1^=61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+e1,Y2^=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4+e2,Y1SD^=3352.909+655.84375x1+249.09375x2+153.96875−577.59375+e3,Y2SD^=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625+e4.Table13 exhibits the detailed results of the simulated responses as obtained from the 33 experimental trials. Moreover, the regression residuals are calculated via difference between the real and estimated output amounts for the objective functions and reported in Table 13.Table 13
Experimental details, predictions, and residuals of fitness functions for Example1.
RunFactorsPFMsResidualsS. No.M2 (x1)M3 (x2)Col. (x3)Seq. (x4)Y1Y2Y1SDY2SDe1e2e3e41−1−1−1169.989525094.724276310.5−3.3780.2501046.591−0.852−11−1150.458318121.958358210.57−17.188−4695.591367.404−0.533−1−11−170.046925115.143319311.980.132634.34313.4660.6784−11−1−168.270324480.632156710.643.9242291.933−1802.78−0.5815−111175.090624478.783110810.547.5981646.058−1414.53−0.3916111164.0820905.797550111.9110.7822682.4471666.7790.76671−1−1166.235523754.198154511.217.06223269.099−1483.1−0.352811−1−149.15517657.079533211.65−0.99677.755650.5280.2169−111−155.507318123.471467010.37−8.685−4080.4992.278−0.68210111−149.154417657.076827210.8−0.84462.5763282.591−0.46511−1−1−1177.084125125.851213611.73.71631.377419.5910.35012111155.935118262.415274010.462.63739.065−1094.22−0.683131−11−149.164317659.954475611.79−6.556−2211.47264.7780.27514−11−1175.146524497.28111211.687.5011679.731−1102.60.58015−111175.093324479.709191910.627.6011646.984−603.534−0.311161−11156.123618323.137486710.5−2.897−2177.141530.966−0.89317−111−150.462818123.472248711.75−13.730−4080.4−1190.720.6981811−1149.22117680.594583111.87−4.230−527.582304.7160.5581911−1155.962318271.136274011.622.51162.962−786.2840.30720−1−11−170.148925151.202225010.630.234670.402−929.534−0.67221−1−11170.040225113.195496111.62−3.1743.5452936.6540.43922−11−1−168.305724493.579157211.123.9602304.88−1797.78−0.101231−1−1−154.054817651.809372311.58−1.818−2204.44−460.284−0.10324111−149.189717670.027481311.62−0.80875.527−176.4090.35525−1−1−1−176.917125071.426644411.826.849605.8023572.4040.34926−1−11169.974125088.895181611.97−3.241−20.755−208.3470.78927−1−1−1−177.064525119.035252511.736.997653.411−346.5970.2592811−1−149.233317684.92533811.72−0.918105.596656.5280.57629000063.453422761.417144911.261.7701416.93−1903.91−0.047301−11−173.789426450.876249911.8118.0696579.451−1992.220.295311−1−1150.998718314.877152611.98−8.175−2170.22−1502.10.418321−1−1−149.172517662.833439910.9−6.701−2193.42215.716−0.783331−11151.017318321.587121011.22−8.003−2178.69−2126.03−0.173Statistical hypothesis testing for the residuals based on applying normal goodness of fit test revealed that there is no trace to reject the null hypothesis (allp values are greater than 0.76). Consequently, the validation of the estimated fitness functions could be accepted statistically.In the present study, an effort has been made to determine the optimal values of objective functions to get the best possible surface quality within the specific test range. Effective optimization methods including goal programming (GP), multichoice goal programming (MCGP), and augmented epsilon-constraint method (AUGMECON) are developed to solve the optimization problem; for this study, it is called as model (23):(23)MinY1^,MinY2^,MinY1SD,^MinY2SD^,s.t.,−1≤xi≤1,i=1,2,3,4.


## 6.2.4. Solution Methods

(1) GP Method. The GP method as a MCDM technique is used for solving the optimization problem; model (23). Thus, modelling the GP technique comes in equations (24)–(29) which is the minimization of deviations from the set of goals. It should be noted that di+ and di− are over achievement and under achievement of the goal i one-to-one, in accordance with GP technique.(24)Min∑i=14di++di−,s.t.,(25)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=49.1544,(26)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=17651.809,(27)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=1108,(28)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=10.37,(29)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−≥0.The present work is an effort to effectively perform the Goal programming method for producing the optimal solution in a multicriteria sequencing problem.(2) MCGP Method. Following the idea of MCGP with utility functions “proposed by Chang [72],” first the utility function is depicted in Figure 11; then, model (23) is formulated as comes in the following equations:Figure 11
The linear and left triangular utility functions for the objectives.
(a)(b)(c)(d)(30)Min∑i=14widi++di−+βifi−,s.t.,(31)61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4−d1++d1−=y1,(32)21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x4−d2++d2−=y2.(33)3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x4−d3++d3−=y3,(34)11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4−d4++d4−=y4,(35)49.1544≤y1≤77.0841,(36)17651.809≤y2≤26450.876,(37)1108≤y3≤8272,(38)10.37≤y4≤11.98,(39)ƛ1≤77.0841−y127.9297,(40)ƛ2≤26450.876−y28799.067,(41)ƛ3≤8272−y37164,(42)ƛ4≤11.98−y41.61,(43)ƛi+fi−=1,i=1,2,3,4,(44)−1≤xi≤1,i=1,2,x3,x4=−1or0or1,di+,di−,fi−,ƛi≥0.(3) AUGMECON Method. In this study, to produce a representation of the Pareto optimal solutions for the problem, the epsilon-constraint approach is used. Following the idea of AUGMECON “proposed by Mavrotas [73],” firstly, payoff table is created as shown in Table 14. Then, the ranges are calculated and shown in Table 15. Finally, model (23) is formulated in equations (45)–(49).Table 14
Payoff table for the model (23).
Y1Y2Y1SDY2SDY151.64851.72563.33360.393Y217908.90017901.30021658.90020527.400Y1SD4411.8204257.8502775.3203178.320Y2SD11.20411.28811.24611.037Table 15
Calculation of the ranges in payoff table for the model (23).
GoalsRangesY111.685Y23757.600Y1SD1636.500Y2SD0.251(45)Min61.68304−7.0973x1−2.861x2−0.0766x3+1.65x4+epsilons23757.600+s31636.500+s40.251,s.t.,(46)Y21000=21344.487−2304.6875x1−1138.4625x2+7.58775x3+314.425x41000,(47)Y1SD1000=3352.909+655.84375x1+249.09375x2+153.96875x3−577.59375x41000,(48)Y2SD=11.30727+0.10625x1−0.125x2−0.084375x3−0.060625x4,(49)−1≤xi≤1,i=1,2,x3,x4=−1or0or1.


## 7. Discussion

The problem is solved based on the solution methods via LINGO software on a PC with MS-Windows 10 Pro, 4.0 GB of RAM, and 4 COMPUTE CORES 2C+2G to achieve the optimal solutions. The results based on GP and MCGP are presented in Table16.Table 16
The results of GP and MCGP optimization methods.
GPx1∗ = 1x2∗ = 1x3∗ = 0x4∗ = 1Goal 1 = 53.375Goal 2 = 18215.756Goal 3 = 3680.253Goal 4 = 11.228MCGPx1∗ = 0x2∗ = 0x3∗ = 1x4∗ = 1Goal 1 = 63.256Goal 2 = 21666.499Goal 3 = 2929.284Goal 4 = 11.217The subsequent results are obtained after solving model (23) with AUGMECON. Regarding this solution method, Pareto optimal solutions are gained. The Pareto optimal solutions are shown in Table 17. In the next step, the ED simulation tool is being attuned to the above results and the above optimal scenarios are being run. Then, according to the defined PFMs (Proposed Simulation Model), the results are shown in Table 18.Table 17
Pareto optimal solutions based on AUGMECON.
Solution no.e2e3e4Y1Y2Y1SDY2SD117901.3003320.8211.20451.64817908.924411.81511.204219153.833320.8211.20451.64817908.924411.81511.204320406.3663866.3211.28851.64817908.924411.81511.204421658.9002775.32011.03751.64817908.924411.81511.204519153.832775.32011.12051.64817908.924411.81511.204619153.832775.32011.20451.64817908.924411.81511.204720406.3662775.32011.12051.64817908.924411.81511.204x1∗ = 1x2∗ = 1x3∗ = 1x4∗ = 0Table 18
Simulation results for optimal scenarios in optimization solution methods.
Methods/PFMsY1Y2Y1SDY2SDGP55.86118238.34269810.72MCGP54.45618252.98473710.57AUGMECON50.81818250.114118311.22There is an acceptable closeness between the solutions in Tables16–18, by comparing the results of CS and RSM. This indicates the accuracy of prediction made by RSM.Consistent with the above results, optimal scenarios have the following characteristics.Supposing that the percentage of part’s performed operations on second and third machines as well is at least 30%:If the sequence type is R-B-P, place the first part in red on second and third machines.If the sequence type is B-R-P, the initial part can also be in purple on the second and third machines.A set of these conditions can provide robust optimization for the examined robotic cell to meet the defined goals.


## 8. Conclusions and Future Research Directions

This study purely tends to investigate the problem of three-machine robotic cell from two view points: when the system produces similar parts which leads to a scheduling-allocation problem and when the robotic cell produces dissimilar parts which leads to a sequencing problem.For the first problem, a mathematical model was proposed. The proposed model was focused basically on minimizing jointly operational cost and cycle time in the robotic manufacturing cell under breakdowns including failures and repairs. In this case, the operating conditions decisions in robotic cells link up with maintenance decisions will improve time and operational costs simultaneously. Consequently, this enhancement increases the production rate. This work is an effort to implement the epsilon-constraint approach for producing the Pareto optimal (nondominated) solutions in a mixed integer nonlinear programming problem in comparison with MOPSO algorithm. Numerical experiments were carried out on several designated instances, and the results were discussed.The solution method for the second problem purely has a tendency to demonstrate robustness by simulation-based optimization to determine the sequencing in a defective three-machine robotic cell underS6 cycle. By finding the best scenarios for the entrance of parts into the robotic cell, a robust simultaneous minimization of the time and operating costs is conceivable.The contribution can be extended to the robotic cell considering robot failures or to the dual-gripper robot instead of single-gripper one. Eliminating flexibility in the robotic cell is also considerable for future research. On the contrary, developing new models with other probability distributions of the time to failure and time to repair for the machines besides the exponential probability distribution is an interest. Moreover, the present study can allow modelling different scenarios to be examined in forthcoming simulations with less expense, time, and resources than experimentation with the real-world system. Thus, enhancement of the manufacturing system earlier or through manufacture will be possible.



---

*Source: 4597827-2020-10-31.xml*