Nilton Cuellar
email: ncuellar@tecgraf.puc-rio.br

Anderson Pereira
email: anderson@tecgraf.puc-rio.br

Ivan Fabio

Mota De Menezes

Americo Cunha
email: americo@ime.uerj.br

Ivan F M Menezes

Non-intrusive polynomial chaos expansion for topology optimization using polygonal meshes

Keywords: topology optimization, stochastic spectral approach, polynomial chaos, Karhunen-Loève expansion, robust optimization, polygonal finite element

Introduction

Due to new requirements of design associated with the most modern engineering applications, mechanical systems with very complex geometrical configurations are becoming increasingly common. In this context, some of the most promising design approaches are based on topology optimization (TO), which seeks to find the best layout for a system, by optimizing the material distribution in a predefined design domain [START_REF] Michell | Lviii. the limits of economy of material in frame-structures[END_REF][START_REF] Bendsøe | Topology Optimization: Theory, Methods, and Applications[END_REF]. The growing popularity of TO solutions is demonstrated by their wide range of application in various fields such as structural mechanics [START_REF] Talischi | Polygonal finite elements for topology optimization: A unifying paradigm[END_REF][START_REF] Talischi | PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes[END_REF][START_REF] Romero | A topology optimization approach applied to laminar flow machine rotor design[END_REF][START_REF] Qizhi | A topology optimization method for geometrically nonlinear structures with meshless analysis and independent density field interpolation[END_REF][START_REF] Dapogny | Geometric constraints for shape and topology optimization in architectural design[END_REF][START_REF] Luo | A multi-material topology optimization approach for wrinkle-free design of cable-suspended membrane structures[END_REF], composite and multi-materials [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Zhang | Multi-material topology optimization with multiple volume constraints: a general approach applied to ground structures with material nonlinearity[END_REF], nanotecnology [START_REF] Nanthakumar | Surface effects on shape and topology optimization of nanostructures[END_REF], fluid mechanics [START_REF] Pereira | Fluid flow topology optimization in polytop: stability and computational implementation[END_REF][START_REF] Duan | Adaptive mesh method for topology optimization of fluid flow[END_REF], fluid-structure interaction [START_REF] Andreasen | Topology optimization of fluid-structure-interaction problems in poroelasticity[END_REF], medicine [START_REF] Park | Design of complex bone internal structure using topology optimization with perimeter control[END_REF], etc.

In general, the physical systems underlying TO applications are subjected to a series of uncertainties (e.g. unknown loads, geometrical imperfections, fluctuations in physical properties, etc) so that, usually, their response is not well predicted by the traditional (deterministic) tools of engineering analysis. For this reason, there is a consensus among computational engineering experts that uncertainties effects must be incorporated into any computational predictive model [START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamics -recent theoretical advances[END_REF]. The modeling and quantification of uncertainties is necessary in order to predict a possible range of variability for the mathematical model response, and to conduct applied tasks, such as analysis and design, in a robust way [START_REF] Banichuk | Structural Optimization with Uncertainties[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. Notice, however, that the majority of the works in TO area, currently available in the literature, are limited to deterministic analyses.

The need for robust design and analysis of uncertainties in topological optimization applications naturally induces the search for computationally efficient frameworks for TO. For this purpose, TO literature started to take uncertainty quantification (UQ) into account over the last decade, as can be seen in several papers addressing the two issues [START_REF] Zhang | Robust shape and topology optimization considering geometric uncertainties with stochastic level set perturbation[END_REF][START_REF] Keshavarzzadeh | Topology optimization under uncertainty via non-intrusive polynomial chaos expansion[END_REF][START_REF] Da Silva | Stress-based topology optimization of continuum structures under uncertainties[END_REF][START_REF] Zhang | Stochastic sampling for deterministic structural topology optimization with many load cases: Density-based and ground structure approaches[END_REF][START_REF] Putek | Shape and topology optimization of a permanent-magnet machine under uncertainties[END_REF][START_REF] Richardson | A unified stochastic framework for robust topology optimization of continuum and truss-like structures[END_REF][START_REF] Zhao | Robust structural topology optimization under random field loading uncertainty[END_REF][START_REF] Dunning | Robust topology optimization: Minimization of expected and variance of compliance[END_REF][START_REF] Jalalpour | Reliability-based topology optimization of trusses with stochastic stiffness[END_REF][START_REF] Tootkaboni | Topology optimization of continuum structures under uncertainty -A polynomial chaos approach[END_REF][START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF][START_REF] Chen | Level set based robust shape and topology optimization under random field uncertainties[END_REF][START_REF] Guest | Structural optimization under uncertain loads and nodal locations[END_REF][START_REF] Kim | Efficient shape optimization under uncertainty using polynomial chaos expansions and local sensitivities[END_REF][START_REF] Wu | Robust topology optimization for structures under interval uncertainty[END_REF].

Some of these works are based on classical techniques for stochastic computation like Monte Carlo (MC) method [START_REF] Chen | Level set based robust shape and topology optimization under random field uncertainties[END_REF][START_REF] Dunning | Robust topology optimization: Minimization of expected and variance of compliance[END_REF] or series expansion [START_REF] Da Silva | Stress-based topology optimization of continuum structures under uncertainties[END_REF][START_REF] Zhao | Robust structural topology optimization under random field loading uncertainty[END_REF][START_REF] Jalalpour | Reliability-based topology optimization of trusses with stochastic stiffness[END_REF][START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF][START_REF] Guest | Structural optimization under uncertain loads and nodal locations[END_REF] which, despite of being very simple in conceptual terms, they are limited by the high computational cost, the former, or very small range of applicability, the latter. These limitations open space for spectral-based approaches [START_REF] Kim | Efficient shape optimization under uncertainty using polynomial chaos expansions and local sensitivities[END_REF][START_REF] Tootkaboni | Topology optimization of continuum structures under uncertainty -A polynomial chaos approach[END_REF][START_REF] Richardson | A unified stochastic framework for robust topology optimization of continuum and truss-like structures[END_REF][START_REF] Zhang | Robust shape and topology optimization considering geometric uncertainties with stochastic level set perturbation[END_REF][START_REF] Putek | Shape and topology optimization of a permanent-magnet machine under uncertainties[END_REF][START_REF] Keshavarzzadeh | Topology optimization under uncertainty via non-intrusive polynomial chaos expansion[END_REF], that use state-of-the-art tools for representing and propagating uncertainties in computational models, like Karhunen-Loève (KL) and generalized polynomial chaos (gPC) expansions. A recent work by Keshavarzzadeh et. al [START_REF] Keshavarzzadeh | Topology optimization under uncertainty via non-intrusive polynomial chaos expansion[END_REF] presents a non-intrusive gPC strategy to propagate uncertainties in topology optimization problems. They use non-intrusive polynomial chaos expansion to evaluate loworder statistics of compliance and volume and the uncertainties are considered in the applied loads and also in the geometry of the problems.

The classical formulation for topology optimization, which corresponds to minimize the structural compliance, is commonly carried out on uniform grids consisting of Lagrangian-type finite elements (e.g., linear quads). However, this choice of discretization, together with density methods, suffer from the well-known numerical instabilities, such as the checkerboard patterns. Unstructured "Voronoi" meshes, generated from an initial set of random points, have been shown to be effective in suppressing checkerboard patterns [START_REF] Talischi | Polygonal finite elements for topology optimization: A unifying paradigm[END_REF]. Moreover, compared to standard Lagrangian-type uniform grids, polygonal elements are more versatile in discretizing complex domains and in reducing mesh dependency in the solutions of topology optimization [START_REF] Talischi | Polygonal finite elements for topology optimization: A unifying paradigm[END_REF][START_REF] Antonietti | On the virtual element method for topology optimization on polygonal meshes: A numerical study[END_REF]. The geometrical flexibility of the polygonal finite elements also make them very attractive for adaptive mesh refinement schemes in topology optimization problems [START_REF] Nguyen-Xuan | A polytree-based adaptive polygonal finite element method for topology optimization[END_REF][START_REF] Hoshina | A simple adaptive mesh refinement scheme for topology optimization using polygonal meshes[END_REF]. The computational code used here was developed based on PolyTop [START_REF] Talischi | PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes[END_REF], a MATLAB code for solving topology optimization problems using either structured or unstructured polygonal meshes in arbitrary two-dimensional domains. The modular structure of PolyTop, where the analysis routine and optimization algorithm are separated from the choice of topology optimization formulation, together with a non-intrusive way of computing the statistical measures, allowed us to implement a robust topology optimization code in a very straightforward way, with only a few modifications in the original PolyTop code.

The aim of this paper is to present a computationally efficient and accurate non-intrusive robust topology optimization approach using polygonal elements. For this purpose, the PolyTop framework by Talischi et al. [START_REF] Talischi | PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes[END_REF], which employs polygonal finite elements in TO, is combined with a consistent methodology for stochastic analysis that uses a non-intrusive gPC strategy to propagate parameters uncertainties through the computational model. This combination generates a framework that is computationally efficient for stochastic simulations, and robust to numerical instabilities typical of TO problems, such as checkerboards, one-node connections, and mesh dependency. The novel approach is used to solve TO problems that seek to minimize an objective function based on the low-order statistical moments of the compliance function of a structure, subjected to uncertainties on the external load, and satisfying volume constraints.

The remainder of this paper is organized as follows. Stochastic spectral methods are introduced in section 2, together with the mathematical formulation and basic steps to obtain KL and gPC expansions. In section 3, the TO problem is briefly described, as well as stochastic procedure to propagate uncertainties within TO, which is called robust topology optimization (RTO). In section 4, numerical examples are presented, and the proposed methodology is compared with a reference result obtained with MC method. Finally, some remarks and suggestions for future work are presented in section 5.

Stochastic spectral methods

Preliminary definitions and notation

Consider a probability space (Θ, F, P ), where Θ is the sample space, F a σ-field over Θ, and P : F → [0, 1] denotes the probability measure. It is assumed that the distribution P X (dx) of any real-valued random variable X in this probability space admits a density x → p X (x) with respect to dx. The set of values where this density is not zero is dubbed the support of X, being denoted by Supp X.

In this probabilistic setting, any realization of random variable X is denoted by X(θ) for θ ∈ Θ, and the mathematical expectation operator is defined by

E {X} = R x P X (dx), (1)
so that the mean value and standard deviation of X are given by µ X = E {X} and

σ X = (E X 2 -E {X} 2 ) 1/2 , respectively. The random variable X is said to be of second-order if E |X| 2 = R |x| 2 P X (dx) < +∞. ( 2
)
The space of all second-order random variables in (Θ, F, P ), denoted by L 2 (Θ, P ), is a Hilbert space [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] equipped with the inner product

•, • Θ : L 2 (Θ, P )×L 2 (Θ, P ) → R such that X, Y Θ = E {X Y } = R 2
x y P X,Y (dx, dy),

where P X,Y (dx, dy) is the joint distribution of the random variables X, Y ∈ L 2 (Θ, P ). This inner product induces a norm

• Θ : L 2 (Θ, P ) → R where X Θ = X, X Θ 1/2 = R |x| 2 P X (dx) 1/2 . ( 4
)
Further ahead it will also be helpful to consider L 2 (D), the set of all real-valued square integrable functions defined on the spatial domain D ⊂ R d , d ≥ 1. This set of functions is also a Hilbert space [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], with inner product

•, • D : L 2 (D)×L 2 (D) → R defined by φ, φ D = D φ(x) φ (x) dx, (5)
for φ, φ ∈ L 2 (D).

Karhunen-Loève expansion

The KL expansion [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF] is one of the most widely used and powerful techniques for analysis and synthesis of random fields, providing a denumerable representation, in terms of the spectral decomposition of the correlation function, for a random field parametrized by a nondenumerable index [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF][START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF]. Let the map (x, θ) ∈ D × Θ → U (x, θ) ∈ R be an arbitrary real-valued random field, indexed by the spatial coordinate vector x ∈ D ⊂ R d , d ≥ 1, denoted in an abbreviated way as U (x) or U . By construction, for a fixed x ∈ D, U (x, •) is a real-valued random variable, and U (•, θ), for a fixed θ ∈ Θ, is a function of x, dubbed realization of the random field.

The correlation of U (x) is the function K U (•, •) : D × D → R defined for any pair of vectors x and x by means of

K U (x, x ) = E U (x) U (x ) . (6)
Suppose that random field U (x) is second-ordered and mean-square continuous, properties respectively defined by

E |U (x)| 2 < +∞, ∀ x ∈ D, (7)
and lim

x →x U (x ) -U (x) 2 Θ = 0. ( 8
)
Under these assumptions, the linear integral operator

K U φ (x) = D K U (x, x ) φ(x ) dx (9)
defines a Hilbert-Schmidt operator [START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF], which has denumerable family of eigenpairs (λn, φn)

+∞ n=1 such that D K U (x, x ) φn(x ) dx = λn φn(x), x ∈ D, (10)
where λn are the eigenvalues and φn the corresponding eigenfunctions of the operator defined by Eq. [START_REF] Betz | Numerical methods for the discretization of random fields by means of the Karhunen-Loève expansion[END_REF]. Besides that, the sequence of eigenvalues is such that

+∞ n=1 λn < +∞ and λ 1 ≥ λ 2 ≥ • • • ≥ λn ≥ • • • → 0;
where Kronecker delta is such that δmn = 1 if m = n and δmn = 0 for m = n. Therefore, applying two standard results of functional analysis [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF], the theorems of Hilbertian basis and orthogonal projection, it is possible to show that the random field U (x) admits a decomposition

U (x) = µ U (x) + +∞ n=1 √ λn φn(x) ξn , (12)
where {ξn} +∞ n=1 is a family of random variables defined by

ξn = 1 √ λn U -µ U , φn D , (13)
which are centered (zero mean) and mutually uncorrelated i.e.

µ ξn = 0, and

E {ξm ξn} = δmn . (14)
A finite dimensional approximation for U (x), denoted by U ν kl (x), is constructed by truncation of the series in Eq. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] i.e.

U ν kl (x) = µ U (x) + ν kl n=1 √ λn φn(x) ξn, (15)
where the integer ν kl is chosen such that

energy(ν kl ) = ν kl n=1 λn +∞ n=1 λn ≥ τ, (16)
for a heuristically chosen threshold τ (e.g. τ = 90%). In practice, as a closed formula for λn is not available in general, energy(ν kl ) is estimated using a finite (but large) number of eigenvalues, instead of an infinite quantity. This procedure is justified in light of the eigenvalues decreasing property. One of the main difficulties to apply KL expansion to discrete random fields is the determination of the eigenvalues and corresponding eigenfunctions of the correlation function. Analytical solutions for the Fredholm integral equation in [START_REF] Beyer | Robust optimization -A comprehensive survey[END_REF] are almost never available. However, for some special cases, such as exponential and Gaussian autocovariance functions, an analytical solution can be obtained by converting the integral equation into a differential equation through successive derivatives [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF].

Several numerical methods can be used to solve the eigenvalue problem of Eq. [START_REF] Beyer | Robust optimization -A comprehensive survey[END_REF], such as the direct method, projection methods, among others [START_REF] Atkinson | The Numerical Solution of Integral Equations of the Second Kind[END_REF][START_REF] Betz | Numerical methods for the discretization of random fields by means of the Karhunen-Loève expansion[END_REF]. In this study, the direct method is employed to transform the Fredholm integral equation into a finite dimensional eigenvalue problem, whose the solution provides an approximation for the desired eigenvalues/eigenvectors of the infinite dimensional problem.

In this numerical procedure, a set of M realizations of the random field U (x) and its mean function µ U (x) are numerically generated1 and grouped into the matrices

U =       U 1 (x 1 ) U 2 (x 1 ) . . . U M (x 1 ) U 1 (x 2 ) U 2 (x 2 ) . . . U M (x 2 ) . . . . . . . . . . . . U 1 (xn) U 2 (xn) . . . U M (xn)       , µ =       µ U (x 1 ) µ U (x 1 ) . . . µ U (x 1 ) µ U (x 2 ) µ U (x 2 ) . . . µ U (x 2 ) . . . . . . . . . . . . µ U (xn) µ U (xn) . . . µ U (xn)       , (17)
which are used to define the zero mean matrix Û = Uµ. Then, the correlation matrix is estimated with the aid of

K Û = 1 M Û Û T , (18)
and the discrete eigenvalue problem

K Û Φ = ΛΦ, (19)
is solved to obtain the matrices Φ and Λ, which present approximations for the first M eigenfunctions/eigenvalues on the columns/main diagonal.

For further details on theoretical and practical aspects of KL expansion the reader is encouraged to see [START_REF] Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Bellizzi | Smooth decomposition of random fields[END_REF][START_REF] Perrin | Karhunen-Loève expansion revisited for vector-valued random fields: Scaling, errors and optimal basis[END_REF][START_REF] Pettersson | Polynomial Chaos Methods for Hyperbolic Partial Differential Equations: Numerical Techniques for Fluid Dynamics Problems in the Presence of Uncertainties[END_REF].

Generalized polynomial chaos expansion

The gPC expansion is a theoretical tool used to construct representations for random fields, with a denumerable or nondenumerable set of index, in terms of a denumerable collection of random variables weigthed by deterministic coefficients [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF][START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF]. It was introduced in the engineering community by R. Ghanem [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Spanos | Stochastic finite element expansion for random media[END_REF][START_REF] Ghanem | Ingredients for a general purpose stochastic finite elements implementation[END_REF] as a tool to compute approximate responses for problems involving random fields with unknown distribution and since the early 2000s, especially after the work of Xiu and Karniadakis [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for stochastic differential equations[END_REF], it has been used in many applications of computational stochastic mechanics, a trend that should increase [START_REF] Stefanou | The stochastic finite element method: Past, present and future[END_REF].

For the sake of theoretical development, consider a second-ordered random variable U : Θ → R which can be written in terms of a (possibly infinity) set of independent random variables ξ

1 (θ), ξ 2 (θ), ξ 3 (θ), • • • . Collecting these independent variables into the random vector ξ(θ) = (ξ 1 (θ), ξ 2 (θ), ξ 3 (θ), • • • ), dubbed the germ, it is possible to rewrite the original variable in the parametric form U = U (ξ).
In this context of second-ordered random variables, gPC expansion theory says that such random variable U admits a spectral representation

U (ξ) = +∞ n=0
un ψn(ξ), [START_REF] Dunning | Robust topology optimization: Minimization of expected and variance of compliance[END_REF] where the set of orthonormal polynomials {ψ 0 , ψ 1 , ψ 2 , • • • } is a basis for L 2 (Θ, P ), and the deterministic coefficients u 0 , u 1 , u 2 , • • • correspond to the coordinates of U in this infinite dimensional base. By definition ψ 0 (ξ) = 1, and due to the orthogonality property ψm, ψn Θ = δmn, one has

µ U = E {U } (21) = E {ψ 0 U } = +∞ n=0 un ψ 0 , ψn Θ = u 0 ,
i.e., the mean of U is equal to the coordinate u 0 of gPC expansion. Besides that, it is easy to see that E {ψn} = 0, n ≥ 1. Thus, it can also be shown that standard deviation of U can be written as

σ 2 U = E U 2 -E {U } 2 (22) = E        +∞ n=1 un ψn(ξ)   2      = +∞ n=1 u 2 n .
The gPC expansion has attracted the attention of many researchers due to its rapid convergence property and its capability to estimate statistical moments, which allows this method to efficiently reduce computational effort in highly nonlinear engineering design applications [START_REF] Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Pettersson | Polynomial Chaos Methods for Hyperbolic Partial Differential Equations: Numerical Techniques for Fluid Dynamics Problems in the Presence of Uncertainties[END_REF].

Of course, for any purpose of numerical computation, it is necessary to parametrize the random variable U with a finite number νrv of independent random variables, i.e., ξ = (ξ 1 , ξ 2 , • • • , ξν rv ), and restrict up to ppc the order of the polynomials in the basis, so that the series in Eq.( 20) is truncated with νpc +1 = (νrv +ppc)!/(νrv! ppc!) terms, giving rise to the approximation

U νpc (ξ) = νpc n=0
un ψn(ξ), [START_REF] Ghanem | Polynomial Chaos: Modeling, Estimation, and Approximation[END_REF] which is mean-square convergent to U when ppc, νrv → +∞.

In practice, νrv is obtained from the stochastic modeling and ppc is specified such that a compromise between accuracy and efficiency can be established. In this work, the gPC polynomial order is selected such that the low-order statistics estimations become invariant for increasing values of ppc.

This spectral expansion can be easily extended to a second-ordered random field U : D × Θ → R, considering this field as a family of random variables U (x, •) ∈ L 2 (Θ, P ), parameterized by the index x, and by letting the deterministic coefficients of the gPC expansion depend on x, i.e. U (x, ξ) = +∞ n=0 un(x) ψn(ξ), [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] where the coordinates are now a deterministic function of x. As in the case of a random variable, the truncation of the series in Eq.( 24) results in the approximation

U νpc (x, ξ) = νpc n=0 un(x) ψn(ξ). (25)
Basically, two approaches for implementing the gPC expansion are available in the literature, one called the intrusive, and another one dubbed non-intrusive. The intrusive approach is based on a stochastic version of Galerkin method [START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF], and consists in modifying the deterministic model in order to take into account the uncertainty propagation. The main disadvantage of this formalism is related to the level of difficulty in modifying the code associated with the deterministic model, particularly if commercial software is being used, due to code access restrictions [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF][START_REF] Ghanem | Polynomial Chaos: Modeling, Estimation, and Approximation[END_REF].

On the other hand, the non-intrusive approach requires no modification in the deterministic model and therefore, can be treated as a black box [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF][START_REF] Ghanem | Polynomial Chaos: Modeling, Estimation, and Approximation[END_REF]. In this case, a probabilistic collocation approach, based on a sparse grid method [START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF], can be used to estimate the coefficients (coordinates) of the expansion. In this work, a different non-intrusive approach based on linear regression is employed, where the random field (variable) of interest is evaluated in a finite set of νgq possible realizations of the germ ξ, and thus the coefficients of the expansion are obtained through u = (Ψ T Ψ ) -1 Ψ T U, the solution of the mean-square problem

Ψ u ≈ U, (26)
where

    ψ 0 (ξ 1 ) • • • ψν pc (ξ 1 ) . . . . . . . . . ψ 0 ξν gq • • • ψν pc ξν gq     Ψ       u 0 u 1 . . . uν pc       u ≈     U (x, ξ 1 ) . . . U (x, ξν gq )     . U (27)
For further information about the basic aspects of gPC expansion the reader is encouraged to see the references [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF][START_REF] Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Pettersson | Polynomial Chaos Methods for Hyperbolic Partial Differential Equations: Numerical Techniques for Fluid Dynamics Problems in the Presence of Uncertainties[END_REF][START_REF] Ghanem | Polynomial Chaos: Modeling, Estimation, and Approximation[END_REF]31], and for more advanced topics [START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF][START_REF] Soize | Computational aspects for constructing realizations of polynomial chaos in high dimension[END_REF][START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF][START_REF] Kundu | Probabilistic optimization of engineering system with prescribed target design in a reduced parameter space[END_REF].

3 Topology optimization

Classical topology optimization

The main objective of TO is to find the optimal distribution of materials, for every point x in a given design domain Ω ⊂ R d , d = 2 or 3, which maximize a certain performance measure subjected to a set of design constraints, i.e., to determine which regions in Ω should not present material (void regions), and obtain the final topology of the structure.

By convention, points where material exists are represented by a density value of 1, otherwise, the density value is 0. Note that, in this way, one has an integerprogramming problem, where the distribution of material is defined by the density map x ∈ Ω → ρ(x) ∈ {0, 1}, for

ρ(x) = 1 if x is structural member, 0 if x is void. ( 28
)
In this paper the optimization problem seeks to minimize an objective function defined by the continuum structure compliance, denoted here by c, subject to a constraint on the final volume of the structure and satisfying the equilibrium equations for a linear Hookean solid material. This formulation, which is equivalent to maximize the structural stiffness subject to the same constraints (see reference [START_REF] Bendsøe | Topology Optimization: Theory, Methods, and Applications[END_REF] for further details), can be stated as

min ρ c u(ρ), ρ = Ω 1 2 σ : dΩ, s.t. v(ρ) = Ω ρ(x) dΩ ≤ v S , (29)
where σ and represent the tensors of stress and strain, respectively, v S is a specified upper bound on the optimized structure volume, and the map x ∈ Ω → u(ρ(x)) ∈ R d is the continum structure displacement, parametrized by ρ and implicitly defined by the elasticity equations

∇ • σ(u) = 0, σ(u) = σ(u) T , (u) = 1 2 ∇u + ∇u T , σ(u) = C(ρ) : (u), (30)
and the boundary conditions

σ(u) • n = t in Γ N , u = 0 in Γ D , (31)
where Γ D is the partition of ∂Ω on which the displacements are prescribed, Γ N is the complimentary partition of ∂Ω on which tractions t are prescribed such that

Γ D ∪ Γ N = ∂Ω and Γ D ∩ Γ N = ∅ and C(ρ)
is the 4th order stiffness tensor that depends on the density function ρ. As posed, finding ρ and C(ρ) becomes a large integer programming problem, which can be impractical to solve. Thus, we recast ρ as a continuous scalar field, ρ(x) ∈ [0, 1]. In order to recover the binary nature of the problem, the SIMP [START_REF] Bendsøe | Topology Optimization: Theory, Methods, and Applications[END_REF] model is employed and the stiffness tensor can be expressed as

C(ρ) = ε + (1 -ε)ρ p C 0 , (32)
where p > 1 is the penalty parameter, C 0 is the elasticity tensor of the constituent material, and 0 < ε 1 is a positive parameter ensuring well-posedness of the governing equations.

In terms of computational implementation, the domain is splited into N elementary regions, i.e., Ω = ∪ N e=1 Ω e , and the finite element method is employed for the solution of the elasticity equations. Thus, the following finite dimensional version of the optimization problem is considered

min ρ C(ρ) = F T U(ρ), (33)
s.t. V (ρ) = N e=1 ρe |Ω e | -v S ≤ 0, (34)
with

K(ρ) U(ρ) = F, (35)
where In order to use gradient-based optimization algorithms to solve the above formulation, gradients of the objective function as well as the volume constraint function are needed. The sensitivity of the objective function C(ρ) with respect to the design variables ρ, is expressed in component form as

ρ = (ρ 1 , ρ 2 , • • • , ρ N ) is a discretized
∂C ∂ρe = -U T ∂K ∂ρe U. ( 36
)
The gradient of the volume constraint function V (ρ) with respect to the design variable ρe is given as

∂V ∂ρe = |Ω e |. (37)
During the solution of a TO problem it is very common to deal with numerical anomalies, such as checkerboards, which are traditionally treated through the use of higher-order elements or filters [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, meshdependencies and local minima[END_REF][START_REF] Bruggi | On the solution of the checkerboard problem in mixed-FEM topology optimization[END_REF]. However, Talischi et al. [START_REF] Talischi | PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes[END_REF] have shown that the use of the PolyTop framework, which employs polygonal finite elements, can naturally address the checkerboard problem. Besides that, this approach also allows flexibility in the optimization strategy to be used, being compatible with the classical approaches based on the optimality criteria (OC) [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] and the method of moving asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF].

An overview of the classical TO procedure, used in PolyTop framework to obtain an optimal design, is illustrated in Figure 1. The sensitivity analysis step described in this schematic is explained in details in section 3.5.

Robust design optimization

Robust optimization, also known as robust design optimization (RDO), is a mathematical procedure that simultaneously addresses optimization and robustness analysis, obtaining an optimal design that is less susceptible to variabilities (uncertainties) in the system parameters. In contrast to conventional optimization, that is deterministic, RDO considers parameters that are random so that it consists of a stochastic problem. The general overview of RDO is explained in Figure 2, which shows a computational model where the input is subjected to uncertainties - that can be in material or geometrical properties, loadings, boundary conditions etc. -and, therefore, the model response has a certain probability distribution. This distribution is used to compute some kind of statistical response of the system, which is conveniently used to update the model input, in order to reduce the output uncertainty.

In this context, it is essential to understand the mathematical definition of robustness, i.e., the choice of the robustness measure that is generally expressed by the combination of statistical properties of the objective function. Several definitions of measures of robustness have been proposed in literature [START_REF] Beyer | Robust optimization -A comprehensive survey[END_REF][START_REF] Birge | Introduction to Stochastic Programming[END_REF][START_REF] Doltsinis | Robust design of structures using optimization methods[END_REF][START_REF] Shin | Computing trade-offs in robust design: Perspectives of the mean squared error[END_REF] and the weighted sum of both the mean and the standard deviation of the objective function is often considered. The tradeoff between these two statistical measures gives rise to a final design that is less sensitive to parameters variations, i.e., a kind of robust design.

Robust topology optimization

In order to increase the optimal design robustness, the concept of robust optimization described in section 3.2 can be applied to TO. This possibility is addressed in this paper where variabilities in the external loading acting on the structure of interest are taken into account. Thus, the force vector and the compliance function become random objects, more precisely, a random vector F(θ) and a random variable C(ρ, θ), both defined on the probability space (Θ, F, P ).

For the sake of computational implementation, these objects are parametrized by a set of νrv suitable random variables that are lumped into the germ ξ(θ) = ξ 1 (θ), ξ 2 (θ), • • • , ξν rv (θ) so that force vector and compliance can be expressed as F(ξ) and C(ρ, ξ).

A straightforward measure of structural performance (robustness of the objective function) in RDO framework is given by the mean of the compliance µ C(ρ,ξ) . However, the final design may still be sensitive to the fluctuation due to external loading uncertainties and this may give rise to the need for a more robust design [START_REF] Dunning | Robust topology optimization: Minimization of expected and variance of compliance[END_REF]. Then, the standard deviation of the compliance σ C(ρ,ξ) is also introduced into the formulation of the structural performance measure C, so that it is a linear combination between the mean and standard deviation of the random compliance C(ρ, ξ). By combining these two statistics it is possible to improve the design by minimizing the variability of the structural performance, satisfying the volume constraint.

This procedure, called robust topology optimization (RTO), can be mathematically formulated as

min ρ C(ρ) = µ C(ρ,ξ) + w σ C(ρ,ξ) , (38)
s.t. V (ρ) = N e=1 ρe |Ω e | ≤ v S ,
which depends on the weight w ≥ 0 and on the random displacement map U(ρ, ξ), implicitly defined by the random equilibrium equation

K(ρ) U(ρ, ξ) = F(ξ). (39)
The RTO problem defined in [START_REF] Park | Design of complex bone internal structure using topology optimization with perimeter control[END_REF] can be solved by considering non-intrusive methods for stochastic computation. The basic idea of non-intrusive methods is to use a set of deterministic model evaluations to construct an approximation of the desired (random) output response. The deterministic evaluations are obtained for a finite set of realizations of parameter ξ with the aid of a deterministic solver (e.g. finite element code), that is used as a black box. Thus, non-intrusive methods offer a very simple way to propagate uncertainties in complex models, such as structural optimization, where only deterministic solvers are available. In this study the focus is on two non-intrusive techniques, namely, MC simulation [START_REF] Kroese | Stochastic finite elements of discretely parameterized random systems on domains with boundary uncertainty[END_REF][START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF] and gPC expansion [START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF][START_REF] Ghanem | Polynomial Chaos: Modeling, Estimation, and Approximation[END_REF].

Low-order statistics for compliance

Monte Carlo (MC) method is one of the simplest crude techniques for stochastic simulation and may be used to construct mean-square consistent and unbiased estimations (approximations) -see [START_REF] Kroese | Stochastic finite elements of discretely parameterized random systems on domains with boundary uncertainty[END_REF] for details -for µ C = µ C(ρ,ξ) and σ C = σ C(ρ,ξ) , respectively defined by

µ C = 1 νmc νmc n=1 Cn, (40)
and

σ C =   1 νmc -1 νmc n=1 (Cn -µ C ) 2   1/2 , ( 41
)
where Cn = C(ρ, ξn), ξn is the n-th realization of the germ ξ and νmc denotes the number of MC realizations.

Despite its simplicity, the slow convergence rate of MC method (∼ 1/ √ νmc)

usually makes it a very expensive stochastic solver in terms of computational cost, particularly for TO problems, where a large number of deterministic model resolutions needs to be obtained to achieve an adequate response characterization. For this reason, a gPC procedure for low-order statistics estimation is also considered in this work. Using the gPC approach, an spectral representation of the compliance function can be written as

C(ρ, ξ) = +∞ n=0
cn(ρ) ψn(ξ), [START_REF] Putek | Shape and topology optimization of a permanent-magnet machine under uncertainties[END_REF] in a way that, because of properties E {ψ 0 } = 1 and E {ψn} = 0, n ≥ 1, the mean value of C(ρ, ξ) writes as

µ C = E C(ρ, ξ) = c 0 (ρ), (43)
where an approximation for the PCE coefficient c 0 is obtained from the linear regression [START_REF] Hoshina | A simple adaptive mesh refinement scheme for topology optimization using polygonal meshes[END_REF]. This procedure induces a Gaussian quadrature estimation of µ C , defined by the estimator

µ C = νgq j=1 W j C j , (44)
where C j = C(ρ, ξ j ) corresponds to the evaluation of the compliance at the Gauss points and the quadrature weights W j = Ψ † 1j are given by the first line entries of Ψ † = (Ψ T Ψ ) -1 Ψ T , the pseudoinverse of the νgq × νpc regression matrix Ψ .

By definition, the standard deviation of compliance is written as

σ C = E C(ρ, ξ) 2 -E C(ρ, ξ) 2 1/2 , (45)
so that a procedure similar to that used to construct the estimator of Eq.( 44) can be adopted now to propose

σ C =   νgq j=1 W j C 2 j -µ C 2   1/2 , (46)
as an estimator for σ C . The gPC-based estimators defined by Eqs.( 44) and ( 46) provide a very accurate and efficient framework for estimation of the compliance low-order statistics, that demands a small number of deterministic model evaluations, once in practice νgq νmc.

Sensitivity analysis

The partial derivative of the objective function C, defined in the optimization problem [START_REF] Park | Design of complex bone internal structure using topology optimization with perimeter control[END_REF], with respect to the element density function ρe is given by

∂ C ∂ρe = ∂µ C ∂ρe + w ∂σ C ∂ρe , (47)
where the partial derivatives on the right hand side can be approximated, via crude MC, with the aid of the estimators

∂µ C ∂ρe ≈ ∂µ C ∂ρe = 1 νmc νmc n=1 ∂Cn ∂ρe (48)
and

∂σ C ∂ρe ≈ ∂σ C ∂ρe = 1 (νmc -1) σ C      νmc n=1 Cn ∂Cn ∂ρe   -νmc µ C ∂µ C ∂ρe    . (49)
However, from the point of view of computational cost, it is more efficient to obtain these sensitivity coefficients using the gPC estimators, i.e.

∂µ C ∂ρe ≈ ∂µ C ∂ρe = νgq j=1 W j ∂C j ∂ρe (50)
and

∂σ C ∂ρe ≈ ∂σ C ∂ρe = 1 σ C      νgq j=1 W j C j ∂C j ∂ρe   -µ C ∂µ C ∂ρe    , (51)
obtained from Eqs.( 44) and ( 46) by differentiation with respect to ρe.

Algorithm for robust topology optimization

The results obtained from the TO algorithm, i.e., the compliance and sensitivities, are used to compute statistical measures in a non-intrusive way. Therefore, the RTO algorithm, for problems with uncertain loading, can be described as follows:

1. Topology Optimization Data: define finite element model, set optimizer and underlying numerical and control parameters; 2. Stochastic Model: parametrize aleatory objects with a set of independent random variables defined by the germ ξ. Choose an appropriate family of orthogonal polynomials, define weight factor w and gPC order ppc;

Objective Function

for each germ realization ξn perform finite element analysis using Eq.( 39) and compute compliance sensitivities with Eq.( 36); -Compute gPC coefficients from Eq.( 26); -Compute statistical estimates for mean and standard deviation with aid of Eqs.( 44) and ( 46); -Compute the sensitivity of the objective function from Eq.( 47). 4. Constraint Function: compute volume constraint using Eq.( 34) and its sensitivity with Eq.( 37); 5. Update the design variables ρ according to the optimizer. Repeat from step 3 until convergence is achieved;

A flowchart of the proposed RTO algorithm is depicted in Figure 3.

Numerical examples

The effectiveness of the proposed gPC RTO algorithm is addressed in this section by means of a study that considers bidimensional mechanical systems subjected to uncertain loads. The goal is to show that different statistical responses can be obtained when using the proposed gPC RTO design algorithm and a non-robust design strategy, where TO is done first (deterministically) and the propagation of uncertainties is computed later, considering the deterministic optimized topology.

For the sake of accuracy verification, a reference crude Monte Carlo RTO solution is employed. This comparison allows one the verify the accuracy of the statistical measures obtained with the proposed gPC approach. The influence of the weight factor in the robust design is also addressed, as well as the different effects that are observed when a random load is treated as a random variable or a random field.

For the examples presented in the following section, consistent units are used.

Cantilever beam design

This first example consists of a simple cantilever beam subjected to a pair of vertical loads, with uncertain magnitudes, applied at the two right edge corners, as illustrated in Figure 4(a). This problem has been studied by Wu et al. [START_REF] Wu | Robust topology optimization for structures under interval uncertainty[END_REF].

The vertical and horizontal dimensions of the structure are 30 and 60 units of length, respectively. The structure is composed of an isotropic material with Young modulus E 0 = 1 and Poisson ratio ν = 0.3. For the void material an elastic modulus value equal to E min = 10 -9 is employed. The prescribed volume fraction of material is set as 0.3, the filter radius as 1.5, the penalization factor 3, and the design domain is discretized by means of a polygonal mesh with N = 7, 200 finite elements. The nominal (deterministic) configuration for this problem adopts the magnitude of the two vertical forces as

F 1 = F 2 = 1, respectively.
On the other hand, in the stochastic case, magnitudes of the forces are assumed to be uncertain and modeled by independent random variables θ ∈ Θ → F 1 (θ) ∈ R and θ ∈ Θ → F 2 (θ) ∈ R, both defined on a suitable probability space (Θ, F, P ). For the sake of simplicity, but being consistent with the physics of the mechanical problem, it is assumed that these two random variables are uniformly distributed on the same positive support Supp F = [F min , Fmax] ⊂ (0, +∞). Three numerical studies are conducted in this example, where Supp The non-robust TO, obtained using the PolyTop with MMA optimizer, is shown in Figure 4(b). The lack of material on the left side of the cantilever is due to the two forces of equal magnitudes applied in opposite directions. Therefore, the stress in the cantilever is distributed only on the right side of the domain. However, to avoid instability (displacements going to infinity), a minimum value of elastic modulus E min is used.

In order to perform the gPC RTO one needs to define the germ ξ = (F 1 , F 2 ), which is over the region [F min , Fmax] × [F min , Fmax] ⊂ (0, +∞) × (0, +∞), so that νrv = 2. In this case the optimal base for gPC expansion is given by the Legendre polynomials [START_REF] Xiu | Numerical Methods for Stochastic Computations: A Spectral Method Approach[END_REF]. For the three different types of uniform distribution considered, a weight factor w = 1 is employed, together with an expansion of order ppc = 5 (so that 1 + νpc = 21)(this value was heuristically chosen to ensure the stochastic convergence.). A total number of νgq = 36 collocation points is used to generate realizations of ξ = (F 1 , F 2 ). To check the accuracy of the gPC RTO strategy, the same problem is addressed using the MC simulation with νmc = 10 4 scenarios of loading magnitudes, a reference result dubbed MC RTO. In Figure 5 the reader can see the optimum topologies obtained by gPC RTO (top) and MC RTO (bottom), for different support of the random variable F . The topologies shown in Figure 5 are different from the deterministic counterparts in Figure 4(b) -some extra members can be observed on the left side of the structure -for different levels of uncertainties (length of Supp F ). As the level of uncertainty increases, more members appear in the final topology. Finally, the robust designs using MC simulation present equivalent topologies and statistical measures, which demonstrates the accuracy of the proposed gPC RTO approach.

Table 1 compares statistical estimates for the compliance expected value µ C and standard deviation σ C , in the cases of robust and non-robust design. Remember that, in this context, non-robust design means first optimizing the topology via classical (deterministic) TO and then using MC simulation to propagate the loading uncertainties through the mechanical system. A good agreement between robust strategies based on gPC and MC is noted, as well as that the statistical measures for the non-robust design tend to approach infinity, since there is no connection between the left and right sides of the domain. It is also worth noting that, while the MC RTO needs 10 4 evaluations of the compliance function, the gPC RTO only needs 36 evaluations. This difference of three orders of magnitude demonstrates the efficiency of the gPC RTO implementation.

Table 1 Low-order statistics of cantilever beam compliance for robust and non-robust TO strategies.

Supp F gPC RTO MC RTO Non-robust TO µ C σ C µ C σ C µ C σ C [0.

Michell type structure

In this second example RTO is performed on a simple Michell type structure considering three loads, with uncertain directions, applied at the bottom edge of the two dimensional system, as illustrated in Figure 6(a). The length and height of the structure are equal to 120 and 50 units, respectively. The design domain is discretized with a polygonal mesh with N = 12, 000 finite elements, and all other parameters are the same as in the first example. For the deterministic case, magnitudes and directions of the three forces are defined as

F 1 = 1, F 2 = 2, F 3 = 1, and α 1 = α 2 = α 3 = -90 o , respectively.
Meanwhile, on the stochastic case, the forces directions are modeled as the independent and identically distributed random

variables θ ∈ Θ → A 1 (θ) ∈ R, θ ∈ Θ → A 2 (θ) ∈ R and θ ∈ Θ → A 3 (θ) ∈ R.
Three scenarios of probabilistic distribution are analyzed: (i) Normal, (ii) Uniform, and (iii) Gumbel. For the Normal and Gumbel distributions, mean values are assumed to be equal to the nominal values of α 1 , α 2 and α 3 , with all the standard deviations equal to 10 o . In the Uniform Now the germ is defined as ξ = (A 1 , A 2 , A 3 ), thus νrv = 3, and the family of orthogonal polynomials (basis) is chosen according to the germ support. For simplicity, Hermite polynomials are used in the case of Gaussian or Gumbel parameters, while Legendre polynomials are the option when the germ is uniform distributed. Employing gPC RTO with an expansion of order ppc = 5 (thus 1 + νpc = 56) total number of νgq = 216 collocation points and and weight factor value w = 1, one obtains the robust designs shown in Figure 8, where connections at fixed points are created to balance the horizontal components of non-vertical forces. Note that the non-robust design in Figure 6(b) only presents four bars connected at the forces application points, no connection at the joints can be seen. This occurs because the forces are always vertical. However, when gPC RTO design is used, there are connections at the joints, because the angle variability induces horizontal force components.

The statistical results of the robust design compared with the non-robust counterpart can be appreciated in Figure 9 and Table 2, which show the compliance probability densities and their low order statistics, respectively, for the different distributions considered in the force angle. One can observe from these simulation results that the range variability of compliance is reduced, which implies that the robust design is less sensitive to loading uncertainties than its non-robust counterpart. As shown in Figure 8, the final topologies are symmetric for the normal and uniform distributions but is asymmetric for the Gumbel distribution. In the second analysis for this example, the influence of the weighting factor in the robust optimization process is addressed. The aim is to minimize the variability 3, all the designs shown in Figure 10 present significant lower expected compliance and standard deviation values when compared to the non-robust solution. The highest values of standard deviation are obtained using w = 0, since we are minimizing only the expected compliance. For w > 0, both expected compliance and its standard deviation contribute to the objective function and we can observe that as the value of w increases, the standard deviation decreases. Based on the numerical experiments presented in Table 3, we recommend the value w = 3, for practical use, because it leads to the best values of standard deviation with only a slight change in the expected compliance values.

2D bridge structure

This last example corresponds to a simple bridge structure subjected to an uncertain distributed loading at the top edge, as illustrated in Figure 11(a). The nominal load is uniform throughout the structure, with magnitude per unit of length equal to F = 1.

For the stochastic case, the load magnitude per unit of length in each point is assumed to be a Gaussian random field (θ,

x) ∈ Θ × [0, l] → F (θ, x) ∈ R with correlation function K F (x, x ) = σ F , (52)
such that the loads at any pair of points x, x are fully correlated. The mean and standard deviation of the random field F (θ, x) are assumed as µ F = 1 and σ F = 0.3, respectively. For the optimization process, the prescribed volume fraction of material is set as 0.3, the filter radius is set as 3, the penalization factor 3, and the design domain is discretized with a polygonal mesh with N = 10, 000 finite elements. The nonrobust TO of the 2D bridge, performed using PolyTop with MMA optimizer, is shown in Figure 11(b). Furthermore, the first two rows of finite elements on the top of structure are fixed during the optimization process, to ensure that the bridge remains attached to the loading conditions. Allowing the final results to be more realistic. It is observed in Figure 11(b) that the non-robust design leads to a final topology which is similar to the classical case of a 2D bridge under an uniformly distributed load.

For the purpose of numerical computation, the random field F (θ, x) is discretized by means of ξ = (F 1 ), a single Gaussian random variable for which loworder statistics are the same as for the random field. The gPC RTO design is obtained using an expansion of order ppc = 5 (with 1 + νpc = 6) and a total number of νgq = 6 collocation points for Hermite orthogonal polynomials, and weight factor values w ∈ {0, 1, 2, 3}. The final results are shown in Figure 12, where one can observe that some bars connected at the bottom of the bridge are different from those of the non-robust case in Figure 11(b). Moreover, as the value of w increases, the structure presents a more robust physical form, which represents a consistent result, because the standard deviation of the compliance is being forced to be smaller. The corresponding mean value and standard deviation of the compliance function, for the different values of w employed, are given in Table 4. As a second analysis, the random load F (θ, x) is assumed to have the same low-order statistics as before, but an exponentially decaying correlation function where lcorr is a correlation length for the random field. Note that, by this assumption, the loads at any two points x, x in the 2D bridge are partially correlated.

K F (x, x ) = σ F exp - |x -x | lcorr , (53)
If the correlation length is increased, a strong correlation is obtained between the points x, x , so that lcorr = ∞ implies a perfectly correlated random field -the previous case where the field depends on a single random variable. On the other hand, when lcorr = 0, the random field is completely uncorrelatedmany independent random variables are necessary for an accurate computational representation. In order to avoid the two limit cases, lcorr = 120 is chosen.

In terms of computational representation for numerical calculations, the random field F (θ, x) is discretized with the aid of Karhunen-Loève expansion described in section 2.2. The number of terms in this expansion is chosen in a heuristic way, seeking to satisfy the criterion presented in [START_REF] Da Silva | Stress-based topology optimization of continuum structures under uncertainties[END_REF]. A good compromise between accuracy and computational efficiency is obtained with ν kl = 7. Therefore, the germ is ξ = (F 1 , F 2 , • • • , Fν rv ), a set of νrv = 7 independent Gaussian random variables for which low-order statistics are the same as for the random field. Then, for this case we use a total number of νgq = 279936 collocation points.

A comparison between non-robust and gPC RTO design, for w = 1 and the different types of distributed load considered, are shown in Figure 13. The difference between the three obtained topologies is very clear, and can also be appreciated in Table 5, which shows the low-order statistics of the compliance in all cases analyzed. This example clearly shows that the nature of the distributed load has a significant effect on the RTO. Also from Table 5, it is possible to see that the compliance low-order statistics for a 2D bridge under a distributed load, emulated by a partially correlated random field, are smaller than those for a fully correlated field.

Conclusions

In the present paper RTO problem has been formulated and solved by means of an optimization procedure which integrates a classical TO algorithm with a stochastic spectral expansion based on gPC. Monte Carlo simulation is used to verify the accuracy and efficiency of the proposed methodology. This approach is introduced to reduce the variability due to uncertain loadings applied to the mechanical structure of interest. The objective function of the robust problem is defined as the weighted sum of the mean and standard deviation of the compliance, and it can be computed by considering a number of additional load cases. This makes the RTO computationally tractable and accessible by any TO algorithm. Furthermore, the gPC is compatible with RTO for computing the statistical measures of the compliance. The numerical examples presented here show a substantial benefit and exhibits topology changes within their design domains compared with their deterministic counterpart. The optimal topology configurations confirm that the uncertainty parameters might change the deterministically obtained optimal topologies. The proposed methodology allows to obtain approximate outcomes with a much lower computational cost than that associated with Monte Carlo simulation, which makes it attractive, particularly in the context of structural topology optimization. Moreover, when using random load fields, the results show different topologies because the forces are correlated, i.e., each force depends on the other and therefore, their interactions with the structure have significant effects on the robust design. The limitation of the gPC can be observed when a large number of random variables is used to parametrize the stochastic model, since in this case a substantial number of terms is necessary to construct the expansion, and, consequently, the computational cost increases significantly with the dimension. This is often referred to as the curse of dimensionality, and it can be reduced using adaptive techniques such as the adaptive sparse grid.

Figure 2

 2 Figure 2 General overview of the robust design optimization procedure.

Figure 3

 3 Figure 3 General flow chart of the gPC RTO integrated procedure.

60 F 1 F 2 30Figure 4

 60124 Figure 4 Cantilever beam structure: (a) original configuration, (b) non-robust TO design.

F

  is chosen as [F min , Fmax] = [0.95, 1.05], [F min , Fmax] = [0.9, 1.1] and [F min , Fmax] = [0.8, 1.2]. Note that these intervals correspond to symmetrical variabilities of up to 5%, 10% and 20% around the mean values µ F1 = µ F2 = 1, respectively.

  (a) gPC RTO design (b) MC RTO design

Figure 5

 5 Figure 5 Optimized topologies for the cantilever beam using gPC RTO and MC RTO designs, for uniform distributions over the intervals [F min , Fmax] = [0.95, 1.05] (left), [F min , Fmax] = [0.9, 1.1] (center) and [F min , Fmax] = [0.8, 1.2] (right).

Figure 6

 6 Figure 6 Michell type structure: (a) original configuration, (b) non-robust TO design.

Figure 7

 7 Figure 7 Probability distributions for the loads angles: normal (left), uniform (center) and Gumbel (right).

Figure 8

 8 Figure 8 Optimized topologies for the Michell type structure using gPC RTO design, with different probability distributions for load angle: (a) Normal, (b) Uniform and (c) Gumbel.

Figure 9 2 (d) w = 3 Figure 10

 92310 Figure 9 Probability density function of the compliance for the Michell type structure using non-robust and gPC RTO robust design, with different probability distributions for load angle: (a) Normal, (b) Uniform and (c) Gumbel.

Figure 11 2D

 11 Figure 11 2D Bridge structure: (a) original configuration, (b) non-robust RO design.

Figure 12

 12 Figure 12 Robust design for the 2D bridge structure with fully correlated distributed load, for different values of weight w.

Figure 13

 13 Figure 13 Optimized topologies for the 2D bridge structure: (a) non-robust TO with fully correlated load, (b) gPC RTO with fully correlated load, (c) gPC RTO with partially correlated load.

  and the family of functions

	{φn} +∞ n=1 defines an orthonormal Hilbertian basis in L 2 (D), i.e.
	φm, φn D = δmn ,

  version of the density map, |Ω e | denotes the volume of the element e, U is the discrete displacement vector, parametrized by ρ and implicitly defined by the equilibrium equation, F is the global loading vector and K represents the global stiffness matrix, which is also dependent on ρ.

  Overview of the classical topology optimization procedure (adapted from[START_REF] Bendsøe | Topology Optimization: Theory, Methods, and Applications[END_REF]).

	Initialize
	Finite Element Method
	Sensitivity Analysis
	Optimization Algorithm
	Update Design Variables
	Converged ?	no
	yes
	Plot Results /
	Post-processing
	Stop
	Robust Design Optimization Figure 1 Computational Model (Input variable) Model parameters i	(Output variable) Response metric Robust Design	i
			Robustness

Table 2

 2 Low-order statistics of the compliance for the Michell type structure using robust and non-robust TO strategies.

	distribution	gPC RTO	Non-robust TO
		µ C	σ C	µ C	σ C
	Normal	251.6	6.0	314.2	113.3
	Uniform	253.3	5.7	262.7	33.3
	Gumbel	249.1	6.7	312.5	128.9

Table 3

 3 Low-order statistics of the compliance for the Michell type structure using gPC RTO, with different force angles distributions and different weight factors.

					Supp A
			[-95 o , -85 o ] [-100 o , -80 o ] [-110 o , -70 o ]
		w	µ C	σ C	µ C	σ C	µ C	σ C
	gPC RTO	0	241.8	5.2	256.0	10.6	249.0	10.2
	gPC RTO	1	247.2	4.3	253.3	5.7	249.2	6.4
	gPC RTO	2	253.0	3.1	249.7	4.0	250.0	6.0
	gPC RTO	3	248.9	2.4	247.0	2.9	250.1	5.6
	Non-robust TO	-	364.4	14.2	366.1	28.4	373.0	56.9

Table 4

 4 Low-order statistics of the compliance for a 2D bridge using gPC RTO and different weight factors.

	w	µ C	σ C
	0	5.4443 E5 3.0666 E5
	1	5.4481 E5 3.0657 E5
	2	5.4462 E5 3.0646 E5
	3	5.4455 E5 3.0642 E5

Table 5

 5 Low-order statistics of the compliance for a 2D bridge considering different design scenarios.

	Non-robust TO	RTO full corr	RTO partial corr
	µ C	σ C	µ C	σ C	µ C	σ C
	5.449 E5 3.076 E5 5.448 E5 3.066 E5 2.361 E5 9.708 E4

These numerical realizations are defined in a computational mesh x 1 , x

, • • • , xn.

Acknowledgments

NC acknowledges the financial support from the Group of Technology in Computer Graphics (Tecgraf/PUC-Rio), Rio de Janeiro, Brazil. AP and IFMM acknowledge the financial support from the National Council for Scientific and Technological Development (CNPq) under projects 312280/2015-7 and 309708/2015-0, respectively. AP and ACJr are thankful for the support from Carlos Chagas Filho Research Foundation of Rio de Janeiro State (FAPERJ) under grants E-26/203.189/2016, E-26/010.002.178/2015 and E-26/010.000.805/2018. The information provided in this paper is the sole opinion of the authors and does not necessarily reflect the views of the sponsoring agencies.


Camilo Salguero

Sandra P Mosquera

Andrés F Barco

Élise Vareilles

Low Cost Intelligent System for the 2D Biomechanical Analysis of Road Cyclists

Keywords: Augmented Reality, Biomechanical analysis, Intelligent Systems, Motion Capture, Road Bicycle, Software architecture

This paper introduces an intelligent system focused on the biomechanical analysis of road bicycle cyclists. This type of analysis is carried out in specialized medical centers that operate using costly resources and are employed, mainly, for studies on athletes of high performance. The proposed system contrasts with these centers in that it provides the rookie cyclist with an accessible and affordable biomechanical analysis, although not as accurate. The architecture of the system rests in the advances in motion capture and augmented reality libraries. In the paper are discussed the motivations of the research, the internal design of the proposed system and the differences with various systems.

Introduction

Biomechanics is defined as the study of the movement of living beings using the mechanics analysis, which involves forces as creators of movement [START_REF] Knudson | Fundamentals of Biomechanics[END_REF]. The discipline studies the capture, measurement and analysis of human body movement with the aim of avoiding injuries and improve sports performance. Our research is focused on the movement performed by cyclists in the bicycle route category [START_REF] Too | Biomechanics of cycling and factors affecting performance[END_REF][START_REF] Mota | Influence of bike fit in lower limb biomechanics in road cyclists[END_REF], one of the most important sports world-wide.

Frequently beginner cyclists, particularly those of bicycle route category, have a bad or uncomfortable posture when using their bicycles which generally has a negative impact on their performance [START_REF] Too | Biomechanics of cycling and factors affecting performance[END_REF]. This situation is due to unknowing of the correct positioning on the bicycle that may lead to an optimal performance. Normally, beginner cyclists come to good postures only through practice and imitation of positioning of professional or high performance athletes. This imitation of expert cyclists is not straightforward as differences on positioning cannot be observed with the naked eye and much less if the cyclist is riding.

Biomechanical studies with traditional systems carried out in specialized centers are usually done in controlled environments, mainly with controlled illumi-nation, with several cameras and specialized sensors and with non versatile proprietary software [START_REF] John | The economic impact of sports facilities, teams and mega-events[END_REF]. In addition to this, access to a biomechanical evaluation is impossible for many cyclists in developing countries as it is an expensive service. In Colombia, for instance, where road cycling is one of the national sports, it exists only two biomechanical centers (visit [14,[START_REF]Biomec: Bikefitting 3d[END_REF] for details in Spanish) that performs an evaluation that cost around 250 US Dollars (around 750.000 COP), which is not affordable for most beginners cyclists.

Our goal is to develop a versatile and portable biomechanical analysis intelligent systems that contrasts with specialized center systems in that it does not require expensive devices or complex features. Given the low cost characteristic, the system does not have the same precision in its bicycle and cyclist measurements and is not able to implement all the biomechanical analysis techniques that standard systems can perform. However, it does strengthen the versatility with which athletes can improve their positioning on their bicycles

The rest of the document presents the following topics. Related work of systems that currently use low resources are presented in Section 2. The design of an intelligent system for the biomechanical analysis is discussed in Section 3. Afterwards, the advantages and disadvantages of our proposal are shown in Section 4. Finally, conclusions and future work are drawn in Section 5.

Related Work

In this section we present related work of systems for the biomechanical analysis that use of low resources but that are not fully automatic systems, meaning that the human intervention is required, given opportunity to configuration failures and complication.

The first related work is called Kinematic analysis of cyclist pedaling using systems of motion capture in 2D and 3D [START_REF] Del Blanco | Análisis cinemático del pedaleo ciclista mediante sistemas de captura de movimiento en 2d y 3d[END_REF]. This research focused on checking the reliability of the data thrown by a low cost analysis system that works in 2D and compares it with a system that works in 3D that makes use of different expensive resources. The study was carried out through the simultaneous use of both systems (2D and 3D) over a group of 12 cyclists, and the subsequent comparison of the data offered by both systems. Among the results and conclusions, authors did not find important differences between the sensitivity of both systems. On the other hand, as a drawback in the 2D measurement systems, the data obtained for the evaluation was not captured in real time, a characteristic of the 3D systems. The intelligent system that we propose works as well in 2D but it does operates real time, which addresses the disadvantage of the 2D system.

A second work on the subject is found in [START_REF] Márquez | Análisis biomecánico para ciclistas[END_REF] and it is called Biomechanical analysis for cyclists. The research focuses on the analysis of the position of the lower joints of the rider. The authors developed an application for Android systems using the Java programming language and the OpenCV library [START_REF] Pulli | Real-time computer vision with opencv[END_REF]. To make use of their system, the cyclist must locate LED bulbs in his lower joints and record a video of the cyclist performing the pedaling movement on the bicycle. Finally, the system presents graphs with the data of the angles of the cyclist's joints movement. In contrast with this system, our system does not force the uses of LED bulbs but whatever color marker that facilitates the motion capture process. Additionally, the capture is done in real time without the need of an intermediate video.

Another work on biomechanical analysis is found under the name Mobile tracking system and optical tracking integration for mobile mixed reality [START_REF] Bajana | Mobile tracking system and optical tracking integration for mobile mixed reality[END_REF]. This work aims at developing a system that combines optical tracking and augmented reality, operating in a local area network as it communicates with a server responsible for processing all tracking data. The idea behind it is to reduce the responsibilities of the mobile device that have limited processing and storage capabilities. The system is possible thanks to the implementation of a fast communication channel between mobile device and servers.

Apart from the described works, there are commercial systems that allow cyclists to perform a biomechanical analysis of low cost. To the best of the authors knowledge, these applications do not have related academic publications in which the underlying models and techniques are described.

Bike Fast Fit [START_REF]Bike fast fit[END_REF] is an application available for iOS devices, which allows to capture in video the pedaling of the cyclist to subsequently measure the angles presented in the recording. This measurement is done to a certain degree manually, meaning that the user has to choose the frames of the video where it is potentially correct to take the measurement. Also the user must choose the body angle to measure, and then locate the measurement points in the cyclist's joints. The advantage of this system is that it has a good precision in the measurements. Nonetheless, the precision may be affected by the intervention of the user as it can make mistakes when measuring and configuring the system.

The Roadie Bike Fit [START_REF]The roadie bike fit[END_REF] is an application available for Android devices, which offers the user a guide and recommendations for the adjustment of the bicycle. The guide begins with the capture of two images, one of them of the bicycle and another of the cyclist. With these two images the application executed an object recognition algorithm with which it identifies body angles and bicycle measurements. The automation and guidance feature make of this a simple and intuitive application for the user. Nevertheless, the measurements are obtained from static positions of the cyclist, which does not give reliable data of the movements that can only be measured while the cyclist pedals.

3 System's Underlying Design

Solution Scheme

In the face of the problem, the proposed system is able to offer the route cyclist a kinematic analysis and recommendations to adapt the bicycle and body to their correct positions. This is done through the measurement of the body angles involved in the movements. The measurements are performed mostly in the sagittal plane (refer to [START_REF] Ferrer Roca | Comparación de diferentes métodos de ajuste de la bicicleta en ciclistas entrenados. influencia de factores biomecánicos y energéticos[END_REF] for information about different positions to capture movement data). Thus, our system uses a camera perpendicular to the direction of the bicycle, focusing on the cyclist's right profile, as shown in Fig. 1.a. Having this configuration, the system captures the data necessary to execute the biomechanical analysis, i.e., the angles formed by the body parts. As presented in Fig. 1.b, in order to measure the angles, the cyclist must locate markers of any color that are distinctive in the image (high contrast to improve image recognition). The system interaction is as follows. Once the capturing device (e.g., smart-phone camera) is in front of the cyclist, it starts marker recognition on the video stream. Using the frames of the video stream it computes angles (as we will see further in Section 3.2). Then the angles are rendered over a display and finally a recommendation based on an analysis is given.

Architecture

The intelligent system implements the architecture presented in Fig. 2. Essentially, the application executes the biomechanical analysis on a device (such as a tablet, mobile device or laptop) and statistics are stored on a server.

Client's side. The client's components are responsible for the most important tasks of the system, the core of the system; motion capture, data processing (biomechanical analysis) and visualization of results.

-Capture: This component captures the movement from the right profile of the cyclist. It uses only one video capturing device as a smart-phone camera and no sensors, making it low cost. The video stream is presented in a first HTML5 canvas. -Processing: The processing of the captured images are performed on the client's side, making use of the library tracking.js [START_REF] Community | Tracking.js[END_REF]. grates functionalities of image recognition and tracking of colors thus allowing to capture in real-time the position of the color markers that the cyclist is wearing. Having captured the positions of the markers, which are points in the vector space R 2 , the system proceeds with the computation of angles. To do so two vectors are need to be computed (as presented in Fig. 3) by using three pixels coordinates (from the markers) and applying basic vector geometry equations. Once the vectors # » V 1 and # » V 2 are computed, the system computes the vectors' magnitudes (Equations 1 and 2), computes the dot product using its algebraic notation (Equation 3) and finally it computes the angle using the dot product geometry notation (isolating θ in Equation 4).

Pixeles X

| # » V 1| = V 1x 2 + V 1y 2 (1)
| # » V 2| = V 2x 2 + V 2y 2 (2) # » V 1 • # » V 2 = V 1x × V 2x + V 1y × V 2y (3) # » V 1 • # » V 2 = | # » V 1| × | # » V 2| × cos θ (4)
Finally, the biomechanical analysis of the cyclist is done by using cyclist and bicycle measures as well as the computed angles. This analysis is under implementation in the system as it requires the knowledge from experts in the field to be extracted. -Visualization: The visualization is done in a fully augmented reality fashion. This technology integrates information captured from the real world with digital information, matching the two of them to build a new coherent complemented and enriched world [START_REF] Van Krevelen | A survey of augmented reality technologies, applications and limitations[END_REF]. The system does so by rendering information related to the positioning of the body and presenting recommendations of modification tailored to the bicycle. The data resulting from the analysis together with the bicycle adjustment recommendations are rendered in real time in a second HTML5 canvas.

Server's side. In the server' side are implemented the persistent storage for statistics on the systems' use.

-Storage: The server is in charge of storing information about the biomechanical analysis made by cyclists. So far, the system does not manage users as it can be executed without sessions but it is planned to add such functionality. Nevertheless, it stores users data and results for future statistical analysis. -Web Server: The server processes the collected data from users and analysis results in order to obtain use statistics. These statistics refer to the most common positioning errors made by novice cyclist, the number of times a given cyclist uses the tool, measurements of bicycles and execution time. This data is planned to be exploited with data mining algorithms in order to improve recommendations.

Key Differences with Other Systems

There exist some clear differences with respect to the systems in specialized centers and also with commercial systems available for mobile devices. We expose these as advantages and disadvantages.

1. One advantage w.r.t. systems developed as part of research or as a commercial product, is that out system requires very little configuration work on the part of the user. This is because the motion capture of the joints is highly automated, requiring only to select the appropriated color for the markers that the cyclist is wearing.

2. By having augmented reality as a visualization technology, the information related to the analysis is clearly displayed and easily understandable by the cyclist. This contrasts with the results in specialized centers that must be interpreted by experts. 3. The intelligent system offers the cyclist graphic recommendations, with the use of augmented reality, so that it can autonomously adjust the bicycle and thus implement a proper posture. No mechanics or expert intervention. 4. The system is implemented using a library that works in web environments, providing flexibility and a multi-platform development. Other systems either have too complex environment or are tied with a given platform. 5. The system uses only a video capture device and a processing unit (e.g. a smart-phone). Ergo, it is low cost unlike specialized centers.

Our system, on the other hand, presents the following disadvantages.

1. First, the system is affected in its measurements by the processing speed because it is designed to perform real-time analysis on a mobile devices that have limited processing and storage capabilities. It can therefore be slower than traditional systems. 2. Secondly, the angle and position of the camera may influence the reliability of the captured data which may generate inaccurate measurements. 3. Since it is a system that identifies the joints from markers, these markers are recognized by the system using an specific color. This process, made by the library tracking.js, is error prune and depends mostly on the color composition of the image (cyclist, bicycle and surroundings). Not being able to identify the marker implies not joints movement capture. 4. Finally, adequate lighting can be a determinant factor when capturing data.

Unlike the specialized centers where there is a controlled environment, the lighting when using the proposed application may vary vastly in different executions making the analysis difficult or inaccurate.

Concluding Remarks

In this paper we have presented the motivations and internal design of an intelligent system focused on the biomechanical analysis of road bicycle cyclists. In particular, the paper briefly discussed the context of cycling, the intrinsic need of novice athletes for analysis of positioning and movement and related work found in the literature and Web. Our research is possible given the advances in mobile devices (cameras), their access, high-speed WI-FI networks and advances in libraries for motion capture and augmented reality. As such, the system is an important starting point for the novice cyclists as it allows to; do motion capture (thanks to the library tracking.js); show the cyclist movement over the display in real time (thanks to the augmented reality technology), do the biomechanical analysis in real time using the limited processing and storage capabilities of a mobile device, store usage statistics for further analysis and be used in a friendly and intuitive way.

Future work. We have identified three strategic directions for future work.

-First, increase the accuracy of the analysis with the use of sensors. This step is fundamental for improving the analysis and is possible thanks to the inclusion of sensors in the latest mobile devices (mainly smart-phones). -Second, it is important to dwell more on usage statistics. Here, it is planned to make a analysis of the stored data (how many times a user executed the system, measurements, etc.) in the aim of providing good recommendations for the novice cyclist. In short, it is planned to apply data mining. -Lastly, it is planned to involve a knowledge base that helps the positioning of cyclist either by applying expertise knowledge or case-based reasoning.

Fig. 1 .

 1 Fig. 1. a) Camera captures right profile (sagittal plane). b) Angles in the sagittal plane involved in the movement.

Fig. 2 .

 2 Fig. 2. System's architectural design.

Fig. 3 .

 3 Fig. 3. Obtaining vectors from point in vector space R 2 .

  The library inte-

		Server Side
	Web Server (Apache)	User Analysis Data Usage Statistics	Relational Database (SQL)
	Client Side
	Recognition, tracking and		Capture
	processing (Tracking.js)	First canvas	(Monocular)
			Visualization
	First and second
	canvas