# A New Multiobjective Time-Cost Trade-Off for Scheduling Maintenance Problem in a Series-Parallel System

Authors: Leyla Sadat Tavassoli; Reza Massah; Arsalan Montazeri; Mirpouya Mirmozaffari; Guang-Jun Jiang; Hong-Xia Chen
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a modified model of Nondominated Sorting Genetic Algorithm 2 (NSGA-II), which is one of the Multiobjective Evolutionary Algorithms, is proposed. This algorithm is a new model designed to make a trade-off between minimizing the cost of preventive maintenance (PM) and minimizing the time taken to perform this maintenance for a series-parallel system. In this model, the limitations of labor and equipment of the maintenance team and the effects of maintenance issues on manufacturing problems are also considered. In the mathematical model, finding the appropriate objective functions for the maintenance scheduling problem requires all maintenance costs and failure rates to be integrated. Additionally, the effects of production interruption during preventive maintenance are added to objective functions. Furthermore, to make a better performance compared with a regular NSGA-II algorithm, we proposed a modified algorithm with a repository to keep more unacceptable solutions. These solutions can be modified and changed with the proposed mutation algorithm to acceptable solutions. In this algorithm, modified operators, such as simulated binary crossover and polynomial mutation, will improve the algorithm to generate convergence and uniformly distributed solutions with more diverse solutions. Finally, by comparing the experimental solutions with the solutions of two Strength Pareto Evolutionary Algorithm 2 (SPEA2) and regular NSGA-II, MNSGA-II generates more efficient and uniform solutions than the other two algorithms.


---

## Body

## 1. Introduction

In today's industrial world, it is crucial for manufacturing companies to keep the production rates of machines constant. Since cost reduction and profit increase are the main goals of all manufacturing companies, the breakdown of production machines can cause a decrease in production or stop production in some cases, which will reduce the profits of companies eventually. In this situation, the need for a preventive maintenance (PM) system to keep the machines running is essential [1, 2].Calculating the optimal time for PM actions will prevent not only the unexpected breakdown of machines but also save cost as too many maintenance operations could potentially increase the cost of production [3]. It is more complex when two preventive maintenance programs coincide. For example, labor shortage or lack of equipment for PM maintenance makes it difficult to perform the PM process on different machines simultaneously [4]. Since the maintenance team is unable to perform more than one maintenance at a time, the maintenance schedule should be changed so that no more than one maintenance is scheduled at a time [5].Therefore, in cases where the machines’ maintenance schedules overlap, we need to change the optimal time of the parts replacement by modifying the schedule. Hence, our goal is to create a trade-off between the optimal time for preventive maintenance and the limitations of labor and equipment in order to create an executable schedule at the lowest cost.Preventive maintenance (PM) is an old well-known problem in the manufacturing world with a proven effect on production scheduling and output [6, 7]. It had been studied from the last century, starting from the policy development of replacement which was based on the limit for the repair cost to more complex methods such as the joint optimization of predictive maintenance planning and production scheduling [8].There are a few approaches to how PM problems can be solved. The first one is to create policies/PM schedules by maximizing system availability. It has been done by using the different methods in [9–13]. Most of the existing research on PM overlooks the interrelationship between maintenance planning and equipment work schedules. As a result, maintenance planning is often done without considering the interactions between these two activities [14, 15]. The most recent research focuses on combining production scheduling with PM scheduling. Thus, the second approach to solve a PM problem is to create a PM scheduling model that considers the production schedule of the machines as the constraints [16, 17]. Another approach is to build the production scheduling model with the PM constraints, such as in [18]. Furthermore, the joint optimization of preventive maintenance and production scheduling planning gained a lot of interest. The goal is to minimize all the costs (maintenance and production-related) and satisfy the demand while creating the optimal production and preventive maintenance schedule [19–21]. Tardiness cost can also have an impact on the production and PM planning and was added to the objective function by [22, 23]. Additionally, the problem can be extended by considering the effect of maintenance on the quality of the products or simultaneously making the quality decision together with the production and maintenance.Any scheduling problem is known as a complex problem that includes multiple objectives that have to be minimized or maximized at the same time. However, multiple scheduling criteria are rarely considered simultaneously in the literature [24] but, in recent years, it has gained more interest through a trade-off between conflicting objectives such as in [25, 26]. Multiple objectives are typically being used in such problems as the assignment-allocation problem [27] or supply chain network design problems [28]. In production scheduling models, it is very important to find the balance between the time and the total production cost [29, 80]. To buttress this point, [81] uses the multiobjective evolutionary algorithm to minimize the total flow time of jobs and the number of tardy jobs. PM scheduling and rescheduling are also done by minimizing total operational cost (job’s total completion times, maintenance cost, and compression cost) and total completion time deviation simultaneously [30].PM problem is categorized as an NP-hard problem as any scheduling model with a lot of constraints and large datasets. Metaheuristic is a to-go method for many scheduling/assignment problems that provides very good results in a short period of time that is a critical concern in the industry. Metaheuristic methods are a wide class of algorithms that can conduct searching phases based on their stochastic cores [31–41]. They have found their place among viral methods to be utilized for intelligent systems [42–46] and hybrid pattern recognition works [47–55]. These methods can be basic or enhanced methods with other evolutionary bases, and some examples include evolutionary algorithms, ant colony optimization (ACO) [36, 47], memetic and hybrid algorithms [56–60], tabu search, simulated annealing, etc. [47, 61]. Evolutionary algorithms have been found to be very successful in solving multiobjective optimization problems. They also have the potential to be integrated with the hybrid and new machine learning models in the future [62–67]. In the past decade, the nondominated sorting genetic algorithm (NSGA-II) has been one of the most popular and practical evolutionary multiobjective optimization (EMO) algorithms [32, 68–70].There are many modified versions of heuristics algorithms that solve the scheduling and assignment problems, such as [35, 52, 55]. For example, in many studies, modified versions of NSGA-II were utilized in order to create algorithms to find more accurate solutions because the regular NSGA-II algorithms were not able to find it [71–74]. The other modified algorithms, such as SPEA-II and MOPSO, were generated in order to find solutions for particular problems which could not be solved by regular algorithms [36, 39, 75, 76].In this study, we first review scheduling maintenance and cost optimization problems and integrate these two problems with other related problems in a series-parallel system. We convert these problems to a mathematical model and propose a modified-NSGA-II algorithm to get optimal solutions. Since many of the generated solutions are not acceptable due to having overlap with other solutions, this algorithm has a repository to keep these unacceptable solutions, which are changed from unacceptable to acceptable solutions by a mutation operator. Therefore, the performance of solution generation in the algorithm is increased. This mutation operator includes a complicated procedure to change the start times of preventive maintenance in solutions in order to convert unacceptable solutions to acceptable solutions in the repository.In this modified NSGA-II algorithm, we propose to add a repository with a new mutation operator algorithm. Then we analyze the comparison of the results obtained from this algorithm with two regular algorithms.


## 2. Problems

### 2.1. Maintenance Scheduling Problem

Let us consider a factory in which several machines are working. Each piece of equipment has parts that wear out (depreciate) over time and need to be repaired or replaced. Obviously, the replacement time of these parts is different. In addition, the time that the maintenance team spends to replace each of these parts varies. Therefore, considering the number of machines in operation and the several parts that must be repaired or replaced, the possibility of scheduling two different maintenance at the same time is very high. Due to the limited labor and equipment of the maintenance team, it is necessary to change the preventive replacement schedule when two or more maintenance coincides with each other. Potentially, it can help to prevent the accumulation of preventive maintenance over a period of time.If we assume the total number of PM actions in a period of time ism so that the maintenance is arranged in order of its start time from 1 to m, then assuming that n is a number between 0 and m, the connection of the end time of the maintenance n with the start time of maintenance n + 1 can have 3 modes.In the first case, as soon as the maintenancen is completed, the maintenance n + 1 begins. Since no time is wasted by the maintenance team between the two maintenance, it is an ideal situation.The second case is when the period of time between PM actionsn and n + 1 overlaps. Indeed, the maintenance, n + 1, begins before the completion of the maintenance, n. Due to the accumulation of the PM program and the limited capacity of the maintenance team, if the volume of maintenance work in a period of time is more than the capacity of the maintenance team, it becomes impossible to perform and the schedule would have to be modified. To solve this problem, we need to shift the maintenance n slightly backward or the maintenance n + 1 slightly forward so that the overlap is eliminated and the maintenance planning becomes feasible.In the last case, we have a gap or time interval between the end time of the maintenancen and the start time of the maintenance n + 1, which means the maintenance team has to be idle during this period. This case is acceptable for the maintenance team, but it is not ideal due to the gap. Since our goal is to perform all periodic maintenance within a defined time frame, the presence of different gaps can waste time and cost. Hence, to create the optimal schedule, we can eliminate or reduce the gap or time interval between the two PM actions to get closer to the ideal state.


### 2.2. Cost Optimization and Failure Rate

A machine that works in a factory has parts that wear out over time. These parts need to be replaced periodically. Since for most of the parts, the possibility of failure increases over time with aging, the late replacement of each part escalates the probability of failure of that part. Failures may result in damage to other parts and components of the machine, which causes additional costs and inactivity of the machine for a period. On the other hand, the early replacement of defective parts can impose additional costs such as increased utilization of spare parts, the maintenance team’s wages, and inactivity of the machine during the replacement period on the factory. Therefore, finding the exact time to replace defective parts is very important in reducing costs.


### 2.3. Production Problem in a Series-Parallel System

As mentioned earlier, keeping the production rate of machines constant is the main goal of the maintenance team. In a series-parallel system shown in Figure1, stopping a machine due to a failure or stopping a manufacturing process due to PM activities will reduce production. This reduction rate or RS is determined as follows:(1)Rs=RiRT,where Ri is the amount of reduction in production and RT is the total production rate of machines in one of subsystems 1 to N while all the machines in that subsystem are working and it can be formulated as follows:(2)RT=∑i=1nRn.Figure 1
Series-parallel systems.Therefore,RS is a number between 0 and 1.


### 2.4. Trade-Off between Problems

As discussed in Section2.2, early and late replacement of parts can increase costs. The exact time of replacement of each part can be calculated. The problem occurs when the maintenance team has limited working time. Therefore, it will be difficult to execute the schedule accurately. On the other hand, as mentioned in Section 2.3, keeping the production rate constant is also one of the goals of the maintenance unit.Therefore, a trade-off between a production rate and PM costs is needed considering the limitations of labor and equipment of the maintenance team. In order to obtain this trade-off, a mathematical model is required.


## 2.1. Maintenance Scheduling Problem

Let us consider a factory in which several machines are working. Each piece of equipment has parts that wear out (depreciate) over time and need to be repaired or replaced. Obviously, the replacement time of these parts is different. In addition, the time that the maintenance team spends to replace each of these parts varies. Therefore, considering the number of machines in operation and the several parts that must be repaired or replaced, the possibility of scheduling two different maintenance at the same time is very high. Due to the limited labor and equipment of the maintenance team, it is necessary to change the preventive replacement schedule when two or more maintenance coincides with each other. Potentially, it can help to prevent the accumulation of preventive maintenance over a period of time.If we assume the total number of PM actions in a period of time ism so that the maintenance is arranged in order of its start time from 1 to m, then assuming that n is a number between 0 and m, the connection of the end time of the maintenance n with the start time of maintenance n + 1 can have 3 modes.In the first case, as soon as the maintenancen is completed, the maintenance n + 1 begins. Since no time is wasted by the maintenance team between the two maintenance, it is an ideal situation.The second case is when the period of time between PM actionsn and n + 1 overlaps. Indeed, the maintenance, n + 1, begins before the completion of the maintenance, n. Due to the accumulation of the PM program and the limited capacity of the maintenance team, if the volume of maintenance work in a period of time is more than the capacity of the maintenance team, it becomes impossible to perform and the schedule would have to be modified. To solve this problem, we need to shift the maintenance n slightly backward or the maintenance n + 1 slightly forward so that the overlap is eliminated and the maintenance planning becomes feasible.In the last case, we have a gap or time interval between the end time of the maintenancen and the start time of the maintenance n + 1, which means the maintenance team has to be idle during this period. This case is acceptable for the maintenance team, but it is not ideal due to the gap. Since our goal is to perform all periodic maintenance within a defined time frame, the presence of different gaps can waste time and cost. Hence, to create the optimal schedule, we can eliminate or reduce the gap or time interval between the two PM actions to get closer to the ideal state.


## 2.2. Cost Optimization and Failure Rate

A machine that works in a factory has parts that wear out over time. These parts need to be replaced periodically. Since for most of the parts, the possibility of failure increases over time with aging, the late replacement of each part escalates the probability of failure of that part. Failures may result in damage to other parts and components of the machine, which causes additional costs and inactivity of the machine for a period. On the other hand, the early replacement of defective parts can impose additional costs such as increased utilization of spare parts, the maintenance team’s wages, and inactivity of the machine during the replacement period on the factory. Therefore, finding the exact time to replace defective parts is very important in reducing costs.


## 2.3. Production Problem in a Series-Parallel System

As mentioned earlier, keeping the production rate of machines constant is the main goal of the maintenance team. In a series-parallel system shown in Figure1, stopping a machine due to a failure or stopping a manufacturing process due to PM activities will reduce production. This reduction rate or RS is determined as follows:(1)Rs=RiRT,where Ri is the amount of reduction in production and RT is the total production rate of machines in one of subsystems 1 to N while all the machines in that subsystem are working and it can be formulated as follows:(2)RT=∑i=1nRn.Figure 1
Series-parallel systems.Therefore,RS is a number between 0 and 1.


## 2.4. Trade-Off between Problems

As discussed in Section2.2, early and late replacement of parts can increase costs. The exact time of replacement of each part can be calculated. The problem occurs when the maintenance team has limited working time. Therefore, it will be difficult to execute the schedule accurately. On the other hand, as mentioned in Section 2.3, keeping the production rate constant is also one of the goals of the maintenance unit.Therefore, a trade-off between a production rate and PM costs is needed considering the limitations of labor and equipment of the maintenance team. In order to obtain this trade-off, a mathematical model is required.


## 3. Mathematical Model

In this paper, to create a mathematical model, we begin with identifying the problem assumptions and stating the constraints. Then the probability of failure of parts in the machine is calculated and the objective function by specifying the decision variables and the parameters are defined.The assumptions of the problem are based on the following:(i)
The maintenance team is not able to perform a PM action on two machines at the same time(ii)
The daily and hourly efficiency of the maintenance team is constant(iii)
The maintenance team is available to perform PM action during all working hours of the company(iv)
The spare parts are available, and there are no restrictions on access to the required spare parts(v)
All required PM actions are performed within a specified time framewhere Ti,j,k is the time allocated for PM action in the kth turn, to repair the part jth of the machine ith.


### 3.1. Identification of Failure Rate

We need an exponential function to calculate the failure time of parts in machines. The Weibull function [77] is used to check the failure time in the depreciation state. The Weibull function is a function with continuous values. The values of this function are nonnegative real numbers, so this distribution can be used in cases where the random variable is related to longevity. Therefore, to find the failure probability density function, the cumulative distribution function (CDF) would be a good option, which can be computed as follows:(3)ry=1−e−y/λk.The part’s age plays a decisive role in calculating the probability of part failure. Ifr (y) is a function of the failure rate of the part, the probability of failure after time L is obtained from the following equation:(4)Fy=∫0Lrydx,where L is the age of the part at time t, and the age of the part at the beginning of the period is equal to 0 due to PM action.


### 3.2. Identification of Decision Variables and Objective Functions

The decision variable in this problem isLi,j,k, which is a matrix with i rows, j columns, and k dimension. In this matrix, i is the number of machines in operation, j is the number of parts of a machine that need preventive maintenance, and k is the number of times parts that need to be repaired. The number shown in each cell of this matrix is assigned to the start time of the PM action.The duration of PM process for each part of the machine is a fixed number and it is indicated byTi,j. The cost of the repair team for each hour of work on the ith machine and jth part is specified as Pi,j, the cost of PM for the part jth of the machine ith is obtained from the following equation:(5)CPMi,j=Ti,j×Pi,j+Mi,j+Si,j,where Si,j is the price of spare part jth of the machine ith and Mi,j is the cost of one hour of inactivity of the machine ith.Then the cost objective function is formulated as follows:(6)minimizeZ=∑i=1I∑j=1J∑k=1KCPMi,j+∫0Li,j,kri,jydy×Ci,j,where Ci,j is the cost of part j’s breakdown of the machine ith.Since inactivity of the production line due to failure or PM process leads to decrease in the volume of the production, time objective function that can be formulated as follows:(7)minimizeZ=∑i=1I∑j=1J∑k=1K∫0Li,j,kri,jydy×Ris×Ti,j′+Qi,j˙×Ti,j×Ris,where Ris is the cost of inactivity for the machine, which has Ti,j′ hour inactivity for the machine i and the part j and Qi, j is the number of maintenance actions performed on the machine i and the part j.s.t.(8)∑i=1I∑j=1J∑k=1KTi,j,k≤Ttotal,i,j∈N,(9)∑i=1ITi′≤Fi,i∈N,(10)At∈0∪N,∀t∈Ttotal.The total number of hours that the maintenance team is available to perform a PM action in the company is defined byTtotal. Then we determine how many preventive maintenance activities should be done for each part of a single machine per year. It should be noted that each maintenance team has a certain capacity in terms of labor and equipment. Hence, the total time spent in a certain period should not exceed the capacity of the maintenance team. This period is usually considered to be annual—this constraint model as equation (8).The constraint modeled by equation (9) allocates to the inactivity time of machine i in the time interval 0 to Ttotal which cannot be longer than the specified value Fi. The last constraint specifies the number of assignments assigned to the maintenance team at time t, denoted by At.


## 3.1. Identification of Failure Rate

We need an exponential function to calculate the failure time of parts in machines. The Weibull function [77] is used to check the failure time in the depreciation state. The Weibull function is a function with continuous values. The values of this function are nonnegative real numbers, so this distribution can be used in cases where the random variable is related to longevity. Therefore, to find the failure probability density function, the cumulative distribution function (CDF) would be a good option, which can be computed as follows:(3)ry=1−e−y/λk.The part’s age plays a decisive role in calculating the probability of part failure. Ifr (y) is a function of the failure rate of the part, the probability of failure after time L is obtained from the following equation:(4)Fy=∫0Lrydx,where L is the age of the part at time t, and the age of the part at the beginning of the period is equal to 0 due to PM action.


## 3.2. Identification of Decision Variables and Objective Functions

The decision variable in this problem isLi,j,k, which is a matrix with i rows, j columns, and k dimension. In this matrix, i is the number of machines in operation, j is the number of parts of a machine that need preventive maintenance, and k is the number of times parts that need to be repaired. The number shown in each cell of this matrix is assigned to the start time of the PM action.The duration of PM process for each part of the machine is a fixed number and it is indicated byTi,j. The cost of the repair team for each hour of work on the ith machine and jth part is specified as Pi,j, the cost of PM for the part jth of the machine ith is obtained from the following equation:(5)CPMi,j=Ti,j×Pi,j+Mi,j+Si,j,where Si,j is the price of spare part jth of the machine ith and Mi,j is the cost of one hour of inactivity of the machine ith.Then the cost objective function is formulated as follows:(6)minimizeZ=∑i=1I∑j=1J∑k=1KCPMi,j+∫0Li,j,kri,jydy×Ci,j,where Ci,j is the cost of part j’s breakdown of the machine ith.Since inactivity of the production line due to failure or PM process leads to decrease in the volume of the production, time objective function that can be formulated as follows:(7)minimizeZ=∑i=1I∑j=1J∑k=1K∫0Li,j,kri,jydy×Ris×Ti,j′+Qi,j˙×Ti,j×Ris,where Ris is the cost of inactivity for the machine, which has Ti,j′ hour inactivity for the machine i and the part j and Qi, j is the number of maintenance actions performed on the machine i and the part j.s.t.(8)∑i=1I∑j=1J∑k=1KTi,j,k≤Ttotal,i,j∈N,(9)∑i=1ITi′≤Fi,i∈N,(10)At∈0∪N,∀t∈Ttotal.The total number of hours that the maintenance team is available to perform a PM action in the company is defined byTtotal. Then we determine how many preventive maintenance activities should be done for each part of a single machine per year. It should be noted that each maintenance team has a certain capacity in terms of labor and equipment. Hence, the total time spent in a certain period should not exceed the capacity of the maintenance team. This period is usually considered to be annual—this constraint model as equation (8).The constraint modeled by equation (9) allocates to the inactivity time of machine i in the time interval 0 to Ttotal which cannot be longer than the specified value Fi. The last constraint specifies the number of assignments assigned to the maintenance team at time t, denoted by At.


## 4. Methodology

On a production line, we often use a significant number of machines in the process, which affects the time and cost decisions of PM actions. In order to mitigate the computational difficulties of solving the resulting large-scale problem because of the nonlinear objective functions, we apply a modified NSGA-II heuristic algorithm. The modifications in the mutation and crossover operators are applied, which generate convergence and uniformly distributed solutions.


### 4.1. Initialization

Before creating the crossover and mutation operators, maintaining a diverse population of candidate solutions is required. This could subsequently be incorporated into a chromosome matrix. To solve the objective functions, the solutions created for this algorithm are in the form of a three-dimensional matrix. We assign a row that corresponds to each machinei in operation, a column that represents each part j of a machine in need of PM action, and a dimension that stands for time k that parts need to be repaired. The number shown in each element indicates the start time of the PM action for the machine i, the part j in the time k.


### 4.2. Simulated Binary Crossover (SBX)

The simulated binary crossover operator produces children in the vicinity of the parents [78]. In this operator, the parents are generated using the roulette wheel selection. Moreover, the two corresponding elements from both matrices are randomly selected and then for u, a random value between 0 and 1 is generated. Then, the value of β is calculated as follows:(11)β=2u1/η+1,ifu≤0.5,121−u1/η+1,otherwise,where ƞ is the distribution index, which is a real number.Two new children are produced using the following formula. We use the following equations to create new solution candidates:(12)x1=0.51+βxi+1−βxj,x2=0.51−βxi+1+βxj,where x1 and x2 are the children produced and xi and xj are the parents which are 2 elements of the decision variable matrix.


### 4.3. Polynomial Mutation

In polynomial mutation, it is possible to find a new solution at any distance from the parent, but it is more likely to be found around the parent than elsewhere [79]. To find the solution, first, the parameter δ is computed as follows:(13)δ=2r1/η−1−1,ifr<0.5,1−21−r1/η+1,ifr≥0.5,where r is a random number between 0 and 1, and ƞ is a positive real number that can directly control the probability distribution as an external parameter.The mutated child is obtained from equation (14), in which xU and xL are the upper and lower bonds, to extract a new solution from the mutation operator.(14)x2=x1+xU−xLδ.


### 4.4. Developed Mutation Algorithm (Mutation Type 2)

This mutation operator is designed to convert unacceptable solutions to acceptable ones. Unacceptable solutions are solutions that overlap at least once with the other solutions. To solve this problem, we define a mutation operator (Figure2). The mutation operator type 2 is designed by picking the first cell of the first row from the first column of the first aisle of the solutions matrix. If there is an overlap, the time interval of PM action for the cell increases or decreases to avoid any overlap with the other cells. Likewise, all cells in a row are reviewed to prevent overlapping between the cells. After reviewing one row, we will review the next row, and after reviewing all the rows of one aisle, we will review the other aisles until the last cell of the matrix is reviewed. After reviewing the last cell, the solution changes from an unacceptable solution to an acceptable solution.Figure 2
Developed mutation algorithm.


### 4.5. Nondominated Sorting

To sort the solutions, we use the convexity crowding distance method. For every chromosome, we have the following:(15)∑i=1Ifix−Ffimax−fimin,where fix is the ith objective function and fimax and fimin are the highest and lowest values of the ith objective function value, respectively. F is the closest value generated to the considered objective value. Solutions with higher distance crowding are better solutions.


### 4.6. Repository

Since the number of acceptable solutions generated per iteration is very small, if we use standard algorithms to solve this problem, many solutions will be unacceptable due to overlap. If these solutions are eliminated, the speed of the algorithm is greatly reduced, and it takes a long time for the algorithm to find the optimal solutions for the problem. So, to solve this problem, we are looking for a new method to convert unacceptable solutions to acceptable ones.The proposed method is to use a repository to store unacceptable solutions generated by mutation type 1 and crossover operators. Since there is a high probability that some of the solutions generated by mutation and crossover operators are unacceptable solutions because of the time overlap among solution matrix cells, these unacceptable generated solutions are stored in the repository and then after performing mutation operator type 2, they are converted to acceptable solutions. Because the time consumption of mutation operator type 2 is significant, not all solutions are sent to this operator. Therefore, after ranking solutions in the repository, a certain number of the best solutions are sent to mutation operator type 2 and then they are converted into acceptable solutions.The cycle of entering solutions into the repository and leaving solutions from the repository is shown in Figure3. The only repository input is the unacceptable solutions generated by crossover and mutation type 1 operators. Since the capacity of the repository is limited and it is not possible to store all unacceptable solutions generated in the repository, in each iteration, a number of lower-ranked solutions in the repository are removed so that the repository has enough space to store the solutions generated in the next iteration. The number of removed solutions in each iteration depends on the number of solutions entered into the repository and the number of solutions removed from the repository so that the number of solutions stored in the repository does not exceed the capacity of the repository.Figure 3
Repository cycle in MNSGA2 algorithm.The use of a repository in the NSGA2 algorithm also has disadvantages. The first drawback is to increase the computation time for each iteration of the algorithm. Although the presence of a repository in the algorithm causes more optimal solutions to be generated by the algorithm in each iteration, more significant time per iteration in the MNSGA2 algorithm is spent compared to the regular NSGA2 algorithm. It is because of the time required to perform operations in the type 2 mutation operator. A comparison between the time of each iteration between the algorithms is performed in the numerical example and result section. Another disadvantage of using this method is that many of the solutions generated by the type 2 mutation operator may no longer be the optimal solutions. The type 2 mutation operator converts unacceptable solutions into acceptable ones by changing the values of each solution matrix cell that overlaps with other cells in the same solution matrix, but these changes may result in to change the ranking of solutions.


### 4.7. Modified NSGA-II Algorithm

The idea behind the proposed algorithm is to increase the speed of generating optimal solutions. By adding a repository and a developed mutation algorithm (mutation type 2) to the regular NSGA-II algorithm, the modified NSGA-II algorithm is created, which has an extra loop compared to the regular NSGA-II algorithm. In the new algorithm, unacceptable optimal solutions that are in danger of being deleted because of having overlap between one cell in matrix solution and another cell in the same matrix solution are stored in the repository. These solutions are converted to acceptable solutions with the help of the mutation type 2 operator and then added to the candidate solutions. We conduct a modified NSGA-II algorithm that is illustrated in Figure4.Figure 4
Modified-NSGA-II algorithm.


## 4.1. Initialization

Before creating the crossover and mutation operators, maintaining a diverse population of candidate solutions is required. This could subsequently be incorporated into a chromosome matrix. To solve the objective functions, the solutions created for this algorithm are in the form of a three-dimensional matrix. We assign a row that corresponds to each machinei in operation, a column that represents each part j of a machine in need of PM action, and a dimension that stands for time k that parts need to be repaired. The number shown in each element indicates the start time of the PM action for the machine i, the part j in the time k.


## 4.2. Simulated Binary Crossover (SBX)

The simulated binary crossover operator produces children in the vicinity of the parents [78]. In this operator, the parents are generated using the roulette wheel selection. Moreover, the two corresponding elements from both matrices are randomly selected and then for u, a random value between 0 and 1 is generated. Then, the value of β is calculated as follows:(11)β=2u1/η+1,ifu≤0.5,121−u1/η+1,otherwise,where ƞ is the distribution index, which is a real number.Two new children are produced using the following formula. We use the following equations to create new solution candidates:(12)x1=0.51+βxi+1−βxj,x2=0.51−βxi+1+βxj,where x1 and x2 are the children produced and xi and xj are the parents which are 2 elements of the decision variable matrix.


## 4.3. Polynomial Mutation

In polynomial mutation, it is possible to find a new solution at any distance from the parent, but it is more likely to be found around the parent than elsewhere [79]. To find the solution, first, the parameter δ is computed as follows:(13)δ=2r1/η−1−1,ifr<0.5,1−21−r1/η+1,ifr≥0.5,where r is a random number between 0 and 1, and ƞ is a positive real number that can directly control the probability distribution as an external parameter.The mutated child is obtained from equation (14), in which xU and xL are the upper and lower bonds, to extract a new solution from the mutation operator.(14)x2=x1+xU−xLδ.


## 4.4. Developed Mutation Algorithm (Mutation Type 2)

This mutation operator is designed to convert unacceptable solutions to acceptable ones. Unacceptable solutions are solutions that overlap at least once with the other solutions. To solve this problem, we define a mutation operator (Figure2). The mutation operator type 2 is designed by picking the first cell of the first row from the first column of the first aisle of the solutions matrix. If there is an overlap, the time interval of PM action for the cell increases or decreases to avoid any overlap with the other cells. Likewise, all cells in a row are reviewed to prevent overlapping between the cells. After reviewing one row, we will review the next row, and after reviewing all the rows of one aisle, we will review the other aisles until the last cell of the matrix is reviewed. After reviewing the last cell, the solution changes from an unacceptable solution to an acceptable solution.Figure 2
Developed mutation algorithm.


## 4.5. Nondominated Sorting

To sort the solutions, we use the convexity crowding distance method. For every chromosome, we have the following:(15)∑i=1Ifix−Ffimax−fimin,where fix is the ith objective function and fimax and fimin are the highest and lowest values of the ith objective function value, respectively. F is the closest value generated to the considered objective value. Solutions with higher distance crowding are better solutions.


## 4.6. Repository

Since the number of acceptable solutions generated per iteration is very small, if we use standard algorithms to solve this problem, many solutions will be unacceptable due to overlap. If these solutions are eliminated, the speed of the algorithm is greatly reduced, and it takes a long time for the algorithm to find the optimal solutions for the problem. So, to solve this problem, we are looking for a new method to convert unacceptable solutions to acceptable ones.The proposed method is to use a repository to store unacceptable solutions generated by mutation type 1 and crossover operators. Since there is a high probability that some of the solutions generated by mutation and crossover operators are unacceptable solutions because of the time overlap among solution matrix cells, these unacceptable generated solutions are stored in the repository and then after performing mutation operator type 2, they are converted to acceptable solutions. Because the time consumption of mutation operator type 2 is significant, not all solutions are sent to this operator. Therefore, after ranking solutions in the repository, a certain number of the best solutions are sent to mutation operator type 2 and then they are converted into acceptable solutions.The cycle of entering solutions into the repository and leaving solutions from the repository is shown in Figure3. The only repository input is the unacceptable solutions generated by crossover and mutation type 1 operators. Since the capacity of the repository is limited and it is not possible to store all unacceptable solutions generated in the repository, in each iteration, a number of lower-ranked solutions in the repository are removed so that the repository has enough space to store the solutions generated in the next iteration. The number of removed solutions in each iteration depends on the number of solutions entered into the repository and the number of solutions removed from the repository so that the number of solutions stored in the repository does not exceed the capacity of the repository.Figure 3
Repository cycle in MNSGA2 algorithm.The use of a repository in the NSGA2 algorithm also has disadvantages. The first drawback is to increase the computation time for each iteration of the algorithm. Although the presence of a repository in the algorithm causes more optimal solutions to be generated by the algorithm in each iteration, more significant time per iteration in the MNSGA2 algorithm is spent compared to the regular NSGA2 algorithm. It is because of the time required to perform operations in the type 2 mutation operator. A comparison between the time of each iteration between the algorithms is performed in the numerical example and result section. Another disadvantage of using this method is that many of the solutions generated by the type 2 mutation operator may no longer be the optimal solutions. The type 2 mutation operator converts unacceptable solutions into acceptable ones by changing the values of each solution matrix cell that overlaps with other cells in the same solution matrix, but these changes may result in to change the ranking of solutions.


## 4.7. Modified NSGA-II Algorithm

The idea behind the proposed algorithm is to increase the speed of generating optimal solutions. By adding a repository and a developed mutation algorithm (mutation type 2) to the regular NSGA-II algorithm, the modified NSGA-II algorithm is created, which has an extra loop compared to the regular NSGA-II algorithm. In the new algorithm, unacceptable optimal solutions that are in danger of being deleted because of having overlap between one cell in matrix solution and another cell in the same matrix solution are stored in the repository. These solutions are converted to acceptable solutions with the help of the mutation type 2 operator and then added to the candidate solutions. We conduct a modified NSGA-II algorithm that is illustrated in Figure4.Figure 4
Modified-NSGA-II algorithm.


## 5. Numerical Example and Result

We used the real data collected from a company that produced car spare parts for the computational experiments. In this company, there are 40 machines that are producing metal parts in a production line. This production line has 10 stations where there are 4 machines in each station. All machines are milling or turning machine tools that are under preventive maintenance scheduling.These machines have some moving parts such as spindle, tool change, turret, pallet change, and moving table, all of which are depreciated due to their mobility. In addition to the mobility of these parts, machining vibration, corrosion of parts due to use of cooling fluid, and damage of machine parts because of metal chips increase the depreciation of machine parts. Typically, ball bearings, spindles, various parts of the cooling system, tool change, and pallet change are subject to preventive maintenance.


### 5.1. Results for Modified-NSGA-II, NSGA-II, and SPEA-II

We initialized the algorithm by randomly generating a population of 100 solutions. The number of the initial population and the solutions provided by mutation and crossover operators, repository capacity, and the maximum number of unacceptable solutions sent from repository to operator mutation type 2 would be different based on their importance on the problem. For this problem, values of algorithm and repository specifications are demonstrated in Table1.Table 1
Algorithm and repository specifications.
Algorithm specificationNumber of produced populationsInitial population100Mutation (type 1)30Crossover30Repository capacity500Maximum number of unacceptable solutions sent from repository to operator mutation type 2 in each iteration20To assess the performance of our Modified-NSGA-II algorithm and compare it with the two algorithms NSGA-II and SPEA-II, we defined a reference set based on 100 best candidate solutions obtained by the Modified-NSGA-II, NSGA-II, and SPEA-II. Each was replicated 1000 times after every 100 iterations. Table2 illustrates the number of nondominated solutions (NNS) by reference in all three algorithms, and as shown in Table 2, the Modified-NSGA-II algorithm has more nondominated solutions than the other two algorithms. We conducted the experiment using the MATLAB platform on a Windows-based server with 16 GB RAM, i7 CPU, and 1.8GHs.Table 2
The number of NNS for MNSGA-II, NSGA-II, and SPEA-II.
GenerationMNSGA-IINSGA-IISPEA210088242008701300800040076005007700600720070072008007200900730010007300For further assessment, we compared the performance of the Modified-NSGA-II algorithm with the regular NSGA-II and SPEA-II algorithms using three methods: inverted generated distance, distribution metric, and spacing metric. The results of the comparison among the MNSGA-II, the NSGA-II, and the SPEA-II can be seen in Table3.Table 3
Comparison of three algorithms.
AlgorithmIGDDeltaSpacingMNSGA-II24490076846NSGA-II232400015393408SPEA2781400063986110


### 5.2. Inverted Generated Distance

This approach is designed to assess diversity and convergence and can be calculated as follows:(16)IGD=∑v∈P∗dν,PP∗,where P∗ is equal to uniformly distributed points in a true Pareto front, P is a nondominated solution obtained by a selected algorithm, ν is a solution that belongs to P∗, and ⅆ (ν, P) is equal to the minimum Euclidean distance between ν and the point in P. The results shown in Table 3 illustrate that, in general, the modified-NSGA-II algorithm has lower IGD values than the other two algorithms.


### 5.3. Delta Index

In the delta approach, the distribution of the solutions in the Pareto front is examined considering the following equation:(17)Δs=∑i=1s−1di−d¯s−1,where di is the Euclidean distance between two consecutive solutions of Pareto front with s optimal solutions and d¯is the average di. The lower the value of delta is, the more uniformly the Pareto front is distributed. According to the values specified in Table 3, with a slight difference after the SPEA-II algorithm, the modified-NSGA-II algorithm has the lowest value.


### 5.4. Spacing Metrics

The uniformity of the distance between the Pareto front and the reference set is analyzed in the Spacing Metric method. A lower space parameter (SP) value indicates a more uniform distance between the Pareto front and the reference. The SP value is calculated as follows:(18)SPs=1S−1∑i=1sd¯−di2.In this case,di is the minimum distance of solution i from the reference point. d¯is the average of di for s optimal solution set. According to the values specified in Table 3, with a significant difference compared to the other two algorithms, the SPEA2 algorithm has a uniform distance between the reference and Pareto front.


### 5.5. Operation Time

As mentioned in Section4.6, although the MNSGA2 algorithm generates better solutions per iteration compared to the SPEA2 and NSGA2 algorithms due to the existence of a repository, it requires more time to perform the calculations. The more complexity the type 2 mutation operator has, the more time is required to perform calculations. Also, the more machines and parts that need preventive maintenance, the more cells the solution matrix has. In this case, the probability of overlap’s existence is higher, and the probability of unacceptable solutions’ existence increases as well. Also, more calculations in the type 2 mutation operator are required to make changes in the solution matrix and produce an acceptable solution. The increase in the computational volume in the type 2 mutation operator increases the whole operation time of the algorithm for each iteration. Figure 5 represents the computation time for this problem between MNSGA2, NSGA2, and SPEA2 algorithms. In Figure 5, the average operation time of NSGA2, MNSGA2, and SPEA2 algorithms is shown for 10 operations so that each one has 1000 iterations.Figure 5
Operation time of algorithms.


### 5.6. Fitness of Generated Solutions

Another problem that the MNSGA2 algorithm has in comparison with the two NSGA2 and SPEA2 algorithms is the change in the fitness or rank of the solutions after converting the unacceptable solutions into the acceptable ones in the type 2 mutation operator. As mentioned earlier, the type 2 mutation operator converts unacceptable solutions into acceptable ones by changing the values of the solution matrix cells that have overlapped with other cells, but after examining the fitness of the solutions, it becomes clear that the values obtained may no longer be our optimal solutions and may be dominated by other solutions obtained in the Pareto front. These solutions are removed after ranking and do not help speed up the algorithm to find optimal solutions. In this case, 73.6% of the total solutions generated by the mutation type 1 and crossover operators were unacceptable, and 68.7% of the solutions generated by the type 2 mutation operator were dominated by other Pareto front solutions. These results are obtained from the average 10 times operation of the MNSGA2 algorithm, which has 1000 iterations each operation.


## 5.1. Results for Modified-NSGA-II, NSGA-II, and SPEA-II

We initialized the algorithm by randomly generating a population of 100 solutions. The number of the initial population and the solutions provided by mutation and crossover operators, repository capacity, and the maximum number of unacceptable solutions sent from repository to operator mutation type 2 would be different based on their importance on the problem. For this problem, values of algorithm and repository specifications are demonstrated in Table1.Table 1
Algorithm and repository specifications.
Algorithm specificationNumber of produced populationsInitial population100Mutation (type 1)30Crossover30Repository capacity500Maximum number of unacceptable solutions sent from repository to operator mutation type 2 in each iteration20To assess the performance of our Modified-NSGA-II algorithm and compare it with the two algorithms NSGA-II and SPEA-II, we defined a reference set based on 100 best candidate solutions obtained by the Modified-NSGA-II, NSGA-II, and SPEA-II. Each was replicated 1000 times after every 100 iterations. Table2 illustrates the number of nondominated solutions (NNS) by reference in all three algorithms, and as shown in Table 2, the Modified-NSGA-II algorithm has more nondominated solutions than the other two algorithms. We conducted the experiment using the MATLAB platform on a Windows-based server with 16 GB RAM, i7 CPU, and 1.8GHs.Table 2
The number of NNS for MNSGA-II, NSGA-II, and SPEA-II.
GenerationMNSGA-IINSGA-IISPEA210088242008701300800040076005007700600720070072008007200900730010007300For further assessment, we compared the performance of the Modified-NSGA-II algorithm with the regular NSGA-II and SPEA-II algorithms using three methods: inverted generated distance, distribution metric, and spacing metric. The results of the comparison among the MNSGA-II, the NSGA-II, and the SPEA-II can be seen in Table3.Table 3
Comparison of three algorithms.
AlgorithmIGDDeltaSpacingMNSGA-II24490076846NSGA-II232400015393408SPEA2781400063986110


## 5.2. Inverted Generated Distance

This approach is designed to assess diversity and convergence and can be calculated as follows:(16)IGD=∑v∈P∗dν,PP∗,where P∗ is equal to uniformly distributed points in a true Pareto front, P is a nondominated solution obtained by a selected algorithm, ν is a solution that belongs to P∗, and ⅆ (ν, P) is equal to the minimum Euclidean distance between ν and the point in P. The results shown in Table 3 illustrate that, in general, the modified-NSGA-II algorithm has lower IGD values than the other two algorithms.


## 5.3. Delta Index

In the delta approach, the distribution of the solutions in the Pareto front is examined considering the following equation:(17)Δs=∑i=1s−1di−d¯s−1,where di is the Euclidean distance between two consecutive solutions of Pareto front with s optimal solutions and d¯is the average di. The lower the value of delta is, the more uniformly the Pareto front is distributed. According to the values specified in Table 3, with a slight difference after the SPEA-II algorithm, the modified-NSGA-II algorithm has the lowest value.


## 5.4. Spacing Metrics

The uniformity of the distance between the Pareto front and the reference set is analyzed in the Spacing Metric method. A lower space parameter (SP) value indicates a more uniform distance between the Pareto front and the reference. The SP value is calculated as follows:(18)SPs=1S−1∑i=1sd¯−di2.In this case,di is the minimum distance of solution i from the reference point. d¯is the average of di for s optimal solution set. According to the values specified in Table 3, with a significant difference compared to the other two algorithms, the SPEA2 algorithm has a uniform distance between the reference and Pareto front.


## 5.5. Operation Time

As mentioned in Section4.6, although the MNSGA2 algorithm generates better solutions per iteration compared to the SPEA2 and NSGA2 algorithms due to the existence of a repository, it requires more time to perform the calculations. The more complexity the type 2 mutation operator has, the more time is required to perform calculations. Also, the more machines and parts that need preventive maintenance, the more cells the solution matrix has. In this case, the probability of overlap’s existence is higher, and the probability of unacceptable solutions’ existence increases as well. Also, more calculations in the type 2 mutation operator are required to make changes in the solution matrix and produce an acceptable solution. The increase in the computational volume in the type 2 mutation operator increases the whole operation time of the algorithm for each iteration. Figure 5 represents the computation time for this problem between MNSGA2, NSGA2, and SPEA2 algorithms. In Figure 5, the average operation time of NSGA2, MNSGA2, and SPEA2 algorithms is shown for 10 operations so that each one has 1000 iterations.Figure 5
Operation time of algorithms.


## 5.6. Fitness of Generated Solutions

Another problem that the MNSGA2 algorithm has in comparison with the two NSGA2 and SPEA2 algorithms is the change in the fitness or rank of the solutions after converting the unacceptable solutions into the acceptable ones in the type 2 mutation operator. As mentioned earlier, the type 2 mutation operator converts unacceptable solutions into acceptable ones by changing the values of the solution matrix cells that have overlapped with other cells, but after examining the fitness of the solutions, it becomes clear that the values obtained may no longer be our optimal solutions and may be dominated by other solutions obtained in the Pareto front. These solutions are removed after ranking and do not help speed up the algorithm to find optimal solutions. In this case, 73.6% of the total solutions generated by the mutation type 1 and crossover operators were unacceptable, and 68.7% of the solutions generated by the type 2 mutation operator were dominated by other Pareto front solutions. These results are obtained from the average 10 times operation of the MNSGA2 algorithm, which has 1000 iterations each operation.


## 6. Conclusion

In this study, we developed a model to solve the scheduling problem in preventive maintenance for series-parallel systems. In this model, we addressed the problems in PM process scheduling in series-parallel systems and then turned this integrated model into a mathematical model. In this mathematical model, we considered a matrix as a decision variable where each cell represents the time scheduling of the PM process for a specific part of a single machine. In the future, the scale of the mathematical model can be changed by adding more decision variables to obtain more results. Then, we defined the objective functions in order to find the optimal solutions. Moreover, we proposed a modified genetic algorithm with customized mutation and crossover operators that allowed us to find a better solution compared to the existing algorithms. We implemented our experiments on a case study drawn from the real data then analyzed and compared the solutions.Our results confirm the better performance of the MNSGA-II algorithm, which proposes that evolutionary algorithms alone do not necessarily provide more viable optimal solutions and need customized operators for better computational efficiency. In our future works, we will consider more objective functions like the scheduling of special equipment’s use in the PM process, the spare parts supply, and the employment of specially trained personnel by utilizing the modified NSGA-II that presented in this study. Analysis of the scheduling time of PM action provides manufacturing companies with valuable information to determine meticulous scheduling decisions and planning, prevent production interruptions, and increase the amount of production.



---

*Source: 5583125-2021-07-01.xml*