Robust topology optimisation of lattice structures with spatially correlated uncertainties

The uncertainties in material and other properties of structures are usually spatially correlated. We introduce an efficient technique for representing and processing spatially correlated random fields in robust topology optimisation of lattice structures. Robust optimisation considers the statistics of the structural response to obtain a design whose performance is less sensitive to the specific realisation of the random field. We represent Gaussian random fields on lattices by leveraging the established link between random fields and stochastic partial differential equations (SPDEs). It is known that the precision matrix, i.e. the inverse of the covariance matrix, of a random field with Mat\'ern covariance is equal to the finite element stiffness matrix of a possibly fractional PDE with a second-order elliptic operator. We consider the discretisation of the PDE on the lattice to obtain a random field which, by design, considers its geometry and connectivity. The so-obtained random field can be interpreted as a physics-informed prior by the hypothesis that the elliptic SPDE models the physical processes occurring during manufacturing, like heat and mass diffusion. Although the proposed approach is general, we demonstrate its application to lattices modelled as pin-jointed trusses with uncertainties in member Young's moduli. We consider as a cost function the weighted sum of the expectation and standard deviation of the structural compliance. To compute the expectation and standard deviation and their gradients with respect to member cross-sections we use a first-order Taylor series approximation. The cost function and its gradient are computed using only sparse matrix operations. We demonstrate the efficiency of the proposed approach using several lattice examples with isotropic, anisotropic and non-stationary random fields and up to eighty thousand random and optimisation variables.


Introduction 1.Motivation
The inherent random variations in material properties and geometry resulting from manufacturing have adverse implications for the performance and design of optimised structures in general and, more specifically, lattice structures considered in this paper.Lattice structures in the form of microarchitected lattices, or mechanical metamaterials, have recently attracted much attention due to their superior mechanical performance compared to conventional bulk solids (Fleck et al, 2010;Gibson and Ashby, 1999;Zheng et al, 2014).Technological advances in additive manufacturing, especially over the past decade, have made it possible to produce lattice structures of various topologies with geometric features spanning multiple length-scales (Gibson et al, 2021;Schaedler and Carter, 2016;Maconachie et al, 2019).The vast design space afforded by additive manufacturing can be utilised using structural optimisation, in particular by optimising the member crosssections and topology.
In structural optimisation it is essential to consider the uncertainties in material properties (Doltsinis and Kang, 2004;Asadpoure et al, 2011;Da Silva and Cardoso, 2017), geometry (Guest and Igusa, 2008;Chen and Chen, 2011;Jansen et al, 2015), loading (Ben-Tal and Nemirovski, 1997;Kogiso et al, 2008;Zhao and Wang, 2014;Torres et al, 2021) and manufacturing parameters (Wang et al, 2011;Schevenels et al, 2011).A design obtained using conventional deterministic optimisation is usually sensitive to variations in material properties, loading and geometry, often eradicating the optimality of the structure.In robust optimisation, uncertainties are taken into account by considering a cost function consisting of the weighted sum of the expectation and standard deviation of standard cost functions used in deterministic optimisation; see Park et al (2006); Beyer and Sendhoff (2007); Schuëller and Jensen (2008) for an overview.Consequently, robust optimisation yields a structure which is optimal in terms of a chosen cost function and is less sensitive to the uncertain properties of the actual structure.For completeness, we note that alternative approaches for considering uncertainties, including reliability-based structural optimisation, are not discussed further in this paper.

Related research
The variations in material properties and geometry of the lattice structure can be represented as spatial random fields, or processes.We assume in this paper that the respective random fields are approximated as Gaussian fields with a mean and a covariance.The covariance, or kernel, function describes the correlation between the field values at different positions in space, and consequently depends on the specific manufacturing and assembly process used.Especially, in the case of the Young's modulus the covariance function must depend on the manufacturing process and the partial differential equations corresponding to the underlying physical processes.For instance, in additive manufacturing heat is used to melt/sinter powders or to melt filaments to make threedimensional lattices layer by layer (Gibson et al, 2021).Consequently, it is expedient to consider the Laplace operator modelling the heat flow in establishing the covariance function of random fields.Indeed the formal link between PDEs and covariance properties of random fields is widely known in spatial statistics (Whittle, 1954;Rue and Held, 2005;Lindgren et al, 2011Lindgren et al, , 2022)).For instance, Whittle (1954) showed that a Gaussian random field with the widely used Matérn covariance function is equivalent to the solution of a stochastic PDE with a Gaussian white noise as the forcing.The respective PDE with a second-order elliptic operator is of fractional order.In other words, after discretising it with finite elements, the fractional power of the corresponding stiffness matrix is equal to the precision matrix, i.e. the inverse of the covariance matrix, of the random field (Lindgren et al, 2011).The so-obtained covariance matrix can be interpreted as a physicsinformed prior describing the spatial correlation in the random field.
The benefits of the PDE representation of the covariance include that the resulting matrices are usually sparse and that the PDE can be easily generalised to anisotropic and non-stationary (non-homogenous) random fields on bounded and non-Euclidean domains.For lattices, the correlated random field can be obtained by considering the discretisation of the PDE operator on the lattice.That is, the respective global diffusion matrix is assembled from the one-dimensional element diffusion matrices and the global mass matrix from the one-dimensional element mass matrices.Subsequently, the covariance matrix is determined using standard procedures developed for the nonlattice case (Bolin and Kirchner, 2020;Koh and Cirak, 2023).The sketched approach is somewhat similar to the graph kernels recently proposed in statistical machine learning (Borovitskiy et al, 2020(Borovitskiy et al, , 2021;;Nikitin et al, 2022;Bolin et al, 2023).
Robust optimisation is computationally significantly more costly than deterministic optimisation because of the dependence of cost and, possibly, constraint functions on the second-order statistics of the finite element solution.The probability density of the solution is the push-forward of the probability density of the material and other random parameters by the lattice solution operator.Hence, even though the input probability densities are Gaussians, the probability density of the finite element solution is usually non-Gaussian and is analytically intractable.A straightforward but computationally prohibitive approach is to use Monte Carlo (MC) sampling (Haldar and Mahadevan, 2000) to estimate the expectation and variance of the objective and cost functions and their gradients (sensitivities) in each optimisation step.MC sampling alone is hardly scalable beyond a handful of random variables and is consequently of limited practical use.A multitude of approaches has been explored to reduce the number of design variables and the number of required finite element solutions to estimate the statistics of the cost and constraint functions by MC sampling.Perhaps the most common approach is to reduce the dimensionality of random fields by using a stochastic series expansion, including Karhunen-Loéve (Sudret and Der Kiureghian, 2000) or the expansion optimal linear estimator (EOLE) (Li and Der Kiureghian, 1993).In structural optimisation series expansions have been used, e.g., in Schevenels et al (2011); Lazarov et al (2012b,a); Jansen et al (2013); Zhao and Wang (2014).
Alternatively, the number of finite element evaluations can be drastically reduced with a truncated Taylor series expansion of the finite element solution with respect to the random parameters to estimate expectation and variance of cost and constraint functions (Doltsinis and Kang, 2004;Guest and Igusa, 2008;Asadpoure et al, 2011;Lazarov et al, 2012b).Usually, only the first-order linear term in the Taylor series expansion is considered so that the probability density of the structural response is a Gaussian, given that the input densities are Gaussians, yielding the so-called perturbation methods (Kriegesmann and Lüdeker, 2019).The perturbation methods inevitably lead to truncation errors.The coefficient of variation of the random parameters must be small to obtain meaningful approximations.
Instead of the series expansion, other approaches for solving stochastic PDEs (Matthies and Keese, 2005), like the stochastic Galerkin method using a polynomial chaos expansion (PCE) (Ghanem and Spanos, 1991), can be utilised to estimate the second order statistics of the random variables and fields required in robust optimisation (Tootkaboni et al, 2012;Keshavarzzadeh et al, 2016).Although the PCE based approaches are appealing, they are unsuitable for problems with a large number of random variables.
For completeness, methods which replace the finite element model with a pre-trained surrogate model, or a metamodel, are worth mentioning as well (Jin et al, 2003;Santner et al, 2003;Sobester et al, 2008).The surrogate model can be, for instance, a Gaussian process regression model (Williams and Rasmussen, 2006) or a more sophisticated probabilistic machine learning model (Vadeboncoeur et al, 2023) trained on a dataset consisting of finite element solutions.The surrogate model is very fast to evaluate and consequently makes MC sampling practical.More recent techniques to speed up the MC sampling include the use of stochastic gradient-based methods prominent in machine learning (De et al, 2020).

Contributions
In this paper, we introduce a novel approach for the robust topology optimisation of lattice structures with up to several tens of thousands of members, see Figure 1.We focus on correlated random fields, especially their efficient representation and consideration within gradient-based robust optimisation.The lattice structure is modelled as a pin-jointed truss and consists of a large number of unit cells each of which have a small number of members.This approximation is sufficient for technologically relevant stretch-dominated lattices with appealing mechanical properties; see, e.g.Deshpande et al (2001).We optimise the topology of the lattice by taking the cross-sectional areas of the members as design variables (Xiao and Cirak, 2022;Christensen and Klarbring, 2008).Although the proposed approach is general, we demonstrate its application to lattices with uncertainties in member Young's moduli, and the cost function is the weighted sum of the expectation and standard deviation of the compliance.
The probability density of the member Young's moduli is given by a multivariate Gaussian distribution with a covariance matrix which we determine by generalising the SPDE representation of random fields to lattices.Because Young's moduli are associated with the members of the lattice, we consider the adjoint lattice for discretising the stochastic PDE.By interpreting the lattice structure as a graph, its adjoint (or, line) graph is another graph representing the adjacencies between the members of the lattice (Harary and Norman, 1960).The adjoint lattice is defined by placing the vertices of the adjoint graph at the centroids of the members of the original lattice.As mentioned, the precision matrix of the random field is equal to the stiffness matrix of the discretised fractional PDE with a second-order elliptic operator on the adjoint lattice.The precision matrix is sparse when the PDE is of integer order and can be approximated as sparse when non-integer order (Bolin and Kirchner, 2020;Koh and Cirak, 2023).We consider only integer order operators in this paper.
Furthermore, we use a first-order Taylor series approximation to determine the second-order statistics of the random compliance.The mean and standard deviation of the compliance and their gradients with respect to the cross-sectional areas are determined using only sparse matrix operations given that the precision matrix of the random Young's moduli is sparse.Consequently, the proposed approach can be applied to lattices with a very larger number of members and joints.
In addition, we apply a density filter to avoid scattered members (similar to checkerboards in continuum structures) and penalisation to obtain cross-sectional areas within a certain prescribed range (Xiao and Cirak, 2022).

Robust design optimisation
In this section, we review the robust topology optimisation of linear elastic truss structures.That is, the minimisation of the member cross-sectional areas for a prescribed loading and material volume.The Young's moduli of the members is assumed to be random and spatially correlated.

Problem formulation
The equilibrium equation for a pin-jointed truss structure with n e members and n d degrees of freedom is given by where f ∈ R n d and u ∈ R n d are the displacement and force vectors, and K ∈ R n d ×n d is the positive definite symmetric global stiffness matrix.
The stiffness matrix and the displacements depend on member design variables and Young's moduli collected in the vectors s ∈ R ne and r ∈ R ne , respectively.The design vector s consists of the relative (pseudo) density of the members defined as s e = A e − A min A max − A min . (2) Here, A e is the cross-sectional area of the member e and A max is the maximum allowable crosssectional area.Furthermore, A min is a small algorithmic parameter to avoid the ill-conditioning of the stiffness matrix when A e ≈ 0. The crosssectional area and Young's modulus of each member is uniform across its length.The Young's modulus of the members is random and has the multivariate Gaussian density where r ∈ R ne is its expected value and C r ∈ R ne×ne its positive definite covariance matrix.
Both are defined as using the expectation operator E[•].The covariance matrix is usually dense and becomes a diagonal matrix when the components of the random vector are uncorrelated, and in case of Gaussian densities statistically independent.Furthermore, the precision matrix Q r is defined as the inverse of the covariance matrix, i.e.Q r = C −1 r .The robust compliance optimisation problem for the truss structure with random member Young's moduli can now be stated as The compliance is defined as where the dot denotes inner product and the last expression is obtained by considering the total potential of the structure at equilibrium.The cost function F (s) is the sum of the expectation and variance, weighted by the two user-chosen parameters α E , α σ ∈ R + .The choice of the two parameters will be detailed in Section 4 on examples.The constraint (5c) ensures that the volume V (s) = l • a(s) of the optimised structure is below the prescribed volume V max , where l ∈ R ne and a(s) ∈ R ne are the vectors of member lengths and cross-sectional areas, respectively.The last constraint (5d) enforces that the relative density s e of the members is between 0 and 1.

Perturbation approximation
The fast evaluation of the cost function and its gradient are critical in iterative structural optimisation.As specified in (5a) the cost function is the weighted sum of the expectation and standard deviation of the compliance.We use a first-order Taylor series expansion of the displacements to approximate both.The first-order series expansion of the random displacement vector u(s, r) around the expected value r of the random vector is given by The gradient of the solution herein is derived by differentiating the equilibrium equation ( 1) yielding We use in the following the same symbols for the approximate vectors and matrices like their exact counterparts to avoid a proliferation of symbols.
Owing to the linear transformation property of Gaussian vectors, see, e.g.Williams and Rasmussen (2006), the (approximate) probability density of the displacement vector is following (8), ( 9) and (3) also a Gaussian, i.e., p (u(s)) = N (u(s), C u (s)) . (10) According to the series expansion (8), we obtain for the expected value and for the covariance matrix Furthermore, using one more time the linear transformation property of Gaussian vectors, we arrive at the probability density of the compliance with the mean and variance For later derivations, we can can express the variance more compactly as where it has been taken into account that

Sensitivities
We determine the gradient, i.e. sensitivity, of the cost function (5a) using the first-order approximation of the random displacements (8).To this end, it is beneficial to consider the cost function as the sum of the member contributions, i.e., By recourse to the approximate expectation (13b) and the compliance (6), the contribution of member e to the gradient of the expected compliance is given by ∂J e (s) ∂s e = −u e (s, r) • ∂K e (s e , r) ∂s e u e (s, r) .(17) Here, s e denotes the design variable, K e the element stiffness matrix and u e the nodal displacement vector of member e.Similarly, we can obtain from the approximate variance ( 14) the member e's contribution to the gradient of the standard deviation as ∂σ J,e (s) ∂s e = 1 2σ J,e (s) where where l e and A e are the length and cross-sectional area of member e.

Penalisation and filtering
We use in some of the examples included in this paper a penalisation approach akin to the SIMP method (Bendsøe and Sigmund, 1999) to obtain optimised structures with a desired crosssectional area distribution.Specifically, we aim to reduce the number of members with a relative density s e < s * .The prescribed parameter s * will usually depend on the specific manufacturing process used, and herein is set to s * = 0.5.Figure 2 shows the used the piecewise C 1 continuous penalised density se (s e ) obtained by combining a univariate B-spline curve and a line.The B-spline is of polynomial degree four, has six control points, and the begin and end points are interpolating.The choice of different penalisation functions has been explored by Xiao and Cirak (2022).
Furthermore, we regularise the optimisation problem by applying a filter to the design variables to avoid checkerboard-like instabilities known from topology optimisation for continuum structures.The used cone kernel function evaluated at the centroidal coordinates c i and c j of the members i and j reads where R is the prescribed filter length and ∥ • ∥ the Euclidean distance.With the obtained filtering weights the filtered design variable of a member e is given by ŝe As proposed in Xiao and Cirak (2022) this definition includes the member length l i rendering the cost function gradient a measure of the strain energy density (elastic energy per unit volume) rather than the total compliance.To apply filtering to the gradients, the vector of design variables s in the cost function ( 16) and the constraint ( 21) is replaced with its filtered counterpart ŝ.Subsequently, their gradients are determined using the chain rule, that is, 3 Spatial random fields We now turn to the computation of the multivariate Gaussian probability density of the member Young's moduli.After reviewing the stochastic PDE representation of Matérn fields and its finite element discretisation on bounded Euclidean domains, we present its generalisation to lattice structures.Our focus is on obtaining the precision matrix corresponding to a generalised Matérn field with prescribed hyperparameters, such as length-scale, smoothness, variance and anisotropy.

Matérn fields on unbounded Euclidean domains
We consider in R d , d ∈ {1, 2, 3}, a zero-mean Gaussian process, or a (spatial) Gaussian random field, r(x) ∼ GP (0, c r (x, x ′ )) , with the Matérn covariance function c r (x, x ′ ) defined as with σ ∈ R + being the standard deviation, ν ∈ R + the smoothness parameter, ℓ ∈ R + the lengthscale parameter, Γ the Gamma function and K ν the modified Bessel function of the second kind of order ν.The parameter κ is defined as Figure 3 shows the covariance function for smoothness parameters ν ∈ {0.5, 1.5, 2.5} and realisations of the corresponding Gaussian zero-mean random fields.As apparent, the covariance function and the respective samples become smoother with increasing ν.
According to Whittle (1954) and Lindgren et al (2011), the random field r(x) is the solution of the stochastic partial differential equation (SPDE) where ∆ is the Laplace operator, g(x) a Gaussian white noise process and the remaining two parameters are defined as Although the exponent β can take any value β > d/4, we choose in this paper β ∈ Z >0 to avoid the need to consider fractional derivatives.Following Lindgren et al (2011), for a smoothness parameter β > 1 the SPDE is solved using the recursion i.e. r(x) ≡ r (β) (x) with k = 2, . . .β. Without going into details, a straightforward finite element discretisation of the PDEs in the recursion yields ) where M is the (lumped) mass matrix, A the stiffness matrix and g a Gaussian white noise vector with the density g ∼ N (0, M ) ; (34) see Koh and Cirak (2023).
Owing to the linear transformation property of Gaussian random vectors, the solution vector r is a Gaussian and has the probability density with the precision matrix (36) We refer to Appendix A for further details on this derivation and the computation of the precision matrix.

Matérn fields on lattices
The stochastic PDE representation of random fields is instrumental in generalising Matérn fields to lattices.The naive approach of evaluating the covariance function ( 27) by introducing the coordinates of lattice points does not take into account the topology of the lattice, and simply replacing the distance with the geodesic, or shortest-path distance, does not always lead to a valid, i.e. positive definite, covariance matrix (Gneiting, 1998).However, starting from the finite element discretised random field (35) and replacing the mass and stiffness matrices with the corresponding matrices of the lattice structure does indeed lead to a positive definite covariance matrix.
As noted, Young's modulus of a member is assumed to be constant along its length, so that the members attached to a joint may have different Young's moduli.This suggests using the adjoint lattice instead of the actual lattice for determining the precision matrix of the random vector r.
That is, we interpret the lattice as a graph G(V, E) where the joints and members of the lattice correspond to the vertices V and edges E, respectively.In the adjoint graph L(G), defining the adjoint lattice, the edges of G become the vertices of L(G) (Harary and Norman, 1960).The adjoint lattice Fig. 4: Two lattices (solid black) and their adjoint lattices (dashed red).The vertices of the adjoint lattice are placed at the centroid of the members of the original lattice.
joints are placed at the centroids of the lattice members.In Figure 4, two example lattices and their adjoints are depicted.
It is straightforward to determine the mass and stiffness matrices of the so-defined adjoint lattice.For instance, for a lattice chain in R 1 the lumped mass and stiffness matrices of a member e are given by where h e is the length of the member.After assembling the element matrices as in standard finite elements, the obtained matrices are introduced into (36) to obtain the precision matrix for the adjoint lattice.As usual, the mass and stiffness matrices of a member embedded in R 2 or R 3 are obtained by applying a coordinate transformation to (37).
Figure 5 depicts the isocontours of random fields sampled from SPDE-defined Gaussian probability densities.As apparent, the obtained random fields depend on the geometry and topology of the lattice structure and its unit cell connectivity.The two lattices in Figures 5a and 5b have the same unit cell connectivity with two diagonal members.Comparing the isocontours in Figures 5a and 5b, the four holes in Figure 5b lead to a change in the overall appearance of the random field.Specifically, the random field values at two opposing ends of a hole are very different despite their proximity in terms of Euclidean distance.This is unsurprising given that we determine the random field by solving the system of equations (33) containing the lattice stiffness and mass matrices and a Gaussian white noise forcing.We chose as boundary conditions for the nodes adjacent to holes homogenous Neumann and the nodes at the outer boundary as infinite domain; see Roininen et al (2014); Khristenko et al (2019) for discussions on the choice of boundary conditions in case of Euclidean domains.
Furthermore, Figure 5c shows a lattice structure with the same geometry and topology but with each unit cell having only one diagonal.This leads again to a change in the overall appearance of the random field in comparison to the lattice structure with two diagonal members in each cell.The observed dependency of the random field on the geometry, topology and unit cell connectivity agrees with our hypothesis that the SPDE models the diffusion processes, i.e. heat or mass diffusion, during manufacturing.Consequently, the obtained random fields represent a sound physics-informed prior for properties of lattices manufactured, e.g., by additive manufacturing or casting.

Non-standard Matérn fields on lattices
The SPDE representation of random fields makes it possible to easily model non-stationary and anisotropic random fields without altering the introduced overall solution process.Specifically, a non-stationary random field is modelled by choosing a spatially varying length-scale ℓ(x) yielding the SPDE parameters To model anisotropic random fields, e.g.resulting from the particulars of the manufacturing process used, we replace the first equation (33a) in the recursion (33) by Here, D is a diagonal matrix with the components where t e ∈ R d is the unit tangent to member e of the lattice structure, n ∈ R d is a prescribed unit vector describing the direction of anisotropy, and d || and d ⊥ ∈ R + are two parameters.In practice, the anisotropy vector n and the two parameters d || and d ⊥ will depend on the specific manufacturing process and protocol used.
Figure 6 shows non-stationary random fields sampled from SPDE-defined Gaussian probability densities.The length-scale of the random field is uniform in Figures 6a and 6c (ℓ = 3 and ℓ = 20, respectively).In Figure 6b the nonuniform length-scale increases linearly from ℓ = 3 at the left end to ℓ = 20 at the right end.The increase of the length-scale is easily discernible from the isocontours of the random field.

Examples
We proceed to demonstrate the utility and efficiency of the proposed robust optimisation approach with selected examples.In all the examples, the random imperfection concerns the Young's moduli of the members.We consider the cost function (5a), which is the weighted sum of the expectation and standard deviation of the structural compliance.We choose the respective weights in dependence on a single prescribed cost function parameter α as The parameter α is chosen such that 0 ≤ α ≤ 1.The normalisation constant J * is the cost function value obtained by optimising with α = 1, and the normalisation constant σ * J is the cost function value obtained with α = 0. Note that the robust optimisation problem can also be interpreted as a multicriteria optimisation so that the obtained solution for a given α represents a Pareto point (Beyer and Sendhoff, 2007).We solve the discrete optimisation problem in all examples using the Method of Moving Asymptotes (MMA) (Svanberg, 1987).

Verification example
To verify the correctness of our derivations and implementation, we consider the lattice structure depicted in Figure 7, which was previously investigated by Asadpoure et al (2011).The lattice consists of 2 × 4 unit cells of size 0.75 × 1.There are in total 15 joints and 38 members.The joints at the left end are fixed, and a single external force of magnitude 1 is applied at the centre of the right end.The prescribed maximum volume is V max = 0.5 and the cross-sectional member areas are constrained not to exceed A max = 1.The Young's moduli of the members are random, spatially uncorrelated and have the probability density where I ne is an all-ones vector and I ne×ne an identity matrix.We note that the corresponding finite element solution has a non-diagonal covariance matrix as given in ( 12).When optimising only with respect to the expected compliance, i.e. by choosing α = 1, in the optimised structure shown in Figure 8a only the members along the centre of the structure are present while all the others are missing.Note that this case is equivalent to conventional deterministic optimisation.Evidently, the obtained structure cannot be considered as robust given that there is only a single load path, and the compliance depends strongly on the random Young's moduli of the four members.To improve robustness, we take into account the standard deviation of the compliance by choosing a smaller α value.As expected, the optimised structure for α = 0.3 depicted in Figure 8b has many more members, a slightly higher compliance and a significantly smaller standard deviation as desired.Irrespective of its practicality, it is feasible only to minimise the standard deviation of the compliance while neglecting its expected value by choosing α = 0.The respective optimised structure is shown in Figure 8c.It bears emphasis that the material volume is the same in all the three optimised structures in Figure 8. Finally, we confirm that the obtained optimised structures are in close agreement with the results reported in Asadpoure et al (2011), verifying the correctness of our implementation.

Tension strip
In this example, we optimise the two-dimensional tension strip depicted in Figure 9 and study the effect of the spatial properties of the random Young's moduli.The precision matrix of the respective probability density function is determined using the SPDE formulation of random fields and its extension to lattice structures introduced in Section 3.2.The width and height of the tension strip are 30 × 20.The joints along the top boundary are fixed, and all the remaining boundary joints are unconstrained.A single load of magnitude 1 is applied at the centre on the bottom edge.Each unit cell is of size 1 × 1 and has two diagonal members.The structure has 2450 members, and the mean of Young's modulus for all members is E = 100.Furthermore, the total allowable volume is V max = 6 and the cross-sectional member areas are constrained not Fig. 8: Verification example.Optimised lattice structures for the cost function parameter values α ∈ {1, 0.3, 0}.Optimising only with respect to the expectation of the compliance, i.e. α = 1, yields an optimised structure consisting of only a colinear chain of four members.Including the standard deviation of the compliance by choosing α < 1 the members become more dispersed.For α = 0 the cost function consists only of the standard deviation of the compliance.
to exceed A max = 1.The radius of the density filter is chosen as R = 2.5, i.e. 2.5 times the unit cell size.

Isotropic, stationary covariance
First, we investigate the effect of the covariance length-scale ℓ ∈ {0.5, 5.0, 10.0} and the cost function weight parameter α ∈ {1, 0.5, 0} on the optimisation results.The covariance smoothness parameter and the standard deviation are ν = 1.5 and σ = 10.For illustration purposes, in Figure 10, three samples from the considered random fields with different length-scales are shown.Note that E(x) − E is equal to the solution of the SPDE (33) with a Gaussian white noise forcing.We use this fact to sample efficiently large random fields like the ones depicted.In Figure 10, the members are coloured according to their Young's The nine optimised lattice structures for different combinations of length-scale parameters ℓ and cost function weight parameters α are depicted in Figure 11.When the cost function is only the expectation of the compliance, i.e. α = 1, the optimal structure and the respective expected compliance J do not depend on the length-scale.However, an increase in ℓ still leads to a rise in the standard deviation of the compliance σ J .The cost function depends on the standard deviation when choosing α < 1.As seen in the middle and right columns in Figure 11, this leads to designs with a more spread-out member layout.The extent of the spreading is inversely proportional to the length-scale ℓ.It is intuitively clear that a decrease in the length-scale leads to increasingly uncorrelated Young's moduli, in which case the best robust optimisation strategy consists of distributing the load to as many as possible members while satisfying all the constraints.

Non-stationary, anisotropic covariance
We now optimise a tension strip with a random Young's modulus with a non-stationary isotropic covariance with a length-scale ℓ(x) varying linearly from the top to the bottom.The obtained optimised structures for two different lengthscale distributions are compared in Figure 12.In Figure 12a, the length-scale increases linearly from the bottom to the top according to ℓ(y) = 0.5(0.95y+ 1), so that the length-scale at the top is ℓ(20) = 10 and at the bottom ℓ(0) = 0.5.In contrast, in Figure 12b, the lengthscale increases linearly from the top to the bottom according to ℓ(y) = 0.5(−0.95y+ 20), leading to ℓ(20) = 0.5 and ℓ(0) = 10.The other covariance parameters σ = 10 and ν = 1.5 are the same for both structures.The Young's moduli vary more rapidly in shorter length-scale regions, prompting the algorithm to spread the material among more members.At last, we consider a tension strip with two different anisotropic random Young's modulus distributions, see Figure 13.The prescribed unit vector n gives the direction of the anisotropy.In the first case (Figure 13a) with d || = 1 and d ⊥ = 5 the members perpendicular to n have higher uncertainty.The robust optimisation leads to a narrower structural system as the load is carried primarily by the vertical members with a lower uncertainty.In contrast, in the second case (Figure 13b) with d || = 5 and d ⊥ = 1 the load is carried mainly by the diagonal members due to the high uncertainty associated with the vertical members.

Cantilever
We consider next a cantilever lattice structure of size 20 × 10 with unit cells of size 0.5 × 0.5, see Figure 14.Each unit cell has two diagonal members, so there are 861 joints and 3260 members in total.A vertical load of magnitude 100 is applied at the mid-height on the right end of the lattice while the joints along the left end are fixed.The uncorrelated member Young's moduli have a constant mean value of E = 7 • 10 7 and a standard deviation of σ = 7 • 10 6 .The prescribed maximum volume is V max = 1.6, and the cross-sectional areas are constrained not to exceed A max = 5.1 × 10 −3 .The density filter radius is chosen as R = 1.
In this example, we penalise the relative member densities se (s e ) with the penalisation function introduced in Figure 2 to obtain a manufacturable structure with cross-sectional areas larger than a prescribed threshold.The optimised structures in Figure 15 correspond to cost function weight parameter values α ∈ {1, 0.5, 0.1} and computations with and without penalisation.Selected histograms of the respective cross-sectional areas are given in Figure 16.It is evident that in the case of no penalisation, there are a large number of, likely impossible to manufacture, members with very small cross-sectional areas.
Furthermore, although the shape of the optimised structures is almost insensitive to α, the mean and standard deviation of the compliance are sensitive to the choice of α.

Engine bracket
As a final example, we demonstrate the optimisation of a bracket lattice structure (loosely motivated by the GE jet engine bracket challenge) to demonstrate the application of the proposed approach to large problems.We simplify the original engine bracket geometry while retaining its essential features, see Figures 1 and 18.In particular, the four holes for fixing the horizontal base plate and the hole for applying the load to the vertical plate are omitted.The dimensions of the bracket are given in Figure 19.The lattice consists of body-centred cubic (BCC) unit cells of size 5 × 5, leading to 82124 members and 12971 joints.
In the finite element model, all the joints on the reentrant surfaces of the base plate are fixed, and a unit load is applied to the 62 joints at the top left corner of the hole in the vertical plate.As illustrated in Figure 18, each unit load has a slope of 45 deg with respect to the horizontal base plate.The prescribed maximum volume is V max = 15200, and the cross-sectional areas are constrained not to exceed A max = 1.0.The relative member densities se (s e ) are penalised using the penalisation function introduced in Figure 2. Furthermore, the mean Young's modulus of the members is E = 100, and for the SPDE-defined covariance matrix the length-scale is ℓ = 10, the smoothness parameter is ν = 1.5, and the standard deviation is σ = 0.1E.
Figure 20 shows the optimised lattice structures for three different cost function parameter values α ∈ {1, 0.2, 0}.As discussed, a decrease in the value of α leads to increased consideration of the standard deviation of the compliance in the cost function.For the present structure, a decrease in α results in a concentration of material around the top ends of the four legs of the bracket.For instance, reducing from α = 1 to α = 0.2 yields a reduction of the standard deviation from 6.19 to 3.67 and an increase of the mean by 5%.At last, in terms of computational efficiency, in each iteration step, it is sufficient to factorise the stiffness matrix of the truss structure only once.Subsequently, the derivative of the displacement vector with respect to the member Young's moduli is computed using back substitution, cf.(20).For each of the 82124 members, it is necessary to apply one back substitution which may represent a performance bottleneck in the case of large problems.However, each member's derivative of the displacement vector can be independently computed so that this step is easily parallelisable.In robust optimisation with α = 0.2, the CPU runtime for one iteration is 18106.19s (1783.77s wall time), taking in total 127 iterations.For comparison, the run-time for a coarse lattice with 3226 members is 16.56 s (3.44 s wall time), and for a lattice with 10676 members it is 200.72 s (38.53 s wall time).The optimisation results for the two coarse meshes are not included in the paper.All experiments were run on an Intel®Xeon®Silver 4116 CPU @ 2.10 GHz (12 cores, 24 logical processors).

Conclusions
We introduced an efficient approach for topology optimisation of large lattice structures with spatially correlated random material properties.Although we considered only members with random Young's moduli, it is straightforward to adapt the proposed approach for other random fields, including the member cross-sectional areas, load distribution and joint positions, by following Guest and Igusa (2008).We obtain the spatially correlated random fields by generalising the SPDE formulation of Matérn fields to lattice structures.In the SPDE formulation, the random field is defined as the solution of a possibly fractional PDE with a second-order elliptic operator and a Gaussian white noise as the forcing.The properties of the resulting random field are governed by the parameters of the SPDE, i.e. the  length-scale, the exponent and the standard deviation of the white noise.We generalise the SPDE formulation to lattices by replacing the stiffness and mass matrices of the finite element discretised SPDE with the respective matrices of the lattice structure.The resulting system matrix is the precision matrix of the Gaussian random field, which is usually sparse or can be approximated as sparse in the case of fractional SPDEs.We hypothesise that the SPDE models the unaccounted physical diffusion processes during the manufacturing of the lattice structure and, hence, is an effective approach to constructing physics-informed random fields.This interpretation of the SPDE formulation of random fields helps generalise them to non-stationary and anisotropic random fields.For robust topology optimisation by optimising the cross-sectional areas of the members, we use a standard gradient-based approach and compute the required mean and standard deviation of the compliance and constraint functions using a first-order Taylor series expansion.The efficiency of the series expansion combined with the sparsity of the precision matrix of the random field allows us to optimise large lattice structures with up to eighty thousand members with random Young's moduli.These computations are only feasible because the expectation and standard deviation of the cost and the constraint functions and their gradients are computed using sparse matrix operations.
We note that it is straightforward to apply the proposed robust optimisation approach to other types of lattice structures, including lattice-skin systems composed of a lattice core and a thinshell skin (Xiao et al, 2019;Xiao and Cirak, 2022) or structures consisting of beam members connected by joints that can transfer both forces and moments (Yin et al, 2020).Furthermore, we assumed in the present paper that the parameters of the SPDE are given.They will depend on the particulars of the manufacturing process, like the deposition rate and laser power in additive manufacturing by selective laser sintering (Gibson et al, 2021).Relating the SPDE parameters to such directly observable manufacturing process parameters is of great practical relevance for designing more robust structures.

Fig. 1 :
Fig. 1: Robust topology optimisation of a bracket consisting of body-centred cubic unit cells.The bracket is fixed at the four corners of the base plate and an external load is applied at the horizontal shaft attached to the vertical plate.For further details on this example see Section 4.4.

Fig. 2 :
Fig. 2: Penalisation function mapping the relative density s e to the penalised relative density se .The piecewise penalisation function consists of a B-spline curve and a line, which are smoothly connected at s * e = 0.5.The dashed polygon and its corners represent the control polygon of the B-spline curve.The penalisation ensures that densities s e = 0 and s e > s * are preferred.

Fig. 5 :
Fig.5: Random fields sampled from SPDEdefined Gaussian probability densities.The size of each lattice is 50 × 30 and the size of each unit cell is 1 × 1.The lattice structures in (a) and (b) have unit cells with two diagonals and the lattice structure in (c) has unit cells with a single diagonal.In all case the length-scale parameter is ℓ = 20, the smoothness parameter is ν = 2.0, the standard deviation is σ = 0.1r, and the mean is r = 100.

Fig. 6 :
Fig. 6: Stationary and non-stationary random fields sampled from SPDE-defined Gaussian probability densities.The size of each lattice is 50 ×30, the size of each unit cell is 1 × 1 and each unit cell has two diagonals.In (a) and (c) the length-scale is uniform throughout the domain and in (b) it increases linearly from ℓ = 3 at the left boundary to ℓ = 20 at the right boundary.The smoothness is ν = 2.0, the standard deviation is σ = 0.1r, and the mean is r = 100.
Fig. 11: Tension strip with isotropic, stationary covariance.Optimised structures for different combinations of length-scale ℓ and cost function parameter α.For each structure the minimum expected compliance J and standard deviation of the compliance σ J are given at the bottom left of each subfigure.

Fig. 15 :Fig. 16 :
Fig. 15: Cantilever.Comparison of optimised structures with and without penalisation of the relative member densities for different cost function weight parameter values.

Fig. 17 :
Fig. 17: Cantilever.Pareto fronts of the multiobjective robust optimisation problem obtained by varying the cost function weight parameter α.The two B-spline penalised curves correspond to different penalisations of the relative member densities, i.e. according to Figure 2 and a (not shown) less curved function.

Fig. 20 :
Fig. 20: Engine bracket.Comparison of the optimised engine bracket for three different cost function weight parameters.For α = 1 (left column) the cost function consists only of the mean and for α = 0 (right column) only of the standard deviation of the compliance.