A sequential two-stage approach based on variational Bayesian inference for reliability-redundancy allocation

Abstract

The reliability-related design has been a crucial chain in complex and reliability-critical engineering systems. It serves as a preventive countermeasure to catastrophic failures caused by uncertainties. To keep up with this rapidly developing field, this paper presents a novel metaheuristic, based on the variational Bayesian inference (VBI) to efficiently solve the optimal reliability design. Specifically, the reliability-redundancy allocation problem (RRAP). The proposed metaheuristic starts from a primary population, then leverages VBI to fully excavate the information of feasible individuals from the previous generation, in order to produce the next-generation population. This process is iterated until the solution converges to the optimal decision scheme. In addition, we set up an automatic stratification strategy, so that the new individuals can approach the optimal solution faster. Furthermore, we divide RRAP into a reliability optimization problem (ROP) and a redundancy allocation problem (RAP). This not only reduces the dimension of decision variables, but also speeds up the convergence. ROP and RAP are solved sequentially and iteratively until the preset stopping condition is satisfied. The case studies showcase that the proposed approach can obtain the optimal or near-optimal solution within a reasonable period.

Introduction

Nowadays, the advancement of high-tech industrial technology has increased the requirement of system functionalities. As feedback, engineering systems have experienced an ever-increasing level of sophistication. Moreover, the ever-harsher work environment together with the increasing awareness of high reliability, creates challenging optimal reliability design for complicated and expensive systems, forcing it to be a hot spot of intensive concern in both academia and industry.^1,2 System reliability can be enhanced by providing redundancy at the component level and/or increasing component reliabilities.^3–5 This, without any doubt, increases the system budget. Thus, it requires a smooth trade-off between budget constraint and system reliability to obtain a cost-effective design. The problem that maximizes the system reliability through redundancy or/and component reliability choice is dubbed as the reliability redundancy allocation problem. It can be further categorized into three streams according to the decision variables’ types, that is, the redundancy allocation problem (RAP, a pure integer nonlinear programming), reliability optimization problem (ROP, a continuous nonlinear programming), and reliability-redundancy allocation problem (RRAP, a mixed-integer nonlinear programming). They are non-deterministic polynomial hard (NP-hard) problems. The topic this paper is concerned with pertains to RRAP.

In the light of designers’ requirements, RRAP can be formulated either to maximize the system reliability under resource constraints or to minimize the total cost under the minimum demand on system reliability. For brevity, we only utilize the former to illustrate the proposed approach. Putting more clearly, the general mathematical formulation of RRAP is:

$$\begin{array}{cc}\multicolumn{1}{c}{\underset{(\mathit{x},\mathit{r})}{max}:}& R(\mathit{x},\mathit{r})\hfill \\ \multicolumn{1}{c}{s.t.}& g_i(\mathit{x},\mathit{r})\le b_i(i=1,2,\cdots ,m)\hfill \\ x_j\in N_{+},\begin{array}{c}\hfill \end{array}{x}_j^L\le x_j\le x_j^U\begin{array}{c}\hfill \end{array}(j=1,2,\cdots ,n)\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}}& r_j^L\le r_j\le r_j^U\begin{array}{c}\hfill \end{array}(j=1,2,\cdots ,n)\hfill \end{array}$$

(1)

where $R(\mathit{x},\mathit{r})$ is the system reliability with regard to the redundancy level $\mathit{x}$ and component reliability $\mathit{r}$ ; $g_i(\mathit{x},\mathit{r})$ for $i=1,2,\cdots ,m$ is the constraint function (e.g., procurement cost, physical weight, volume, etc.), $m$ is the number of constraints, and $b_i$ is the corresponding threshold; $\mathit{x}=\{x_1,x_2,\cdots ,x_n\}$ and $\mathit{r}=\{r_1,r_2,\cdots ,r_n\}$ are decision variables, $n$ is the number of subsystems; $x_i/r_i$ is the redundancy level/reliability choice of subsystem $i$ ; $x_j^L/r_j^L$ and $x_j^U/r_j^U$ are respectively the lower and upper bounds of $x_j/r_j$ for $j=1,2,\cdots ,n$ ; $N_{+}$ is the set of positive integers.

The solution of equation (1) includes two parts, that is, component (reliability) choices and redundancy levels. RRAP is classified as the hardest problem in the optimal reliability design community. This is because the involved decision variables are mixed-integer and the system reliability function is nonlinear, non-separable, and non-convex. Hence, most of the canonical mathematical methods have not been successful in finding an optimal or near-optimal solution for RRAP. In contrast, solving RRAP via metaheuristics has attracted increasing attention in recent years.

In the past few decades, the genetic algorithm (GA) has attracted the widest attention.^6,7 GA is a numerical simulation approach coming from mimicking the natural selection and genetic mechanism. GA converts the process solving the optimization problem into a process similar to the crossover and mutation of chromosome genes in biological evolution. When solving complex combinatorial optimization problems, GA can usually rapidly obtain good optimization results. The GA-based metaheuristics for solving RRAP, for example, Ardakan et al. ⁸ investigated RRAP of hybrid-redundancy systems, and a modified version of GA was developed to solve the given optimization problem. Kim and Kim ⁹ investigated RRAP considering the optimal redundancy strategy (either active or cold-standby), and applied the parallel GA to solve the formed optimization problem. Particle swarm optimization (PSO) is another prevalent approach for solving RRAP. PSO is a random search algorithm that was constructed by imitating the foraging behavior of birds.^10,11 PSO is initialized as a group of random particles (random solutions) within the feasible domain, and the optimal solution is found through subsequent iterations. In each iteration, the particles update themselves by tracking two optimal solutions, that is, the one found by the particle itself (individual best or personal best) and the other one found by the whole population (global best). For the PSO-based metaheuristics to solve RRAP, Wu et al. ¹² proposed an improved PSO algorithm involving two position updating strategies to solve RRAP, in which a mutation operator was introduced after position updating. This not only prevents the algorithm from trapping into the local optimum, but also enhances its space developing ability. Ouyang et al. ¹³ introduced heterogeneous components into RRAP, and the formulated optimization problem was resolved by an improved PSO algorithm with stochastic perturbation nature.

Artificial bee colony algorithm (ABC) is an optimization method proposed via imitating the behavior of bees. It does not need to learn the special information of the problem, but only compare the advantages and disadvantages of the obtained information. Through the local optimization behavior of the individual artificial bee, the global optimum is finally emerged among these local optimums. The application of ABC algorithm for solving RRAP, see Refs.^14–16Cuckoo search algorithm (CS), originated from imitating the brood parasitism of cuckoo, is a population-based intelligent optimization approach. It mainly depends on two strategies, namely, the brood parasitism of cuckoo and the Levy flight mechanism. ¹⁷ The main advantages of CS algorithm include: (1) small number of parameters, (2) easy implementation, (3) random search path, and (4) strong optimization ability. The applications of CS in RRAP, see Refs.^18–20 The other metaheuristics for solving RRAP are, for example, the artificial fish swarm algorithm, ²¹ biogeography optimization approach, ²² and stochastic fractal search algorithm.^23–25 In addition, many studies have devoted their vigour to hybrid metaheuristics, in hope of obtaining more efficient algorithms than the single metaheuristic, for example, the combination method of PSO and simplified swarm optimization, ²⁶ the hybrid Jaya algorithm based on a time-variant acceleration coefficient, ²⁷ and the bilevel evolutionary algorithm based on quadratic approximations. ²⁸

Although fruitful achievements have been made in RRAP, obviously, the research on RRAP will not stop here. It remains an exciting topic with many new research challenges and opportunities. Motivated by this fact, a novel optimization algorithm based on the variational Bayesian inference (VBI) is proposed, aiming to search for the optimal decision scheme more quickly and globally. VBI is one of the variational methods in statistical inference, which can asymptotically obtain the posterior distribution of latent variables from a given variational family.^29–31 VBI leverages the mean field theory (MFT) to expand the posterior joint distribution of latent variables into the product of the marginal distributions of latent variables, in order to obtain the calculation framework. It can be utilized as a low computational alternative to the Markov chain Monte Carlo (MCMC), and has been applied to some machine learning algorithms, such as the learning of variational autoencoder.

The proposed optimization approach starts from an initial population randomly generated from the design region. Then, the next-generation population is produced by VBI, through leveraging the information provided by feasible individuals of the previous-generation population. This process is iterated until the optimal decision scheme is found. It is worth mentioning that the proposed algorithm does not utilize the information of all feasible individuals in each generation to produce new individuals, but first layers feasible individuals according to their objective function values by setting an intermediate threshold, and carries out the optimization algorithm along the direction of the reduction of the objective function value (for the minimization problem), in order to search for the optimal solution faster. The proposed optimization approach is a numerical-simulation-based technique, which is suitable for optimization problems with discontinuous, discrete, and multiple feasible regions. To apply the proposed optimization algorithm to solve RRAP, we first divide RRAP into ROP together with RAP, and solve these two formulated problems sequentially and iteratively until the stopping condition is satisfied. Such segmentation strategy transforms the original $2n$ -dimensional optimization problem into two $n$ -dimensional problems, thus reduces the complexity of the optimization problem. When $n$ is large, the proposed algorithm is more efficient to obtain the optimal or near-optimal solution.

This paper is structured as follows. The problem formulation is briefly expounded in Section 2. Section 3 concretely presents the proposed optimization algorithm for solving RRAP. The performance of the proposed approach is validated by several examples in Section 4. Conclusions are drawn in Section 5.

Problem formulation

Conventionally, we first transform the RRAP in equation (1) (a maximization problem) into a minimization problem:

$$\begin{array}{cc}\multicolumn{1}{c}{\underset{\mathit{z}}{min}:}& f(\mathit{z})=-R(\mathit{z})\hfill \\ \multicolumn{1}{c}{s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}}& g_i(\mathit{z})\le b_i\begin{array}{c}\hfill \end{array}(i=1,2,\cdots ,m)\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}}& z_j^L\le z_j\le z_j^U\begin{array}{c}\hfill \end{array}(j=1,2,\cdots ,n)\hfill \end{array}$$

(2)

where $\mathit{z}=\{\mathit{x},\mathit{r}\}$ is the set of decision variables; $z_j^L$ and $z_j^U$ are the lower and upper bounds of $z_j$ respectively, for $j=1,2,\cdots ,n$ ; $f(\mathit{z})=-R(\mathit{z})$ denotes the new objective function.

We then suppose that the decision variables in $\mathit{z}$ are all continuous (real variables). Let $\{{\mathit{z}}_1,{\mathit{z}}_2,\cdots ,{\mathit{z}}_N\}$ be $N$ initial individuals of decision variables $\mathit{z}$ , which are randomly generated from the design space. The infeasible individuals are filtered out by utilizing the given constraints $g_i(\mathit{z})\le b_i$ ( $i=1,2,\cdots ,m$ ). The feasible individuals are retained, denoted as ${\mathbf{S}}_f=\{{\mathit{z}}_1,{\mathit{z}}_2,\cdots ,{\mathit{z}}_{N_f}\}$ where $N_f$ is the number of feasible individuals.

To search for the optimal solution of design variables, the intuitive measure is to produce the next-generation better feasible individuals based on the information provided by current feasible individuals, and repeat this process to mine the feasible domain until the optimal solution comes to light. This is a general procedure of the greedy optimization algorithm. Undoubtedly, the quality of the next-generation population severely affects the convergence speed of the optimization process and the precision of the final solution. In other words, the strategy of producing the next-generation population is really crucial. The most classic strategy to generate the next-generation population should be the MCMC, but MCMC has two fatal defects: (1) the generation of the latter individual is reliant upon previous individuals, resulting in the generated individuals being related rather than independent; (2) the decision of acceptance or rejection of new individuals relies on the proposal distribution and a uniformly distributed variable. MCMC is a random approach, making it difficult to converge to a stable state. Hence, this paper adopts VBI instead of MCMC to produce the next-generation population. VBI is a deterministic approach, which is more stable than MCMC.

Proposed approach for solving RRAP

To facilitate understanding, we first list the general execution process of the proposed “crude” optimization algorithm. After showcasing the key techniques of the proposed optimization approach, we design an “elaborate” optimization algorithm for exploring the optimal decision scheme of RRAP.

Proposed non-directional optimization approach

This section first presents the general whole algorithm of the proposed “crude” optimization algorithm, as shown in Algorithm 1. This “crude” optimization algorithm is referred to as the non-directional approach, and the corresponding reason will be concretely explained later. After that, the key factors of the proposed approach are revealed, in order to pave the way for the subsequent design of an improved optimization algorithm.

Algorithm 1. Non-directional approach for solving the optimal reliability design
I. Initialization
1.	Let $k=0$ be the number of iterations, and ${\mathbf{S}}_f$ be an empty matrix.
II. First-generation population
2.	Generate the primary population ${\mathbf{S}}^{(k)}=\{{\mathit{z}}_1^{(k)},{\mathit{z}}_2^{(k)},\cdots ,{\mathit{z}}_N^{(k)}\}$ in the design domain where $N$ stands for the number of individuals.
3.	Screen the feasible individuals ${\mathbf{S}}_f^{(k)}=\{{\mathit{z}}_1^{(k)},{\mathit{z}}_2^{(k)},\cdots ,{\mathit{z}}_{N_f^{(k)}}^{(k)}\}$ in ${\mathbf{S}}^{(k)}$ by the given constraints, where $N_f^{(k)}$ is the number of feasible individuals.
4.	Let ${\mathbf{S}}_f={\mathbf{S}}_f\bigcup {\mathbf{S}}_f^{(k)}$ .
5.	Compute the objective function values $\{f({\mathit{z}}_1^{(k)}),f({\mathit{z}}_2^{(k)}),\cdots ,f({\mathit{z}}_{N_f^{(k)}}^{(k)})\}$ of feasible individuals in ${\mathbf{S}}_f^{(k)}$ .
6.	Obtain the current optimal objective function value $f^{*(k)}=min\{f({\mathit{z}}_1^{(k)}),f({\mathit{z}}_2^{(k)}),\cdots ,f({\mathit{z}}_{N_f^{(k)}}^{(k)})\}$ and the corresponding optimal solution ${\mathit{z}}^{*(k)}=\arg min\{f({\mathit{z}}_1^{(k)}),f({\mathit{z}}_2^{(k)}),\cdots ,f({\mathit{z}}_{N_f^{(k)}}^{(k)})\}$ .
III. New-generation population
	While (1)
7.	Produce the new-generation individuals ${\mathbf{S}}^{(k+1)}=\{{\mathit{z}}_1^{(k+1)},{\mathit{z}}_2^{(k+1)},\cdots ,{\mathit{z}}_N^{(k+1)}\}$ .
8.	Sift out the feasible individuals ${\mathbf{S}}_f^{(k+1)}=\{{\mathit{z}}_1^{(k+1)},{\mathit{z}}_2^{(k+1)},\cdots ,{\mathit{z}}_{N_f^{(k+1)}}^{(k+1)}\}$ from ${\mathbf{S}}^{(k+1)}$ .
9.	Let ${\mathbf{S}}_f={\mathbf{S}}_f\bigcup {\mathbf{S}}_f^{(k+1)}$ .
10.	Compute the objective function values $\{f({\mathit{z}}_1^{(k+1)}),f({\mathit{z}}_2^{(k+1)}),\cdots ,f({\mathit{z}}_{N_f^{(k+1)}}^{(k+1)})\}$ for individuals in ${\mathbf{S}}_f^{(k+1)}$ .
11.	Obtain the current optimal objective function value $f^{*(k+1)}=min\{f^{*(k)},f({\mathit{z}}_1^{(k+1)}),f({\mathit{z}}_2^{(k+1)}),\cdots ,f({\mathit{z}}_{N_f^{(k+1)}}^{(k+1)})\}$ and the corresponding optimal solution ${\mathit{z}}^{*(k+1)}=\arg min\{f^{*(k)},f({\mathit{z}}_1^{(k+1)}),f({\mathit{z}}_2^{(k+1)}),\cdots ,f({\mathit{z}}_{N_f^{(k+1)}}^{(k+1)})\}$ .
12.	If $|f^{*(k)}-f^{*(k+1)}|/|f^{*(k+1)}|\le {\epsilon }_1$ (where ${\epsilon }_1$ is a small positive real number, e.g., ${10}^{-3}$ )  Break End if
13.	Let $k=k+1$ .
	End while
IV. Output the optimal decision scheme
14.	Output the optimal solution ${\mathit{z}}^{*}={\mathit{z}}^{*(k+1)}$ and the corresponding objective function value $f^{*}=f^{*(k+1)}$ .

From Algorithm 1, we can see that the key of the proposed approach is how to establish a generator to produce the next-generation population based on the information about the existing feasible individuals. This is to increase the number of feasible individuals and procure individuals that are better than those of previous generations. To address this problem, the first idea that comes to mind is to fit a conditional probability density function (PDF) $P(\mathit{z}|\mathbb{F})$ of the decision variables $\mathit{z}$ , where $\mathbb{F}$ represents the feasible domain, and generate new individuals following this conditional PDF. Theoretically, the distribution $P(\mathit{z}|\mathbb{F})$ can be expressed as:

$$P(\mathit{z}|\mathbb{F})=\frac{I_{\mathbb{F}}(\mathit{z})P(\mathit{z})}{{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{I}_{\mathbb{F}}(\mathit{z})P(\mathit{z})\mathrm{d}\mathit{z}}$$

(3)

where $I_{\mathbb{F}}(\mathit{z})$ is the indicator function of the feasible domain $\mathbb{F}$ , such that $I_{\mathbb{F}}(\mathit{z})=1$ if $\mathit{z}\in \mathbb{F}$ , $I_{\mathbb{F}}(\mathit{z})=0$ otherwise. $P(\mathit{z})$ is the prior distribution of decision variables $\mathit{z}$ , which can simply be assumed to be uniform in $[{\mathit{z}}^L,{\mathit{z}}^U]$ where ${\mathit{z}}^L=\{z_1^L,z_2^L,\cdots ,z_n^L\}$ and ${\mathit{z}}^U=\{z_1^U,z_2^U,\cdots ,z_n^U\}$ are the lower bound and upper bound vectors of $\mathit{z}=\{z_1,z_2,\cdots ,z_n\}$ , respectively. The denominator ${\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{I}_{\mathbb{F}}(\mathit{z})P(\mathit{z})\mathrm{d}\mathit{z}$ is a normalized coefficient such that $P(\mathit{z}|\mathbb{F})$ is a PDF, that is, ${\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}P(\mathit{z}|\mathbb{F})\mathrm{d}\mathit{z}=1$ .

However, the conditional PDF $P(\mathit{z}|\mathbb{F})$ in equation (3) is not available in practice, since the feasible region $\mathbb{F}$ cannot be straightforwardly obtained in advance, that is, $I_{\mathbb{F}}(\mathit{z})$ cannot be precisely known. Therefore, we can only utilize the information of the feasible region that has been mined (i.e., ${\mathbf{S}}_f^{(k)}$ ) to establish current conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ . Undoubtedly, we can only construct a PDF as close as possible to the real conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ , and produce new individuals that are approximately following $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ . In conclusion, the current problem is how to establish an optimal asymptotical model $Q_{\mathit{z}}(\mathit{z})$ of the conditional distribution $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ so that the individuals generated from $Q_{\mathit{z}}(\mathit{z})$ are similar to those generated from $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ . In this paper, this problem is solved by VBI (see Section 3.2), that is, VBI is utilized as a generator integrating the information of current feasible individuals to produce the next-generation population.

The process of the proposed non-directional approach for solving the optimal reliability design is depicted in Figure 1, where MCS denotes the Monte Carlo simulation. In Figure 1, MCS is first used to generate the primary-generation individuals (blue points) in the design space. These individuals are divided into infeasible individuals (circular points) that do not meet the constraints and feasible individuals (triangular points) that satisfy the constraints. By calculating the objective function values of these feasible individuals, we can obtain the current minimum of the objective function ( $f^{*(1)}$ ) as well as the current optimal design ( ${\mathit{z}}^{*(1)}$ ). Then, based on the information provided by current feasible individuals, the new-generation individuals (green points) are produced via VBI. Similarly, we can choose feasible individuals in terms of the given constraints. This will expand the feasible domain. We can find a better solution ( ${\mathit{z}}^{*(2)}$ ) that achieves smaller objective function value ( $f^{*(2)}$ ), compared with the optimal solution ( ${\mathit{z}}^{*(1)}$ ) obtained in the last iteration. In the same vein, the third and fourth-generation individuals (red points and pink points) can be generated by VBI based on the acquired information of the existing feasible individuals. We can then screen out the feasible individuals and find the current optimal individuals. This process continues till the algorithm converges. For the case in Figure 1, the optimization algorithm achieves the convergence condition through an initialization step and three iterations. The obtained optimal solutions of these four steps are ${\mathit{z}}^{*(1)}\to {\mathit{z}}^{*(2)}\to {\mathit{z}}^{*(3)}\to {\mathit{z}}^{*(4)}$ , and the corresponding objective function values are $f^{*(1)}\to f^{*(2)}\to f^{*(3)}\to f^{*(4)}$ . We can draw the conclusion that, with the increased number of iterations, the number of feasible individuals increase, and the mining of feasible regions is more thorough. As a result, the obtained solution becomes better (the objective function value becomes smaller).

Figure 1.

Illustration of the VBI-based non-directional approach for searching for the optimal design.

Variational Bayesian inference

This section first expounds the basic theory and general algorithm of VBI. Then, some details are illustrated using an example.

Basic theory of VBI

The KL-divergence of the asymptotical model $Q_{\mathit{z}}(\mathit{z})$ and conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ is:

$$\mathrm{KL}(Q_{\mathit{z}}||P)={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log \left(\frac{Q_{\mathit{z}}(\mathit{z})}{P(\mathit{z}|{\mathbf{S}}_f^{(k)})}\right)\mathrm{d}\mathit{z}$$

(4)

Then, the optimal asymptotical model $Q_{\mathit{z}}(\mathit{z})$ for the conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ can be obtained by minimizing equation (4).

Furthermore, we convert equation (4) into:

$$\begin{array}{cc}\multicolumn{1}{c}{}& \mathrm{KL}(Q_{\mathit{z}}||P)=-{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log \left(\frac{P(\mathit{z}|{\mathbf{S}}_f^{(k)})}{Q_{\mathit{z}}(\mathit{z})}\right)\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& =-{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log \left(\frac{P(\mathit{z},{\mathbf{S}}_f^{(k)})}{Q_{\mathit{z}}(\mathit{z})P({\mathbf{S}}_f^{(k)})}\right)\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})(\log [Q_{\mathit{z}}(\mathit{z})]+\log [P({\mathbf{S}}_f^{(k)})])\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& -{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})]\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}=\log [P({\mathbf{S}}_f^{(k)})]+{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [Q_{\mathit{z}}(\mathit{z})]\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& -{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})]\mathrm{d}\mathit{z}\hfill \end{array}$$

(5)

where $P({\mathbf{S}}_f^{(k)})$ is the PDF of ${\mathbf{S}}_f^{(k)}$ , while $P(\mathit{z},{\mathbf{S}}_f^{(k)})$ denotes the joint PDF of $\mathit{z}$ and ${\mathbf{S}}_f^{(k)}$ .

Let $L(Q_{\mathit{z}})={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})]\mathrm{d}\mathit{z}-{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [Q_{\mathit{z}}(\mathit{z})]\mathrm{d}\mathit{z}$ , then equation (5) can be reformulated as:

$$\mathrm{KL}(Q_{\mathit{z}}||P)+L(Q_{\mathit{z}})=\log [P({\mathbf{S}}_f^{(k)})]$$

(6)

Since $\log [P({\mathbf{S}}_f^{(k)})]$ is independent of $\mathit{z}$ , it can be viewed as a constant. To minimize $\mathrm{KL}(Q_{\mathit{z}}||P)$ , we just need to maximize $L(Q_{\mathit{z}})$ . Furthermore, since $\mathrm{KL}(Q_{\mathit{z}}||P)\ge 0$ , we can obtain $\log [P({\mathbf{S}}_f^{(k)})]\ge L(Q_{\mathit{z}})$ . Thus, $L(Q_{\mathit{z}})$ can be viewed as the lower bound of $\log [P({\mathbf{S}}_f^{(k)})]$ , referred to as the evidence lower bound (ELBO). In conclusion, we can obtain the optimal asymptotical model $Q_{\mathit{z}}(\mathit{z})$ through maximizing the ELBO of $\log [P({\mathbf{S}}_f^{(k)})]$ , that is:

$$\begin{array}{cc}\multicolumn{1}{c}{}& max:L(Q_{\mathit{z}})={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})]\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& -{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [Q_{\mathit{z}}(\mathit{z})]\mathrm{d}\mathit{z}\hfill \end{array}$$

(7)

Following the MFT, we can obtain $Q_{\mathit{z}}(\mathit{z})=\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)$ where $Q_i(z_i)$ is the marginal distribution of $z_i$ . Then, ${\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [Q_{\mathit{z}}(\mathit{z})]\mathrm{d}\mathit{z}$ can be reformulated as:

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [Q_{\mathit{z}}(\mathit{z})]\mathrm{d}\mathit{z}={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)\log (\underset{j=1}{\overset{n}{\Pi }}{Q}_j(z_j))\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& ={\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)\sum _{j=1}^n\log [Q_j(z_j)]\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^n{\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)\log [Q_j(z_j)]\mathrm{d}\mathit{z}\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^n({\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j(z_j)]\mathrm{d}{z}_j)({\int }_{{\mathit{z}}_{-j}^L}^{{\mathit{z}}_{-j}^U}\underset{i=1,i\ne j}{\overset{n}{\Pi }}{Q}_i(z_i)\mathrm{d}{\mathit{z}}_{-j})\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^n({\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j(z_j)]\mathrm{d}{z}_j)\hfill \end{array}$$

(8)

where ${\mathit{z}}_{-j}=\{z_1,z_2,\cdots ,z_{j-1},z_{j+1},\cdots ,z_n\}$ .

In addition, ${\int }_{{\mathit{z}}^L}^{{\mathit{z}}^U}{Q}_{\mathit{z}}(\mathit{z})\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})]\mathrm{d}\mathit{z}$ can be reorganized as:

$$\begin{array}{l}{\int }_{z^L}^{z^U}{Q}_z(z)\log [P(z,S_f^{(k)})]dz={\int }_{z^L}^{z^U}\prod _{i=1}^n{Q}_i(z_i)\log [P(z,S_f^{(k)})]dz\\ ={\int }_{z^L}^{z^U}{Q}_j(z_j)\left(\prod _{i=1,i\ne j}^n{Q}_i(z_i)\log [P(z,S_f^{(k)})]d{z}_{-j}\right)d{z}_j\\ ={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j){\mathbb{E}}_{i\ne j}\left(\log [P(z,S_f^{(k)})]\right)d{z}_j={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [e^{{\mathbb{E}}_{i\ne j}\left(\log [P(z,S_f^{(k)})]\right)}]d{z}_j\\ ={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log \left[\frac{e^{{\mathbb{E}}_{i\ne j}\left(\log [P(z,S_f^{(k)})]\right)}}{{\int }_{z_j^L}^{z_j^U}{e}^{{\mathbb{E}}_{i\ne j}\left(\log [P(z,S_f^{(k)})]\right)}d{z}_j}\right]d{z}_j+{\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log \left[{\int }_{z_j^L}^{z_j^U}{e}^{{\mathbb{E}}_{i\ne j}\left(\log [P(z,S_f^{(k)})]\right)}d{z}_j\right]d{z}_j\\ ={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j^{*}(z_j)]d{z}_j+\xi \end{array}$$

(9)

where $Q_j^{*}(z_j)=\frac{e^{{\mathbb{E}}_{i\ne j}(\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})])}}{{\int }_{z_j^L}^{z_j^U}{e}^{{\mathbb{E}}_{i\ne j}(\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})])}\mathrm{d}{z}_j}$ , and $\xi ={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log \left[{\int }_{z_j^L}^{z_j^U}{e}^{{\mathbb{E}}_{i\ne j}(\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})])}\mathrm{d}{z}_j\right]\mathrm{d}{z}_j$ is a constant.

Combining equations (8) and (9), the ELBO $L(Q_{\mathit{z}})$ can be re-expressed as:

$$\begin{array}{cc}\multicolumn{1}{c}{L(Q_{\mathit{z}})}& ={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j^{*}(z_j)]\mathrm{d}{z}_j\hfill \\ \multicolumn{1}{c}{}& +\xi -\sum _{j=1}^n({\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j(z_j)]\mathrm{d}{z}_j)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}=({\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j^{*}(z_j)]\mathrm{d}{z}_j-{\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log [Q_j(z_j)]\mathrm{d}{z}_j)\hfill \\ \multicolumn{1}{c}{}& -\sum _{i=1,i\ne j}^n({\int }_{z_i^L}^{z_i^U}{Q}_i(z_i)\log [Q_i(z_i)]\mathrm{d}{z}_i)+\xi \hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}={\int }_{z_j^L}^{z_j^U}{Q}_j(z_j)\log \left[\frac{Q_j^{*}(z_j)}{Q_j(z_j)}\right]\mathrm{d}{z}_j\hfill \\ \multicolumn{1}{c}{}& -\sum _{i=1,i\ne j}^n({\int }_{z_i^L}^{z_i^U}{Q}_i(z_i)\log [Q_i(z_i)]\mathrm{d}{z}_i)+\xi \hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}=-\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))+\underset{i=1,i\ne j}{\overset{n}{\Pi }}H(Q_i(z_i))+\xi \hfill \end{array}$$

(10)

where $H(Q_i(z_i))=-{\int }_{z_i^L}^{z_i^U}{Q}_i(z_i)\log [Q_i(z_i)]\mathrm{d}{z}_i$ stands for the information entropy.

Since $\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))\ge 0$ and $H(Q_i(z_i))\ge 0$ , to maximize $L(Q_{\mathit{z}})$ we just need to minimize $\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))$ , that is, VBI transforms the minimization of $\mathrm{KL}(Q_{\mathit{z}}||P)$ into the following minimization problems:

$$\{\begin{array}{c}min:\mathrm{KL}(Q_1(z_1)||Q_1^{*}(z_1))\hfill \\ \multicolumn{1}{c}{min:\mathrm{KL}(Q_2(z_2)||Q_2^{*}(z_2))}\\ \multicolumn{1}{c}{\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}⋮\hfill \end{array}}\\ \multicolumn{1}{c}{min:\mathrm{KL}(Q_n(z_n)||Q_n^{*}(z_n))}\end{array}$$

(11)

Because $\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))\ge 0$ , to minimize $\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))$ we can set $\mathrm{KL}(Q_j(z_j)||Q_j^{*}(z_j))=0$ . Then, we can obtain:

$$Q_j(z_j)=Q_j^{*}(z_j)=\frac{e^{{\mathbb{E}}_{i\ne j}(\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})])}}{{\int }_{z_j^L}^{z_j^U}{e}^{{\mathbb{E}}_{i\ne j}(\log [P(\mathit{z},{\mathbf{S}}_f^{(k)})])}\mathrm{d}{z}_j}(j=1,2,\cdots ,n)$$

(12)

Finally, the asymptotical model $Q_{\mathit{z}}(\mathit{z})$ for conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ can be obtained by $Q_{\mathit{z}}(\mathit{z})=\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)$ . The general process of VBI constructing the asymptotical model $Q_{\mathit{z}}(\mathit{z})$ for conditional PDF $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ is as Algorithm 2.

Algorithm 2. Construct the asymptotical model of conditional PDF via VBI
1.	Allocate a prior distribution for $P(\mathit{z}|{\mathbf{S}}_f^{(k)})$ .
2.	Construct the joint PDF $P(\mathit{z},{\mathbf{S}}_f^{(k)})$ by the prior distribution and $P({\mathbf{S}}_f^{(k)})$ .
3.	Determine the form of $Q_{\mathit{z}}(\mathit{z})$ based on the form of $P(\mathit{z},{\mathbf{S}}_f^{(k)})$ .
4.	Obtain $Q_{\mathit{z}}(\mathit{z})=\underset{i=1}{\overset{n}{\Pi }}{Q}_i(z_i)$ by solving $Q_i(z_i)$ $(i=1,2,\cdots ,n)$ .

Case study

We take the rosenbrock function (Banana function) as an example to demonstrate the process of VBI approximating a PDF. The 2D banana function is as follows:

$$f(\mathit{z})=(\eta -z_1{)}^2+\delta (z_2-z_1^2{)}^2$$

(13)

where $\eta$ and $\delta$ are two parameters.

It is seen that the global minimum of $f(\mathit{z})$ is achieved at $\mathit{z}=(\eta ,{\eta }^2)$ , and the corresponding minimum is $f(\mathit{z})=0$ . Without the loss of generality, we set $\eta =1$ and $\delta =100$ . Then, the banana function of interest is $f(\mathit{z})=(1-z_1{)}^2+100(z_2-z_1^2{)}^2$ . The contour map of this banana function is shown in Figure 2, where the global minimum is located within a long, narrow, and parabolic valley. It is not difficult to find the valley. However, it is difficult to reach the global minimum.

Figure 2.

Contour map of the 2D banana function.

Now, we need to explore this special valley to find the minimum of the 2D banana function. We suppose that the information of some samples located in or around this valley has been learned. Let ${\mathbf{S}}_{\mathit{z}}=\{{\mathit{z}}_1,{\mathit{z}}_2,\cdots ,{\mathit{z}}_{N_1}\}$ be the set consisting of these learned samples where $N_1$ denotes the size of ${\mathbf{S}}_{\mathit{z}}$ . We want to fit an approximate model $Q_{\mathit{z}}(\mathit{z})$ of conditional PDF $P(\mathit{z}|{\mathbf{S}}_{\mathit{z}})$ in order to mine the information of this valley.

First of all, we need to allocate a prior distribution to the conditional PDF $P(\mathit{z}|{\mathbf{S}}_{\mathit{z}})$ . Actually, both uniform and normal prior distributions are popular choices. The statistical parameters (e.g., mean, variance, etc.) of a prior distribution could be tuned according to the problem at hand, in order to achieve a good performance when faced with different target distributions. Here, we consider six scenarios, as shown in Table 1, where $\mathbf{\mu }$ and $\mathbf{\sigma }$ represent the mean vector and standard deviation vector respectively. Then, the joint PDF $P(\mathit{z},{\mathbf{S}}_{\mathit{z}})$ is the product of the prior distribution and $P({\mathbf{S}}_{\mathit{z}})$ . In addition, the form of $Q_{\mathit{z}}(\mathit{z})$ can be determined in terms of the form of $P(\mathit{z},{\mathbf{S}}_{\mathit{z}})$ .

After that, the optimal asymptotical model $Q_{\mathit{z}}(\mathit{z})$ of $P(\mathit{z}|{\mathbf{S}}_{\mathit{z}})$ can be obtained by equation (11). As mentioned, $Q_{\mathit{z}}(\mathit{z})$ can be expressed as $Q_{\mathit{z}}(\mathit{z})=Q_1(z_1)Q_2(z_2)$ . Moreover, $Q_1(z_1)$ and $Q_2(z_2)$ can be obtained by equation (12). It is seen that we need to calculate two univariate integrals to obtain the marginal PDFs $Q_1(z_1)$ and $Q_2(z_2)$ , and this can be done by numerical integration approaches or other methods. This paper adopts the Bayesian quadrature to complete this task. Finally, the optimal asymptotical model $Q_{\mathit{z}}(\mathit{z})$ is obtained. The marginal PDFs $Q_1(z_1)$ and $Q_2(z_2)$ under these listed scenarios are depicted in Figure 3. Meanwhile, the approximate valley zones obtained in terms of these PDFs are also demonstrated in Figure 3.

Figure 3.

The approximate valley zones of the 2D banana function: (a) Scenario 1, (b) Scenario 2, (c) Scenario 3, (d) Scenario 4, (e) Scenario 5, and (f) Scenario 6.

It can be concluded from Figure 3 that, different prior distributions will obtain different posterior distributions, causing the approximate valley zones to be different from each other. Moreover, even if the form of the prior distribution is same, the parameter settings also have an impact on the posterior distribution. The ELBO estimates and corresponding standard deviations (Sd. (ELBO)) and the estimated minimum values (Min) of the banana function are presented in Table 1. It should be pointed out that, the uncertainty of the ELBO estimate comes from the Bayesian quadrature. It is seen that the minimum of the banana function is related to the prior distribution. For this case study, the normal prior with mean values {0.5, 0.5} and standard deviations {5, 5} (i.e., Scenario 2) is the best choice. Furthermore, we consider a special situation, as follows: if the obtained minimum is smaller than ${10}^{-5}$ , its value will be set to 0, then the minimum values under Scenarios 1, 2, 4, 5, 6 can be regarded as same. That is, under this particular hypothesis, the final result is immunized to the prior distribution if the allocated prior is not outrageous. Of course, if the allocated prior distribution is in a bad state, such as the Scenario 3, the final result will deviate far from the optimal solution.

Table 1.

Results of the case study under different prior distributions.

Scenario	Prior		ELBO	Sd. (ELBO)	Min
	Distribution	Parameters
1	Normal	$\mathbf{\mu }=\{0,0\}$ , $\mathbf{\sigma }=\{3,3\}$	−2.2755	$5.5423\times {10}^{-4}$	$8.4327\times {10}^{-9}$
2	Normal	$\mathbf{\mu }=\{0.5,0.5\}$ , $\mathbf{\sigma }=\{5,5\}$	−3.0263	$5.5024\times {10}^{-4}$	$2.2808\times {10}^{-11}$
3	Normal	$\mathbf{\mu }=\{1,1\}$ , $\mathbf{\sigma }=\{0.1,0.1\}$	−0.0471	$5.5002\times {10}^{-4}$	0.2147
4	Uniform	$[-1,1{]}^2$	0.7814	$3.9277\times {10}^{-4}$	$1.6864\times {10}^{-7}$
5	Uniform	$[-2,2{]}^2$	−0.5295	$4.6333\times {10}^{-4}$	$1.2384\times {10}^{-6}$
6	Uniform	$[-3,3{]}^2$	−1.2750	$4.7326\times {10}^{-4}$	$1.6230\times {10}^{-8}$

Proposed directional optimization approach

Until now, we have presented the general process and key techniques of the proposed optimization approach for solving optimal reliability design. As shown in Figure 1, the first-generation individuals are generated randomly in the design region, then the feasible individuals are filtered. On the premise of the information provided by current feasible individuals, the next-generation individuals are produced by VBI, so the explored feasible areas are expanded and the number of feasible individuals are increased. The optimal solution can be procured through repeating this process until the termination condition is met. However, in this approach, although the number of feasible individuals increase iteratively, these feasible individuals do not necessarily increase along the direction of the reduction of the objective function value. In other words, the next-generation population is produced based only on the information provided by current feasible individuals. This increases the number of feasible individuals disorderly and aimlessly, and the optimal individual of the new-generation population may not be as good as that of the previous-generation population. It can be seen from Figure 1 that the feasible individuals of the next generation are not necessarily better than those of the previous generations (the objective function values of the new feasible individuals are not necessarily smaller than those of the previous feasible individuals). This fact may incur in slow convergence of the optimization algorithm. For the above reasons, this approach is dubbed as a non-directional method.

General whole algorithm of the proposed directional optimization approach

In response to the above issues, we will make an improvement to the non-directional approach, so that the improved approach could produce more feasible individuals along the direction of reducing the value of the objective function to speed up the convergence of the algorithm. The basic measure of the improved approach is to set a series of intermediate events, forcing new individuals to be produced along the direction that the objective function value is smaller than the threshold of the intermediate event. Then, the new-generation individuals will fall (as far as possible) in the area whose objective function value is smaller than the preset threshold. These intermediate events are automatically determined instead of being set in advance, so the improved approach is very easy to be implemented without any additional constraints. The specific execution steps of the proposed improved optimization algorithm are listed in Algorithm 3.

Algorithm 3. Directional approach for solving the optimal reliability design
I. Initialization
1.	Let $k=0$ be the number of iterations, and set a probability $p_0$ satisfying $0<p_0<1$ .
II. First-generation population
2.	Generate the primary population ${\mathbf{S}}^{(k)}=\{{\mathit{z}}_1^{(k)},{\mathit{z}}_2^{(k)},\cdots ,{\mathit{z}}_N^{(k)}\}$ randomly in the design domain.
3.	Screen out the infeasible individuals, and retain the feasible individuals ${\mathbf{S}}_f^{(k)}=\{{\mathit{z}}_1^{(k)},{\mathit{z}}_2^{(k)},\cdots ,{\mathit{z}}_{N_f^{(k)}}^{(k)}\}$ in ${\mathbf{S}}^{(k)}$ by the given constraints.
4.	Calculate the objective function values of feasible individuals in ${\mathbf{S}}_f^{(k)}$ , and arrange the objective function values in descending order, that is, $f({\mathit{z}}_{[1]}^{(k)})\ge f({\mathit{z}}_{[2]}^{(k)})\ge \cdots \ge f({\mathit{z}}_{[N_f^{(k)}]}^{(k)})$ .
5.	Obtain the current optimal objective function value $f^{*(k)}=f({\mathit{z}}_{[N_f^{(k)}]}^{(k)})$ and optimal solution ${\mathit{z}}^{*(k)}={\mathit{z}}_{[N_f^{(k)}]}^{(k)}$ .
III. New-generation population
	While (1)
6.	Let the $(1-p_0)N_f^{(k)}$ th objective function value ${\zeta }_k=f({\mathit{z}}_{[(1-p_0)N_f^{(k)}]}^{(k)})$ be the threshold of the event $\mathbf{D}\mathbf{=}\{\mathit{z}|f(\mathit{z})\le {\zeta }_k,\mathit{z}\in {\mathbf{S}}_f^{(k)}\}$ .
7.	Use VBI to obtain the asymptotical model $Q_{\mathit{z}}(\mathit{z})$ of conditional distribution $P(\mathit{z}|\mathbf{D})$ .
8.	Produce the new population ${\mathbf{S}}^{(k+1)}=\{{\mathit{z}}_1^{(k+1)},{\mathit{z}}_2^{(k+1)},\cdots ,{\mathit{z}}_N^{(k+1)}\}$ by $Q_{\mathit{z}}(\mathit{z})$ .
9.	Screen the feasible individuals ${\mathbf{S}}_f^{(k+1)}=\{{\mathit{z}}_1^{(k+1)},{\mathit{z}}_2^{(k+1)},\cdots ,{\mathit{z}}_{N_f^{(k+1)}}^{(k+1)}\}$ in ${\mathbf{S}}^{(k+1)}$ .
10.	Calculate the objective function values of feasible solutions in ${\mathbf{S}}_f^{(k+1)}$ , and arrange the objective function values in descending order, that is, $f({\mathit{z}}_{[1]}^{(k+1)})\ge f({\mathit{z}}_{[2]}^{(k+1)})\ge \cdots \ge f({\mathit{z}}_{[N_f^{(k+1)}]}^{(k+1)})$ .
11.	Obtain the current optimal objective function value $f^{*(k+1)}=f({\mathit{z}}_{[N_f^{(k+1)}]}^{(k+1)})$ and optimal solution ${\mathit{z}}^{*(k+1)}={\mathit{z}}_{[N_f^{(k+1)}]}^{(k+1)}$ .
12.	If $|f^{*(k)}-f^{*(k+1)}|/|f^{*(k+1)}|\le {\epsilon }_1$   Break End if
13.	Let $k=k+1$ .
	End while
IV. Output the optimal decision scheme
14.	Output the optimal solution ${\mathit{z}}^{*}={\mathit{z}}^{*(k+1)}$ and corresponding objective function value $f^{*}=f^{*(k+1)}$ .

Comparing Algorithm 1 and Algorithm 3, we can observe that the improved algorithm has a filtering operation (Step 6). In this step, a fraction $(1-p_0)$ of worse individuals are abandoned while the fraction $p_0$ of better individuals are kept. The conceptual illustration of the proposed directional optimization approach is demonstrated in Figure 4. It is seen that the first step of the directional approach is the same as that of the non-directional approach, from which we get the first optimal solution ${\mathit{z}}^{*(1)}$ and its objective function value $f^{*(1)}$ . Now we are going to produce the second-generation individuals. First, we divide the current feasible individuals into two groups based on their objective function values, namely, the group whose objective function value is larger than ${\zeta }_1$ and the other group whose objective function value is smaller than ${\zeta }_1$ . Only the second-group feasible individuals are retained, that is, $\mathbf{D}\mathbf{=}\{\mathit{z}|f(\mathit{z})\le {\zeta }_1,\mathit{z}\in {\mathbf{S}}_f^{(1)}\}$ . Based on the information provided by $\mathbf{D}$ , VBI is applied to generate the new-generation individuals (green points). Pick the feasible individuals (green triangle points) and compute their objective function values. We can see that the objective function values of new feasible individuals are smaller than ${\zeta }_1$ . This fact makes a distinction between the non-directional and directional approaches. After that, the updated optimal solution ${\mathit{z}}^{*(2)}$ and its objective function value $f^{*(2)}$ are obtained. In the same vein, we can acquire the third and fourth-generation individuals, and update the optimal solution ${\mathit{z}}^{*(2)}\to {\mathit{z}}^{*(3)}\to {\mathit{z}}^{*(4)}$ and objective function value $f^{*(2)}\to f^{*(3)}\to f^{*(4)}$ . Finally, it is seen that the newly produced feasible individuals tend to be distributed along the direction of reduction of the objective function value. Moreover, in each generation, the newly produced individuals possess objective function values smaller than the current threshold ${\zeta }_k$ . This implies that the convergence speed of the directional approach is generally higher than that of the non-directional approach.

Figure 4.

Illustration of the VBI-based directional approach for searching for the optimal design.

Mapping function for integer decision variables

The proposed directional approach for solving the optimal reliability design involving real decision variables has been concretely presented. Here, we assume that the optimization problem only contains integer decision variables, that is, $\mathit{z}$ is a set of integer design variables. The strategy of handling integer decision variables is as follows. First, the integer decision variables $\mathit{z}$ are treated as real variables in the design domain $[{\mathit{z}}^L,{\mathit{z}}^U]$ . Then, we utilize the proposed approach to generate the population of real decision variables $\mathit{z}$ . After that, a mapping function is leveraged to map new individuals from the real number space into the integer space. For simplicity, the rounding approach is applied to round the new population to generate new individuals that satisfy the integer constraint.

Until now, the proposed optimization algorithm has been introduced in detail. To facilitate understanding, the flowchart of the proposed optimization algorithm for solving the optimal reliability design is depicted in Figure 5.

Figure 5.

Flowchart of the proposed optimization algorithm.

Proposed two-stage approach for solving RRAP

As mentioned, the solution of RRAP contains two parts, namely, the reliability choices as well as the redundancy levels. Of course, the proposed optimization approach can simultaneously handle the integer and real design variables, but such a hybrid measure has two defects: the computational effort of handling the integer and real decision variables simultaneously is huge. This is because the sum of dimensions of these two types of design variables is high, leading to the longer process time in producing new individuals; moreover, this high dimension also reduces the accuracy of VBI, resulting in the produced new individuals to be sub-optimal.

In response to this fact, we propose a two-stage approach to efficiently and precisely solve RRAP. The first stage is to solve ROP, in order to determine the reliability of components when the redundancy level is fixed. The second stage is to solve RAP, aiming to determine the redundancy level of subsystems when the reliability choice is known. These two stages are conducted sequentially and iteratively until the termination criterion is satisfied. The concrete algorithm of the proposed two-stage approach for solving RRAP is presented in Algorithm 4.

Algorithm 4. Proposed two-stage approach for solving RRAP
I. Initialization
1. Set the initial design solution ${\mathit{z}}^{(0)}=({\mathit{x}}^{(0)},{\mathit{r}}^{(0)})$ , and let $k=0$ be the number of iterations.
II. Optimization
While (1)
2. Let $k=k+1$ .
3. Solve the following ROP by the proposed directional approach to obtain the optimal solution ${\mathit{r}}^{(k)}$ .              $\begin{array}{cc}\multicolumn{1}{c}{}& \underset{\mathit{r}}{min}:f(\mathit{r})=-R({\mathit{x}}^{(k-1)},\mathit{r})\hfill \\ \multicolumn{1}{c}{}& s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{g}_i({\mathit{x}}^{(k-1)},\mathit{r})\le b_i\begin{array}{c}\hfill \end{array}(i=1,2,\cdots ,m)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{r}_j^L\le r_j\le r_j^U\begin{array}{c}\hfill \end{array}(j=1,2,\cdots ,n)\hfill \end{array}$
4. Solve the following RAP by the proposed directional approach to obtain the optimal solution ${\mathit{x}}^{(k)}$ .             $\begin{array}{cc}\multicolumn{1}{c}{}& \underset{\mathit{x}}{min}:f(\mathit{x})=-R(\mathit{x},{\mathit{r}}^{(k)})\hfill \\ \multicolumn{1}{c}{}& s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{g}_i(\mathit{x},{\mathit{r}}^{(k)})\le b_i\begin{array}{c}\hfill \end{array}(i=1,2,\cdots ,m)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{x}_j\in N_{+},\begin{array}{c}\hfill \end{array}{x}_j^L\le x_j\le x_j^U\begin{array}{c}\hfill \end{array}(j=1,2,\cdots ,n)\hfill \end{array}$
5. Update the optimal solution ${\mathit{z}}^{(k)}=({\mathit{x}}^{(k)},{\mathit{r}}^{(k)})$ .
6. Check the stopping condition.   If $|f({\mathit{z}}^{(k)})-f({\mathit{z}}^{(k-1)})|/|f({\mathit{z}}^{(k)})|<{\epsilon }_2$ (where ${\epsilon }_2\to 0^{+}$ is a real number)    Break   End if
End while
III. Output
7. Output the optimal decision scheme ${\mathit{z}}_{\mathrm{opt}}=({\mathit{x}}^{(k)},{\mathit{r}}^{(k)})$ and the system reliability $R_{\mathrm{opt}}=-f({\mathit{z}}_{\mathrm{opt}})$ .

Discussion and remarks

• (1) The proposed optimization algorithm starts from the string set of individuals instead of a single initial point. Thus, it is not easy for it to be caught into the local optimum, and is conducive to global optimum.

• (2) The proposed optimization algorithm processes multiple individuals in a population, simultaneously. This not only reduces the risk of falling into the local optimum, but also facilitates parallelization.

• (3) The proposed optimization algorithm does not need the information of the design space or any auxiliary information in advance, but only utilizes the given constraints to screen feasible individuals, and uses the objective function to evaluate individuals to guide its search direction. The objective function or constraint function is not constrained by any assumption or premise. These characteristics greatly expand the scope of application of the proposed optimization algorithm.

• (4) The proposed optimization algorithm generates a new population based on a group of excellent individuals rather than the most prominent individual. This avoids the emergence of super individual and eventually circumvents the local optimum.

• (5) In the proposed optimization algorithm, each generation only retains the $p_0$ proportion of better individuals and discards the $(1-p_0)$ proportion of poor individuals. Generally, the larger the $p_0$ , the more are the iterations required and consequently, the slower is the convergence speed. But this can comprehensively explore the feasible region and retain more excellent individuals to produce new individuals. The smaller $p_0$ means that fewer iterations are needed and the convergence speed is faster. However, this may lead to incomplete exploration of the feasible region, producing inferior individuals, and eventually leading to the explored optimal solution far away from the real optimal one.

• (6) To get more excellent individuals at each generation, we do not require that the individuals strictly satisfy the constraints, but retain those individuals that meet the constraints and those that are close to the constraint boundary. This aims to avoid missing the individuals that contribute greatly to the accuracy of the optimization algorithm. Therefore, the “feasible” individuals of the $k$ th generation can be modified as:

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mathbf{S}}_f^{(k)}=\{{\mathit{z}}^{(k)}|g_i({\mathit{z}}^{(k)})\le b_i,\forall i=1,2,\cdots ,m\}\bigcup \hfill \\ \multicolumn{1}{c}{}& \{{\mathit{z}}^{(k)}||g_i({\mathit{z}}^{(k)})-b_i|\le {\gamma }_i,\forall i=1,2,\cdots ,m\}\hfill \end{array}$$

(14)

where ${\gamma }_i\to 0^{+}$ is a preset small real number. The first term denotes that the feasible individuals solidly satisfy the constraints, and the second term stands for the individuals that softly satisfy the constraints, that is, the individuals distributed in the vicinity of the constraint boundary.

• (7) The proposed two-stage approach converts the solution of RRAP into a process of sequentially solving ROP and RAP. In each iteration, only the reliability choice or redundancy level needs to be optimized. This reduces the dimension of design variables that need to be optimized and raises the optimization speed.

• (8) This paper uses the uniform prior distribution in VBI, because the decision variables $\mathit{z}$ are generally regarded as uniform variables in the design space $[{\mathit{z}}^L,{\mathit{z}}^U]$ . The boundary of uniform prior distribution can be determined by the existing sample information. In addition, if we have acquired some information of a variable in advance, for example, its unconditional distribution is normal. Then setting the prior distribution to normal could usually get better results.

Examples

To showcase the performance of the proposed optimization approach for solving RRAP, three examples, including a series system with 10 subsystems, an oil transportation system, and a 15-unit system, are presented in this section. The CPU type which this paper used is “11th Gen Intel(R) Core(TM) i7-11800H @ 2.30 GHz.”

A series system with 10 subsystems

A series system with 10 subsystems, modified from Refs.,^12,17 is first used to validate the feasibility of the proposed algorithm. The optimization problem is to determine the reliability choice and redundancy level of each subsystem, as follows:

$$\underset{(\mathit{x},\mathit{r})}{max}:\begin{array}{c}\hfill \end{array}R(\mathit{x},\mathit{r})=\underset{i=1}{\overset{10}{\Pi }}[1-{(1-r_i)}^{x_i}]$$

(15)

$$s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{g}_1(\mathit{x},\mathit{r})=\sum _{i=1}^{10}{\upsilon }_i{x}_i^2-V\le 0$$

(16)

$$g_2(\mathit{x},\mathit{r})=\sum _{i=1}^{10}{\beta }_i{x}_i|\cos (2\pi r_i-\pi )|-C\le 0$$

(17)

$$g_3(x,r)=\sum _{i=1}^{10}{\omega }_i{x}_i{e}^{x_i/4}-W\le 0$$

(18)

$$0.9\le r_i\le 0.95,\begin{array}{c}\hfill \end{array}1\le x_i\le 5,\begin{array}{c}\hfill \end{array}{x}_i\in N_{+}\begin{array}{c}\hfill \end{array}(i=1,2,\cdots ,10)$$

(19)

where the objective function equation (15) is the reliability of this series system. The first constraint equation (16) is the volume constraint related to the redundancy level. The second constraint equation (17) is the budget constraint that is related to the redundancy level and reliability choice. The third constraint equation (18) is the weight constraint. $\mathbf{\upsilon }=\{{\upsilon }_i,i=1,2,\cdots ,n\}$ , $\mathbf{\beta }=\{{\beta }_i,i=1,2,\cdots ,n\}$ , and $\mathbf{\omega }=\{{\omega }_i,i=1,2,\cdots ,n\}$ are the parameters of constraints equations (16) to (18), respectively. Here, $n=10$ , $\mathbf{\upsilon }=\{1,2,3,3,2,3,4,3,2,4\}$ , $\mathbf{\beta }=\{5.0,4.5,3.5,3.5,3.5,5.0,2.5,3.5,2.5,2.5\}$ , and $\mathbf{\omega }=\{6,5,5,6,4,5,4,5,4,6\}$ . The thresholds of these three constraints are $V=120$ , $C=80$ , and $W=180$ . In this paper, equations (15), (16), and (18) are directly taken from the literature. But the budget constraint is modified, which considers not only the redundancy level $\mathit{x}$ , but also the reliability choice $\mathit{r}$ . Last but not least, we have to point out that the constraints in this paper are dimensionless, which can be viewed as constraints after standardization.

For comparison, we first list the results obtained by GA under different population sizes ( $N$ ) and different crossover fractions ( $cf$ ) in Tables 2 and 3, respectively. It is seen that the population size $N$ has small impact on the optimal solution of this example, whereas the crossover fraction $\mathrm{cf}$ has large influence on the optimal solution. Moreover, when $\mathrm{cf}\ge 0.5$ the effect of crossover fraction $\mathrm{cf}$ on optimal solution is decreased. Meanwhile, we can also observe that the population size and crossover fraction all have evident impact on the computational efficiency.

Table 2.

Optimal solutions of example 1 obtained by GA under different population sizes ( $\mathrm{cf}=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.94998	0.94998	0.94992	0.94950	0.94991
$r_2$	0.94999	0.94992	0.94998	0.94919	0.95000
$r_3$	0.94996	0.94994	0.94998	0.94949	0.94992
$r_4$	0.94995	0.94994	0.94991	0.94941	0.94996
$r_5$	0.94998	0.94987	0.94996	0.94916	0.94982
$r_6$	0.94996	0.95000	0.94992	0.94976	0.94992
$r_7$	0.95000	0.94998	0.94998	0.94966	0.94995
$r_8$	0.94997	0.94995	0.94964	0.94914	0.94995
$r_9$	0.94823	0.94996	0.94991	0.94938	0.94999
$r_{10}$	0.94992	0.94999	0.94995	0.94984	0.94998
$x_1~x_{10}$	2,2,2,2,2,2,2,2,3,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,2,2,2,2,3,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2
$R(\mathit{x},\mathit{r})$	0.97756	0.97757	0.97755	0.97714	0.97756
$g_1(\mathit{x},\mathit{r})$	−2	−2	−2	−2	−2
$g_2(\mathit{x},\mathit{r})$	−9.1766	−8.2028	−9.1584	−8.2790	−8.2055
$g_3(\mathit{x},\mathit{r})$	−2.9136	−2.9136	−2.9136	−2.9136	−2.9136
Time (s)	5.32	10.70	17.42	7.36	20.14

Table 3.

Optimal solutions of example 1 obtained by GA under different choices of $\mathrm{cf}$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.90053	0.94655	0.94998	0.94997	0.94997
$r_2$	0.92217	0.94440	0.94997	0.94994	0.94999
$r_3$	0.92731	0.93699	0.94998	0.94996	0.94996
$r_4$	0.93438	0.94019	0.94998	0.94996	0.94998
$r_5$	0.91820	0.93400	0.94973	0.94997	0.94960
$r_6$	0.92184	0.93347	0.94998	0.94996	0.94994
$r_7$	0.91134	0.94482	0.94993	0.94995	0.94996
$r_8$	0.92773	0.94321	0.94995	0.94990	0.94998
$r_9$	0.92967	0.93881	0.94992	0.94976	0.94999
$r_{10}$	0.92020	0.93539	0.94997	0.94997	0.94996
$x_1~x_{10}$	2,2,2,2,3,2,2,2,3,1	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,2,2,2,2,3,2	2,2,2,2,3,2,2,2,2,2
$R(\mathit{x},\mathit{r})$	0.87902	0.9681	0.97757	0.97756	0.97757
$g_1(\mathit{x},\mathit{r})$	−4	−2	−2	−2	−2
$g_2(\mathit{x},\mathit{r})$	−13.7815	−9.8930	−8.2050	−9.1559	−8.2071
$g_3(\mathit{x},\mathit{r})$	−2.7799	−2.9136	−2.9136	−2.9136	−2.9136
Time (s)	3.78	11.92	37.74	31.52	11.38

Then, to illustrate the performance of the proposed approach for exploring the optimal design of this system, the results obtained by the proposed approach under different population sizes $N$ and different choices of $p_0$ are presented in Tables 4 and 5, respectively. It is worth pointing out that the solutions obtained by the proposed approach are the average results coming by repeatedly and independently running the proposed algorithm for 20 times. We can conclude that the population size has evident impact on the computational efficiency, but has little impact on the final optimal solution. In addition, the choice of $p_0$ has relatively large impact on the optimal solution and computational efficiency. When $p_0<0.5$ the optimal solution seems to be suboptimal. But when $p_0\ge 0.5$ , $p_0$ has small influence on the computational accuracy. Moreover, we can conclude that the proposed approach tends to obtain better solutions compared with GA in terms of the system reliability, although the proposed approach needs longer CPU time. That is to say, the system obtained by the proposed approach tends to be more reliable.

Table 4.

Optimal solutions of example 1 obtained by the proposed approach under different population sizes ( $p_0=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.94970	0.95000	0.95000	0.95000	0.94981
$r_2$	0.95000	0.95000	0.95000	0.95000	0.95000
$r_3$	0.94999	0.95000	0.95000	0.95000	0.95000
$r_4$	0.94990	0.95000	0.95000	0.95000	0.95000
$r_5$	0.95000	0.95000	0.94982	0.95000	0.95000
$r_6$	0.95000	0.95000	0.94997	0.95000	0.95000
$r_7$	0.95000	0.95000	0.95000	0.95000	0.95000
$r_8$	0.94999	0.95000	0.95000	0.94997	0.95000
$r_9$	0.95000	0.95000	0.95000	0.94977	0.95000
$r_{10}$	0.95000	0.94989	0.95000	0.95000	0.94996
$x_1~x_{10}$	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,2,2,2,2,3,2	2,2,2,2,2,2,2,2,3,2
$R(\mathit{x},\mathit{r})$	0.97756	0.97759	0.97760	0.97760	0.97758
$g_1(\mathit{x},\mathit{r})$	−2	−2	−2	−2	−2
$g_2(\mathit{x},\mathit{r})$	−8.2026	−8.1963	−8.1995	−9.1501	−9.1504
$g_3(\mathit{x},\mathit{r})$	−2.9136	−2.9136	−2.9136	−2.9136	−2.9136
Time (s)	49.89	51.71	58.15	56.32	54.26

Table 5.

Optimal solutions of example 1 obtained by the proposed approach under different choices of $p_0$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.94966	0.95000	0.95000	0.94998	0.95000
$r_2$	0.95000	0.94975	0.95000	0.95000	0.95000
$r_3$	0.94991	0.95000	0.95000	0.95000	0.95000
$r_4$	0.94994	0.94975	0.95000	0.95000	0.95000
$r_5$	0.93691	0.95000	0.95000	0.95000	0.95000
$r_6$	0.95000	0.95000	0.94997	0.94991	0.95000
$r_7$	0.95000	0.95000	0.95000	0.94998	0.95000
$r_8$	0.95000	0.95000	0.95000	0.95000	0.95000
$r_9$	0.94893	0.95000	0.95000	0.95000	0.95000
$r_{10}$	0.94971	0.94921	0.94989	0.94995	0.95000
$x_1~x_{10}$	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,2,2,2,2,3,2	2,2,2,2,3,2,2,2,2,2	2,2,2,2,3,2,2,2,2,2
$R(\mathit{x},\mathit{r})$	0.97729	0.97747	0.97759	0.97758	0.97761
$g_1(\mathit{x},\mathit{r})$	−2	−2	−2	−2	−2
$g_2(\mathit{x},\mathit{r})$	−8.5176	−8.2109	−9.1480	−8.1983	−8.1953
$g_3(\mathit{x},\mathit{r})$	−2.9136	−2.9136	−2.9136	−2.9136	−2.9136
Time (s)	54.03	59.24	56.60	58.94	52.99

The possible mechanism of population size $N$ affecting the calculation time is as follows. When the number of per-generation individuals are small, although the calculation efficiency of each generation is high, such small-size population cannot fully explore the feasible region. Therefore, it needs more iterations to reach the optimal solution, which eventually leads to the reduction of calculation efficiency. When the number of individuals in each generation are large, although the computational efficiency of each generation decreases, such large-size population fully excavates the information of the feasible region. So it needs less iterations to find the optimal solution and the overall efficiency of the optimization algorithm is high.

The possible principle of the parameter $p_0$ affecting the computational efficiency is as follows. When $p_0$ is small, only small-fraction of better individuals are retained to generate new individuals. This theoretically accelerates the convergence speed of the proposed algorithm (the objective function value declines faster), but these smaller-part better individuals may lead to the insufficient use of feasible individuals, making the search direction deviate from the optimal one, ultimately affecting the efficiency of the whole algorithm. When $p_0$ is large, each generation retains more number of excellent individuals. Theoretically, the objective function value declines slowly and it needs longer time to procure the optimal solution. However, this can make full use of the information of feasible individuals. It thus needs less iterations to achieve the stopping condition.

An oil transportation system

Consider an oil transportation system shown in Figure 6, which is abstracted from practical applications. This system needs to transport the oil from the starting city to the destination city. This system possesses multiple transportation routes, and there are multiple transfer stations built along each route. The oil continues to be transported after passing through transfer stations till it is transported to the destination. Now the oil supplier needs to establish these transfer stations at the necessary 10 positions, that is, there are 10 subsystems in this oil transportation system. It is supposed that the reliability of this system is only reliant upon the reliabilities of transfer stations $\{R_1,R_2,\cdots ,R_{10}\}$ , that is, other facilities such as the pipes of this system are very reliable (with the reliability approximate to 1). Meanwhile, the behavior of a transfer station does not affect that of others.

Figure 6.

Schematic view of example 2.

Then, the reliability of this system can be obtained according to its structure, as follows:

$$\begin{array}{cc}\multicolumn{1}{c}{}& R_s=R_1{R}_2{R}_3+R_2{R}_3{R}_4+R_7{R}_9{R}_{10}+R_8{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_4+R_1{R}_2{R}_5{R}_6+R_2{R}_4{R}_5{R}_6-R_7{R}_8{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_5{R}_6\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}-R_2{R}_3{R}_4{R}_5{R}_6-R_1{R}_2{R}_4{R}_5{R}_6+R_1{R}_2{R}_3{R}_4{R}_5{R}_6-R_1{R}_2{R}_3{R}_7{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_8{R}_9{R}_{10}-R_2{R}_3{R}_4{R}_7{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}-R_2{R}_3{R}_4{R}_8{R}_9{R}_{10}-R_1{R}_2{R}_5{R}_6{R}_7{R}_9{R}_{10}-R_2{R}_4{R}_5{R}_6{R}_7{R}_9{R}_{10}+R_1{R}_2{R}_3{R}_4{R}_7{R}_9{R}_{10}-R_1{R}_2{R}_5{R}_6{R}_8{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}-R_2{R}_4{R}_5{R}_6{R}_8{R}_9{R}_{10}+R_1{R}_2{R}_3{R}_4{R}_8{R}_9{R}_{10}+R_1{R}_2{R}_3{R}_7{R}_8{R}_9{R}_{10}+R_2{R}_3{R}_4{R}_7{R}_8{R}_9{R}_{10}+R_1{R}_2{R}_3{R}_5{R}_6{R}_7{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+R_2{R}_3{R}_4{R}_5{R}_6{R}_7{R}_9{R}_{10}+R_1{R}_2{R}_4{R}_5{R}_6{R}_7{R}_9{r}_{10}+R_1{R}_2{R}_3{R}_5{R}_6{R}_8{R}_9{R}_{10}+R_2{R}_3{R}_4{R}_5{R}_6{R}_8{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+R_1{R}_2{R}_4{R}_5{R}_6{R}_8{R}_9{R}_{10}+R_1{R}_2{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}+R_2{R}_4{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_4{R}_7{R}_8{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}-R_1{R}_2{R}_3{R}_4{R}_5{R}_6{R}_7{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_4{R}_5{R}_6{R}_8{R}_9{R}_{10}-R_1{R}_2{R}_3{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}-R_2{R}_3{R}_4{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}-R_1{R}_2{R}_4{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}+R_1{R}_2{R}_3{R}_4{R}_5{R}_6{R}_7{R}_8{R}_9{R}_{10}\hfill \end{array}$$

(20)

where $R_i=1-(1-r_i{)}^{x_i}$ ( $i=1,2,\cdots ,10$ ) is the reliability of subsystem $i$ , in which $r_i$ and $x_i$ are the reliability choice and redundancy level of subsystem $i$ . The design task of this system is to maximize the system reliability $R_s$ through allocating the component reliability and corresponding redundancy level, under the volume constraint equation (21), weight constraint equation (22), and budget constraint equation (23), as follows:

$\underset{(\mathit{x},\mathit{r})}{max}:\begin{array}{c}\hfill \end{array}$ equation (20)

$$s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{g}_1(\mathit{x},\mathit{r})=\sum _{i=1}^{10}{\upsilon }_i{x}_i^2-V\le 0$$

(21)

$$g_2(x,r)=\sum _{i=1}^{10}{\omega }_i{e}^{x_i/4}\log (x_i+1)-W\le 0$$

(22)

$$g_3(\mathit{x},\mathit{r})=\sum _{i=1}^{10}{\beta }_i{x}_i|\cos (2\pi r_i-\pi )|-C\le 0$$

(23)

$$0.9\le r_i\le 0.95,\begin{array}{c}\hfill \end{array}1\le x_i\le 3,\begin{array}{c}\hfill \end{array}{x}_i\in N_{+},\begin{array}{c}\hfill \end{array}i=1,2,\cdots ,10$$

(24)

Here, $\mathbf{\upsilon }=\{0.11,0.20,0.13,0.16,0.15,0.14,0.10,0.14,0.16\}$ , $\mathbf{\beta }=\{0.92,0.92,0.91,0.93,0.92,0.92,0.92,0.90,0.94,0.93\}$ , and $\mathbf{\omega }=\{0.37,0.13,0.25,0.15,0.15,0.18,0.24,0.18,0.31,0.37\}$ . The thresholds of these constraints are $V=5$ , $W=5$ , and $C=10$ .

Similarly, we first apply GA to obtain comparison results. The optimal solutions obtained by GA under different population sizes and crossover fractions are listed in Tables 6 and 7, respectively. It can be concluded that the population size affects the final optimal solution as well as the computational efficiency. With the increased size of population, GA needs longer time to procure the optimal solution. Meanwhile, the crossover fraction also severely affects the computational precision and efficiency. When the crossover fraction is smaller than 0.5, the obtained optimal solution seems to be suboptimal.

Table 6.

Optimal solutions of example 2 obtained by GA under different population sizes ( $\mathrm{cf}=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.90033	0.95000	0.94999	0.94998	0.90112
$r_2$	0.94985	0.95000	0.94998	0.95000	0.94985
$r_3$	0.94197	0.95000	0.94990	0.95000	0.94077
$r_4$	0.90009	0.95000	0.95000	0.95000	0.90016
$r_5$	0.90036	0.90000	0.90039	0.90000	0.90049
$r_6$	0.90042	0.90000	0.90002	0.90000	0.90055
$r_7$	0.94994	0.90000	0.90436	0.90001	0.94890
$r_8$	0.94987	0.90000	0.90013	0.90001	0.94984
$r_9$	0.94963	0.95000	0.94998	0.94999	0.94995
$r_{10}$	0.94996	0.95000	0.94999	0.95000	0.94964
$x_1~x_{10}$	1,1,1,1,1,1,1,1,2,2	1,2,2,1,1,1,1,1,1,1	1,2,2,1,1,1,1,1,1,1	1,2,2,1,1,1,1,1,1,1	1,1,1,1,1,1,1,1,2,2
$R_s$	0.999473	0.999418	0.999419	0.999418	0.999469
$g_1(\mathit{x},\mathit{r})$	−2.6700	−2.5800	−2.5800	−2.5800	−2.6700
$g_2(\mathit{x},\mathit{r})$	−2.2998	−2.5762	−2.5762	−2.5762	−2.2998
$g_3(\mathit{x},\mathit{r})$	0.0003	−0.0202	−0.0043	−0.0202	−0.0003
Time (s)	2.35	4.77	7.45	8.44	12.17

Table 7.

Optimal solutions of example 2 obtained by GA under different choices of $\mathrm{cf}$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.95000	0.95000	0.90007	0.95000	0.90079
$r_2$	0.90000	0.94999	0.95000	0.94999	0.94994
$r_3$	0.95000	0.90072	0.94351	0.94963	0.93694
$r_4$	0.95000	0.94999	0.90002	0.94998	0.90049
$r_5$	0.95000	0.94998	0.90003	0.90010	0.90170
$r_6$	0.95000	0.94995	0.90000	0.90002	0.90098
$r_7$	0.90436	0.90042	0.94993	0.90499	0.94971
$r_8$	0.94998	0.90104	0.94981	0.90011	0.94957
$r_9$	0.95000	0.94999	0.94994	0.95000	0.94985
$r_{10}$	0.95000	0.94999	0.94996	0.94998	0.94992
$x_1~x_{10}$	1,3,1,1,1,1,1,1,1,1	1,2,2,1,1,1,1,1,1,1	1,1,1,1,1,1,1,1,2,2	1,2,2,1,1,1,1,1,1,1	1,1,1,1,1,1,1,1,2,2
$R_s$	0.999149	0.999367	0.999478	0.999419	0.999468
$g_1(\mathit{x},\mathit{r})$	−1.9700	−2.5800	2.6700	−2.5800	−2.6700
$g_2(\mathit{x},\mathit{r})$	−2.6604	−2.5762	−2.2998	−2.5762	−2.2998
$g_3(\mathit{x},\mathit{r})$	0.0010	−0.0078	0.0010	−0.0042	−0.0012
Time (s)	7.93	10.45	14.53	11.76	14.87

The optimal solutions obtained by the proposed approach under different population sizes $N$ and choices of $p_0$ are listed in Tables 8 and 9, respectively. It can be seen that the population size influences the computational precision and speed, but the influence on computational efficiency is relatively small when $N=1000~5000$ . We can also say that $p_0$ has a large effect on the computational precision, but has small impact on computation time. Moreover, when $p_0$ is proper (e.g., $p_0\ge 0.5$ ), the effect of $p_0$ on the computational accuracy is small. Last but not least, although the CPU time consumed by the proposed approach is longer than that consumed by GA, the proposed approach enhances the reliability of the elaborately designed system on the premise that the slight violation of constraint is acceptable.

Table 8.

Optimal solutions of example 2 obtained by the proposed approach under different population sizes ( $p_0=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.90000	0.90000	0.90000	0.90000	0.90000
$r_2$	0.95000	0.94985	0.95000	0.95000	0.94999
$r_3$	0.94010	0.94999	0.94999	0.94267	0.95000
$r_4$	0.90386	0.90001	0.90000	0.90000	0.90667
$r_5$	0.90006	0.91720	0.90026	0.90168	0.90002
$r_6$	0.90000	0.90000	0.90000	0.90016	0.90740
$r_7$	0.95000	0.95000	0.95000	0.95000	0.94995
$r_8$	0.94996	0.94999	0.94986	0.95000	0.95000
$r_9$	0.95000	0.95000	0.95000	0.94996	0.95000
$r_{10}$	0.95000	0.95000	0.94912	0.95000	0.95000
$x_1~x_{10}$	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2
$R_s$	0.999477	0.999492	0.999482	0.999479	0.999495
$g_1(\mathit{x},\mathit{r})$	−2.6700	−2.6700	−2.6700	−2.6700	−2.6700
$g_2(\mathit{x},\mathit{r})$	−2.2998	−2.2998	−2.2998	2.2998	−2.2998
$g_3(\mathit{x},\mathit{r})$	0.0076	0.0673	0.0110	0.0058	0.0601
Time (s)	46.77	43.76	43.36	42.07	41.98

Table 9.

Optimal solutions of example 2 obtained by the proposed approach under different choices of $p_0$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.94958	0.95000	0.90000	0.90000	0.90117
$r_2$	0.95000	0.95000	0.95000	0.95000	0.95000
$r_3$	0.95000	0.94913	0.95000	0.95000	0.94949
$r_4$	0.94930	0.94992	0.90000	0.90001	0.90310
$r_5$	0.94119	0.95000	0.90000	0.90272	0.90000
$r_6$	0.95000	0.95000	0.90000	0.90053	0.90000
$r_7$	0.90000	0.90035	0.94876	0.94750	0.95000
$r_8$	0.90397	0.93474	0.95000	0.95000	0.95000
$r_9$	0.94966	0.90000	0.95000	0.95000	0.95000
$r_{10}$	0.90001	0.95000	0.95000	0.95000	0.95000
$x_1~x_{10}$	1,2,1,1,1,1,1,1,1,2	1,2,1,1,1,1,1,1,2,1	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2	1,1,1,1,1,1,1,1,2,2
$R_s$	0.999290	0.999348	0.999484	0.999481	0.999490
$g_1(\mathit{x},\mathit{r})$	−2.4900	−2.5500	−2.6700	−2.6700	−2.6700
$g_2(\mathit{x},\mathit{r})$	−2.4656	−2.5209	−2.2998	−2.2998	−2.2998
$g_3(\mathit{x},\mathit{r})$	−0.0107	0.0994	0.0113	0.0200	0.0272
Time (s)	41.06	44.56	42.62	44.79	44.72

A 15-unit system

Consider a 15-unit system shown in Figure 7, which is modified from Refs.^17,32 The system reliability is:

$$\begin{array}{cc}\multicolumn{1}{c}{}& R_s=R_1{R}_2{R}_3{R}_4{R}_5{R}_6+R_9{R}_{10}{R}_{11}{R}_{12}{R}_{13}{R}_{14}{R}_{15}(Q_1+R_1{Q}_2+R_1{R}_2{Q}_3+R_1{R}_2{R}_3{Q}_4+R_1{R}_2{R}_3{R}_4{Q}_5+R_1{R}_2{R}_3{R}_4{R}_5{Q}_6)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+R_4{R}_5{R}_6{R}_7{R}_9{R}_{10}(Q_{11}+R_{11}{Q}_{12}+R_{11}{R}_{12}{Q}_{13}+R_{11}{R}_{12}{R}_{13}{Q}_{14}+R_{11}{R}_{12}{R}_{13}{R}_{14}{Q}_{15})(Q_1+R_1{Q}_2)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+\{(Q_1+R_1{Q}_2)(Q_3+R_3{Q}_4+R_3{R}_4{Q}_7)+R_1{R}_2{Q}_7(Q_3+R_3{Q}_4)\}(Q_{13}+R_{13}{Q}_{14}+R_{13}{R}_{14}{Q}_{15})R_5{R}_6{R}_8{R}_9{R}_{10}{R}_{11}{R}_{12}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+R_1{R}_2{R}_5{R}_6{R}_7{R}_8{R}_{11}{R}_{12}(R_9{R}_{10}+Q_9+R_9{Q}_{10})(Q_3+R_3{Q}_4)(Q_{13}+R_{13}{Q}_{14}+R_{13}{R}_{14}{Q}_{15})+(Q_5+R_5{Q}_6)\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\times \{(Q_7+R_7{Q}_{11}+R_7{R}_{11}{Q}_{12})(Q_9+R_9{Q}_{10})+R_9{R}_{10}(Q_{11}+R_{11}{Q}_{12})\}{R}_1{R}_2{R}_3{R}_4{R}_8{R}_{13}{R}_{14}{R}_{15}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}+R_1{R}_2{R}_7{R}_{11}{R}_{12}{R}_{13}{R}_{14}{R}_{15}(Q_9+R_9{Q}_{10})(Q_3+R_3{Q}_4+R_3{R}_4{Q}_5+R_3{R}_4{R}_5{Q}_6)+R_3{R}_4{R}_7{R}_8{R}_9{R}_{10}{R}_{13}{R}_{14}{R}_{15}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\times (Q_1+R_1{Q}_2)(Q_{11}+R_{11}{Q}_{12})(Q_5+R_5{Q}_6)\hfill \end{array}$$

(25)

where $R_i=1-(1-r_i{)}^{x_i}$ ( $i=1,2,\cdots ,15$ ) is the reliability of the $i$ th subsystem; $r_i$ and $x_i$ are the component reliability and redundancy level of subsystem $i$ , respectively; $Q_i=1-R_i$ is the unreliability (or failure probability) of subsystem $i$ .

Figure 7.

Structure of the 15-unit system.

The optimization model of this 15-unit system is:

$\underset{(\mathit{x},\mathit{r})}{max}:\begin{array}{c}\hfill \end{array}$ equation (25)

$$s.t.\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}\begin{array}{c}\hfill \end{array}{g}_1(\mathit{x},\mathit{r})=\sum _{i=1}^{15}{\upsilon }_i{x}_i-V\le 0$$

(26)

$$g_2(\mathit{x},\mathit{r})=\sum _{i=1}^{15}{\beta }_i{x}_i|\cos (2\pi r_i-\pi )|-C\le 0$$

(27)

$$0.8\le r_i\le 0.85,\begin{array}{c}\hfill \end{array}1\le x_i\le 5,\begin{array}{c}\hfill \end{array}{x}_i\in N_{+}\begin{array}{c}\hfill \end{array}(i=1,2,\cdots ,15)$$

(28)

where $\mathbf{\upsilon }=\{{\upsilon }_i,i=1,2,\cdots ,15\}$ and $\mathbf{\beta }=\{{\beta }_i,i=1,2,\cdots ,15\}$ are parameters of constraint functions, as follows:

$$\upsilon =\left\{14.66,32.21,26.65,41.47,17.52,48.21,14.93,25.81,32.18,29.33,40.89,28.67,14.14,34.96,12.12\right\}$$

(29)

$$\beta =\left\{1.18,2.65,2.31,2.19,2.86,2.46,1.86\text{,}2.70,2.49,1.88,2.00,1.83,2.34,1.78,2.56\right\}$$

(30)

Moreover, $V=1000$ and $C=40$ . The system reliability function comes from Refs.^17,32 The volume constraint equation (26) is established by mimicking that in Valian et al., ¹⁷ and the budget constraint equation (27) is a newly added one since the original paper did not consider the budget constraint.

We first apply GA to solve this optimal reliability design problem, and the optimal solutions under different population sizes and choices of crossover fractions are listed in Tables 10 and 11, respectively. In this example, the population size affects the final optimal solution and CPU time. In addition, the crossover fraction severely influences the final optimal solution, but has little effect on the CPU time.

Table 10.

Optimal solutions of example 3 obtained by GA under different population sizes ( $\mathrm{cf}=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.82441	0.82782	0.83536	0.84067	0.82972
$r_2$	0.80701	0.84756	0.81834	0.83773	0.82730
$r_3$	0.82426	0.84312	0.81039	0.83741	0.81358
$r_4$	0.82934	0.82941	0.83834	0.84028	0.83334
$r_5$	0.84687	0.83864	0.81518	0.84135	0.83632
$r_6$	0.84967	0.82174	0.82207	0.84397	0.84667
$r_7$	0.82437	0.81379	0.84644	0.83003	0.81496
$r_8$	0.82341	0.81777	0.82623	0.84418	0.80848
$r_9$	0.84646	0.83205	0.84555	0.82698	0.84435
$r_{10}$	0.83173	0.81644	0.81501	0.83256	0.82704
$r_{11}$	0.82833	0.82144	0.82723	0.84062	0.83453
$r_{12}$	0.82093	0.83711	0.83434	0.84532	0.82901
$r_{13}$	0.80828	0.84661	0.81707	0.82306	0.82577
$r_{14}$	0.82815	0.83586	0.84691	0.82771	0.84898
$r_{15}$	0.81908	0.82888	0.81974	0.82623	0.82443
$x_1~x_{15}$	4,3,3,3,4,3,3,1,2,2,2,2,2,1,2	3,2,3,2,3,3,3,2,3,2,3,2,2,1,3	3,2,2,2,3,2,2,2,2,2,3,3,3,3,3	2,2,3,3,2,2,2,1,3,3,2,2,3,3,3	2,2,2,2,2,1,2,2,3,3,3,3,4,3,4
$R_s$	0.997636	0.995825	0.996593	0.997872	0.997883
$g_1(\mathit{x},\mathit{r})$	−5.4400	−0.3000	−9.5400	−7.4600	−2.1600
$g_2(\mathit{x},\mathit{r})$	−1.4250	−0.2549	−2.7238	−0.0153	−0.1641
Time (s)	0.51	0.36	0.64	1.12	0.60

Table 11.

Optimal solutions of example 3 obtained by GA under different choices of $\mathrm{cf}$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.84544	0.82078	0.82960	0.83836	0.84167
$r_2$	0.82777	0.81922	0.83160	0.83486	0.82534
$r_3$	0.81224	0.80860	0.81480	0.83236	0.83489
$r_4$	0.82802	0.84069	0.84079	0.80649	0.83574
$r_5$	0.82149	0.82173	0.80883	0.82370	0.80508
$r_6$	0.81890	0.84947	0.80589	0.84436	0.84766
$r_7$	0.84149	0.83417	0.83799	0.81636	0.82922
$r_8$	0.80772	0.84403	0.82246	0.84254	0.84579
$r_9$	0.81899	0.81484	0.81779	0.81938	0.82101
$r_{10}$	0.83499	0.83405	0.82357	0.82400	0.84149
$r_{11}$	0.82422	0.84003	0.82356	0.83861	0.82837
$r_{12}$	0.81203	0.82801	0.83716	0.84267	0.83649
$r_{13}$	0.81366	0.81981	0.84850	0.83793	0.80269
$r_{14}$	0.83624	0.84115	0.80842	0.84838	0.84106
$r_{15}$	0.80367	0.84025	0.80798	0.83800	0.81521
$x_1~x_{15}$	4,2,3,2,3,2,2,2,2,3,2,3,3,2,4	4,2,3,2,3,2,2,2,3,3,2,2,2,3,2,3	3,3,2,2,4,3,2,2,2,2,2,3,2,2,3	3,3,3,3,4,2,2,1,2,2,2,2,2,3,3	4,3,3,2,4,3,2,2,2,2,2,2,2,2,2
$R_s$	0.996077	0.996801	0.996550	0.997717	0.997329
$g_1(\mathit{x},\mathit{r})$	−2.6300	−11.2400	−1.5900	−1.2000	−1.0700
$g_2(\mathit{x},\mathit{r})$	−3.6115	−0.6735	−4.1024	−0.0335	−1.6052
Time (s)	0.31	0.42	0.43	0.66	0.46

The optimal solutions obtained by the proposed approach under different population sizes and choices of $p_0$ are demonstrated in Tables 12 and 13, respectively. The following conclusions can be drawn: (1) the population size has relatively small effect on the computational efficiency and final optimal solution; (2) the choice of $p_0$ severely influences the optimal solution and CPU time; (3) although the proposed approach consumes longer CPU time than GA, it is more likely to produce higher reliable systems.

Table 12.

Optimal solutions of example 3 obtained by the proposed approach under different population sizes ( $p_0=0.8$ ).

	1000	2000	3000	4000	5000
$r_1$	0.85000	0.85000	0.85000	0.84706	0.85000
$r_2$	0.85000	0.85000	0.80342	0.85000	0.84996
$r_3$	0.84999	0.84842	0.84999	0.84819	0.84488
$r_4$	0.84994	0.85000	0.84932	0.85000	0.84526
$r_5$	0.80000	0.84544	0.85000	0.80000	0.80000
$r_6$	0.84965	0.85000	0.80003	0.84934	0.84997
$r_7$	0.81226	0.84350	0.80004	0.80427	0.83930
$r_8$	0.80199	0.80000	0.83157	0.84133	0.80000
$r_9$	0.85000	0.80121	0.85000	0.85000	0.84980
$r_{10}$	0.85000	0.84824	0.83997	0.80370	0.84919
$r_{11}$	0.84995	0.85000	0.85000	0.85000	0.80074
$r_{12}$	0.85000	0.80000	0.85000	0.85000	0.83744
$r_{13}$	0.80086	0.83673	0.84882	0.84996	0.84938
$r_{14}$	0.80001	0.84960	0.85000	0.81360	0.85000
$r_{15}$	0.80082	0.80000	0.80000	0.80000	0.80000
$x_1~x_{15}$	3,3,3,3,4,3,1,1,2,2,2,2,2,2,3	3,3,3,3,3,3,1,1,2,2,2,2,2,2,3	3,3,2,2,2,2,2,1,2,2,3,3,3,3,4	3,3,3,3,4,2,2,1,2,2,2,2,2,3,3	3,3,3,3,4,3,1,1,2,2,2,2,2,2,3
$R_s$	0.998696	0.998472	0.998434	0.998155	0.998692
$g_1(\mathit{x},\mathit{r})$	−2.8800	−20.4000	−8.5400	−1.2000	−2.8800
$g_2(\mathit{x},\mathit{r})$	−1.2587	−0.2037	−0.5864	−0.4332	−0.4757
Time (s)	85.40	83.70	85.39	88.70	88.53

Table 13.

Optimal solutions of example 3 obtained by the proposed approach under different choices of $p_0$ ( $N=5000$ ).

	0.1	0.3	0.5	0.7	0.9
$r_1$	0.85000	0.84999	0.85000	0.85000	0.84999
$r_2$	0.85000	0.85000	0.85000	0.84692	0.81507
$r_3$	0.83859	0.83118	0.84905	0.85000	0.84120
$r_4$	0.85000	0.85000	0.85000	0.85000	0.85000
$r_5$	0.85000	0.84826	0.80003	0.80000	0.80332
$r_6$	0.84866	0.84697	0.84957	0.85000	0.85000
$r_7$	0.80000	0.80000	0.85000	0.80087	0.80000
$r_8$	0.80101	0.80000	0.84990	0.80000	0.85000
$r_9$	0.85000	0.85000	0.80017	0.85000	0.84950
$r_{10}$	0.85000	0.83251	0.85000	0.85000	0.84420
$r_{11}$	0.82410	0.80000	0.84888	0.85000	0.85000
$r_{12}$	0.80000	0.85000	0.84310	0.84999	0.85000
$r_{13}$	0.85000	0.84260	0.81504	0.85000	0.80000
$r_{14}$	0.85000	0.84761	0.84659	0.80247	0.85000
$r_{15}$	0.80000	0.80000	0.80097	0.80000	0.83642
$x_1~x_{15}$	3,2,3,3,3,3,1,1,2,2,2,3,2,2,3	3,3,3,3,3,3,1,1,2,2,2,2,2,2,3	3,3,3,3,3,3,1,1,2,2,2,2,3,2,3	3,3,3,3,4,2,2,1,2,2,2,2,2,3,3	3,3,3,3,4,2,2,1,2,2,2,2,3,2,3
$R_s$	0.997777	0.998431	0.998609	0.998199	0.998173
$g_1(\mathit{x},\mathit{r})$	−23.9400	−20.4000	−6.2600	−1.2000	−22.0200
$g_2(\mathit{x},\mathit{r})$	−0.6789	−0.2301	−0.3498	−0.5139	−0.0048
Time (s)	87.63	77.76	78.62	82.10	76.29

Conclusions

Two issues are addressed in this paper, that is, proposing an effective optimization algorithm and applying it to solve RRAP. The proposed optimization approach is initialized as a small size population, then new populations are generated by VBI through sufficiently utilizing the information of the existing feasible individuals. Moreover, an automatic stratifying strategy is proposed to accelerate the speed of convergence. Thus, the new population will be generated iteratively based on informative feasible individuals, rather than all feasible individuals, until the optimal solution comes to light. The proposed optimization approach can avoid being caught into local optimum. It does not make any additional assumptions or preconditions for the optimization problem, so it can be easily applied to various kinds of optimization problems. After that, a sequential two-stage approach based on the proposed optimization algorithm is developed to solve RRAP, which divides RRAP into two parts, i.e., ROP and RAP, and solves these two problems sequentially and iteratively. This reduces the number of decision variables, thus decreasing the computation complexity of solving RRAP. The listed examples showcase that the proposed approach tends to obtain more reliable systems than GA, although it may need longer CPU time.