A simulation-based optimization approach for free distributed repairable multi-state availability-redundancy allocation problems

Contents lists available at ScienceDirect Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress A simulation-based optimization approach for free distributed repairable multi-state availability-redundancy allocation problems Ahmad Attar, Sadigh Raissi, Kaveh Khalili-Damghani ⁎ School of Industrial Engineering, South Tehran Branch, Islamic Azad University, Tehran, Iran ARTICLE INFO Keywords: Availability-redundancy allocation Simulation based optimization (SBO) Multi-state systems Distribution-free failure/repair time Non-dominated sorting genetic algorithm Repairable systems ABSTRACT A simulation-based optimization (SBO) method is proposed to handle multi-objective joint availabilityredundancy allocation problem (JARAP). Here, there is no emphasis on probability distributions of time to failures and repair times for multi-state multi-component series-parallel configuration under active, cold and hot standby strategies. Under such conditions, estimation of availability is not a trivial task. First, an efficient computer simulation model is proposed to estimate the availability of the aforementioned system. Then, the estimated availability values are used in a repetitive manner as parameter of a two-objective joint availabilityredundancy allocation optimization model through SBO mechanism. The optimization model is then solved using two well-known multi-objective evolutionary computation algorithms, i.e., non-dominated sorting genetic algorithm (NSGA-II), and Strength Pareto Evolutionary Algorithm (SPEA2). The proposed SBO approach is tested using non-exponential numerical example with multi-state repairable components. The results are presented and discussed through different demand scenarios under cold and hot standby strategies. Furthermore, performance of NSGA-II and SPEA2 are statistically compared regarding multi-objective accuracy, and diversity metrics. 1. Introduction Designing reliable systems is a vital optimization task [1]. Many authors have addressed such problem with different assumptions. In general, three categories are recognized for reliability optimization. They are: “reliability allocation [2–5]”,“redundancy allocation [6–11]”, and “joint reliability-redundancy allocation [12]”. The most common system structure for reliability optimization problems is series-parallel which can represent a wide verity of real systems [2,4,6,10,13–15]. In a binary-state series-parallel system, each component has two possible states: perfect functioning and complete failure [11]. However, real systems, called multi-state, can usually take some medial states [1–4,13,14,16,17]. In this paper, we focus on extending the multi-state joint availability-redundancy allocation problem (JARAP) proposed by [13] to a more practical environment. In practice, the redundancy strategy can be either active or standby [18,19]. However, almost all of the JARAP models in the literature have considered active redundancy [20]. Due to simplicity, the majority of analytic models has assumed the exponential distribution as the probability density function (PDF) of failure and repair time [13]. Bowles [23] pointed out the consequences of using such simplifications and concluded that such assumptions can lead to significantly unacceptable reliability estimations. The Monte Carlo simulation (MCS) and discrete event simulation (DES) are the two possible replacements that countenance all forms of distributions and can handle realistic situations [24,25,16]. The multi-objective approaches of JARAP were also mostly limited to the weighting technique and the traditional utility function. Yet, it is widely-known that these aggregation methods decline the chance of finding multiple non-dominated (Pareto) solutions [21]. High quality Pareto-optimal solutions with outstanding grade levels of diversity can be evolved using multi-objective evolutionary algorithms (MOEA) such as the non-dominated sorting genetic algorithm (NSGA-II) [10,22].In the present study, the JARAP proposed by [13] is extended regarding the following issues: (1) admittance of either of cold standby and hot standby strategies; (2) introducing an accurate time-efficient objectoriented DES; (3) providing a multi-objective model for the JARAP and considering availability and total cost as objective functions; and (4) customizing NSGA-II as an efficient solution algorithm and combining it with the proposed DES using advanced programming elements. The sections of the paper are organized as follows. Literature survey is presented and classified in Section 2.InSection 3, the proposed availability estimation method is presented. In Section 4, a twoobjective mathematical model is developed. The proposed SBO methhttp://dx.doi.org/10.1016/j.ress.2016.09.006 Received 10 August 2015; Received in revised form 11 July 2016; Accepted 18 September 2016 ⁎ Correspondence to: Department of Industrial Engineering, Faculty of Industrial Engineering, South-Tehran Branch, Islamic Azad University, Tehran, Iran. E-mail addresses: [email protected] (A. Attar), [email protected] (S. Raissi), [email protected],[email protected] (K. Khalili-Damghani). Reliability Engineering and System Safety 157 (2017) 177–191 0951-8320/ © 2016 Elsevier Ltd. All rights reserved. Available online 20 September 2016 crossmark odology is presented in Section 5. The benchmark instances and statistical analysis are presented in Section 6. And finally, concluding remarks and future research directions are provided in Section 7. 2. Literature review 2.1. Reliability-redundancy allocation models In general, the JARAP has gained far less attention in comparison with the simple redundancy allocation approach. We classify the existing models of the joint approach based on four fundamental attributes: (a) objective function; (b) components characteristics; (c) redundancy strategy; and (d) solution methodology. Components usually receive three characteristics that are graphically shown in Fig. 1 as the characteristics’triangle; namely, (1) states (binary or multi-state), (2) repair conditions, and (3) failure/repair distribution function. Misra and Ljubojevic [12] introduced the first JARAP model for a series-parallel system. The JARAPs with non-repairable binary components were also studied in [6,20,26–33] and. Elegbede and Adjallah [34] assumed repair capable binary components with exponential function. Tian et al. [14] modeled the JARAP for engineering systems where each component could take more than two states. Determining the relationship between component cost and its state distribution was still rather difficult in real practice [13]. Tian et al. [13] proposed a practical approach to cover systems in which components were both multi-state and repairable. Recently, Hamadani and Khorshidi [35] and Khorshidi and Nikfalazar [36] added time value of money to the model proposed in [13] and considered cost and availability as objective functions. The repairable components, was surveyed in [13,35,36]. To our best knowledge in the JARAP literature, the only model with cold-standby strategy is the one proposed by Ardakan and Hamadani [20]. Nevertheless, it still has very simplified component characteristics; i.e., binary-state non-repairable components with exponential distribution. 2.2. Availability estimation techniques Applying the Markov chain is a common task in availability estimation and redundancy allocation problems [17]. More advanced methodologies such as universal generating function (UGF) were also used [13,35]. The main deficiency of such analytic approach turns back a property that is usually characterized as "memory-less-ness": the probability distribution of the next state depends only on the current state and not on the sequence of events that preceded it. Therefore it is limited to use of the exponential probability distribution for time to failure (TTF) and time to repair (TTR). If these conditions are not ruling, simulation-based methodology is frequently recommended [24]. Simulation methods are time consuming. In order to resolve this problem, different simulation approaches were introduced in [5,16,25]. Marseguerra et al. [25] proposed a SBO approach called “drop-bydrop”based on Monte Carlo simulation (MCS) theory and Genetic algorithm (GA). This approach was mainly established upon this fact that, in GA good chromosomes usually appear in successive generations (Elitism). Lins and Droguett [5] offered a stepwise DES method for estimating the availability of repairable binary systems. Apart from the mentioned issues, the above declared methods still take at least few seconds to estimate the reliability/availability of a medium-size system [16]. 2.3. Solution methodologies Both JARAP and RAP are known to be NP-hard problems [39]. There are several records of successful applications of meta-heuristic methods such as GA, immune system (IS) [30], particle swarm optimization (PSO), cuckoo search (CS) [32], imperialist competitive algorithm (ICA), non-dominated sorting genetic algorithm (NSGA-II) [9], niched Pareto genetic algorithm (NPGA) [41] to solve JARAPs and RAPs considering non-repairable binary-state components. In general, a good multi-objective algorithm should provide: (i) high quality non-dominated solutions (i.e., accurate estimation of the Pareto front), and (ii) proper diversity on the Pareto front [42]. NSGAII, which was introduced by Deb et al. [22], is one of the most successful algorithms for solving reliability optimization [10,43,9,19,44,45,46]. 2.4. Summary and motivations The literature of the JARAP is summarized in Table 1. Here, studies are listed in chronological order and the classification is made based on the attributes and the characteristics that were illustrated in Fig. 1. As seen in Table 1, the majority of the JARAPs contemplate availability as objective function. Multi-objective variants of JARAP are not frequent. From component characteristics standpoint, only few researches have considered multi-state or repairable components and none of the past studies in this field have considered a non-exponential PDF for the failure/repair times. Moreover, except for one recent research (i.e., [20]), all of the JARAPs were restricted to active redundancy strategy. Notably, none of the MOEAs were employed for solving multi-state repairable models of JARAP, The last row in Table 1 demonstrates the assumptions and characteristics of the this study. 3. The proposed methodology for availability estimation In order to estimate the availability, the following assumptions are made in this study: •System consist of multi-state components and repair allowed independently by a single repair team. •The state transitions are freely distributed and can follow any given probability density function. •Components might have different performance rates (outputs) in different states. •Switching acts perfectly that means the reliability of the switch is assumed to be equal to 100%. The notations used in this study are presented in Table 2. Fig. 1. An overview of attributes of a joint reliability-redundancy optimization problem. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 178 Table 1 An overview of the past works in joint reliability-redundancy allocation problem. Research paper Objective Component characteristics Redundancy strategy Solution methodology Algorithm a States Distribution Reparability Cost Availability Multiple Binary-state Multi-state Exponential All forms Non-Repairable Repairable Active Standby Exact Heuristic Meta-heuristic [12] ☑☑☑☑ ☑ ☑ [26] ☑☑ ☑ ☑ ☑ ☑ [34] ☑☑ ☑ ☑ ☑ ☑ GA [33] ☑☑☑☑ ☑ ☑ [27] ☑☑ ☑ ☑ ☑ ☑ [6] ☑☑☑☑ ☑☑ DP [9] ☑☑ ☑ ☑ ☑ ☑ NSGA-II [28] ☑☑ ☑ ☑ ☑ ☑IS [14] ☑☑☑ ☑ ☑ ☑GA [41] ☑☑ ☑ ☑ ☑ ☑ NPGA [13] ☑☑☑☑☑☑GA [29] ☑☑ ☑ ☑ ☑ ☑PSO [30] ☑☑ ☑ ☑ ☑ ☑IS [40] ☑☑ ☑ ☑ ☑ ☑PSO [35] ☑☑☑ ☑☑ ☑GA [20] ☑☑ ☑ ☑ ☑ ☑GA [36] ☑☑☑ ☑☑ ☑GA, ICA [32] ☑☑ ☑ ☑ ☑ ☑CS Present study ☑☑ ☑ ☑☑ ☑NSGA-II a GA: Genetic algorithm; DP: Dynamic programming; NSGA-II: Non-dominated sorting genetic algorithm; IS: Immune system; NPGA: Niched Pareto genetic algorithm; PSO: Particle swarm optimization; CS: Cuckoo search; ICA: Imperialist competitive algorithm. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 179 Since availability calculation method is assumed as the core of a JARAP model, this part is described first. Different simulation software is available for estimating the system availability with non-exponential transition PDFs. In this study “Enterprise Dynamics™(ED)”is applied for object oriented computer simulation modeling. Features and successful applications of ED simulation software can be found in recent studies such as [47,48]. 3.1. Design and development of multi-state server atom In order to achieve better results four simulation tricks are used using designing new ED atoms. This atom models a multi-state component. There are two major events, i.e., Failure and Repair, that randomly happen in multi-state components. Failure event is the event through which the component goes to another state with a lower performance rate. Improving the component condition and restoring its performance rate is called repair event. Also as a common assumption, transitions are considered to happen among neighbor states. Therefore, for state j, we denote the time-to-failure and time-torepair distributions by TTF j and TTR j respectively, where jtakes values from 1 (perfect functioning) to m(complete failure). A typical events for a three-state component is illustrated in Fig. 2. We imitate this process in our Multi-state Server atom by defining two events in the event handler of this object as depicted in Fig. 3. When a failure event occurs, the first flowchart (i.e., Fig. 3a) is executed. This flowchart degrades the component condition by increasing the state indicator j by 1, and schedules a repair event to occur Ri time units later, where Ri is randomly generated based on the TTR distribution of the current state of the component. Next, the flowchart schedules a failure event to occur in F i time units, where F i is randomly generated based on the TTF distribution of the component in the current state. Since the component faces no more failures in the complete failure state (i.e., jS= i ), scheduling the new failure event is restricted to the operating states (i.e., jS< i ). When a repair event occurs, the second flowchart is executed (i.e., Fig. 3b). This flowchart, first, improves the component conditions by decreasing j by 1. Then, it schedules a failure event and a repair event to occur in F i and Ri time units, where F i and Ri are defined similar to the failure event. Since the component condition cannot be improved more than perfect functioning (i.e., j= 1 ), the repair event is only scheduled for the degraded states (i.e., j> 1 ). Note that, we did not need to set such condition in the failure flowchart for scheduling of repair events. Thus, it is assure that repair events are happened in state jsuch that jS 2 ≤≤ i . Moreover, when a component is set on cold-standby mode, its condition and the associated performance rate should not be degraded. Therefore, in the designed atom, once a component becomes standby the failure schedule is paused and the occurrence of the failure event is postponed until the component returns active. For the same reason, scheduling the failure event in the repair flowchart is only considered for active components. In order to make sure that each component has at most one failure event and one repair event scheduled at a time, the previously scheduled events are eliminated at the beginning of flowcharts of Fig. 3. The triggering process is illustrated in Fig. 3c, in which the failure and repair events are initialized based on the present conditions of the component. One of the advantages of ED software and the designed atom is the resemblance of its 2D layout to the common reliability diagrams. For instance, an example series-parallel system is presented in Fig. 4 using ordinary reliability diagram and the associated simulation model in ED software. As seen in Fig. 4, the 2D representation of the model in ED software is very similar to the ordinary reliability diagram, which is usually used to represent the series-parallel systems. This feature improves the readability of the simulation models layout. 3.2. Design and development of subsystem atom The intentions of creating this atom are: to link parallel components, to aggregate the total performance rates/outputs of each subsystem, and to determine its condition (up or down) regarding the desired demand. In order to connect each Multi-state Server to the associated subsystem, all needed is to connect its central channel to an input channel on the subsystem. Fig. 5 represents the 2D preview of the three proposed atoms. As seen in Fig. 5b, the Subsystem atom displays the real-time performance rate (PR) of the subsystem, while its color reveals the current status of the subsystem: green for up and red for down, where down means the PR has dropped under the desired demand level. 3.3. Design and development of switcher and availability calculator atom Switching between components of a subsystem and calculating the Table 2 Definitions of the used notations. N Number of subsystems in the system S i Number of states available for the components in subsystem i V i Number of components versions for subsystem i K i s Number of subsystem-level actions in subsystem i K i c Number of component-level actions in subsystem i D The desired demand level for the system C i 0 The overhead cost of installing components in subsystem i Cj iThe unit cost of purchasing component version jfor subsystem i V C jl i The unit cost of applying action l on component version jin subsystem i F C j l i The overhead (fixed) cost of applying action l on component version j in subsystem i CiThe overall cost of subsystem i Cco m iThe total component cost of subsystem i Cact iThe total action cost of subsystem i R j m i The performance rate of component version jin subsystem iin state m AF l m The effect of action l on the failure distribution of the component in state m AR l m The effect of action l on the repair distribution of the component in state m G q i Action group q of subsystem i Q iNumber of action groups defined in subsystem i At( ) Instantaneous availability of the system at time t At() A Achieved availability of the system at time t A s Steady-state availability of the system x ij An integer variable for the number of components of version j installed in subsystem i y c jl i A binary variable representing whether component-level action l is applied to components of version jin subsystem i y sliA binary variable representing whether subsystem-level action l is applied to subsystem i XThe redundancy vector Y The action vector CX Y(, ) The total cost of the system AX Y(, ) The total availability of the system Fig. 2. The schematic overview of the occurrence of failure and repair events over time in a multi-state component with three states. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 180 availability of the system need access to all subsystems. For this reason, an integrated atom is designed and developed for handling these two tasks concurrently. The switching process acts immediately by considering perfect switching. Therefore, the subsystem does not fail unless all of its components fail [37,38]. In real systems, sensors and signaling mechanisms are usually used to make the instant switching possible [18,38]. Here, a similar signaling mechanism is defined using “create event”function in ED software. In this mechanism, whenever a component changes its state, a signal is sent toward the switcher. The switcher checks the subsystem and performs the switching if needed. Then, the new subsystem status (up or down) and the related PR is determined by the Subsystem atom. Afterward, based on the new PRs and status of all subsystems, the overall system status and PR is determined in the “Switcher and avail. calc.”atom. Using the proposed signaling mechanism, the real-time estimation of availability with the maximum accuracy and low computational efforts is proposed. Instantaneous availability (i.e., At( ) ) is probability that the system is up at time t. Estimation of availability for long time is commonly used [1] [13,35]. Eq. (1) is used to calculate instantaneous availability [49,50]. AAt= lim ( ) st→∞ (1) Estimation of availability using Eq. (1) is quite easy for Markovbased models ([13,17,50]). However, when a simulation model is used, obtaining an acceptable estimation for At( ) in the steady-state condition requires performing long simulation runs. Hence, most of the simulation-based studies have neglected the steady-state availability and focused on small t values. In this study, another type of availability Fig. 3. The flowcharts of the failure event and repair event and the flowcharts of the triggering process of the designed Multi-state Server atom. Fig. 4. The reliability diagram and the associated ED object-oriented simulation model for an example series-parallel system with two multi-state subsystems and hot standby strategy. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 181 is used to overcome this limitation “Achieved Availability”, which was defined for real running systems [49], is estimated using Eq. (2). AUp time Total Time = A (2) Since the simulation model imitates the real system behavior, the A A value is calculated and reported in the Switcher and avail. calc. atom. The 2D representation designed for this atom is shown in Fig. 5c. The long-run achieved availability of a system is equal to its steadystate availability [49]. 3.4. Model validation and speed checking The validity of the proposed simulation procedure is checked using a benchmark test problem adopted from [13]. The system in this test problem consists of two subsystems, where the 1st subsystem has three component versions and the 2nd subsystem has four component versions. For each subsystem, the number of redundant components is given in Table 3. The desired demand is 1000 and the reported availability for this example is 95.39%, which was calculated in [13] using the Markov and the UGF methods. After parameter setting, 30 simulation runs based on 10,000 observations period are carried out and the first 10% time missed as warm-up to avoid any biased estimation. It is worthy that each run did not take more than one second long due to the designed especial atoms, which is considerably more time-efficient than the CPU-time reported by other simulation methods in the literature [16]. The results of the availability estimations in these 30 runs are summarized in Table 4.As Table 4 shows, the average simulated availability (i.e., 95.41%) is very close to the exact value reported in [13]. Since the behavior of the estimated availability values (Table 4) follows a normal distribution (P-value=0.196 for Kolmogorov-Smirnov test), we undertake a two-tailed one-sample t-test for comparing these results with those of [13]. This statistical test checks the null hypothesis that the mean value of the simulated availabilities is equal to the one reported by [13] (H 0 :µ=95.39%, against the alternative hypothesis H 1 :µ≠95.39%. The described hypothesis test resulted in p-value of 0.623, which indicates that there is no evidence to reject H 0and the simulation results are statistically acceptable and comparable with those reported in [13]. Although the output of the proposed simulation method was very promising in this JARAP benchmark test problem, another Markovian numerical example with standby components is used to check the accuracy of the cold-standby part of the proposed simulation method. Consider a system with cold-standby redundancy strategy which has two identical repairable components connected in parallel with λ and μ as the failure and repair rates, respectively. The Markov state transition diagram of this system is given in Fig. 6. The first and second letters in each state indicate the conditions of the first and second component in that state; where A,S, and Rare reserved for active, standby, and under repair (failed) conditions, respectively. For instance, in state 2 the first component is under repair and the second one is active [17]. The instantaneous availability of the system is given by Eq. (3) where Pt( ) 5is the probability of having both of the components under repair at time t .APt=1− ( ) t() 5 (3) (a) Multi-State Server (b) Subsystem (c) Switcher and availability calculato r Fig. 5. An overview of 2D representations of the designed atoms in ED. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.). Table 3 Benchmark example adopted from [13]. Subsystem Version Redundancy Actions 1 1 5 {8} 2 8 {8} 3 1 {8} 2 1 4 {4} 2 2 {4} 3 2 {3,4} 4 2 {4} Table 4 Summary of simulation results and statistical tests of the benchmark problem for validating the proposed simulation model under a hot standby JRRAP. Summary of simulation results and statistical tests I. Estimated Availability Values Sample size Average Std. Dev. 30 95.41% 0.003 II. Normality test K-S statistic p-value 0.132 0.196 III. t-test T statistic p-value 0.498 0.623 Fig. 6. State transition diagram of a standby system with two parallel components. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 182 The steady-state availability of the system is estimated using both the proposed simulation method and the Markov method (using Eqs. (3), and (1), respectively) and the results are reported in Table 5. The results of cold-standby test problem are also compared through a t-test. The results of statistical analysis are given in Table 5. It is observed that the p-values of the test is significantly high, which confirms the equality of the results of the two methods. Based on the results of Tables 4 and 5, it can be concluded that the proposed simulation method is accurate enough to be used for estimating the availability of series-parallel systems under both cold and hot standby strategies. 4. Mathematical optimization model Multi-objective model (4)–(9) is proposed for the JARAP CX Y AX Y { Min Total Cost = ( , ), Max Total Availability = ( , )} (4) ∑yc iN jV qQ≤1; ∀(1≤≤,1≤≤,1≤ ≤ ) lG jl iii ∈q i(5) ∑ys i N q Q≤1 ; ∀ (1≤ ≤ , 1≤ ≤ ) lG l ii ∈q i(6) x ziNjV∈ ∪{0} ; ∀ (1≤ ≤ , 1≤ ≤ ) ij i +(7) y ciNjVlK∈{0,1}; ∀(1≤≤,1≤≤,1≤≤ ) jl iii c(8) y siNKlKK∈{0,1} ; ∀(1≤ ≤ , +1≤ ≤ + ) l i i c i c i s(9) where, N represents the total number of subsystems, V i stands for the number of component versions available for subsystemi. Moreover, for subsystemi, number of component-level actions and number of subsystem-level actions are denoted by K i c and K i s , respectively. The design variable (i.e., redundancy Xand action Y ) are defined in Eqs. (10)–(13). XXX X=( , ,…, ) N12 (10) Xxx x= ( , ,…, ) iii iv12 i (11) Y YY Y=( , ,… ) N12 (12) Yyc ycyc yc yc ycys ys = ( ,…, , ,…, ,…, ,…, , , …, ) ii K ii K i V i VK i K i KK i 11121 2 1 +1 + icicii icic icis(13) Here, Xiand Y idenote the redundancy and action vectors of the ith subsystem, x i j is a non-negative integer variable which symbolizes the number of components of type jin subsystem i. For subsystemi, y cj l iand y s l i are binary variables that represent performing component-level actions and subsystem-level actions, respectively. The qth group of mutually exclusive actions in subsystem iis denoted by G q i , and Q i is the total number of groups defined in subsystem i. The term “technical and organizational actions”comprises managerial decisions that can increase the TTF or reduce the TTR of the components. For instance, providing advanced training courses for the machine operators to deal with a specific version of machines may lead to better usage and operations or earlier fault detections. On the other hand, providing advanced repair courses for the repair crew, doubling the repair team or even considering extra working shifts for the maintenance department can improve the repair characteristics of the components. In addition, proper spare parts supply may play an important role in reducing the repair time. So, adopting a certain spare part inventory control system can also be considered as a technical and organizational action. Eq. (4) presents the conflictive objective functions of the model (i.e., minimizing the total cost of the system and maximizing the total availability of the system). Technical and organizational actions might be “independent”or “mutually exclusive”[13]. For instance, doubling the maintenance crew and increasing the crew by 50% are mutually exclusive actions [13]. Constraints (5) and (6) restrict the model to choose at most one action from among all mutually exclusive actions. Finally, constraints (7)–(9) define the domain of decision variables. Total cost of a series-parallel system is the summation of subsystems’ costs which is calculated using Eq. (14). Subsystem cost (Eq. (15)) includes the costs associated with components ( C com i ) and costs of applying technical and organizational actions ( C ac t i ). ∑ CX C()= i N i =1 (14) CC C=+ icom i act i(15) Component cost of Eq. (16) is consisted of the related fixed cost of installing components in the subsystem (Ci 0) and the installation/ purchase cost of each component. Action cost of subsystem iis given in Eq. (17). The first part of Eq. (17) represents the fixed (FCjl i) and variable ( V Cjl i) costs of applying component-level actions and the second part is related to the subsystem-level actions. ∑ CC xC=+ com ii j V ij j i 0 =1 i (16) ∑∑ C yc FC x VC ys FC=(+)+ act i l K jl i jl iij jl i lK KK l i jl i =1 = +1 + ic ic icis (17) 5. Solution methodology NSGA-II [22] is adopted to solve the proposed optimization model (4)–(9). The mechanism of NSGA-II is graphically illustrated in Fig. 7. Fig. 7 shows that solutions of each front are sorted based on a diversity measure called crowding distance, which determines how far each solution is from its neighbors in the front. Further details can be found in [22]. 5.1. Solution representation and handling the constraints Defining suitable solution representation and handling the constraints of the model (4)–(9) are essential steps in applying NSGA-II. A Table 5 Summary of comparison of the steady-state availability of the proposed simulation model and Markov method for a cold standby example for the validation purpose under different scenarios and the related statistical tests results. λ =0. 3 λ =0. 5 μ =0. 5 μ =0. 8 μ =0. 3 μ =0. 8 A S (Markov Model) 0.8989 0.9514 0.6575 0.8927 A S (Simulation Model) 0.8995 0.9514 0.6580 0.8931 P-value (t-test) 0.371 0.986 0.577 0.599 Eliminated Pt Ot F1 F3 F2 Non-dominated sorting Sorting by Crowding distance and selection Pt+1 Fig. 7. NSGA-II mechanism adopted form [22]. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 183 fully customized chromosome is proposed to codify the solution of the model (4)–(9). Each chromosome has N parts, where N is the number of subsystems. The proposed chromosome is presented in Fig. 8.As seen in Fig. 8, each part contains the redundancy levels (i.e., x i j ) and the component-level actions (i.e., y cj l i) of each component version along with the subsystem-level actions (i.e., y s l i ) of the subsystem. Constraints (7)–(9) can easily be considered in the structure of the chromosome through defining the sufficient values for each gene (i.e., non-negative integer and binary for the redundancy and action selection variables, respectively). But, handling constraints (5) and (6) for the mutually exclusive actions is not a trivial task. Different methods like imposing penalty costs for constraint violations can be used for applying constraints (5) and (6) [8]. An easy to implement transformation mechanism that keeps all of the solutions produced by the algorithm in the feasible area is proposed. In this mechanism, the selected action in each group is determined by the summation of the binary variables associated with the actions in that group. Considering u q as the summation of binary variables associated with the actions in group q, the selected action of the group will be the u q th member of this group. In case u=0 q, none of the actions in this group will be selected. For instance, assume a group with 4 actions and let y 1 to y 4 be their associated decision variables. For such a group, the transformation mechanism is shown in Fig. 9 for all possible situations. Using the proposed mechanism, the acceptable chromosome is extracted just before calculating the objective values while the raw chromosome is used for the crossover and mutation operations. This assures that the transformations do not affect the outcome of these operations and the diversity of the individuals is retained for the next population. 5.2. Combining the simulation model with NSGA-II ED Simulation software offers an object-oriented modeling environment and supports event-based simulation programming that helps building fast simulation models. Meta-heuristic algorithms are efficiently implemented in MATLAB™software. NSGA-II is coded in MATLAB software. The main challenge is connecting and interaction between ED and MATLAB soft-wares. We take advantage of the ActiveX standard to transfer commands and data between ED and MATLAB. MATLAB software is programmed to control the main steps of interaction every time that we have to estimate the availability of a new chromosome, i.e. hundreds of times in each iteration. Fig. 10 illustrates the combination of the simulation model and the optimization algorithm. At the beginning of the proposed NSGA-II procedure, the initial population is generated randomly. Then objective function values are calculated for each individual; the availability of the system is estimated using the simulation model described in Section 3. Next, the offspring set ( O t) is formed by applying single point crossover and swap mutation operators, where the parents for each individual in O t are selected from P t using tournament selection scheme. By combining O tand P t , R t is created and then the new population is selected from the existing individuals in R t based on the non-dominated sorting and crowding distance concepts. The algorithm keeps repeating these steps until the stopping criteria are met. A maximum number of iterations is defined as stopping criterion. The framework of the proposed SBO and its transactions with the mathematical model are graphically summarized in Fig. 11. As seen in Fig. 11, the proposed SBO is designed modularly. So, each module of the proposed SBO can alternatively be replaced with another self-created or existing methods. For instance, another optimization algorithm can be replaced with NSGA-II. 6. Numerical results and discussion The numerical example discussed in [13] and [35] is customized for non-exponential failure and repair PDFs under standby situations. Like the original example in [13] and [35], a system is assumed with two subsystems, three component versions for the first subsystem and four component versions for the second subsystem. The number of states available for the components in subsystem 1 and 2 are assumed to be 3 and 2, respectively (as assumed in [13]). The fixed costs (Ci 0) are also assumed 50 and 60 for the first and second subsystems, respectively. The first subsystem has 8 actions, where actions 1–5 are componentlevel actions and actions 6–8 are subsystem-level actions. These actions are divided into the following four distinct groups (i.e., Q = 4 1): G ={1,2} 1 1, G ={3,4} 2 1, G ={5} 3 1, G ={6,7,8} 4 1.Table 6 contains the unit cost of the components (Cj i ), their performance rates (Rj m i), and the failure and repair PDFs in each state. The effect of the actions on component versions 1–3 for subsystem Fig. 8. Structure of designed chromosome for solution representation in NSGA-II. y1y2y3y4 0 i y In the raw chromosome In the transformed chromosome 0000 1000 0100 0010 0001 2 i y 3 i y 4 i y 1 i y Fig. 9. A schematic view of the raw and transformed decision variables for an example group of actions with four members. Fig. 10. The flowchart of the simulation-optimization algorithm. A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 184 1 are given in Table 7. Note that, in this example, the technical and organizational actions only influence the expected values of the state transition distributions. So, the distribution type and other specifications of the components are kept unchanged. Here, FCjl iand V Cjl iare the fixed and variable costs of applying the actions, respectively. For state m ,AF lm and ARlm represent the effect of action lon the expected values of the failure and repair PDFs, respectively. For subsystem 2, three component-level actions and one subsystem-level action are assumed and all of the actions are independent which means Q = 4 2and G ={1} 1 2, G ={2} 2 2 , G ={3} 3 2, G ={4} 4 2 . Table 8 presents the performance rates, failure and repair PDFs, and the associated unit costs. Similar to the first subsystem, the effects of the actions on the versions of components in the second subsystem are given in the Table 9. This example is investigated under six different scenarios. Three demand levels (i.e., 600, 1000 and 1200 units) as well as cold and hot standby strategies are considered. The best values for the parameters of NSGA-II are determined experimentally. The total iteration number is set equal to 200, the tournament size is set equal to 5, and population size is set equal to 100. The crossover and mutation rates are set equal to 0.8 and 0.3, respectively. Fig. 12 presents the estimated Pareto fronts of the numerical example for all six scenarios based on system availability and total cost of the system. As observed in the Pareto fronts of Fig. 12, solutions are perfectly spread over the solution space which reveals that the proposed algorithm was successful in maintaining the diversity. Density of the solutions in the optimal Pareto fronts is acceptable. A wide range solutions is also included for each scenario, from availability around 0 to the highest possible availability, which provides a wide interval of total cost for the decision maker to choose from. A sensitivity analysis is accomplished on standby strategy and the demand level. Since analyzing the Pareto fronts graphically is not Fig. 11. The conceptual framework of the proposed hybrid approach. Table 6 Parameters of components of subsystem 1. Versio n jCj 1 R j 1 1 R j 2 1 Time-to-Failure (TTF) Time-to-Repair (TTR) State 1 State 2 State 2 State 3 1 2 060 0 30 0 Weibull (12, 1. 5 ) a Weibull (15, 1. 8 ) Weibull (14, 0. 8 ) Weibull (19, 0. 30 ) 2 25 1000 5 0 0 Weibull (11, 1. 9 ) Weibull (14, 1. 7 ) Weibull (13, 0. 7 ) Weibull (19, 0. 35 ) 3 40 1200 60 0 Weibull (17, 2. 2 ) Weibull (20, 2. 1 ) Weibull (13, 0. 5 ) Weibull (23, 0. 40 ) a Weibull (λ, k); λ: Scale parameter, k: Shape parameter. Table 7 The effects of different actions over failure and repair characteristics of all component versions for subsystem 1. Version Action type l F C j l i V C jl i AF l 1 AF l 2 AR l 2 AR l3 j=1 Component-level 1 0.5 5 1.1 1 1 1 2 2 7.5 1.25 1 1 1 3 4 15.5 1.25 1.1 1 1 4 0 20 1.4 1.25 1 1 51021110.7 Subsystem-level 6 32 0 1 1 0.7 0.8 7 40 0 1 1 0.7 0.7 8 53 0 1 1 0.35 0.5 j=2 Component-level 1 0.5 5 1 1 1 1 2 2.5 7.5 1 1 1 1 3 4.5 15.5 1.25 1 1 1 4 0 20 1.4 1.1 1 1 51021110.6 Subsystem-level 6 32 0 1 1 1 0.8 7 40 0 1 1 0.8 0.6 8 53 0 1 1 0.5 0.4 j=3 Component-level 1 0.5 5 1 1 1 1 2 2 7.5 1 1 1 1 3 5 15.5 1 1 1 1 4 0 21 1 1.1 1 1 5 10 2.5 1 1 1 0.7 Subsystem-level 6 32 0 1 1 0.8 1 7 40 0 1 1 0.7 0.7 8 53 0 1 1 0.3 0.5 A. Attar et al. Reliability Engineering and System Safety 157 (2017) 177–191 185