Topology optimization with worst-case handling of material uncertainties

TL;DR

This paper develops a topology optimization method that accounts for material uncertainties in a worst-case scenario, using a minimax formulation with convex relaxation and regularization schemes, demonstrated on additive manufacturing and degradation cases.

Contribution

It introduces a novel worst-case approach for topology optimization under material uncertainties, with a convex relaxation and regularization schemes for efficient solution.

Findings

The method effectively handles material uncertainties in practical examples.

Barrier regularization improves the solution process compared to Tikhonov regularization.

A continuation scheme mitigates over-conservativeness for large uncertainties.

Abstract

In this article a topology optimization method is developed, which is aware of material uncertainties. The uncertainties are handled in a worst-case sense, i.e. the worst possible material distribution over a given uncertainty set is taken into account for each topology. The worst-case approach leads to a minimax problem, which is analyzed throughout the paper. A conservative convex relaxation for the inner maximization problem is suggested, which allows to treat the minimax problem by minimization of an optimal value function. A Tikhonov type and a barrier regularization scheme are developed, which render the resulting minimization problem continuously differentiable. The barrier regularization scheme turns out to be more suitable for the practical solution of the problem, as it can be closely linked to a highly efficient interior point approach used for the evaluation of the optimal value function and its gradient. Based on this, the outer minimization problem can be approached by a gradient based optimization solver like the method of moving asymptotes. Examples from additive manufacturing as well as material degradation are examined, demonstrating the practical applicability as well as the efficiency of the suggested method. Finally, the impact of the convex relaxation of the inner problem is investigated. It is observed that for large material deviations, the robust solution may become too conservative. As a remedy, a RAMP-type continuation scheme for the inner problem is suggested and numerically tested.

Tables

Table 1. Table 1: From left to right, the total uncertainty bound, the compliance of the worst-case topology with homogeneous material E 0 subscript 𝐸 0 E_{0} ( 𝜹 ≡ 0 𝜹 0 \bm{\delta}\equiv 0 ), the compliance of the nominal topology with worst-case material distribution and the compliance of the worst-case topology with worst-case material distribution are displayed. All compliance values are specified as increase in comparison to the compliance of the nominal topology with nominal material, i.e. 𝜹 ≡ 0 𝜹 0 \bm{\delta}\equiv 0 on Ω Ω \Omega . The chosen material parameters are E 0 = 1 subscript 𝐸 0 1 E_{0}=1 and E D = 0.7 subscript 𝐸 𝐷 0.7 E_{D}=0.7 in all cases.

total uncertainty	w.-c. topo, $𝜹\equiv 0$	nom. topo, worst-case $𝜹$	w.-c. topo, worst-case $𝜹$
$D=0.03$	+0.182%	+7.69%	+7.25%
$D=0.06$	+0.317%	+11.5%	+10.8%
$D=0.10$	+0.383%	+15.6%	+14.8%

Table 2. Table 2: Compliance increase in comparison to nominal topology with “average material” ( 𝜹 ≡ 0.4 𝜹 0.4 \bm{\delta}\equiv 0.4 ); see Table  1 for a detailed explanation of columns 2 to 4. Material parameters are again chosen as E 0 = 1 subscript 𝐸 0 1 E_{0}=1 and E D = 0.7 subscript 𝐸 𝐷 0.7 E_{D}=0.7 .

total uncertainty	w.-c. topo, worst-case $𝜹$	w.-c. topo, $𝜹\equiv 0.4$	nom. topo, worst-case $𝜹$
$D=0.01$	+0.050%	+3.07%	+2.93%
$D=0.02$	+0.108%	+3.86%	+3.63%
$D=0.03$	+0.170%	+4.32%	+3.98%
$D=0.04$	+0.195%	+4.58%	+4.18%
$D=0.06$	+0.304%	+4.90%	+4.28%
$D=0.10$	+0.331%	+5.08%	+4.29%

Table 3. Table 3: Compliance increases in comparison to nominal topology with fully stiff material E 0 subscript 𝐸 0 E_{0} ( 𝜹 ≡ 0 𝜹 0 \bm{\delta}\equiv 0 ); for a detailed explanation of columns 2 to 4, see Table  1 . Material parameters are now chosen as E 0 = 1 subscript 𝐸 0 1 E_{0}=1 and E D = 0.01 subscript 𝐸 𝐷 0.01 E_{D}=0.01 corresponding to a serious degradation scenario.

total uncertainty	w.-c. topo, worst-case $𝜹$	w.-c. topo, $𝜹\equiv 0$	nom. topo, worst-case $𝜹$
$D=0.0002$	+1.19%	+17.4%	+13.6%
$D=0.0003$	+1.96%	+25.4%	+18.5%
$D=0.0004$	+2.80%	+34.3%	+22.7%
$D=0.0005$	+3.38%	+40.9%	+26.3%
$D=0.0006$	+3.87%	+48.0%	+29.7%
$D=0.001$	+5.71%	+80.1%	+42.5%
$D=0.002$	+7.47%	+148%	+71.7%
$D=0.003$	+8.01%	+216%	+99.8%
$D=0.004$	+8.67%	+283%	+127%
$D=0.008$	+10.5%	+532%	+230%

Table 4. Table 4: Compliance increases in comparison to nominal topology with only stiffest material ( 𝜹 ≡ 0 𝜹 0 \bm{\delta}\equiv 0 ); see Table  1 for a detailed explanation of columns 2 to 4. Material parameters are chosen to be E 0 = 1 subscript 𝐸 0 1 E_{0}=1 and E D = 0.01 subscript 𝐸 𝐷 0.01 E_{D}=0.01 . Worst-case compliances were computed in three different ways: with a RAMP continuation scheme in 10 steps from q = ( E 0 − E D ) / E D 𝑞 subscript 𝐸 0 subscript 𝐸 𝐷 subscript 𝐸 𝐷 q=(E_{0}-E_{D})/E_{D} to q = 0 𝑞 0 q=0 (“contin”, in bold), direct solution without continuation for q = 0 𝑞 0 q=0 (“direct”) and inverse interpolation (“inverse”) for comparison.

total uncertainty	w.-c. topo, $𝜹\equiv 0$	nom. topo, worst-case $𝜹$			w.-c. topo, worst-case $𝜹$
		contin	direct	inverse	contin	direct	inverse
$D=0.0005$	+0.28%	+32.6%	+17.6%	+42.3%	+22.4%	+15.0%	+37.7%
$D=0.001$	+5.75%	+67.0%	+70.5%	+81.6%	+29.6%	+28.7%	+43.1%
$D=0.002$	+8.78%	+139%	+77.9%	+158%	+59.0%	+52.4%	+73.2%
$D=0.004$	+8.62%	+286%	+229%	+306%	+119%	+61.2%	+135%
$D=0.006$	+7.30%	+425%	+281%	+444%	+180%	+138%	+207%
$D=0.008$	+7.55%	+554%	+380%	+575%	+239%	+187%	+270%
$D=0.01$	+7.19%	+675%	+440%	+699%	+296%	+223%	+330%
$D=0.02$	+12.5%	+1203%	+442%	+1234%	+539%	+449%	+597%

Equations

$$(F\ast \mathbf{\rho }{)}_e=\frac{{\sum }_{j\in V(e)}{\rho }_j{\int }_{j}F(\mathbf{x}-{\mathbf{c}}_e)\text{dx}}{{\sum }_{j\in V(e)}{\int }_{j}F(\mathbf{x}-{\mathbf{c}}_e)\text{dx}}=:\frac{{\sum }_{j\in V(e)}{\rho }_j{F}_{ej}}{{\sum }_{j\in V(e)}{F}_{ej}}.$$

$$\mathbf{K}(\mathbf{\rho },\mathbf{\delta })\mathbf{u}=\mathbf{f}$$

$${\mathbf{K}}_e(\mathbf{\rho },\mathbf{\delta })=\sum _{l=1}^{n_{ig}}{\mathbf{B}}_{el}^{\top }{\mathbf{C}}_e\left(\mathbf{\rho },\mathbf{\delta }\right){\mathbf{B}}_{el},$$

$$\hat{\mathbf{C}}(\alpha )=(1-\alpha ){\mathbf{C}}^0+\alpha {\mathbf{C}}^1$$

$${\hat{\mathbf{C}}}_e(\mathbf{\delta })={\left(\frac{1-{\delta }_e}{E_0}+\frac{{\delta }_e}{E_D}\right)}^{-1}\mathbf{C}$$

$$\mathbf{C}=\frac{1}{(1+\nu )(1-2\nu )}\left(\begin{array}{cccc}1-\nu & \nu & \nu & 0\\ \nu & 1-\nu & \nu & 0\\ 0& 0& 0& (1-2\nu )/2\end{array}\right)$$

$${\mathbf{C}}_e(\mathbf{\rho },\mathbf{\delta })={\tilde{\rho }}_e^p{\hat{\mathbf{C}}}_e(\mathbf{\delta })$$

$$\displaystyle\left\{\begin{array}[]{l}\min\limits_{\bm{\rho}\in U_{\text{ad}}}\max\limits_{\bm{\delta}\in\Delta_{\text{ad}}}\bm{f}^{\mathrm{T}}\bm{u}(\bm{\rho},\bm{\delta})\\[8.53581pt] \,\mbox{s.t.:}\ \bm{K}(\bm{\rho},\bm{\delta})\bm{u}=\bm{f}.\\ \end{array}\right.$$

$$\displaystyle U_{\text{ad}}=\bigg{\{}\bm{\rho}\in\mathbb{R}^{n}:$$

$$\displaystyle\sum_{e=1}^{n}v_{e}\rho_{e}/|\Omega|\leq V\bigg{\}},$$

$$\displaystyle\Delta_{\text{ad}}=\bigg{\{}\bm{\delta}\in\mathbb{R}^{n}:$$

$$\displaystyle\sum_{e=1}^{n}v_{e}\delta_{e}/|\Omega|=D\bigg{\}},$$

$$\displaystyle\Delta_{\text{ad}}(\bm{\rho})=\bigg{\{}\bm{\delta}\in\mathbb{R}^{n}:$$

$$\displaystyle\sum_{e=1}^{n}v_{e}\tilde{\rho}_{e}^{p}\delta_{e}/|\Omega|=D\bigg{\}},$$

$$\displaystyle\Delta_{\text{ad}}(\bm{\rho})=\bigg{\{}\bm{\delta}\in[0,1]^{n}:$$

$$\displaystyle\sum_{e=1}^{n}\frac{v_{e}(\delta_{e}-0.4)^{2}}{|\Omega|}\leq D_{2}\bigg{\}},$$

$$\displaystyle\Delta_{\text{ad}}(\bm{\rho})=\bigg{\{}\bm{\delta}\in[0,1]^{n}:$$

$$\displaystyle\sum_{e=1}^{n}\frac{v_{e}(\bm{\tilde{\rho}}_{e}^{p}\delta_{e}-0.4)^{2}}{|\Omega|}\leq D_{2}\bigg{\}},$$

$$\underset{\mathbf{\rho }\in U_{\text{ad}}}{\min }\underset{\mathbf{\delta }\in {\Delta }_{\text{ad}},\mathbf{u}\in {\mathbb{R}}^m}{\max }2{\mathbf{f}}^{\top }\mathbf{u}-\sum _{e=1}^n{\mathbf{u}}^{\top }{\mathbf{K}}_e(\mathbf{\rho },\mathbf{\delta })\mathbf{u}.$$

$$\underset{\mathbf{\rho }\in U_{\text{ad}}}{\min }\underset{\mathbf{\delta }\in {\Delta }_{\text{ad}},\mathbf{u}\in {\mathbb{R}}^m}{\max }{\mathcal{J}}_{ϵ}(\mathbf{\rho },\mathbf{\delta },\mathbf{u}).$$

$${\mathcal{J}}_{ϵ}(\mathbf{\rho },\mathbf{\delta },\mathbf{u})=2{\mathbf{f}}^{\top }\mathbf{u}-\sum _{e=1}^n{\mathbf{u}}^{\top }{\mathbf{K}}_e(\mathbf{\rho },\mathbf{\delta })\mathbf{u}-\frac{ϵ}{2}\Vert \mathbf{\delta }{\Vert }^2,$$

$${\hat{\mathbf{K}}}_e=\sum _{l=1}^{n_{ig}}{\mathbf{B}}_{el}^{\top }\mathbf{C}{\mathbf{B}}_{el}.$$

$$\mathbf{K}=\sum _{e=1}^n{\mathbf{K}}_e=\sum _{e=1}^n{\tilde{\rho }}_e^{p}E({\delta }_e){\hat{\mathbf{K}}}_e.$$

$$\frac{{\partial }^2}{\partial {\delta }_i\partial {\delta }_j}\mathcal{J}$$

$$\frac{{\partial }^2}{\partial {\delta }_i\partial \mathbf{u}}\mathcal{J}$$

$$\frac{{\partial }^2}{(\partial \mathbf{u}{)}^2}\mathcal{J}$$

$${\left(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right)}^{\top }$$

$$={\mathbf{x}}^{\top }\frac{{\partial }^2{\mathcal{J}}_{ϵ}}{{\left(\partial \mathbf{\delta }\right)}^2}\mathbf{x}+2{\mathbf{y}}^{\top }\frac{{\partial }^2{\mathcal{J}}_{ϵ}}{\partial \mathbf{\delta }\partial \mathbf{u}}\mathbf{x}+{\mathbf{y}}^{\top }\frac{{\partial }^2{\mathcal{J}}_{ϵ}}{(\partial \mathbf{u}{)}^2}\mathbf{y}$$

$$=-\sum _{e=1}^n{x}_e^2{\tilde{\rho }}_e^p(E({\delta }_e){)}^3{\left(\frac{1}{E_D}-\frac{1}{E_0}\right)}^2{\mathbf{u}}^{\top }{\hat{\mathbf{K}}}_e\mathbf{u}$$

$$+2\sum _{e=1}^n{x}_e{\tilde{\rho }}_e^p(E({\delta }_e){)}^2\left(\frac{1}{E_D}-\frac{1}{E_0}\right){\mathbf{y}}^{\top }{\hat{\mathbf{K}}}_e\mathbf{u}$$

$$-{\mathbf{y}}^{\top }\mathbf{K}\mathbf{y}-ϵ\Vert \mathbf{x}{\Vert }^2$$

$$=-\sum _{e=1}^n{\tilde{\rho }}_e^{p}E({\delta }_e)\left\{{\mathbf{v}}^{\top }{\hat{\mathbf{K}}}_e\mathbf{v}-2{\mathbf{v}}^{\top }{\hat{\mathbf{K}}}_e\mathbf{y}+{\mathbf{y}}^{\top }{\hat{\mathbf{K}}}_e\mathbf{y}\right\}$$

$$-ϵ\Vert \mathbf{x}{\Vert }^2$$

$$=-\sum _{e=1}^n{\tilde{\rho }}_e^{p}E({\delta }_e)\left\{{\left[\mathbf{y}-\mathbf{v}\right]}^{\top }{\hat{\mathbf{K}}}_e\left[\mathbf{y}-\mathbf{v}\right]\right\}-ϵ\Vert \mathbf{x}{\Vert }^2$$

$$\{\begin{array}{ll}\le 0& ϵ=0\\ <0& ϵ>0\end{array}\forall \left(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right)\in {\mathbb{R}}^{n+m}$$

$$F_{ϵ}(\mathbf{\rho }):=\underset{\mathbf{\delta }\in {\Delta }_{\text{ad}},\mathbf{u}\in {\mathbb{R}}^m}{\max }{\mathcal{J}}_{ϵ}(\mathbf{\rho },\mathbf{\delta },\mathbf{u})$$

$$\underset{\mathbf{\delta }\in {\Delta }_{\text{ad}},\mathbf{u}\in {\mathbb{R}}^m}{\max }{\mathcal{J}}_{ϵ}(\mathbf{\rho },\mathbf{\delta },\mathbf{u})$$

$$E_e^{\text{eff}}=E_e({\delta }_e)=E_D+\frac{1-{\delta }_e}{1+q{\delta }_e}(E_0-E_D).$$

Full text

1 Introduction

In the last decade, structural optimization has seen a vast increase in publications considering uncertainties in the optimization problem. Uncertainties in structural optimization can arise from a multitude of different sources, such as modelling errors, inaccurate loading scenarios, manufacturing precision or defective material. In this work, we will focus on the latter case for compliance problems in linear elasticity.

Local material failure can be caused by unexpectedly high loads, leading to yielding and a local loss in stiffness. Other sources of local material degradation may be corrosion, fatigue failure or accidents. Furthermore, the fabrication process may lead to variances in the achieved material properties. In the linear elasticity setting, all of these material defects may be modeled by either a reduced material stiffness or a removal of material accounting for a complete failure of material.

In the literature, a number of publications can be found on topology optimization under consideration of uncertainties in geometry or material properties. Topological uncertainties are considered in Wang et al (2011) for an etching fabrication process through uniform over- and underetching. In Chen et al (2010), uncertainties in the loads or in the material parameters are represented with random fields using a Karhunen-Loeve expansion for a level-set method. In Lazarov et al (2012), the material properties and the topology are subject to random field uncertainties using a collocation method for density-based topology optimization. Both of the latter methods are able to yield robustly optimized structures, but suffer from a seriously increasing computational complexity if the distribution of the uncertain degradation is allowed to vary not only on a broad scale over the body. Furthermore, even if the uncertainty distribution is sufficiently well known and can be represented properly using the random field uncertainties, the methods may only give a probability for the actual structural performance. In this work, on the other hand, we would like to obtain an actual bound for the structural performance, even if the allowed uncertainties in the material properties are distributed in the worst possible way. One method to obtain worst-case structures accounting for total material failure is given in Jansen et al (2014). Here, the design domain is subdivided into finite element (FE) patches and, in each iteration, the material is removed for the FE patch in a worst-case manner. However, one FE simulation is needed per considered FE patch, so a fine resolution of material failure or allowing for the failure of a combination of multiple of the considered FE patches would lead to a prohibitively high computational complexity. In Bendsøe and Díaz (1998), a topology design problem in a continuum with worst-case material defects is solved using a relaxation to a two-scale problem with material mixtures of damaged and undamaged material. However, this is not done without the authors stating that a mathematically consistent model was not pursued. More publications studying worst-case optimization models can be found for truss optimization problems. In Calafiore and Dabbene (2008), robust truss optimization problems for uncertainties of polytopic shape in load and Young’s modulus are reformulated as convex semi-definite programs. The complexity of the reformulated problem, however, increases drastically with the number of vertices of the uncertainty polytope. The method thus only allows for the solution of problems with relatively low-dimensional uncertainty characteristics. In Achtziger and Bendsøe (1999), a bilevel truss optimization is discussed for minimum compliance with worst case material degradation. In fact, it is shown that the robust truss topology stays the same for the used material degradation model and pure cross-sectional bar area optimization without upper bounds.

Contrary to the discussed publications, the aim of this work is to provide a mathematically rigorous model which is capable of solving continuum topology optimization problems with a worst-case approach for uncertain material properties, covering local failure as well as manufacturing uncertainties arising, e.g., from powderbed based additive manufacturing technologies.

The structure of the paper is as follows: In section 2 we develop a minimax problem, which is well defined in the sense that the inner maximization problem can be solved to global optimality. In section 3 differentiability properties for the minimax problem are discussed. Moreover two regularization approaches are suggested. In both cases the marginal function associated with the minimax problem can be shown to be continuously differentiable. In section 4 a complete algorithm for the solution of the regularized minimax problem is established. Finally, in section 5 numerical examples are presented for which uncertainties arise from additive manufacturing as well as local damage.

2 Problem description

2.1 State problem

The worst-case problem is discussed solely on the discretized level, i.e. we assume that the design domain $\Omega\subset\mathbb{R}^{2}$ is discretized in $n$ elements $\Omega_{i}$, on which material properties and densities are constant. Furthermore we introduce a discrete convolution operator as follows: for an element $e$ we define the set $V(e)$ of neighbour elements whose distance is less than a given radius $R>0$, and $c_{e}$ the coordinates of the center of element $e$. Then, for a discretized material density distribution $\rho\in\mathbb{R}^{n}$, the filtered value on element $e$ becomes

[TABLE]

We abbreviate the filtered pseudo-density variable by $\bm{\tilde{\rho}}:=F*\bm{{\rho}}$. The state problem then takes the form

[TABLE]

where $\bm{f}$ is the numerically integrated load form, $\bm{u}$ the displacement vector and $\bm{K}$ is the global stiffness matrix assembled by the local stiffness matrices $\bm{K}=\sum_{e=1}^{n}\bm{K}_{e}$. The local stiffness matrices are

[TABLE]

where $\bm{B}_{il}$, $l=1,\ldots,n_{ig}$ denote standard strain displacement matrices and $n_{ig}$ the number of integration points. Finally, the precise form of the material property function $\bm{C}_{e}\left(\bm{\rho},\bm{\delta}\right)$, which depends on the pseudo-density $\bm{{\rho}}\in\mathbb{R}^{n}$ as well as on an uncertainty vector $\bm{{\delta}}\in\mathbb{R}^{n}$ is discussed in Section 2.2.

2.2 Material model

In additive manufacturing, not only the averaged material properties are anisotropic. The variations in the material properties discussed in this section can exhibit a different behaviour in different directions as well. Commonly, the variation of the elastic modulus in the build direction of a part is expected to be slightly larger than the variation within the build layer (see e. g. vdi (2015); Wegner and Witt (2012)). Nevertheless, we choose to use an isotropic model to keep our final model numerically manageable. The chosen formulation is based on a publication about truss topology design with degradation effects (Achtziger and Bendsøe (1999)), in which the degradation is considered as a decrease of material stiffness. For two elasticity tensors $\bm{C^{0}},\bm{C^{1}}$ the effective material property is computed by interpolation as

[TABLE]

with $\alpha\in[0,1]$. This simplified model was used for an elastic continuum by Francfort and Marigo (1993) and can be traced back to Kachanov (1958) and Lemaitre (1978). Achtziger and Bendsøe first used this model for damage in an optimization problem for maximum compliance in a truss structure in Achtziger and Bendsøe (1995). They revisited the formulation in Achtziger et al (1998) and were able to obtain lower and upper damage bounds by studying inner and outer approximations. This same inverse interpolation was used in Achtziger and Bendsøe (1999) for a basic truss topology design problem. Applied to a continuum setting, it reads

[TABLE]

where $e$ denotes the finite element number, $E_{0}>E_{D}>0$ the largest and smallest Young’s moduli allowed and

[TABLE]

a unit stiffness tensor with the Poisson’s ratio $\nu$ of the material. The parameter ${\delta}_{e}$ is used to interpolate between the least stiff and the stiffest material and is constant on each finite element $e$. This parameter represents the uncertainty in the material properties. A major advantage of this formulation will be the concavity of the compliance functional with respect to $\bm{\delta}$ in the optimization problem introduced in the next subsection. Furthermore, the material stiffness is strictly decreasing in ${\delta}_{e}$, so ${\delta}_{e}$ can be interpreted as a material degradation parameter.

As usual in the classic SIMP approach, the topology itself is represented through a pseudo-density variable $\bm{\rho}$. Consequently, the formula for the material property function becomes

[TABLE]

where $\hat{\bm{C}}_{e}$ is defined as in (5), $\tilde{\rho}$ denotes the pseudo-density after application of the convolution and $p>1$ is a penalty parameter.

2.3 Optimization problem

We start out with the formulation of a minimax problem, where the compliance is minimized w. r. t. the topology $\bm{\rho}$ and maximized w. r. t. the material uncertainties $\bm{\delta}$:

[TABLE]

Here, $U_{\text{ad}}$ and $\Delta_{\text{ad}}$ are the admissible design sets containing bound and volume constraints, $\bm{K}(\bm{\rho},\bm{\delta})$ the symmetric and positive definite stiffness matrix, $\bm{u}$ the vector containing the nodal displacements and $\bm{f}$ the external load vector. The stiffness matrix $\bm{K}$ is assembled analogously to (2) with $C_{e}$ as in (7), if not stated differently. The admissible set of pseudo-densities is given as

[TABLE]

where $\rho_{\text{min}}$ is a small positive number, $v_{e}$ is the volume of element $\Omega_{e}$ and $V$ is an upper bound on the admissible volume. Further, we use the uncertainty sets

[TABLE]

or alternatively

[TABLE]

with $D$ being an upper bound on the total material uncertainty. Note that for a given binary distribution $\bm{\tilde{\rho}}$ the two definitions of $\Delta_{\text{ad}}$ will yield similar solutions of the inner maximization problem, as material deviations will naturally be concentrated on elements where $\tilde{\rho}_{e}=1$. The first definition avoids the dependence on the upper level optimization variable $\bm{\rho}$, while the second definition is more rigorous from an application point of view.

Remark 1

In the application context of additive manufacturing (see section 5.2), we will further study the impact of additional uncertainty sets of the following form:

[TABLE]

or similarly

[TABLE]

For an interpretation of these sets we refer to section 5.2. From a theoretical point of view, for fixed $\bm{\rho}$, either of the sets is compact and convex and thus no special attention is required throughout this section.

Now, in order to analyse the structure of the minimax problem we rephrase problem (10) equivalently as:

[TABLE]

Along with (16) we study the problem

[TABLE]

where

[TABLE]

is a regularized variant of the objective function with $\epsilon>0$ being a small positive number causing the last term to act as a Tikhonov type regularization. We are now able to state the following important result.

Theorem 2.1

The functional $\mathcal{J}_{\epsilon}$ is concave w. r. t. $(\bm{\delta},\bm{u})$ for $\epsilon=0$ and strictly concave for $\epsilon>0$.

Proof

For the convenience of the reader, we first introduce the local stiffness matrices for a unit material tensor from (6) as

[TABLE]

With $E(\delta):=\left(\frac{1-\delta}{E_{0}}+\frac{\delta}{E_{D}}\right)^{-1}$ we then have

[TABLE]

and for $i,j=1,\ldots n$ it follows from straightforward calculations that

[TABLE]

Now we are able to show the (strict) concavity by demonstrating the negative (semi-)definiteness of the Hessian of $\mathcal{J}_{\epsilon}$:

[TABLE]

where $\bm{v}=x_{e}E(\delta_{e})\left(\frac{1}{E_{D}}-\frac{1}{E_{0}}\right)\bm{u}$. The latter inequality holds because of the local stiffness matrices $\hat{\bm{K}}_{e}$ being positive semi-definite, $\tilde{\rho}_{e}^{p}E(\delta_{e})>0$ and $\sum_{e=1}^{n}\tilde{\rho}_{e}^{p}E(\delta_{e})\hat{\bm{K}}_{e}$ positive definite.

A direct consequence of Theorem 2.1 is that the inner maximization problem in (10) has – for either of the choices of the uncertainty set outlined above – a unique solution for any feasible choice of of the pseudo-density $\bm{\rho}$ provided that the regularization parameter $\epsilon$ is positive. Moreover the computation of this maximizer is numerically tractable. In other words, the marginal function

[TABLE]

is well defined and can be efficiently evaluated. Moreover, as a maximum over infinitely many linear functions $F_{epsilon}$ is even continuous w.r.t. the design variable $\bm{\rho}$.

2.4 Discussion of the inverse interpolation model

As outlined previously the inverse interpolation formula (5) has the strong advantage that the resulting adversarial problem

[TABLE]

is concave and can thus be solved to global optimality with a reasonable effort. On the other hand, from an application point of view using the standard interpolation model (4) appears to be much more natural.

Motivated by this, throughout this section, we would like to discuss the impact of using the interpolation model (5) as an approximation of the linear interpolation model (4) as well as a connection of the two models through a continuation scheme for the uncertainty variable $\bm{\delta}$.

We start with a closer look at formula (5). As we have seen in the previous section, the properties of (5) were essential for the compliance to be concave in $\bm{\delta}$. In fact, this was one of the goals Stolpe and Svanberg wanted to achieve with the suggestion of an alternative penalization model in topology optimization, see Stolpe and Svanberg (2001). This interpolation scheme later became better known by the acronym RAMP for for Rational Approximation of Material Properties. Using the RAMP model to prescribe the effect of material uncertainty on the effective material stiffness in an $e$, we may write:

[TABLE]

Here, the parameter $q\geq 0$ serves is the RAMP penalization parameter. It follows directly from Stolpe and Svanberg (2001) that using this interpolation scheme in a linear elastic framework, the compliance is concave for all $q\leq(E_{0}-E_{D})/E_{D}$. Moreover, it can be easily checked that in the case $q=(E_{0}-E_{D})/E_{D}$ the interpolation scheme (27) coincides with our inverse interpolation scheme (5), i.e. we can represent (5) in a RAMP fashion.

Examples of the effective material stiffness computed by (5) for different values of ${\delta}$ and $E_{D}$ are shown in Figure 1. It is seen that for an increasing difference between $E_{0}$ and $E_{D}$, the penalization parameter $q=(E_{0}-E_{D})/E_{D}$ leads to increasing penalization of intermediate values strictly between 0 and 1. However, as the compliance is maximized w.r.t. $\bm{\delta}$ (see (26)), this means rewarding intermediate values if a linear constraint is used to bound the total material uncertainty as in (12) or (13).

Now, in order to compare the effect of the linear and the inverse interpolation scheme using the same linear total uncertainty bound, we can sketch the effective material properties obtained in both cases. Figure 2 demonstrates this for a two element setup with a total uncertainty bound $\delta_{1}+\delta_{2}\leq 1$ and material data $E_{0}=1$ and $E_{D}=0.7$.

It is clearly seen that the admissible set obtained using the inverse interpolation model (magenta set) is a relatively tight outer approximation of the admissible set obtained using the linear interpolation model (blue set) for this choice of material parameters (see Figure 2d). This is true in general, if the material deviations characterized by the difference $(E_{0}-E_{D})/E_{0}$ are sufficiently close to 1. If this is the case, applying model (5) as an outer approximation of model (4) in the framework of the robust optimization model (10) will lead to only slightly more conservative topological solutions. This appears to be a small price, taking the numerical tractability into account. It is noted that $(E_{0}-E_{D})/E_{0}$ is comparably large all additive manufacturing scenarios discussed in section 5.

On the other hand Figure 3 reveals that for smaller and smaller $E_{D}$, it becomes more and more likely that the worst-case effective material (i.e. the maximizer of the compliance function), is close to the black spots highlighted in Figure 3. If this happens, the solution of the adversarial problem (26) can be quite far away from the original uncertainty sets (blue sets in Figure 3) and thus lead to an over-conservative solution. As a consequence, from an application point of view it might be more attractive to work directly with the linear interpolation model. As the latter can not be expected to be solved to global optimality for a meaningful size of the problem, we suggest the following heuristic based on the RAMP interpretation of model (5): we start with RAMP parameter $q=(E_{0}-E_{D})/E_{D}$ – corresponding precisely to the concave variant of the adversarial problem – and from there slowly decrease the RAMP parameter towards 0, corresponding precisely to the linear interpolation rule. It should be stressed that there is no guarantee that this procedure leads to the global maximizer of (26) with linear interpolation model, however it provides a lower bound for the worst case compliance of a given topology. As the inverse model at the same time provides a rigorous upper bound for the same, at least this information can be used to estimate how conservative a computed robust solution is. Moreover numerical experiments in section 5.3 will reveal that much more reasonable results are obtained compared to a direct application of the uncertainty model (4).

We finally note that one may also think about inner approximations (green sets in Figure 3) instead of outer approximations, but then the true worst-case might be missed, when solving the adversarial problem and thus the term worst-case looses its meaning. This is why we did not further pursue this direction in this paper.

3 Obtaining derivative information

In Section 2.3 we have seen that the marginal function $F_{\epsilon}(\bm{\rho})$ is well-defined and continuous. Taking additionally the compactness of the set $U_{\text{ad}}$ (see (11)) into account it is clear that the problem

[TABLE]

which is fully equivalent to the regularized problem (17), has a solution. As the objective $F_{\epsilon}$ depends only on $\bm{\rho}$ in an explicit manner, we may consider using standard, gradient-based optimzation solvers for the numerical solution of (28). This is our motivation to analyse differentiability properties of $F_{\epsilon}$ in the sequel.

3.1 The $\epsilon$-regularized minimax approach

In order to investigate the differentiability properties of $F_{\epsilon}$ we make use of results from literature formulated for marginal functions of the type

[TABLE]

where $\mathpzc{j}:\mathbb{R}^{l}\times\mathbb{R}^{k}\to\bm{\mathbb{R}}$ and $Y(\cdot)$ is a point-to-set map from $\mathbb{R}^{k}$ to $\mathbb{R}^{l}$. Moreover we assume that the set $Y(x)$ (which is in the context of our paper the uncertainty set) is given in the following form:

[TABLE]

where $G:\mathbb{R}^{l}\times\mathbb{R}^{k}\to\mathbb{R}^{s},(x,y)\mapsto(g_{1}(x,y),\ldots,g_{k}(x,y))$ and $H:\mathbb{R}^{l}\times\mathbb{R}^{k}\to\mathbb{R}^{q},(x,y)\mapsto(h_{1}(x,y),\ldots,h_{q}(x,y))$ are continuously differentiable in both variables.

Next we state two well known constraint qualifications for the maximization problem in (29).

Definition 1 (Mangasarian-Fromowitz constraint

qualification)

A feasible point $\bm{\bar{y}}\in Y(\bm{\bar{x}})$ is said to satisfy the Mangasarian-Fromowitz constraint qualification (MFCQ) if the following holds:

There exists a direction $\bm{d}\in\mathbb{R}^{k}$ such that $\nabla_{\bm{y}}H(\bm{\bar{x}},\bm{\bar{y}})^{\top}\bm{d}=\bm{0}$ and $\nabla_{\bm{y}}g_{i}(\bm{\bar{x}},\bm{\bar{y}})^{\top}\bm{d}<0$ for all $i\in I(\bm{\bar{x}},\bm{\bar{y}})=\left\{i\in\{1,\ldots,s\}\mid g_{i}(\bm{\bar{x}},\bm{\bar{y}})=0\right\}$. 2. 2.

The gradients of the equality constraints $\nabla_{\bm{y}}h_{j},j=1,\ldots,q$ are linearly independent at $(\bm{\bar{x}},\bm{\bar{y}})$.

Definition 2 (Linear independence constraint qualification)

A feasible point $\bm{\bar{y}}\in Y(\bm{\bar{x}})$ is said to satisfy the linear independence constraint qualification (LICQ), if the system of gradients $\nabla_{\bm{y}}h_{i}(\bm{\bar{x}},\bm{\bar{y}}),\ i=1,\ldots,q$, $\nabla_{\bm{y}}g_{j}(\bm{\bar{x}},\bm{\bar{y}}),\ j\in I(\bm{\bar{x}},\bm{\bar{y}})$ is linearly independent.

Note that the LICQ is stronger and in fact implies MFCQ. Furthermore it is well known that the LICQ immediately guarantees the uniqueness of Lagrangian multipliers associated with the constraints $G(\bm{x},\bm{y})\leq\bm{0},H(\bm{x},\bm{y})=\bm{0}$. Moreover we need the following definition.

Definition 3

The set $Y(\bar{x})$ is called uniformly compact near $\bm{\bar{x}}$ if there is a neighborhood $N(\bm{\bar{x}})$ of $\bm{\bar{x}}$, such that the closure of $\bigcup_{\bm{x}\in N(\bm{\bar{x}})}Y(\bm{x})$ is compact.

After these preliminaries, we can state the following result taken from Gauvin and Dubeau (1982):

Theorem 3.1

If $Y(\bm{\bar{x}})$ is nonempty, $Y(\bm{\bar{x}})$ is uniformly compact near $\bm{\bar{x}}$ and the MFCQ is satisfied at every point $\bm{\bar{y}}\in M(\bm{\bar{x}})$, then the following inclusion holds true for the generalized gradient in the sense of Clarke:

[TABLE]

Here $M(\bm{\bar{x}})$ is the set of maximizers for a given $\bm{\bar{x}}$, i. e.

[TABLE]

and $K(\bm{\bar{x}},\bm{\bar{y}})$ is the set of optimal Lagrange multipliers

[TABLE]

If, in addition, the LICQ holds at every point $\bm{\bar{y}}\in M(\bm{\bar{x}})$

[TABLE]

where $(\bm{\bar{\lambda}},\bm{\bar{\mu}})$ are the unique multipliers from $K(\bm{\bar{x}},\bm{\bar{y}})$.

Finally, the following theorem is taken from Clarke (1975):

Theorem 3.2

The following are equivalent:

$\partial\psi(\bm{x})=\left\{\bm{\zeta}\right\}$, i. e. the generalized gradient reduces to a singleton. 2. 2.

$\nabla\psi(\bm{x})$* exists, $\nabla\psi(\bm{x})=\bm{\zeta}$, and $\nabla\psi$ is continuous at $\bm{x}$ relative to the set on which it is defined.*

Based on these general results, we can state the differentiability result for $F_{\epsilon}$:

Theorem 3.3

Let

[TABLE]

with $\epsilon>0$ and $\Delta_{\text{ad}}(\bm{\rho})$ is given by formula (13). Let further $(\bm{\bar{\delta}},\bm{\bar{u}})$ denote the unique maximizer with $F_{\epsilon}(\bm{\bar{\rho}})=\mathcal{J}_{\epsilon}(\bm{\bar{\rho}},\bm{\bar{\delta}},\bm{\bar{u}})$ for $\bm{\bar{\rho}}\in U_{\text{ad}}$ and assume that the LICQ condition is satisfied at every point $(\bm{\bar{\delta}},\bm{\bar{u}})\in M(\bm{\bar{\rho}})$.

Then the following holds: $F_{\epsilon}$ is continuously differentiable at $\bm{\bar{\rho}}$ and

[TABLE]

Remark 2

Although we can not prove that LICQ holds in general, we would like to remark that LICQ is expected to hold in any practically relevant situation. To see this, we note that the only way that LICQ will not hold at some point $\bm{\bar{\rho}}$ is that the worst-case solution $\bm{\bar{\delta}}$ is fully binary and $\sum_{e:\bar{\delta}_{e}=1}v_{e}\tilde{\bar{\rho}}_{e}/|\Omega|=D$. This is very unlikely for any $D<|\Omega|$. Furthermore, as $\bm{\bar{\delta}}$ is mostly large where $\bm{\tilde{\rho}}^{p}$ has value $1$, $D$ may be chosen such that $\sum_{e:\bar{\delta}_{e}=1}v_{e}/|\Omega|$ cannot equal $D$ (e. g. for a uniform finite element size s. t. $n\cdot D\notin\mathbb{N}$). Finally, in practice, no setting could be found which resulted in a 0-1 distribution of $\bm{\delta}$. We furthermore remark that analogous results can be proven, if instead of formula (13), one of the formulae (12), (14) or (15) are used.

Proof (Theorem 3.3)

We want use Theorem 3.1. We have $Y(\bm{\rho})={\Delta_{\text{ad}}(\bm{\rho})\times\mathbb{R}^{m}}$ and $\mathpzc{j}=\mathcal{J}_{\epsilon}(\bm{\rho},\bm{\delta},\bm{u})$. While $Y(\bm{\rho})$ is nonempty and $\Delta_{\text{ad}}(\bm{\rho})$ is uniformly compact, $\mathbb{R}^{m}$ is evidently unbounded. However, as in the optimum we have $\bm{u}=\bm{K}^{-1}(\bm{\rho},\bm{\delta})\bm{f}$ with $\bm{K}$ uniformly positive definite, the optimal $\bm{u}$ is indeed uniformly bounded. Thus, we can modify the admissible set of $\bm{u}$ to include artificial lower and upper bounds without changing the solution of the optimization problem. Now, as we assume LICQ to hold, from Theorem 2.1 and Theorem 3.2 we directly obtain the continuous differentiability of $F_{\epsilon}$. Finally, the gradient formula (33) follows directly from the general formula (31).

Remark 3

In principal, we can now solve problem (29) by our favorite gradient based optimization solver, where in each iteration we solve a maximization problem to evaluate $F_{\epsilon}$ and apply formula (33) to evaluate the associated gradient. However, this will only work well in practice, if the inner maximization problem is solved with high accuracy in each iteration. Moreover, in order to obtain a good approximation of the gradient, also tight numerical approximations of the Lagrangian multipliers associated with the constraints of the inner problem have to be accessible. Taking these observations as well as the sparsity structure of the inner problem into account it is advisable to use a second order method for the evaluation of $F_{\epsilon}$. As a matter of fact, interior point methods belong to the most efficient second order optimization methods, in particular if the problem to be solved belongs to the class of convex minimization problems. As the maximization problem required to evaluate $F_{\epsilon}$ can be easily turned into a problem of this kind, we suggest to use an interior point method also here. Still the problem remains that the precise evaluation of $F_{\epsilon}$ can be numerically quite demanding. As a remedy, in the following section, we suggest an alternative regularization scheme which is more tightly connected to interior point methods and thus allows in general for a more efficient solution of the regularized inner maximization problem.

3.2 Regularization using a barrier method

In this section, we study the approach of using a barrier term for the regularization of the minimax problem. This will allow to show differentiability without using the minimax theorems of the previous section and, more importantly, turn out to be superior w. r. t. computational efficiency.

We consider (16) as a specific bilevel optimization problem as follows:

[TABLE]

Furthermore, we write the uncertainty set in the generic form

[TABLE]

where

[TABLE]

for the set (12) and

[TABLE]

for (13). Note that again we may alternatively use the uncertainty sets (14) or (15) with only minor modifications.

The idea is now to replace the lower-level (or inner) optimization problem by its optimality conditions and solve the upper-level (or outer) optimization problem as a constrained problem with these optimality conditions given as constraints. We remark that for fixed ${\bm{\rho}}$ the uncertainty sets above are affine linear and thus no constraint qualification is required. Thus we may write the optimality conditions of the lower-level optimization problem as follows:

[TABLE]

However it turns out that the resulting problem is difficult to solve in practice mainly due to the presence of the complementarity constraints on the lower level. Thus, we try to get rid of those in the sequel. First of all we find

[TABLE]

Using this it is seen from (36) that the constraint for the total uncertainty will be always at the bound and can hence simply replaced by an equality constraint without changing the solution of the problem. This leaves us with the bound constraints on $\bm{\delta}$, which are simply approximated using appropriate barrier terms in the objective. We obtain the problem

[TABLE]

where $\mu^{*}$ is a small positive barrier parameter. Even though the application of this barrier method does indeed change the solution of the optimization problem, the influence will be minor if $\mu^{*}$ is chosen small enough. Moreover, the additional term has another benefit. Using similar arguments as in the proof of Theorem 2.1, we can show for the regularized function

[TABLE]

that

[TABLE]

for all $\left(\bm{x}^{\top},\bm{y}^{\top}\right)\neq\bm{0}$ and thus strict concavity of $\mathcal{J}_{\mu^{*}}$ in the joint variable $(\bm{\delta},\bm{u})$ for every $\mu^{*}>0$.

This implies that the constraint in (LABEL:eq:problem-bilevel) can be equivalently stated through optimality conditions, which read as

[TABLE]

Now, having eliminated all complementarity constraints, we may formulate a new bilevel optimization problem as

[TABLE]

and its implicit counterpart as

[TABLE]

Here $\bm{\delta}(\bm{\rho})$ and $\bm{u}(\bm{\rho})$ are the solution operators associated with (44). Note that $\bm{\delta}(\bm{\rho})$ and $\bm{u}(\bm{\rho})$ are well defined as the system (44) has a unique solution for each feasible $\bm{\rho}$ thanks to the strict concavity of (3.2).

For the latter problem, we can now use a standard adjoint approach for the calculation of the gradient. In order to do so, we use the Lagrangian corresponding to the lower-level problem

[TABLE]

in order to rewrite the KKT conditions (44) as

[TABLE]

Now, using the adjoint calculus for problem (LABEL:eq:problem-nested2), we get:

[TABLE]

where $\left(\bm{\tau}_{\bm{\delta}},\bm{\tau}_{\bm{u}},\tau_{\lambda}\right)$ are the solutions of the adjoint equations

[TABLE]

Finally, if we slightly modify problem (LABEL:eq:problem-nested2) by replacing the upper level objective $J_{0}$ by $J_{\mu^{*}}$ and use again formula (47) and (48), the right hand side simplifies to

[TABLE]

With this modified right hand side, we can directly find $\left(\bm{\tau}_{\bm{\delta}},\bm{\tau}_{\bm{u}},\tau_{\lambda}\right)=\left(\bm{0},\bm{0},-\lambda\right)$ to be a solution of system (50). Thus, no adjoint problem has to be solved and we obtain the derivative formula:

[TABLE]

It is seen that this formula has precisely the same structure as the derivative formula for $F_{\epsilon}$ in Section 3.1. This is not too surprising, as, by changing the objective as above, (LABEL:eq:problem-nested2) is essentially turned into a minimax problem and thus the generic theory for the minimax-problems in Section 3.1 also applies to this situation.

4 Algorithmic solution

4.1 Lower-level problem solution

Even though we have seen that gradients of the upper-level problem can be easily computed using, e. g., the minimax theorems in Section 3.1, the lower-level problem still has to be solved very precisely in each outer iteration for a sufficiently accurate gradient. Using, for instance, the standard topology optimization algorithm MMA for the solution of the inner problem will hence lead to an very high number of total iterations required. On the other hand, the sparse Hessian of the lower-level maximization problem and its strict concavity lend themselves to the use of second-order optimization methods and their much better speed of convergence. For this, we take a quick look at interior-point methods.

Interior-point methods replace inequalities with exactly the same barrier term that was used in Section 3.2 and iteratively solve the resulting problems for a decreasing sequence of barrier parameters $\mu\searrow 0$. For the sake of completeness, We demonstrate this concept for a simple model problem with objective $f:\mathbb{R}^{n}\to\mathbb{R}$, and a constraint $g:\mathbb{R}^{n}\to\mathbb{R}$: The original problem

[TABLE]

is converted into a sequence of barrier problems parametrized by the barrier parameter $\mu>0$:

[TABLE]

Each of the latter problems can be now approached by solving the optimality conditions

[TABLE]

In order to get a full picture of the error made in comparison to the KKT system of the original optimization problem, the term $(\mu/x_{e})_{e=1}^{n}\geq\bm{0}$ is interpreted as a Lagrange vector $\bm{z}\geq\bm{0}$ of the inequality constraints. Substituting $\bm{z}:=(\mu/x_{e})_{e=1}^{n}$ in the above optimality conditions then gives the primal-dual equations

[TABLE]

Note that for $\mu=0$ these are the KKT conditions of the original optimization problem including the complementarity constraints. Hence, the last line in the equation set is sometimes called perturbed complementarity and $\mu>0$ is the desired complementarity error. This system, sometimes also referred to as central path equations is now solved repeatedly for decreasing $\mu$.

For the problem considered in this work, the barrier objective obtained by the interior point procedure is exactly the cost function $\mathcal{J}_{\mu}$, see (3.2). Hence, we can use existing interior point solvers and simply stop the decreasing sequence of $\mu$ at $\mu^{*}$. Moreover, this automatically implies a numerical advantage of the regularization with a barrier term compared with the $\epsilon$-regularization of Section 3.1: The solution of the lower-level problem is faster because the decreasing $\mu$ sequence can be stopped consistently sooner. Furthermore, we have seen that $\mu$ has a clear interpretation as complementarity error, while the influence of $\epsilon$ on optimization results is much less clear.

4.2 Algorithm

We combine two different optimization algorithms for the upper- and lower-level problems. We solve the lower-level optimization problem in $(\bm{\delta},\bm{u})$ with the second order interior point optimizer IPOPT, see Wächter and Biegler (2006), and then compute the gradients w. r. t. $\bm{\rho}$ using formula (51). For the upper-level problem in $\bm{\rho}$, we use the Method of Moving Asymptotes (MMA, see Svanberg (2002)). This procedure is repeated until convergence is achieved.

Note that line 1 is optional, however doing at least a couple of “cheaper” non-robust optimization steps will result in being closer to an optimal topology as a starting value and the convergence speed of the subsequent robust iterations can be improved. For the results in this work, the (pure) topology optimization is done until convergence, as it is used to obtain a reference compliance value in Section 5.

4.2.1 Software settings

As a default, IPOPT starts with a smoothing parameter of mu_init=0.1, which is subsequently decreased during the optimization process unto a final value of mu_target. The final smoothing parameter mu_target corresponds to $\mu^{*}$ in our formulation and we set it to mu_target=$10^{-7}$. Note that the choice of $\mu^{*}$ depends on the scaling of the inner cost function. If it is chosen too large compared to the cost function, the barrier term is large enough to push $\bm{\delta}$ significantly towards intermediate values, thus noticeably influencing the optimization results. If it is chosen too small, the intended strict concavity and regularization effects get less pronounced. More importantly, the number of iterations needed for solving the lower-level problem as well as the computation times decrease when choosing $\mu^{*}$ larger, as this determines how early IPOPT’s usual decrease of $\mu$ to 0 may be stopped. As the initial barrier parameter, in the first iteration we use the default value mu_init=0.1. After the first iteration, when we can reuse the previous iterates of $\bm{\delta}$, $\bm{u}$ and IPOPT’s Lagrangian parameters as a starting value, we start IPOPT directly with mu_init=$\mu^{*}$.

As termination criterion, IPOPT uses essentially the maximum norm of the KKT error for the central path equations associated with mu_target. The termination criterion is set to a very strict value of tol=$10^{-10}$. Furthermore, the terminal maximum constraint violation and complementarity error can be set seperately as well. We use the values of constr_viol_tol=$10^{-10}$ and compl_inf_tol=$0.0001$, respectively. The strict termination criteria ensure that the optimality assumptions for the computation of the upper-level gradients are fulfilled sufficiently accurate.

As the linear system solver within IPOPT, we use the HSL Mathematical Software Library HSL (2013).

5 Numerical results

5.1 Numerical example setup

As exemplary setup we use a loaded cantilever fixed on the left side with a rectangular design domain of the dimension $[0,2]\times[0,1]$ as depicted in Figure 4. The cantilever is loaded with a line load at the lower right bottom starting at $x=1.9$ of total magnitude $0.3$ in direction $-y$. The domain is discretized using $300\times 150$ uniform rectangular finite elements. We use a linearly decaying density filter with a radius of $0.045$ corresponding to the length of $4.5$ elements. As resource constraints, we use a total relative material volume of 50% ($V=0.5$). The penalty parameter is set to $p=4$ and the lower bound of the topology variable to $\rho_{\text{min}}=0.01$. The maximum number of iterations is $500$.

The starting value for the robust optimization is for all results the solution of the respective SIMP topology optimization problem without material uncertainties, i.e. ($\bm{\delta}\equiv 0$). This reference problem and the resulting topology will be denoted as nominal problem and nominal topology, respectively. For this nominal topology, we compute the compliance using the material with Young’s modulus $E_{0}$ as a nominal compliance. This nominal “best-case compliance” is then compared to the worst-case compliance of the nominal topology, as well as the best- and worst-case compliances of the robustly optimized topology.

5.2 Additive manufacturing

The complex process in additive manufacturing can lead to an inhomogeneous microstructure of the material, e.g. by microscopic pores. This introduces variances in the local effective material properties. This variance will be modeled as a uncertainty in the Young’s modulus, which can also be observed in measurements. The measured maximum and minimum Young’s modulus will then be used in the worst-case framework detailed in this work. Specifically, the uncertain material properties are allowed to deviate by 30%, i. e. the normalized values in the optimization problem are $E_{D}=0.7$ and $E_{0}=1$. These values are based on tensile testing results for AlSi10Mg (vdi, 2015). Note that the variations in Young’s modulus will change for different processes, materials, machines and process parameters. E. g. for reused polyamide powder, the elastic modulus may decrease rapidly up to half the original value after several uses (Wegner and Witt, 2012) while for the maraging steel 18Ni-300 the variations may be closer to 14% (see Hermann Becker et al (2016)). As the amount of material variations occurring is even harder to estimate, we use a number of different volumes for the degeneration variable ($D=0.03,\,D=0,06$ and $D=0.10$). However, we restrict ourselves in the visualization to a single representative.

Before stating the results, we can estimate the maximum and minimum effect that these numbers can theoretically have on the compliance. As a change in material has a linear influence on the compliance, a decrease in material stiffness from $1$ to $0.7$ on the whole domain would amount to $1/E_{D}=1/0.7\approx 143\%$ times the compliance compared to using only the stiffest material. Similarly, we can distribute the fraction of degenerated material uniformly on all elements to get a lower bound of $((1-D)/E0+D/E_{D})$ for $\bm{\delta}\equiv D$. However, as we can assume an approximate 0-1-distribution for the topology, a better estimate would be to distribute the degeneration variable only on the elements with $\bm{\rho}=1$. Neglecting that the topology variable is not actually binary, this gives us $\bm{\delta}\equiv D/V$ and a (slightly too large) lower bound estimate of the minimum effect of the degradation of $((1-D/V)/E0+(D/V)/E_{D})$. For the total uncertainty bounds of $D=0.03,\,D=0,06$ and $D=0.1$, this gives a compliance increase of about $2.57\%,5.14\%$ and $8.57\%$, respectively.

The compliance values obtained from numerical results can be found in Table 1 and are denoted as the increase w. r. t. the compliance of the nominal topology (standard SIMP result) evaluated with material $E_{0}$ only. As expected, the worst-case compliance values greatly surpass those of the lower bound estimates. Furthermore, the compliance of the robust topology evaluated with $\bm{\delta}\equiv 0$ slightly deteriorates, while the worst-case compliance improves about twice as much, if we pass from the nominal to the robustly optimized topology. Looking at the visualization of some results in Figure 5(a), we first observe that the robust optimization does not change the nominal SIMP topology a lot. Only on the left side we see that the structure gets slightly thicker along the support. It also becomes apparent that the worst case distribution of the degradation variable $\bm{\delta}$ is almost binary.

Next we test the impact of changing the uncertainty set. The motivation for this is that a worst-case behaviour in the sense that the material properties are systematically weakened is not really expected in additive manufacturing practice. Looking e. g. at vdi (2015), we rather expect the Young’s modulus to fluctuate around a mean value. An uncertainty model reflecting this behaviour has already been defined in Section 2 and is repeated here for the convenience of the reader:

[TABLE]

This choice of $\Delta_{\text{ad}}$ causes that the material defects to fluctuate around an expected value between the stiffest and the least stiff material. Note that by the constraint $(v_{e}\bm{\tilde{\rho}}_{e}^{p}\delta_{e})/|\Omega|=0.4$, for $E_{0}=1$ and $E_{D}=0.7$ we get an average Young’s modulus of $(0.6/E_{0}+0.4/E_{D})^{-1}\approx 0.8537$, which is very close to the true mean value of $0.85$. Furthermore, the quadratic constraint $\sum_{e=1}^{n}(v_{e}(\bm{\tilde{\rho}}_{e}^{p}\delta_{e}-0.4)^{2})/|\Omega|\leq D$ models that values of $\bm{\delta}$ close to the average of $0.4$ (and accordingly Young’s moduli close to the expected value) are comparatively cheap, while values further away are significantly more expensive.

For this setting, we use the compliance of the nominal topology with the average material generated by $\bm{\delta}\equiv 0.4$ on the whole domain as reference compliance corresponding to the mean stiffness expected from tensile testing. The worst-case compliance is expected to be much closer to the reference value now, as the allowed variation from the reference is halved. In addition the choice of the uncertainty set implies that weak material at one location has to be compensated in equal amounts by stiffer material in other places. We consider different bounds on the total material variations defined now by the quadratic constraint. Specifically we use $D\in\{0.01,0.02,0.03,0.04,0.06,0.1\}$. The relative compliance values obtained by numerical experiments for this setting can be found in Table 2; selected results are visualized in Figure 6. We note that the compliance values are now rather close to each other. For $\bm{\delta}\equiv 0.4$ the robust topology is only slightly worse than the nominal topology while still consistently outperforming the latter by about twice that margin for the worst-case material distribution. Looking at the images, the worst-case material distribution looks a lot more consistent with the continuous variations of the Young’s modulus that we see in tensile testing. While the worst-case material distribution still shows material with the lowest possible stiffness in the most critical areas of the cantilever, now material stiffnesses in the whole allowed range can be found, especially for the lower total uncertainty bounds, e.g. $D=0.003$.

In summary, the presented results suggest that in additive manufacturing the influence of the uncertain material parameters on the compliance of topology optimized parts is only minor. Nevertheless, the robust optimization methods are consistently decreasing the influence of the worst-case material distribution.

5.3 Robust topology optimization with material degeneration

As for the additive manufacturing scenarios, only a minor effect of material uncertainties was observed, a second scenario is studied, which allows more easily to investigate the strength of our robust model and algorithm. Specifically, we will use a significantly lower material parameter $E_{D}$ modelling that the material is heavily degraded at local spots. As this situation clearly corresponds to a degradation scenario, we apply the uncertainty model (13) again. We start our investigations with values $E_{0}=1$ and $E_{D}=0.01$. For these we consider a large number of different total uncertainty bounds ranging from $D=0.0002$ to $D=0.008$, as the robust topology and compliances are heavily influenced by this bound.

Selected optimal results are visualized in Figure 7(a) and the relative compliance values can be found in Table 3. As in the case of additive manufacturing, the cantilever gets thicker at the support, but this time by a much larger margin. It can be seen quite nicely, how – with increasing material uncertainty – more and more of the support is occupied with “degenarate” material. The compliance values in Table 3 reveal that now the worst-case compliance is decreased immensely for the robust topology compared to the nominal one, while sacrificing comparatively little structural performance for the case without material degeneration. For the nominal SIMP topology, degradation is almost exclusively observed at the support, whereas for the robustly optimized topology degradation starts appearing at the load already for $D=0.0002$ and is getting more pronounced for higher values of $D$.

At this point we would like to recall that – as outlined in Section 2.4 – in the degeneration case, i.e. $E_{D}\ll E_{0}$, the difference between the inverse interpolation model (5) and the linear interpolation model (4) may become substantial. In order to investigate this in practice, we would first like to get a sense of what results for the linear interpolation would look like. Simply replacing model (5) by (4) in the lower-level problem however leads to a non-concave maximization problem, for which it is very hard to find a result close to a global maximizer. In order to facilitate finding at least a good local maximum, we lend from the idea of a continuation of the penalty parameter in standard topology optimization using the interpretation of the inverse interpolation model as the RAMP scheme (27) with $q=(E_{0}-E_{D})/E_{D}$. We start by solving the same problem as before with $q=(E_{0}-E_{D})/E_{D}$ and continuously decrease $q$ unto a value of 0. Note that the lower-level problem fails to be strictly convex for $q<(E_{0}-E_{D})/E_{D}$, but that we are still able to obtain locally optimal solutions using IPOPT. At convergence of the continuation scheme, i. e. for $q=0$, the interpolation between $E_{D}$ and $E_{0}$ becomes linear. It can be easily checked that the objective becomes convex in this case. As we are maximizing the objective in the lower level problem, this means that have to expect of a worst-case material distribution which is close to binary. On the other hand, for $\bm{\delta}\in\{0,1\}^{n}$, the inverse interpolation model yields the same effective material properties as the linear interpolation model and thus the same compliance. This means, the less intermediate values of $\bm{\delta}$ appear in the worst-case material distribution, when we apply the inverse interpolation model, the closer the gap to the compliance obtained with the linear model is expected to be. Finally, we would like to recall that – for each fixed topology $\bm{\rho}$ – the worst-case compliance computed using the inverse model is an upper bound for the true worst-case compliance (in the sense of a global maximum) using the linear model. Moreover the approximate solution computing using the RAMP continuation scheme constitutes a lower bound for the same value.

We can now apply the RAMP scheme in two different ways. We either can use it to check the conservatism of the result obtained from the solution of the minimax problem based on the inverse interpolation model. In that case, the RAMP scheme is applied only once for the final robustly optimized topology. If the gap between lower and upper bound for the worst case compliance turns out to be small, we can be sure that we did not waste much performance by the effective relaxation of the uncertainty set (cf. again Figure 2 for an illustration). The second possibility is to directly perform the optimization based on the RAMP approach, i. e. in each iteration we essentially attempt to solve the adversarial problem formulated based on the linear interpolation model. It has to be mentioned that – from a theoretical point of view – this is not rigorous as the marginal function might by even discontinuous whenever we fail to find the global maximum of the inner problem. Nevertheless, we follow this second route in the remainder of this section.

Technically, after trying a number of different continuation schemes, we settled on 10 continuation steps from the inverse to the linear model as a reasonable compromise between a smooth continuation and computational performance. We performed numerical experiments with material parameters $E_{0}=1$ and $E_{D}=.01$ as well as total uncertainty bounds varying from $0.0005$ to $0.02$. It is interesting to note that we did not observe any technical difficulties which could be traced back to a potential discontinuity of the marginal function during numerical computations. The results in terms of relative compliances are summarized in Table 4. The associated optimal topologies and worst-case results are illustrated in Figure 8(a). In Table 4, for either of the worst-case computations we state three different compliance values: The worst-case compliance obtained using the inverse model (which serves as an upper bound here), the result obtained using the RAMP continuation scheme and the result obtained from attempting to solve the worst-case problem in a brute force way for the linear model. While we observe that there is only one exceptional case ($D=.001$ for the nominal topology), for which the worst-case compliance obtained by the direct application of the linear model is better than the one obtained by the continuation scheme, in all other cases the continuation scheme turns out to be by far superior. It should be noted that for the direct solution of the optimization problem with linear interpolation and larger uncertainty bounds IPOPT sometimes struggled to obtain a solution with values of $\bm{\delta}$ close to 0, explaining the extreme difference in comparison with the continuation scheme. It is also worth to note that for the nominal topology for which a close to binary worst-case material distribution has been already obtained based on the inverse model, the compliance calculated by the continuation scheme comes very close to the upper bound, which indicates that we found a quite good approximation of the true worst case also in the linear interpolation model. Even more surprising is that the situation for the final worst topologies is not much different: also here the gap between lower and upper bounds becomes surprisingly tight. The corresponding compliance values differ for larger choices of $D$ by only 10-20%. Besides this, the worst-case compliance values shown in Table 4 again look rather impressive: in almost all cases the compliance loss can be consistently more than halved for the robust topology compared with the nominal one, while again the compliance of the worst case topology is only little larger than the nominal topology for the best case material parameters corresponding to $\bm{\delta}\equiv 0$.

Looking at Figure 8(a) we first see again that for an increasing bound on total material degradation the robust topology becomes thicker along the support and more and more degraded material appears where the force is applied. On the other hand, in sharp contrast to the results obtained using the inverse interpolation model, now also for the robustly optimized topologies a binary degradation pattern is observed. As outlined above, this is expected due to the convex nature of the maximization problem using the linear interpolation scheme.

5.4 Computational efficiency

All computations were done on a machine with an Intel Xeon CPU E3-1280 v2 with four cores with 3.60 GHz and 32 GB of memory. Despite the very strict optimality and constraint tolerance of $10^{-10}$ used in IPOPT, the second order algorithm and the reuse of the previous iterates brought down the iteration numbers for each solution of the lower-level problem to a reasonable range. An individual lower level iteration took about one second for $300\times 150$ finite elements. For the additive manufacturing example with the average-based uncertainty set, the first solution of the lower-level problem took 21-22 interior point iterations, each solution after that only 4-5 and occasionally 6 interior point iterations. For the material degradation examples with $E_{D}=0.01$, the first outer iteration needed between 13 and 59 lower-level iterations, which decreased to no more than 10 iterations after the first 3-8 outer iteration steps. For the RAMP continuation, IPOPT struggled to solve the lower-level problem as efficiently. Without a noticeable change between outer iterations, the number for the RAMP continuation increased with growing total uncertainty bound from about 300 iterations ($D=0.0005$) to 1150 ($D=0.02$). It should be mentioned however that the warm start strategy was much less successful in this case, leaving room for improvements in the future.

For larger problem sizes, the iteration numbers and computation times are expected to scale quite well, as the main complexity comes from the solution of a strictly convex problem with a second-order method. We tested this for the examples with $E_{D}=0.01$ on a grid of $600\times 300$ finite elements and found a quite consistent scaling for different total uncertainty bounds. With 4 times the finite element number, average times per outer iteration rose by a factor of 4.8, while the average number of iterations for the lower-level problem solution actually decreased by a factor of 0.82.

6 Conclusion

We proposed a model and an algorithm for the robust optimal compliance design with unpredictable material variations and applied it successfully to a number of numerical examples. For additive manufacturing examples the effect of the uncertain material parameters turned out to be minor even in the worst case. This is an interesting result by itself. On the other hand, the robust topology optimization accounting for material degradation showed very promising results. A rigorous worst-case degeneration formulation could be solved with a still manageable increase in computational cost compared to the nominal solution. Moreover, the robustly optimized topologies showed a much better worst-case performance than the nominal ones.

One partly open point is however that the exact nature of robustness of the optimal result and its performance is not entirely clear. While the worst-case performance for the robust topologies improved strongly, it consistently remained far away from that of the nominal solution. This might be partly due to the nature of the studied examples, however another reason might be the fact that the convexification we applied for the lower level problem leads to conservative solutions. To get a feeling about the latter, a RAMP continuation strategy was suggested, which lead to a binary worst-case material degeneration. By this, the gap could be tightened. On the other hand there is no guarantee, that the worst case was identified in this case. Moreover, the continuation strategy came at a substantially increased computational demand.

Future research directions include studying more elaborate schemes for the RAMP continuation and possibly different solvers for the lower-level problem in the non-convex RAMP case. For additive manufacturing, different cost functionals might be considered which are more sensitive to relatively small changes in the material parameters. One difficulty in this case would be however that even with the inverse interpolation model the inner maximization problem might fail to be convex. As a consequence alternative methods from global optimization would be required to identify the worst-case material distribution.

Acknowledgements.

The authors thank the German Research Foundation (DFG) for funding this research work within Collaborative Research Centre 814, subproject C2.