# A New Biobjective Model to Optimize Integrated Redundancy Allocation and Reliability-Centered Maintenance Problems in a System Using Metaheuristics

Authors: Shima MohammadZadeh Dogahe; Seyed Jafar Sadjadi
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A novel integrated model is proposed to optimize the redundancy allocation problem (RAP) and the reliability-centered maintenance (RCM) simultaneously. A system of both repairable and nonrepairable components has been considered. In this system, electronic components are nonrepairable while mechanical components are mostly repairable. For nonrepairable components, a redundancy allocation problem is dealt with to determine optimal redundancy strategy and number of redundant components to be implemented in each subsystem. In addition, a maintenance scheduling problem is considered for repairable components in order to identify the best maintenance policy and optimize system reliability. Both active and cold standby redundancy strategies have been taken into account for electronic components. Also, net present value of the secondary cost including operational and maintenance costs has been calculated. The problem is formulated as a biobjective mathematical programming model aiming to reach a tradeoff between system reliability and cost. Three metaheuristic algorithms are employed to solve the proposed model: Nondominated Sorting Genetic Algorithm (NSGA-II), Multiobjective Particle Swarm Optimization (MOPSO), and Multiobjective Firefly Algorithm (MOFA). Several test problems are solved using the mentioned algorithms to test efficiency and effectiveness of the solution approaches and obtained results are analyzed.


---

## Body

## 1. Introduction and Literature Review

In general, reliability is defined as ability of a system to meet required performance standards under specified conditions during a determined time horizon. It has a significant effect on manufacturing cost, company’s fame, production efficiency and environment, and so forth. There are two main approaches to enhance system reliability: implementing a proper maintenance policy and using effective redundancy strategies. Applying these approaches leads to increases in system reliability along with increasing costs of other resources. Thus, reaching a tradeoff between system reliability, cost, volume, weight, and so forth is significantly important [1].There are two types of maintenance policies: corrective maintenance (CM) and preventive maintenance (PM). In corrective maintenance, the system is repaired or replaced after failure. However, prescheduled periodic maintenance actions are taken in preventive maintenance. It is obvious that preventive maintenance prevents major failures that impose high costs on the system. In recent years, many authors conducted variety of research on preventive maintenance scheduling problems. A mathematical model has been proposed by Goel and Gupta [2] for a multistate system with repair and replacement policy. Goel et al. [3] presented a formulation for design, production, and maintenance planning to incorporate the reliability allocation problem at the design stage. Then, a simultaneous optimization framework has been employed to solve the proposed model. Tsai et al. [4] studied the preventive maintenance scheduling problem for a multicomponent system by assuming three maintenance actions: mechanical service, repair, and replacement. Periodic preventive maintenance actions have been considered in order to maximize availability of the system. Mohanta et al. [5] have employed both GA and hybrid GA/SA techniques to optimize maintenance scheduling for a power plant and compared the obtained results by the algorithms.Martorell et al. [6] proposed a multiobjective model to optimize the maintenance scheduling problem by integrating human and material resources. Reliability, availability, maintainability, and cost were considered as objectives of the model and the Genetic Algorithm was used to solve the problem. Another multiobjective maintenance model for a series-parallel system was investigated by Certa et al. [7]. They implemented periodic PM policy and considered maintenance cost and makespan as objectives of the problem. An effective Pareto optimal frontier approach was applied to solve the multiobjective problem. Moghaddass et al. [8] focused on finding an optimal tradeoff between design of a repairable multistate system with binary components and its maintenance strategy. Also, they considered both active and standby redundancy strategies. Doostparast et al. [9] developed a reliability-based approach to optimize preventive maintenance scheduling in coherent systems. A system is called coherent when its performance is related to all of the components. In other words, for coherent systems each component is relevant and system structure function is monotone and nondecreasing [10]. They studied periodic PM performance in three types of coherent systems and used a Simulated Annealing (SA) algorithm to solve the problem trying to minimize total costs along with meeting the minimum predetermined reliability level. Several other studies dealing with preventive maintenance scheduling problem have been done by [11–13] during the last two decades.However, selecting a proper maintenance policy is not all we can do to maximize system reliability. Identifying and implementing the best redundancy strategy is another way to optimize system reliability. One of the famous problems in field of reliability optimization is redundancy allocation problem (RAP). Redundant components are incorporated into the system to back up different parts of the system and prevent system breakdown under different redundancy strategies. There are two main redundancy strategies: (1) active redundancy, in which all redundant components are implemented in a parallel structure together from time zero and only one component is required to work at any given time, and (2) standby redundancy, in which a sequential order is determined for using the redundant components at component failure time. Three variants of the standby redundancy strategy are called cold, warm, and hot. Each strategy can be implemented in a different part of a system. RAPs are proved to be NP-hard by Chern [14]. Therefore, metaheuristic algorithms have been widely used in the literature to solve such problems. Literature of the redundancy allocation problem (RAP) could be reviewed from several points of view. In this paper, literature of the problem has been accurately reviewed by assuming three main characteristics: objectives of the problem, applied solution algorithms, and considered redundancy strategies.Coit [15] studied the redundancy optimization problem using integer programming approach and logarithm function to develop an equivalent formulation of the problem and obtained high quality solutions. He assumed nonconstant component hazard functions, Erlang distribution for component time-to-failure, imperfect switching, and multiple component choices for each subsystem. Cold standby redundancy strategy is considered for a nonrepairable series-parallel system. Zhao and Liu [16] proposed a stochastic model for the redundancy optimization problem aiming at maximizing system lifetime or system reliability by considering both active and standby redundancy strategies. They used stochastic simulation, Genetic Algorithm (GA), and Neural Network (NN) to develop a hybrid intelligent algorithm to solve the problem. Liang and Smith [17] proposed an Ant Colony Optimization (ACO) Algorithm to solve the redundancy allocation problem (RAP) for a series-parallel system. The objective is to maximize system reliability when active redundancy strategy has been implemented in the system. Restrictions are set on system cost and system weight, in addition. They found the ACO algorithm to be very effective and efficient for solving NP-hard reliability design problems because it brings GA flexibility, robustness, and ease of implementation along with improving its random behavior. Tavakkoli-Moghaddam et al. [18] studied RAP for a series-parallel system by considering both active and standby redundancy strategies. They formulated the problem as a nonlinear integer programming model and used a Genetic Algorithm to solve the NP-hard problem and maximize system reliability. Sadjadi and Soltani [19] proposed a heuristic and a hybrid GA for the RAP in a series-parallel system to maximize its reliability. Parameters of the proposed hybrid GA are calibrated using Taguchi’s robust design method to enhance efficiency and effectiveness of the algorithm. Solving numerical examples indicated that the proposed heuristic method is time-efficient and produces comparable solutions to the hybrid GA in terms of quality. Kumar et al. [20] studied a multiobjective multilevel RAP and proposed multiobjective hierarchical Genetic Algorithms to solve two numerical examples. They integrated the hierarchical genotype encoding scheme with two multiobjective Genetic Algorithms. Beji et al. [21] proposed a hybrid metaheuristic algorithm based on Particle Swarm Optimization (PSO) and local search for RAP in a series-parallel system and tried to maximize system reliability.A large number of studies on redundancy allocation problem have been conducted after 2010. Among those who studied multiobjective RAP (MORAP), Zio and Bazzo [22] used a level diagram analysis of Pareto solutions to assist the decision maker in selecting his/her preferred system design in terms of reliability and availability. Soylu and Ulusoy [23] applied UTADIS sorting procedure to categorize the Pareto solutions obtained by augmented epsilon constraint method into preference ordered classes. They considered maximization of the minimum system reliability along with minimization of the overall system cost and weight as objectives. Safari [24] studied a MORAP by considering system reliability and overall system cost as objectives and both active and standby strategies as candidate redundancy strategies. He used a Nondominated Sorting Genetic Algorithm (NSGA-II) to solve the multiobjective RAP. Khalili-Damghani and Amiri [25] applied three solution methods, epsilon constraint, multistart partial bound enumeration algorithm, and Data Envelopment Analysis (DEA) to optimize system reliability, cost, and weight in a RAP for a series-parallel system. Chambari et al. [26] studied a biobjective RAP trying to maximize system reliability and minimize overall cost along with making decision about using active and/or standby redundancy strategies for a system with nonrepairable components. They proposed two metaheuristics, NSGA-II and MOPSO, to solve the problem. However, Zoulfaghari et al. [27] considered system reliability and availability as objectives of a RAP with both repairable and nonrepairable components and proposed a mixed integer nonlinearprogramming (MINLP) model for the problem. Cao et al. [28] used a decomposition-based exact approach to solve a multiobjective RAP of a mixed system trying to optimize system reliability, cost, and weight. Garg and Sharma [29] studied a multiobjective RAP with nonstochastic uncertain parameters by considering system reliability and cost as objectives. They formulated the problem as a fuzzy multiobjective optimization problem.Firefly Algorithm as a new metaheuristic optimization method was introduced in 2008 by Yang [30]. dos Santos Coelho et al. introduced a modified FA approach combined with chaotic sequences to optimize reliability-redundancy problem [31]. Many other authors proposed metaheuristic algorithms for RAP. Sadjadi and Soltani [32] developed a heuristic method and a honey bee mating algorithm to solve the large-scale RAP. Hsieh and Yeh [33] applied a penalty guided bee colony algorithm for RAP in a series-parallel system. Several other metaheuristic algorithms have been proposed by other researchers in [34–38].In this paper, a novel mathematical model of a system of repairable and nonrepairable components is formulated. The model contains two objectives: firstly, it aims to select a proper redundancy strategy for nonrepairable part of the system and secondly, it offers a maintenance policy for repairable part of the system. Minimizing net present value of total cost and maximizing system reliability are objectives of the problem. In addition, different types of redundancy strategies, repair, and replacement actions are considered in order to model the problem as realistic as possible. Other practical constraints such as available budget for purchasing redundant components, volume, weight, and maximum allowed failure rate in each inspection period are taken into account. Due to NP-hardness of the problem, the authors tried to employ metaheuristic methods to solve proposed model. Three common solution approaches called NSGA-II, MOPSO, and MOFA were selected based on the Vanoye and Parra classification. Ruiz-Vanoye and Díaz-Parra [39] classified metaheuristics into three groups: metaheuristics based on gene transfer (like Genetic Algorithms), metaheuristics based on interactions among individual insects (e.g., Ant Colony, Honey Bees, and Firefly Algorithms) and metaheuristics based on biological aspects of alive beings (such as Simulated Annealing, Tabu Search, and Particle Swarm Optimization Algorithms).Remainder of the paper is organized as follows. In Section2, mathematical formulation of the problem is proposed followed by detailed explanation of objective functions and constraints. Three metaheuristic algorithms are presented in Section 3 to solve the proposed model. A set of numerical examples have been solved using the metaheuristics in Section 4. Then, obtained results are indicated and computational analysis is carried out. Finally, a summary of the paper and conclusions have been presented in Section 5.


## 2. Problem Formulation

In this section, a new integrated mathematical model is proposed for redundancy allocation and reliability-centered maintenance problems. Objective of the reliability problems could be one or a set of the following objectives: maximizing system reliability and minimizing cost, weight, and volume of the system. In this paper, system reliability and costs including maintenance and operational costs are considered as objectives.In most articles, the system under study includes either repairable or nonrepairable components. However, systems usually consist of repairable and nonrepairable components simultaneously in real world [27]. Generally, components of electronic devices are not repairable and should be replaced by new ones after failure. However, components of mechanical systems are usually repairable and repairing or replacing the broken component after failure brings the system back to the normal condition.In this paper, a system of electronic and mechanical components has been considered. Figure1 represents configuration of the system. Two approaches are applied to achieve the highest possible system reliability: (1) maximizing reliability of each component by using a diverse set of high quality and reliable redundant parallelized components (heterogeneous redundancy); (2) choosing optimal maintenance policies. It is obvious that the first approach can be used for nonrepairable electronic components and the second one is applied on repairable mechanical components.Figure 1
A system of mechanical and electronic subsystems.Assumptions. The following assumptions have been taken into account in proposing mathematical formulation of the problem under study:(i)
The system is comprised of two subsystems in series: mechanical components and electronic components (Figure1). Mechanical components are repairable while electronic components are nonrepairable.(ii)
Selecting optimal maintenance policy for mechanical components is considered in order to maximize system reliability.(iii)
Selecting a proper redundancy strategy, active or cold standby, and determining the number of redundant components in the electronic section is taken into account aimed at improving system reliability. Selecting active redundancy strategy adds operational costs to the system cost while selecting cold standby redundancy strategy threatens system performance by imperfect switching.(iv)
It is possible to use different type of components with different initial and operational costs and failure rates for the electronic subsystem.(v)
Since unstable market and economic conditions may have serious effects on results, inflation rate and time value of money are considered in computations.(vi)
Required resources such as financial resources, human resources, volume, and weight are known deterministically.(vii)
Time-to-failure distributions of components are independent.(viii)
A fixed amount of budget is available at time zero to purchase electronic components (initial cost).(ix)
Secondary cost of the system is calculated by taking operational costs of the electronic subsystem and maintenance costs of the mechanical subsystem into account during the system running period (mission time).(x)
The system mission time is finite.(xi)
Repair and replacement times and restoration times are calculated as the system downtime cost.(xii)
In cold standby strategy, redundant components do not fail before their activation. In addition, failure rates in active redundancy strategy are larger than cold standby strategy because active redundant components are exposed to the operational stresses. (i.e.,λij1≥λij2)(xiii)
Failure detection mechanism and switching are imperfect.


### 2.1. Mathematical Model

The constraints of the problem can be formulated as follows:(1)∑k=12∑j=1Si∑i=1Ewij·xijk≤W,(2)∑k=12∑j=1Si∑i=1Evij·xijk≤V,(3)∑k=12∑j=1Sixijk≤Nmaxi∀i=1,…,E,(4)xijk≤Myik∀i=1,…,E,j=1,…,Si,k=1,2,(5)∑k=12yik=1∀i=1,…,E,(6)∑l=1Lxmrlt+∑u=1Srxrrut≤1∀r=1,…,M,t=1,…,mT,(7)∑r=1M∑l=1Lhrl·xmrlt+∑u=1Srhru·xrrut≤Ho∀t=1,…,mT,(8)λrt=λrt-1+ηr1-∑l=1Lxmrlt+∑u=1Srxrrut+∑l=1Lxmrltλ-rl+∑u=1Srxrrutλ̿ru∀r=1,…,M,t=1,…,mT,(9)λrt≤λmax.r∀r=1,…,M,t=1,…,mT,(10)∑i=1E∑j=1Si∑k=12xijkcij≤B0,(11)xmrlt∈0,1,xrrut∈0,1,yik∈0,1,xijk∈Z+.Constraints (1) and (2) ensure that total weight and volume of the electronic components are lower than the maximum allowed amounts. Constraint (3) represents the maximum allowed number of redundant components in each subsystem. Constraints (4) and (5) ensure that each subsystem uses only one redundancy strategy. Constraint (6) states that components of the mechanical section can be repaired or replaced only at inspection points. Constraint (7) ensures that the number of required operators to perform maintenance actions does not exceed the number of available operators. Constraint (8) calculates failure rates of mechanical components in each period. Failure rate is increased by ηr if no repair or replacement action is taken on the component at period t and is changed into λ-rl or λ̿ru if the component is repaired or replaced at period t, respectively. Constraint (9) ensures that failure rates of mechanical components in each period do not exceed the maximum allowed amount. Constraint (10) represents the maximum available budget for purchasing electronic components at the beginning of the mission. Finally, constraint (11) represents domain restrictions of the decision variables.


### 2.2. Objective Functions

The proposed model contains two objectives:(1)
Maximizing system reliability in each period by selecting optimal redundancy strategy and maintenance policy.(2)
Minimizing secondary cost of the system including operational costs of the electronic subsystem and maintenance costs of the mechanical subsystem during the system mission.In the following of this section, detailed explanations on the introduced objectives have been proposed.


#### 2.2.1. System Reliability

The system under study consists of two electronic and mechanical subsystems connected in series according to Figure1. Thus, system reliability is equal to RS(t)=RE(t)×RM(t), where RE(t) and RM(t) are reliability of electronic and mechanical subsystem, respectively.Electronic Subsystem. The electronic part consists of E subsystems connected in series. Each subsystem can use either active or cold standby redundancy strategy to maximize reliability of the subsystem. As discussed previously, electronic components are nonrepairable and failed components can only be replaced by redundant components until all redundant components are used. Generally, cold standby redundancy leads to higher system reliability in comparison with active redundancy [15].LetH subsystems use active redundancy and C (C=E-H) subsystems use cold standby redundancy. Thus, reliability of the electronic subsystem is calculated according to the following equation:(12)REt=RHottRColdt.RHott is reliability of the subsystem under active redundancy strategy that is calculated by considering it as a subsystem with a number of redundant components and calculated as follows: (13)RHott=∏i∈H1-1-rijtxijk,where rijt is the reliability of component j in the subsystem i at period t. Also, reliability of a system under cold standby redundancy is calculated as (14)(14)RColdt=∏i∈Crijt+∑n=1xijk-1∫oTρiufijnurijT-udu,where fijnt is the density function of nth failure for type j of the subsystem i. Also, ρiu refers to switch reliability. Coit [15] calculated a lower bound for RColdt and proposed an approximation of the system reliability by considering imperfect switching using the following equation: (15)R~Coldt=∏i∈Crijt+ρit∑n=1xijk-1∫oTfijnurijT-udu.According to Coit’s approximation [26] and since most electronic devices have exponential failure distributions and system working time is equal to the sum of working times of components, time-to-failure distribution function of each component is Erlang (λ,k); thus, reliability of the electronic subsystem is calculated as follows:(16)REt=∏i∈H1-∏j=1Si1-e-λijk·t∑l=0kij-1λijktll!xijk∏i∈C1-∏j=1Si1-e-λijt∑l=0kij-1λijktll!+ρit∑l=kijkijxijk-1λijktll!.Mechanical Subsystem. The mechanical subsystem includes M devices connected in series. Therefore, all devices should be in the working condition to have an active subsystem during the mission time. In order to prevent system shutdown, scheduled maintenance actions should be taken. Reliability of the mechanical subsystem is calculated by multiplying reliability of all components as follows:(17)RMt=∏r=1Mprt,where pr(t) indicates reliability of the component r at period t. Each time period is divided into m equal intervals for inspection. Thus, we have m×T inspections during the time horizon. All components are inspected in each interval and one of the following maintenance actions is performed on each component: (1) inspection, (2) preventive repair, and (3) preventive replacement [9]. Each action has a special effect on reliability of the component with a different cost dependent on the amount of required resources. Lifetime of the mechanical devices has Weibull α,1/λ distribution function. Thus, reliability of the mechanical subsystem is calculated as follows:(18)RMt=∏r=1Me-λrttα.In order to calculate reliability of the system in each period, a few points should be considered: (1) each component deteriorates after activation with an increase in its failure rate; (2) failure rate of a component is reduced or fixed by performing maintenance actions [3].


#### 2.2.2. System Cost

Total cost of the system (CT) could be divided into two parts: initial costs and secondary costs. Initial costs (CI) include purchasing costs of the electronic devices while secondary costs (CS) refer to operational costs of the electronic subsystem (CSO) [40] plus inspection and maintenance costs of the mechanical subsystem (CSM). The objective is to minimize secondary costs (CS=CSO+CSM) according to restrictions on the initial budget for purchasing costs and other constraints.Initial Cost. Initial cost is calculated based on the purchase cost and number of the components implemented in electronic subsystem at the beginning of running the system(19)CI=∑i=1E∑j=1Si∑k=12xijkcij.Secondary Cost. As mentioned before, secondary cost is sum of the operational costs of the electronic components and maintenance and inspection costs of the mechanical components. Operational cost is calculated by multiplying the number of operating components by unit operational cost (c´ij). Number of operating components in each period is determined according to the redundancy strategy used in the system. If active redundancy is used, operational cost is calculated for all working components because all redundant components are operating from time zero in this strategy. However under cold standby redundancy, operational cost is only calculated for the operating redundant component because one component is required to be operating in this strategy.It should be noted that, in active redundancy, expected value for failure of each component in each period (λijk) is deducted from the total amount. Using compound interest rate, operational costs of different periods are converted to the present time and net present value of the total operational cost is calculated using the following equation:(20)CSO=∑t=1mT∑i=1E∑j=1Si∑k=12cij′xijk-t-1λijkyik+1-yik·1+It.Maintenance costs of the mechanical components consist of repair/replacement costs, plus system downtime cost due to performing maintenance actions. Repairing the component adds a repair cost while replacing it adds a purchasing cost to the secondary cost. System downtime cost is calculated based on the selected maintenance. Present value of maintenance costs is calculated using (21) and added to the secondary costs(21)CSM=∑t=1mT∑r=1M∑l=1Lcrl1+crl3xmrlt+∑u=1Srcru2+cru4xrrut1+It.


## 2.1. Mathematical Model

The constraints of the problem can be formulated as follows:(1)∑k=12∑j=1Si∑i=1Ewij·xijk≤W,(2)∑k=12∑j=1Si∑i=1Evij·xijk≤V,(3)∑k=12∑j=1Sixijk≤Nmaxi∀i=1,…,E,(4)xijk≤Myik∀i=1,…,E,j=1,…,Si,k=1,2,(5)∑k=12yik=1∀i=1,…,E,(6)∑l=1Lxmrlt+∑u=1Srxrrut≤1∀r=1,…,M,t=1,…,mT,(7)∑r=1M∑l=1Lhrl·xmrlt+∑u=1Srhru·xrrut≤Ho∀t=1,…,mT,(8)λrt=λrt-1+ηr1-∑l=1Lxmrlt+∑u=1Srxrrut+∑l=1Lxmrltλ-rl+∑u=1Srxrrutλ̿ru∀r=1,…,M,t=1,…,mT,(9)λrt≤λmax.r∀r=1,…,M,t=1,…,mT,(10)∑i=1E∑j=1Si∑k=12xijkcij≤B0,(11)xmrlt∈0,1,xrrut∈0,1,yik∈0,1,xijk∈Z+.Constraints (1) and (2) ensure that total weight and volume of the electronic components are lower than the maximum allowed amounts. Constraint (3) represents the maximum allowed number of redundant components in each subsystem. Constraints (4) and (5) ensure that each subsystem uses only one redundancy strategy. Constraint (6) states that components of the mechanical section can be repaired or replaced only at inspection points. Constraint (7) ensures that the number of required operators to perform maintenance actions does not exceed the number of available operators. Constraint (8) calculates failure rates of mechanical components in each period. Failure rate is increased by ηr if no repair or replacement action is taken on the component at period t and is changed into λ-rl or λ̿ru if the component is repaired or replaced at period t, respectively. Constraint (9) ensures that failure rates of mechanical components in each period do not exceed the maximum allowed amount. Constraint (10) represents the maximum available budget for purchasing electronic components at the beginning of the mission. Finally, constraint (11) represents domain restrictions of the decision variables.


## 2.2. Objective Functions

The proposed model contains two objectives:(1)
Maximizing system reliability in each period by selecting optimal redundancy strategy and maintenance policy.(2)
Minimizing secondary cost of the system including operational costs of the electronic subsystem and maintenance costs of the mechanical subsystem during the system mission.In the following of this section, detailed explanations on the introduced objectives have been proposed.


### 2.2.1. System Reliability

The system under study consists of two electronic and mechanical subsystems connected in series according to Figure1. Thus, system reliability is equal to RS(t)=RE(t)×RM(t), where RE(t) and RM(t) are reliability of electronic and mechanical subsystem, respectively.Electronic Subsystem. The electronic part consists of E subsystems connected in series. Each subsystem can use either active or cold standby redundancy strategy to maximize reliability of the subsystem. As discussed previously, electronic components are nonrepairable and failed components can only be replaced by redundant components until all redundant components are used. Generally, cold standby redundancy leads to higher system reliability in comparison with active redundancy [15].LetH subsystems use active redundancy and C (C=E-H) subsystems use cold standby redundancy. Thus, reliability of the electronic subsystem is calculated according to the following equation:(12)REt=RHottRColdt.RHott is reliability of the subsystem under active redundancy strategy that is calculated by considering it as a subsystem with a number of redundant components and calculated as follows: (13)RHott=∏i∈H1-1-rijtxijk,where rijt is the reliability of component j in the subsystem i at period t. Also, reliability of a system under cold standby redundancy is calculated as (14)(14)RColdt=∏i∈Crijt+∑n=1xijk-1∫oTρiufijnurijT-udu,where fijnt is the density function of nth failure for type j of the subsystem i. Also, ρiu refers to switch reliability. Coit [15] calculated a lower bound for RColdt and proposed an approximation of the system reliability by considering imperfect switching using the following equation: (15)R~Coldt=∏i∈Crijt+ρit∑n=1xijk-1∫oTfijnurijT-udu.According to Coit’s approximation [26] and since most electronic devices have exponential failure distributions and system working time is equal to the sum of working times of components, time-to-failure distribution function of each component is Erlang (λ,k); thus, reliability of the electronic subsystem is calculated as follows:(16)REt=∏i∈H1-∏j=1Si1-e-λijk·t∑l=0kij-1λijktll!xijk∏i∈C1-∏j=1Si1-e-λijt∑l=0kij-1λijktll!+ρit∑l=kijkijxijk-1λijktll!.Mechanical Subsystem. The mechanical subsystem includes M devices connected in series. Therefore, all devices should be in the working condition to have an active subsystem during the mission time. In order to prevent system shutdown, scheduled maintenance actions should be taken. Reliability of the mechanical subsystem is calculated by multiplying reliability of all components as follows:(17)RMt=∏r=1Mprt,where pr(t) indicates reliability of the component r at period t. Each time period is divided into m equal intervals for inspection. Thus, we have m×T inspections during the time horizon. All components are inspected in each interval and one of the following maintenance actions is performed on each component: (1) inspection, (2) preventive repair, and (3) preventive replacement [9]. Each action has a special effect on reliability of the component with a different cost dependent on the amount of required resources. Lifetime of the mechanical devices has Weibull α,1/λ distribution function. Thus, reliability of the mechanical subsystem is calculated as follows:(18)RMt=∏r=1Me-λrttα.In order to calculate reliability of the system in each period, a few points should be considered: (1) each component deteriorates after activation with an increase in its failure rate; (2) failure rate of a component is reduced or fixed by performing maintenance actions [3].


### 2.2.2. System Cost

Total cost of the system (CT) could be divided into two parts: initial costs and secondary costs. Initial costs (CI) include purchasing costs of the electronic devices while secondary costs (CS) refer to operational costs of the electronic subsystem (CSO) [40] plus inspection and maintenance costs of the mechanical subsystem (CSM). The objective is to minimize secondary costs (CS=CSO+CSM) according to restrictions on the initial budget for purchasing costs and other constraints.Initial Cost. Initial cost is calculated based on the purchase cost and number of the components implemented in electronic subsystem at the beginning of running the system(19)CI=∑i=1E∑j=1Si∑k=12xijkcij.Secondary Cost. As mentioned before, secondary cost is sum of the operational costs of the electronic components and maintenance and inspection costs of the mechanical components. Operational cost is calculated by multiplying the number of operating components by unit operational cost (c´ij). Number of operating components in each period is determined according to the redundancy strategy used in the system. If active redundancy is used, operational cost is calculated for all working components because all redundant components are operating from time zero in this strategy. However under cold standby redundancy, operational cost is only calculated for the operating redundant component because one component is required to be operating in this strategy.It should be noted that, in active redundancy, expected value for failure of each component in each period (λijk) is deducted from the total amount. Using compound interest rate, operational costs of different periods are converted to the present time and net present value of the total operational cost is calculated using the following equation:(20)CSO=∑t=1mT∑i=1E∑j=1Si∑k=12cij′xijk-t-1λijkyik+1-yik·1+It.Maintenance costs of the mechanical components consist of repair/replacement costs, plus system downtime cost due to performing maintenance actions. Repairing the component adds a repair cost while replacing it adds a purchasing cost to the secondary cost. System downtime cost is calculated based on the selected maintenance. Present value of maintenance costs is calculated using (21) and added to the secondary costs(21)CSM=∑t=1mT∑r=1M∑l=1Lcrl1+crl3xmrlt+∑u=1Srcru2+cru4xrrut1+It.


## 2.2.1. System Reliability

The system under study consists of two electronic and mechanical subsystems connected in series according to Figure1. Thus, system reliability is equal to RS(t)=RE(t)×RM(t), where RE(t) and RM(t) are reliability of electronic and mechanical subsystem, respectively.Electronic Subsystem. The electronic part consists of E subsystems connected in series. Each subsystem can use either active or cold standby redundancy strategy to maximize reliability of the subsystem. As discussed previously, electronic components are nonrepairable and failed components can only be replaced by redundant components until all redundant components are used. Generally, cold standby redundancy leads to higher system reliability in comparison with active redundancy [15].LetH subsystems use active redundancy and C (C=E-H) subsystems use cold standby redundancy. Thus, reliability of the electronic subsystem is calculated according to the following equation:(12)REt=RHottRColdt.RHott is reliability of the subsystem under active redundancy strategy that is calculated by considering it as a subsystem with a number of redundant components and calculated as follows: (13)RHott=∏i∈H1-1-rijtxijk,where rijt is the reliability of component j in the subsystem i at period t. Also, reliability of a system under cold standby redundancy is calculated as (14)(14)RColdt=∏i∈Crijt+∑n=1xijk-1∫oTρiufijnurijT-udu,where fijnt is the density function of nth failure for type j of the subsystem i. Also, ρiu refers to switch reliability. Coit [15] calculated a lower bound for RColdt and proposed an approximation of the system reliability by considering imperfect switching using the following equation: (15)R~Coldt=∏i∈Crijt+ρit∑n=1xijk-1∫oTfijnurijT-udu.According to Coit’s approximation [26] and since most electronic devices have exponential failure distributions and system working time is equal to the sum of working times of components, time-to-failure distribution function of each component is Erlang (λ,k); thus, reliability of the electronic subsystem is calculated as follows:(16)REt=∏i∈H1-∏j=1Si1-e-λijk·t∑l=0kij-1λijktll!xijk∏i∈C1-∏j=1Si1-e-λijt∑l=0kij-1λijktll!+ρit∑l=kijkijxijk-1λijktll!.Mechanical Subsystem. The mechanical subsystem includes M devices connected in series. Therefore, all devices should be in the working condition to have an active subsystem during the mission time. In order to prevent system shutdown, scheduled maintenance actions should be taken. Reliability of the mechanical subsystem is calculated by multiplying reliability of all components as follows:(17)RMt=∏r=1Mprt,where pr(t) indicates reliability of the component r at period t. Each time period is divided into m equal intervals for inspection. Thus, we have m×T inspections during the time horizon. All components are inspected in each interval and one of the following maintenance actions is performed on each component: (1) inspection, (2) preventive repair, and (3) preventive replacement [9]. Each action has a special effect on reliability of the component with a different cost dependent on the amount of required resources. Lifetime of the mechanical devices has Weibull α,1/λ distribution function. Thus, reliability of the mechanical subsystem is calculated as follows:(18)RMt=∏r=1Me-λrttα.In order to calculate reliability of the system in each period, a few points should be considered: (1) each component deteriorates after activation with an increase in its failure rate; (2) failure rate of a component is reduced or fixed by performing maintenance actions [3].


## 2.2.2. System Cost

Total cost of the system (CT) could be divided into two parts: initial costs and secondary costs. Initial costs (CI) include purchasing costs of the electronic devices while secondary costs (CS) refer to operational costs of the electronic subsystem (CSO) [40] plus inspection and maintenance costs of the mechanical subsystem (CSM). The objective is to minimize secondary costs (CS=CSO+CSM) according to restrictions on the initial budget for purchasing costs and other constraints.Initial Cost. Initial cost is calculated based on the purchase cost and number of the components implemented in electronic subsystem at the beginning of running the system(19)CI=∑i=1E∑j=1Si∑k=12xijkcij.Secondary Cost. As mentioned before, secondary cost is sum of the operational costs of the electronic components and maintenance and inspection costs of the mechanical components. Operational cost is calculated by multiplying the number of operating components by unit operational cost (c´ij). Number of operating components in each period is determined according to the redundancy strategy used in the system. If active redundancy is used, operational cost is calculated for all working components because all redundant components are operating from time zero in this strategy. However under cold standby redundancy, operational cost is only calculated for the operating redundant component because one component is required to be operating in this strategy.It should be noted that, in active redundancy, expected value for failure of each component in each period (λijk) is deducted from the total amount. Using compound interest rate, operational costs of different periods are converted to the present time and net present value of the total operational cost is calculated using the following equation:(20)CSO=∑t=1mT∑i=1E∑j=1Si∑k=12cij′xijk-t-1λijkyik+1-yik·1+It.Maintenance costs of the mechanical components consist of repair/replacement costs, plus system downtime cost due to performing maintenance actions. Repairing the component adds a repair cost while replacing it adds a purchasing cost to the secondary cost. System downtime cost is calculated based on the selected maintenance. Present value of maintenance costs is calculated using (21) and added to the secondary costs(21)CSM=∑t=1mT∑r=1M∑l=1Lcrl1+crl3xmrlt+∑u=1Srcru2+cru4xrrut1+It.


## 3. Solution Approaches

As mentioned in Introduction and Literature Review, the solution methods were selected based on the Ruiz-Vanoye and Díaz-Parra [39] metaheuristics classification. In this section, the mechanisms of three algorithms, NSGA-II, MOPSO, and MOFA, are explained; then in Section 4, some numerical examples based on the proposed model are generated and are solved using these methods. Finally, obtained results by algorithms have been compared and analyzed.The proposed model for the problem contains two conflicting objectives. We try to make a tradeoff between these objectives to achieve a desired level of optimality for each objective. One of the common approaches to solve multiobjective problems is the weighted-sum method that converts the problem into a single objective problem by making a weighted linear combination (WLC) of objectives. Although it is a very popular method due to its simplicity and ease of implementation, it has some major disadvantages such as determining weight of each objective and lack of information about it. Thus, another method called Pareto set has been developed. Pareto set method produces set of solutions within the feasible region of the problem that dominate other feasible solutions. The nondominated solution sets are called Pareto optimal solutions and other inferior solutions are called dominated solutions. The decision maker selects the final set of Pareto optimal solutions according to his/her preferences and considered criteria. In this paper, three metaheuristic algorithms have been employed to produce Pareto optimal solutions.


### 3.1. Nondominated Sorting Genetic Algorithm (NSGA-II)

The NSGA-II is developed for solving multiobjective problems by adding two operators to the classic Genetic Algorithm (GA) to find Pareto optimal sets instead of finding the unique optimal solution [41]. The additional operators are (1) ranking operator which assigns a rank to each member of the generation based on nondominated sorting and (2) diversity operator which increases diversity of the produced solutions with equal ranks. In the following of this section, procedure of the proposed NSGA-II has been presented.In our study, the proposed chromosome has two parts. The first part represents the electronic section and the second part represents the mechanical section. The electronic section containsE rows, which refer to the number of electronic subsystems, and Si+1 columns where Si columns are allocated to different types of components and the last column represents the selected redundancy strategy. If active redundancy is selected, the related element in the last column is equal to 1 while selecting cold standby redundancy makes this element equal to 2. Element (i,j) indicates the number of implemented redundant components for component type j of subsystem i. In other words, if element (2,3) is equal to 1 it means that one redundant component is implemented for component type 3 in subsystem 2 (Table 1).Table 1
Electronic section matrix.
i
j
1
2
⋯
S
i
Redundancy strategy
1
2
0
⋯
1
1
2
1
1
⋯
2
2
⋮
⋮
⋮
⋯
⋮
⋮
E
0
2
⋯
1
1The proposed chromosome for the mechanical section is indicated by two matrixes withM rows representing mechanical subsystems and mT columns representing inspection intervals. Elements of the first matrix take three values: 0, 1, or 2. Number 0 indicates that no repair or replacement action is required for that subsystem. Number 1 refers to replacement and number 2 refers to repair action. The second matrix determines type of the replacement and repair actions.Crossover Operator. The proposed matrix for the electronic section is divided into two parts in order to perform crossover. The first part is related to selecting a redundancy strategy for each subsystem, that is, the last column, while the second part is related to determining the number of redundant components for components of each subsystem, that is, (i,j) elements. The chromosome is cut from a randomly selected point to diversify the selected redundancy strategies. Then, the resulting two slices from parent chromosomes are interchanged. This process changes the selected redundancy strategy for some subsystems. For the second part of the chromosome, a vertical cut is applied from a randomly selected point. Then, the resulting two slices are interchanged. Figure 2 indicates the applied crossover operation.Figure 2
Crossover in Genetic Algorithm.Mutation Operator. Mutation process for the first chromosome, that is, electronic section, is performed as explained in the following. First, a subsystem is selected randomly. Next, one of the allowed components is selected for the selected subsystem. Then, a random number between 1 and Nmaxi is assigned to the selected subsystem and component. Finally, redundancy strategy is selected randomly. For the mechanical section, mutation process is started by selecting a subsystem and a period randomly. Then, selected maintenance policy and type of repair/replacement actions are changed.


### 3.2. Particle Swarm Algorithm

The Multiobjective Particle Swarm Optimization (MOPSO) is a metaheuristic algorithm capable of producing high quality nondominated Pareto optimal solutions with high diversities for multiobjective problems. The MOPSO is widely used by researchers due to its simplicity and successful performance in continuous optimization problems. The idea of this algorithm is inspired by a swarm of birds looking for food [42]. In this algorithm, each feasible solution is indicated as a particle with known velocity and fitness value. Particles move in the search space and obtained results are classified based on a fitness criterion at the end of each time interval. Particles are changed into particles with higher fitness values gradually based on the following model:(22)V→it+1=W×V→it+C1r1P→i-X→it+C2r2P→g-X→iti=1,2,…,m,X→it+1=X→it+V→it+1i=1,2,…,m,where W, C1, C2, and t represent inertia weight, cognitive factor, social factor, and iteration number, respectively. Also r1 and r2 are random numbers in [0,1] and V→i and X→i are velocity and position vector of ith particles, respectively. Multiobjective approach of this algorithm (MOPSO) was proposed by Coello Coello and Lechuga in 2002 [43]. The improved version of this algorithm, in which constraint-handling mechanism and a mutation operator have been considered, was presented in 2004 [44].In MOPSO, all objective functions are calculated and evaluated for each particle and the nondominated solutions (based on Pareto set concept) found by the particles are stored in a repository. The size of repository is limited and is set by decision maker. In addition, the search space is divided into hypercubes with a fitness value. Fitness value is inversely proportional to the number of particles it contains [45]. Then, the selection method is used to choose a hypercube and the best nondominated position (called leader) and finally, the velocities and positions of the particles are updated. This process is repeated for a certain number of iterations.


### 3.3. Firefly Algorithm

Firefly Algorithm (FA) is the last method applied in this study. It was introduced by Yang in 2010 [30] as a new approach for solving complex and continuous problems. The algorithm is inspired by the flashing behavior of fireflies to attract each other. Considering the following three rules, Firefly Algorithm introduced by Yang:(1)
All fireflies are unisex and attracting a firefly by another one is regardless of its sex.(2)
Attractiveness is proportional to the brightness and both of these features will decrease with increasing distance. Less bright fireflies are always attracted to the brighter one and if there is no brighter one, the move will be randomly.(3)
The brightness of fireflies is defined according to the objective function (like fitness function in Genetic Algorithm).Variation of light intensity and formulation of the attractiveness of the FA are two important issues. Brightness of each firefly at place x is defined as I(x)∝f(x) and attractiveness β(r) is defined with respect to distance of the firefly i from firefly j that is calculated as following equation:(23)βr=β0e-γrmm≥1,where β0 is attractiveness at r=0 and γ is fixed light absorption coefficient. Based on the Cartesian distance, the distance between the firefly i and firefly j at xi and xj is obtained as rij=xi-xj=∑k=1dxi,k-xj,k2 where xi,k is the kth component of xi in spatial coordinate. Calculating distances is not limited only to the Euclidean coordinate, but given the type of the problem it can be defined differently, for example, as time interval. Moving of firefly i towards more attractive firefly j is calculated according to (24)xi+1=xi+β0e-γrij2xj-xi+α0ϵi.The second term of (24) is movement because of attraction and the third term is the random movement where α0 is a randomization parameter and ϵi is vector of random numbers with a Gaussian or uniform distribution. Although the basic design of this method was for continuous problems, recent studies have shown that this approach is also very efficient in discrete problems [46].Single objective form of Firefly Algorithm (FA) was developed to Multiobjective Firefly Algorithm (MOFA) in 2013 by Yang [47]. In this approach, nondominated solutions are detected based on the objective function value. The main steps of the MOFA can be summarized as shown in Algorithm 1.Algorithm 1
Define objective functions:f1x,f2x,…,fkx where x=(x1,x2,…,xd)T
Generate initial population of firefliesxi=(i=1,…,n)
Formulate light intensityI
Define absorption coefficientγ; randomization parameter α0,ϵi; vector of random numbers and maximum iteration
Whilet < Maximum Iteration
fori,j=1 : n (all n fireflies)
Evaluate approximationsPFi and PFj to the Pareto Front
ifPFj dominates PFi
Move fireflyi towards j using (24)
if pervious position doesn’t dominate new one
New position replaced with old one
end if
end if
end for
Update and pass non dominated solution to next generation
Updatet=t+1
End whileSince that, variables of our problem are binary and integer and MOPSO algorithm and MOFA find solutions in continuous space; round function is used to convert real number to the integer, and to change detected solutions to the binary style, sigmoid function is applied.


## 3.1. Nondominated Sorting Genetic Algorithm (NSGA-II)

The NSGA-II is developed for solving multiobjective problems by adding two operators to the classic Genetic Algorithm (GA) to find Pareto optimal sets instead of finding the unique optimal solution [41]. The additional operators are (1) ranking operator which assigns a rank to each member of the generation based on nondominated sorting and (2) diversity operator which increases diversity of the produced solutions with equal ranks. In the following of this section, procedure of the proposed NSGA-II has been presented.In our study, the proposed chromosome has two parts. The first part represents the electronic section and the second part represents the mechanical section. The electronic section containsE rows, which refer to the number of electronic subsystems, and Si+1 columns where Si columns are allocated to different types of components and the last column represents the selected redundancy strategy. If active redundancy is selected, the related element in the last column is equal to 1 while selecting cold standby redundancy makes this element equal to 2. Element (i,j) indicates the number of implemented redundant components for component type j of subsystem i. In other words, if element (2,3) is equal to 1 it means that one redundant component is implemented for component type 3 in subsystem 2 (Table 1).Table 1
Electronic section matrix.
i
j
1
2
⋯
S
i
Redundancy strategy
1
2
0
⋯
1
1
2
1
1
⋯
2
2
⋮
⋮
⋮
⋯
⋮
⋮
E
0
2
⋯
1
1The proposed chromosome for the mechanical section is indicated by two matrixes withM rows representing mechanical subsystems and mT columns representing inspection intervals. Elements of the first matrix take three values: 0, 1, or 2. Number 0 indicates that no repair or replacement action is required for that subsystem. Number 1 refers to replacement and number 2 refers to repair action. The second matrix determines type of the replacement and repair actions.Crossover Operator. The proposed matrix for the electronic section is divided into two parts in order to perform crossover. The first part is related to selecting a redundancy strategy for each subsystem, that is, the last column, while the second part is related to determining the number of redundant components for components of each subsystem, that is, (i,j) elements. The chromosome is cut from a randomly selected point to diversify the selected redundancy strategies. Then, the resulting two slices from parent chromosomes are interchanged. This process changes the selected redundancy strategy for some subsystems. For the second part of the chromosome, a vertical cut is applied from a randomly selected point. Then, the resulting two slices are interchanged. Figure 2 indicates the applied crossover operation.Figure 2
Crossover in Genetic Algorithm.Mutation Operator. Mutation process for the first chromosome, that is, electronic section, is performed as explained in the following. First, a subsystem is selected randomly. Next, one of the allowed components is selected for the selected subsystem. Then, a random number between 1 and Nmaxi is assigned to the selected subsystem and component. Finally, redundancy strategy is selected randomly. For the mechanical section, mutation process is started by selecting a subsystem and a period randomly. Then, selected maintenance policy and type of repair/replacement actions are changed.


## 3.2. Particle Swarm Algorithm

The Multiobjective Particle Swarm Optimization (MOPSO) is a metaheuristic algorithm capable of producing high quality nondominated Pareto optimal solutions with high diversities for multiobjective problems. The MOPSO is widely used by researchers due to its simplicity and successful performance in continuous optimization problems. The idea of this algorithm is inspired by a swarm of birds looking for food [42]. In this algorithm, each feasible solution is indicated as a particle with known velocity and fitness value. Particles move in the search space and obtained results are classified based on a fitness criterion at the end of each time interval. Particles are changed into particles with higher fitness values gradually based on the following model:(22)V→it+1=W×V→it+C1r1P→i-X→it+C2r2P→g-X→iti=1,2,…,m,X→it+1=X→it+V→it+1i=1,2,…,m,where W, C1, C2, and t represent inertia weight, cognitive factor, social factor, and iteration number, respectively. Also r1 and r2 are random numbers in [0,1] and V→i and X→i are velocity and position vector of ith particles, respectively. Multiobjective approach of this algorithm (MOPSO) was proposed by Coello Coello and Lechuga in 2002 [43]. The improved version of this algorithm, in which constraint-handling mechanism and a mutation operator have been considered, was presented in 2004 [44].In MOPSO, all objective functions are calculated and evaluated for each particle and the nondominated solutions (based on Pareto set concept) found by the particles are stored in a repository. The size of repository is limited and is set by decision maker. In addition, the search space is divided into hypercubes with a fitness value. Fitness value is inversely proportional to the number of particles it contains [45]. Then, the selection method is used to choose a hypercube and the best nondominated position (called leader) and finally, the velocities and positions of the particles are updated. This process is repeated for a certain number of iterations.


## 3.3. Firefly Algorithm

Firefly Algorithm (FA) is the last method applied in this study. It was introduced by Yang in 2010 [30] as a new approach for solving complex and continuous problems. The algorithm is inspired by the flashing behavior of fireflies to attract each other. Considering the following three rules, Firefly Algorithm introduced by Yang:(1)
All fireflies are unisex and attracting a firefly by another one is regardless of its sex.(2)
Attractiveness is proportional to the brightness and both of these features will decrease with increasing distance. Less bright fireflies are always attracted to the brighter one and if there is no brighter one, the move will be randomly.(3)
The brightness of fireflies is defined according to the objective function (like fitness function in Genetic Algorithm).Variation of light intensity and formulation of the attractiveness of the FA are two important issues. Brightness of each firefly at place x is defined as I(x)∝f(x) and attractiveness β(r) is defined with respect to distance of the firefly i from firefly j that is calculated as following equation:(23)βr=β0e-γrmm≥1,where β0 is attractiveness at r=0 and γ is fixed light absorption coefficient. Based on the Cartesian distance, the distance between the firefly i and firefly j at xi and xj is obtained as rij=xi-xj=∑k=1dxi,k-xj,k2 where xi,k is the kth component of xi in spatial coordinate. Calculating distances is not limited only to the Euclidean coordinate, but given the type of the problem it can be defined differently, for example, as time interval. Moving of firefly i towards more attractive firefly j is calculated according to (24)xi+1=xi+β0e-γrij2xj-xi+α0ϵi.The second term of (24) is movement because of attraction and the third term is the random movement where α0 is a randomization parameter and ϵi is vector of random numbers with a Gaussian or uniform distribution. Although the basic design of this method was for continuous problems, recent studies have shown that this approach is also very efficient in discrete problems [46].Single objective form of Firefly Algorithm (FA) was developed to Multiobjective Firefly Algorithm (MOFA) in 2013 by Yang [47]. In this approach, nondominated solutions are detected based on the objective function value. The main steps of the MOFA can be summarized as shown in Algorithm 1.Algorithm 1
Define objective functions:f1x,f2x,…,fkx where x=(x1,x2,…,xd)T
Generate initial population of firefliesxi=(i=1,…,n)
Formulate light intensityI
Define absorption coefficientγ; randomization parameter α0,ϵi; vector of random numbers and maximum iteration
Whilet < Maximum Iteration
fori,j=1 : n (all n fireflies)
Evaluate approximationsPFi and PFj to the Pareto Front
ifPFj dominates PFi
Move fireflyi towards j using (24)
if pervious position doesn’t dominate new one
New position replaced with old one
end if
end if
end for
Update and pass non dominated solution to next generation
Updatet=t+1
End whileSince that, variables of our problem are binary and integer and MOPSO algorithm and MOFA find solutions in continuous space; round function is used to convert real number to the integer, and to change detected solutions to the binary style, sigmoid function is applied.


## 4. Numerical Examples

In this section, assigning proper values to the parameters of algorithms and using comparison metrics for evaluating solution methods are discussed first. Then, three different sets of test problems (small, medium, and large size) are tackled and solved using the chosen solution methods. Finally, obtained results by each algorithm have been compared and the obtained results for an example are explained.


### 4.1. Setting Parameters of the Algorithms

Setting proper values for the control parameters of metaheuristic algorithms has a significant effect on their desirable performance. Well-tuned parameters empower the algorithms in producing better solutions within shorter computation times. Thus, setting proper values for control parameters is a critical task [48]. A few additional techniques are applied to tune parameters. One of the common methods is response surface methodology (RSM). RSM is a mathematical and statistical technique that examines the relationship between one or more response variables and the set of parameters (input variables) influencing them. Using this method, the levels of parameters that optimize the response variables are identified. In the first step, main parameters and range of them are determined. In this study, we used related literature to identify these parameters and their range. Table 2 shows range of the main parameters of each algorithm.Table 2
Range of the main parameters.
NSGA-II
MOPSO
MOFA
Parameter
Range
Parameter
Range
Parameter
Range
Pop. size (Npop)
20–100
Pop. size (Npop)
20–100
Pop. size (Npop)
10–100
Max iteration number (Maxit)
50–200
Max iteration number (Maxit)
50–200
Max iteration number (Maxit)
50–200
Cross rate (Cr)
0.5–0.9
Inertia weight (W)
0.4–0.9
Randomization parameter(∝0)
0.1–0.9
Mutation rate (Mr)
0.01–0.3
Cognitive factor(C1)
1-2
Fixed light absorption coefficient(γ)
1–3
Social factor(C2)
1-2In the next step, the response variables should be determined. Three performance metrics are chosen as response variables which are CPU time, number of nondominated solutions (NNS), and Diversification Metric (DM). These metrics were selected based on the two features, the convergence speed and diversity of the detected solutions. More details on the performance metrics can be found in Section4.2.Central composite design (CCD) with 6 center points is applied for the experiments. Experiments are run by Design Expert 9. According to the number of input variables and type of the design, different number of experiments should be run. For instance, in case of the four parameters and 6 center points’ design, 46 experiments are required. After performing the experiments, analysis of variance (ANOVA) is applied to fit an adequate model to the experimental data. Last step is setting goals for responses to generate optimal condition (optimal level of the parameters). Here, we aim to minimize CPU time and maximize NNS and DM.Figure3 shows counter plot of RSM results for MOFA approach. Desirability displays the amount of goals that have been met. It can be concluded from Figure 3 that the highest desirability is obtained when Npop level is medium and Maxit is high. RSM is applied for all three methods and the optimum values for the algorithms parameters are presented in Table 3.Table 3
Optimum values of the algorithms parameters.
NSGA-II
MOPSO
MOFA
Pop. size
100
Pop. size
100
Pop. size
25
Max iteration number
50
Max iteration number
50
Max iteration number
200
Cross rate
0.9
W
0.9
∝
0
0.1
Mutation  rate
0.3
C
1
1
γ
1
C
2
2Figure 3
Counter plot for MOFA: desirability versus Maxit andNpop.Other parameters are set according to the literatures as follows:MOPSO:r1,r2=rand[0,1], repository size = 500 and MOFA: ϵi=[rand-1/2], β0=1.


### 4.2. Performance Measures

According to Deb et al. [49], there are two main features that must be taken into account in order to evaluate performance of metaheuristic algorithms in solving multiobjective optimization problems: (1) convergence to the Pareto optimal set and (2) diversity of the produced Pareto optimal set. Numerous metrics have been proposed for quantifying these features in the literature. Yu and Gen [50] and Zitzler and Thiele [51] proposed some criteria such as number of nondominated solutions (NNS), Error Ratio (ER), and Generational Distance (GD) for measuring accuracy of algorithms. Higher values for NNS and lower values for GD and ER are preferred and show that the algorithm under analysis performed better in providing a set of Pareto optimal solutions. In addition, Diversification Metric (DM), Spacing Metric (SM), and Maximum Spread Metrics (MS) have been proposed to check ability of the algorithms in providing a diverse set of Pareto optimal solutions. In this study we utilized DM and MS to compare algorithms performance which are calculated using(25)DM=∑i=1Nmax⁡xi-yi1/2,MS=∑m=1Mmaxi=1N⁡fmi-mini=1N⁡fmi21/2,where xi-yi is the Euclidean distance between of the nondominated solution xi and the nondominated solution yi and N denote number of the Pareto solutions [52]. The MS presents the distance between the boundary solutions in the obtained nondominated solutions (NDS) where fmi is the mth objective function value for ith Pareto solution.


### 4.3. Numerical Examples

In order to examine performance of the algorithms, three sets of test problems with different sizes (small, medium, and large) are simulated. The dimensions of these sets are shown in Table4. The properties of the test problems were generated using random functions. All three algorithms are programmed in MATLAB and the test problems have been solved on a PC with 4 GB RAM/1.80 GHz CPU. Using the algorithm’s parameters in Table 3, twenty different problems in each set are solved by NSGA-II, MOPSO, and MOFA. All problems are run for five times and the mean of the obtained values for the metrics is considered for algorithm’s performance comparison. In each run, 500 solutions are generated (according to the number of population and maximum iterations) that included infeasible, dominated, and similar solutions. Detected solutions by the algorithms are checked for feasibility firstly and if all the constraints are satisfied, the generated solution is accepted; otherwise another solution would be found. This process has been coded in MATLAB. Mean and standard deviation values of CPU time, NNS, DM, and MS for each case are presented in Table 5.Table 4
Test problems dimensions.
Case
Number of electronic components;i
Number of mechanical components;r
Number of component types in the electronic subsystem;j
Number of component types in the mechanical subsystem;u
Small size
∗
rand
[
3
,
8
]
rand
[
3
,
7
]
rand
[
2,4
]
rand
[
2,4
]
Medium size
rand
[
8
,
13
]
rand
[
7,11
]
rand
[
2,4
]
rand
[
2,4
]
Large size
rand
[
13
,
18
]
rand
[
11,15
]
rand
[
2,4
]
rand
[
2,4
]
∗This function returns an integer number between lower and upper bounds randomly.Table 5
Mean and standard deviation of the metrics for the different sizes of the test problems.
Problem size
Time (seconds)
NNS
DM
MS
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Small
1.91(0.26)
0.78 (0.23)
3.50(1.16)
22.9 (13.54)
12.3(3.45)
3.3(1.26)
55.35(30.22)
60.69 (20.54)
22.60(11.72)
0.01493(0.0257)
0.07999 (0.0542)
0.04230(0.0376)
Medium
2.59(0.21)
1.45 (0.46)
6.99(2.05)
12.85 (5.60)
11.90(3.80)
3.85(1.14)
39.02(17.32)
65.79 (23.50)
28.65(12.71)
0.0238(0.0367)
0.1562 (0.2004)
0.0482(0.0303)
Large
31.49(6.91)
22.35 (9.36)
91.25(25.48)
15.60 (6.22)
12.84(4.19)
4.11(1.49)
51.35(25.77)
97.44 (30.93)
34.35(13.81)
0.00419(0.0048)
0.14625 (0.0285)
0.07408(0.0418)
∗Numbers in the parentheses are standard deviation.Four metrics are considered to evaluate algorithms performance (Table5). These are CPU time, NNS, DM, and MS. The lower values for time and the higher values for the rest of the metrics are desirable. The first rank of the CPU time metric in all cases (small, medium, and large sized instances) belongs to MOPSO approach. It means that MOPSO is more time-efficient. Second rank belongs to NSGA-II and MOFA represents a very long CPU time in comparison with the other competing methods.It also can be concluded from Table5 that the average of NNS in the NSGA-II method is higher than other methods. This means that the NSGA-II can generate more nondominated solutions on average. However, MOPSO has provided close values of this metric to NSGA-II in medium and large sized cases.Obtained values for the DM and MS indicate higher diversity solutions of the produced nondominated solutions by MOPSO in all three cases. This means that the MOPSO method has a wider spread.Our model has two conflicting objectives, reliability and cost. As the reliability of the system goes higher its cost increases too. We try to find an optimal tradeoff between these objectives. To compare quality of solutions generated by three approaches, mean and lower and upper bounds of objectives are reported in Tables6 and 7. It should be noted that the average values of the parameters had been calculated in each case.Table 6
Boundary values for reliability.
Method
Small size problems
Medium size problems
Large size problems
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
0.991548
0.945431
0.96181
0.985369
0.914826
0.938127
0.99101
0.86809
0.90353
Lower bound
0.981583
0.908162
0.939269
0.970327
0.860268
0.909683
0.98771
0.80074
0.86255
Upper bound
0.996518
0.988152
0.980975
0.994157
0.974125
0.957834
0.99190
0.94698
0.93663Table 7
Boundary values for cost.
Method
Small size problems
Medium size problems
Large size problems
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
1713.687
1734.326
1870.389
2944.94
3035.20
3189.41
4725.122
4900.882
5166.778
Lower bound
1652.382
1632.584
1790.569
2891.92
2883.30
3070.47
4662.533
4586.47
5015.414
Upper bound
1848.888
2014.234
1987.349
3065.41
3356.64
3346.01
4889.001
5529.174
5369.501As it can be concluded from Tables6 and 7, MOPSO method provides wider ranges of the objective functions but the quality of the solutions which are detected by NSGA-II is better; in other words, solutions with higher system reliability and lower values of cost were generated by NSGA-II.In remainder of this section, an example is considered based on the proposed model and the results are explained. Tables11 and 12 provide data of failure rate, Weibull distribution parameters, repair, replacement, and downtime costs, and required human resources (number of operators) for each maintenance action in the mechanical section. The electronic section includes 14 subsystems connected in series in which each subsystem can have 3 or 4 types of components. Also, the mechanical section includes 11 repairable and replaceable subsystems. Table 13 presents upper bounds for some problem parameters.This example is solved using three approaches for ten times and the mean and standard deviation values for the performance metrics and objective functions boundary are reported in Tables8 and 9.Table 8
Mean and standard deviation values for performance metrics.
Algorithm
NSGA-II
MOPSO
MOFA
Mean
Std Dev
Mean
Std Dev
Mean
Std Dev
Time
2.74
0.12
1.49
0.06
6.84
0.49
NNS
20.40
5.19
13.20
3.65
3.70
1.06
DM
44.96
10.79
91.46
15.11
30.17
12.21
MS
0.0068
0.0165
0.1436
0.0132
0.0508
0.0218Table 9
Objectives boundary values for the sample.
Method
Reliability
Cost
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
0.999421
0.87153
0.90719
3848.62
3900.48
4285.62
Lower bound
0.999825
0.802065
0.841807
3760.333
3485.712
4044.209
Upper bound
0.999833
0.977786
0.951062
4027.349
4488.861
4585.62As concluded before, MOPSO method provides better values for the metrics (time, NNS, DM, and MS) and wider boundary values for the objectives but as Table9 shows, NSGA-II approach detects Pareto solutions lead to higher system reliability and lower cost.The proposed model in this paper aims at finding proper maintenance policies and effective redundancy strategies. Table10 represents one of the nondominated solutions that are found by NSGA-II for the mentioned example. The selected maintenance policies and redundancy strategies including number and type of the redundant components and the implemented redundancy strategy in the electronic section along with the selected repair and replacement actions in the mechanical section are shown in Table 10. The first four columns in the electronic section represent the number of redundant components and the last column indicates the implemented redundancy strategy in each subsystem, where letter A means active strategy and letter C indicates cold strategy. In the mechanical subsystem, Xr and Xm indicate replacement and repair action in each period, respectively. Xr and Xm are followed by a number which determines the type of replacement/repair action. For instance, Xr2 means replacement type two is required and element 00 indicates that no repair or replacement action has been taken in that period. Figure 4 shows Pareto solutions obtained by NSGA-II and MOPSO methods.Table 10
Selected maintenance actions and redundancy strategies by the NSGA-II.
Electronic section
Mechanical section
i
j
r
t
1
2
3
4
Strategy
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
0
0
1
0
A
1
00
00
00
X
r
1
00
X
r
2
00
X
r
2
00
00
X
r
2
X
m
1
00
X
r
1
00
2
1
1
0
0
C
2
00
00
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
X
r
2
X
r
1
X
m
2
00
3
1
0
0
0
A
3
00
X
r
1
00
X
r
1
00
00
X
r
1
X
r
2
00
00
X
r
2
X
r
2
X
m
2
X
r
2
00
4
0
1
0
0
C
4
00
00
00
X
r
2
00
00
X
m
2
00
00
X
m
2
00
X
r
2
X
r
2
X
m
2
00
5
0
1
0
0
A
5
00
X
r
1
00
00
X
r
2
00
00
X
r
1
00
00
X
r
1
00
X
r
2
X
r
1
00
6
0
1
0
1
A
6
00
X
r
2
00
X
r
2
X
r
2
00
00
X
r
2
00
00
X
r
1
00
X
m
2
X
m
2
00
7
0
0
1
0
A
7
00
X
r
2
00
00
X
m
2
00
00
X
r
2
00
00
X
r
1
00
X
r
2
X
m
2
X
r
2
8
0
1
0
0
A
8
00
X
r
2
00
00
X
r
2
X
r
2
00
00
X
m
2
00
00
X
r
2
00
X
r
2
00
9
1
0
1
0
A
9
00
X
r
2
X
m
2
X
r
2
00
00
X
r
2
00
00
X
r
1
00
X
r
1
00
X
r
2
00
10
1
0
0
0
C
10
00
X
r
2
00
X
r
2
00
00
X
r
1
00
00
X
r
1
00
X
r
2
X
r
1
X
r
2
00
11
0
0
1
0
C
11
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
X
r
1
X
m
2
00
12
2
0
0
1
A
13
0
1
0
0
C
14
1
0
0
0
ATable 11
Component data for electronic subsystem.
(a)
i
Choice 1 (j=1)
Choice 2 (j=2)
Choice 3 (j=3)
Choice 4 (j=4)
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
1
0.00532
0.004
2
0.000726
0.0012
1
0.00499
0.004
2
0.00818
0.007
3
2
0.00818
0.007
3
0.000619
0.0001
1
0.00431
0.003
2
—
—
—
3
0.0133
0.012
3
0.0110
0.011
3
0.0124
0.011
3
0.00466
0.004
2
4
0.00741
0.006
2
0.0124
0.012
3
0.00683
0.006
2
—
—
—
5
0.00619
0.005
1
0.00413
0.004
2
0.00818
0.007
3
—
—
—
6
0.00436
0.003
3
0.00567
0.005
3
0.00268
0.002
2
0.000408
0.0004
1
7
0.0105
0.01
3
0.00466
0.004
2
0.00394
0.003
2
—
—
—
8
0.0105
0.01
3
0.00105
0.0006
1
0.0105
0.01
3
—
—
—
9
0.00268
0.002
2
0.000101
0.00005
1
0.00408
0.003
1
0.000943
0.00005
1
10
0.0141
0.013
3
0.00683
0.006
2
0.00105
0.0005
1
—
—
—
11
0.00394
0.003
2
0.00355
0.003
2
0.00314
0.002
2
—
—
—
12
0.00236
0.001
1
0.00769
0.007
2
0.0133
0.012
3
0.0110
0.01
3
13
0.00215
0.001
2
0.00536
0.005
3
0.00665
0.006
3
—
—
—
14
0.0110
0.001
3
0.00834
0.003
1
0.00355
0.003
2
0.00436
0.004
3
(b)
i
Choice 1 (j=1)
Choice 2 (j=2)
Choice 3 (j=3)
Choice 4 (j=4)
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
1
10
1
3
3
10
1
4
4
20
2
2
2
20
2
5
5
2
20
2
8
8
10
1
10
10
10
1
9
9
—
—
—
—
3
20
2
7
7
30
3
5
5
10
1
6
6
40
4
4
4
4
30
3
5
5
40
4
6
6
50
5
4
4
—
—
—
—
5
20
2
4
4
20
2
3
3
30
3
5
5
—
—
—
—
6
30
3
5
5
30
3
4
4
20
2
5
5
20
2
4
4
7
40
4
7
7
40
4
8
8
50
5
9
9
—
—
—
—
8
30
3
4
4
50
5
7
7
60
6
6
6
—
—
—
—
9
20
2
8
8
30
3
9
9
40
4
7
7
30
3
8
8
10
40
4
6
6
40
4
5
5
50
5
6
6
—
—
—
—
11
30
3
5
5
40
4
6
6
50
5
6
6
—
—
—
—
12
20
2
4
4
30
3
5
5
40
4
6
6
50
5
7
7
13
20
2
5
5
30
3
5
5
20
2
6
6
—
—
—
—
14
40
4
4
6
40
4
6
7
50
5
6
6
60
6
5
9Table 12
Component data for mechanical subsystem.
Parameters
r
1
2
3
4
5
6
7
8
9
10
11
λ
r
t
—
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004167
0.0004167
0.0004167
λ
-
r
l
l
=
1
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.0005834
0.0005834
0.0005834
l
=
2
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.0008784
0.0008784
0.0008784
λ
̿
r
u
u
=
1
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004167
0.0004167
0.0004167
u
=
2
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005167
0.0005167
0.0005167
α
r
—
2.5
2.5
2.5
2.5
2.5
2.5
2.5
2.6
2.6
2.6
2.4
c
r
l
1
l
=
1
2
1.5
2.5
3
1.8
3.5
3
2
2.5
1.5
2
l
=
2
1
1
1.5
1.5
1
2
1.5
1
1.5
0.5
1
c
r
u
2
u
=
1
4
3
5
6
4
7
6
4
3.5
3
4.5
u
=
2
3
1.5
3
2
2
4
3
1.5
2
1.5
2.5
h
r
l
l
=
1
1
1
2
1
2
1
1
2
1
1
2
l
=
2
2
1
1
1
1
2
1
1
2
1
1
h
r
u
u
=
1
1
2
1
2
2
1
2
2
1
2
1
u
=
2
1
2
2
1
2
2
1
2
2
1
2
c
r
l
3
l
=
1
2
2
2
2
2
2
2
2
2
2
2
l
=
2
1.5
1.5
1.5
1
1
1.5
1.5
1
1
1.5
1.5
c
r
u
4
u
=
1
8
7
6
8
9
7
6
5
7
6
9
u
=
2
3
2
5
4
5
5
2
3
3
2
3Table 13
Upper bound of parameters.
Parameter
T
M
N
max
⁡
i
λ
max.
r
B
0
H
o
W
V
η
r
I
ρ
i
(
t
)
Value
5
3
4
0.0009
1000
20
200
200
0.001
0.03
0.99Figure 4
Reliability-cost Pareto solutions obtained by (a) NSGA-II and (b) MOPSO.
(a)
(b)


## 4.1. Setting Parameters of the Algorithms

Setting proper values for the control parameters of metaheuristic algorithms has a significant effect on their desirable performance. Well-tuned parameters empower the algorithms in producing better solutions within shorter computation times. Thus, setting proper values for control parameters is a critical task [48]. A few additional techniques are applied to tune parameters. One of the common methods is response surface methodology (RSM). RSM is a mathematical and statistical technique that examines the relationship between one or more response variables and the set of parameters (input variables) influencing them. Using this method, the levels of parameters that optimize the response variables are identified. In the first step, main parameters and range of them are determined. In this study, we used related literature to identify these parameters and their range. Table 2 shows range of the main parameters of each algorithm.Table 2
Range of the main parameters.
NSGA-II
MOPSO
MOFA
Parameter
Range
Parameter
Range
Parameter
Range
Pop. size (Npop)
20–100
Pop. size (Npop)
20–100
Pop. size (Npop)
10–100
Max iteration number (Maxit)
50–200
Max iteration number (Maxit)
50–200
Max iteration number (Maxit)
50–200
Cross rate (Cr)
0.5–0.9
Inertia weight (W)
0.4–0.9
Randomization parameter(∝0)
0.1–0.9
Mutation rate (Mr)
0.01–0.3
Cognitive factor(C1)
1-2
Fixed light absorption coefficient(γ)
1–3
Social factor(C2)
1-2In the next step, the response variables should be determined. Three performance metrics are chosen as response variables which are CPU time, number of nondominated solutions (NNS), and Diversification Metric (DM). These metrics were selected based on the two features, the convergence speed and diversity of the detected solutions. More details on the performance metrics can be found in Section4.2.Central composite design (CCD) with 6 center points is applied for the experiments. Experiments are run by Design Expert 9. According to the number of input variables and type of the design, different number of experiments should be run. For instance, in case of the four parameters and 6 center points’ design, 46 experiments are required. After performing the experiments, analysis of variance (ANOVA) is applied to fit an adequate model to the experimental data. Last step is setting goals for responses to generate optimal condition (optimal level of the parameters). Here, we aim to minimize CPU time and maximize NNS and DM.Figure3 shows counter plot of RSM results for MOFA approach. Desirability displays the amount of goals that have been met. It can be concluded from Figure 3 that the highest desirability is obtained when Npop level is medium and Maxit is high. RSM is applied for all three methods and the optimum values for the algorithms parameters are presented in Table 3.Table 3
Optimum values of the algorithms parameters.
NSGA-II
MOPSO
MOFA
Pop. size
100
Pop. size
100
Pop. size
25
Max iteration number
50
Max iteration number
50
Max iteration number
200
Cross rate
0.9
W
0.9
∝
0
0.1
Mutation  rate
0.3
C
1
1
γ
1
C
2
2Figure 3
Counter plot for MOFA: desirability versus Maxit andNpop.Other parameters are set according to the literatures as follows:MOPSO:r1,r2=rand[0,1], repository size = 500 and MOFA: ϵi=[rand-1/2], β0=1.


## 4.2. Performance Measures

According to Deb et al. [49], there are two main features that must be taken into account in order to evaluate performance of metaheuristic algorithms in solving multiobjective optimization problems: (1) convergence to the Pareto optimal set and (2) diversity of the produced Pareto optimal set. Numerous metrics have been proposed for quantifying these features in the literature. Yu and Gen [50] and Zitzler and Thiele [51] proposed some criteria such as number of nondominated solutions (NNS), Error Ratio (ER), and Generational Distance (GD) for measuring accuracy of algorithms. Higher values for NNS and lower values for GD and ER are preferred and show that the algorithm under analysis performed better in providing a set of Pareto optimal solutions. In addition, Diversification Metric (DM), Spacing Metric (SM), and Maximum Spread Metrics (MS) have been proposed to check ability of the algorithms in providing a diverse set of Pareto optimal solutions. In this study we utilized DM and MS to compare algorithms performance which are calculated using(25)DM=∑i=1Nmax⁡xi-yi1/2,MS=∑m=1Mmaxi=1N⁡fmi-mini=1N⁡fmi21/2,where xi-yi is the Euclidean distance between of the nondominated solution xi and the nondominated solution yi and N denote number of the Pareto solutions [52]. The MS presents the distance between the boundary solutions in the obtained nondominated solutions (NDS) where fmi is the mth objective function value for ith Pareto solution.


## 4.3. Numerical Examples

In order to examine performance of the algorithms, three sets of test problems with different sizes (small, medium, and large) are simulated. The dimensions of these sets are shown in Table4. The properties of the test problems were generated using random functions. All three algorithms are programmed in MATLAB and the test problems have been solved on a PC with 4 GB RAM/1.80 GHz CPU. Using the algorithm’s parameters in Table 3, twenty different problems in each set are solved by NSGA-II, MOPSO, and MOFA. All problems are run for five times and the mean of the obtained values for the metrics is considered for algorithm’s performance comparison. In each run, 500 solutions are generated (according to the number of population and maximum iterations) that included infeasible, dominated, and similar solutions. Detected solutions by the algorithms are checked for feasibility firstly and if all the constraints are satisfied, the generated solution is accepted; otherwise another solution would be found. This process has been coded in MATLAB. Mean and standard deviation values of CPU time, NNS, DM, and MS for each case are presented in Table 5.Table 4
Test problems dimensions.
Case
Number of electronic components;i
Number of mechanical components;r
Number of component types in the electronic subsystem;j
Number of component types in the mechanical subsystem;u
Small size
∗
rand
[
3
,
8
]
rand
[
3
,
7
]
rand
[
2,4
]
rand
[
2,4
]
Medium size
rand
[
8
,
13
]
rand
[
7,11
]
rand
[
2,4
]
rand
[
2,4
]
Large size
rand
[
13
,
18
]
rand
[
11,15
]
rand
[
2,4
]
rand
[
2,4
]
∗This function returns an integer number between lower and upper bounds randomly.Table 5
Mean and standard deviation of the metrics for the different sizes of the test problems.
Problem size
Time (seconds)
NNS
DM
MS
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Small
1.91(0.26)
0.78 (0.23)
3.50(1.16)
22.9 (13.54)
12.3(3.45)
3.3(1.26)
55.35(30.22)
60.69 (20.54)
22.60(11.72)
0.01493(0.0257)
0.07999 (0.0542)
0.04230(0.0376)
Medium
2.59(0.21)
1.45 (0.46)
6.99(2.05)
12.85 (5.60)
11.90(3.80)
3.85(1.14)
39.02(17.32)
65.79 (23.50)
28.65(12.71)
0.0238(0.0367)
0.1562 (0.2004)
0.0482(0.0303)
Large
31.49(6.91)
22.35 (9.36)
91.25(25.48)
15.60 (6.22)
12.84(4.19)
4.11(1.49)
51.35(25.77)
97.44 (30.93)
34.35(13.81)
0.00419(0.0048)
0.14625 (0.0285)
0.07408(0.0418)
∗Numbers in the parentheses are standard deviation.Four metrics are considered to evaluate algorithms performance (Table5). These are CPU time, NNS, DM, and MS. The lower values for time and the higher values for the rest of the metrics are desirable. The first rank of the CPU time metric in all cases (small, medium, and large sized instances) belongs to MOPSO approach. It means that MOPSO is more time-efficient. Second rank belongs to NSGA-II and MOFA represents a very long CPU time in comparison with the other competing methods.It also can be concluded from Table5 that the average of NNS in the NSGA-II method is higher than other methods. This means that the NSGA-II can generate more nondominated solutions on average. However, MOPSO has provided close values of this metric to NSGA-II in medium and large sized cases.Obtained values for the DM and MS indicate higher diversity solutions of the produced nondominated solutions by MOPSO in all three cases. This means that the MOPSO method has a wider spread.Our model has two conflicting objectives, reliability and cost. As the reliability of the system goes higher its cost increases too. We try to find an optimal tradeoff between these objectives. To compare quality of solutions generated by three approaches, mean and lower and upper bounds of objectives are reported in Tables6 and 7. It should be noted that the average values of the parameters had been calculated in each case.Table 6
Boundary values for reliability.
Method
Small size problems
Medium size problems
Large size problems
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
0.991548
0.945431
0.96181
0.985369
0.914826
0.938127
0.99101
0.86809
0.90353
Lower bound
0.981583
0.908162
0.939269
0.970327
0.860268
0.909683
0.98771
0.80074
0.86255
Upper bound
0.996518
0.988152
0.980975
0.994157
0.974125
0.957834
0.99190
0.94698
0.93663Table 7
Boundary values for cost.
Method
Small size problems
Medium size problems
Large size problems
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
1713.687
1734.326
1870.389
2944.94
3035.20
3189.41
4725.122
4900.882
5166.778
Lower bound
1652.382
1632.584
1790.569
2891.92
2883.30
3070.47
4662.533
4586.47
5015.414
Upper bound
1848.888
2014.234
1987.349
3065.41
3356.64
3346.01
4889.001
5529.174
5369.501As it can be concluded from Tables6 and 7, MOPSO method provides wider ranges of the objective functions but the quality of the solutions which are detected by NSGA-II is better; in other words, solutions with higher system reliability and lower values of cost were generated by NSGA-II.In remainder of this section, an example is considered based on the proposed model and the results are explained. Tables11 and 12 provide data of failure rate, Weibull distribution parameters, repair, replacement, and downtime costs, and required human resources (number of operators) for each maintenance action in the mechanical section. The electronic section includes 14 subsystems connected in series in which each subsystem can have 3 or 4 types of components. Also, the mechanical section includes 11 repairable and replaceable subsystems. Table 13 presents upper bounds for some problem parameters.This example is solved using three approaches for ten times and the mean and standard deviation values for the performance metrics and objective functions boundary are reported in Tables8 and 9.Table 8
Mean and standard deviation values for performance metrics.
Algorithm
NSGA-II
MOPSO
MOFA
Mean
Std Dev
Mean
Std Dev
Mean
Std Dev
Time
2.74
0.12
1.49
0.06
6.84
0.49
NNS
20.40
5.19
13.20
3.65
3.70
1.06
DM
44.96
10.79
91.46
15.11
30.17
12.21
MS
0.0068
0.0165
0.1436
0.0132
0.0508
0.0218Table 9
Objectives boundary values for the sample.
Method
Reliability
Cost
NSGA-II
MOPSO
MOFA
NSGA-II
MOPSO
MOFA
Mean
0.999421
0.87153
0.90719
3848.62
3900.48
4285.62
Lower bound
0.999825
0.802065
0.841807
3760.333
3485.712
4044.209
Upper bound
0.999833
0.977786
0.951062
4027.349
4488.861
4585.62As concluded before, MOPSO method provides better values for the metrics (time, NNS, DM, and MS) and wider boundary values for the objectives but as Table9 shows, NSGA-II approach detects Pareto solutions lead to higher system reliability and lower cost.The proposed model in this paper aims at finding proper maintenance policies and effective redundancy strategies. Table10 represents one of the nondominated solutions that are found by NSGA-II for the mentioned example. The selected maintenance policies and redundancy strategies including number and type of the redundant components and the implemented redundancy strategy in the electronic section along with the selected repair and replacement actions in the mechanical section are shown in Table 10. The first four columns in the electronic section represent the number of redundant components and the last column indicates the implemented redundancy strategy in each subsystem, where letter A means active strategy and letter C indicates cold strategy. In the mechanical subsystem, Xr and Xm indicate replacement and repair action in each period, respectively. Xr and Xm are followed by a number which determines the type of replacement/repair action. For instance, Xr2 means replacement type two is required and element 00 indicates that no repair or replacement action has been taken in that period. Figure 4 shows Pareto solutions obtained by NSGA-II and MOPSO methods.Table 10
Selected maintenance actions and redundancy strategies by the NSGA-II.
Electronic section
Mechanical section
i
j
r
t
1
2
3
4
Strategy
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
0
0
1
0
A
1
00
00
00
X
r
1
00
X
r
2
00
X
r
2
00
00
X
r
2
X
m
1
00
X
r
1
00
2
1
1
0
0
C
2
00
00
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
X
r
2
X
r
1
X
m
2
00
3
1
0
0
0
A
3
00
X
r
1
00
X
r
1
00
00
X
r
1
X
r
2
00
00
X
r
2
X
r
2
X
m
2
X
r
2
00
4
0
1
0
0
C
4
00
00
00
X
r
2
00
00
X
m
2
00
00
X
m
2
00
X
r
2
X
r
2
X
m
2
00
5
0
1
0
0
A
5
00
X
r
1
00
00
X
r
2
00
00
X
r
1
00
00
X
r
1
00
X
r
2
X
r
1
00
6
0
1
0
1
A
6
00
X
r
2
00
X
r
2
X
r
2
00
00
X
r
2
00
00
X
r
1
00
X
m
2
X
m
2
00
7
0
0
1
0
A
7
00
X
r
2
00
00
X
m
2
00
00
X
r
2
00
00
X
r
1
00
X
r
2
X
m
2
X
r
2
8
0
1
0
0
A
8
00
X
r
2
00
00
X
r
2
X
r
2
00
00
X
m
2
00
00
X
r
2
00
X
r
2
00
9
1
0
1
0
A
9
00
X
r
2
X
m
2
X
r
2
00
00
X
r
2
00
00
X
r
1
00
X
r
1
00
X
r
2
00
10
1
0
0
0
C
10
00
X
r
2
00
X
r
2
00
00
X
r
1
00
00
X
r
1
00
X
r
2
X
r
1
X
r
2
00
11
0
0
1
0
C
11
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
00
X
r
2
00
X
r
1
X
m
2
00
12
2
0
0
1
A
13
0
1
0
0
C
14
1
0
0
0
ATable 11
Component data for electronic subsystem.
(a)
i
Choice 1 (j=1)
Choice 2 (j=2)
Choice 3 (j=3)
Choice 4 (j=4)
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
λ
i
j
1
λ
i
j
2
k
i
j
1
0.00532
0.004
2
0.000726
0.0012
1
0.00499
0.004
2
0.00818
0.007
3
2
0.00818
0.007
3
0.000619
0.0001
1
0.00431
0.003
2
—
—
—
3
0.0133
0.012
3
0.0110
0.011
3
0.0124
0.011
3
0.00466
0.004
2
4
0.00741
0.006
2
0.0124
0.012
3
0.00683
0.006
2
—
—
—
5
0.00619
0.005
1
0.00413
0.004
2
0.00818
0.007
3
—
—
—
6
0.00436
0.003
3
0.00567
0.005
3
0.00268
0.002
2
0.000408
0.0004
1
7
0.0105
0.01
3
0.00466
0.004
2
0.00394
0.003
2
—
—
—
8
0.0105
0.01
3
0.00105
0.0006
1
0.0105
0.01
3
—
—
—
9
0.00268
0.002
2
0.000101
0.00005
1
0.00408
0.003
1
0.000943
0.00005
1
10
0.0141
0.013
3
0.00683
0.006
2
0.00105
0.0005
1
—
—
—
11
0.00394
0.003
2
0.00355
0.003
2
0.00314
0.002
2
—
—
—
12
0.00236
0.001
1
0.00769
0.007
2
0.0133
0.012
3
0.0110
0.01
3
13
0.00215
0.001
2
0.00536
0.005
3
0.00665
0.006
3
—
—
—
14
0.0110
0.001
3
0.00834
0.003
1
0.00355
0.003
2
0.00436
0.004
3
(b)
i
Choice 1 (j=1)
Choice 2 (j=2)
Choice 3 (j=3)
Choice 4 (j=4)
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
c
i
j
c
i
j
′
v
i
j
w
i
j
1
10
1
3
3
10
1
4
4
20
2
2
2
20
2
5
5
2
20
2
8
8
10
1
10
10
10
1
9
9
—
—
—
—
3
20
2
7
7
30
3
5
5
10
1
6
6
40
4
4
4
4
30
3
5
5
40
4
6
6
50
5
4
4
—
—
—
—
5
20
2
4
4
20
2
3
3
30
3
5
5
—
—
—
—
6
30
3
5
5
30
3
4
4
20
2
5
5
20
2
4
4
7
40
4
7
7
40
4
8
8
50
5
9
9
—
—
—
—
8
30
3
4
4
50
5
7
7
60
6
6
6
—
—
—
—
9
20
2
8
8
30
3
9
9
40
4
7
7
30
3
8
8
10
40
4
6
6
40
4
5
5
50
5
6
6
—
—
—
—
11
30
3
5
5
40
4
6
6
50
5
6
6
—
—
—
—
12
20
2
4
4
30
3
5
5
40
4
6
6
50
5
7
7
13
20
2
5
5
30
3
5
5
20
2
6
6
—
—
—
—
14
40
4
4
6
40
4
6
7
50
5
6
6
60
6
5
9Table 12
Component data for mechanical subsystem.
Parameters
r
1
2
3
4
5
6
7
8
9
10
11
λ
r
t
—
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004167
0.0004167
0.0004167
λ
-
r
l
l
=
1
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.00056
0.0005834
0.0005834
0.0005834
l
=
2
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.00085
0.0008784
0.0008784
0.0008784
λ
̿
r
u
u
=
1
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004167
0.0004167
0.0004167
u
=
2
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005
0.0005167
0.0005167
0.0005167
α
r
—
2.5
2.5
2.5
2.5
2.5
2.5
2.5
2.6
2.6
2.6
2.4
c
r
l
1
l
=
1
2
1.5
2.5
3
1.8
3.5
3
2
2.5
1.5
2
l
=
2
1
1
1.5
1.5
1
2
1.5
1
1.5
0.5
1
c
r
u
2
u
=
1
4
3
5
6
4
7
6
4
3.5
3
4.5
u
=
2
3
1.5
3
2
2
4
3
1.5
2
1.5
2.5
h
r
l
l
=
1
1
1
2
1
2
1
1
2
1
1
2
l
=
2
2
1
1
1
1
2
1
1
2
1
1
h
r
u
u
=
1
1
2
1
2
2
1
2
2
1
2
1
u
=
2
1
2
2
1
2
2
1
2
2
1
2
c
r
l
3
l
=
1
2
2
2
2
2
2
2
2
2
2
2
l
=
2
1.5
1.5
1.5
1
1
1.5
1.5
1
1
1.5
1.5
c
r
u
4
u
=
1
8
7
6
8
9
7
6
5
7
6
9
u
=
2
3
2
5
4
5
5
2
3
3
2
3Table 13
Upper bound of parameters.
Parameter
T
M
N
max
⁡
i
λ
max.
r
B
0
H
o
W
V
η
r
I
ρ
i
(
t
)
Value
5
3
4
0.0009
1000
20
200
200
0.001
0.03
0.99Figure 4
Reliability-cost Pareto solutions obtained by (a) NSGA-II and (b) MOPSO.
(a)
(b)


## 5. Conclusion

In this paper, a biobjective reliability model by integrating redundancy allocation problem (RAP) and reliability-centered maintenance (RCM) problem for a system of nonrepairable electronic components and repairable mechanical components has been proposed. Objectives of the problem are maximizing the system reliability and minimizing the system operational and maintenance costs. In order to improve system reliability, active and cold standby redundancy strategies and periodic maintenance actions are considered for electronic section and mechanical section, respectively. Total system cost includes initial costs for purchasing the required equipment and secondary costs such as operational costs of the electronic section and maintenance costs of the mechanical section. Initial costs are taken into account by setting a budget constraint and operational costs are considered as the second objective of the problem. Three metaheuristic algorithms, NSGA-II, MOPSO, and MOFA, are used to solve the proposed model. Different sets of test problems (small, medium, and large size) have been generated for evaluating solution methods. Obtained results indicate that MOPSO algorithm requires less time to produce Pareto optimal sets with high diversities according to the DM and MS. NSGA-II outperforms MOPSO and MOFA in terms of generating more nondominated solutions (NNS) with better values for system reliability and cost. Finally an example was solved and the obtained results were explained.



---

*Source: 396864-2015-08-02.xml*