# A Globally Optimal Robust Design Method for Complex Systems

Authors: Yue Chen; Jian Shi; Xiao-jian Yi
Journal: Complexity
 (2020)

---

## Abstract

The uncertainty of the engineering system increases with the growing complexity of the engineering system; therefore, the tolerance to the uncertainty is essential. In the design phase, the output performance should reach the design criterion, even under large variations of design parameters. The tolerance to design parameter variations may be measured by the size of a solution space in which the output performance is guaranteed to deliver the required performance. In order to decouple dimensions, a maximum solution hyperbox, expressed by intervals with respect to each design parameter, is sought. The proposed approach combines the metaheuristic algorithm with the DIRECT algorithm where the former is used to seek the maximum size of hyperbox, and the latter is used as a checking technique that guarantees the obtained hyperbox is indeed a solution hyperbox. There are three advantages of the proposed approach. First, it is a global search and has a considerable high possibility to produce the globally maximum solution hyperbox. Second, it can be used for both analytically known and black-box performance functions. Third, it guarantees that any point selected within the obtained hyperbox satisfies the performance criterion as long as the performance function is continuous. Furthermore, the proposed approach is illustrated by numerical examples and real examples of complex systems. Results show that the proposed approach outperforms the GHZ and CES-IA methods in the literature.


---

## Body

## 1. Introduction

With the rapid development of technology in recent years, the systems applied in various engineering areas such as electronics, communication, and networking have become more and more complex. The increasing complexity of the system provides a new source for the uncertainty. Uncertainty arises because some design parameters cannot yet be specified exactly or they may be changed over the course of development [1, 2]. However, the higher requirements on systems place great importance on analyzing the uncertainty.Traditional optimization techniques seek an optimum in the design space. Typically, without considering the uncertainty, the optimum design is frequently pushed to the constraint boundary of the design. Consequently, this type of optimum design may be nonrobust and sensitive to parameter variabilities. Some authors even believe that optimization is actually just the opposite of robustness [3]. As reliability is nevertheless mandatory, engineers have to look for robust designs in order to avoid unexpected deviations from the nominal performance [4]. To this end, more advanced methods have been developed to include uncertainties of the parameters and robustness criteria in the optimization. These include(i)
Methods for selecting a design that is insensitive to the variations without eliminating possible variations of design parameters, i.e.,robust design optimization (RDO), see [5–7](ii)
Methods for identifying optimum designs which are characterized by a low probability of failure, i.e., reliability-based design optimization (RBDO), see [1, 2, 8](iii)
Methods for computing the effect of variability of input parameters on variability of objective function value, i.e.,sensitivity analysis (SA), see [9]Unfortunately, RDO, RBDO, and SA suffer from certain disadvantages which exclude them from some important applications. RBDO and RDO deal with models in the case where the variability of design parameters is given. If, however, the variability of the design parameters is not known completely, it has to be estimated, which is not always possible. When applying SA, information on how to improve a nonrobust solution is limited, what design parameter needs to be adjusted and what value it should assume is unknown.In this paper, we aim to infer the maximum uncertainty or variability of design parameters that we can tolerate without violating the required performance. That is, we focus on finding a maximum box-shaped solution space rather than a single optimum solution. The box-shaped solution space, representing a hyperbox, can be expressed by intervals for each design parameter and has the following benefits:(i)
The performance function delivers the required performance as long as the design parameters lie in the box-shaped solution space(ii)
For a design to be good, the choice of the value of one parameter within its assigned interval does not depend on the values of other parameters as long as they are within their respective intervals(iii)
The intervals may be used to assess the robustness and sensitivity of the uncertain design parameters by the widths of the associated intervals.(iv)
The intervals are independent of each other and may be combined with intervals of other disciplines, enabling distributed design in separate disciplinesTwo categories of approaches that solve similar problems could be identified in the literature. The first category is based on data-mining and machine-learning techniques [10–13]. The stochastic approach proposed in [10] combines query and online learning, which probes a box-shaped candidate solution space by stochastic sampling and then readjusts its boundaries in order to remove designs with insufficient performance and explore more design space that has not been probed before. The quality of the results and the efficiency of this stochastic approach are studied in detail in [11]. This stochastic approach, however, produces a box-shaped solution space that may contain some bad designs and may not have the best size because the Monte Carlo sampling is used to estimate the locations of good and bad designs. The second category is based on an analytical technique, called interval arithmetic [14]. The algorithm presented in [15], for example, applies interval arithmetic within an iterative optimization scheme to compute the feasibility of candidate boxes. Interval arithmetic, however, limits the applicability of the algorithm, as the accuracy of the results depends on the problem and cannot be assessed for general cases. Another drawback of the algorithm in [15] is that the performance function needs to be analytically known and cannot be treated as a black box [16].The maximum box-shaped solution space must both fulfill the “quantitative” viewpoint (i.e., its volume is maximum) and “feasible” viewpoint (i.e., it is indeed a solution space). To address these problems, we propose an innovative approach which combines the metaheuristic algorithm with the DIRECT algorithm. More specifically, we use the metaheuristic algorithm to search the maximum hyperbox and use the DIRECT algorithm to guarantee that all designs in the obtained hyperbox satisfy the performance criterion (i.e., the obtained hyperbox is indeed a solution hyperbox). Actually, in the proposed approach, metaheuristic algorithms do not rely on mathematical properties to be applied, and DIRECT algorithm only requires evaluations of the performance function. Therefore, our approach in essence treats the performance function as a black box and can be applied to various types of practical engineering problems for complex systems, such as reliability allocation, availability allocation, and maintenance analysis. As a real example of the complex system, the reliability allocation problem of the power-shift steering transmission control system (PSSTCS) of a heavy vehicle is analyzed by the proposed approach in this paper.The proposed approach is similar to the traditional RDO, as the variations of design parameters are considered. However, it differs from the traditional method in which the optimization is to seek the intervals of a permissible design range rather than a single design, and interval boundaries are used as degrees of freedom rather than design parameters.The paper is organized as follows. In Section2, motivated by a problem from engineering practice, the optimization problem is formulated for identifying the maximum hyperbox which guarantees the performance requirement. Section 3 introduces the DIRECT algorithm. Section 4 presents in detail the proposed approach. Section 5 discusses the numerical performance of the proposed approach. In Section 6, engineering cases are investigated. Section 7 presents some concluding remarks.


## 2. Motivation and Problem Statement

### 2.1. Example Problem

We consider the voltage divider shown in Figure1 in [15]. The terminal voltage Vab is expressed as Vab=V0R/R+Ra and the resistance Rin=R+Ra.Figure 1
Voltage divider example.Good designs fulfilling the design goals should satisfy(1)4.5V≤Vab≤5.5V,80Ω≤Rin≤120Ω.Figure2 shows the complete solution space (gray region). However, the classical optimum that minimizes Rin and maximizes Vab is not robust. The immediate neighborhood of the optimum includes designs that violate condition (1) and thus fails to meet at least one design goal.Figure 2
The complete solution space (gray region) and the solution box with the largest volume (region delimited by red lines).The classical optimum would not be a suitable choice because of the lack of robustness. Using the complete solution space instead would also be impractical, as the tolerable range forRa would depend on the choice of R. For decoupling the design parameters, hyperboxes contained in the complete solution space, expressed by intervals for each design parameter, are desirable.The robustness of uncertain design parameters [17] is defined as the widths of the associated intervals whereby the performance criterion (1) is met. Maximum robustness is provided if the intervals are as large as possible. The choice of intervals is not unique, and the size of the interval for Ra depends on the choice of the interval for R. The approach proposed here is to maximize the volume of the solution hyperbox, and the resulting unique intervals are shown in Figure 2.


### 2.2. Problem Statement

Letx=x1,x2,…,xm represent the design parameters or design points or designs, where m is the number of dimensions. The output performance at x is given by(2)y=fx,with f being the performance function. If the analytical form of f is not known, it is considered as a black box, and y may have to be computed numerically. In typical optimization problems, the performance is optimized, that is, a design x is sought such that y obtains an extreme value. In the approach presented here, the performance has to be sufficient by satisfying the performance criterion:(3)fx≥fc,with the threshold value fc. The performance criterion is the requirement of performance. For example, if the performance function is the reliability, the performance criterion is that the reliability should be above the given threshold (e.g., 0.99).Design points satisfying (3) are called good, otherwise are called bad. The design spaceDxl,xu is assumed to be continuous and defined by the lower and upper bounds of the design parameters, given by xl and xu, that is, Dxl,xu=xl,xu=x:xil≤xi≤xiu,i=1,…,m. The complete solution space is the set of all good design points in the design space. The shape of the complete solution space depends on the output performance and the performance criterion given by (2) and (3), respectively. In general, boundaries of the range of one parameter depend on the values of other parameters. In order to obtain interval boundaries of each parameter that are independent of other parameters, only subspaces that are hyperboxes are considered. The hyperbox D=Dxlow,xup can be expressed by(4)Dxlow,xup=defxlow,xup=x:xilow≤xi≤xiup,i=1,…,m,that is, the lower bound and upper bound of the hyperbox Dxlow,xup are xlow and xup, respectively. Moreover, let xbound denote the bounds of the hyperbox, namely,(5)xbound=defxlow,xup=x1low,…,xmlow,x1up,…,xmup.The volume of the hyperboxDxlow,xup is thus given by(6)vxlow,xup=def∏i=1mxiup−xilow.A hyperbox which only includes good design points is calleda solution hyperbox. Simply, the hyperbox which meets the condition fx≥fc for all x∈Dxlow,xup is the solution hyperbox. Furthermore, fminxlow,xup is defined as the global minimization of the performance function on the hyperbox, that is, fminxlow,xup=defminx∈Dxlow,xupfx.With these preparations at hand, the problem of searching for a maximum solution hyperbox can be formulated as the following constrained optimization problem:(7)maxxlow,xupvxlow,xupsubject tominx∈Dxlow,xupfx≥fcxl≤xlow≤xup≤xu.Iffx is analytically known and minx∈xlow,xupfx (i.e., fminxlow,xup) can be explicitly expressed by xlow and xup, the optimization problem (7) can be solved by traditional optimization algorithms. Unfortunately, in most of the engineering problems, the following two cases are more common. First, the performance function fx is not analytically known; therefore, it is produced by numerical simulations. Second, minx∈xlow,xupfx cannot be explicitly expressed by xlow and xup. To cope with these two cases, this paper aims to present a new approach based on performance function samples (i.e., the performance function can be analytically known or black box). There are two key difficulties for solving this optimization problem. First is how to maximize the size of the hyperbox. Second is how to guarantee the obtained hyperbox is a solution hyperbox.Remark 1.
If the performance criterion isfl≤fx≤fu, the optimization problem formulation (7) sets fx=defminfx−fl,−fx+fu and fc=def0.Remark 2.
Typically, engineering problems have multiple performance functionsfix≥fic,i=1,…,q. In the optimization problem formulation (7), minx∈Dfx≥fc is replaced by minx∈Dfix≥fic,i=1,…,q.


## 2.1. Example Problem

We consider the voltage divider shown in Figure1 in [15]. The terminal voltage Vab is expressed as Vab=V0R/R+Ra and the resistance Rin=R+Ra.Figure 1
Voltage divider example.Good designs fulfilling the design goals should satisfy(1)4.5V≤Vab≤5.5V,80Ω≤Rin≤120Ω.Figure2 shows the complete solution space (gray region). However, the classical optimum that minimizes Rin and maximizes Vab is not robust. The immediate neighborhood of the optimum includes designs that violate condition (1) and thus fails to meet at least one design goal.Figure 2
The complete solution space (gray region) and the solution box with the largest volume (region delimited by red lines).The classical optimum would not be a suitable choice because of the lack of robustness. Using the complete solution space instead would also be impractical, as the tolerable range forRa would depend on the choice of R. For decoupling the design parameters, hyperboxes contained in the complete solution space, expressed by intervals for each design parameter, are desirable.The robustness of uncertain design parameters [17] is defined as the widths of the associated intervals whereby the performance criterion (1) is met. Maximum robustness is provided if the intervals are as large as possible. The choice of intervals is not unique, and the size of the interval for Ra depends on the choice of the interval for R. The approach proposed here is to maximize the volume of the solution hyperbox, and the resulting unique intervals are shown in Figure 2.


## 2.2. Problem Statement

Letx=x1,x2,…,xm represent the design parameters or design points or designs, where m is the number of dimensions. The output performance at x is given by(2)y=fx,with f being the performance function. If the analytical form of f is not known, it is considered as a black box, and y may have to be computed numerically. In typical optimization problems, the performance is optimized, that is, a design x is sought such that y obtains an extreme value. In the approach presented here, the performance has to be sufficient by satisfying the performance criterion:(3)fx≥fc,with the threshold value fc. The performance criterion is the requirement of performance. For example, if the performance function is the reliability, the performance criterion is that the reliability should be above the given threshold (e.g., 0.99).Design points satisfying (3) are called good, otherwise are called bad. The design spaceDxl,xu is assumed to be continuous and defined by the lower and upper bounds of the design parameters, given by xl and xu, that is, Dxl,xu=xl,xu=x:xil≤xi≤xiu,i=1,…,m. The complete solution space is the set of all good design points in the design space. The shape of the complete solution space depends on the output performance and the performance criterion given by (2) and (3), respectively. In general, boundaries of the range of one parameter depend on the values of other parameters. In order to obtain interval boundaries of each parameter that are independent of other parameters, only subspaces that are hyperboxes are considered. The hyperbox D=Dxlow,xup can be expressed by(4)Dxlow,xup=defxlow,xup=x:xilow≤xi≤xiup,i=1,…,m,that is, the lower bound and upper bound of the hyperbox Dxlow,xup are xlow and xup, respectively. Moreover, let xbound denote the bounds of the hyperbox, namely,(5)xbound=defxlow,xup=x1low,…,xmlow,x1up,…,xmup.The volume of the hyperboxDxlow,xup is thus given by(6)vxlow,xup=def∏i=1mxiup−xilow.A hyperbox which only includes good design points is calleda solution hyperbox. Simply, the hyperbox which meets the condition fx≥fc for all x∈Dxlow,xup is the solution hyperbox. Furthermore, fminxlow,xup is defined as the global minimization of the performance function on the hyperbox, that is, fminxlow,xup=defminx∈Dxlow,xupfx.With these preparations at hand, the problem of searching for a maximum solution hyperbox can be formulated as the following constrained optimization problem:(7)maxxlow,xupvxlow,xupsubject tominx∈Dxlow,xupfx≥fcxl≤xlow≤xup≤xu.Iffx is analytically known and minx∈xlow,xupfx (i.e., fminxlow,xup) can be explicitly expressed by xlow and xup, the optimization problem (7) can be solved by traditional optimization algorithms. Unfortunately, in most of the engineering problems, the following two cases are more common. First, the performance function fx is not analytically known; therefore, it is produced by numerical simulations. Second, minx∈xlow,xupfx cannot be explicitly expressed by xlow and xup. To cope with these two cases, this paper aims to present a new approach based on performance function samples (i.e., the performance function can be analytically known or black box). There are two key difficulties for solving this optimization problem. First is how to maximize the size of the hyperbox. Second is how to guarantee the obtained hyperbox is a solution hyperbox.Remark 1.
If the performance criterion isfl≤fx≤fu, the optimization problem formulation (7) sets fx=defminfx−fl,−fx+fu and fc=def0.Remark 2.
Typically, engineering problems have multiple performance functionsfix≥fic,i=1,…,q. In the optimization problem formulation (7), minx∈Dfx≥fc is replaced by minx∈Dfix≥fic,i=1,…,q.


## 3. The DIRECT Algorithm

We consider the global optimization problem of the form(8)minx∈Dxlow,xupfx.As mentioned in Section 2, Dxlow,xup=xlow,xup=x:xilow≤xi≤xiup,i=1,…,m is an m-dimensional hyperbox (hyper-rectangle), and fx is a performance function.In the literature, there are numerous methods for dealing with problem (8) (see, e.g., [18–27]). DIRECT (DIviding RECTangles) is one of the most known partitioning-based algorithms that balance local and global search in an attempt to find the global minimizer efficiently [24–27]. The main procedure of the DIRECT algorithm involves partitioning potentially optimal (the most promising) hyper-rectangles and evaluating the performance function at the centers of these hyper-rectangles [28–31].DIRECT’s performance is independent of the scaling ofD and the problem is typically scaled as follows:(9)Tx=Ax−xlow,A=diag1x1up−x1low,…,1xmup−xmlow.Thus, the bound-constrained problem (8) can be reduced to the following problem:(10)minx∈D˜fA−1x+xlow,over the unit hypercube D˜=0,1m. That is, the search space is transformed to the m-dimensional unit hypercube.


### 3.1. Selection Scheme

To select the potentially optimal hyper-rectangles, DIRECT assesses the goodness based on the lower bound estimates for the performance function over each hyper-rectangle. The requirement of potential optimality is stated formally in Definition1.Definition 1.
(potentially optimal hyper-rectangle). Suppose that we have a partition of the unit hypercube inton hyper-rectangles. Let ci denote the center point of the ith hyper-rectangle, di denote the distance from the center point to the vertices, ε>0 be a positive constant, and fmin be the current best function value. A hyper-rectangle j is said to be potentially optimal if there exist some K˜>0 such that(11)fA−1cj+xlow−K˜dj≤fA−1ci+xlow−K˜di,for alli=1,…,n,(12)fA−1cj+xlow−K˜dj≤fmin−εfmin.
The hyper-rectanglej is potentially optimal if the lower Lipschitz bound for the performance function computed by the left-hand side of (11) is the smallest one with some positive constant K˜ among all hyper-rectangles. In (12), the parameter ε is used to protect from an excessive refinement of the local minimum [24, 31].
DIRECT identifies potentially optimal hyper-rectangles in at least two different ways: using modified Graham’s scan algorithm [32] or the rule described by Lemma 2.3 in [33]. Usually this does not impose significant difference; thus, in this paper, we use the modified Graham’s scan algorithm.


### 3.2. Division and Sampling Scheme

We start by sampling the pointsc±δei, i=1,…,m, where c is the center point of the hypercube, δ is one-third the side length of the hypercube, and ei is the ith unit vector (i.e., a vector with a one in the ith position and zeros elsewhere). The strategy in [31] is used to divide the hypercube.More precisely, we adopt the following rule for subdividing a hypercube. Let(13)ωi=minfA−1c+δei+xlow,fA−1c−δei+xlowbe the best of the function values sampled along dimension i. Then, we start by splitting along the dimension with the smallest ω value. Once this is done, we split the rectangle containing c into thirds along the dimension with the next smallest ω value. Continue in this way until we have split on all long dimensions. By dividing only along the long dimensions, we ensure that the rectangles shrink on every dimension. A formal description of the rectangle division procedure is given in Algorithm 1.Algorithm 1:Procedure for dividing rectangles.
(1)
Identify the setI of dimensions with the maximum side length. Let δ equal one-third of this maximum side length.
(2)
Sample the function at the pointsc±δei for all i∈I, where c is the center of rectangle and ei is the ith unit vector.
(3)
Divide the rectangle containingc into thirds along the dimensions in I, starting with the dimension with the lowest value of ωi=minfA−1c+δei+xlow,fA−1c−δei+xlow, and continuing to the dimension with the next smallest ωi until we have split on all dimensions i∈I.


### 3.3. The DIRECT Algorithm Flowchart

We now have all the ingredients for the DIRECT algorithm. We initialize the search by sampling at the center of the unit hypercube. Each iteration begins with identifying the set of potentially optimal hyper-rectangles as described in Section3.1. These hyper-rectangles are then sampled and subdivided as described in Section 3.2. The process continues until a prespecified iteration limit is reached. A formal statement of the DIRECT algorithm is given in Algorithm 2.Algorithm 2:The DIRECT algorithm.
(0) Parameter setting: set the maximum number of iterationsKmax.
(1) Normalize the search space to be the unit hypercube. Letc1 be the center point of this hypercube and evaluate fA−1c1+xlow. Set fmin=fA−1c1+xlow, n=1, and k=1.
(2) Identify the setS of potentially optimal rectangles.
(3) Select any rectanglej∈S.
(4) Use Algorithm1 to determine where to sample within rectangle j and how to divide the rectangle into subrectangles. Update fmin, set n=n+Δn, where Δn is the number of new points sampled.
(5) SetS=S−j; If S≠∅ go to Step 3.
(6) Setk=k+1. If k=Kmax, stop and output the minimum value fmin; Otherwise, go to Step 2.In order to better illustrate how the DIRECT algorithm works, Figure3shows the first three iterations of the algorithm on the two-dimensional Branin function in [34]. For each iteration, the first row shows the set of potentially optimal rectangles (gray colored), and the second row shows how these rectangles are sampled (red points) and subdivided.Figure 3
Three iterations of the DIRECT algorithm on Branin’s function: gray colored hyper-rectangles highlight potentially optimal ones. They are selected in the first phase of the current iteration and divided in the second phase of the current iteration. Red points are sample points where the performance function is evaluated. (a) Iteration 1, (b) Iteration 2, and (c) Iteration 3.
(a)
(b)
(c)


### 3.4. Convergence

The DIRECT algorithm is guaranteed to converge to the globally optimal function value if the performance function is continuous or at least continuous in the neighborhood of a global optimum [31]. This follows the fact that, as the number of iterations goes to infinity, the set of points sampled by DIRECT forms a dense subset of the unit hypercube. That is, given any point x in the unit hypercube and any η>0, DIRECT will eventually sample a point within the neighborhood of x with radius η. It is proved that DIRECT with a large enough iteration number can converge with a high accuracy [31, 35–37].


## 3.1. Selection Scheme

To select the potentially optimal hyper-rectangles, DIRECT assesses the goodness based on the lower bound estimates for the performance function over each hyper-rectangle. The requirement of potential optimality is stated formally in Definition1.Definition 1.
(potentially optimal hyper-rectangle). Suppose that we have a partition of the unit hypercube inton hyper-rectangles. Let ci denote the center point of the ith hyper-rectangle, di denote the distance from the center point to the vertices, ε>0 be a positive constant, and fmin be the current best function value. A hyper-rectangle j is said to be potentially optimal if there exist some K˜>0 such that(11)fA−1cj+xlow−K˜dj≤fA−1ci+xlow−K˜di,for alli=1,…,n,(12)fA−1cj+xlow−K˜dj≤fmin−εfmin.
The hyper-rectanglej is potentially optimal if the lower Lipschitz bound for the performance function computed by the left-hand side of (11) is the smallest one with some positive constant K˜ among all hyper-rectangles. In (12), the parameter ε is used to protect from an excessive refinement of the local minimum [24, 31].
DIRECT identifies potentially optimal hyper-rectangles in at least two different ways: using modified Graham’s scan algorithm [32] or the rule described by Lemma 2.3 in [33]. Usually this does not impose significant difference; thus, in this paper, we use the modified Graham’s scan algorithm.


## 3.2. Division and Sampling Scheme

We start by sampling the pointsc±δei, i=1,…,m, where c is the center point of the hypercube, δ is one-third the side length of the hypercube, and ei is the ith unit vector (i.e., a vector with a one in the ith position and zeros elsewhere). The strategy in [31] is used to divide the hypercube.More precisely, we adopt the following rule for subdividing a hypercube. Let(13)ωi=minfA−1c+δei+xlow,fA−1c−δei+xlowbe the best of the function values sampled along dimension i. Then, we start by splitting along the dimension with the smallest ω value. Once this is done, we split the rectangle containing c into thirds along the dimension with the next smallest ω value. Continue in this way until we have split on all long dimensions. By dividing only along the long dimensions, we ensure that the rectangles shrink on every dimension. A formal description of the rectangle division procedure is given in Algorithm 1.Algorithm 1:Procedure for dividing rectangles.
(1)
Identify the setI of dimensions with the maximum side length. Let δ equal one-third of this maximum side length.
(2)
Sample the function at the pointsc±δei for all i∈I, where c is the center of rectangle and ei is the ith unit vector.
(3)
Divide the rectangle containingc into thirds along the dimensions in I, starting with the dimension with the lowest value of ωi=minfA−1c+δei+xlow,fA−1c−δei+xlow, and continuing to the dimension with the next smallest ωi until we have split on all dimensions i∈I.


## 3.3. The DIRECT Algorithm Flowchart

We now have all the ingredients for the DIRECT algorithm. We initialize the search by sampling at the center of the unit hypercube. Each iteration begins with identifying the set of potentially optimal hyper-rectangles as described in Section3.1. These hyper-rectangles are then sampled and subdivided as described in Section 3.2. The process continues until a prespecified iteration limit is reached. A formal statement of the DIRECT algorithm is given in Algorithm 2.Algorithm 2:The DIRECT algorithm.
(0) Parameter setting: set the maximum number of iterationsKmax.
(1) Normalize the search space to be the unit hypercube. Letc1 be the center point of this hypercube and evaluate fA−1c1+xlow. Set fmin=fA−1c1+xlow, n=1, and k=1.
(2) Identify the setS of potentially optimal rectangles.
(3) Select any rectanglej∈S.
(4) Use Algorithm1 to determine where to sample within rectangle j and how to divide the rectangle into subrectangles. Update fmin, set n=n+Δn, where Δn is the number of new points sampled.
(5) SetS=S−j; If S≠∅ go to Step 3.
(6) Setk=k+1. If k=Kmax, stop and output the minimum value fmin; Otherwise, go to Step 2.In order to better illustrate how the DIRECT algorithm works, Figure3shows the first three iterations of the algorithm on the two-dimensional Branin function in [34]. For each iteration, the first row shows the set of potentially optimal rectangles (gray colored), and the second row shows how these rectangles are sampled (red points) and subdivided.Figure 3
Three iterations of the DIRECT algorithm on Branin’s function: gray colored hyper-rectangles highlight potentially optimal ones. They are selected in the first phase of the current iteration and divided in the second phase of the current iteration. Red points are sample points where the performance function is evaluated. (a) Iteration 1, (b) Iteration 2, and (c) Iteration 3.
(a)
(b)
(c)


## 3.4. Convergence

The DIRECT algorithm is guaranteed to converge to the globally optimal function value if the performance function is continuous or at least continuous in the neighborhood of a global optimum [31]. This follows the fact that, as the number of iterations goes to infinity, the set of points sampled by DIRECT forms a dense subset of the unit hypercube. That is, given any point x in the unit hypercube and any η>0, DIRECT will eventually sample a point within the neighborhood of x with radius η. It is proved that DIRECT with a large enough iteration number can converge with a high accuracy [31, 35–37].


## 4. The Proposed Approach

As metaheuristic algorithms are general and do not rely on mathematical properties for application and the DIRECT algorithm converges to the global optimum with a high accuracy, an innovative approach which combines metaheuristic algorithms with the DIRECT algorithm is proposed in this paper to tackle the two challenges mentioned in Section2.2. First, the metaheuristic algorithm is used for seeking the largest hyperbox. Second, the DIRECT algorithm is used as a checking technique to ensure all designs in the obtained hyperbox satisfy the performance criterion.More specifically, in this paper, two methods are proposed: one combines simulated annealing with the DIRECT algorithm and the other combines the distributed covariance matrix adaptation evolution strategy with the DIRECT algorithm. They will be referred to as the simulated annealing DIRECT algorithm (SA-DIRECT) and the distributed covariance matrix adaptation evolution strategy DIRECT algorithm (d CMA-ES-DIRECT), respectively. Of course, one can use other metaheuristic algorithms, but these two metaheuristic algorithms are derivative-free and insensitive to initial solution and global search.


### 4.1. Simulated Annealing DIRECT Algorithm (SA-DIRECT)

Simulated annealing (SA) proposed by Kirkpatrick et al. [38] is based on the similarity between the solid annealing process and solving global optimization problems. It is a generic probabilistic meta-algorithm for the global optimization problem. SA has been applied to a wide range of problems especially in those cases where traditional optimization techniques have shown poor performances or simply have failed [39].Simulated annealing DIRECT (SA-DIRECT) algorithm is illustrated in Algorithm3. Initially, the SA-DIRECT sets the design space and builds an initial solution for problem (7). More specifically, we generate a solution x0low,x0up, where x0low and x0up are stochastically generated from the range xl,xu. Then, we evaluate minx∈Dx0low,x0upfx by the DIRECT algorithm in order to check whether this generated solution meets constraints. If the constraint conditions, i.e., xl≤x0low≤x0up≤xu andminx∈Dx0low,x0upfx≥fc, are met, we calculate the associated value of the fitness function (i.e., the volume expressed by equation (6)). Otherwise, we discard this generated solution and generate a new initial solution.Algorithm 3: The proposed simulated annealing DIRECT (SA-DIRECT) algorithm.
(0) Parameter setting:
Set the initial temperatureT0, the cooling ratio β, the maximum number of the inner iterations lmax and the maximum number of the outer iterations tmax.
(1) Initialization:
(1.1) Set the lower and upper bounds of design space,xl and xu, and the threshold level fc.
(1.2) Generate randomlyx0low and x0up from the range xl,xu.
(1.3) Evaluatefmin=minx∈Dx0low,x0upfx by Algorithm 2 with Kmax=⌊50ln2m⌋.
(1.4) If all the constraints are satisfied. i.e.,xl≤x0low≤x0up≤xu and fmin≥fc, then calculate the volume by equation (6) and store it in v0, i.e., v0=vx0low,x0up; Otherwise, go to Step 1.2.
(1.5) Setxlow=xclow=x0low, xup=xcup=x0up, vb=vc=v0, and T=T0.
(2) Fort=1,…,tmaxdo:
(2.1) Setl=1.
(2.2) Whilel≤lmaxdo:
(2.2.1) Generate a new solutionxllow,xlup in the neighborhood of xclow,xcup, i.e., xllow=xclow+xu−xl⊙2rl−13, xlup=xcup+xu−xl⊙2ru−13, where rl and ru are randomly and uniformly generated from range 0,1m, and ⊙ denotes componentwise multiplication.
(2.2.2) Evaluatefmin=minx∈Dxllow,xlupfx by Algorithm 2 with the parameter setting Kmax=⌊100ln2m⌋.
(2.2.3) If constraintsxl≤xllow≤xlup≤xu and fmin≥fc are all met, then evaluate the volume vl by equation (6), i.e., vl=vxllow,xlup; Otherwise go to Step 2.2.1.
(2.2.4) EvaluateΔv=vl−vc, if Δv≥0, then set xclow=xllow and xcup=xlup; If vl−vb>0, then set xlow=xllow, xup=xlup.
(2.2.5) IfΔv<0, then generate a random number r in 0,1 and evaluate P=expΔv/T; If P>r, then set xclow=xllow, xcup=xlup, vc=vl.
(2.2.6) Setl=l+1.
(2.3) SetT=βT.
(3) Verification:
(3.1) Evaluatefmin=minx∈Dxlow,xupfx by Algorithm 2 with the parameter setting Kmax=⌊1000ln2m⌋.
(3.2) Iffmin≥fc, output the best solution xlow,xup; Otherwise go to Step 1.Next, during the process of optimization, in each iteration, a neighbor solution to the current solutionxllow,xlup is generated according to a predefined neighborhood structure, and minx∈Dxllow,xlupfx is evaluated by the DIRECT algorithm. If the condition minx∈Dxllow,xlupfx≥fc is not satisfied, this neighbor solution is rejected. Otherwise, the volume vxllow,xlup is calculated by equation (6). The improving move (the volume of the neighbor is larger than that of the current solution) is always accepted, whilst worse neighbor is accepted with a certain probability determined by the Boltzmann probability [40], expΔv/T>r, where Δv is the difference between the volume of the current solution and the generated neighbor, and r is randomly generated from the range 0,1. Moreover, T is a parameter (called the temperature) which periodically lowers during the search process according to the temperature updating rule. The temperature updating rule (as adopted in [41]) is T=βtT0, where T0 is the initial temperature, β is the cooling ratio, and t is the number of times the temperature has been lowered. The halting criterion of SA is “reaching the maximum number of times the temperature could be lowered.”Finally, the SA-DIRECT verifies whether the best solution meets the constraints by the DIRECT algorithm with a relatively large maximum number of iterations. If the constraints are not met, the SA-DIRECT goes back to Step 1. Otherwise, the SA-DIRECT outputs the solution.In this paper, the maximum number of the iterations at each temperaturelmax=200 and the maximum number of times the temperature could be lowered tmax=1000.


### 4.2. Distributed Covariance Matrix Adaptation Evolution Strategy DIRECT Algorithm (dCMA-ES-DIRECT)

The evolution strategy (ES) is developed as a powerful tool for numerical optimization tasks [42]. Covariant matrix adaptation evolution strategy (CMA-ES) acts as an improved robust form of evolution strategy [43]. The main feature of the CMA-ES is the ability of being invariant to landscape transformations and scaling modulation. The CMA-ES is also invariant to applications of rotation, reflection, and translation, besides maintaining order and monotonicity [43]. It offers no discrepancy in behavior towards varied nature of functions and can be easily generalized. Complexity of algorithm is largely reduced with update schemes of CMA-ES, and thus it offers an extremely prospective mode of maximization in fitting function landscapes [44, 45]. The implementation details of CMA-ES are given in the appendix.The CMA-ES is powerful and performs well. However, better results can be obtained by distributed covariant matrix adaptation evolution strategy (dCMA-ES) with multiple subpopulations and proper migration strategy [46]. Particularly, the population model of dCMA-ES divides the large population into multiple small demes. These demes evolve independently of each other for a certain number of generations, and then a number of individuals are migrated from one deme to another. The dCMA-ES preserves diversity in the population through multiple demes, while increasing the selection pressure through periodic migration [47].The migration operator includes four parameters: (a) migration topology that defines the topology of the connections between demes, (b) the migration rate (the fraction of the population that migrates) that controls how many individuals migrate, (c) migration interval that affects the frequency of migrations, and (d) migration policy that selects emigrants and replaces existing individuals with incoming migrants [48].In this paper, the ring topology is selected as the migration scheme [47]. In other words, individuals are transferred between directionally adjacent demes. The migration policy is that the best individuals (the larger volume value) are selected as migrants to replace the worst individuals at the receiving demes. The migration interval is set to 20 generations to permit the demes to partially converge prior to migration. The migration rate is set to 5 percent to provide a reasonable selection pressure. The final remaining experimental parameters directly related to the dCMA-ES are the deme size and deme count.The distributed covariant matrix adaptation evolution strategy DIRECT (dCMA-ES-DIRECT) algorithm is illustrated in Algorithm 4. The dCMA-ES-DIRECT begins with generating the initial mean, covariance matrix, and step size. Subsequently, the new individual is randomly generated by sampling from a multivariate normal distribution and is checked whether the constraints are satisfied by the DIRECT algorithm. If the constraints are not satisfied, resampling is performed until the generated individual becomes good. The fitness, i.e., the volume of the hyperbox in this paper, is calculated. The next step involves storing the individual with the maximum volume. Then, the migration operator is performed. Next, the update schemes by equations (A.2)–(A.6) in the appendix are applied to equip the individuals for the next generation.Algorithm 4: The proposed distributed covariant matrix adaptation evolution strategy DIRECT (dCMA-ES-DIRECT) algorithm.
(0)Parameter setting:
Set the default strategy parameters according to Appendix A.
(1)Initialization:
(1.1) Set the lower and upper bounds of the design space,xl and xu, the best volume vb=0, the maximum number of the outer iterations gmax, and the threshold level fc.
(1.2) Set the deme sizeds, the deme count dc, the migration rate migr, the migration interval migi.
(1.3) Set the initial covariance matrixC0j, the step-size σ0j, the evolution paths pσj and pcj, j=1,…,dc (details in the appendix).
(1.4) Generate the initial meanm0j uniformly and randomly from the range xl,xl,xu,xu, j=1,…,dc.
(2) Forg=0,…,gmax do:
(2.1) Forj=1,…,dc do:
(2.1.1) Fori=1,…,ds do:
(2.1.1.1) Setk=j−1×ds+i, then generate a random number xkbound from normal distribution Nmgj,σgj2Cgj. Based on equation (5), xklow=xk1bound,…,xkmbound and xkup=xk,m+1bound,…,xk,2mbound.
(2.1.1.2) Evaluatefmin,k=minx∈Dxklow,xkupfx by Algorithm 2 with parameter setting Kmax=⌊100ln2m⌋.
(2.1.1.3) If all constraints are satisfied, i.e.,xl≤xklow≤xkup≤xu, and fmin,k≥fc, then calculate the volume by equation (6), i.e., vk=vxklow,xkup; Otherwise go to Step 2.1.1.1.
(2.1.2) Setk1=j−1×ds+1 and k2=j×ds. Sort ds individuals xk1bound,…,xk2bound based on the value of the volume from largest to smallest. If vs1>vb, set vb=vs1, xlow=xs1low and xup=xs1up, where si denotes the index of ith ranked individual.
(2.2)Migrate:
(2.2.1) Ifg+1 is divisible by migi do:
(2.2.1.1) Forj=1,…,dc do:
(2.2.1.1.1) Setk1=j−1×ds+1 and k2=j×ds. If j=dc, then set k3=1 and k4=ds; Otherwise set k3=j×ds+1 and k4=j+1×ds.
(2.2.1.1.2) Select⌊migr×ds⌋ best individual out of xk1bound,…,xk2bound to replace ⌊migr×ds⌋ worst individual out of xk3bound,…,xk4bound.
(2.3) Forj=1,…,dc do:
(2.3.1)Recombination:
(2.3.1.1) Updatemg+1,j, according to equation (A.3) in the appendix.
(2.3.2)Step-size control
(2.3.2.1) Updatepσ,j, according to equation (A.4) in the appendix.
(2.3.2.2) Updateσg+1,j, according to equation (A.5) in the appendix.
(2.3.3)Covariance matrix adaptation
(2.3.3.1) Updatepc,j, according to equation (A.6) in the appendix.
(2.3.3.2) UpdateCg+1,j, according to equation (A.7) in the appendix.
(3)Verification:
(3.1) Evaluatefmin=minx∈Dxlow,xupfx by Algorithm 2 with the parameter setting Kmax=⌊1000ln2m⌋.
(3.2) Iffmin≥fc, output the individual xlow,xup; Otherwise go to Step 1.Finally, thedCMA-ES-DIRECT verifies whether the best individual meets the constraints by the DIRECT algorithm with a relatively large maximum number of iterations. If the constraints are not met, the dCMA-ES-DIRECT goes back to Step 1. Otherwise, the dCMA-ES-DIRECT outputs the best individual xlow,xup.


### 4.3. Characteristics of the Proposed Approach

As these two metaheuristic algorithms are insensitive to initial solution, derivative-free and global search [49], the DIRECT algorithm with a large enough iteration number converges to the global optimum with a high accuracy [31]; the proposed methods which combine metaheuristic algorithms with the DIRECT algorithm have the following good characteristics:(1)
The quality of the hyperbox obtained finally is not affected by the initial solution, except that the computational time may increase with improper starting designs.(2)
The proposed methods have the advantage of not getting trapped in local optimum and have great possibility to reach the globally maximum solution hyperbox.(3)
Because of the discrete nature of the performance function evaluations in the DIRECT algorithm, the proposed methods can be used for both analytically known and black-box performance functions.(4)
The proposed methods guarantee that any point selected within the hyperbox obtained finally satisfies the performance criterion as long as the performance function is continuous.In the majority of practical engineering problems, the continuity of the performance function is easily satisfied, no matter it is analytically known or a black box. Therefore, the proposed approach has strong applicability and is valuable to practical engineering applications.From Algorithms 3 and 4, we see that the DIRECT algorithm with a relatively small maximum number of iterations (i.e.,Kmax=⌊100ln2m⌋) is adopted in the iterative optimization phase, while the DIRECT algorithm with a relatively large maximum number of iterations (i.e., Kmax=⌊1000ln2m⌋) is adopted in the verification phase. This setting reduces the computation time while ensuring that any point selected within the obtained hyperbox satisfies the performance criterion.Due to the stochastic nature of metaheuristic algorithms, it is important to show the robustness of the proposed approaches. Therefore, as adopted in [15], the proposed approach runs 20 times, and the best solution among the 20 runs is used in practice. In addition, the boxplot containing the minimum, the maximum, the median, and the first and third quartiles of the 20 runs is shown.


## 4.1. Simulated Annealing DIRECT Algorithm (SA-DIRECT)

Simulated annealing (SA) proposed by Kirkpatrick et al. [38] is based on the similarity between the solid annealing process and solving global optimization problems. It is a generic probabilistic meta-algorithm for the global optimization problem. SA has been applied to a wide range of problems especially in those cases where traditional optimization techniques have shown poor performances or simply have failed [39].Simulated annealing DIRECT (SA-DIRECT) algorithm is illustrated in Algorithm3. Initially, the SA-DIRECT sets the design space and builds an initial solution for problem (7). More specifically, we generate a solution x0low,x0up, where x0low and x0up are stochastically generated from the range xl,xu. Then, we evaluate minx∈Dx0low,x0upfx by the DIRECT algorithm in order to check whether this generated solution meets constraints. If the constraint conditions, i.e., xl≤x0low≤x0up≤xu andminx∈Dx0low,x0upfx≥fc, are met, we calculate the associated value of the fitness function (i.e., the volume expressed by equation (6)). Otherwise, we discard this generated solution and generate a new initial solution.Algorithm 3: The proposed simulated annealing DIRECT (SA-DIRECT) algorithm.
(0) Parameter setting:
Set the initial temperatureT0, the cooling ratio β, the maximum number of the inner iterations lmax and the maximum number of the outer iterations tmax.
(1) Initialization:
(1.1) Set the lower and upper bounds of design space,xl and xu, and the threshold level fc.
(1.2) Generate randomlyx0low and x0up from the range xl,xu.
(1.3) Evaluatefmin=minx∈Dx0low,x0upfx by Algorithm 2 with Kmax=⌊50ln2m⌋.
(1.4) If all the constraints are satisfied. i.e.,xl≤x0low≤x0up≤xu and fmin≥fc, then calculate the volume by equation (6) and store it in v0, i.e., v0=vx0low,x0up; Otherwise, go to Step 1.2.
(1.5) Setxlow=xclow=x0low, xup=xcup=x0up, vb=vc=v0, and T=T0.
(2) Fort=1,…,tmaxdo:
(2.1) Setl=1.
(2.2) Whilel≤lmaxdo:
(2.2.1) Generate a new solutionxllow,xlup in the neighborhood of xclow,xcup, i.e., xllow=xclow+xu−xl⊙2rl−13, xlup=xcup+xu−xl⊙2ru−13, where rl and ru are randomly and uniformly generated from range 0,1m, and ⊙ denotes componentwise multiplication.
(2.2.2) Evaluatefmin=minx∈Dxllow,xlupfx by Algorithm 2 with the parameter setting Kmax=⌊100ln2m⌋.
(2.2.3) If constraintsxl≤xllow≤xlup≤xu and fmin≥fc are all met, then evaluate the volume vl by equation (6), i.e., vl=vxllow,xlup; Otherwise go to Step 2.2.1.
(2.2.4) EvaluateΔv=vl−vc, if Δv≥0, then set xclow=xllow and xcup=xlup; If vl−vb>0, then set xlow=xllow, xup=xlup.
(2.2.5) IfΔv<0, then generate a random number r in 0,1 and evaluate P=expΔv/T; If P>r, then set xclow=xllow, xcup=xlup, vc=vl.
(2.2.6) Setl=l+1.
(2.3) SetT=βT.
(3) Verification:
(3.1) Evaluatefmin=minx∈Dxlow,xupfx by Algorithm 2 with the parameter setting Kmax=⌊1000ln2m⌋.
(3.2) Iffmin≥fc, output the best solution xlow,xup; Otherwise go to Step 1.Next, during the process of optimization, in each iteration, a neighbor solution to the current solutionxllow,xlup is generated according to a predefined neighborhood structure, and minx∈Dxllow,xlupfx is evaluated by the DIRECT algorithm. If the condition minx∈Dxllow,xlupfx≥fc is not satisfied, this neighbor solution is rejected. Otherwise, the volume vxllow,xlup is calculated by equation (6). The improving move (the volume of the neighbor is larger than that of the current solution) is always accepted, whilst worse neighbor is accepted with a certain probability determined by the Boltzmann probability [40], expΔv/T>r, where Δv is the difference between the volume of the current solution and the generated neighbor, and r is randomly generated from the range 0,1. Moreover, T is a parameter (called the temperature) which periodically lowers during the search process according to the temperature updating rule. The temperature updating rule (as adopted in [41]) is T=βtT0, where T0 is the initial temperature, β is the cooling ratio, and t is the number of times the temperature has been lowered. The halting criterion of SA is “reaching the maximum number of times the temperature could be lowered.”Finally, the SA-DIRECT verifies whether the best solution meets the constraints by the DIRECT algorithm with a relatively large maximum number of iterations. If the constraints are not met, the SA-DIRECT goes back to Step 1. Otherwise, the SA-DIRECT outputs the solution.In this paper, the maximum number of the iterations at each temperaturelmax=200 and the maximum number of times the temperature could be lowered tmax=1000.


## 4.2. Distributed Covariance Matrix Adaptation Evolution Strategy DIRECT Algorithm (dCMA-ES-DIRECT)

The evolution strategy (ES) is developed as a powerful tool for numerical optimization tasks [42]. Covariant matrix adaptation evolution strategy (CMA-ES) acts as an improved robust form of evolution strategy [43]. The main feature of the CMA-ES is the ability of being invariant to landscape transformations and scaling modulation. The CMA-ES is also invariant to applications of rotation, reflection, and translation, besides maintaining order and monotonicity [43]. It offers no discrepancy in behavior towards varied nature of functions and can be easily generalized. Complexity of algorithm is largely reduced with update schemes of CMA-ES, and thus it offers an extremely prospective mode of maximization in fitting function landscapes [44, 45]. The implementation details of CMA-ES are given in the appendix.The CMA-ES is powerful and performs well. However, better results can be obtained by distributed covariant matrix adaptation evolution strategy (dCMA-ES) with multiple subpopulations and proper migration strategy [46]. Particularly, the population model of dCMA-ES divides the large population into multiple small demes. These demes evolve independently of each other for a certain number of generations, and then a number of individuals are migrated from one deme to another. The dCMA-ES preserves diversity in the population through multiple demes, while increasing the selection pressure through periodic migration [47].The migration operator includes four parameters: (a) migration topology that defines the topology of the connections between demes, (b) the migration rate (the fraction of the population that migrates) that controls how many individuals migrate, (c) migration interval that affects the frequency of migrations, and (d) migration policy that selects emigrants and replaces existing individuals with incoming migrants [48].In this paper, the ring topology is selected as the migration scheme [47]. In other words, individuals are transferred between directionally adjacent demes. The migration policy is that the best individuals (the larger volume value) are selected as migrants to replace the worst individuals at the receiving demes. The migration interval is set to 20 generations to permit the demes to partially converge prior to migration. The migration rate is set to 5 percent to provide a reasonable selection pressure. The final remaining experimental parameters directly related to the dCMA-ES are the deme size and deme count.The distributed covariant matrix adaptation evolution strategy DIRECT (dCMA-ES-DIRECT) algorithm is illustrated in Algorithm 4. The dCMA-ES-DIRECT begins with generating the initial mean, covariance matrix, and step size. Subsequently, the new individual is randomly generated by sampling from a multivariate normal distribution and is checked whether the constraints are satisfied by the DIRECT algorithm. If the constraints are not satisfied, resampling is performed until the generated individual becomes good. The fitness, i.e., the volume of the hyperbox in this paper, is calculated. The next step involves storing the individual with the maximum volume. Then, the migration operator is performed. Next, the update schemes by equations (A.2)–(A.6) in the appendix are applied to equip the individuals for the next generation.Algorithm 4: The proposed distributed covariant matrix adaptation evolution strategy DIRECT (dCMA-ES-DIRECT) algorithm.
(0)Parameter setting:
Set the default strategy parameters according to Appendix A.
(1)Initialization:
(1.1) Set the lower and upper bounds of the design space,xl and xu, the best volume vb=0, the maximum number of the outer iterations gmax, and the threshold level fc.
(1.2) Set the deme sizeds, the deme count dc, the migration rate migr, the migration interval migi.
(1.3) Set the initial covariance matrixC0j, the step-size σ0j, the evolution paths pσj and pcj, j=1,…,dc (details in the appendix).
(1.4) Generate the initial meanm0j uniformly and randomly from the range xl,xl,xu,xu, j=1,…,dc.
(2) Forg=0,…,gmax do:
(2.1) Forj=1,…,dc do:
(2.1.1) Fori=1,…,ds do:
(2.1.1.1) Setk=j−1×ds+i, then generate a random number xkbound from normal distribution Nmgj,σgj2Cgj. Based on equation (5), xklow=xk1bound,…,xkmbound and xkup=xk,m+1bound,…,xk,2mbound.
(2.1.1.2) Evaluatefmin,k=minx∈Dxklow,xkupfx by Algorithm 2 with parameter setting Kmax=⌊100ln2m⌋.
(2.1.1.3) If all constraints are satisfied, i.e.,xl≤xklow≤xkup≤xu, and fmin,k≥fc, then calculate the volume by equation (6), i.e., vk=vxklow,xkup; Otherwise go to Step 2.1.1.1.
(2.1.2) Setk1=j−1×ds+1 and k2=j×ds. Sort ds individuals xk1bound,…,xk2bound based on the value of the volume from largest to smallest. If vs1>vb, set vb=vs1, xlow=xs1low and xup=xs1up, where si denotes the index of ith ranked individual.
(2.2)Migrate:
(2.2.1) Ifg+1 is divisible by migi do:
(2.2.1.1) Forj=1,…,dc do:
(2.2.1.1.1) Setk1=j−1×ds+1 and k2=j×ds. If j=dc, then set k3=1 and k4=ds; Otherwise set k3=j×ds+1 and k4=j+1×ds.
(2.2.1.1.2) Select⌊migr×ds⌋ best individual out of xk1bound,…,xk2bound to replace ⌊migr×ds⌋ worst individual out of xk3bound,…,xk4bound.
(2.3) Forj=1,…,dc do:
(2.3.1)Recombination:
(2.3.1.1) Updatemg+1,j, according to equation (A.3) in the appendix.
(2.3.2)Step-size control
(2.3.2.1) Updatepσ,j, according to equation (A.4) in the appendix.
(2.3.2.2) Updateσg+1,j, according to equation (A.5) in the appendix.
(2.3.3)Covariance matrix adaptation
(2.3.3.1) Updatepc,j, according to equation (A.6) in the appendix.
(2.3.3.2) UpdateCg+1,j, according to equation (A.7) in the appendix.
(3)Verification:
(3.1) Evaluatefmin=minx∈Dxlow,xupfx by Algorithm 2 with the parameter setting Kmax=⌊1000ln2m⌋.
(3.2) Iffmin≥fc, output the individual xlow,xup; Otherwise go to Step 1.Finally, thedCMA-ES-DIRECT verifies whether the best individual meets the constraints by the DIRECT algorithm with a relatively large maximum number of iterations. If the constraints are not met, the dCMA-ES-DIRECT goes back to Step 1. Otherwise, the dCMA-ES-DIRECT outputs the best individual xlow,xup.


## 4.3. Characteristics of the Proposed Approach

As these two metaheuristic algorithms are insensitive to initial solution, derivative-free and global search [49], the DIRECT algorithm with a large enough iteration number converges to the global optimum with a high accuracy [31]; the proposed methods which combine metaheuristic algorithms with the DIRECT algorithm have the following good characteristics:(1)
The quality of the hyperbox obtained finally is not affected by the initial solution, except that the computational time may increase with improper starting designs.(2)
The proposed methods have the advantage of not getting trapped in local optimum and have great possibility to reach the globally maximum solution hyperbox.(3)
Because of the discrete nature of the performance function evaluations in the DIRECT algorithm, the proposed methods can be used for both analytically known and black-box performance functions.(4)
The proposed methods guarantee that any point selected within the hyperbox obtained finally satisfies the performance criterion as long as the performance function is continuous.In the majority of practical engineering problems, the continuity of the performance function is easily satisfied, no matter it is analytically known or a black box. Therefore, the proposed approach has strong applicability and is valuable to practical engineering applications.From Algorithms 3 and 4, we see that the DIRECT algorithm with a relatively small maximum number of iterations (i.e.,Kmax=⌊100ln2m⌋) is adopted in the iterative optimization phase, while the DIRECT algorithm with a relatively large maximum number of iterations (i.e., Kmax=⌊1000ln2m⌋) is adopted in the verification phase. This setting reduces the computation time while ensuring that any point selected within the obtained hyperbox satisfies the performance criterion.Due to the stochastic nature of metaheuristic algorithms, it is important to show the robustness of the proposed approaches. Therefore, as adopted in [15], the proposed approach runs 20 times, and the best solution among the 20 runs is used in practice. In addition, the boxplot containing the minimum, the maximum, the median, and the first and third quartiles of the 20 runs is shown.


## 5. Numerical Examples and Comparisons

A stochastic approach (we denote this method by “GHZ” hereafter) that computes the solution hyperbox has been discussed in [11]. An approach (we denote this method by “CES-IA” hereafter) based on the use of cellular evolutionary strategies (CES) and interval arithmetic (IA) has been proposed in [15].In order to compare the proposed SA-DIRECT anddCMA-ES-DIRECT methods with the existing GHZ and CES-IA methods, the following two numerical examples are considered:(i)
Example 1 studies multiple performance functions. Since each performance function is monotone, the exact solution of the optimization problem (7) exists and can be considered globally optimal.(ii)
Example 2 studies a performance function which has multiple local minima. Unfortunately, the exact solution of the optimization problem (7) does not exist.In addition, for comparison purpose, the solutions of thedCMA-ES-IA method which combines the distributed covariance matrix adaptation evolution strategy (dCMA-ES) with interval arithmetic (IA) is also shown.


### 5.1. Example 1: Multiple Performance Functions

Example 1 comes from [11], including multiple performance functions, where xl=0,0 and xu=4,4 are the lower and upper bounds of the design space, respectively, and(14)f1xf2xf3xf4xf5xf6x=−18x1−14x2−417x1−217x212x1−12x212x1+13x213x1+23x2−x1+32x2,f1cf2cf3cf4cf5cf6c=−1−1−111−1.Sincefix is monotone, as the definition fmin,ixlow,xup=defminDxlow,xupfix, we have(15)fmin,1xlow,xupfmin,2xlow,xupfmin,3xlow,xupfmin,4xlow,xupfmin,5xlow,xupfmin,6xlow,xup=−18x1up−14x2up−417x1up−217x2up12x1low−12x2up12x1low+13x2low13x1low+23x2low−x1up+32x2low.That is,fmin,ixlow,xup has an explicit expression, for i=1,…,6.Therefore, the optimization problem (7) can be restated as follows:(16)maxxlow,xup∏i=12xiup−xilowsubject tofmin,ixlow,xup≥fic,i=1,…,6.Obviously, (16) is an optimization problem with inequality constraints. The exact solution can be found by the Lagrange multiplier method [50] and has the value tabulated in the second column of Table 1.Table 1
.The exact, GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT solutions of Example 1.
Exact
GHZ
CES-IA
dCMA-ES-IA
SA-DIRECT
dCMA-ES-DIRECT
T
0
=
100
T
0
=
150
T
0
=
150
T
0
=
200
d
s
=
8
d
s
=
8
β
=
0.94
β
=
0.94
β
=
0.98
β
=
0.98
d
c
=
1
d
c
=
3
x
1
low
1.1800
1.1740
1.1804
1.1826
1.1805
1.1876
1.1744
1.1957
1.1790
1.1799
x
1
up
2.8500
2.7350
2.8196
2.8392
2.8432
2.8298
2.8547
2.8111
2.8473
2.8453
x
2
low
1.2300
1.2350
1.2298
1.2262
1.2294
1.2206
1.2384
1.2081
1.2316
1.2302
x
2
up
2.5800
2.6330
2.5902
2.5802
2.5770
2.5850
2.5727
2.5943
2.5763
2.5773
Volume
2.2545
2.1823
2.2396
2.2432
2.2408
2.2406
2.2419
2.2398
2.2436
2.2436
Error (%)
—
3.20
0.66
0.50
0.60
0.62
0.56
0.65
0.48
0.48The result in [11] obtained by the GHZ method is listed in the third column of Table 1. The CES-IA method adopts the same parameter setting as in [15] and the best solution among its 20 runs is shown in the fourth column of Table 1. The best solution among 20 runs of the dCMA-ES-IA method is listed in the fifth column of Table 1. The row entitled “Volume” contains the volume of the obtained hyperbox, and the row entitled “Error” contains the relative error between the volume and the exact volume.In the proposed SA-DIRECT method, the initial temperatureT0 is set to 100, 150, and 200, and the cooling ratio β is set to 0.94 and 0.98. In the proposed dCMA-ES-DIRECT method, two sets of deme size and count values are chosen: 8 and 1 and 8 and 3, respectively. The best solutions among their respective 20 runs of the SA-DIRECT and CMA-ES-DIRECT methods are presented in Table 1. A visualization of results of these methods in Table 1 is shown in Figure 4. Moreover, we magnify the region A in Figure 4 and show it in Figure 5.Figure 4
Example 1: the gray region is the complete solution space. The regions are obtained by the following methods: the exact (red), GHZ (green), CES-IA (yellow),dCMA-ES-IA (cyan), SA-DIRECT with T0=100andβ=0.94 (blue), SA-DIRECT with T0=150andβ=0.94 (gold), SA-DIRECT with T0=150andβ=0.98 (pink), SA-DIRECT with T0=200andβ=0.98 (orange), dCMA-ES-DIRECT with ds=8anddc=1 (magenta), and dCMA-ES-DIRECT with ds=8anddc=3 (black).Figure 5
Magnification of region A in Figure4.According to Table1, the volumes of the obtained hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are approximately the same as the exact value, whereas the volume of the obtained hyperbox by the GHZ method has a somewhat large error. Moreover, we can see that the SA-DIRECT method is insensitive to parameter variations of simulated annealing.From Figure4, we observe that the obtained hyperboxes by the SA-DIRECT and dCMA-ES-DIRECT methods are in close proximity to the exact solution. From Figure 5, it can be seen that the green lines exceed the gray region, which means that the hyperbox obtained by the GHZ method contains some design points which are not within the complete solution space, and therefore is not a solution hyperbox. However, the hyperboxes obtained by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods all locate within the complete solution space and therefore are solution hyperboxes.Consequently, the hyperboxes obtained by the CES-IA,dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods not only are close to the exact hyperbox but also guarantee that any point selected within those hyperboxes satisfies the performance criterion.Furthermore, we sketch the boxplots using the volume of the hyperbox of 20 runs, which are presented in Figure6. From Figure 6, we observe that, in each case of parameter setting, the volumes of the obtained hyperboxes are all above 2.2000 after eliminating outliers, and the range of the volumes is relatively small; therefore, the two proposed approaches are robust.Figure 6
Example 1: the boxplots of the volume of the obtained hyperbox. (a) SA-DIRECT method and (b)dCMA-ES-DIRECT method.
(a)
(b)


### 5.2. Numerical Example 2: Multiple Local Minima

The maximization of the solution hyperbox is considered in case of the Michalewics function:(17)fx=−∑i=12sinxisinixi2π20.Stripinis et al. [51] have pointed out that the Michalewics function has two local minima and one global minimum. In this example, the lower and upper bounds of design space are xl=1,1 and xu=π,π, respectively, and the threshold value is fc=−1.5.Unfortunately, the bound-constrained optimization problemminDxlow,xupfx has no explicit expressions; therefore, we cannot obtain the exact solution for the optimization problem (7) in this example. However, minDxlow,xupfx can be solved by the FMINCON function of MATLAB. For comparison purpose, an approach (we denote this method by “SA-FMINCON” hereafter) which combines simulated annealing and FMINCON function is presented to solve problem (7).The best solutions of their respective 20 runs of the GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON methods are shown in Table 2. A visualization of these results is shown in Figure 7. According to Table 2 and Figure 7, the hyperbox obtained by the SA-FMINCON method is the maximum and the hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods are much larger than those by the GHZ, CES-IA, and dCMA-ES-IA methods. Besides, we see that the SA-DIRECT method is relatively insensitive to the initial parameters of simulated annealing.Table 2
The best GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON solutions of Example 2.
GHZ
CES-IA
dCMA-ES-IA
SA-DIRECT
dCMA-ES-DIRECT
SA-FMINCON
T
0
=
100
T
0
=
150
T
0
=
150
T
0
=
200
d
s
=
8
d
s
=
8
β
=
0.94
β
=
0.94
β
=
0.98
β
=
0.98
d
c
=
1
d
c
=
3
x
1
low
1.5095
1.2388
1.0001
1.0002
1.0007
1.0007
1.0007
1.0000
1.0000
1.0009
x
1
up
2.6406
2.2659
1.9619
3.1403
3.1406
3.1414
3.1410
3.1416
3.1415
3.1412
x
2
low
1.6520
1.8628
1.8756
1.6622
1.6625
1.6633
1.6625
1.6620
1.6620
1.0037
x
2
up
3.1415
2.5521
2.5673
3.1410
3.1414
3.1414
3.1408
3.1416
3.1416
3.1415
Volume
1.6849
0.7080
0.7426
3.1649
3.1648
3.1647
3.1640
3.1687
3.1685
4.5755Figure 7
Example 2: the regions obtained by the following methods: the GHZ (green), CES-IA (yellow),dCMA-ES-IA (cyan), SA-DIRECT with T0=100andβ=0.94 (blue), dCMA-ES-DIRECT with ds=8anddc=1 (magenta), and SA-FMINCON (red).The values of the performance functionfx on the obtained hyperboxes by the GHZ, CES-IA, dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON methods are shown in Figure 8. We can see that the values of the performance function fx on the hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are all above the performance criterion while some values of the performance function fx on those by the GHZ and SA-FMINCON methods are below the performance criterion. This indicates that the obtained hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods only include good design points while those by the GHZ and SA-FMINCON methods include some bad design points.Figure 8
The values of the performance function f(x) on the hyperboxes obtained by the following methods: (a) GHZ, (b) CES-IA, (c) dCMA-ES-IA, (d) SA-DIRECT with T0 = 100, β = 0.94 , (e) SA-DIRECT with T0 = 150, β = 0.94, (f) SA-DIRECT with T0 = 150, β = 0.98, (g) SA-DIRECT with T0 = 200, β = 0.98, (h) dCMA-ES-DIRECT with ds = 8, dc = 1, (i) dCMA-ES-DIRECT with ds = 8, dc = 3, and (j) SA-fmincon.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The boxplots of the volume of the hyperbox of 20 runs are given in Figure9. From Figure 9, we see that, in each case of parameter setting, the volumes of the obtained hyperboxes are all above 3.1400 after eliminating outliers, and the range of the volumes is relatively small; therefore, the robustness of the proposed approaches are verified.Figure 9
Example 2: the boxplots of the volume of the obtained hyperbox. (a) SA-DIRECT method and (b)dCMA-ES-DIRECT method.
(a)
(b)


### 5.3. Discussion

The CES-IA,dCMA-ES-IA, and SA-FMINCON methods require an analytically known performance function, while GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods can be used for both analytically known and black-box performance functions.The hyperboxes obtained by the CES-IA anddCMA-ES-IA methods in Example 1 are very close to those obtained by the proposed SA-DIRECT and dCMA-ES-DIRECT methods while the hyperboxes obtained by the CES-IA and dCMA-ES-IA methods in Example 2 are much smaller than those obtained by the proposed SA-DIRECT and dCMA-ES-DIRECT methods. This is mainly caused by the dependency problem of the interval arithmetic. In the CES-IA and dCMA-ES-IA methods, the range of the performance function inside the generated hyperbox is calculated by the interval arithmetic. In Example 1, each design parameter appears only once in the performance function fix; thus, the interval arithmetic can determine the range of the performance function very accurately. However, in Example 2, the performance function fx is unable to be described by the output interval directly, instead the Taylor interval extension [52] is used. In this case, the design parameters occur several times in the calculation, and each occurrence is taken independently. Consequently, the range of the performance function inside the generated hyperbox is overestimated and the actually solution hyperbox is discarded. Rocco et al. [15] have pointed out that the main drawback of the interval arithmetic adopted in the CES-IA method is that results can be overestimated due to the dependency problem, and overestimation cannot be estimated and could cause that no solution hyperbox is found.Numerical examples show that the hyperboxes obtained by the GHZ method may include some bad design points. This is mainly because the GHZ method uses Monte Carlo sampling to estimate locations of good and bad designs. Monte Carlo sampling cannot guarantee that each design point in the obtained hyperbox is good for even a large number of samples, especially when the performance function has multiple local minima. In addition, numerical examples show that the hyperboxes by the GHZ method are smaller than those by the proposed SA-DIRECT anddCMA-ES-DIRECT methods. This is because the GHZ method removes bad sample points sequentially in only one order. The order, however, has an effect on the result. Zimmermann and von Hoessle [10] have concluded that the GHZ method is not necessarily globally optimal, even not locally optimal.In Example 2, although the hyperbox obtained by the SA-FMINCON method is larger than those by the proposed SA-DIRECT anddCMA-ES-DIRECT methods, it includes some bad design points, and therefore is not a solution hyperbox. This demonstrates that it is improper to use the FMINCON function of MATLAB to solve the optimization problem (7) especially when the performance function has multiple local minima.From the “quantitative” viewpoint, the solution hyperboxes obtained by the proposed SA-DIRECT anddCMA-ES-DIRECT methods are nearly the same as the exact result in Example 1, and the solution hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods are much larger than those obtained by the CES-IA and GHZ methods in Example 2 without exact result. From the “feasible” viewpoint, no matter Example 1 or Example 2, any design point selected within the solution hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods satisfies the performance criterion. Therefore, numerical examples prove that the proposed global search is effective.Moreover, from Figures6 and 9, we observe that the SA-DIRECT and dCMA-ES-DIRECT methods are robust. As the SA-DIRECT method with T0=150andβ=0.98 has a relatively large maximum value and a relatively small interquartile range, we therefore set T0=150 and β=0.98 in the remainder of paper. Similarly, the d-CMA-ES-DIRECT method with ds=3+⌊3+ln2m⌋(i.e., ds=8 in Examples 1 and 2) and dc=3 has a relatively large maximum value and a relatively small interquartile range; we therefore set ds=3+⌊3+ln2m⌋ and dc=3 in the rest of the paper.


## 5.1. Example 1: Multiple Performance Functions

Example 1 comes from [11], including multiple performance functions, where xl=0,0 and xu=4,4 are the lower and upper bounds of the design space, respectively, and(14)f1xf2xf3xf4xf5xf6x=−18x1−14x2−417x1−217x212x1−12x212x1+13x213x1+23x2−x1+32x2,f1cf2cf3cf4cf5cf6c=−1−1−111−1.Sincefix is monotone, as the definition fmin,ixlow,xup=defminDxlow,xupfix, we have(15)fmin,1xlow,xupfmin,2xlow,xupfmin,3xlow,xupfmin,4xlow,xupfmin,5xlow,xupfmin,6xlow,xup=−18x1up−14x2up−417x1up−217x2up12x1low−12x2up12x1low+13x2low13x1low+23x2low−x1up+32x2low.That is,fmin,ixlow,xup has an explicit expression, for i=1,…,6.Therefore, the optimization problem (7) can be restated as follows:(16)maxxlow,xup∏i=12xiup−xilowsubject tofmin,ixlow,xup≥fic,i=1,…,6.Obviously, (16) is an optimization problem with inequality constraints. The exact solution can be found by the Lagrange multiplier method [50] and has the value tabulated in the second column of Table 1.Table 1
.The exact, GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT solutions of Example 1.
Exact
GHZ
CES-IA
dCMA-ES-IA
SA-DIRECT
dCMA-ES-DIRECT
T
0
=
100
T
0
=
150
T
0
=
150
T
0
=
200
d
s
=
8
d
s
=
8
β
=
0.94
β
=
0.94
β
=
0.98
β
=
0.98
d
c
=
1
d
c
=
3
x
1
low
1.1800
1.1740
1.1804
1.1826
1.1805
1.1876
1.1744
1.1957
1.1790
1.1799
x
1
up
2.8500
2.7350
2.8196
2.8392
2.8432
2.8298
2.8547
2.8111
2.8473
2.8453
x
2
low
1.2300
1.2350
1.2298
1.2262
1.2294
1.2206
1.2384
1.2081
1.2316
1.2302
x
2
up
2.5800
2.6330
2.5902
2.5802
2.5770
2.5850
2.5727
2.5943
2.5763
2.5773
Volume
2.2545
2.1823
2.2396
2.2432
2.2408
2.2406
2.2419
2.2398
2.2436
2.2436
Error (%)
—
3.20
0.66
0.50
0.60
0.62
0.56
0.65
0.48
0.48The result in [11] obtained by the GHZ method is listed in the third column of Table 1. The CES-IA method adopts the same parameter setting as in [15] and the best solution among its 20 runs is shown in the fourth column of Table 1. The best solution among 20 runs of the dCMA-ES-IA method is listed in the fifth column of Table 1. The row entitled “Volume” contains the volume of the obtained hyperbox, and the row entitled “Error” contains the relative error between the volume and the exact volume.In the proposed SA-DIRECT method, the initial temperatureT0 is set to 100, 150, and 200, and the cooling ratio β is set to 0.94 and 0.98. In the proposed dCMA-ES-DIRECT method, two sets of deme size and count values are chosen: 8 and 1 and 8 and 3, respectively. The best solutions among their respective 20 runs of the SA-DIRECT and CMA-ES-DIRECT methods are presented in Table 1. A visualization of results of these methods in Table 1 is shown in Figure 4. Moreover, we magnify the region A in Figure 4 and show it in Figure 5.Figure 4
Example 1: the gray region is the complete solution space. The regions are obtained by the following methods: the exact (red), GHZ (green), CES-IA (yellow),dCMA-ES-IA (cyan), SA-DIRECT with T0=100andβ=0.94 (blue), SA-DIRECT with T0=150andβ=0.94 (gold), SA-DIRECT with T0=150andβ=0.98 (pink), SA-DIRECT with T0=200andβ=0.98 (orange), dCMA-ES-DIRECT with ds=8anddc=1 (magenta), and dCMA-ES-DIRECT with ds=8anddc=3 (black).Figure 5
Magnification of region A in Figure4.According to Table1, the volumes of the obtained hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are approximately the same as the exact value, whereas the volume of the obtained hyperbox by the GHZ method has a somewhat large error. Moreover, we can see that the SA-DIRECT method is insensitive to parameter variations of simulated annealing.From Figure4, we observe that the obtained hyperboxes by the SA-DIRECT and dCMA-ES-DIRECT methods are in close proximity to the exact solution. From Figure 5, it can be seen that the green lines exceed the gray region, which means that the hyperbox obtained by the GHZ method contains some design points which are not within the complete solution space, and therefore is not a solution hyperbox. However, the hyperboxes obtained by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods all locate within the complete solution space and therefore are solution hyperboxes.Consequently, the hyperboxes obtained by the CES-IA,dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods not only are close to the exact hyperbox but also guarantee that any point selected within those hyperboxes satisfies the performance criterion.Furthermore, we sketch the boxplots using the volume of the hyperbox of 20 runs, which are presented in Figure6. From Figure 6, we observe that, in each case of parameter setting, the volumes of the obtained hyperboxes are all above 2.2000 after eliminating outliers, and the range of the volumes is relatively small; therefore, the two proposed approaches are robust.Figure 6
Example 1: the boxplots of the volume of the obtained hyperbox. (a) SA-DIRECT method and (b)dCMA-ES-DIRECT method.
(a)
(b)


## 5.2. Numerical Example 2: Multiple Local Minima

The maximization of the solution hyperbox is considered in case of the Michalewics function:(17)fx=−∑i=12sinxisinixi2π20.Stripinis et al. [51] have pointed out that the Michalewics function has two local minima and one global minimum. In this example, the lower and upper bounds of design space are xl=1,1 and xu=π,π, respectively, and the threshold value is fc=−1.5.Unfortunately, the bound-constrained optimization problemminDxlow,xupfx has no explicit expressions; therefore, we cannot obtain the exact solution for the optimization problem (7) in this example. However, minDxlow,xupfx can be solved by the FMINCON function of MATLAB. For comparison purpose, an approach (we denote this method by “SA-FMINCON” hereafter) which combines simulated annealing and FMINCON function is presented to solve problem (7).The best solutions of their respective 20 runs of the GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON methods are shown in Table 2. A visualization of these results is shown in Figure 7. According to Table 2 and Figure 7, the hyperbox obtained by the SA-FMINCON method is the maximum and the hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods are much larger than those by the GHZ, CES-IA, and dCMA-ES-IA methods. Besides, we see that the SA-DIRECT method is relatively insensitive to the initial parameters of simulated annealing.Table 2
The best GHZ, CES-IA,dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON solutions of Example 2.
GHZ
CES-IA
dCMA-ES-IA
SA-DIRECT
dCMA-ES-DIRECT
SA-FMINCON
T
0
=
100
T
0
=
150
T
0
=
150
T
0
=
200
d
s
=
8
d
s
=
8
β
=
0.94
β
=
0.94
β
=
0.98
β
=
0.98
d
c
=
1
d
c
=
3
x
1
low
1.5095
1.2388
1.0001
1.0002
1.0007
1.0007
1.0007
1.0000
1.0000
1.0009
x
1
up
2.6406
2.2659
1.9619
3.1403
3.1406
3.1414
3.1410
3.1416
3.1415
3.1412
x
2
low
1.6520
1.8628
1.8756
1.6622
1.6625
1.6633
1.6625
1.6620
1.6620
1.0037
x
2
up
3.1415
2.5521
2.5673
3.1410
3.1414
3.1414
3.1408
3.1416
3.1416
3.1415
Volume
1.6849
0.7080
0.7426
3.1649
3.1648
3.1647
3.1640
3.1687
3.1685
4.5755Figure 7
Example 2: the regions obtained by the following methods: the GHZ (green), CES-IA (yellow),dCMA-ES-IA (cyan), SA-DIRECT with T0=100andβ=0.94 (blue), dCMA-ES-DIRECT with ds=8anddc=1 (magenta), and SA-FMINCON (red).The values of the performance functionfx on the obtained hyperboxes by the GHZ, CES-IA, dCMA-ES-IA, SA-DIRECT, dCMA-ES-DIRECT, and SA-FMINCON methods are shown in Figure 8. We can see that the values of the performance function fx on the hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are all above the performance criterion while some values of the performance function fx on those by the GHZ and SA-FMINCON methods are below the performance criterion. This indicates that the obtained hyperboxes by the CES-IA, dCMA-ES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods only include good design points while those by the GHZ and SA-FMINCON methods include some bad design points.Figure 8
The values of the performance function f(x) on the hyperboxes obtained by the following methods: (a) GHZ, (b) CES-IA, (c) dCMA-ES-IA, (d) SA-DIRECT with T0 = 100, β = 0.94 , (e) SA-DIRECT with T0 = 150, β = 0.94, (f) SA-DIRECT with T0 = 150, β = 0.98, (g) SA-DIRECT with T0 = 200, β = 0.98, (h) dCMA-ES-DIRECT with ds = 8, dc = 1, (i) dCMA-ES-DIRECT with ds = 8, dc = 3, and (j) SA-fmincon.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The boxplots of the volume of the hyperbox of 20 runs are given in Figure9. From Figure 9, we see that, in each case of parameter setting, the volumes of the obtained hyperboxes are all above 3.1400 after eliminating outliers, and the range of the volumes is relatively small; therefore, the robustness of the proposed approaches are verified.Figure 9
Example 2: the boxplots of the volume of the obtained hyperbox. (a) SA-DIRECT method and (b)dCMA-ES-DIRECT method.
(a)
(b)


## 5.3. Discussion

The CES-IA,dCMA-ES-IA, and SA-FMINCON methods require an analytically known performance function, while GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods can be used for both analytically known and black-box performance functions.The hyperboxes obtained by the CES-IA anddCMA-ES-IA methods in Example 1 are very close to those obtained by the proposed SA-DIRECT and dCMA-ES-DIRECT methods while the hyperboxes obtained by the CES-IA and dCMA-ES-IA methods in Example 2 are much smaller than those obtained by the proposed SA-DIRECT and dCMA-ES-DIRECT methods. This is mainly caused by the dependency problem of the interval arithmetic. In the CES-IA and dCMA-ES-IA methods, the range of the performance function inside the generated hyperbox is calculated by the interval arithmetic. In Example 1, each design parameter appears only once in the performance function fix; thus, the interval arithmetic can determine the range of the performance function very accurately. However, in Example 2, the performance function fx is unable to be described by the output interval directly, instead the Taylor interval extension [52] is used. In this case, the design parameters occur several times in the calculation, and each occurrence is taken independently. Consequently, the range of the performance function inside the generated hyperbox is overestimated and the actually solution hyperbox is discarded. Rocco et al. [15] have pointed out that the main drawback of the interval arithmetic adopted in the CES-IA method is that results can be overestimated due to the dependency problem, and overestimation cannot be estimated and could cause that no solution hyperbox is found.Numerical examples show that the hyperboxes obtained by the GHZ method may include some bad design points. This is mainly because the GHZ method uses Monte Carlo sampling to estimate locations of good and bad designs. Monte Carlo sampling cannot guarantee that each design point in the obtained hyperbox is good for even a large number of samples, especially when the performance function has multiple local minima. In addition, numerical examples show that the hyperboxes by the GHZ method are smaller than those by the proposed SA-DIRECT anddCMA-ES-DIRECT methods. This is because the GHZ method removes bad sample points sequentially in only one order. The order, however, has an effect on the result. Zimmermann and von Hoessle [10] have concluded that the GHZ method is not necessarily globally optimal, even not locally optimal.In Example 2, although the hyperbox obtained by the SA-FMINCON method is larger than those by the proposed SA-DIRECT anddCMA-ES-DIRECT methods, it includes some bad design points, and therefore is not a solution hyperbox. This demonstrates that it is improper to use the FMINCON function of MATLAB to solve the optimization problem (7) especially when the performance function has multiple local minima.From the “quantitative” viewpoint, the solution hyperboxes obtained by the proposed SA-DIRECT anddCMA-ES-DIRECT methods are nearly the same as the exact result in Example 1, and the solution hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods are much larger than those obtained by the CES-IA and GHZ methods in Example 2 without exact result. From the “feasible” viewpoint, no matter Example 1 or Example 2, any design point selected within the solution hyperboxes obtained by the SA-DIRECT and dCMA-ES-DIRECT methods satisfies the performance criterion. Therefore, numerical examples prove that the proposed global search is effective.Moreover, from Figures6 and 9, we observe that the SA-DIRECT and dCMA-ES-DIRECT methods are robust. As the SA-DIRECT method with T0=150andβ=0.98 has a relatively large maximum value and a relatively small interquartile range, we therefore set T0=150 and β=0.98 in the remainder of paper. Similarly, the d-CMA-ES-DIRECT method with ds=3+⌊3+ln2m⌋(i.e., ds=8 in Examples 1 and 2) and dc=3 has a relatively large maximum value and a relatively small interquartile range; we therefore set ds=3+⌊3+ln2m⌋ and dc=3 in the rest of the paper.


## 6. Case Studies

The first case is the life-support system in a space capsule which has been studied by Rocco et al. [15] and is reinvestigated by the proposed approach. The second case is the power-shift steering transmission control system (PSSTCS) which is a typical mechatronics complex control system with a price of approximately 500,000 USD. Therefore, it is importantly significant to improve the robustness of the PSSTCS.


### 6.1. Life-Support System in a Space Capsule

#### 6.1.1. Life-Support System in a Space Capsule: Reliability Constraint

The complex system (life-support system in a space capsule) in [15], as shown in Figure 10, consists of four components. The reliability of the system is given by(18)Rs=fR=1−R3R¯1R¯42−R¯31−R21−R¯1R¯42,where R=R1,R2,R3,R4 and R¯i=1−Ri for i=1,…,4.Figure 10
Example of the reliability system.Rocco et al. [15] have applied the CES-IA method to find a largest symmetric hyperbox using a specified point as symmetry center, whereby any point selected within this hyperbox satisfies the performance criterion. Particularly, the center point is Rc=0.90,0.90,0.90,0.90, the lower and upper bounds of the design space are Rl=0 and Ru=1, respectively, and the performance criterion is Rs≥0.99. The best solution among the 20 runs of the CES-IA method shown in [15] is given in the fourth column of Table 3 in terms of the intervals of design parameters.Table 3
Comparison I: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8000, 1.0000]
[0.8058, 0.9942]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
R
2
[0.8245, 0.9755]
[0.8068, 0.9932]
[0.8220, 0.9770]
[0.8241, 0.9759]
[0.8246, 0.9754]
R
3
[0.8000, 1.0000]
[0.8065, 0.9935]
[0.8010, 0.9980]
[0.8001, 0.9999]
[0.8000, 1.0000]
R
4
[0.8000, 1.0000]
[0.8059, 0.9941]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
Volume
1.2134e − 03
1.2354e − 03
1.1786e − 03
1.2114e − 03
1.2120e − 03
Error (%)
—
1.81
2.90
0.16
0.12For comparison purpose, the optimization problem considered in this paper is the same as that considered in [15], and it is stated similarly to the optimization problem (7):(19)maxRup∏i=142Riup−Ric,subject tominR∈Rlow,RupfR≥0.99,Rc≤Rup≤Ru,where Rlow=2Rc−Rup.The best solutions among the 20 runs of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown in the third, fifth, and sixth columns of Table 3, respectively. To show how close to the global solution are the results obtained by the proposed methods, the second column of Table 3 presents the exact solution obtained by the Lagrange multiplier method. To show whether the obtained hyperboxes satisfy the performance criterion, the intervals of the performance function Rs on the hyperboxes obtained by the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are, respectively, calculated and listed in Table 4.Table 4
The intervals ofRs on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes, denoted by Rs,bounded.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9900, 1.0000]
[0.9892, 1.0000]
[0.9900, 0.9999]
[0.9900, 1.0000]
[0.9900, 1.0000]Based on Tables3 and 4, the SA-DIRECT and dCMA-ES-DIRECT methods outperform the GHZ and CES-IA methods in which the volume errors are smaller and the minimum values of the system reliability on their obtained hyperboxes are above the reliability criterion (0.99). Although the hyperbox obtained by the GHZ method is maximum, from the third column of Table 4, we can see that it includes some design points whose corresponding reliabilities are below the reliability criterion. Although the hyperbox obtained by the CES-IA method only includes good design points, its volume is smaller than those by the SA-DIRECT and dCMA-ES-DIRECT methods.Besides, the average volumes are, respectively,1.2111×10−3 and 1.2115×10−3 in 20 runs of the SA-DIRECT and dCMA-ES-DIRECT methods, and the standard deviations are, respectively, 1.2184×10−6 and 1.2745×10−6. This implies that the proposed methods are robust.


#### 6.1.2. Life-Support System in a Space Capsule: Additional Cost Constraint

Normally, the design problem also seeks to minimize the cost functionCs:(20)Cs=2∑i=14KiRiαi.In [15], K1=100, K2=100, K3=100, K4=150, and αi=1 for all i. Using the above values and the previous SA-DIRECT ranges of Ri, the corresponding range of Cs is 724.8847,895.1153. Using the exact results of Ri, the corresponding range of Cs is 724.8300,895.1000.Rocco et al. [15] studied the problem to derive the ranges for each Ri subject to Rs≥0.99 and 750≤Cs≤850. We solve the same problem by the proposed methods for comparison. The best solution shown in [15] by the CES-IA method is given in the fourth column of Table 5.Table 5
Comparison II: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8500, 0.9500]
[0.8549, 0.9451]
[0.8482, 0.9517]
[0.8502, 0.9498]
[0.8502, 0.9498]
R
2
[0.8500, 0.9500]
[0.8553, 0.9447]
[0.8519, 0.9480]
[0.8498, 0.9502]
[0.8498, 0.9502]
R
3
[0.8500, 0.9500]
[0.8538, 0.9462]
[0.8522, 0.9477]
[0.8503, 0.9497]
[0.8503, 0.9497]
R
4
[0.8667, 0.9333]
[0.8551, 0.9449]
[0.8652, 0.9347]
[0.8664, 0.9336]
[0.8664, 0.9336]
Volume
6.6667e-05
6.6977e-05
6.6016e-05
6.6656e-05
6.6661e-05
Error (%)
_______
0.46
0.97
0.02
0.02Since the cost functionCs is monotone, the exact solution can be obtained by the Lagrange multiplier method and is shown in the second column of Table 5. The best solutions among the 20 runs of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods are also shown in the third, fifth, and sixth columns of Table 5, respectively. The intervals of Rs and Cs on the obtained hyperbox, denoted by Rs,bounded and Cs,bounded, respectively, are shown in Table 6. From Tables 5 and 6, we can see that the proposed methods are better than the GHZ and CES-IA methods because the volume errors are smaller and the performance criterion is satisfied in the obtained hyperboxes.Table 6
The intervals forRs and Cs, denoted by Rs,bounded and Cs,bounded, respectively, based on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9955, 0.9998]
[0.9958, 0.9998]
[0.9956, 0.9998]
[0.9954, 0.9999]
[0.9958, 0.9999]
C
s
,
bounded
[770.0100, 849.9900]
[769.3231, 850.6769]
[770.0810, 849.9189]
[770.0000, 849.9999]
[770.0000, 849.9999]The relative errors between the exact solution and the best solutions of the proposed SA-DIRECT anddCMA-ES-DIRECT methods are both 0.02. This indicates that the results are globally optimal.


### 6.2. The Power-Shift Steering Transmission Control System (PSSTCS)

The power-shift steering transmission control system (PSSTCS) is a key complex system with multicharacteristics of heavy vehicle to achieve the control of the steering, speed changing, fan driving, and lubricating. The PSSTCS is composed of a hydraulic oil supply system, an integration pump-motor system, a fan control system, an electronic control system, and a hydraulic control system. The hydraulic oil supply system consists of a fill oil and constant pressure system of pressure oil tank, a pump group, and a fill oil system of transmission control and fan control. The function constitutes and the structure principle drawing of PSSTCS are shown in Figures11 and 12, respectively. The PSSTCS is a nonmonotonic coherent system; therefore, its reliability function is not analytically known and is evaluated by a graphical inductive analysis method based on the principles of the decision tree, i.e., the goal-oriented (GO) reliability assessment method, as illustrated in the appendix.Figure 11
Function constitutes of PSSTCS.Figure 12
Structure principle drawing of PSSTCS.As shown in Table7, the number of design parameters is 86. The optimization problem (7) is given as follows:(21)maxRup,Rlow∑i=186logRiup−Rilowsubject tominR∈Rlow,RupfR≥fcRl≤Rlow≤Rup≤Ru,where Rl=0.9990, Ru=1.0000, fc=0.9900, and f is the reliability function and is evaluated by the goal-oriented (GO) reliability assessment method. Note that the volume here is replaced by the log-volume (logarithmic transformation of the volume) in order to calculate conveniently.Table 7
The GHZ, SA-DIRECT, anddCMA-ES-DIRECT solutions.
No. (variable)
Unit
GHZ
SA-DIRECT
dCMA-ES-DIRECT
1
Oil pan
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
2
Filter LF1
[0.9996, 0.9998]
[0.9992, 1.0000]
[0.9994, 1.0000]
3
Filter LF1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
4
Pump P1
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
5
Pump P1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
6
Power of P1, P2, and P3
[0.9996, 0.9998]
[0.9996, 1.0000]
[0.9995, 1.0000]
7
Filter LF2
[0.9995, 0.9997]
[0.9993, 1.0000]
[0.9993, 1.0000]
8
Filter LF2
[0.9995, 0.9998]
[0.9995, 1.0000]
[0.9994, 1.0000]
9
Bypass-valve LF2B
[0.9995, 0.9997]
[0.9992, 1.0000]
[0.9993, 1.0000]
10
Pressure oil tank
[0.9996, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
11
Pump P3
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
12
Pump P2
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
13
One-way valve CV1
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
14
Filter LF3
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9997, 1.0000]
15
One-way valve CV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
16
Bypass-valve LF3B
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
17
Constant pressure valve RV1
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9996, 1.0000]
18
Power P4 and P5
[0.9997, 0.9999]
[0.9999, 1.0000]
[0.9998, 1.0000]
19
Pump P4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
20
Filter LF4
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
21
Constant pressure valve RV2
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
22
Steering wheel signal
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9995, 1.0000]
23
Closed-loop feedback link
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
24
Valve body of SDV1
[0.9996, 0.9998]
[0.9997, 1.0000]
[0.9997, 1.0000]
25
Piston valve block
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9996, 1.0000]
26
Hydraulic cylinder DLU
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9995, 1.0000]
27
Swash plate pump of P5
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]
28
Two-way variable displacement pump P5
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
29
Overflow valve OV1
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
30
One-way valve CV4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9995, 1.0000]
31
Overflow valve OV2
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
32
One-way valve CV3
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
33
Group valve RVG
[0.9997, 0.9999]
[0.9998, 1.0000]
[0.9996, 1.0000]
34
Reversing motor M
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
35
Power of electronic control system
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
36
Control panel
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
37
Panel switch
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
38
Handle signal
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
39
State signal sensor
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
40
Switch D1
[0.9996, 0.9999]
[0.9994, 1.0000]
[0.9996, 1.0000]
41
Switch D2
[0.9996, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
42
Pilot valve PV2
[0.9997, 0.9999]
[0.9998, 1.0000]
[0.9998, 1.0000]
43
Liquid viscous clutch cylinder of left fan
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
44
Active friction plate of left pan
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
45
Friction plate of left pan
[0.9997, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
46
Power of liquid viscous clutch
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]
47
Liquid viscous clutch cylinder of right fan
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9996, 1.0000]
48
Active friction plate of right fan
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
49
Friction plate of right fan
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9996, 1.0000]
50
Overflow valve OV3
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
51
Throttle valve TV1
[0.9995, 0.9998]
[0.9992, 1.0000]
[0.9993, 1.0000]
52
Throttle valve TV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
53
Overflow valve RV3
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9997, 1.0000]
54
Overflow valve RV4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
55
Signal of SDV
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
56
Hand valve SDV
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
57
Electric controllable valve DV1
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
58
Manual hydraulic control valve SV1
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
59
Liquid control valve HV1
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
60
Constant pressure throttle valve EV1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
61
Oil cylinder E1
[0.9995, 0.9998]
[0.9995, 1.0000]
[0.9993, 1.0000]
62
Electric controllable valve DV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
63
Manual hydraulic control valve SV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9995, 1.0000]
64
Liquid control valve HV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
65
Constant pressure throttle valve EV2
[0.9992, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
66
Oil cylinder F1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
67
Electric controllable valve DV3
[0.9996, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
68
Manual hydraulic control valve SV3
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
69
Liquid control valve HV3
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
70
Constant pressure throttle valve EV3
[0.9992, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
71
Oil cylinder E2
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
72
Electric controllable valve DV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
73
Manual hydraulic control valve SV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9995, 1.0000]
74
Liquid control valve HV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
75
Oil cylinder F2
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
76
Electric controllable valve DV5
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
77
Manual hydraulic control valve SV5
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
78
Liquid control valve HV5
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9994, 1.0000]
79
Constant pressure throttle valve EV5
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
80
Oil cylinder E3
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
81
Electric controllable valve DV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
82
Manual hydraulic control valve SV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
83
Liquid control valve HV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
84
Constant pressure throttle valve EV6
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
85
Oil cylinder F3
[0.9995, 0.9998]
[0.9992, 1.0000]
[0.9993, 1.0000]
86
Pilot valve PV1
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]Since the reliability function of the PSSCTS is not analytically known, the CES-IA method is not suitable for this complex system. The best solutions among the 20 runs of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown in the third, fourth, and fifth columns of Table 7, respectively. The log-volumes of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes are listed in Table 8. We see that the log-volumes of the SA-DIRECT and dCMA-ES-DIRECT hyperboxes are all much larger than that of the GHZ hyperbox.Table 8
The log-volumes of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT hyperboxes.
Method
GHZ
SA-DIRECT
dCMA-ES-DIRECT
Log-volume
−686.1204
−646.3142
−646.8860To show whether the obtained hyperboxes satisfy the performance criterion,N sample design points are randomly chosen from each of the obtained hyperbox by using the Latin hypercube sampling; then, the rate of good sample design points (the performance criterion satisfaction probability for the hyperbox) is calculated as follows:(22)rate=1N∑i=1N#fRi≥fc,where #⋅ is the indicator function, i.e., #fRi≥fc=1,fRi≥fc0,fRi<fc.The performance criterion satisfaction probabilities of the obtained hyperboxes by the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown under different sample sizes in Figure 13. It can be observed that the performance criterion satisfaction probabilities for the SA-DIRECT and dCMA-ES-DIRECT hyperboxes are all 1 while those for the GHZ hyperbox are below 1. This implies that the obtained hyperboxes by the SA-DIRECT and dCMA-ES-DIRECT methods are solution hyperboxes. Therefore, our approach is superior to the GHZ method.Figure 13
The performance criterion satisfaction probabilities for the obtained hyperboxes. (a) GHZ method, (b) SA-DIRECT method, and (c)dCMA-ES-DIRECT method.
(a)
(b)
(c)


## 6.1. Life-Support System in a Space Capsule

### 6.1.1. Life-Support System in a Space Capsule: Reliability Constraint

The complex system (life-support system in a space capsule) in [15], as shown in Figure 10, consists of four components. The reliability of the system is given by(18)Rs=fR=1−R3R¯1R¯42−R¯31−R21−R¯1R¯42,where R=R1,R2,R3,R4 and R¯i=1−Ri for i=1,…,4.Figure 10
Example of the reliability system.Rocco et al. [15] have applied the CES-IA method to find a largest symmetric hyperbox using a specified point as symmetry center, whereby any point selected within this hyperbox satisfies the performance criterion. Particularly, the center point is Rc=0.90,0.90,0.90,0.90, the lower and upper bounds of the design space are Rl=0 and Ru=1, respectively, and the performance criterion is Rs≥0.99. The best solution among the 20 runs of the CES-IA method shown in [15] is given in the fourth column of Table 3 in terms of the intervals of design parameters.Table 3
Comparison I: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8000, 1.0000]
[0.8058, 0.9942]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
R
2
[0.8245, 0.9755]
[0.8068, 0.9932]
[0.8220, 0.9770]
[0.8241, 0.9759]
[0.8246, 0.9754]
R
3
[0.8000, 1.0000]
[0.8065, 0.9935]
[0.8010, 0.9980]
[0.8001, 0.9999]
[0.8000, 1.0000]
R
4
[0.8000, 1.0000]
[0.8059, 0.9941]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
Volume
1.2134e − 03
1.2354e − 03
1.1786e − 03
1.2114e − 03
1.2120e − 03
Error (%)
—
1.81
2.90
0.16
0.12For comparison purpose, the optimization problem considered in this paper is the same as that considered in [15], and it is stated similarly to the optimization problem (7):(19)maxRup∏i=142Riup−Ric,subject tominR∈Rlow,RupfR≥0.99,Rc≤Rup≤Ru,where Rlow=2Rc−Rup.The best solutions among the 20 runs of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown in the third, fifth, and sixth columns of Table 3, respectively. To show how close to the global solution are the results obtained by the proposed methods, the second column of Table 3 presents the exact solution obtained by the Lagrange multiplier method. To show whether the obtained hyperboxes satisfy the performance criterion, the intervals of the performance function Rs on the hyperboxes obtained by the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are, respectively, calculated and listed in Table 4.Table 4
The intervals ofRs on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes, denoted by Rs,bounded.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9900, 1.0000]
[0.9892, 1.0000]
[0.9900, 0.9999]
[0.9900, 1.0000]
[0.9900, 1.0000]Based on Tables3 and 4, the SA-DIRECT and dCMA-ES-DIRECT methods outperform the GHZ and CES-IA methods in which the volume errors are smaller and the minimum values of the system reliability on their obtained hyperboxes are above the reliability criterion (0.99). Although the hyperbox obtained by the GHZ method is maximum, from the third column of Table 4, we can see that it includes some design points whose corresponding reliabilities are below the reliability criterion. Although the hyperbox obtained by the CES-IA method only includes good design points, its volume is smaller than those by the SA-DIRECT and dCMA-ES-DIRECT methods.Besides, the average volumes are, respectively,1.2111×10−3 and 1.2115×10−3 in 20 runs of the SA-DIRECT and dCMA-ES-DIRECT methods, and the standard deviations are, respectively, 1.2184×10−6 and 1.2745×10−6. This implies that the proposed methods are robust.


### 6.1.2. Life-Support System in a Space Capsule: Additional Cost Constraint

Normally, the design problem also seeks to minimize the cost functionCs:(20)Cs=2∑i=14KiRiαi.In [15], K1=100, K2=100, K3=100, K4=150, and αi=1 for all i. Using the above values and the previous SA-DIRECT ranges of Ri, the corresponding range of Cs is 724.8847,895.1153. Using the exact results of Ri, the corresponding range of Cs is 724.8300,895.1000.Rocco et al. [15] studied the problem to derive the ranges for each Ri subject to Rs≥0.99 and 750≤Cs≤850. We solve the same problem by the proposed methods for comparison. The best solution shown in [15] by the CES-IA method is given in the fourth column of Table 5.Table 5
Comparison II: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8500, 0.9500]
[0.8549, 0.9451]
[0.8482, 0.9517]
[0.8502, 0.9498]
[0.8502, 0.9498]
R
2
[0.8500, 0.9500]
[0.8553, 0.9447]
[0.8519, 0.9480]
[0.8498, 0.9502]
[0.8498, 0.9502]
R
3
[0.8500, 0.9500]
[0.8538, 0.9462]
[0.8522, 0.9477]
[0.8503, 0.9497]
[0.8503, 0.9497]
R
4
[0.8667, 0.9333]
[0.8551, 0.9449]
[0.8652, 0.9347]
[0.8664, 0.9336]
[0.8664, 0.9336]
Volume
6.6667e-05
6.6977e-05
6.6016e-05
6.6656e-05
6.6661e-05
Error (%)
_______
0.46
0.97
0.02
0.02Since the cost functionCs is monotone, the exact solution can be obtained by the Lagrange multiplier method and is shown in the second column of Table 5. The best solutions among the 20 runs of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods are also shown in the third, fifth, and sixth columns of Table 5, respectively. The intervals of Rs and Cs on the obtained hyperbox, denoted by Rs,bounded and Cs,bounded, respectively, are shown in Table 6. From Tables 5 and 6, we can see that the proposed methods are better than the GHZ and CES-IA methods because the volume errors are smaller and the performance criterion is satisfied in the obtained hyperboxes.Table 6
The intervals forRs and Cs, denoted by Rs,bounded and Cs,bounded, respectively, based on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9955, 0.9998]
[0.9958, 0.9998]
[0.9956, 0.9998]
[0.9954, 0.9999]
[0.9958, 0.9999]
C
s
,
bounded
[770.0100, 849.9900]
[769.3231, 850.6769]
[770.0810, 849.9189]
[770.0000, 849.9999]
[770.0000, 849.9999]The relative errors between the exact solution and the best solutions of the proposed SA-DIRECT anddCMA-ES-DIRECT methods are both 0.02. This indicates that the results are globally optimal.


## 6.1.1. Life-Support System in a Space Capsule: Reliability Constraint

The complex system (life-support system in a space capsule) in [15], as shown in Figure 10, consists of four components. The reliability of the system is given by(18)Rs=fR=1−R3R¯1R¯42−R¯31−R21−R¯1R¯42,where R=R1,R2,R3,R4 and R¯i=1−Ri for i=1,…,4.Figure 10
Example of the reliability system.Rocco et al. [15] have applied the CES-IA method to find a largest symmetric hyperbox using a specified point as symmetry center, whereby any point selected within this hyperbox satisfies the performance criterion. Particularly, the center point is Rc=0.90,0.90,0.90,0.90, the lower and upper bounds of the design space are Rl=0 and Ru=1, respectively, and the performance criterion is Rs≥0.99. The best solution among the 20 runs of the CES-IA method shown in [15] is given in the fourth column of Table 3 in terms of the intervals of design parameters.Table 3
Comparison I: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8000, 1.0000]
[0.8058, 0.9942]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
R
2
[0.8245, 0.9755]
[0.8068, 0.9932]
[0.8220, 0.9770]
[0.8241, 0.9759]
[0.8246, 0.9754]
R
3
[0.8000, 1.0000]
[0.8065, 0.9935]
[0.8010, 0.9980]
[0.8001, 0.9999]
[0.8000, 1.0000]
R
4
[0.8000, 1.0000]
[0.8059, 0.9941]
[0.8010, 0.9980]
[0.8000, 1.0000]
[0.8000, 1.0000]
Volume
1.2134e − 03
1.2354e − 03
1.1786e − 03
1.2114e − 03
1.2120e − 03
Error (%)
—
1.81
2.90
0.16
0.12For comparison purpose, the optimization problem considered in this paper is the same as that considered in [15], and it is stated similarly to the optimization problem (7):(19)maxRup∏i=142Riup−Ric,subject tominR∈Rlow,RupfR≥0.99,Rc≤Rup≤Ru,where Rlow=2Rc−Rup.The best solutions among the 20 runs of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown in the third, fifth, and sixth columns of Table 3, respectively. To show how close to the global solution are the results obtained by the proposed methods, the second column of Table 3 presents the exact solution obtained by the Lagrange multiplier method. To show whether the obtained hyperboxes satisfy the performance criterion, the intervals of the performance function Rs on the hyperboxes obtained by the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT methods are, respectively, calculated and listed in Table 4.Table 4
The intervals ofRs on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes, denoted by Rs,bounded.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9900, 1.0000]
[0.9892, 1.0000]
[0.9900, 0.9999]
[0.9900, 1.0000]
[0.9900, 1.0000]Based on Tables3 and 4, the SA-DIRECT and dCMA-ES-DIRECT methods outperform the GHZ and CES-IA methods in which the volume errors are smaller and the minimum values of the system reliability on their obtained hyperboxes are above the reliability criterion (0.99). Although the hyperbox obtained by the GHZ method is maximum, from the third column of Table 4, we can see that it includes some design points whose corresponding reliabilities are below the reliability criterion. Although the hyperbox obtained by the CES-IA method only includes good design points, its volume is smaller than those by the SA-DIRECT and dCMA-ES-DIRECT methods.Besides, the average volumes are, respectively,1.2111×10−3 and 1.2115×10−3 in 20 runs of the SA-DIRECT and dCMA-ES-DIRECT methods, and the standard deviations are, respectively, 1.2184×10−6 and 1.2745×10−6. This implies that the proposed methods are robust.


## 6.1.2. Life-Support System in a Space Capsule: Additional Cost Constraint

Normally, the design problem also seeks to minimize the cost functionCs:(20)Cs=2∑i=14KiRiαi.In [15], K1=100, K2=100, K3=100, K4=150, and αi=1 for all i. Using the above values and the previous SA-DIRECT ranges of Ri, the corresponding range of Cs is 724.8847,895.1153. Using the exact results of Ri, the corresponding range of Cs is 724.8300,895.1000.Rocco et al. [15] studied the problem to derive the ranges for each Ri subject to Rs≥0.99 and 750≤Cs≤850. We solve the same problem by the proposed methods for comparison. The best solution shown in [15] by the CES-IA method is given in the fourth column of Table 5.Table 5
Comparison II: the exact, GHZ, CES-IA, SA-DIRECT, anddCMA-ES-DIRECT solutions.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
1
[0.8500, 0.9500]
[0.8549, 0.9451]
[0.8482, 0.9517]
[0.8502, 0.9498]
[0.8502, 0.9498]
R
2
[0.8500, 0.9500]
[0.8553, 0.9447]
[0.8519, 0.9480]
[0.8498, 0.9502]
[0.8498, 0.9502]
R
3
[0.8500, 0.9500]
[0.8538, 0.9462]
[0.8522, 0.9477]
[0.8503, 0.9497]
[0.8503, 0.9497]
R
4
[0.8667, 0.9333]
[0.8551, 0.9449]
[0.8652, 0.9347]
[0.8664, 0.9336]
[0.8664, 0.9336]
Volume
6.6667e-05
6.6977e-05
6.6016e-05
6.6656e-05
6.6661e-05
Error (%)
_______
0.46
0.97
0.02
0.02Since the cost functionCs is monotone, the exact solution can be obtained by the Lagrange multiplier method and is shown in the second column of Table 5. The best solutions among the 20 runs of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods are also shown in the third, fifth, and sixth columns of Table 5, respectively. The intervals of Rs and Cs on the obtained hyperbox, denoted by Rs,bounded and Cs,bounded, respectively, are shown in Table 6. From Tables 5 and 6, we can see that the proposed methods are better than the GHZ and CES-IA methods because the volume errors are smaller and the performance criterion is satisfied in the obtained hyperboxes.Table 6
The intervals forRs and Cs, denoted by Rs,bounded and Cs,bounded, respectively, based on the exact, GHZ, CES-IA, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes.
Exact
GHZ
CES-IA
SA-DIRECT
dCMA-ES-DIRECT
R
s
,
bounded
[0.9955, 0.9998]
[0.9958, 0.9998]
[0.9956, 0.9998]
[0.9954, 0.9999]
[0.9958, 0.9999]
C
s
,
bounded
[770.0100, 849.9900]
[769.3231, 850.6769]
[770.0810, 849.9189]
[770.0000, 849.9999]
[770.0000, 849.9999]The relative errors between the exact solution and the best solutions of the proposed SA-DIRECT anddCMA-ES-DIRECT methods are both 0.02. This indicates that the results are globally optimal.


## 6.2. The Power-Shift Steering Transmission Control System (PSSTCS)

The power-shift steering transmission control system (PSSTCS) is a key complex system with multicharacteristics of heavy vehicle to achieve the control of the steering, speed changing, fan driving, and lubricating. The PSSTCS is composed of a hydraulic oil supply system, an integration pump-motor system, a fan control system, an electronic control system, and a hydraulic control system. The hydraulic oil supply system consists of a fill oil and constant pressure system of pressure oil tank, a pump group, and a fill oil system of transmission control and fan control. The function constitutes and the structure principle drawing of PSSTCS are shown in Figures11 and 12, respectively. The PSSTCS is a nonmonotonic coherent system; therefore, its reliability function is not analytically known and is evaluated by a graphical inductive analysis method based on the principles of the decision tree, i.e., the goal-oriented (GO) reliability assessment method, as illustrated in the appendix.Figure 11
Function constitutes of PSSTCS.Figure 12
Structure principle drawing of PSSTCS.As shown in Table7, the number of design parameters is 86. The optimization problem (7) is given as follows:(21)maxRup,Rlow∑i=186logRiup−Rilowsubject tominR∈Rlow,RupfR≥fcRl≤Rlow≤Rup≤Ru,where Rl=0.9990, Ru=1.0000, fc=0.9900, and f is the reliability function and is evaluated by the goal-oriented (GO) reliability assessment method. Note that the volume here is replaced by the log-volume (logarithmic transformation of the volume) in order to calculate conveniently.Table 7
The GHZ, SA-DIRECT, anddCMA-ES-DIRECT solutions.
No. (variable)
Unit
GHZ
SA-DIRECT
dCMA-ES-DIRECT
1
Oil pan
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
2
Filter LF1
[0.9996, 0.9998]
[0.9992, 1.0000]
[0.9994, 1.0000]
3
Filter LF1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
4
Pump P1
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
5
Pump P1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
6
Power of P1, P2, and P3
[0.9996, 0.9998]
[0.9996, 1.0000]
[0.9995, 1.0000]
7
Filter LF2
[0.9995, 0.9997]
[0.9993, 1.0000]
[0.9993, 1.0000]
8
Filter LF2
[0.9995, 0.9998]
[0.9995, 1.0000]
[0.9994, 1.0000]
9
Bypass-valve LF2B
[0.9995, 0.9997]
[0.9992, 1.0000]
[0.9993, 1.0000]
10
Pressure oil tank
[0.9996, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
11
Pump P3
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
12
Pump P2
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
13
One-way valve CV1
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
14
Filter LF3
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9997, 1.0000]
15
One-way valve CV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
16
Bypass-valve LF3B
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
17
Constant pressure valve RV1
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9996, 1.0000]
18
Power P4 and P5
[0.9997, 0.9999]
[0.9999, 1.0000]
[0.9998, 1.0000]
19
Pump P4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
20
Filter LF4
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
21
Constant pressure valve RV2
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
22
Steering wheel signal
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9995, 1.0000]
23
Closed-loop feedback link
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
24
Valve body of SDV1
[0.9996, 0.9998]
[0.9997, 1.0000]
[0.9997, 1.0000]
25
Piston valve block
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9996, 1.0000]
26
Hydraulic cylinder DLU
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9995, 1.0000]
27
Swash plate pump of P5
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]
28
Two-way variable displacement pump P5
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
29
Overflow valve OV1
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
30
One-way valve CV4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9995, 1.0000]
31
Overflow valve OV2
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
32
One-way valve CV3
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
33
Group valve RVG
[0.9997, 0.9999]
[0.9998, 1.0000]
[0.9996, 1.0000]
34
Reversing motor M
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
35
Power of electronic control system
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
36
Control panel
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
37
Panel switch
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
38
Handle signal
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
39
State signal sensor
[0.9990, 1.0000]
[0.9990, 1.0000]
[0.9990, 1.0000]
40
Switch D1
[0.9996, 0.9999]
[0.9994, 1.0000]
[0.9996, 1.0000]
41
Switch D2
[0.9996, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
42
Pilot valve PV2
[0.9997, 0.9999]
[0.9998, 1.0000]
[0.9998, 1.0000]
43
Liquid viscous clutch cylinder of left fan
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
44
Active friction plate of left pan
[0.9996, 0.9999]
[0.9996, 1.0000]
[0.9997, 1.0000]
45
Friction plate of left pan
[0.9997, 0.9999]
[0.9993, 1.0000]
[0.9997, 1.0000]
46
Power of liquid viscous clutch
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]
47
Liquid viscous clutch cylinder of right fan
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9996, 1.0000]
48
Active friction plate of right fan
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
49
Friction plate of right fan
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9996, 1.0000]
50
Overflow valve OV3
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
51
Throttle valve TV1
[0.9995, 0.9998]
[0.9992, 1.0000]
[0.9993, 1.0000]
52
Throttle valve TV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
53
Overflow valve RV3
[0.9996, 0.9999]
[0.9997, 1.0000]
[0.9997, 1.0000]
54
Overflow valve RV4
[0.9996, 0.9999]
[0.9998, 1.0000]
[0.9997, 1.0000]
55
Signal of SDV
[0.9996, 0.9999]
[0.9999, 1.0000]
[0.9997, 1.0000]
56
Hand valve SDV
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9996, 1.0000]
57
Electric controllable valve DV1
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
58
Manual hydraulic control valve SV1
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
59
Liquid control valve HV1
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
60
Constant pressure throttle valve EV1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
61
Oil cylinder E1
[0.9995, 0.9998]
[0.9995, 1.0000]
[0.9993, 1.0000]
62
Electric controllable valve DV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
63
Manual hydraulic control valve SV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9995, 1.0000]
64
Liquid control valve HV2
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
65
Constant pressure throttle valve EV2
[0.9992, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
66
Oil cylinder F1
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
67
Electric controllable valve DV3
[0.9996, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
68
Manual hydraulic control valve SV3
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
69
Liquid control valve HV3
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
70
Constant pressure throttle valve EV3
[0.9992, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
71
Oil cylinder E2
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
72
Electric controllable valve DV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
73
Manual hydraulic control valve SV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9995, 1.0000]
74
Liquid control valve HV4
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
75
Oil cylinder F2
[0.9995, 0.9998]
[0.9990, 1.0000]
[0.9993, 1.0000]
76
Electric controllable valve DV5
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9993, 1.0000]
77
Manual hydraulic control valve SV5
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
78
Liquid control valve HV5
[0.9995, 0.9998]
[0.9991, 1.0000]
[0.9994, 1.0000]
79
Constant pressure throttle valve EV5
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
80
Oil cylinder E3
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
81
Electric controllable valve DV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
82
Manual hydraulic control valve SV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
83
Liquid control valve HV6
[0.9995, 0.9998]
[0.9993, 1.0000]
[0.9993, 1.0000]
84
Constant pressure throttle valve EV6
[0.9995, 0.9998]
[0.9994, 1.0000]
[0.9993, 1.0000]
85
Oil cylinder F3
[0.9995, 0.9998]
[0.9992, 1.0000]
[0.9993, 1.0000]
86
Pilot valve PV1
[0.9996, 0.9999]
[0.9995, 1.0000]
[0.9997, 1.0000]Since the reliability function of the PSSCTS is not analytically known, the CES-IA method is not suitable for this complex system. The best solutions among the 20 runs of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown in the third, fourth, and fifth columns of Table 7, respectively. The log-volumes of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT hyperboxes are listed in Table 8. We see that the log-volumes of the SA-DIRECT and dCMA-ES-DIRECT hyperboxes are all much larger than that of the GHZ hyperbox.Table 8
The log-volumes of the GHZ, SA-DIRECT, anddCMA-ES-DIRECT hyperboxes.
Method
GHZ
SA-DIRECT
dCMA-ES-DIRECT
Log-volume
−686.1204
−646.3142
−646.8860To show whether the obtained hyperboxes satisfy the performance criterion,N sample design points are randomly chosen from each of the obtained hyperbox by using the Latin hypercube sampling; then, the rate of good sample design points (the performance criterion satisfaction probability for the hyperbox) is calculated as follows:(22)rate=1N∑i=1N#fRi≥fc,where #⋅ is the indicator function, i.e., #fRi≥fc=1,fRi≥fc0,fRi<fc.The performance criterion satisfaction probabilities of the obtained hyperboxes by the GHZ, SA-DIRECT, anddCMA-ES-DIRECT methods are shown under different sample sizes in Figure 13. It can be observed that the performance criterion satisfaction probabilities for the SA-DIRECT and dCMA-ES-DIRECT hyperboxes are all 1 while those for the GHZ hyperbox are below 1. This implies that the obtained hyperboxes by the SA-DIRECT and dCMA-ES-DIRECT methods are solution hyperboxes. Therefore, our approach is superior to the GHZ method.Figure 13
The performance criterion satisfaction probabilities for the obtained hyperboxes. (a) GHZ method, (b) SA-DIRECT method, and (c)dCMA-ES-DIRECT method.
(a)
(b)
(c)


## 7. Conclusions

To improve the design robustness, rather than optimizing one single design point, the approach presented in this paper maximizes the region in the design space where all design points meet the required performance goal. Moreover, regions expressed by hyperboxes are considered in order to decouple parameters.To this end, an innovative global approach that combines metaheuristic algorithms and the DIRECT algorithm is proposed to seek a maximum solution hyperbox. The metaheuristic algorithm is used to obtain a maximum hyperbox and the DIRECT algorithm is used as a checking technique to guarantee that the obtained hyperbox is a solution hyperbox. More specifically, the SA-DIRECT anddCMA-ES-DIRECT methods are presented in detail. The results of studies on complex numerical examples and engineering cases have shown that these two methods have better performance than the GHZ and CES-IA methods. Since the DIRECT algorithm only requires evaluating the values of the performance function at the sample points, the proposed approach can be used for both analytically known and black-box performance functions. The performance function is continuous in most engineering problems; therefore, the proposed approach is a powerful tool that engineering designers can use to obtain a maximum solution hyperbox.The optimality of the proposed approach mainly depends on the properties of the metaheuristic algorithms. The two metaheuristic algorithms we adopt are global search and have great possibility to reach the globally optimal solution. Our work can be further improved provided that more advanced global search algorithms are developed.



---

*Source: 3697161-2020-05-11.xml*