A set strategy approach for multidisciplinary robust design optimization under interval uncertainty

Abstract

Uncertainties widely exist in complex engineering systems. Robust design is one of the most used method for designing under uncertainty and has been gaining more attention. For the wide range of uncertainties, this article proposes a multidisciplinary robust design optimization method based on the set strategy. In this method, a robust design model that utilizes the maximum variation analysis is developed for uncertainty analysis. Then, a set strategy–based approach is employed to build a system optimization model, which is used to coordinate the coupling variables between full autonomy subsystems and obtains a new design space. Finally, the system robust optimal solution and the optimal robust design space are obtained through the sequential optimization, which provide a direction for the subsequent analysis. Two mathematics examples and the speed reducer design problem are taken to verify the validity and accuracy of the proposed method. A practical engineering problem, namely, air cooling battery thermal management system design problem, is successfully solved by the proposed method.

Introduction

Multidisciplinary design optimization (MDO) as the key technology of future cutting-edge engineering design can significantly shorten the design cycle through parallel design, which is the main development trend of modern complex engineering systems design.^1–5 In addition, there are numerous uncertainties in the whole life cycle of complex systems, such as randomness of loads, geometrical dimensions caused by machining errors, external conditions of operating conditions, and uncertainty of establishing simulation models. These uncertainties and their coupling effects between different disciplines lead to the fact that the final design is difficult to achieve in an optimal state and cannot be applied directly to engineering practice. MDO under uncertainty has become one of the most challenging problems.^6–12

To effectively reduce the impact of uncertainties on the quality of the product, the robust design has aroused the concern of many scholars. Taguchi proposed the robust design ideas in the 1970s for the first time. ¹³ After years of development, some sound robust design methods are maturing and are applied to solve MDO problems under uncertainty in various engineering fields.^14,15 Du and Chen ¹⁶ utilized efficient methods for uncertainty analysis in MDO and proposed a multidisciplinary robust design optimization (MRDO) procedure based on system uncertainty analysis (SUA) and concurrent subsystem uncertainty analysis (CSSUA) methods. ¹⁷ Mourelatos and Liang ¹⁸ reviewed the reassessing reliability and robustness in engineering design and proposed a new unifying formulation. Liu et al. ¹⁹ proposed a probabilistic analytical target cascading (PATC) formulation and addressed several implementation issues, for instance, the representation of probabilistic design targets, matching responses and linking variables under uncertainty, and coordination strategies. Xu et al. ²⁰ employed a time-varying sensitivity method to establish an MRDO model. Wang et al. ²¹ proposed a unified framework of combining reliability with robustness for integrated optimization under uncertainty. Giassi et al. ²² described an innovative multidisciplinary robust design approach which engaged concurrent work and high precision optimization results. Zaman and Mahadevan ²³ proposed formulations and algorithms for the perspective of system robustness design optimization of multidisciplinary systems under uncertainty and developed a single-loop approach for the design optimization which does not require any coupled multidisciplinary uncertainty propagation analysis. Hu et al. ²⁴ presented a new approximation-assisted multi-objective collaborative robust optimization approach under interval uncertainty. More excellent research about MRDO under uncertainty can be found in Li and Azarm, ²⁵ and Xia et al. ²⁶

MDO methods can be divided into single-level and multilevel methods. Each single-level MDO method only has one optimizer, and the engineering system does not need to be decomposed—for example, multidisciplinary feasible (MDF), ²⁷ individual disciplinary feasible (IDF), ²⁸ simultaneous analysis and design (SAND)^29,30 methods. As to multilevel MDO methods, each level has an optimizer, such as collaborative optimization (CO), ³¹ concurrent subspace optimization (CSSO), ³² bi-level integrated system synthesis (BLISS) ³³ methods, and analytical target cascading (ATC). ³⁴

The above MDO methods only discussed the point design in which only an optimal scheme is obtained in each design cycle. For a typical design problem, after the initial analysis, the engineers may choose the most promising option. If the promising option cannot meet the requirements, engineers will go back to the initial analysis and then get a new solution, through this design process until a satisfactory solution is found. This lengthy design process greatly limits the efficiency of design optimization.

For the design of complex systems, due to the performance, cost of the product and customer requirements may change. A single-point design is difficult to adapt to changing requirements. Another option is the set-based design, which provides flexibility in the design process and gradually reduces design space. The set-based design is a design methodology that seeks to offer improvements over the traditional point-based design approach.^35–37

Madhavan et al. ³⁸ proposed a set-based multiscale and multidisciplinary design method to help distributed designers manage interdependencies by exchanging targets and Pareto sets of solutions. Avigad and Moshaiov ³⁹ provided a novel approach to perform concept selection in multi-objective design. A computing visualization system for multidisciplinary design and analysis named DataView has been developed by Wang. ⁴⁰ Wang and Terpenny ⁴¹ combined set-based design and fuzzy design trade-off strategy into evolutionary synthesis for addressing the preliminary engineering design problems. Hannapel ⁴² developed a set-based MDO algorithm. It has been successfully applied to a ship design analysis. Nahm and Ishikawa ⁴³ proposed a novel space-based design methodology to obtain a ranged set solution that can satisfy changing sets of performance requirements for preliminary engineering design. Malak et al. ⁴⁴ studied a set theory method to effectively depict the inaccuracy of design and used it in the conceptual design of automobile transmission.

In a set-based design, designers begin with a broad design set values for design variables, and then gradually reduce the sets to obtain more information. Therefore, a new MRDO method based on set strategy is proposed in this article. This method can not only improve the efficiency of robust design but also obtain the optimal robust design space for a complex multidisciplinary system. The remainder of this article is structured as follows. In “Technical background” section, reviews of robust design, MRDO, and the set strategy are introduced. The proposed MRDO method based on set strategy is elaborated in “The developed method” section. Four examples are given in “Tests and results” section. “Conclusion” section concludes the article.

Technical background

In this section, several definitions and terminologies are provided. First, a robust design method is discussed. Second, an MDO formulation with a two-subsystem is given. Finally, a set-based strategy is discussed.

Robust design model for maximum variation analysis

From the content of robust design, the robust design can be divided into two aspects: the robustness of the objective function and the robustness of the constraints. The robustness of objective function discusses how to reduce the variation of the objective function. The robustness of constraints studies how to ensure that the design solution falls within the range of feasible space under the impact of uncertainties.

A general formulation of single objective optimization problems with uncertainty is given

$$\begin{array}{c}minf(\mathbf{x},\tilde{\mathbf{p}})\hfill \\ s.t.g(\mathbf{x},\tilde{\mathbf{p}})⩽0\hfill \\ \mathbf{x}\in [{\mathbf{x}}^{\mathbf{L}},{\mathbf{x}}^{\mathbf{U}}]\hfill \end{array}$$

(1)

where the x is a vector that contains N-element design variables, that is, x = [x₁, …, x_N]. p = [p₁, …, p_M] is a vector of M-element uncertainty parameters. x^L and x^U are the lower and upper bounds of x, respectively.

It is difficult to obtain the accurate distribution of uncertainty in the practical engineering design. The uncertainties in this article are quantified by a series of intervals, namely, the uncertainty of a design scheme x, which can be expressed by a mapping from x’s tolerance region in p space, that is

$$\tilde{\mathbf{p}}\in [{\mathbf{p}}_0-\mathbf{\Delta }\mathbf{p},{\mathbf{p}}_0+\mathbf{\Delta }\mathbf{p}]$$

where $\tilde{\mathbf{p}}$ is the uncertainty parameter vector with the variation range Δp around its nominal value p₀. It is important to note that the region of uncertainty given here is symmetric about the nominal value, which is not necessary. The proposed method in the article can solve the optimization problem with uncertainty in both variables and parameters (x and p). For convenience, the article uses $\tilde{\mathbf{p}}$ to represent uncertainty including both variables and parameters (x and p), while the exams only contain uncertainty in variables.

We assume that an initial design scheme (x₀, p₀) has been identified, due to the existence of uncertainties in design variables and design parameters. The corresponding change of constraint can be expressed as follows

$$\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})=g({\mathbf{x}}_0,\tilde{\mathbf{p}})-g({\mathbf{x}}_0,{\mathbf{p}}_0)$$

(2)

where the $\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ is the variation of the inequality constraint, $g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ is a variable constraint, and $g({\mathbf{x}}_0,{\mathbf{p}}_0)$ is the original constraint.

$\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ is not a fixed value; its value will accompany a different Δp and have a different value. But we can obtain the maximum $\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ , when Δp changes in the hypercube space. Because Δp is the only variable of $\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ , $\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})$ does not contain any constraint. It can be seen as a simple unconstrained optimization problem.

A robust indicator of constraint is used to judge the feasibility of a design scheme. The inequality constraints of equation (1) can be rewritten as

$$g({\mathbf{x}}_0,\tilde{\mathbf{p}})+\mathrm{\Delta }g({\mathbf{x}}_0,\tilde{\mathbf{p}})⩽0$$

(3)

Similarly, the variation of the equality constraint and objective function can be expressed as follows

$$\mathrm{\Delta }h({\mathbf{x}}_0,\tilde{\mathbf{p}})=h({\mathbf{x}}_0,\tilde{\mathbf{p}})-h({\mathbf{x}}_0,{\mathbf{p}}_0)$$

(4)

$$\mathrm{\Delta }f({\mathbf{x}}_0,\tilde{\mathbf{p}})=f({\mathbf{x}}_0,\tilde{\mathbf{p}})-f({\mathbf{x}}_0,{\mathbf{p}}_0)$$

(5)

The robustness of the objective function refers to how to determine the value of the design variables so that the variation of the objective function is as small as possible due to uncertainties. In practical engineering design, designers may not seek the most robust design of the objective function but require the objective function to be optimal within a prescribed range. Therefore, the objective function’s robustness index can be converted into constraints which are added to the original optimization problem. A robust design optimization model is given as follows

$$\begin{array}{c}minf(\mathbf{x},{\mathbf{p}}_0)\hfill \\ s.t.g_i(\mathbf{x},{\mathbf{p}}_0)+\mathrm{\Delta }{g}_i⩽0,i=1,\cdots ,n\hfill \\ \mathrm{\Delta }{h}_j⩽\mathrm{\Delta }{h}^B,j=1,\cdots ,m\hfill \\ \mathrm{\Delta }f⩽\mathrm{\Delta }{f}^B\hfill \\ \mathbf{x}\in [{\mathbf{x}}^{\mathbf{L}},{\mathbf{x}}^{\mathbf{U}}]\hfill \end{array}$$

(6)

where f is the objective function, g is the vector of inequality constraints, and Δg, Δh, and Δf represent the maximum variation vector of inequality constraint, equality constraint, and the objective function, respectively. Δh ^B and Δf ^B are the bound values given by designers. The n denotes the number of inequality constraints and m denotes the number of equality constraints. x^L and x^U are the lower and upper bounds of x, respectively. Δg, Δh, and Δf are calculated by the following formulas

$$\begin{array}{c}\mathrm{\Delta }f=max\left\{|f({\mathbf{x}}_0,\tilde{\mathbf{p}})-f({\mathbf{x}}_0,{\mathbf{p}}_0)|\right\}\}\hfill \\ \mathrm{\Delta }{g}_i=max\{g_i({\mathbf{x}}_0,\tilde{\mathbf{p}})-g_i({\mathbf{x}}_0,{\mathbf{p}}_0)\}\hfill \\ \mathrm{\Delta }{h}_j=max\{h_j({\mathbf{x}}_0,\tilde{\mathbf{p}})-h_j({\mathbf{x}}_0,{\mathbf{p}}_0)\}\hfill \\ \tilde{\mathbf{p}}\in [{\mathbf{p}}_0-\mathrm{\Delta }\mathbf{p},{\mathbf{p}}_0+\mathrm{\Delta }\mathbf{p}]\hfill \end{array}$$

(7)

In the robust design of the objective function, if the goal is not to find the most robust solution of the objective function, it is necessary to prescribe the maximum variation bound of the objective function. However, Δh and Δf cannot be set arbitrarily. Because the values are set too large, it means that the equality constraints and the objective function can be in a large range of arbitrary fluctuations, so that the robust design will lose its original meaning. If the setting is too small, it is very likely that there is no solution, that is, in any solution design, solution of the equation variations or the objective function variations is greater than the set value, so it is extremely important to set a reasonable Δh and Δf. For the sake of expression, the changes of all uncertainties in the following text are expressed in terms of p. The interval uncertainty robust design optimization model is established (Figure 1); it is a two-level optimization mathematical model—the top box of the figure is a conventional optimization problem and the below box of the figure is an optimization problem to solve the robustness. The simplified solving process can be expressed with three steps. (1) A design result is obtained in the conventional optimization, and it is passed to the lower level optimization. (2) Design results are obtained in step (1) as a constant in the lower level optimization. Therefore, maximum variation values of equality and inequality are solved by genetic algorithm (GA). (3) The maximum variation values in step (2) are taken into the upper level optimization; a new design result is solved and it is passed to the lower level optimization. Looping is done until the upper level optimization converges to the optimal solution.

Figure 1.

Robust design model under uncertainty.

MDO formulation

The general formulation of MDO problem is defined by

$$\begin{array}{c}minf(\mathbf{x},\mathbf{y})\hfill \\ s.t.g(\mathbf{x},\mathbf{y})⩽0\hfill \\ \mathbf{y}=Y(\mathbf{x},{\mathbf{y}}^c)\hfill \\ \mathbf{x}\in [{\mathbf{x}}^{\mathbf{L}},{\mathbf{x}}^{\mathbf{U}}]\hfill \end{array}$$

(8)

where g and f represent the constraints and the objective function of system, respectively. Y is the multidisciplinary analysis (MDA) function to obtain coupled state variable y.

Take a two-subsystem problem as an example (Figure 2).

Figure 2.

MDO model of two subsystems.

Subsystems 1 and 2 are defined as follows

$$\begin{array}{c}min{f}_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{12},{\mathbf{y}}_{21}^{\mathrm{c}})\hfill \\ s.t.{\mathbf{y}}_{12}=Y_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{21},{\mathbf{y}}_{21}^{\mathrm{c}})\hfill \\ g_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{12},{\mathbf{y}}_{21}^{\mathrm{c}})⩽0\hfill \end{array}$$

(9)

$$\begin{array}{c}min{f}_{\mathrm{sub}(2)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{21},{\mathbf{y}}_{12}^{\mathrm{c}})\hfill \\ s.t.{\mathbf{y}}_{21}=Y_{\mathrm{sub}(2)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{12},{\mathbf{y}}_{12}^{\mathrm{c}})\hfill \\ g_{\mathrm{sub}(2)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{21},{\mathbf{y}}_{12}^{\mathrm{c}})⩽0\hfill \end{array}$$

(10)

where g_sub(i) and f_sub(i) represent the constraints and objective function of ith subsystem, respectively. Y_sub(i) denotes the analysis function associated with coupled state variable y . y _ij is the coupled state variable vector which is the output from subsystem i and input to subsystem j. x_sub(i) is a vector of subsystem variables. y _ij ^c is the constant vector of the coupled state variable y _ij , which is used to coordinate the coupling between subsystems. For the subsystem 1, the value of the coupling variable is output to the subsystem 2 via the MDA. However, before the MDA of the subsystem 1, the coupling variable y₂₁ must be known and the coupling variable y₂₁ is obtained with the coupling variable y₁₂ as a prerequisite. This reflects the coupling relationships between disciplines; the key to solving MDO problems is to coordinate the coupling relationships between subsystems.

Set strategy–based design optimization

To improve the efficiency of optimization, we adopt a set strategy to develop a novel MRDO method. The developed method is different from the previous MDO method, which can not only effectively decouple but also rapidly reduce the optimization search region and shorten the repetitive work of the optimization process. For ease of understanding, Figure 3 describes the difference between a point-based design and a set-based design.

Figure 3.

Point-based design and set-based design.

Figure 3 illustrates boundary change of a simple design problem, with only two design variables. The light gray shade part of the figure shows the original design space. The left part of the figure describes the traditional point-based design. First, an initial point is given; then based on a gradient method, the search direction is constructed; or through an intelligent algorithm, a better solution is obtained. Finally, the optimal value of the objective function can be acquired by iterating sequentially. The right part of the figure shows the new design space, which is the smaller, purple shaded area. As described in “Introduction” section, one of the set-based design principles is to describe the design variables by changing the set of variables in the design process. The boundary between the design variables of the original MDO problem is defined in equation (8). Since the boundaries of design variables are getting smaller in the optimization process, this manner gives designers a way to decide. The purpose of the new MDO algorithm is to use a set-based approach to design and motivate design optimization. This means the design space is changed through the optimization process. In this way, we can avoid repeated analysis work and get the best design space. The results provide designers with the direction for subsequent analysis and decision making.

If a multidisciplinary system includes N subsystems, the objective functions of these subsystems may vary in order of magnitude. In addition, some target functions may vary greatly in the overall design space, while other changes are less. Therefore, the evaluation and comparison of different objective functions are contained. At first, those objective functions need to be normalized.

The main steps to reduce the design space are as follows:

• First, normalize objective function value of subsystems.

The objective function values of subsystem i are normalized as follows

$$\begin{array}{c}\multicolumn{1}{c}{f_{\mathrm{sub}(i)}^{\mathrm{Nom}}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}})}\\ \multicolumn{1}{c}{=\frac{f_{\mathrm{sub}(i)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}})-f_{\mathrm{sub}(i)}^{min}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}})}{f_{\mathrm{sub}(i)}^{max}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}})-f_{\mathrm{sub}(i)}^{min}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}})}}\end{array}$$

(11)

where $f_{\mathrm{sub}(i)}^{\mathrm{Nom}}()$ is the objective function value after normalization, and $f_{\mathrm{sub}(i)}^{max}$ and $f_{\mathrm{sub}(i)}^{min}$ are the maximum and minimum value of the subsystem objective function, respectively. y_sub(i) is the coupling state variable; it contains the state variables of the subsystem i and the coupling variables of other subsystem to subsystem i ( ${\mathbf{y}}_i^c$ ).

After normalization, we can find that the values of $f_{\mathrm{sub}(i)}^{\mathrm{Nom}}()$ are in the range of [0, 1]. When the value is smaller, the subsystem is closer to the minimum value.

• Second, define and calculate subsystem space reduction rate.

Due to the complex coupling of MDO, a series of attempts have been made to reduce the search space. A simple and easy to calculate method is given. For instance, in the first cycle, for subsystem i, subsystem space reduction rate is defined and calculated by

$${\mathrm{\Psi }}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{cycle}}=\frac{|{\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{(}\mathbf{*},\mathbf{cycle}\mathbf{)}}-{\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{start}}|}{{\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{U}}-{\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{L}}}$$

(12)

where ${\mathrm{\Psi }}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{cycle}}$ is the subsystem i space reduction rate, ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{(}\mathbf{*},\mathbf{cycle}\mathbf{)}}$ is the optimal value of the current cycle of subsystem i, and ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{U}}$ and ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{L}}$ are the upper and lower bounds of the subsystem variable x, respectively. In the first cycle, ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{start}}$ is the initial point of design. In the next cycle, ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{start}}$ will be replaced by the optimal solution ${\mathbf{x}}_{\mathrm{sys}}^{\mathbf{(}\mathbf{*},\mathbf{cycle}\mathbf{)}}$ of the system.

• Third, define and calculate system space reduction rate.

To ensure the accuracy of the solution, we should ensure that the search region contains the optimal solution. System space reduction rate is defined and calculated by

$${\mathrm{\Psi }}_{\mathbf{sys}}^{\mathbf{cycle}}=max({\mathrm{\Psi }}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}^{\mathbf{cycle}},\cdots ,{\mathrm{\Psi }}_{\mathbf{sub}\mathbf{(}\mathbf{i}\mathbf{)}}^{\mathbf{cycle}},\cdots ,{\mathrm{\Psi }}_{\mathbf{sub}\mathbf{(}\mathbf{N}\mathbf{)}}^{\mathbf{cycle}})$$

(13)

• Finally, get new reduced design space.

The new search space in system level optimization is

$${\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle}+1)}={\mathbf{x}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}-({\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{U},\mathbf{cycle})}-{\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle})}){\mathbf{\psi }}_{\mathbf{sys}}^{\mathbf{cycle}}$$

(14)

$${\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{U},\mathbf{cycle}+1)}={\mathbf{x}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}+({\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{U},\mathbf{cycle})}-{\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle})}){\mathbf{\psi }}_{\mathbf{sys}}^{\mathbf{cycle}}$$

(15)

The developed method

In this section, the developed approach is introduced thoroughly. First, after initial setting, a robust design is obtained for each subsystem. Second, the set strategy is utilized to search the robust space and solve the coupling variables during the optimization. Finally, the detailed steps and the head flowchart of the presented method are given.

MRDO under uncertainty

A complex mechanical system consists of different subsystems. Each subsystem has a specific function. Coupled relationships exist between different subsystems, which makes it difficult to evaluate the effect of uncertainties. The uncertainties of one subsystem not only affect its own subsystem but also affect other subsystems due to the relationships between these subsystems. This leads to the difficulty of solving multidisciplinary problems under uncertainty. This article develops a MRDO method by using a set-based strategy. The specific steps (take a two-subsystems as an example) are detailed below.

The following optimization framework is established for subsystem 1 considering variation of objective function, constraint condition, and the coupling state variable

$$\begin{array}{c}min{f}_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{12},{\mathbf{y}}_{21}^{\mathbf{c}})+\mathrm{\Delta }{f}_{\mathrm{sub}(1)}\hfill \\ s.t.g_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{12},{\mathbf{y}}_{21}^{\mathbf{c}})+\mathrm{\Delta }{g}_{\mathrm{sub}(1)}⩽0\hfill \\ \mathrm{\Delta }{\mathbf{y}}_{12}⩽\mathrm{\Delta }{{\mathbf{y}}_{12}}^{\mathbf{B}}\hfill \\ \mathrm{\Delta }{f}_{\mathrm{sub}(1)}⩽\mathrm{\Delta }{f}_{\mathrm{sub}(1)}^B\hfill \\ {\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}\in [{\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}^{\mathbf{L}},{\mathbf{x}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}^{\mathbf{U}}]\hfill \end{array}$$

(16.1)

where

$$\begin{array}{c}\mathrm{\Delta }{f}_{\mathrm{sub}(1)}=max\hfill \\ \multicolumn{1}{c}{\left\{|f_{\mathrm{sub}(1)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}})-f_{\mathrm{sub}(1)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}})|\right\}}\\ \mathrm{\Delta }{g}_{\mathrm{sub}(1)}=max\hfill \\ \multicolumn{1}{c}{\{g_{\mathrm{sub}(1)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}})-g_{\mathrm{sub}(1)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}})\}}\\ \mathbf{\Delta }{\mathbf{y}}_{12}=\Vert {\mathbf{y}}_{12}\mathbf{-}{\tilde{\mathbf{y}}}_{12}\Vert ;{\tilde{\mathbf{y}}}_{12}\mathbf{=}{Y}_{\mathrm{sub}(1)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\mathbf{y}}_{21}^{\mathbf{c}})\hfill \\ {\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}\in [{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}-\mathrm{\Delta }{\mathbf{p}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}+\mathrm{\Delta }{\mathbf{p}}_{\mathbf{sub}\mathbf{(}1\mathbf{)}}]\hfill \end{array}$$

(16.2)

where $f_{\mathrm{sub}(1)}()$ represents the objective function of subsystem 1. $\mathrm{\Delta }{y}_{\mathrm{ij}}^B$ is the bound value given by designers. In subsystem 1, the design variable vector is x_sub(1) and y₁₂, while $y_{21}^c$ is the constant value. ${\tilde{\mathbf{p}}}_{\mathbf{sub}\left(1\right)}$ is the uncertainty vector with the variation range Δp_sub(1) around its nominal value $\overline{\mathbf{p}}$ . $\tilde{\mathbf{y}}$ is the uncertainty vector result from $\tilde{\mathbf{p}}$ . Any proper solution strategy can be employed to solve equation (16). In general, gradient-based methods can only guarantee obtaining a local optimal design solution unless the functions used in the problem have special properties such as convexity over convex feasible regions. The non-gradient algorithm has a good global convergence effect. Note that the below problem is an unconstrained maximizing optimization problem. Considering the trade-off between computational efficiency and accuracy, the sequential quadratic programming (SQP) method will be employed to address the upper level optimization and the GA is used to solve the lower level optimization in equation (16) in this article.

A similar structure for subsystem 2 is established in equation (17) and the same optimization strategy can be adopted

$$\begin{array}{c}min{f}_{\mathrm{sub}(2)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{21},{\mathbf{y}}_{12}^{\mathbf{c}})+\mathrm{\Delta }{f}_{\mathrm{sub}(2)}\hfill \\ s.t.g_{\mathrm{sub}(1)}({\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{21},{\mathbf{y}}_{12}^{\mathbf{c}})+\mathrm{\Delta }{g}_{\mathrm{sub}(2)}⩽0\hfill \\ \mathbf{\Delta }{\mathbf{y}}_{21}⩽\mathbf{\Delta }{{\mathbf{y}}_{21}}^{\mathbf{B}}\hfill \\ \mathrm{\Delta }{f}_{\mathrm{sub}(2)}⩽\mathrm{\Delta }{f}_{\mathrm{sub}(2)}^B\hfill \\ {\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}\in [{\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}^{\mathbf{L}},{\mathbf{x}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}^{\mathbf{U}}]\hfill \end{array}$$

(17.1)

where

$$\begin{array}{c}\mathrm{\Delta }{f}_{\mathrm{sub}(2)}=max\hfill \\ \multicolumn{1}{c}{\left\{|f_{\mathrm{sub}(2)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}})-f_{\mathrm{sub}(2)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}})|\right\}}\\ \mathrm{\Delta }{g}_{\mathrm{sub}(2)}=max\hfill \\ \multicolumn{1}{c}{\{g_{\mathrm{sub}(2)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}})-g_{\mathrm{sub}(2)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}})\}}\\ \mathbf{\Delta }{\mathbf{y}}_{21}=\Vert {\mathbf{y}}_{21}\mathbf{-}{\tilde{\mathbf{y}}}_{21}\Vert ;{\tilde{\mathbf{y}}}_{21}\mathbf{=}{Y}_{\mathrm{sub}(2)}({\overline{\mathbf{x}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\mathbf{y}}_{12}^{\mathbf{c}})\hfill \\ {\tilde{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}\in [{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}-\mathrm{\Delta }{\mathbf{p}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}},{\overline{\mathbf{p}}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}+\mathrm{\Delta }{\mathbf{p}}_{\mathbf{sub}\mathbf{(}2\mathbf{)}}]\hfill \end{array}$$

(17.2)

The system-level optimization is to coordinate the coupling variables from the full autonomy subsystems optimization. The system stage optimization model is given in equation (18). Furthermore, a two-level optimization model framework is developed, and the MRDO model can be established in Figure 4

$$\begin{array}{c}min{f}_{\mathrm{sys}}({\mathbf{x}}_{\mathbf{sys}},{\overline{\mathbf{p}}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}})\hfill \\ s.t.g_j({\mathbf{x}}_{\mathbf{sys}},{\overline{\mathbf{p}}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}})+\mathrm{\Delta }{g}_j⩽0,j=1,\cdots ,m\hfill \\ \mathrm{\Delta }{f}_{\mathrm{sub}(i)}⩽\mathrm{\Delta }{f}_{\mathrm{sub}(i)}^B\hfill \\ \mathrm{\Delta }{h}_k⩽\mathrm{\Delta }{h}_k^B,k=1,\cdots ,l\hfill \\ \Vert {\mathbf{y}}_{12}-{\mathbf{y}}_{12}^{\mathbf{*}}\Vert ⩽\epsilon \hfill \\ \Vert {\mathbf{y}}_{21}-{\mathbf{y}}_{21}^{\mathbf{*}}\Vert ⩽\epsilon \hfill \\ \mathbf{x}\in [{\mathbf{x}}^{\mathbf{L}},{\mathbf{x}}^{\mathbf{U}}]\hfill \end{array}$$

(18)

where f_sys() represents the objective function of system. Several strategies can be used to handle equation (18), for instance, MDF, IDF, SAND, ATC, and so on. In this article, we use the SAND^29,30 strategy to address multidisciplinary coupling problem. The main advantage of SAND method is that it does not require a lot of system analysis and iterative process, which reduces the computational cost. The subsystem does not carry out systems analysis (SA) and subsystem optimization, and only calculates the maximum variation and value of the current point, then returns to system objective function.

Figure 4.

The flowchart of MRDO.

Set strategy MDO procedure

In order to coordinate the coupling between subsystems, the improvement of system-level optimization model is as follows

$$\begin{array}{c}\min f_{\mathrm{sys}}({\mathbf{x}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}})=f_N+f_{\mathrm{\Psi }}\hfill \\ f_N=\sum _i^n{e}^{(f_{\mathrm{sub}(i)}^{\mathrm{Nom}}({\mathbf{x}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}}))},f_{\mathrm{\Psi }}=e^{\Pi {\mathbf{\psi }}_{\mathbf{sys}}}\hfill \\ s.t.g_j({\mathbf{x}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}})⩽0,j=1,\cdots ,m\hfill \\ h_k({\mathbf{x}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sys}})=0,k=1,\cdots ,l\hfill \\ {\mathbf{y}}_{\mathbf{sys}}=Y_{\mathrm{sys}}({\mathbf{x}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}},{\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{c}})\hfill \\ {\mathbf{x}}^{\mathbf{L}}⩽\mathbf{x}⩽{\mathbf{x}}^{\mathbf{U}}\hfill \end{array}$$

(19)

where the f_N represents the effect of the subsystem objective function, which is designed to make each subsystem relatively small. The closer the system is to the solution of the subsystem, the closer f_N is to zero. Therefore, to minimize f_N is to minimize the objective function of the subsystem. The f_ψ is to narrow the size of the design space as much as possible.

The steps and flowchart of the proposed method

We now summarize the implementation procedure of set-based MRDO. The stopping criteria of set-based MRDO method are as follows. (1) The objective approaches stability: the convergence error of the objective function between two consecutive cycles is less than a given value ε. (2) The cycle times exceed the maximum cycle, which means that the solution fails.

Main steps of proposed method are given as follows:

• Step 1: Give initial points and bounds: ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{start}}$ , ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{L}}$ , ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{U}}$ , y_sub(i), p_sub(i), ${\mathbf{x}}_{\mathbf{sys}}^{\mathbf{start}}$ , y_sys, p_sys, Δh ^B , Δy ^B , and Δf ^B .

• Step 2: Set $f_{\mathrm{sys}}^{\mathrm{last}}$  = 1, $f_{\mathrm{sub}(i)}^{\mathrm{last}}$  = 1, convergence error = 1, cycle = 1, and let the maximum cycle, that is, Maxcycle, be 50.

• Step 3: Calculate all subsystem maximum $f_{\mathrm{sub}(i)}^{max}$ and minimum $f_{\mathrm{sub}(i)}^{min}$ .

• Step 4: Full autonomy is given to each subsystem to obtain robust optimal solution, ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{(*,\mathbf{cycle})}$ and ${\mathbf{y}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{(*,\mathbf{cycle})}$ , by handling the optimization model defined in equations (16) and (17). (This step can be seen in the dotted red box of Figure 5.)

• Step 5: According to ${\mathbf{x}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}$ , ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{L}}$ and ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{U}}$ , and the initial point of current cycle, calculate the ${\mathbf{\psi }}_{\mathbf{sub}(i)}^{(*,\mathbf{cycle})}$ by using equation (12).

• Step 6: In order to search the optimal solution as much as possible, we need to obtain the maximum ${\mathbf{\psi }}_{\mathbf{sub}(i)}^{(\mathbf{cycle})}$ (equation (13)).

• Step 7: Let the initial point of system stage optimization be ${\mathbf{x}}_{\mathbf{sys}}^{\mathbf{start}}=({\mathbf{x}}_{\mathbf{sys}}^{\mathbf{start}},{\mathbf{\psi }}_{\mathbf{sys}}^{\mathbf{cycle}},{\mathbf{y}}_{\mathbf{sys}},{\mathbf{p}}_{\mathbf{sys}})$ . Construct system optimization model and get system optimal solution ${\mathbf{x}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}$ and ${\mathbf{y}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}$ (equations (18–19)).

• Step 8: According to ${\mathbf{x}}_{\mathbf{sys}}^{(*,\mathbf{cycle})}$ , update ${\mathbf{x}}_{\mathbf{sub}\mathbf{(}i\mathbf{)}}^{\mathbf{start}}$ . Narrow the search region by adopting equations (14) and (15). Furthermore, it is necessary to ensure that the new interval bound does not exceed the initial bound, do (if ${\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle}+1)}⩽{\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle}},{\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle}+1)}=$ ${\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{L},\mathbf{cycle})}$ , similar to ${\mathbf{x}}_{\mathbf{sys}}^{(\mathbf{U},\mathbf{cycle}+1)}$ ).

• Step 9: Let cycle = cycle + 1, implement next iteration.

• Step 10: If the cycle < Maxcycle, and the convergence error ≥ ε (convergence error = abs $(f_{\mathrm{sys}}^{(*,\mathrm{cycle})}-f_{\mathrm{sys}}^{(*,\mathrm{cycle}-1)}$ ), in this article, we set the ε = 0.0001), then go back to step 4.

• Step 11: Otherwise, terminate the cycle process, output current optimal value.

Figure 5.

(a) The architecture of set-based MRDO and (b) the flowchart of set-based MRDO.

Based on the above discussions, for a two-subsystems MRDO problem, the architecture of set-based MRDO is shown in Figure 5(a) and a simplified flowchart of set-based MRDO is given in Figure 5(b).

Tests and results

In this section, two numerical examples and a speed reducer design problem are adopted to verify the accuracy and applicability of the developed approach. In addition, an air cooling battery thermal management system (BTMS) design problem is successfully solved by the proposed method. This means that the proposed method is the ability for addressing the practical engineering problem.

Test 1: Mathematical example 1

First, a simple math example was tested. The structure of problem is given in Figure 6. The simple problem was defined with only two design variables so that the results of the proposed method can easily be visualized in a two-dimensional plot. The two design variables are the x₁ and the x₂. The initial interval of x₁ is [–5, 5] and x₂ is [0, 6]. There are two subsystems: f_sub(1) and f_sub(2). Details for the mathematical model are provided as follows

$$\begin{array}{c}min{f}_1={(x_1^2+x_2-11)}^2\hfill \\ min{f}_2={(x_2^2+x_1-7)}^2\hfill \\ s.t.g_1=\frac{x_1^2{x}_2}{20-1}⩽0\hfill \\ g_2=\frac{x_1^2}{10-x_2+4}⩽0\hfill \end{array}$$

(20)

Figure 6.

The structure of math example 1.

In this problem, the initial design point is [0,0]. We assume an uncontrollable variation in design variables x₁ and x₂. The two cases are calculated separately. The variations of Δx₁ and Δx₂ are given to be continuous and within ±0.05 from the nominal value in Case 2, and the Case 3 variation of Δx₁ and Δx₂ is ±0.1. Case 1 is a nominal design without uncertainty. The variation bounds in Δf₁ and Δf₂ are given to be 12 and 120, respectively. The set strategy–based MRDO method was applied to the example problem, and the results are shown in Table 1. The results were plotted in Figures 7 and 8 which include the three cases of individual subsystem optima. We can see that the individual subsystem optima are in different boundaries of the design region. With the increase of uncertainty, the robust solution is becoming increasingly deviated from the original boundary. The yellow region is the feasible region. The new boundaries for the design region are shown as purple rectangle in the upper left corner (Case 2: x₁ = [–3.4839, –0.6107], x₂ = [2.9510, 6], Case 3: x₁ = [–3.3469, –0.603], x₂ = [3.006, 6]). For Case 3, the calculation times of subsystem 1 and subsystem 2 are 578 and 5206, correspondingly. Figure 9 shows the convergence error of the system-level cycle. It can also be seen that the accuracy and efficiency of the proposed method take only three cycles.

Table 1.

Optimization results.

		x 1	x 2	f 1 *	f 2 *
Case 1	Sub 1	−1.8257	6	2.7778	–
	Sub 2	−0.6191	4.0383	–	75.4987
Case 2	Sub 1	−1.7682	6	4.401	–
	Sub 2	−0.6107	4.0937	–	91.8203
Case 3	Sub 1	−1.7107	6	6.2780	–
	Sub 2	−0.6025	4.1493	–	111.3963

Figure 7.

The optimal solution of test 1 in individual subsystem.

Figure 8.

The optimal solution of the individual subsystem 2.

Figure 9.

The convergence error of optimization model.

In order to demonstrate the accuracy of the proposed method, a multilevel coordinate search (MCS) method using uniform distribution disordering the variation over the tolerance region was employed. Two cases’ robust optimal designs were used to verify the verification for the constraint robustness. The comparisons of constraint function robustness are shown in Figure 10 (1000 points are shown). In Figure 10, when there are uncertainties, the constraint conditions that the robust design is feasible and the nominal design scheme is hard to meet the design requirements.

Figure 10.

The comparisons of constraint function robustness for g₁ and g₂.

We defined the reduction ratio of design region ψ as follows

$$\psi =\left(\frac{1-\left(O_{region}-N_{region}\right)}{O_{region}}\right)\times 100\%$$

(21)

For Case 3, the ψ = (1 – 8.215 / 60) × 100% = 86.31%. The design space has been reduced by more than 86%, which provides more valuable information for the designer’s subsequent analysis and decision. In addition, if it is smaller than new design region, then the optimal solution of the subsystem will not be included, and then the key information will be lost.

Test 2: Mathematical example 2

A classic MDO problem including two disciplines was used for verification. ²⁸ Each discipline has a coupled state variable; nonlinear coupled relationship exists between two subsystems (Figure 11). The optimization model is expressed as follows

$$\begin{array}{c}minf=x_1^2+x_2+y_1+e^{-y_2}\hfill \\ s.t.g_1=1-\frac{y_1}{3.16}⩽0\hfill \\ g_2=\frac{y_2}{24}-1⩽0\hfill \\ h_1=x_1+x_2+x_3^2-0.2y_2-y_1\hfill \\ h_2=\sqrt{y_1}+x_2+x_3-y_2\hfill \\ 0⩽x_1,x_2⩽10,-10⩽x_3⩽10\hfill \end{array}$$

(22)

Figure 11.

The MDO model of math problem 2.

To illustrate this question more intuitively, we divide the problem into two subsystems, as shown in the Figure 11. The variations of Δx = [Δx₁, Δx₂, Δx₃] is given to be continuous and within ±0.2 from the nominal value in Case 2, while Case 1 is without uncertainty. In this example, the deterministic solution is x* = [0.0001, 0, 1.9776, 3.16, 3.7553]. The set strategy–based MRDO method was employed to address the example problem. We can obtain a system solution x* = [0.0007, 0.0000, 2.0328, 3.3600, 3.8659]. (Note that this solution is not the robust solution.) For nominal design, the optimal solution for subsystem 1 is x_sub(1) = [0.0007, 3.9104, 0] and subsystem 2 is x_sub(2) = [9.4091, 0, 9.9374] and coupled state variable y = [ $y_{12}^c$ , $y_{21}^c$ ] = [3.16, 13.1240]. For robust design, the optimal robust solution for subsystem 1 is x_sub(1) = [0.0005, 4.9840,0] and subsystem 2 is x_sub(2) = [8.3615, 0, 9.8418], and y = [ $y_{12}^c$ , $y_{21}^c$ ] = [3.36, 3.8659]. The results are shown in Table 2. The results in Table 2 are plotted in Figure 12. For Case 2, the calculation times of Sys, Sub 1, and Sub 2 are 1362, 96, and 256, correspondingly. The convergence error of the system-level cycle is given in Figure 13. It can be seen that the accuracy and efficiency of the proposed method take only three cycles. MCS results of g₁ and g₂ are presented in Figure 14. In Figure 13, when there are uncertainties, the constraint conditions that the robust design is feasible and the g₁ of nominal design scheme is hard to meet the design requirements.

Table 2.

Optimization results.

		x 1	x 2	x 3
Origin	x l	0	0	−10
	x u	10	10	10
New	x l	0.0005	0	−5.7760
	x u	8.3615	4.9840	9.8417

Figure 12.

The optimal solution of two cases in individual subsystem.

Figure 13.

The convergence error of optimization model.

Figure 14.

The comparisons of constraint function robustness for g₁ and g₂.

In Table 2, the “Origin” row describes the initial lower and upper bounds of design variables. The “New” row describes the design results using the proposed method. Due to the existence of the coupled state variables, we cannot obtain the optimal solution of subsystem in individual optimization. The new boundaries for the design space are shown as light green cube, and the subsystem optima within the new design space are shown as pink markers. In this problem, the ψ = (1 – 650.86 / 2000) × 100% = 67.46%. The design region has been reduced by more than 67%, that is, this method can greatly decrease the rework and reanalysis. It has broad prospects for comprehensive design of complex systems.

Test 3: The speed reducer design

The speed reducer design problem is a standard MDO test. ⁴⁵ Figure 15 illustrates the speed reducer diagram. The design objective is to minimize the volume of the structure (objective function). The main constrains are the bending stress and contact stress of the gear tooth, the torsional deformation, and stress requirements of the shaft

$$\begin{array}{c}minf=0.7854x_1{x}_2^2(3.333x_3^2+14.9334x_3-43.0934)\hfill \\ \multicolumn{1}{c}{+7.477x_6^3+}\\ 0.7854x_4{x}_6^2-1.5079x_1{x}_6^2+7.477x_7^3\hfill \\ \multicolumn{1}{c}{+0.7854x_5{x}_7^2-1.5079x_1{x}_7^2}\\ s.t.g_1=\frac{27}{x_1{x}_2^2{x}_3}-1⩽0\hfill \\ g_2=\frac{397.5}{x_1{x}_2^2{x}_3^2}-1⩽0\hfill \\ g_3=\frac{1.93x_4^2}{x_2{x}_3{x}_6^4}-1⩽0\hfill \\ g_4=\frac{1.93x_5^2}{x_2{x}_3{x}_7^4}-1⩽0\hfill \\ g_5=\frac{10}{x_6^3}\sqrt{{\left(\frac{745x_4}{x_2{x}_3}\right)}^2+1.69\times {10}^7}-1100⩽0\hfill \\ g_6=\frac{10}{x_7^3}\sqrt{{\left(\frac{745x_5}{x_2{x}_3}\right)}^2+1.575\times {10}^8}-850⩽0\hfill \\ g_7=x_2{x}_3-40⩽0\hfill \\ g_8=5-\frac{x_1}{x_2}⩽0\hfill \\ g_9=\frac{x_1}{x_2}-12⩽0\hfill \end{array}$$

$$\begin{array}{c}{g}_{10}=\frac{1.5x_6+1.9}{x_4}-1⩽0\hfill \\ g_{11}=\frac{1.5x_7+1.9}{x_5}-1⩽0\hfill \end{array}$$

(23)

where x₁ is gear tooth width, x₂ is the teeth module, x₃ is the number of teeth of the pinion, x₄ is the distance between the bearing 1, x₅ is the distance between bearing 2, x₆ is the diameter of the shaft 1, x₇ is the diameter of the shaft 2, g₁ and g₂ are constraints of gear tooth bending stress and contact stress, g₃∼g₈ are constraints of the shaft deformation, stress, and g₉∼g₁₁ are geometric constraints. Assume that the x₂ and x₆ have interval uncertainty (Δx₂, Δx₆) = (0.01, 0.1). We divide the reducer design model into three subsystems, namely, subsystem 1 (gear), subsystem 2 (shaft 1), and subsystem 3 (shaft 2), respectively (Figure 16). The acceptable objective variation range of $\mathrm{\Delta }{f_1}^{\mathrm{B}}$ and $\mathrm{\Delta }{f_2}^{\mathrm{B}}$ is presumed 150 and 100, respectively. In subsystem 3, there is no interval uncertainty, so $\mathrm{\Delta }{f_3}^{\mathrm{B}}$ does not need to be given. In the speed reducer design example, the deterministic solution is x* = [3.5, 0.7, 17, 7.3, 7.715, 3.35, 5.287]. The set-based MRDO method is applied to this problem. For nominal design, the optimal solution for subsystem 1 is ${\mathbf{x}}_{\mathbf{sub}(1)}^{*}$  = [3.5, 0.7, 17], subsystem 2 is ${\mathbf{x}}_{\mathbf{sub}(2)}^{*}$  = [3.6, 0.8, 28, 7.3, 3.345], and subsystem 3 is ${\mathbf{x}}_{\mathbf{sub}(3)}^{*}$  = [3.6, 0.8, 28, 7.714, 5.286]; the optimal objective value is f* = [f₁*, f₂*, f₃*] = [1581.3, 283.3, 1121.8].

Figure 15.

The structural diagram of speed reducer.

Figure 16.

The MDO model of speed reducer.

For robust design, the optimal robust solution for subsystem 1 is ${\mathbf{x}}_{\mathbf{sub}(1)}^{*}$  = [3.55, 0.7, 17], subsystem 2 is ${\mathbf{x}}_{\mathbf{sub}(2)}^{*}$  = [3.6, 0.8, 28, 7.3, 3.445], and subsystem 3 is ${\mathbf{x}}_{\mathbf{sub}(3)}^{*}$  = [3.6, 0.8, 28, 7.714, 5.286]; the optimal robust objective value is f* = [f₁*, f₂*, f₃*] = [1603.9, 309.3, 1121.8].

The results are shown in Table 3. For Case 2, the calculation times of Sys, Sub 1, Sub 2, and Sub 3 are 3025, 391, 731, and 46, correspondingly. The convergence error of the system-level cycle is shown in Figure 17. It also only takes three cycles to get the optimal which is the accuracy and efficiency of the proposed method. The “O” row gives the initial bounds of design variables. The “N” row represents the design results using the proposed method. Because there are more than three design variables, we cannot be intuitive to see the design space changes. So the “Δ” row describes the reduction ratio of each design variables (by equation (21)). It is can be seen from the “Δ” row that every design variables have been reduced by more than 87% except x₂ and x₃. The biggest of the multidisciplinary compromise is in f₂ (robust design solution is 9% more than nominal design). The system optima for f₁ and f₃ are relatively close to the subsystem optimal solution. As similar to the discussion in “Test 1: Mathematical example 1” and “Test 2: Mathematical example 2” sections, MCS method is applied. The nominal design and robust design were used to verify the confirmation for the constraint robustness. There are no uncertainties in subsystem 3, so we do not analyze the constraints of g₄, g₆, and g₁₁. The comparisons of constraint function robustness are shown in Table 4 and Figure 18 (1000 points are shown). In Figure 17, when there are uncertainties, the constraint conditions that the robust design is feasible and the nominal design scheme for g₅ and g₈ is difficult to meet the design requirements.

Table 3.

Optimization results (unit: mm).

	x 1	x 2	x 3	x 4	x 5	x 6	x 7
O	[2.6, 3.6]	[0.7, 0.8]	[17, 28]	[7.3, 8.3]	[7.3, 8.3]	[2.9, 3.9]	[5.0, 5.5]
N	[3.47, 3.6]	[0.7, 0.8]	[17, 28]	[7.30, 7.34]	[7.68, 7.75]	[3.40, 3.50]	[5.27, 5.31]
Δ	87%	0	0	96%	93%	90%	92%

Figure 17.

The convergence error of optimization model for speed reducer.

Table 4.

The probability of constraint function robustness for subsystem 2 and subsystem 3.

	g 1	g 2	g 3	g 5	g 7	g 8	g 9	g 10
Nominal	1	1	1	0.512	1	0.497	1	1
Robust	1	1	1	1	1	1	1	1

Figure 18.

The comparisons of constraint function robustness for g₅ and g₈.

Test 4: Air cooling BTMS design

An air-cooled module containing eight cells had been well defined by National Renewable Energy Laboratory in their previous work. ⁴⁶ It is discharged and charged under normal driving conditions. Air cooling BTMS design involves five design variables that are d₁, d₂, d₃, d₄, and v. Schematic diagram of the battery module is shown in Figure 19. More details can be seen in Fan et al. ⁴⁶ and Li et al. ⁴⁷

Figure 19.

Schematic diagram of the battery module.

Some notations are defined in this work. TD is the maximum temperature difference of eight cells (calculated by mean temperature) and TSD is the temperature standard deviation. Computational domain of the battery module with air cooling is shown in Figure 20.

Figure 20.

Computational domain of the battery module with air cooling.

Without loss of generality, the mathematical model of multidisciplinary optimization is given with all-in-one formulation as follows

$$\{\begin{array}{c}\mathrm{find}\mathbf{x}=[d_1,d_2,d_3,d_4,v]\hfill \\ \multicolumn{1}{c}{minf=\mathrm{TD}}\\ \multicolumn{1}{c}{s.t.h_1=\mathrm{TD}-6-g_1;}\\ \multicolumn{1}{c}{h_2=\mathrm{TSD}-1.5-g_2;}\\ \multicolumn{1}{c}{1⩽d_1⩽4,1⩽d_2⩽4}\\ \multicolumn{1}{c}{1⩽d_3⩽4,1⩽d_4⩽4}\\ \multicolumn{1}{c}{0.002⩽v⩽0.02}\end{array}$$

(24)

In this article, to reduce the cost of the expensive computational fluid dynamics (CFD) process and have high approximation accuracy, Kriging model is applied for surrogate model and built similar to Li et al. ⁴⁷ and Stein. ⁴⁸ The maximum temperature difference and temperature standard deviation surrogate models are represented by equations (25 and 26)

$$\mathrm{TD}={\tilde{f}}_{\mathrm{TD}}(d_1,d_2,d_3,d_4,v)$$

(25)

$$\mathrm{TSD}={\tilde{f}}_{\mathrm{TSD}}(d_1,d_2,d_3,d_4,v)$$

(26)

In this work, the air cooling battery module is analyzed in ANSYS by incorporating the basic information as mentioned in Zhang et al. ⁴⁹ In this example, large weight is assigned to the battery temperature difference. The design variables x with Δx = [±0.05, ±0.05, ±0.05, ±0.05, ±0.001] are considered to be uncertainty variables. The results are shown in Table 5. Temperature distributions (unit: K) within the cells for the initial scheme of battery module, nominal optimum of battery module, and robust optimum of battery module are shown in Figure 21(a)–(c), correspondingly. Compared to the initial scheme, the robust scheme has better cooling performance, indicating that the developed method is effective.

Table 5.

Optimization results of air cooling BTMS.

	Notation	Range	Initial	DO	RO
Design variables	d 1 (mm)	[1.00, 4.00]	4.0000	1.6459	3.3500
	d 2 (mm)	[1.00, 4.00]	4.0000	1.7491	1.1150
	d 3 (mm)	[1.00, 4.00]	4.0000	1.8622	2.6245
	d 4 (mm)	[1.00, 4.00]	4.0000	2.9487	3.7854
	v (kg/s)	[0. 002, 0.02]	0.0120	0.015367	0.018709
Objective	TD (K)	–	7.2399	4.1787	5.2654
Constraint	TSD (K)	–	2.7553	0.9543	1.0465

BTMS: battery thermal management system.

Figure 21.

Temperature distributions (unit: K) within the cells for the (a) initial scheme of battery module, (b) nominal optimum of battery module, and (c) robust optimum of battery module.

MCS method is applied to verify the accuracy of the developed method. The comparisons of constraint function robustness for g₁ and g₂ are presented in Figure 22. Only 50 points are shown, since the simulation is time-consuming. Note that this operation does not change the nature of the problem. In Figure 22, when there are no uncertainties, the nominal design scheme and the robust design scheme for g₁ are both feasible. However, in Figure 22, when there are uncertainties, the constraint conditions that the robust design is still feasible and the nominal design scheme for g₂ is hard to meet the design requirements.

Figure 22.

The comparisons of constraint function robustness for g₁ and g₂.

Conclusion

Coupled relationships exist between different subsystems, which makes it difficult to evaluate the effect of uncertainties. The uncertainties of one subsystem not only affect its own subsystem but also affect other subsystems due to the relationships between these subsystems. This leads to the difficulty of solving multidisciplinary problems under uncertainty. In this article, first, two design methodologies based on point design and set design are compared. Then, a two-level optimization model framework is developed; the lower level is to obtain the optimal robust solution for each subsystem by adopting full autonomy optimization and the system-level optimization is to coordinate the coupling variables from the full autonomy subsystem optimization. Furthermore, we make a try to integrate the set strategy into the MRDO method. Finally, two mathematics examples and the speed reducer design problem verify the validity and accuracy of the proposed method. A practical engineering problem, namely air cooling BTMS design problem, is successfully solved by the proposed method. This means that the proposed method is the ability for addressing the practical engineering problem.

Moreover, the main advantages of the developed method can be summarized as follows. (1) It can effectively coordinate all subsystems and get the global optimal. (2) It obviously narrows the design space and provides feasible direction for the designer’s subsequent design and decision. (3) It is relatively simple and easy to be implemented and the solving process is accurate and stable.

In addition, although a set-based design found many merits for complex systems, it has a lot of work to be done, for instance, if there are time-varying uncertainties in the system. As part of future work, the developed method with surrogate model should be developed.