A Globally Optimal Robust Design Method for Complex Systems

,


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. e 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. ese 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][6][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.
at is, we focus on finding a maximum box-shaped solution space rather than a single optimum solution.
e box-shaped solution space, representing a hyperbox, can be expressed by intervals for each design parameter and has the following benefits: (i) e 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) e intervals may be used to assess the robustness and sensitivity of the uncertain design parameters by the widths of the associated intervals. (iv) e intervals are independent of each other and may be combined with intervals of other disciplines, enabling distributed design in separate disciplines Two categories of approaches that solve similar problems could be identified in the literature. e first category is based on data-mining and machine-learning techniques [10][11][12][13]. e 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. e quality of the results and the efficiency of this stochastic approach are studied in detail in [11].
is stochastic approach, however, produces a boxshaped 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.
e second category is based on an analytical technique, called interval arithmetic [14]. e 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]. e 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 DI-RECT 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. erefore, 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. e 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. e paper is organized as follows. In Section 2, 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.

Motivation and Problem Statement
2.1. Example Problem. We consider the voltage divider shown in Figure 1 in [15]. e terminal voltage V ab is expressed as V ab � V 0 R/(R + R a ) and the resistance Good designs fulfilling the design goals should satisfy (1) Figure 2 shows the complete solution space (gray region). However, the classical optimum that minimizes R in and maximizes V ab is not robust. e immediate neighborhood of the optimum includes designs that violate condition (1) and thus fails to meet at least one design goal. e 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 for R a 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.
2 Complexity e 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. e choice of intervals is not unique, and the size of the interval for R a depends on the choice of the interval for R. e approach proposed here is to maximize the volume of the solution hyperbox, and the resulting unique intervals are shown in Figure 2.
. , x m ) represent the design parameters or design points or designs, where m is the number of dimensions. e output performance at x is given by 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: with the threshold value f c . e 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.
e design spaceD(x l , x u ) is assumed to be continuous and defined by the lower and upper bounds of the design parameters, given by x l and x u , that is, ..,m . e complete solution space is the set of all good design points in the design space. e 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. e hyperbox D � D(x low ,x up ) can be expressed by that is, the lower bound and upper bound of the hyperbox D(x low , x up ) are x low and x up , respectively. Moreover, let x bound denote the bounds of the hyperbox, namely,   Complexity A hyperbox which only includes good design points is called a solution hyperbox. Simply, the hyperbox which meets the condition f(x) ≥ f c for all x ∈ D(x low , x up ) is the solution hyperbox. Furthermore, f min (x low , x up ) is defined as the global minimization of the performance function on the hyperbox, that is, With these preparations at hand, the problem of searching for a maximum solution hyperbox can be formulated as the following constrained optimization problem: If f(x) is analytically known and min x∈[x low ,x up ] f(x) (i.e., f min (x low , x up )) can be explicitly expressed by x low and x up , 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 f(x) is not analytically known; therefore, it is produced by numerical simulations. Second, min x∈[x low ,x up ] f(x) cannot be explicitly expressed by x low and x up . 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). ere 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 2.
Typically, engineering problems have multiple In the optimization problem formulation (7),

The DIRECT Algorithm
We consider the global optimization problem of the form As mentioned in Section 2, D(x low , is an m-dimensional hyperbox (hyper-rectangle), and f(x) is a performance function.
DIRECT's performance is independent of the scaling of D and the problem is typically scaled as follows: us, the bound-constrained problem (8) can be reduced to the following problem: over the unit hypercube D � [0, 1] m . at is, the search space is transformed to the m-dimensional unit hypercube.

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. e requirement of potential optimality is stated formally in Definition 1.
Definition 1 (potentially optimal hyper-rectangle). Suppose that we have a partition of the unit hypercube into n hyperrectangles. Let c i denote the center point of the ith hyperrectangle, d i denote the distance from the center point to the vertices, ε > 0 be a positive constant, and f min be the current best function value. A hyper-rectangle j is said to be potentially optimal if there exist some K > 0 such that e hyper-rectangle j 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.
(1) Identify the set I of dimensions with the maximum side length. Let δ equal one-third of this maximum side length.
(2) Sample the function at the points c ± δe i for all i ∈ I, where c is the center of rectangle and e i is the ith unit vector.
(3) Divide the rectangle containing c into thirds along the dimensions in I, starting with the dimension with the lowest value of , and continuing to the dimension with the next smallest ω i until we have split on all dimensions i ∈ I.
ALGORITHM 1: Procedure for dividing rectangles.  Complexity position and zeros elsewhere). e strategy in [31] is used to divide the hypercube.
More precisely, we adopt the following rule for subdividing a hypercube. Let (13) be the best of the function values sampled along dimension i. en, 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.

e 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 Section 3.1. ese hyper-rectangles are then sampled and subdivided as described in Section 3.2. e process continues until a prespecified iteration limit is reached. A formal statement of the DIRECT algorithm is given in Algorithm 2.
In order to better illustrate how the DIRECT algorithm works, Figure 3shows 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.

Convergence.
e 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]. is 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. at 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][36][37]].

The Proposed Approach
As metaheuristic algorithms are general and do not rely on mathematical properties for application and the DI-RECT 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 Section 2.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. ey 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.

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 Algorithm 3. Initially, the SA-DIRECT sets the design space and builds an initial solution for problem (7). More specifically, we generate a solution ( 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.
Next, during the process of optimization, in each iteration, a neighbor solution to the current solution (x low l , x up l ) is generated according to a predefined neighborhood structure, and min x∈D(x low is calculated by equation (6). e 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. e temperature updating rule (as adopted in [41]) is T � β t T 0 , where T 0 is the initial temperature, β is the cooling ratio, and t is the number of times the temperature has been lowered. e 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 temperature l max � 200 and the maximum number of times the temperature could be lowered t max � 1000.

Distributed Covariance Matrix Adaptation Evolution Strategy DIRECT Algorithm (dCMA-ES-DIRECT).
e 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]. e main feature of the CMA-ES is the ability of being invariant to landscape transformations and scaling modulation. e 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]. e implementation details of CMA-ES are given in the appendix. e 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. ese 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. e dCMA-ES preserves diversity in the population through multiple demes, while increasing the selection pressure through periodic migration [47]. e 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. e migration policy is that the best individuals (the larger volume value) are selected as migrants to replace the worst individuals at the receiving demes. e migration interval is set to 20 generations to permit the demes to partially converge prior to migration. e migration rate is set to 5 percent to provide a reasonable selection pressure. e final remaining experimental parameters directly related to the dCMA-ES are the deme size and deme count. e distributed covariant matrix adaptation evolution strategy DIRECT (dCMA-ES-DIRECT) algorithm is illustrated x u and f min ≥ f c , then calculate the volume by equation (6) and store it where r l and r u are randomly and uniformly generated from range [0, 1] m , and ⊙ denotes componentwise multiplication.
l ≤ x u and f min ≥ f c are all met, then evaluate the volume v l by equation (6) in Algorithm 4. e 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. e fitness, i.e., the volume of the hyperbox in this paper, is calculated. e next step involves storing the individual with the maximum volume. en, 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. Finally, the dCMA-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 (x low , x up ).

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) e quality of the hyperbox obtained finally is not affected by the initial solution, except that the computational time may increase with improper starting designs. (1.1) Set the lower and upper bounds of the design space, x l and x u , the best volume v b � 0, the maximum number of the outer iterations g max , and the threshold level f c . (1.2) Set the deme size d s , the deme count d c , the migration rate mig r , the migration interval mig i . (1.3) Set the initial covariance matrix C 0j , the step-size σ 0j , the evolution paths p σj and p cj , j � 1, . . . , d c (details in the appendix).
(1.4) Generate the initial mean m 0j uniformly and randomly from the range (2) For g � 0, . . . , g max do: (2.1) For j � 1, . . . , d c do:   8 Complexity proposed methods can be used for both analytically known and black-box performance functions. (4) e 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. erefore, 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., K max � ⌊100 ln 2m⌋) is adopted in the iterative optimization phase, while the DIRECT algorithm with a relatively large maximum number of iterations (i.e., K max � ⌊1000 ln 2m⌋) is adopted in the verification phase.
is 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.
erefore, 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.

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 and dCMA-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 the dCMA-ES-IA method which combines the distributed covariance matrix adaptation evolution strategy (dCMA-ES) with interval arithmetic (IA) is also shown.

Example 1: Multiple Performance
Functions. Example 1 comes from [11], including multiple performance functions, where x l � (0, 0) and x u � (4, 4) are the lower and upper bounds of the design space, respectively, and Since f i (x) is monotone, as the definition at is, f min,i (x low , x up ) has an explicit expression, for i � 1, . . . , 6. erefore, the optimization problem (7) can be restated as follows: Obviously, (16) is an optimization problem with inequality constraints. e exact solution can be found by the Lagrange multiplier method [50] and has the value tabulated in the second column of Table 1. e result in [11] obtained by the GHZ method is listed in the third column of Table 1. e 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. e best solution among 20 runs of the dCMA-ES-IA method is listed in the fifth column of Table 1. e 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 temperature T 0 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. e 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.
According to Table 1, 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 Figure 4, 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 Figure 6. 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.

Numerical Example 2: Multiple Local Minima.
e maximization of the solution hyperbox is considered in case of the Michalewics function: 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 x l � (1, 1) and x u � (π, π), respectively, and the threshold value is f c � − 1.5.
Unfortunately, the bound-constrained optimization problem min D(x low ,x up ) f(x) has no explicit expressions; therefore, we cannot obtain the exact solution for the optimization problem (7) 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.   Figure 5: Magnification of region A in Figure 4.   e boxplots of the volume of the hyperbox of 20 runs are given in Figure 9. From Figure 9, we see that, in each case       DIRECT methods. is 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 f i (x); thus, the interval arithmetic can determine the range of the performance function very accurately. However, in Example 2, the performance function f(x) 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.
is 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 and dCMA-ES-DIRECT methods.
is is because the GHZ method removes bad sample points sequentially in only one order. e 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 and dCMA-ES-DIRECT methods, it includes some bad design points, and therefore is not a solution hyperbox. is 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 and dCMA-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. erefore, numerical examples prove that the proposed global search is effective.
Moreover, from Figures 6 and 9, we observe that the SA-DIRECT and dCMA-ES-DIRECT methods are robust. As the SA-DIRECT method with T 0 � 150 and β � 0.98 has a relatively large maximum value and a relatively small interquartile range, we therefore set T 0 � 150 and β � 0.98 in the remainder of paper. Similarly, the d-CMA-ES-DIRECT method with d s � 3 + ⌊3 + ln 2 m⌋(i.e., d s � 8 in Examples 1 and 2) and d c � 3 has a relatively large maximum value and a relatively small interquartile range; we therefore set d s � 3 + ⌊3 + ln 2 m⌋ and d c � 3 in the rest of the paper.

Case Studies
e 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. e 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. erefore, it is importantly significant to improve the robustness of the PSSTCS.

Life-Support System in a Space Capsule: Reliability
Constraint. e complex system (life-support system in a space capsule) in [15], as shown in Figure 10, consists of four components. e reliability of the system is given by 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 R c � (0.90, 0.90, 0.90, 0.90), the lower and upper bounds of the design space are R l � 0 and R u � 1, respectively, and the performance criterion is R s ≥ 0.99. e 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.
For 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): where R low � 2R c − R up . e best solutions among the 20 runs of the GHZ, SA-DIRECT, and dCMA-ES-DIRECTmethods 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 R s 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.
Based on Tables 3 and 4, the SA-DIRECTand 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. is implies that the proposed methods are robust.

Life-Support System in a Space Capsule: Additional
Cost Constraint. Normally, the design problem also seeks to minimize the cost function C s : In [15], K 1 � 100, K 2 � 100, K 3 � 100, K 4 � 150, and α i � 1 for all i. Using the above values and the previous SA-DIRECT ranges of R i , the corresponding range of C s is Rocco et al. [15] studied the problem to derive the ranges for each R i subject to R s ≥ 0.99 and 750 ≤ C s ≤ 850. We solve the same problem by the proposed methods for comparison. e best solution shown in [15] by the CES-IA method is given in the fourth column of Table 5.
Since the cost function C s is monotone, the exact solution can be obtained by the Lagrange multiplier method and is shown in the second column of Table 5. e 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. e intervals of R s and C s on the obtained hyperbox, denoted by R s,bounded and C s,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. e relative errors between the exact solution and the best solutions of the proposed SA-DIRECT and dCMA-ES-DIRECT methods are both 0.02.
is indicates that the results are globally optimal.

e Power-Shift Steering Transmission Control System (PSSTCS).
e 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. e 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. e 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. e function constitutes and the structure principle drawing of PSSTCS are shown in Figures 11 and 12, respectively. e 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.
As shown in Table 7, the number of design parameters is 86. e optimization problem (7) is given as follows: where R l � 0.9990, R u � 1.0000, f c � 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.
Since the reliability function of the PSSCTS is not analytically known, the CES-IA method is not suitable for this complex system. e best solutions among the 20 runs of the GHZ, SA-DIRECT, and dCMA-ES-DIRECT methods are shown in the third, fourth, and fifth columns of Table 7, respectively. e logvolumes 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.
To 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:  Figure 10: Example of the reliability system.

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. e 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 and dCMA-ES-DIRECT methods are presented in detail. e 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. e 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. e optimality of the proposed approach mainly depends on the properties of the metaheuristic algorithms. e 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. e default strategy parameters are μ, w i�1,...,u , c σ , d σ , c c , μ cov , and c cov . Default strategy parameters values are given in [53] and shown in Table 9. Hansen [53] has pointed out these default parameters are in particular chosen to be a robust setting, and thus are not recommended to be changed.

B. The Goal-Oriented (GO) Reliability
Assessment Method e goal-oriented (GO) reliability assessment method for the power-shift steering transmission control system (PSSTCS) is given in [54] in detail. e function GO operators, logical GO operators, and auxiliary GO operators are selected to describe the unit itself, its logical relationships, and auxiliary GO operation in the PSSTCS, respectively, as presented in Tables 10 and 11. e GO model of the PSSTCS is developed by using the signal flows to connect the above GO operators, as shown in Figure 14.
e signal flow 129 is the system's reliability output.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.