# Robust Topology Optimization Based on Stochastic Collocation Methods under Loading Uncertainties

Authors: Qinghai Zhao; Xiaokai Chen; Zheng-Dong Ma; Yi Lin
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A robust topology optimization (RTO) approach with consideration of loading uncertainties is developed in this paper. The stochastic collocation method combined with full tensor product grid and Smolyak sparse grid transforms the robust formulation into a weighted multiple loading deterministic problem at the collocation points. The proposed approach is amenable to implementation in existing commercial topology optimization software package and thus feasible to practical engineering problems. Numerical examples of two- and three-dimensional topology optimization problems are provided to demonstrate the proposed RTO approach and its applications. The optimal topologies obtained from deterministic and robust topology optimization designs under tensor product grid and sparse grid with different levels are compared with one another to investigate the pros and cons of optimization algorithm on final topologies, and an extensive Monte Carlo simulation is also performed to verify the proposed approach.


---

## Body

## 1. Introduction

Structural optimization is an intrinsic element of engineering design for the purpose of improving the structural performance and reducing costs. The major ingredients of structural optimization are, in general, divided into three levels, namely, sizing, shape, and topology optimization. The aim of classic topology optimization is to obtain an optimal material or design parameter distribution at a given structural domain under nominal (i.e., deterministic) material properties, geometry, and loading conditions. Until now, a range of topology optimization strategies has emerged, including the homogenization method, introduced by Bendsøe and Kikuchi [1], the density-based method [2], bilinear evolutionary structural optimization (BESO) [3], and level set method [4]. Details of various representative schemes can be found in comprehensive reviews and text books [5–7].Traditionally, structural topology optimizations are conducted in a deterministic manner, known as deterministic topology optimization (DTO), where the design is determined without taking into account various sources of uncertainties [6]. However, uncertainties are unavoidably observed in real-world applications due to insufficient knowledge, manufacturing errors, changeable environment, and so forth. This may lead to the vulnerable optimum structure or infeasible topologies due to the fluctuation of the structure performance. Therefore, there is a strongly increasing requirement to take the effect of uncertainty into consideration for optimal topologies in structural design.Recent years have witnessed various statements accounting for uncertainty in topology optimization process. In reliability-based topology optimization (RBTO), optimal designs with acceptable failure probabilities have been achieved, where uncertainty is manifested in the probability constraint functions. Many works dealing with applications of RBTO in structural design have appeared over the last decades. One of pioneering works can be traced back to Kharmanda et al. [8], where reliability analysis is decoupled from optimization procedures and carried out at the beginning of the optimization loops, followed by equivalent DTO. Different optimization strategies have been increasingly applied for reliable designs in various mechanic and multiphysical fields. To cite a few, optimization algorithms via the density-based method have been demonstrated for microelectromechanical systems (MEMS) [9], geometrically nonlinear structures [10], and an automotive control arm [11]. Optimization problems via the BESO have been reported for a vehicle’s hood reinforcement [12] and electrothermal-compliant mechanisms [13]. Topology optimization of compliant mechanisms with loads, material properties, and member geometries uncertainties is demonstrated for the level set method [14].Distinct from the abovementioned formulations considering the uncertainties in structural topology optimization problems, robust topology optimization (RTO) has recently emerged to optimize the objective performance while simultaneously minimizing its sensitivity with respect to uncertainties. When the uncertainty is assigned unknown-but-bounded convex models, the worst case approach is followed in general. De Gournay et al. [15] propose the so-called worst case or robust optimal design problem for minimal compliance based on the level set method. Amir et al. [16] demonstrate a fair comparison of a worst case formulation and a stochastic formulation in robust topology optimization procedure. Guo et al. [17] present structural topology optimization considering the uncertainty of boundary variations via level set approach, where the objective function is modeled by the compliance and fundamental frequency of structure enduring the worst case perturbation.In the context of uncertainty parameters defining probability distributions, the probabilistic description of the performance function is usually characterized by its statistical quantities such as mean value and standard deviation. In this scenario, a fundamental issue associated with these RTO procedures is to calculate the statistical moments accurately and efficiently. The Monte Carlo simulations (MCS) are the most widely used sampling-based methods in robust design due to their accuracy and easy implementation. However, it is time-consuming. Schevenels et al. [18] present a robust topology optimization approach for the design of macro-, micro-, or nanostructures, accounting for spatially varying manufacturing errors, where the statistical moments and their sensitivities are estimated by means of the MCS. Jansen et al. [19] propose a robust approach to topology optimization taking into account geometric imperfections, where a sampling method (i.e., 100 Monte Carlo samples) is used to estimate these statistics that are defined as a weighted sum objective function. A new efficient and accurate approach to robust structural topology optimization under loading uncertainty is developed by Zhao and Wang [20]. The Monte Carlo method and matrix decomposition are employed for concentrated loads, and the Karhunen-Loeve expansion and its orthogonal properties of random variables are used for distributed loads. Alternatively, the stochastic collocation method has been gaining more attention in robust design. The idea is to calculate the statistical moments by direct numerical integration. Chen et al. [21] employ the univariate dimension-reduction (UDR) method combined with Gauss type quadrature sampling for calculating statistical moments with consideration of random field uncertainty in loading and material properties. Lazarov et al. [22] introduce the stochastic collocation methods in topology optimization for mechanical systems with material and geometric uncertainties. The robust designs are obtained by utilizing well-developed deterministic solvers.Efficient and accurate generation of statistic moments is a central problem for robust design, especially for structural topology optimization with infinite-dimensional property. In this paper, we present a methodology for robust topology optimization problems with the stochastic collocation method to integrate the statistic moments considering loading uncertainty. Then, the robust topology optimization formulation is transformed into a weighted multiple loading deterministic problem at the collocation points. For the construction of the stochastic collocation method, the collocation points with corresponding weights are generated by full tensor product grid and Smolyak sparse grid, respectively. Integrations based on full tensor product suffers from thecurse of dimensionality, while the sparse grid method originating from the Smolyak algorithm will dramatically reduce the number of collocation points without losing much accuracy. In addition, existing commercial software packages are capable of application in robust topology optimization in these settings.The paper is organized as follows: A brief description of the deterministic and robust topology optimization formulation is presented in Section2. In Section 3, the stochastic collocation methods including tensor product grid and sparse grid are demonstrated. The computational results for deterministic and robust topology optimization together with several numerical examples are illustrated in Section 4. Conclusion and future works are discussed in Section 5.


## 2. Topology Optimization Formulation

### 2.1. Deterministic Topology Optimization Formulation

The generalized topology optimization problem as shown in Figure1 is that the predetermined design domain Ω, with respect to the given loads and boundary conditions, is represented by the material distribution function so as to optimize given objective function f (e.g., compliance) under limited volume constraint.Figure 1
Generalized topology optimization problem.The structural topology optimization considering the general elasticity problem is reformulated in the uniform context of the principle of virtual displacement. The optimization problem for deterministic design can be written in a weak form as(1)
min
⁡
l
u
=
∫
Ω
ε
T
u
D
ρ
ε
v
d
Ω
s.t.
a
u
,
v
=
l
v
,
∀
v
∈
K
∫
Ω
ρ
d
Ω
≤
V
∗
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where K is the space of kinematically admissible displacement fields. u is the equilibrium displacement field, and v is the arbitrary virtual displacement in the space K. V
∗ is the prescribed volume constraint. The energy bilinear form and the load linear form are defined as follows:(2)
a
u
,
v
=
∫
Ω
ε
T
u
D
ρ
ε
v
d
Ω
,
l
v
=
∫
Ω
F
·
v
d
Ω
+
∫
Γ
t
t
·
v
d
Γ
,where ε is the strain tensor. F and t are the body forces and the boundary tractions on the boundary Γ
t, respectively. Considering linear isotropic materials, Young’s modulus is a function of design variables ρ given by the density-based method [6] as follows: (3)
E
=
E
ρ
=
ρ
p
E
0
,where E
0 is Young’s modulus of solid material and p is the penalization power to ensure solid-void solutions.By using FEM discretization for the linear elasticity, the optimization problem can be rewritten as follows:(4)
min
⁡
f
ρ
=
F
T
U
=
U
T
K
U
=
∑
e
=
1
N
e
E
e
x
e
u
e
T
k
0
u
e
s.t.
V
=
∑
e
∈
N
ρ
e
v
e
≤
V
∗
K
ρ
U
=
F
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where f
(
ρ
) is the compliance of the continuum structure, U is the displacement vector, and N is a set of finite elements. u
e and k
0 are the element displacement and the element stiffness matrix with unit Young’s modulus, respectively, v
e is the element volume, and the density of each element is represented by one density variable ρ
e. Therefore, the standard FEM global stiffness matrix K is expressed as(5)
K
ρ
=
∑
e
∈
N
K
e
ρ
e
=
∑
e
∈
N
∫
Ω
e
B
T
D
ρ
e
B
d
Ω
,where K
e
(
ρ
e
) is the element stiffness matrix, B is the strain-displacement matrix of shape function derivatives, Ω
e is element domain, and D
(
ρ
e
) is the element constitutive matrix, defined in 2D as(6)
D
ρ
e
=
E
ρ
e
1
-
ν
2
1
ν
0
ν
1
0
0
0
1
-
ν
2
,where ν denotes Poisson’s ratio. In general, the optimal solution obtained by direct implementation of the method is prone to problems with checkerboards and mesh dependency. A large number of restriction schemes have been proposed to ensure manufacturability and mesh independency, such as perimeter control method, mesh-independent filtering methods, and Heaviside projection. A comprehensive comparison of the existing restriction schemes is provided in the paper of Sigmund [23].


### 2.2. Robust Topology Optimization Formulation

The topology optimization of linear elastic structures subject to loading uncertainties is considered in this paper. Due to the variation of uncertain parametersξ, the performance function f
(
ρ
,
ξ
) has probability distribution description. Thus, the robustness of the design objective can be achieved by optimizing the mean performance and minimizing the standard deviation of the performance simultaneously, as shown in Figure 2. In Figure 2, point A denotes the deterministic optimum not considering the uncertainty, while point B denotes the robust optimum. The performance of deterministic optimum is better than that of the robust optimum. However, its variation is wider than the robust optimum.Figure 2
The concept of robust design model.Therefore, a main challenge in performing RTO for realistic problems is the implementation of an algorithm, which inherently represents a multiobjective optimization problem. One efficient way is to combine all objective functions into a single weighted function. In mathematical terms, a general robust topology optimization task can be stated as(7)
find:
ρ
min:
μ
f
ρ
,
ξ
+
κ
σ
f
ρ
,
ξ
s.t.:
V
ρ
V
∗
-
1
≤
0
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where κ is the weight to be determined by the designer that satisfies κ
⩾
0. Note that the weights have to be set a priori, which make a trade-off decision between an optimal mean performance and minimization of standard deviation of the performance to obtain a compromise solution. μ
f and σ
f are the mean value and the standard deviation of the performance function, respectively. V
∗ is the material resource constraint. The design variables ρ are bounded by their lower and upper bounds ρ
m
i
n and ρ
m
a
x, respectively. In a probabilistic setting, the mean value μ
f and the standard deviation σ
f can be expressed as(8)
μ
f
ρ
,
ξ
=
E
f
ρ
,
ξ
=
∫
Ω
f
ξ
p
ξ
d
ξ
,
σ
f
ρ
,
ξ
=
E
f
ρ
,
ξ
-
μ
f
ρ
,
ξ
2
=
∫
Ω
f
ξ
-
μ
f
2
p
ξ
d
ξ
1
/
2
,where E is the expectation operator and p
(
ξ
) is the joint probability density function (PDF) of the random variables. In practice, it is very difficult and even impossible to calculate the mean value and the variation of the performance function through the multidimensional integration in (8). Such difficulties have motivated the development of various numerical methods such as simulation-based methods, local expansion-based methods, and the stochastic collocation methods. The reader is referred to Lee and Chen [24] for a more detailed review.


## 2.1. Deterministic Topology Optimization Formulation

The generalized topology optimization problem as shown in Figure1 is that the predetermined design domain Ω, with respect to the given loads and boundary conditions, is represented by the material distribution function so as to optimize given objective function f (e.g., compliance) under limited volume constraint.Figure 1
Generalized topology optimization problem.The structural topology optimization considering the general elasticity problem is reformulated in the uniform context of the principle of virtual displacement. The optimization problem for deterministic design can be written in a weak form as(1)
min
⁡
l
u
=
∫
Ω
ε
T
u
D
ρ
ε
v
d
Ω
s.t.
a
u
,
v
=
l
v
,
∀
v
∈
K
∫
Ω
ρ
d
Ω
≤
V
∗
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where K is the space of kinematically admissible displacement fields. u is the equilibrium displacement field, and v is the arbitrary virtual displacement in the space K. V
∗ is the prescribed volume constraint. The energy bilinear form and the load linear form are defined as follows:(2)
a
u
,
v
=
∫
Ω
ε
T
u
D
ρ
ε
v
d
Ω
,
l
v
=
∫
Ω
F
·
v
d
Ω
+
∫
Γ
t
t
·
v
d
Γ
,where ε is the strain tensor. F and t are the body forces and the boundary tractions on the boundary Γ
t, respectively. Considering linear isotropic materials, Young’s modulus is a function of design variables ρ given by the density-based method [6] as follows: (3)
E
=
E
ρ
=
ρ
p
E
0
,where E
0 is Young’s modulus of solid material and p is the penalization power to ensure solid-void solutions.By using FEM discretization for the linear elasticity, the optimization problem can be rewritten as follows:(4)
min
⁡
f
ρ
=
F
T
U
=
U
T
K
U
=
∑
e
=
1
N
e
E
e
x
e
u
e
T
k
0
u
e
s.t.
V
=
∑
e
∈
N
ρ
e
v
e
≤
V
∗
K
ρ
U
=
F
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where f
(
ρ
) is the compliance of the continuum structure, U is the displacement vector, and N is a set of finite elements. u
e and k
0 are the element displacement and the element stiffness matrix with unit Young’s modulus, respectively, v
e is the element volume, and the density of each element is represented by one density variable ρ
e. Therefore, the standard FEM global stiffness matrix K is expressed as(5)
K
ρ
=
∑
e
∈
N
K
e
ρ
e
=
∑
e
∈
N
∫
Ω
e
B
T
D
ρ
e
B
d
Ω
,where K
e
(
ρ
e
) is the element stiffness matrix, B is the strain-displacement matrix of shape function derivatives, Ω
e is element domain, and D
(
ρ
e
) is the element constitutive matrix, defined in 2D as(6)
D
ρ
e
=
E
ρ
e
1
-
ν
2
1
ν
0
ν
1
0
0
0
1
-
ν
2
,where ν denotes Poisson’s ratio. In general, the optimal solution obtained by direct implementation of the method is prone to problems with checkerboards and mesh dependency. A large number of restriction schemes have been proposed to ensure manufacturability and mesh independency, such as perimeter control method, mesh-independent filtering methods, and Heaviside projection. A comprehensive comparison of the existing restriction schemes is provided in the paper of Sigmund [23].


## 2.2. Robust Topology Optimization Formulation

The topology optimization of linear elastic structures subject to loading uncertainties is considered in this paper. Due to the variation of uncertain parametersξ, the performance function f
(
ρ
,
ξ
) has probability distribution description. Thus, the robustness of the design objective can be achieved by optimizing the mean performance and minimizing the standard deviation of the performance simultaneously, as shown in Figure 2. In Figure 2, point A denotes the deterministic optimum not considering the uncertainty, while point B denotes the robust optimum. The performance of deterministic optimum is better than that of the robust optimum. However, its variation is wider than the robust optimum.Figure 2
The concept of robust design model.Therefore, a main challenge in performing RTO for realistic problems is the implementation of an algorithm, which inherently represents a multiobjective optimization problem. One efficient way is to combine all objective functions into a single weighted function. In mathematical terms, a general robust topology optimization task can be stated as(7)
find:
ρ
min:
μ
f
ρ
,
ξ
+
κ
σ
f
ρ
,
ξ
s.t.:
V
ρ
V
∗
-
1
≤
0
ρ
min
⁡
≤
ρ
≤
ρ
max
⁡
,where κ is the weight to be determined by the designer that satisfies κ
⩾
0. Note that the weights have to be set a priori, which make a trade-off decision between an optimal mean performance and minimization of standard deviation of the performance to obtain a compromise solution. μ
f and σ
f are the mean value and the standard deviation of the performance function, respectively. V
∗ is the material resource constraint. The design variables ρ are bounded by their lower and upper bounds ρ
m
i
n and ρ
m
a
x, respectively. In a probabilistic setting, the mean value μ
f and the standard deviation σ
f can be expressed as(8)
μ
f
ρ
,
ξ
=
E
f
ρ
,
ξ
=
∫
Ω
f
ξ
p
ξ
d
ξ
,
σ
f
ρ
,
ξ
=
E
f
ρ
,
ξ
-
μ
f
ρ
,
ξ
2
=
∫
Ω
f
ξ
-
μ
f
2
p
ξ
d
ξ
1
/
2
,where E is the expectation operator and p
(
ξ
) is the joint probability density function (PDF) of the random variables. In practice, it is very difficult and even impossible to calculate the mean value and the variation of the performance function through the multidimensional integration in (8). Such difficulties have motivated the development of various numerical methods such as simulation-based methods, local expansion-based methods, and the stochastic collocation methods. The reader is referred to Lee and Chen [24] for a more detailed review.


## 3. Stochastic Collocation Method

The main idea of stochastic collocation method is to sample the quantity of interest at specific collocation points in the random space, where the sample points are deterministic and are associated with quadrature formula for the evaluation of integral statistics such as mean and standard deviation. Examples of such quadrature rules include the Newton-Cotes (midpoint, rectangle, and trapezoidal) rule, the Gauss quadrature (Legendre, Chebyshev, Laguerre, Hermite, Jacobi, Kronrod, and Patterson) rule, and the Clenshaw-Curtis rule [25–27]. In the case of a one-dimensional random space, these points are quadrature abscissas selected according to the probability measure of the random variable. Tensor product grid and sparse grid constructions are then utilized to generate the multidimensional abscissas.


### 3.1. Tensor Product Grid

Without loss of generality, consider the class of functionf
(
ξ
) over the d-dimensional unit hypercube Ω
∶
=
-
1,1
d with numerical integration representation(9)
I
d
f
=
∫
Ω
f
ξ
d
ξ
.In the one-dimensional case, a sequence of univariate quadrature operators by m
l
1-point quadrature formulas with level l
∈
N is(10)
Q
l
1
f
=
∑
i
=
1
m
l
1
f
ξ
l
i
w
l
i
,where ξ
l
i and w
l
i are the collocation points and weights, respectively. To obtain a quadrature formula for the multivariate case d
>
1, the full tensor product formula is defined as(11)
Q
l
1
1
⊗
⋯
⊗
Q
l
d
1
f
=
∑
l
1
=
1
m
l
1
1
⋯
∑
l
d
=
1
m
l
d
1
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
·
w
l
1
i
1
⊗
⋯
⊗
w
l
d
i
d
.Then, the mean and standard deviation of the objective performance are derived from the following equations:(12)
μ
f
=
E
f
=
∫
Ω
f
ξ
p
ξ
d
ξ
=
∑
l
1
=
1
m
l
1
1
w
l
1
i
1
⋯
∑
l
d
=
1
m
l
d
1
w
l
d
i
d
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
,
σ
f
=
E
f
-
E
f
2
=
∫
Ω
f
ξ
-
μ
f
2
p
ξ
d
ξ
1
/
2
=
∑
l
1
=
1
m
l
1
1
w
l
1
i
1
⋯
∑
l
d
=
1
m
l
d
1
w
l
d
i
d
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
-
μ
f
2
1
/
2
.Clearly, the tensor product grid requires M
=
(
m
l
1
1
,
…
,
m
l
d
1
) collocation points, which are sampled on the full grid. Assuming the same level accuracy in each dimension—that is, m
l
1
1
=
⋯
=
m
l
d
1
=
m—the total number of points is M
=
m
d. Thus, the computational burden increases significantly with the dimension d. This is often referred to as thecurse of dimensionality [28].


### 3.2. Sparse Grid

The sparse grid method is a special discretization technique, which can be traced back to the Smolyak algorithm [29]. It is based on hierarchical basis, a representation of a discrete function space which is equivalent to the conventional nodal basis, and a sparse tensor product construction [30, 31]. The readers are referred to the review paper by Xiong et al. [32] and the references therein for more information. The one-dimensional difference quadrature formula is defined as (13)
Δ
k
1
f
=
Q
k
1
-
Q
k
-
1
1
f
,
with
Q
0
1
f
=
0
.Then the Smolyak quadrature formula for d-dimensional functions f with level l
∈
N is given by (14)
Q
l
d
f
=
∑
k
≤
l
+
d
Δ
k
1
1
⊗
⋯
⊗
Δ
k
d
1
f
.Note that |
k
| denotes the summation of the multi-indices (
|
k
|
=
k
1
+
⋯
+
k
d
). Alternatively, the above formula can be written as(15)
Q
l
d
f
=
∑
l
+
1
≤
k
≤
l
+
d
-
1
l
+
d
-
k
·
d
-
1
l
+
d
-
k
·
Q
k
1
1
⊗
⋯
⊗
Q
k
d
1
f
.To compute Q
l
d
(
f
), a specific tensor product rule of sparse grid is needed, defined as (16)
U
l
d
=
⋃
l
+
1
≤
k
≤
l
+
d
U
k
1
1
⊗
⋯
⊗
U
k
d
1
,where U
k
1
1 denotes the set of collocation points for one-dimensional l-level accuracy. With the Smolyak algorithm, the weight w
i corresponding to the ith collocation point ξ
i is defined as (17)
w
i
=
-
1
l
+
d
-
k
d
-
1
l
+
d
-
k
w
k
1
i
1
⊗
⋯
⊗
w
k
d
i
d
.Then, the mean and standard deviation of the objective performance by sparse grid method can be computed by(18)
μ
f
=
E
f
=
∑
i
=
1
n
l
d
w
k
f
ξ
k
,
σ
f
=
E
f
-
E
f
2
=
∑
i
=
1
n
l
d
w
k
f
ξ
k
-
μ
f
2
1
/
2
,where m
l
d is the total number of points in sparse grid that is determined as(19)
m
l
d
=
⋃
l
+
1
≤
k
≤
l
+
d
m
k
1
1
⊗
⋯
⊗
m
k
d
1
.Suppose that the same quadrature rule is chosen under the dimension d with same level in each dimension, and m represents the collocation points in each dimension and the tensor-product grid scales as m
d, whereas the sparse grid scales as m
l
o
g
d, significantly mitigating thecurse of dimensionality. Accurate quadrature can be created by using higher values of the level l; however, this also leads to more sampling points.The sparse grid method discretizes the stochastic space in hierarchical structure based on nested collocation nodes, such as the Chebyshev nodes or Gauss-Patterson nodes, leading to Clenshaw-Curtis rule or Gauss-Patterson rule, respectively. Here, the Clenshaw-Curtis type sparse grid designH
C
C, with nonequidistant nodes, is constructed [33–35]. The collocation points ξ
k
i are defined as(20)
ξ
k
i
=
-
cos
⁡
π
i
-
1
m
k
-
1
,
for
i
=
1
,
…
,
m
k
,
if
m
k
>
1
,
0
,
if
m
k
=
1
.The sequence m
k is chosen as (21)
m
k
=
1
,
if
k
=
1
,
2
k
-
1
+
1
,
if
k
>
1and the corresponding weights are given by(22)
w
k
i
=
1
m
k
m
k
-
2
,
for
i
=
1
,
m
k
,
if
m
k
>
1
,
2
m
k
-
1
1
-
cos
⁡
π
i
-
1
m
k
m
k
-
2
-
2
∑
j
=
1
m
k
-
3
/
2
1
4
j
2
-
1
·
cos
⁡
2
π
j
i
-
1
m
k
-
1
,
for
i
=
2
,
…
,
m
k
-
1
,
if
m
k
>
1
,
2
,
if
m
k
=
1
.Examples of two-dimensional collocation points of full tensor product grid and Smolyak sparse grid based on Clenshaw-Curtis rule with levell
=
3 are shown in Figure 3. It is clearly observed that, compared with full tensor product grid, the sparse grid dramatically reduces the number of collocation points, while preserving a high level of accuracy.Figure 3
Two-dimensional collocation points and weights by Clenshaw-Curtis rule in full tensor product grid (a) and Smolyak sparse grid (b) with levell
=
1,2, and 3.
(a)
Tensor product grid
(b)
Sparse grid


## 3.1. Tensor Product Grid

Without loss of generality, consider the class of functionf
(
ξ
) over the d-dimensional unit hypercube Ω
∶
=
-
1,1
d with numerical integration representation(9)
I
d
f
=
∫
Ω
f
ξ
d
ξ
.In the one-dimensional case, a sequence of univariate quadrature operators by m
l
1-point quadrature formulas with level l
∈
N is(10)
Q
l
1
f
=
∑
i
=
1
m
l
1
f
ξ
l
i
w
l
i
,where ξ
l
i and w
l
i are the collocation points and weights, respectively. To obtain a quadrature formula for the multivariate case d
>
1, the full tensor product formula is defined as(11)
Q
l
1
1
⊗
⋯
⊗
Q
l
d
1
f
=
∑
l
1
=
1
m
l
1
1
⋯
∑
l
d
=
1
m
l
d
1
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
·
w
l
1
i
1
⊗
⋯
⊗
w
l
d
i
d
.Then, the mean and standard deviation of the objective performance are derived from the following equations:(12)
μ
f
=
E
f
=
∫
Ω
f
ξ
p
ξ
d
ξ
=
∑
l
1
=
1
m
l
1
1
w
l
1
i
1
⋯
∑
l
d
=
1
m
l
d
1
w
l
d
i
d
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
,
σ
f
=
E
f
-
E
f
2
=
∫
Ω
f
ξ
-
μ
f
2
p
ξ
d
ξ
1
/
2
=
∑
l
1
=
1
m
l
1
1
w
l
1
i
1
⋯
∑
l
d
=
1
m
l
d
1
w
l
d
i
d
f
ξ
l
1
i
1
,
…
,
ξ
l
d
i
d
-
μ
f
2
1
/
2
.Clearly, the tensor product grid requires M
=
(
m
l
1
1
,
…
,
m
l
d
1
) collocation points, which are sampled on the full grid. Assuming the same level accuracy in each dimension—that is, m
l
1
1
=
⋯
=
m
l
d
1
=
m—the total number of points is M
=
m
d. Thus, the computational burden increases significantly with the dimension d. This is often referred to as thecurse of dimensionality [28].


## 3.2. Sparse Grid

The sparse grid method is a special discretization technique, which can be traced back to the Smolyak algorithm [29]. It is based on hierarchical basis, a representation of a discrete function space which is equivalent to the conventional nodal basis, and a sparse tensor product construction [30, 31]. The readers are referred to the review paper by Xiong et al. [32] and the references therein for more information. The one-dimensional difference quadrature formula is defined as (13)
Δ
k
1
f
=
Q
k
1
-
Q
k
-
1
1
f
,
with
Q
0
1
f
=
0
.Then the Smolyak quadrature formula for d-dimensional functions f with level l
∈
N is given by (14)
Q
l
d
f
=
∑
k
≤
l
+
d
Δ
k
1
1
⊗
⋯
⊗
Δ
k
d
1
f
.Note that |
k
| denotes the summation of the multi-indices (
|
k
|
=
k
1
+
⋯
+
k
d
). Alternatively, the above formula can be written as(15)
Q
l
d
f
=
∑
l
+
1
≤
k
≤
l
+
d
-
1
l
+
d
-
k
·
d
-
1
l
+
d
-
k
·
Q
k
1
1
⊗
⋯
⊗
Q
k
d
1
f
.To compute Q
l
d
(
f
), a specific tensor product rule of sparse grid is needed, defined as (16)
U
l
d
=
⋃
l
+
1
≤
k
≤
l
+
d
U
k
1
1
⊗
⋯
⊗
U
k
d
1
,where U
k
1
1 denotes the set of collocation points for one-dimensional l-level accuracy. With the Smolyak algorithm, the weight w
i corresponding to the ith collocation point ξ
i is defined as (17)
w
i
=
-
1
l
+
d
-
k
d
-
1
l
+
d
-
k
w
k
1
i
1
⊗
⋯
⊗
w
k
d
i
d
.Then, the mean and standard deviation of the objective performance by sparse grid method can be computed by(18)
μ
f
=
E
f
=
∑
i
=
1
n
l
d
w
k
f
ξ
k
,
σ
f
=
E
f
-
E
f
2
=
∑
i
=
1
n
l
d
w
k
f
ξ
k
-
μ
f
2
1
/
2
,where m
l
d is the total number of points in sparse grid that is determined as(19)
m
l
d
=
⋃
l
+
1
≤
k
≤
l
+
d
m
k
1
1
⊗
⋯
⊗
m
k
d
1
.Suppose that the same quadrature rule is chosen under the dimension d with same level in each dimension, and m represents the collocation points in each dimension and the tensor-product grid scales as m
d, whereas the sparse grid scales as m
l
o
g
d, significantly mitigating thecurse of dimensionality. Accurate quadrature can be created by using higher values of the level l; however, this also leads to more sampling points.The sparse grid method discretizes the stochastic space in hierarchical structure based on nested collocation nodes, such as the Chebyshev nodes or Gauss-Patterson nodes, leading to Clenshaw-Curtis rule or Gauss-Patterson rule, respectively. Here, the Clenshaw-Curtis type sparse grid designH
C
C, with nonequidistant nodes, is constructed [33–35]. The collocation points ξ
k
i are defined as(20)
ξ
k
i
=
-
cos
⁡
π
i
-
1
m
k
-
1
,
for
i
=
1
,
…
,
m
k
,
if
m
k
>
1
,
0
,
if
m
k
=
1
.The sequence m
k is chosen as (21)
m
k
=
1
,
if
k
=
1
,
2
k
-
1
+
1
,
if
k
>
1and the corresponding weights are given by(22)
w
k
i
=
1
m
k
m
k
-
2
,
for
i
=
1
,
m
k
,
if
m
k
>
1
,
2
m
k
-
1
1
-
cos
⁡
π
i
-
1
m
k
m
k
-
2
-
2
∑
j
=
1
m
k
-
3
/
2
1
4
j
2
-
1
·
cos
⁡
2
π
j
i
-
1
m
k
-
1
,
for
i
=
2
,
…
,
m
k
-
1
,
if
m
k
>
1
,
2
,
if
m
k
=
1
.Examples of two-dimensional collocation points of full tensor product grid and Smolyak sparse grid based on Clenshaw-Curtis rule with levell
=
3 are shown in Figure 3. It is clearly observed that, compared with full tensor product grid, the sparse grid dramatically reduces the number of collocation points, while preserving a high level of accuracy.Figure 3
Two-dimensional collocation points and weights by Clenshaw-Curtis rule in full tensor product grid (a) and Smolyak sparse grid (b) with levell
=
1,2, and 3.
(a)
Tensor product grid
(b)
Sparse grid


## 4. RTO Algorithm and Numerical Examples

### 4.1. Optimization Algorithm

The implementation of optimization algorithm for RTO is based on the commercial software package Optistruct 11.0 of Altair Hyperworks, where the solid isotropic material with penalization is employed as material interpolation model and convex programming method as optimization solver. In addition, the coordinates and corresponding weights of the collocation points for random variables are evaluated in an external environment of Optistruct for auxiliary conditions. The optimization algorithm can be summarized as follows:(1)
State problem discretization and initialization, including initial design domainΩ and boundary conditions Г.
(2)
Identify the random loading uncertaintyξ as well as the related loading angle and magnitude probability distributions.
(3)
Compute collocation pointsξ
k
i and the corresponding weights w
k
i by tensor product gird and sparse grid, according to the dimension d and level l information.
(4)
Define multiple loading cases through mapping the loading uncertainty with the collocation points.
(5)
Establish the meanμ
f and standard deviation σ
f of the performance function through specifying the user-defined equations, according to the prescribed format.
(6)
Implement the weighted function by combining the mean and standard deviation functions with predefined weightκ, and choose the function as the objective function to be minimized.
(7)
Define the manufacturing restriction schemes, optimization strategy, and volume fraction constraint.
(8)
Run the optimization loop until the convergence criterion is satisfied.


### 4.2. Numerical Examples

In this section, we will illustrate the proposed RTO approach using two numerical example problems—a 2D simple supported beam and a 3D L-shaped structure—to demonstrate the effectiveness and efficiency for tensor product grid and sparse grid. A simple Monte-Carlo simulation is also performed to verify the accuracy of the statistical moment calculated by the proposed approach. The procedure of topology optimization is conducted on a Windows 7 workstation (Intel Xeon (R) E5440, 2.83 GHz, 4.00 GB RAM, 8 cores).


#### 4.2.1. A 2D Simply Supported Beam

As the first example, the design domain, the boundary, and loading conditions of the simply supported beam are illustrated in Figure4. The dimensions of the beam are L
=
90 mm and H
=
30 mm and the thickness is T
=
1 mm. The isotropic material is assumed to have Young’s modulus E
0 of 1 MPa and Poisson’s ratio ν of 0.30. A single random load case is considered with which the angle and magnitude are characterized by two independent random variables. The angle is assumed to follow a uniform distribution with the interval of [
-
3
π
/
4
,
-
π
/
4
]. The magnitude follows a Gaussian distribution with mean of 1 and standard deviation of 0.3. In robust design, the optimization problem is to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint, where the weight κ is set to be 1. The design domain is discretized by 2700 (90 × 30) four-node quadrilateral elements.Figure 4
Design domain, boundary conditions, and loading conditions of a 2D simply supported beam.The optimal topologies obtained by deterministic and the proposed robust topology optimization approaches are presented in Figures5 and 6, and the corresponding results are compared in Table 1, including the required number of points, statistical information of the performance function, volume fraction constraint, iterations, and computing time. Here, full tensor product grid and sparse grid under level l
=
1,2, and 3 are investigated in this section. Extensive 1000 Monte Carlo simulations are also performed to verify the accuracy of the proposed robust approach.Table 1
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
22.36
12.03
0.25
54
35 s
Monte Carlo simulation
1000
21.37
11.13
0.25
55
302 s
Robust
Tensor product grid (l
=
1)
9
22.26
11.85
0.25
51
59 s
Tensor product grid (l
=
2)
25
21.87
11.66
0.25
54
90 s
Tensor product grid (l
=
3)
81
21.77
11.34
0.25
54
147 s
Sparse grid (l
=
1)
5
22.23
11.84
0.25
53
59 s
Sparse grid (l
=
2)
13
21.92
11.70
0.25
54
68 s
Sparse grid (l
=
3)
29
21.52
11.27
0.25
54
71 sFigure 5
Optimal designs of the simply supported beam.
(a)
Deterministic design
(b)
Monte Carlo simulationFigure 6
Robust designs of the simple supported beam.
(a)
Tensor product formula(
l
=
1
)
(b)
Sparse grid method(
l
=
1
)
(c)
Tensor product formula(
l
=
2
)
(d)
Sparse grid method(
l
=
2
)
(e)
Tensor product formula(
l
=
3
)
(f)
Sparse grid method(
l
=
3
)For comparison purposes, the deterministic design is performed using the mean values of the load and the robust design is conducted with multiple load cases due to randomness of load angle and amplitude. The DTO and RTO resulted in significantly different topologies, as shown in Figure5. The remarkable difference lies in the robust design exhibiting an asymmetric layout compared with the deterministic design characterized by symmetric configuration. It can be attributed to the fact that the deterministic design suffers from symmetric vertical constraints at the lower left and the lower right corner. However, in robust design, the structure possesses an asymmetrical degree of freedom constraint in horizontal direction under the condition of horizontal component of load.As the results listed in Table1 show, the robust designs present more robust topologies than the deterministic approach, according to the results summary of the statistical information. This emphasizes the importance of considering the uncertainties in load for structure design. It is also noted that the corresponding computation cost for robust designs is larger than for the deterministic design since multiple load cases are required for calculating the mean and deviation values of the performance function. The results of Monte Carlo simulations in Table 1 show that the proposed robust approach possesses improved efficiency without losing much accuracy.From the comparison of the results, with the improvement of levell accuracy, the number of collocation points in tensor product grid and sparse grid corresponding increases. As the results listed in Table 1 show, the mean and standard deviation of the compliance for sparse grid and tensor product grid are gradually close to those for Monte Carlo simulation under the increase of level l. This indicates that more robust designs are clearly obtained under higher level l. Meanwhile, the higher computation cost is required due to the increase in the number of evaluations of performance function responses. This means robustness enhancement is inevitably accompanied with the sacrifice of convergence time. It can also be seen from the results listed in Table 1 that, under the same level l, sparse grid uses fewer collocation points than those of tensor product grid to produce similar even more robust results. At the same time, compared with sparse grid, the computation burden for the exponential dependent calculation of tensor product grid can be overcome to a certain extent. This means sparse grid, in a manner, performs better than the tensor product grid in compromise between robustness and cost in robust design.


#### 4.2.2. A 3D L-Shaped Structure

The design domain and loading conditions of this example are illustrated in Figure7. The geometrical dimensions of the structure are L
=
120 mm and T
=
30 mm. Young’s modulus and Poisson’s ratio of the fully solid material are 2.1 × 105 MPa and 0.30, respectively. Two concentrated forces, F
Y and F
Z, are applied at central point of right bottom surface simultaneously. The magnitudes of the forces are characterized by two uncorrelated random variables, which follow Gaussian distributions with mean of 1000 N and 30% variance from the mean values. The optimization problem is still to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint. The design domain is divided into 1.2 × 104 eight-node hexahedron elements.Figure 7
Design domain of the 3D L-shaped structure.The optimal topologies obtained by deterministic, Monte Carlo simulation and the robust approaches with density value above 0.3 thresholds are presented in Figures8 and 9, respectively. The corresponding optimization results are summarized in Table 2, including the required number of points, statistical moments, volume constraint, iterations, and computing time. Meanwhile, full tensor product grid and sparse grid are investigated under level l
=
1,2, and 3.Table 2
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
241.08
131.42
0.06
44
7 min 14 s
Monte Carlo simulation
1000
182.20
111.34
0.06
51
28 min 27 s
Robust
Tensor product formula (l
=
1)
9
183.95
121.71
0.06
46
8 min 19 s
Tensor product formula (l
=
2)
25
182.62
119.30
0.06
49
9 min 44 s
Tensor product formula (l
=
3)
81
182.91
113.72
0.06
49
17 min 26 s
Sparse grid method (l
=
1)
5
183.20
122.66
0.06
47
8 min 18 s
Sparse grid method (l
=
2)
13
182.71
119.24
0.06
50
8 min 47 s
Sparse grid method (l
=
3)
29
182.93
114.02
0.06
50
10 min 25 sFigure 8
Optimal topologies of the L-shaped structure.
(a)
Deterministic design
(b)
Monte Carlo simulation designFigure 9
Optimal designs of the L-shaped structure by robust topology optimization.
(a)
Tensor product grid(
l
=
1
)
(b)
Sparse grid(
l
=
1
)
(c)
Tensor product grid(
l
=
2
)
(d)
Sparse grid(
l
=
2
)
(e)
Tensor product grid(
l
=
3
)
(f)
Sparse grid(
l
=
3
)Comparing the deterministic design with the proposed robust designs, obvious distinctions can be clearly observed in optimal topologies as shown in Figures8 and 9. From the viewpoint of geometry, the robust designs represent more rational material layouts, which lead to more robust structural configuration. According to Table 2, the final statistical moments of the robust designs are superior to that of deterministic design. This indicates that the robust design is less sensitive to the variations of load magnitude. It can also be seen from the results listed in Table 2 that as the level l is increasing, the mean and standard deviation of the structure compliance for tensor product grid and sparse grid are closer to those for Monte Carlo simulation in which more robust designs are obtained, but it is also revealed that the corresponding computing time is increasing. This means robustness and cost are conflicting criteria in robust topology optimization design. For comparison of the results of tensor product grid and sparse grid, these indicate that sparse grid, in a manner, performs better than the tensor product grid comprehensively considering the accuracy and efficiency.


## 4.1. Optimization Algorithm

The implementation of optimization algorithm for RTO is based on the commercial software package Optistruct 11.0 of Altair Hyperworks, where the solid isotropic material with penalization is employed as material interpolation model and convex programming method as optimization solver. In addition, the coordinates and corresponding weights of the collocation points for random variables are evaluated in an external environment of Optistruct for auxiliary conditions. The optimization algorithm can be summarized as follows:(1)
State problem discretization and initialization, including initial design domainΩ and boundary conditions Г.
(2)
Identify the random loading uncertaintyξ as well as the related loading angle and magnitude probability distributions.
(3)
Compute collocation pointsξ
k
i and the corresponding weights w
k
i by tensor product gird and sparse grid, according to the dimension d and level l information.
(4)
Define multiple loading cases through mapping the loading uncertainty with the collocation points.
(5)
Establish the meanμ
f and standard deviation σ
f of the performance function through specifying the user-defined equations, according to the prescribed format.
(6)
Implement the weighted function by combining the mean and standard deviation functions with predefined weightκ, and choose the function as the objective function to be minimized.
(7)
Define the manufacturing restriction schemes, optimization strategy, and volume fraction constraint.
(8)
Run the optimization loop until the convergence criterion is satisfied.


## 4.2. Numerical Examples

In this section, we will illustrate the proposed RTO approach using two numerical example problems—a 2D simple supported beam and a 3D L-shaped structure—to demonstrate the effectiveness and efficiency for tensor product grid and sparse grid. A simple Monte-Carlo simulation is also performed to verify the accuracy of the statistical moment calculated by the proposed approach. The procedure of topology optimization is conducted on a Windows 7 workstation (Intel Xeon (R) E5440, 2.83 GHz, 4.00 GB RAM, 8 cores).


### 4.2.1. A 2D Simply Supported Beam

As the first example, the design domain, the boundary, and loading conditions of the simply supported beam are illustrated in Figure4. The dimensions of the beam are L
=
90 mm and H
=
30 mm and the thickness is T
=
1 mm. The isotropic material is assumed to have Young’s modulus E
0 of 1 MPa and Poisson’s ratio ν of 0.30. A single random load case is considered with which the angle and magnitude are characterized by two independent random variables. The angle is assumed to follow a uniform distribution with the interval of [
-
3
π
/
4
,
-
π
/
4
]. The magnitude follows a Gaussian distribution with mean of 1 and standard deviation of 0.3. In robust design, the optimization problem is to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint, where the weight κ is set to be 1. The design domain is discretized by 2700 (90 × 30) four-node quadrilateral elements.Figure 4
Design domain, boundary conditions, and loading conditions of a 2D simply supported beam.The optimal topologies obtained by deterministic and the proposed robust topology optimization approaches are presented in Figures5 and 6, and the corresponding results are compared in Table 1, including the required number of points, statistical information of the performance function, volume fraction constraint, iterations, and computing time. Here, full tensor product grid and sparse grid under level l
=
1,2, and 3 are investigated in this section. Extensive 1000 Monte Carlo simulations are also performed to verify the accuracy of the proposed robust approach.Table 1
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
22.36
12.03
0.25
54
35 s
Monte Carlo simulation
1000
21.37
11.13
0.25
55
302 s
Robust
Tensor product grid (l
=
1)
9
22.26
11.85
0.25
51
59 s
Tensor product grid (l
=
2)
25
21.87
11.66
0.25
54
90 s
Tensor product grid (l
=
3)
81
21.77
11.34
0.25
54
147 s
Sparse grid (l
=
1)
5
22.23
11.84
0.25
53
59 s
Sparse grid (l
=
2)
13
21.92
11.70
0.25
54
68 s
Sparse grid (l
=
3)
29
21.52
11.27
0.25
54
71 sFigure 5
Optimal designs of the simply supported beam.
(a)
Deterministic design
(b)
Monte Carlo simulationFigure 6
Robust designs of the simple supported beam.
(a)
Tensor product formula(
l
=
1
)
(b)
Sparse grid method(
l
=
1
)
(c)
Tensor product formula(
l
=
2
)
(d)
Sparse grid method(
l
=
2
)
(e)
Tensor product formula(
l
=
3
)
(f)
Sparse grid method(
l
=
3
)For comparison purposes, the deterministic design is performed using the mean values of the load and the robust design is conducted with multiple load cases due to randomness of load angle and amplitude. The DTO and RTO resulted in significantly different topologies, as shown in Figure5. The remarkable difference lies in the robust design exhibiting an asymmetric layout compared with the deterministic design characterized by symmetric configuration. It can be attributed to the fact that the deterministic design suffers from symmetric vertical constraints at the lower left and the lower right corner. However, in robust design, the structure possesses an asymmetrical degree of freedom constraint in horizontal direction under the condition of horizontal component of load.As the results listed in Table1 show, the robust designs present more robust topologies than the deterministic approach, according to the results summary of the statistical information. This emphasizes the importance of considering the uncertainties in load for structure design. It is also noted that the corresponding computation cost for robust designs is larger than for the deterministic design since multiple load cases are required for calculating the mean and deviation values of the performance function. The results of Monte Carlo simulations in Table 1 show that the proposed robust approach possesses improved efficiency without losing much accuracy.From the comparison of the results, with the improvement of levell accuracy, the number of collocation points in tensor product grid and sparse grid corresponding increases. As the results listed in Table 1 show, the mean and standard deviation of the compliance for sparse grid and tensor product grid are gradually close to those for Monte Carlo simulation under the increase of level l. This indicates that more robust designs are clearly obtained under higher level l. Meanwhile, the higher computation cost is required due to the increase in the number of evaluations of performance function responses. This means robustness enhancement is inevitably accompanied with the sacrifice of convergence time. It can also be seen from the results listed in Table 1 that, under the same level l, sparse grid uses fewer collocation points than those of tensor product grid to produce similar even more robust results. At the same time, compared with sparse grid, the computation burden for the exponential dependent calculation of tensor product grid can be overcome to a certain extent. This means sparse grid, in a manner, performs better than the tensor product grid in compromise between robustness and cost in robust design.


### 4.2.2. A 3D L-Shaped Structure

The design domain and loading conditions of this example are illustrated in Figure7. The geometrical dimensions of the structure are L
=
120 mm and T
=
30 mm. Young’s modulus and Poisson’s ratio of the fully solid material are 2.1 × 105 MPa and 0.30, respectively. Two concentrated forces, F
Y and F
Z, are applied at central point of right bottom surface simultaneously. The magnitudes of the forces are characterized by two uncorrelated random variables, which follow Gaussian distributions with mean of 1000 N and 30% variance from the mean values. The optimization problem is still to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint. The design domain is divided into 1.2 × 104 eight-node hexahedron elements.Figure 7
Design domain of the 3D L-shaped structure.The optimal topologies obtained by deterministic, Monte Carlo simulation and the robust approaches with density value above 0.3 thresholds are presented in Figures8 and 9, respectively. The corresponding optimization results are summarized in Table 2, including the required number of points, statistical moments, volume constraint, iterations, and computing time. Meanwhile, full tensor product grid and sparse grid are investigated under level l
=
1,2, and 3.Table 2
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
241.08
131.42
0.06
44
7 min 14 s
Monte Carlo simulation
1000
182.20
111.34
0.06
51
28 min 27 s
Robust
Tensor product formula (l
=
1)
9
183.95
121.71
0.06
46
8 min 19 s
Tensor product formula (l
=
2)
25
182.62
119.30
0.06
49
9 min 44 s
Tensor product formula (l
=
3)
81
182.91
113.72
0.06
49
17 min 26 s
Sparse grid method (l
=
1)
5
183.20
122.66
0.06
47
8 min 18 s
Sparse grid method (l
=
2)
13
182.71
119.24
0.06
50
8 min 47 s
Sparse grid method (l
=
3)
29
182.93
114.02
0.06
50
10 min 25 sFigure 8
Optimal topologies of the L-shaped structure.
(a)
Deterministic design
(b)
Monte Carlo simulation designFigure 9
Optimal designs of the L-shaped structure by robust topology optimization.
(a)
Tensor product grid(
l
=
1
)
(b)
Sparse grid(
l
=
1
)
(c)
Tensor product grid(
l
=
2
)
(d)
Sparse grid(
l
=
2
)
(e)
Tensor product grid(
l
=
3
)
(f)
Sparse grid(
l
=
3
)Comparing the deterministic design with the proposed robust designs, obvious distinctions can be clearly observed in optimal topologies as shown in Figures8 and 9. From the viewpoint of geometry, the robust designs represent more rational material layouts, which lead to more robust structural configuration. According to Table 2, the final statistical moments of the robust designs are superior to that of deterministic design. This indicates that the robust design is less sensitive to the variations of load magnitude. It can also be seen from the results listed in Table 2 that as the level l is increasing, the mean and standard deviation of the structure compliance for tensor product grid and sparse grid are closer to those for Monte Carlo simulation in which more robust designs are obtained, but it is also revealed that the corresponding computing time is increasing. This means robustness and cost are conflicting criteria in robust topology optimization design. For comparison of the results of tensor product grid and sparse grid, these indicate that sparse grid, in a manner, performs better than the tensor product grid comprehensively considering the accuracy and efficiency.


## 4.2.1. A 2D Simply Supported Beam

As the first example, the design domain, the boundary, and loading conditions of the simply supported beam are illustrated in Figure4. The dimensions of the beam are L
=
90 mm and H
=
30 mm and the thickness is T
=
1 mm. The isotropic material is assumed to have Young’s modulus E
0 of 1 MPa and Poisson’s ratio ν of 0.30. A single random load case is considered with which the angle and magnitude are characterized by two independent random variables. The angle is assumed to follow a uniform distribution with the interval of [
-
3
π
/
4
,
-
π
/
4
]. The magnitude follows a Gaussian distribution with mean of 1 and standard deviation of 0.3. In robust design, the optimization problem is to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint, where the weight κ is set to be 1. The design domain is discretized by 2700 (90 × 30) four-node quadrilateral elements.Figure 4
Design domain, boundary conditions, and loading conditions of a 2D simply supported beam.The optimal topologies obtained by deterministic and the proposed robust topology optimization approaches are presented in Figures5 and 6, and the corresponding results are compared in Table 1, including the required number of points, statistical information of the performance function, volume fraction constraint, iterations, and computing time. Here, full tensor product grid and sparse grid under level l
=
1,2, and 3 are investigated in this section. Extensive 1000 Monte Carlo simulations are also performed to verify the accuracy of the proposed robust approach.Table 1
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
22.36
12.03
0.25
54
35 s
Monte Carlo simulation
1000
21.37
11.13
0.25
55
302 s
Robust
Tensor product grid (l
=
1)
9
22.26
11.85
0.25
51
59 s
Tensor product grid (l
=
2)
25
21.87
11.66
0.25
54
90 s
Tensor product grid (l
=
3)
81
21.77
11.34
0.25
54
147 s
Sparse grid (l
=
1)
5
22.23
11.84
0.25
53
59 s
Sparse grid (l
=
2)
13
21.92
11.70
0.25
54
68 s
Sparse grid (l
=
3)
29
21.52
11.27
0.25
54
71 sFigure 5
Optimal designs of the simply supported beam.
(a)
Deterministic design
(b)
Monte Carlo simulationFigure 6
Robust designs of the simple supported beam.
(a)
Tensor product formula(
l
=
1
)
(b)
Sparse grid method(
l
=
1
)
(c)
Tensor product formula(
l
=
2
)
(d)
Sparse grid method(
l
=
2
)
(e)
Tensor product formula(
l
=
3
)
(f)
Sparse grid method(
l
=
3
)For comparison purposes, the deterministic design is performed using the mean values of the load and the robust design is conducted with multiple load cases due to randomness of load angle and amplitude. The DTO and RTO resulted in significantly different topologies, as shown in Figure5. The remarkable difference lies in the robust design exhibiting an asymmetric layout compared with the deterministic design characterized by symmetric configuration. It can be attributed to the fact that the deterministic design suffers from symmetric vertical constraints at the lower left and the lower right corner. However, in robust design, the structure possesses an asymmetrical degree of freedom constraint in horizontal direction under the condition of horizontal component of load.As the results listed in Table1 show, the robust designs present more robust topologies than the deterministic approach, according to the results summary of the statistical information. This emphasizes the importance of considering the uncertainties in load for structure design. It is also noted that the corresponding computation cost for robust designs is larger than for the deterministic design since multiple load cases are required for calculating the mean and deviation values of the performance function. The results of Monte Carlo simulations in Table 1 show that the proposed robust approach possesses improved efficiency without losing much accuracy.From the comparison of the results, with the improvement of levell accuracy, the number of collocation points in tensor product grid and sparse grid corresponding increases. As the results listed in Table 1 show, the mean and standard deviation of the compliance for sparse grid and tensor product grid are gradually close to those for Monte Carlo simulation under the increase of level l. This indicates that more robust designs are clearly obtained under higher level l. Meanwhile, the higher computation cost is required due to the increase in the number of evaluations of performance function responses. This means robustness enhancement is inevitably accompanied with the sacrifice of convergence time. It can also be seen from the results listed in Table 1 that, under the same level l, sparse grid uses fewer collocation points than those of tensor product grid to produce similar even more robust results. At the same time, compared with sparse grid, the computation burden for the exponential dependent calculation of tensor product grid can be overcome to a certain extent. This means sparse grid, in a manner, performs better than the tensor product grid in compromise between robustness and cost in robust design.


## 4.2.2. A 3D L-Shaped Structure

The design domain and loading conditions of this example are illustrated in Figure7. The geometrical dimensions of the structure are L
=
120 mm and T
=
30 mm. Young’s modulus and Poisson’s ratio of the fully solid material are 2.1 × 105 MPa and 0.30, respectively. Two concentrated forces, F
Y and F
Z, are applied at central point of right bottom surface simultaneously. The magnitudes of the forces are characterized by two uncorrelated random variables, which follow Gaussian distributions with mean of 1000 N and 30% variance from the mean values. The optimization problem is still to minimize the weighted summation of the mean μ
f and standard deviation σ
f of the compliance of the structure under the volume constraint. The design domain is divided into 1.2 × 104 eight-node hexahedron elements.Figure 7
Design domain of the 3D L-shaped structure.The optimal topologies obtained by deterministic, Monte Carlo simulation and the robust approaches with density value above 0.3 thresholds are presented in Figures8 and 9, respectively. The corresponding optimization results are summarized in Table 2, including the required number of points, statistical moments, volume constraint, iterations, and computing time. Meanwhile, full tensor product grid and sparse grid are investigated under level l
=
1,2, and 3.Table 2
Comparison of robust and deterministic topology optimization design results.
Tests
Number of points
Meanμ
f
Standard deviationσ
f
V
/
V
∗
Iterations
Time
Deterministic
1
241.08
131.42
0.06
44
7 min 14 s
Monte Carlo simulation
1000
182.20
111.34
0.06
51
28 min 27 s
Robust
Tensor product formula (l
=
1)
9
183.95
121.71
0.06
46
8 min 19 s
Tensor product formula (l
=
2)
25
182.62
119.30
0.06
49
9 min 44 s
Tensor product formula (l
=
3)
81
182.91
113.72
0.06
49
17 min 26 s
Sparse grid method (l
=
1)
5
183.20
122.66
0.06
47
8 min 18 s
Sparse grid method (l
=
2)
13
182.71
119.24
0.06
50
8 min 47 s
Sparse grid method (l
=
3)
29
182.93
114.02
0.06
50
10 min 25 sFigure 8
Optimal topologies of the L-shaped structure.
(a)
Deterministic design
(b)
Monte Carlo simulation designFigure 9
Optimal designs of the L-shaped structure by robust topology optimization.
(a)
Tensor product grid(
l
=
1
)
(b)
Sparse grid(
l
=
1
)
(c)
Tensor product grid(
l
=
2
)
(d)
Sparse grid(
l
=
2
)
(e)
Tensor product grid(
l
=
3
)
(f)
Sparse grid(
l
=
3
)Comparing the deterministic design with the proposed robust designs, obvious distinctions can be clearly observed in optimal topologies as shown in Figures8 and 9. From the viewpoint of geometry, the robust designs represent more rational material layouts, which lead to more robust structural configuration. According to Table 2, the final statistical moments of the robust designs are superior to that of deterministic design. This indicates that the robust design is less sensitive to the variations of load magnitude. It can also be seen from the results listed in Table 2 that as the level l is increasing, the mean and standard deviation of the structure compliance for tensor product grid and sparse grid are closer to those for Monte Carlo simulation in which more robust designs are obtained, but it is also revealed that the corresponding computing time is increasing. This means robustness and cost are conflicting criteria in robust topology optimization design. For comparison of the results of tensor product grid and sparse grid, these indicate that sparse grid, in a manner, performs better than the tensor product grid comprehensively considering the accuracy and efficiency.


## 5. Conclusion

This paper integrates the stochastic collocation methods into structural topology optimization for solving the RTO problem with random loading uncertainty. Specially, tensor product grid and sparse grid are employed to transform the RTO into a simplified weighted multiple loading deterministic problem at the collocation points. The proposed robust approaches are implemented, discussed, and tested on two numerical examples by implementation in existing commercial topology optimization software package. An extensive Monte Carlo simulation is also performed to verify the accuracy and validation of the proposed approach. From the comparison of the results, the following conclusions are summarized:(
1
) The proposed RTO approach possesses more robust topologies than DTO and may exhibit manifest distinct topology configuration. It is revealed that robust topology optimization is essential when involving the uncertainty parameters in optimization design. The stochastic collocation methods are compatible with RTO for calculating statistical moments of the objective performance response and are readily implemented in existing commercial topology optimization software package.(
2
) For tensor product grid and sparse grid, more robust results are clearly obtained under higher level setting. Meanwhile, the required computation cost is more expensive. This indicates that robustness enhancement is inevitably accompanied with the sacrifice of computational expense. Sparse grid demands fewer collocation points than tensor product grid to generate similar even more robust results under the same level. At the same time, the computation burden for tensor product grid can be overcome to a certain extent. It is demonstrated that sparse grid, in a manner, performs better than the tensor product grid in compromise between robustness and computation cost.(
3
) It is important to notice that, compared to deterministic topology optimization, the application of robust topology optimization to practical structural engineering is still a challenge task because of the high RTO computation cost. Future research is still recommended to implement the proposed RTO approach for optimization problems subject to non-Gaussian type random loading uncertainty.



---

*Source: 580980-2015-08-30.xml*