1 Introduction
System reliability analysis and optimization are important to utilize available resources and part types efficiently and to develop the preferred or optimal system design architecture. The roots of the mathematical treatment of optimization problems can be traced even in Antiquity (Žilinskas and Zhigljavsky, 2016). In recent years, reliability and availability have expanded their influence in various industries and fields, thus these concepts serve as essential quality elements in many systems. Allocating redundant elements in the subsystems has been recognized as an effective means to meet the system reliability or availability requirement. The RAP problem is the problem of finding an optimal allocation of redundant components subject to a set of resource constraints (Caserta and Voß, 2015). It is one of the best-developed problems in reliability engineering studies. The RAP, which involves choosing appropriate elements and placing them redundantly to form an optimal system structure with high reliability and low cost, has received much attention in the literature. This essential problem has many applications in the real-world systems, such as production systems design, etc. The optimal reliability design aims to determine a system structure that achieves higher levels of reliability by exchanging the existing components with more reliable components or/and using redundant components.
Researchers have studied this problem from many different perspectives (Kuo and Wan, 2007; Yeh and Hsieh, 2011). The RAP is commonly considered in a multi-criteria decision-making (MCDM) environment, which has many applications in engineering problems (Keshavarz Ghorabaee, 2016; Keshavarz Ghorabaee et al., 2014). MCDM problems are generally divided into two classes: multi-objective decision-making (MODM) methods and multi-attribute decision-making (MADM) methods. MODM methods are generally used for dealing with the RAP. The traditional objectives for the RAP are maximizing the reliability and minimizing the cost of the system. Both of these objectives are increased by including more components. This trade-off requires the problem to be evaluated in the multi-objective context. In multi-objective problems, a set of non-dominated Pareto optimal solutions are obtained instead of a single optimal solution (Radziukynienė and Žilinskas, 2014). Data perception is frequently a complex problem, especially when data point to a complicated phenomenon described by many parameters, i.e. multidimensional data are analysed (Dzemyda et al., 2013). The estimation of intrinsic dimensionality of high-dimensional data still remains a challenging issue (Karbauskaitė and Dzemyda, 2015, 2016). Various methods are used to determine the Pareto optimal solutions of a multi-objective problem. The ε-constraint is one of the MODM methods that can be applied to find them (Soylu and Kapan Ulusoy, 2011). One of early program realizations of MODM methods is the system MOP (Dzemyda and Šaltenis, 1994).
System availability, a concept closely related to reliability, refers to the scale of measuring the reliability of a repairable system. The repairable system indicates a system that can be repaired to operate normally in the event of any failure. When the time-to-repair (TTR) is not negligible in relation to the operational time, system’s reliability is measured by its availability (Høyland and Rausand, 2004).
RAP is a challenging subject which has attracted the attention of many authors. Generally, in the RAP there are two strategies for using the redundant components: active and standby (Ardakan et al., 2015).
Generally, there are two types of the redundancy allocation problems. In the first type, we deal with discrete component choices with predefined characteristics (reliability, cost, weight, etc.). The aim of solving the problems in this type is choosing the components and the corresponding redundancy levels. In the second type, component reliability or a distribution parameter is treated as a design variable, and component cost is a known increasing function of component reliability (Coit and Liu, 2000). The majority of works in the first type have focused on the reliability of non-repairable systems. Using the metaheuristics is the most common approach for dealing with this problem because the RAP is an NP-hard problem and the computational time of optimal algorithms for the NP-hard problems is exponentially increased by enlarging the size of the problem (Amiri et al., 2014; Chern, 1992). Gupta et al. (2009) considered the problem of constrained redundancy allocation of the series system with interval valued reliability of components. They formulated the problem as an unconstrained integer programming problem with interval coefficients by penalty function technique and solved it by an advanced genetic algorithm (GA). Beji et al. (2010) proposed a hybrid algorithm based on particle swarm optimization and local search algorithm for the RAP. In addition, they introduced an adaptive penalty function for encouraging the algorithm to explore the feasible and near feasible region. Zhang et al. (2014) proposed a practical approach, combining bare-bones particle swarm optimization and sensitivity-based clustering for solving multi-objective reliability redundancy allocation problems. Wang and Li (2014) advanced a meta-heuristic approach called particle swarm optimization and applied it to obtain near-optimal solutions of the RAP in the multi-state systems with bridge topology. Keshavarz Ghorabaee et al. (2015) considered a RAP related to a system of s independent k-out-of-n subsystems in series. Maximization of the system reliability and minimization of the system cost were the objectives of the problem. They proposed four multi-objective genetic algorithms to handle this problem. It should be said that the current study is also related to the first type of the redundancy allocation problems.
Ardakan and Hamadani (2014) applied the mixed-integer non-linear optimization of reliability–redundancy allocation problem. Later, Abouei Ardakan et al. (2016) presented a new interpretation and formulation of the reliability–RAP and compared solution with traditional approaches.
Most of the studies for availability maximization have pertained to the second class of problems. There are few researches on availability of repairable systems that belong to the first type of problems. Castro and Cavalca (2003) presented an availability optimization problem of an engineering system assembled in a series configuration which has the redundancy of units and teams of maintenance as optimization parameters and developed a genetic algorithm to solve this problem. Zoulfaghari et al. (2014) studied a bi-objective RAP for a system with mixed repairable and non-repairable components. They proposed a new mixed integer nonlinear programming (MINLP) model to analyse this problem and used an efficient genetic algorithm to solve it. Lins and Droguett (2011) proposed a multi-objective genetic algorithm (GA) coupled with discrete event simulation to solve redundancy allocation problems in systems subjected to imperfect repairs. They validated the multi-objective GA via examples with analytical solutions and showed its superior performance when compared to a multi-objective ant colony algorithm (ACO).
The majority of these studies have considered systems that involve components with constant failure and repair rates. In these cases, time-to-failure (TTF) and time-to-repair (TTR) of the components have an exponential distribution. The hypotheses of constant failure rates are rarely met in real situations and may reduce the estimation accuracy of the entire system reliability or availability. In the classical system reliability theory, problems with such features are usually handled by a Markovian model. If we consider a system with components that have other distributions (e.g. Gamma, Normal, Weibull, etc.) for TTF and TTR, we cannot usually use an analytical method to obtain the system availability. In this study, we propose a new approach to design a system with the components that could have any distribution for their TTF and TTR.
The flexibility given by the simulation enables us to interpret many real aspects in problem modelling. Simulation is especially useful in situations where an analytical treatment is not utilizable (Lins and Droguett, 2009). A simulation model is a causal model of a real system. A meta-model is an approximation of the input/output (I/O) transformation that is implied by the simulation model. There are different types of meta-models such as polynomial regression models, splines, neural networks, etc. (Kleijnen and Sargent, 2000). Many researchers have developed meta-model-based approaches for optimization engineering problems. Aytug et al. (1996) proposed a method optimizing the number of kanbans in a pull production system by using simulation meta-modelling. Through meta-modelling, they determined the relationship between the number of kanbans and the average time to fill a customer order. Yang and Tseng (2002) introduced a new approach to optimize throughput and cycle time performance of integrated circuit ink-marking machines based on a simulation meta-model, a hybrid response surface method, and lexicographical goal programming approach. Noguera and Watson (2006) developed a general simulation meta-model in a particular company of chemicals industry in order to understand better how plant design parameters could be optimized in order to maximize plant throughput under certain environmental conditions and with certain asset investment constraints. Wang and Shan (2007) reviewed the meta-model-based techniques in supporting design optimization, including model approximation, design space exploration, problem formulation, and solving various types of optimization problems. Azadeh et al. (2010) presented an integrated approach based on simulation, meta-modelling using DOE and goal programming to solve the job shop scheduling problem with multiple objectives. Zakerifar et al. (2011) described the application of Kriging meta-modelling in multiple-objective simulation optimization. They utilized a simulation model of an $(s,S)$ inventory system to demonstrate the capabilities of Kriging meta-modelling as a simulation tool. Amiri and Mohtashami (2012) proposed a multi-objective formulation of the buffer allocation problem in unreliable production lines. They used a factorial design which has to build a meta-model for estimating production rate based on simulation. A genetic algorithm was utilized for solving the model and determining the optimal (or near optimal) size of each buffer storage.
In this research, we consider a bi-objective RAP in a series-parallel repairable system. Maximization of the system steady-state availability and minimization of the system cost are the objectives of the problem, and the system is constrained by a predefined weight. Although the TTF and TTR of the components could have any distribution in the proposed approach, we consider an example where the components’ TTF and TTR follow Gamma distribution. The Gamma distribution is considered because it is a flexible one and it can be used to model increasing, decreasing, and constant failure rates, similarly to the Weibull distribution. Moreover, it can be used to approximate several component failure time distributions (Amari, 2012). Also, this distribution has an important application for modelling the distribution of the TTF of a component, subjected to shocks whose arrivals follow a homogeneous Poisson process with intensity lambda. If the component is subjected to partial damage or degradation by each shock and fails completely at the k-th shock, the distribution of the time to failure of the component is given by the Gamma distribution (Modarres et al., 2009).
In the proposed approach, the steady-state availability of each subsystem in the series-parallel system is evaluated with an individual meta-model. The design of experiment, simulation, and the stepwise regression are used to build meta-models for calculating an approximation of the steady-state availability of each subsystem. In other words, steady-state availability of each subsystem is the response of the experiments and redundancy levels are the factors of them. Besides cost objective function, created meta-models are considered in the mathematical model using a max–min approach. The augmented ε-constraint method is utilized for obtaining the Pareto (near-Pareto) optimal solutions to this problem.
The rest of this paper is organized as follows. In Section 2, we show the general and max–min formulations of the RAP for the series-parallel systems. In Section 3, we describe the proposed approach in detail. In Section 4, we use an illustrative example to represent the performance of the proposed approach. The conclusions are discussed in Section 5.
2 Preliminaries
2.1 Redundancy Allocation Problem (RAP) Formulation
We consider the design of a system formed by a predefined number of subsystems in series (s), which may have several components in parallel as shown in Fig. 1. The objectives of this RAP are maximizing the system steady-state availability (A) and minimizing the system cost (C). There are ${m_{i}}$ $(i=1,2,\dots ,s)$ component types (choices) with different parameters for each subsystem. The TTF and TTR of each component could have any distribution such as the Gamma, Normal, Weibull, etc. Each subsystem may have a minimum number of allowed components $({n_{i,\min }})$. Moreover, the system is constrained by a maximum value of weight (${W_{\max }})$. Based on this problem, we can define two formulations: the general and the max–min formulation. We first present the general formulation in the following because of the vast application of it. However, the max–min formulation presented in the next subsection is used in the process of the proposed approach. The max–min formulation enables us to design the experiments of each subsystem individually, and this may lead to a reduction in the number of experiments.
2.1.1 The General Mathematical Formulation
The general mathematical modelling for the RAP in a repairable system can be shown as follows:
2.1.2 Max–Min Problem Formulation
A series system fails if any of its components fail; therefore, improving the availability of the least available component will readily impact on system availability. Although system availability is not directly maximized in this formulation, the resulting solutions lead to high system availability (Ramirez-Marquez et al., 2004).
2.1.3 Notation
| i | Subsystem index, |
| j | Component type (choice) index, |
| s | Number of subsystems, |
| ${A_{s}}$ | System steady-state availability, |
| ${C_{s}}$ | System cost, |
| ${W_{s}}$ | System weight, |
| ${A_{ij}}$ | Steady-state availability of jth component for subsystem i, |
| ${A_{i}}$ | Steady-state availability for subsystem i, |
| ${x_{ij}}$ | Quantity of jth component for subsystem i, |
| ${c_{ij}}$ | Cost of jth component for subsystem i, |
| ${w_{ij}}$ | Weight of jth component for subsystem i, |
| ${W_{\max }}$ | Maximum value of weight, |
| ${m_{i}}$ | Number of available components choice for subsystem i, |
| ${n_{i}}$ | Number of components used in subsystem i, |
| ${n_{i,\min }}$ | Lower bound for ${n_{i}}$. |
2.2 ε-Constraint Method
ε-constraint method is one of the most famous and efficient approaches in comparison with traditional weighting approaches. It has too many advantages in the case of having one main objective function in the multi-objective decision-making problems. There are two major kinds of this method, traditional and augmented. Although the augmented method is used in the process of the proposed approach, we also present the traditional method in this section to introduce the basics of the augmented method.
2.2.1 Traditional ε-Constraint Method
Consider a multi-objective optimization problem with maximization objectives:
where ${f_{i}}(X)$ is the ith objective function, p is the number of objective functions, S is the solution space, and X is the solution vector. With the traditional ε-constraint method, one objective is optimized and the other objectives are added into constraint space for guaranteeing that basic requirements are satisfied. As a result, the aforementioned model can be re-established as follows:
where ${\varepsilon _{1}},{\varepsilon _{2}},\dots ,{\varepsilon _{p-1}}$ are the satisfaction levels which stipulate the minimum requirements on the constrained objectives. Solutions can be obtained by parametrical variations of satisfaction levels ${\varepsilon _{1}},{\varepsilon _{2}},\dots ,{\varepsilon _{p-1}}$ in the right-hand side of constraints (Deb, 2001). Mavrotas (2009) showed that the results of this method are usually not efficient. In order to circumvent the defect of the traditional ε-constraint method, he proposed the augmented ε-constraint method.
(3)
\[ \begin{array}{l}\operatorname{Max}\big\{{f_{1}}(X),{f_{2}}(X),\dots ,{f_{p}}(X)\big\},\\ {} \text{s.t.}\hspace{2.5pt}X\in S,\end{array}\](4)
\[ \begin{array}{l}\operatorname{Max}{f_{p}}(X),\\ {} \text{s.t.}\hspace{2.5pt}{f_{1}}(X)\geqslant {\varepsilon _{1}},\\ {} {f_{2}}(X)\geqslant {\varepsilon _{2}},\\ {} \hspace{1em}\vdots \\ {} {f_{p-1}}(X)\geqslant {\varepsilon _{p-1}},\\ {} X\in S,\end{array}\]2.2.2 Augmented ε-Constraint Method
The augmented ε-constraint method transforms inequality constraints of constrained objectives into equality constraints by introducing non-negative slack variables or surplus variables and then augments the objective function with the weighted sum of these slack or surplus variables (Mavrotas, 2009). The above model can be reformulated as the following:
where δ is an adequately small number usually between ${10^{-3}}$ and ${10^{-6}}$ and ${r_{i}}$, $i\in [1,p-1]$, is the range of the ith objective.
(5)
\[ \begin{array}{l}\operatorname{Max}{f_{p}}(X)+\delta \bigg(\displaystyle \frac{{s_{1}}}{{r_{1}}}+\displaystyle \frac{{s_{2}}}{{r_{2}}}+\cdots +\displaystyle \frac{{s_{p-1}}}{{r_{p-1}}}\bigg),\\ {} \text{s.t.}\hspace{2.5pt}{f_{1}}(X)-{s_{1}}={\varepsilon _{1}},\\ {} {f_{2}}(X)-{s_{2}}={\varepsilon _{2}},\\ {} \hspace{1em}\vdots \\ {} {f_{p-1}}(X)-{s_{p-1}}={\varepsilon _{p-1}},\\ {} X\in S\hspace{1em}\text{and}\hspace{1em}{s_{i}}\in {R^{+}},\hspace{1em}i\in [1,p-1],\end{array}\]3 The Proposed Approach
In this section, we propose an approach to deal with RAP in a repairable series-parallel system with s subsystems. The framework of using the proposed approach is depicted in Fig. 2.
3.1 Design of Experiments (DOE)
The Design of Experiments (DOE) was developed in the early 1920s by Sir Ronald Fisher at the Rothamsted Agricultural Field Research Station in London, England. His initial experiments were concerned with determining the effect of various fertilizers on different plots of land (Roy, 2001). Since then DOE has been widely accepted and applied in biological and agricultural fields. Several successful applications of DOE have been presented by many US and European manufacturers over the last years (Antony, 2003).
In performing a design experiment, we intentionally make changes to the input process or variables (factors) in order to observe corresponding changes in the output (response). In this study, we consider the redundancy levels as the factors and the steady-state availability of each subsystem as the response of the experiments.
3.1.1 Estimate the Range of the Factors
For the design of experiments, we need the range of the factors. In this problem, we have not any information about lower and upper bound of the factors (${x_{ij}}$). Therefore, we must estimate these values for experimental design. To deal with this issue, we propose a heuristic way based on right-hand side (RHS) and technological coefficients of the weight constraint. Let ${x^{L}}$ and ${x^{U}}$ denote the approximated lower bound and upper bound of factors (${x_{ij}}$), respectively. The following equations are used for estimating ${x^{L}}$ and ${x^{U}}$:
where
In the experimental layout, ${x^{L}}$ and ${x^{U}}$ are replaced by −1 and 1 (coded values), respectively. It should be noted that the average of ${x^{L}}$ and ${x^{U}}$ is the central point and replaced by zero (midrange). In other words, the actual values of the factors are in the range of ${x^{L}}$ to ${x^{U}}$ whereas, the coded values of the factors are in the range of −1 to 1. Let ${x^{\prime }_{ij}}$ denote the coded value of ${x_{ij}}$ (actual value). The coded values and actual values of the factors can be converted to each other by the following equation:
(7)
\[ {x^{U}}=\left\{\begin{array}{l@{\hskip4.0pt}l}2\big\lfloor \frac{\omega }{2}\big\rfloor +1\hspace{1em}& \text{if}\hspace{2.5pt}{x^{L}}\hspace{2.5pt}\text{is odd},\\ {} 2\big\lfloor \frac{\omega }{2}+1\big\rfloor \hspace{1em}& \text{if}\hspace{2.5pt}{x^{L}}\hspace{2.5pt}\text{is even},\end{array}\right.\](8)
\[ \omega =\frac{{W_{\max }}}{\Big(\frac{{\textstyle\textstyle\sum _{i=1}^{s}}{\textstyle\textstyle\sum _{j=1}^{{m_{i}}}}{w_{ij}}}{{\textstyle\textstyle\sum _{i=1}^{s}}{m_{i}}}\Big){n_{i,\min }}s}.\]3.1.2 Central Composite Design (CCD)
The most popular response surface design is the central composite design. It is useful for building a second-order model for the response variable without needing to use a complete three-level factorial experiment (Leksakul and Limcharoen, 2014). It combines a two-level fractional factorial and two other kinds of points (Lin et al., 2012):
A central composite design is usually utilized for Badkar et al. (2012):
Face centred design is a type of central composite design. In this design, the axial points are at the centre of each face of the factorial space, and the design requires three levels of each factor (Ghoreishi, 2006). We use this type of central composite design for designing experiments in this study. The design of experiments for each subsystem is performed individually; therefore, the total number of experiments is easily reduced. According to this approach, we must build a meta-model for each subsystem.
-
– Centre points, for which all the factor values are at the zero (or midrange) value. This point is often replicated in order to improve the precision of the experiment;
-
– Axial (or star) points, for which all but one factor are set at zero (midrange) and that one factor is set at outer (axial) values.
3.2 Simulation
Simulation is the imitation of the operation of a real-world process or system over time (Banks et al., 2010). Simulation can be used to show the eventual real effects of alternative conditions and courses of action. Simulation is also used when the real system cannot be engaged, because it may not be accessible, or it may be dangerous or unacceptable to engage, or it is being designed but not yet built, or it may simply not exist (Sokolowski and Banks, 2011). The aim of the RAP is designing a system with high level of reliability or availability. In this paper, we consider a repairable system that TTF and TTR of components can follow any distribution (such as the Gamma, Normal, Weibull, etc.). In this situation, we cannot use an analytical method to obtain the steady-state availability of the system. Therefore, simulation is utilized for computing steady-state availability. We have a simulation response for each experiment in any subsystem design. These responses are used to build a meta-model for each subsystem.
3.3 Stepwise Regression to Build Meta-Models
Stepwise regression is a type of multiple linear regressions that can select the best-fitted combination of independent variables with forward-adding and backward-deleting variables (Ralston and Wilf, 1960). Stepwise regression is used when there is little theory to guide the selection of terms in a model, and the modeller wants to use whatever seems to provide a good fit (Draper and Smith, 1981). The stepping procedure begins as an initial model definition, with a stepped forward addition of a variable to the previous model. The critical F-value is then used to check the eligibility of the added variable. With a new variable added, the previous variables in the model may lose their predictive ability. Thus, stepping criteria are used to check the significance of all the included variables. If the variable is insignificant, then the backward method is used to delete it. Forward-adding and backward-deleting are repeated until no variable is added or removed. The stepping procedure is eliminated when the optimized model is established (Chen et al., 2013). We perform stepwise regression to build a meta-model for each subsystem. In this research, the considered p-value threshold for both forward and backward direction of stepwise regression is 0.25. Subsequently, ${\hat{A}_{i}}$ stands for the steady-state availability meta-model of the ith subsystem.
3.4 Reformulating the Problem and Obtaining the Solutions
The meta-models of each subsystem (${\hat{A}_{i}}$), the augmented ε-constraint method and the max–min formulation are used to reformulate the considered redundancy allocation problem. The new formulation of the problem is shown as follows:
where s is the surplus variable, ${r_{c}}$ is the range of ${C_{s}}$ (system cost) and ε is the satisfaction level. We have to vary the satisfaction level (ε) to find the Pareto (near-Pareto) optimal solutions. Following equation is utilized for varying ε value:
where ${C_{L}}$ is the lower bound of system cost (${C_{s}}$) and α is a real number that can be changed in the range of 0 to 1.
This mathematical formulation is not applicable to obtain the steady-state availability of the system because the meta-models of subsystems (${\hat{A}_{i}}$) are not limited to the range $[0,1]$. However, solving this mathematical model yields the configuration of the system and the redundancy levels that maximize the steady-state availability and minimize the cost of the system.
4 Illustrative Example
Table 1
Distribution parameters of components for given example.
| Subsystem $(i)$ | Choice 1 $(j=1)$ | Choice 2 $(j=2)$ | Choice 3 $(j=3)$ | Choice 4 $(j=4)$ | ||||||||
| ${\alpha _{ij}}$ | ${\beta _{ij}^{f}}$ | ${\beta _{ij}^{r}}$ | ${\alpha _{ij}}$ | ${\beta _{ij}^{f}}$ | ${\beta _{ij}^{r}}$ | ${\alpha _{ij}}$ | ${\beta _{ij}^{f}}$ | ${\beta _{ij}^{r}}$ | ${\alpha _{ij}}$ | ${\beta _{ij}^{f}}$ | ${\beta _{ij}^{r}}$ | |
| 1 | 2.2 | 0.00639 | 0.00173 | 1.3 | 0.01233 | 0.00323 | 1.7 | 0.00517 | 0.00074 | 1 | 0.0105 | 0.00508 |
| 2 | 1.9 | 0.00575 | 0.00151 | 1.1 | 0.01111 | 0.00156 | – | – | – | – | – | – |
| 3 | 1.13 | 0.01103 | 0.00222 | 1.35 | 0.01262 | 0.00281 | 1.75 | 0.00531 | 0.00135 | – | – | – |
| 4 | 1.47 | 0.01365 | 0.00318 | 1.55 | 0.01057 | 0.00249 | 1.68 | 0.00554 | 0.00131 | – | – | – |
| 5 | 2.3 | 0.00644 | 0.00163 | 1.12 | 0.01278 | 0.00296 | – | – | – | – | – | – |
Table 2
Costs and weights of components for given example.
| Subsystem $(i)$ | Choice 1 $(j=1)$ | Choice 2 $(j=2)$ | Choice 3 $(j=3)$ | Choice 4 $(j=4)$ | ||||
| ${c_{ij}}$ | ${w_{ij}}$ | ${c_{ij}}$ | ${w_{ij}}$ | ${c_{ij}}$ | ${w_{ij}}$ | ${c_{ij}}$ | ${w_{ij}}$ | |
| 1 | 100 | 11 | 79 | 7 | 58 | 7 | 90 | 5 |
| 2 | 93 | 25 | 95 | 10 | – | – | – | – |
| 3 | 65 | 8 | 42 | 9 | 69 | 10 | – | – |
| 4 | 41 | 9 | 28 | 10 | 53 | 13 | – | – |
| 5 | 92 | 8 | 83 | 7 | – | – | – | – |
We use an example of redundancy allocation problems in repairable systems to demonstrate the performance of the proposed approach. The example is a series-parallel system that is designed with 5 subsystems. For each subsystem there are two, three or four component choices (types). TTF and TTR of each component follow a 2-parameter Gamma distribution with shape parameter α and scale parameter β (the mean of distribution is $\alpha /\beta $). The distribution parameters of the components are provided in Table 1, and components’ cost and weight are provided in Table 2. We use an equivalent shape parameter for both TTF and TTR of any component. In Table 1, ${\beta ^{f}}$ denotes the scale parameter of TTF and ${\beta ^{r}}$ refers to the scale parameter of TTR. A minimum number of components within a subsystem i (${n_{i,\min }}$) is 3 for all subsystems. Maximization of the system steady-state availability and minimization of the system cost are the objectives of the example, and the maximum allowed system weight is 500.
4.1 DOE and Simulation
Table 3
Design of experiments and simulation responses for subsystem 1.
| Experiment number | Factors | Response | |||
| ${x^{\prime }_{11}}$ | ${x^{\prime }_{12}}$ | ${x^{\prime }_{13}}$ | ${x^{\prime }_{14}}$ | ${A_{1}}$ | |
| 1 | −1 | −1 | 1 | 1 | 0.8837 |
| 2 | 0 | 0 | 0 | −1 | 0.7978 |
| 3 | −1 | 1 | −1 | −1 | 0.7854 |
| 4 | −1 | 0 | 0 | 0 | 0.8298 |
| 5 | 1 | 1 | −1 | 1 | 0.9423 |
| 6 | 1 | 1 | −1 | −1 | 0.8694 |
| 7 | 0 | 0 | 0 | 0 | 0.8745 |
| 8 | −1 | 1 | −1 | 1 | 0.9009 |
| 9 | 1 | −1 | −1 | −1 | 0.7797 |
| 10 | 0 | 0 | 1 | 0 | 0.8904 |
| 11 | −1 | −1 | −1 | −1 | 0.6480 |
| 12 | −1 | 1 | 1 | −1 | 0.8253 |
| 13 | 0 | 0 | −1 | 0 | 0.8475 |
| 14 | 0 | 0 | 0 | 0 | 0.8797 |
| 15 | −1 | −1 | −1 | 1 | 0.8183 |
| 16 | −1 | −1 | 1 | −1 | 0.7222 |
| 17 | 1 | 0 | 0 | 0 | 0.9118 |
| 18 | 0 | 0 | 0 | 1 | 0.9090 |
| 19 | 1 | 1 | 1 | 1 | 0.9516 |
| 20 | 1 | −1 | 1 | −1 | 0.8067 |
| 21 | 1 | 1 | 1 | −1 | 0.8861 |
| 22 | −1 | 1 | 1 | 1 | 0.9291 |
| 23 | 0 | −1 | 0 | 0 | 0.8264 |
| 24 | 1 | −1 | −1 | 1 | 0.8973 |
| 25 | 0 | 1 | 0 | 0 | 0.8920 |
| 26 | 1 | −1 | 1 | 1 | 0.9129 |
As mentioned in previous sections, we have to estimate a lower and an upper bound for factors (${x_{ij}}$). The Eqs. (6), (7), and (8) are used for estimation of ${x^{L}}$ and ${x^{U}}$. The results are shown as follows:
According to Eq. (9), we can represent the relation between the coded variable and actual variable:
Next, a face centred design is used to design the experiments of each subsystem. The response of each experiment (steady-state availability of subsystem) is obtained by simulation. We use Enterprise Dynamics 7 (Falcon version) software to simulate each subsystem experiments. The detailed information of the simulation with this software is not presented because the simulation step of the proposed approach could be performed by any simulation software. The design of experiments (with coded factors) and the responses of experiments for each subsystem are shown in Tables 3 to 7. It should be noted that the responses shown in these tables are calculated by averaging the responses of simulation in 10 runs. As can be seen in these tables, using lower bound (${x^{\prime }_{ij}}=-1$) for all factors leads to the lowest response (steady-state availability of subsystem). On the other hand, the highest steady-state availability of subsystem is the result of using upper levels (${x^{\prime }_{ij}}=1$) for all factors. This relation between the levels of factors and the response shows the good performance of estimation, DOE, and simulation.
Table 4
Design of experiments and simulation responses for subsystem 2.
| Experiment number | Factors | Response | |
| ${x^{\prime }_{21}}$ | ${x^{\prime }_{22}}$ | ${A_{2}}$ | |
| 1 | −1 | −1 | 0.3074 |
| 2 | −1 | 0 | 0.4005 |
| 3 | 1 | 0 | 0.6304 |
| 4 | 1 | 1 | 0.6800 |
| 5 | 1 | −1 | 0.5672 |
| 6 | 0 | 0 | 0.5204 |
| 7 | 0 | −1 | 0.4499 |
| 8 | −1 | 1 | 0.4717 |
| 9 | 0 | 1 | 0.5858 |
| 10 | 0 | 0 | 0.5123 |
Table 5
Design of experiments and simulation responses for subsystem 3.
| Experiment number | Factors | Response | ||
| ${x^{\prime }_{31}}$ | ${x^{\prime }_{32}}$ | ${x^{\prime }_{33}}$ | ${A_{3}}$ | |
| 1 | 0 | 0 | −1 | 0.6370 |
| 2 | −1 | 1 | −1 | 0.6333 |
| 3 | 0 | 1 | 0 | 0.7655 |
| 4 | 1 | 1 | −1 | 0.7337 |
| 5 | 0 | 0 | 0 | 0.7146 |
| 6 | 1 | 0 | 0 | 0.7361 |
| 7 | 1 | −1 | −1 | 0.6301 |
| 8 | 1 | −1 | 1 | 0.7602 |
| 9 | −1 | −1 | −1 | 0.4479 |
| 10 | 0 | 0 | 0 | 0.6915 |
| 11 | 0 | 0 | 1 | 0.7606 |
| 12 | −1 | −1 | 1 | 0.6790 |
| 13 | 1 | 1 | 1 | 0.8395 |
| 14 | −1 | 1 | 1 | 0.7594 |
| 15 | −1 | 0 | 0 | 0.6459 |
| 16 | 0 | −1 | 0 | 0.6340 |
Table 6
Design of experiments and simulation responses for subsystem 4.
| Experiment number | Factors | Response | ||
| ${x^{\prime }_{41}}$ | ${x^{\prime }_{42}}$ | ${x^{\prime }_{43}}$ | ${A_{4}}$ | |
| 1 | 0 | 0 | 0 | 0.7220 |
| 2 | 1 | −1 | 1 | 0.7767 |
| 3 | −1 | −1 | −1 | 0.4761 |
| 4 | 0 | 0 | −1 | 0.6589 |
| 5 | −1 | 1 | −1 | 0.6761 |
| 6 | 0 | 1 | 0 | 0.7703 |
| 7 | 1 | 1 | −1 | 0.7811 |
| 8 | 0 | 0 | 0 | 0.6870 |
| 9 | 1 | 0 | 0 | 0.7641 |
| 10 | 1 | −1 | −1 | 0.6640 |
| 11 | 1 | 1 | 1 | 0.8458 |
| 12 | −1 | 1 | 1 | 0.7568 |
| 13 | 0 | 0 | 1 | 0.7747 |
| 14 | −1 | −1 | 1 | 0.6387 |
| 15 | −1 | 0 | 0 | 0.6608 |
| 16 | 0 | −1 | 0 | 0.6733 |
4.2 Stepwise Regression and Building Meta-Models
Table 8
Analysis of variance for fitted model of subsystem 1.
| Source | DF | Sum of squares | Mean square | F Ratio | p-value |
| Model | 12 | 0.12123 | 0.01010 | 230.184 | <0.0001 |
| Lack of fit | 12 | 0.00055 | 0.00004 | 3.4335 | 0.4007 |
| Error | 13 | 0.00057 | 0.00004 | ||
| C. Total | 25 | 0.12180 | |||
| RSquare = 0.9953 | |||||
Table 9
Analysis of variance for fitted model of subsystem 2.
| Source | DF | Sum of squares | Mean square | F Ratio | p-value |
| Model | 4 | 0.11045 | 0.02761 | 937.5039 | <0.0001 |
| Lack Of Fit | 4 | 0.00011 | 0.00002 | 0.8835 | 0.6527 |
| Error | 5 | 0.00014 | 0.00002 | ||
| C. Total | 9 | 0.11060 | |||
| RSquare = 0.9986 | |||||
Table 10
Analysis of variance for fitted model of subsystem 3.
| Source | DF | Sum of squares | Mean square | F Ratio | p-value |
| Model | 7 | 0.11907 | 0.01701 | 79.097 | <.0001 |
| Lack Of Fit | 7 | 0.00145 | 0.00020 | 0.7728 | 0.7073 |
| Error | 8 | 0.00172 | 0.00021 | ||
| C. Total | 15 | 0.12079 | |||
| RSquare = 0.9857 | |||||
Table 11
Analysis of variance for fitted model of subsystem 4.
| Source | DF | Sum of squares | Mean square | F Ratio | p-value |
| Model | 7 | 0.10900 | 0.01557 | 64.9566 | <0.0001 |
| Lack Of Fit | 7 | 0.00130 | 0.00018 | 0.3065 | 0.8862 |
| Error | 8 | 0.00191 | 0.00024 | ||
| C. Total | 15 | 0.11092 | |||
| RSquare = 0.9827 | |||||
Table 12
Analysis of variance for fitted model of subsystem 5.
| Source | DF | Sum of squares | Mean square | F Ratio | p-value |
| Model | 4 | 0.10925 | 0.02731 | 110.56 | <0.0001 |
| Lack Of Fit | 4 | 0.00087 | 0.00021 | 0.6042 | 0.7323 |
| Error | 5 | 0.00123 | 0.00024 | ||
| C. Total | 9 | 0.11049 | |||
| RSquare = 0. 9888 | |||||
According to results of the previous step, we can use the stepwise regression to build each subsystem’s meta-model. Analysis of variance for the fitted model of the response of subsystem 1 is shown in Table 8. If the p-value of a statistical model is smaller than significance level ($\alpha =0.05$), it can be concluded that the model is significant at the α level. We can see that the p-value of the statistical model of subsystem 1 has a very small value (less than 0.0001). It can be considered as evidence that the statistical model is significant. The p-value for the Lack of Fit test is greater than 0.05 and represents an insignificant lack of fit. RSquare estimates the proportion of variation in the response that can be attributed to the model rather than to random error. It can be seen in Table 8 that RSquare of the fitted model is close to 1 (RSquare = 0.9953). An RSquare closer to 1 indicates a better fit. Tables 9 to 12 show the analysis of variance for fitted models of subsystems 2 to 5, respectively. As can be seen in these tables, all the p-values of the statistical models are less than 0.0001, all the p-values of the Lack of Fit tests are greater than 0.05 and all the RSquares are greater than 0.95. With respect to these results, it can be said that all statistical models are significant and well fitted in the α significance level.
Table 13
The estimates of subsystem 1 meta-model.
| Term | Estimate | Std error | t Ratio | p-value |
| Intercept | 0.87134 | 0.00251 | 346.52 | <0.0001 |
| ${x^{\prime }_{11}}$ | 0.03417 | 0.00156 | 21.88 | <0.0001 |
| ${x^{\prime }_{12}}$ | 0.03816 | 0.00156 | 24.44 | <0.0001 |
| ${x^{\prime }_{13}}$ | 0.01773 | 0.00156 | 11.36 | <0.0001 |
| ${x^{\prime }_{14}}$ | 0.05692 | 0.00156 | 36.46 | <0.0001 |
| ${x^{\prime }_{11}}{x^{\prime }_{12}}$ | −0.00724 | 0.00165 | −4.37 | 0.0008 |
| ${x^{\prime }_{11}}{x^{\prime }_{13}}$ | −0.00870 | 0.00165 | −5.25 | 0.0002 |
| ${x^{\prime }_{12}}{x^{\prime }_{13}}$ | −0.00550 | 0.00165 | −3.32 | 0.0055 |
| ${x^{\prime }_{11}}{x^{\prime }_{14}}$ | −0.01181 | 0.00165 | −7.13 | <0.0001 |
| ${x^{\prime }_{12}}{x^{\prime }_{14}}$ | −0.01237 | 0.00165 | −7.47 | <0.0001 |
| ${x^{\prime }_{13}}{x^{\prime }_{14}}$ | −0.00246 | 0.00165 | −1.49 | 0.1609 |
| ${x^{\prime }_{12}}{x^{\prime }_{12}}$ | −0.00922 | 0.00365 | −2.52 | 0.0257 |
| ${x^{\prime }_{14}}{x^{\prime }_{14}}$ | −0.01503 | 0.00365 | −4.11 | 0.0012 |
Table 14
The estimates of subsystem 2 meta-model.
| Term | Estimate | Std error | t Ratio | p-value |
| Intercept | 0.51714 | 0.00271 | 190.57 | <.0001 |
| ${x^{\prime }_{21}}$ | 0.11635 | 0.00221 | 52.51 | <.0001 |
| ${x^{\prime }_{22}}$ | 0.06883 | 0.00221 | 31.07 | <.0001 |
| ${x^{\prime }_{21}}{x^{\prime }_{22}}$ | −0.01286 | 0.00271 | −4.74 | 0.0052 |
| ${x^{\prime }_{21}}{x^{\prime }_{21}}$ | −0.00757 | 0.00350 | −2.16 | 0.0831 |
Table 15
The estimates of subsystem 3 meta-model.
| Term | Estimate | Std error | t Ratio | p-value |
| Intercept | 0.70059 | 0.00598 | 117.02 | <0.0001 |
| ${x^{\prime }_{31}}$ | 0.05340 | 0.00463 | 11.52 | <0.0001 |
| ${x^{\prime }_{32}}$ | 0.05801 | 0.00463 | 12.51 | <0.0001 |
| ${x^{\prime }_{33}}$ | 0.07167 | 0.00463 | 15.46 | <0.0001 |
| ${x^{\prime }_{31}}{x^{\prime }_{32}}$ | −0.01035 | 0.00518 | −2 | 0.0809 |
| ${x^{\prime }_{31}}{x^{\prime }_{33}}$ | −0.01515 | 0.00518 | −2.92 | 0.0192 |
| ${x^{\prime }_{32}}{x^{\prime }_{33}}$ | −0.01617 | 0.00518 | −3.12 | 0.0143 |
| ${x^{\prime }_{31}}{x^{\prime }_{31}}$ | −0.01404 | 0.00757 | −1.85 | 0.1008 |
Table 16
The estimates of subsystem 4 meta-model.
| Term | Estimate | Std error | t Ratio | p-value |
| ${x^{\prime }_{41}}$ | 0.06232 | 0.00489 | 12.73 | <0.0001 |
| ${x^{\prime }_{42}}$ | 0.06014 | 0.00489 | 12.28 | <0.0001 |
| ${x^{\prime }_{43}}$ | 0.05364 | 0.00489 | 10.96 | <0.0001 |
| ${x^{\prime }_{41}}{x^{\prime }_{42}}$ | −0.01650 | 0.00547 | −3.01 | 0.0167 |
| ${x^{\prime }_{41}}{x^{\prime }_{43}}$ | −0.00823 | 0.00547 | −1.5 | 0.1713 |
| ${x^{\prime }_{42}}{x^{\prime }_{43}}$ | −0.01623 | 0.00547 | −2.97 | 0.018 |
| ${x^{\prime }_{41}}{x^{\prime }_{41}}$ | −0.01031 | 0.00799 | −1.29 | 0.2331 |
Table 17
The estimates of subsystem 5 meta-model.
| Term | Estimate | Std error | t Ratio | p-value |
| Intercept | 0.58234 | 0.00786 | 74.1 | <0.0001 |
| ${x^{\prime }_{51}}$ | 0.10056 | 0.00642 | 15.67 | <0.0001 |
| ${x^{\prime }_{52}}$ | 0.08653 | 0.00642 | 13.48 | <0.0001 |
| ${x^{\prime }_{51}}{x^{\prime }_{52}}$ | −0.02305 | 0.00786 | −2.93 | 0.0325 |
| ${x^{\prime }_{51}}{x^{\prime }_{51}}$ | −0.02531 | 0.01015 | −2.49 | 0.0548 |
The meta-model parameter estimations considering the main effects of factors, interactions between them and 2nd power of them, are shown in Tables 13 to 17 for subsystems 1 to 5, respectively. If the p-value of an estimate of a term (coefficient) is smaller than α level, we can say that this coefficient is significant in the considered significance level. As can be seen in these tables, most of the estimates of terms (coefficients) have p-values which are smaller than $\alpha =0.05$. Therefore, the obtained meta-models have significant coefficients in the considered α level. Moreover, the low values of standard errors (Std Error < 0.01) in these tables show the small confidence intervals of coefficients and high accuracy of estimates of terms.
With respect to Eq. (15) and Tables 13 to 17, we can represent the meta-model of each subsystem in actual (non-coded) variables as follows:
(16)
\[\begin{aligned}{}{\hat{A}_{1}}=& 0.87134+0.03417({x_{11}}-2)+0.03816({x_{12}}-2)+0.01773({x_{13}}-2)\\ {} & +0.05692({x_{14}}-2)-0.00724({x_{11}}-2)({x_{12}}-2)\\ {} & -0.00870({x_{11}}-2)({x_{13}}-2)-0.00550({x_{12}}-2)({x_{13}}-2)\\ {} & -0.01181({x_{11}}-2)({x_{14}}-2)-0.01237({x_{12}}-2)({x_{14}}-2)\\ {} & -0.00246({x_{13}}-2)({x_{14}}-2)-0.00922{({x_{12}}-2)^{2}}-0.01503{({x_{14}}-2)^{2}},\end{aligned}\](17)
\[\begin{aligned}{}{\hat{A}_{2}}=& 0.51714+0.11635({x_{21}}-2)+0.06883({x_{22}}-2)\\ {} & -0.01286({x_{21}}-2)({x_{22}}-2)-0.00757{({x_{21}}-2)^{2}},\end{aligned}\](18)
\[\begin{aligned}{}{\hat{A}_{3}}=& 0.70059+0.05340({x_{31}}-2)+0.05801({x_{32}}-2)+0.07167({x_{33}}-2)\\ {} & -0.01035({x_{31}}-2)({x_{32}}-2)-0.01515({x_{31}}-2)({x_{33}}-2)\\ {} & -0.01617({x_{32}}-2)({x_{33}}-2)-0.01404{({x_{31}}-2)^{2}},\end{aligned}\]4.3 Using ε-Constraint Method to Solve the Model
We utilize the formulation that has been presented based on the augmented ε-constraint method and the max–min approach for solving the problem. The mathematical formulation of the example is shown as follows:
where ${\hat{A}_{i}}$ is the meta-model of ith subsystem. We have ${r_{c}}=5022$, ${C_{L}}=912$ and $\delta ={10^{-5}}$. Therefore, the equation for varying ε is obtained as follows:
By varying the value of α we can alter ε in the specified range. The above mathematical model has been solved with different values of α. We have used the global solver of LINGO 14 software for solving this model. The solutions (configurations) that resulted from this approach and the corresponding system costs are presented in Table 18.
(21b)
\[\begin{array}{l}\displaystyle {C_{s}}=100{x_{11}}+79{x_{12}}+58{x_{13}}+90{x_{14}}+93{x_{21}}+95{x_{22}}+65{x_{31}}\\ {} \displaystyle \phantom{{C_{s}}=}+42{x_{32}}+69{x_{33}}+41{x_{41}}+28{x_{42}}+53{x_{43}}+92{x_{51}}+83{x_{52}},\end{array}\]Table 18
Solutions of the mathematical model.
| α | Actual variables | Cost of the system | |||||||||||||
| ${x_{11}}$ | ${x_{12}}$ | ${x_{13}}$ | ${x_{14}}$ | ${x_{21}}$ | ${x_{22}}$ | ${x_{31}}$ | ${x_{32}}$ | ${x_{33}}$ | ${x_{41}}$ | ${x_{42}}$ | ${x_{43}}$ | ${x_{51}}$ | ${x_{52}}$ | ||
| 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 3 | 0 | 0 | 3 | 0 | 0 | 3 | 912 |
| 0.02 | 0 | 1 | 2 | 0 | 3 | 0 | 0 | 4 | 0 | 0 | 4 | 0 | 0 | 3 | 1003 |
| 0.04 | 0 | 0 | 2 | 1 | 3 | 0 | 0 | 4 | 0 | 0 | 4 | 0 | 3 | 0 | 1041 |
| 0.06 | 0 | 0 | 2 | 1 | 4 | 0 | 0 | 4 | 0 | 1 | 3 | 0 | 0 | 4 | 1203 |
| 0.08 | 0 | 2 | 1 | 0 | 4 | 0 | 0 | 5 | 0 | 1 | 4 | 0 | 0 | 4 | 1283 |
| 0.1 | 0 | 2 | 1 | 0 | 4 | 0 | 0 | 5 | 0 | 0 | 6 | 0 | 4 | 0 | 1334 |
| 0.12 | 0 | 3 | 0 | 0 | 5 | 0 | 0 | 6 | 0 | 0 | 6 | 0 | 4 | 0 | 1490 |
| 0.14 | 0 | 3 | 0 | 0 | 5 | 0 | 0 | 6 | 0 | 0 | 7 | 0 | 5 | 0 | 1610 |
| 0.16 | 0 | 3 | 0 | 0 | 5 | 0 | 0 | 7 | 0 | 1 | 6 | 0 | 5 | 0 | 1665 |
| 0.18 | 0 | 0 | 0 | 3 | 5 | 1 | 0 | 7 | 0 | 2 | 5 | 0 | 5 | 0 | 1806 |
| 0.2 | 0 | 0 | 0 | 3 | 6 | 0 | 0 | 7 | 0 | 0 | 8 | 0 | 5 | 1 | 1889 |
| 0.22 | 0 | 0 | 0 | 4 | 6 | 0 | 0 | 7 | 0 | 0 | 8 | 0 | 6 | 0 | 1988 |
| 0.24 | 0 | 0 | 0 | 4 | 6 | 0 | 0 | 8 | 0 | 0 | 9 | 0 | 6 | 0 | 2058 |
| 0.26 | 0 | 0 | 0 | 4 | 7 | 0 | 0 | 8 | 0 | 1 | 8 | 0 | 6 | 0 | 2164 |
| 0.28 | 0 | 0 | 0 | 4 | 7 | 0 | 0 | 9 | 0 | 1 | 8 | 0 | 7 | 0 | 2298 |
| 0.30–0.32 | 0 | 0 | 0 | 5 | 7 | 0 | 0 | 9 | 0 | 1 | 9 | 0 | 7 | 0 | 2416 |
| 0.34 | 0 | 0 | 0 | 5 | 8 | 0 | 0 | 10 | 0 | 0 | 11 | 0 | 7 | 0 | 2566 |
| 0.36 | 0 | 0 | 1 | 5 | 8 | 0 | 0 | 10 | 0 | 0 | 11 | 0 | 8 | 0 | 2716 |
| 0.38 | 0 | 0 | 0 | 6 | 8 | 0 | 0 | 10 | 0 | 3 | 8 | 0 | 8 | 0 | 2787 |
| 0.4 | 0 | 0 | 0 | 6 | 8 | 0 | 0 | 0 | 7 | 0 | 12 | 0 | 8 | 0 | 2839 |
| 0.42 | 0 | 0 | 0 | 6 | 9 | 0 | 0 | 0 | 7 | 6 | 5 | 0 | 8 | 0 | 2982 |
| 0.44 | 0 | 0 | 0 | 7 | 9 | 0 | 0 | 0 | 7 | 6 | 5 | 0 | 8 | 0 | 3072 |
| 0.46–0.54 | 0 | 1 | 2 | 6 | 8 | 0 | 0 | 0 | 7 | 3 | 8 | 0 | 8 | 1 | 3128 |
| 0.56–0.60 | 0 | 7 | 0 | 0 | 0 | 13 | 0 | 0 | 7 | 1 | 0 | 8 | 0 | 9 | 3483 |
| 0.62–64 | 0 | 0 | 1 | 7 | 0 | 14 | 0 | 0 | 8 | 8 | 5 | 0 | 9 | 0 | 3866 |
| 0.66 | 0 | 0 | 0 | 7 | 0 | 15 | 0 | 0 | 8 | 1 | 0 | 9 | 0 | 10 | 3955 |
| 0.68 | 0 | 0 | 0 | 9 | 0 | 15 | 2 | 11 | 0 | 0 | 0 | 8 | 9 | 0 | 4079 |
| 0.70 | 0 | 10 | 0 | 0 | 0 | 15 | 0 | 0 | 9 | 0 | 0 | 8 | 10 | 0 | 4180 |
| 0.72 | 0 | 0 | 0 | 8 | 0 | 15 | 1 | 0 | 9 | 0 | 0 | 9 | 10 | 0 | 4228 |
| 0.74 | 0 | 0 | 0 | 8 | 0 | 15 | 1 | 0 | 9 | 0 | 0 | 9 | 0 | 13 | 4387 |
| 0.76 | 0 | 0 | 1 | 11 | 0 | 15 | 1 | 0 | 9 | 0 | 0 | 9 | 9 | 0 | 4464 |
| 0.78 | 1 | 1 | 1 | 10 | 0 | 15 | 0 | 0 | 9 | 12 | 0 | 0 | 0 | 11 | 4588 |
| 0.80–1 | 0 | 0 | 0 | 14 | 0 | 15 | 0 | 0 | 8 | 0 | 0 | 8 | 11 | 0 | 4673 |
Table 19
The average, standard deviation, minimum and maximum of steady-state availability of the system obtained by simulation according to each solution.
| α | Average | Standard deviation | Minimum | Maximum |
| 0 | 0.0157 | 0.0004 | 0.0150 | 0.0164 |
| 0.02 | 0.0374 | 0.0004 | 0.0367 | 0.0380 |
| 0.04 | 0.0393 | 0.0004 | 0.0386 | 0.0398 |
| 0.06 | 0.0525 | 0.0011 | 0.0508 | 0.0541 |
| 0.08 | 0.0629 | 0.0008 | 0.0618 | 0.0640 |
| 0.1 | 0.0745 | 0.0007 | 0.0733 | 0.0756 |
| 0.12 | 0.1043 | 0.0038 | 0.0982 | 0.1082 |
| 0.14 | 0.1237 | 0.0061 | 0.1133 | 0.1353 |
| 0.16 | 0.1387 | 0.0048 | 0.1299 | 0.1453 |
| 0.18 | 0.2074 | 0.0015 | 0.2055 | 0.2089 |
| 0.2 | 0.2393 | 0.0024 | 0.2344 | 0.2420 |
| 0.22 | 0.2726 | 0.0056 | 0.2612 | 0.2798 |
| 0.24 | 0.2872 | 0.0036 | 0.2827 | 0.2940 |
| 0.26 | 0.3322 | 0.0081 | 0.3216 | 0.3484 |
| 0.28 | 0.3807 | 0.0091 | 0.3676 | 0.3930 |
| 0.3–0.32 | 0.4128 | 0.0077 | 0.4038 | 0.4288 |
| 0.34 | 0.4504 | 0.0095 | 0.4319 | 0.4611 |
| 0.36 | 0.4833 | 0.0121 | 0.4658 | 0.5053 |
| 0.38 | 0.5016 | 0.0086 | 0.4838 | 0.5112 |
| 0.4 | 0.5059 | 0.0105 | 0.4896 | 0.5273 |
| 0.42 | 0.5107 | 0.0097 | 0.4939 | 0.5244 |
| 0.44 | 0.5302 | 0.0083 | 0.5203 | 0.5418 |
| 0.46–0.54 | 0.5870 | 0.0127 | 0.5708 | 0.6090 |
| 0.56–0.60 | 0.6343 | 0.0053 | 0.6279 | 0.6448 |
| 0.62–64 | 0.7928 | 0.0125 | 0.7719 | 0.8121 |
| 0.66 | 0.8166 | 0.0090 | 0.8005 | 0.8307 |
| 0.68 | 0.8457 | 0.0091 | 0.8298 | 0.8605 |
| 0.70 | 0.8683 | 0.0109 | 0.8507 | 0.8830 |
| 0.72 | 0.8852 | 0.0095 | 0.8704 | 0.8977 |
| 0.74 | 0.9003 | 0.0076 | 0.8862 | 0.9113 |
| 0.76 | 0.9126 | 0.0091 | 0.8983 | 0.9255 |
| 0.78 | 0.9337 | 0.0102 | 0.9145 | 0.9472 |
| 0.80–1 | 0.9481 | 0.0086 | 0.9355 | 0.9623 |
As previously mentioned, we cannot obtain the steady-state availability of the system by solving the mathematical model. For verifying the validity of solutions, we simulate each solution 10 times. The average steady-state availability for each solution that obtained from 10 runs of simulation is given in Table 19. To show the stability of simulation, the standard deviation, minimum, and maximum of each solution are also shown in this table. According to this table, the standard deviation of responses of all solutions is smaller than 0.05. Therefore, we can say that the results of our simulation have a good stability. With respect to the average values of these simulation responses and the cost to each solution (last column of Table 18), we can obtain an approximation of the Pareto (near-Pareto) optimal frontier as shown in Fig. 3. To represent the obtained frontier more accurately, the minimum and maximum values of steady-state availability of the system in each solution are shown, in addition to the average points. As can be seen in Fig. 3, there is a trade-off between the steady-state availability and the cost of the system and none of the obtained points (solutions) are dominated by the other points. This result shows the fact that the obtained frontier is close to Pareto optimal frontier and also represents the good performance of the proposed approach.
The decision-maker can choose one of the obtained solutions of this approach based on the Pareto (near-Pareto) optimal frontier and her/his subjective evaluations. It is obvious that desirable thing in solving a RAP is obtaining the optimal configuration of the system, but not necessarily the value of the steady-state availability of the system. However, using a logarithmic transformation of the response (steady-state availability) before building meta-models is suggested for the future researches. The use of a logarithmic transformation function can lead to building more accurate meta-models from responses, which are limited within the range $[0,1]$. Then a back-transformation function can be used after solving the model to obtain the approximate values of steady-state availability.
5 Conclusion
We have proposed a methodology for solving a bi-objective RAP in a series-parallel repairable system. Maximization of steady-state availability and minimization of cost has been assumed as the objectives of the problem. The TTF and TTR of the components of the system can follow any distribution (such as Gamma, Normal, Weibull, etc.). In the proposed methodology, DOE, simulation and the stepwise regression have been used to build meta-models for finding the approximations of steady-state availability of subsystems in the series-parallel system. Obtained meta-models together with the cost function of the system have been utilized for modelling the problem according to a max–min approach. We have developed an augmented ε-constraint method for reformulating and solving the model. An illustrative example has been explained to show the performance of the proposed approach. In the example, the TTF and TTR of components follow the Gamma distribution. All meta-models of this example are statistically significant and have acceptable (close to 1) values of RSquare. We can also see a trade-off between the steady-state availability and the cost of the system in the obtained Pareto (near-Pareto) optimal frontier. Regarding these two facts, we can say that the proposed approach has a good performance.
