# A first step towards the use of Proper General Decomposition method for structural optimization

Adrien Leygue, Erwan Verron

► To cite this version:

Adrien Leygue, Erwan Verron. A first step towards the use of Proper General Decomposition method for structural optimization. Archives of Computational Methods in Engineering, 2010, 17 (4), pp.465-472. 10.1007/s11831-010-9052-3. hal-01004831

HAL Id: hal-01004831

https://hal.science/hal-01004831

Submitted on 11 Jun 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A First Step Towards the Use of Proper General Decomposition Method for Structural Optimization

A. Leygue · E. Verron

A. Leygue · E. Verron
Institut de Recherche en Génie Civil et Mécanique—GeM, UMR
CNRS 6183, Ecole Centrale Nantes, BP 92101, 44321 Nantes
cedex 3, France
e-mail: erwan.verron@ec-nantes.fr

A. Leygue
e-mail: adrien.leygue@ec-nantes.fr

## 1 Introduction

A classical problem of structural optimization consists in determining the optimal set of n parameters $\alpha^{\text{opt}} = (\alpha_i^{\text{opt}})_{i=1,n}$ such that

$$ \alpha^{\text{opt}} = \min_{\alpha} C(\alpha) \qquad (1) $$

where C is a cost function, and $\alpha = (\alpha_i)_{i=1,n}$ are parameters of different types: geometrical and material; the cost function can be written as

$$ C(\alpha) = \mathcal{F}(u(\mathbf{x}, t)|_{\alpha}). \qquad (2) $$

In this equation, $\mathcal{F}$ is a functional of the space $(\mathbf{x})$-time $(t)$ field $\mathbf{u}$ which can be a scalar, a vectorial or a tensorial field, e.g. temperature, displacement, stress. $\mathbf{u}$ is usually the solution of a structural problem, i.e. a set of partial differential equations that can be non-linear; the subscript $|_{\alpha}$ signifies that the parameters $(\alpha_i)_{i=1,n}$ are fixed. For complete reviews on optimization problems, the reader can refer to the recent papers: [3] and the references herein for general optimization methods, and [8] for the special problem of structural topology optimization. This later type of optimization problem is not considered here; we only focus on optimization problems in which simple geometrical quantities are optimized.

The most common optimization strategy consists in considering only the design variables as optimization variables; it is referred to as the conventional formulation in [3], or Nested ANalysis and Design (NAND). The central idea of this approach is to treat the field $\mathbf{u}$ as an implicit function of the design variables $\alpha$. Consequently the evaluations of the cost function and its gradient require the repeated solution of the structural problem, classically by the finite element method. In order to circumvent the implicit nature

**Abstract** In structural optimization, the implicit nature of the cost function with respect to the optimization parameters, i.e. through the solution of the structural problem calculated with fixed values of these parameters, leads to prohibitive computations whatever the adopted formulation.

Consequently, it yields limitations in both the number of parameters and the size of the structural problem. Moreover, some know-how is required to define a relevant structural problem and a well-behaved cost function.

Here, we profit from the ability of the Proper Generalized Decomposition (PGD) method to handle large-dimensionality problems to transform the optimization parameters into variables of an augmented-structural problem which is solved prior to optimization. As a consequence, the cost function becomes explicit with respect to the parameters.

As the augmented-structural problem is solved a priori, it becomes independent from the a posteriori optimization. Obviously, such approach promises numerous advantages, e.g. the solution of the structural problem can be easily analyzed to provide a guide to define the cost function and advanced optimization schemes become numerically tractable because of the easy evaluation of the cost function and its gradients.
of the cost function with respect to the optimization parameters, Schmit et al. proposed to treat simultaneously both the design and the field variables as optimization variables [9]. The structural problem then reduces to equality constraints; therefore no explicit structural analysis is needed. This method referred to as Simultaneous Analysis and Design (SAND) in the context of structural optimization was applied by Fuchs [6, 7] to trusses. A summary of the respective advantages and disadvantages of the above-mentioned methods is presented in Table 2 of [3]. In both methods, difficulties arise from large-scale problems. Obviously, for the conventional method, numerous solutions of the large scale structural problem rapidly become numerically prohibitive. For the SAND method the huge number of equality constraints necessitates specific optimization algorithms. Moreover, large scale problems produce a number of optimization variables that can exceed the capacity of optimization codes and computers. Consequently, one must reduce a priori his ambition by carefully selecting the most relevant optimization parameters and cost function. Therefore, the definition of optimization problems is usually constrained by numerical and algorithmic possibilities.

With the above classical methods, the primary problem of an optimization study consists in the evaluation and minimization of the cost function C and changing the cost function leads to the definition of a new problem. The structural problem (which solution is u) appears only as secondary: in the evaluation of the cost function or as additional constraints. Our objective in this work is to change this paradigm by considering the structural problem as the primary problem and the minimization of the cost function as the secondary one. Our strategy consists in changing the relationship between the design parameters α and the field u from implicit to explicit. In this way, the optimization parameters α are no longer considered as parameters but as additional variables of an augmented structural problem. The resulting optimization problem consists in:

(i) solving the augmented structural problem which solution (exact or approximated) is $\tilde{\mathbf{u}}(\mathbf{x}, t, \alpha)$,

(ii) defining a cost function

$$C(\alpha) = \mathcal{F}(\tilde{\mathbf{u}}(\mathbf{x}, t, \alpha)), \quad (3)$$

and minimizing it to obtain $\alpha^{\text{opt}}$.

This proposal may appear unrealistic. Nevertheless, recent advances in numerical methods can help us to handle such approach: the so-called Proper Orthogonal Decomposition (PGD) method has been successfully used to solve large dimensionality problems [1, 2, 4]. In the present paper, we propose to extend this method to optimization problems by considering geometrical parameters as new dimensions of the structural problem which is solved by the PGD method. We study a simple multilayer thermal problem which additional variables are the thicknesses of the layers.

# 2 Methods

## 2.1 Basics of the PGD Method

Here, we only recall the basics of the method. For more complete information, the reader can refer to [5].

### 2.1.1 General Statements

Consider a differential problem

$$\mathcal{L}(\mathbf{u}) = 0, \quad (4)$$

with both Dirichlet and Neumann boundary conditions. u is the unknown field and depends on a large number of variables $\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_D$ such as $(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_D) \in \Omega_1 \times \Omega_2 \times \dots \times \Omega_D$. $(\Omega_{d_i})_{i=1,D}$ are subsets of $(\mathbb{R}^{d_i})_{i=1,D}$ and then the dimensionality of the problem is $\sum_{j=1}^D d_j$. Note that the variables $(\mathbf{x}_i)_{i=1,D}$ can be space variables, time, temperature, but also material parameters, loading conditions, ...

The PGD method consists in searching an approximated solution of the problem (4) under a separated form

$$\tilde{\mathbf{u}}(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_D) = \sum_{i=1}^{N} F_1^i(\mathbf{x}_1) F_2^i(\mathbf{x}_2) \cdots F_D^i(\mathbf{x}_D). \quad (5)$$

The product $F_1^i(\mathbf{x}_1) F_2^i(\mathbf{x}_2) \cdots F_D^i(\mathbf{x}_D)$ is called the $i^{\text{th}}$ mode of the solution, and the number of modes N is chosen large enough to ensure convergence at a given precision $\epsilon$, i.e.

$$\int_{\Omega} |\mathcal{L}(\mathbf{u})|^2 d\mathbf{x}_1 d\mathbf{x}_2 \cdots d\mathbf{x}_D < \epsilon. \quad (6)$$

Practically the solution (5) is built incrementally, i.e. mode by mode. We consider the weak form of the problem (4)

$$\int_{\Omega_1 \times \Omega_2 \times \cdots \times \Omega_D} \mathbf{u}^* \cdot \mathcal{L}(\mathbf{u}) d\mathbf{x}_1 d\mathbf{x}_2 \cdots d\mathbf{x}_D = 0 \quad \forall \mathbf{u}^* \quad (7)$$

where $\mathbf{u}^*$ satisfies the Dirichlet boundary conditions. We assume that the first $(i-1)$ modes have been determined and we note

$$\tilde{\mathbf{u}}_{i-1}(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_D) = \sum_{k=1}^{i-1} F_1^k(\mathbf{x}_1) F_2^k(\mathbf{x}_2) \cdots F_D^k(\mathbf{x}_D) \quad (8)$$

the previous solution. Then, determining the $i^{\text{th}}$ mode consists in calculating functions $F_1(\mathbf{x}_1), F_2(\mathbf{x}_2), \dots, F_D(\mathbf{x}_D)$ such as

$$|\mathcal{L}(\tilde{\mathbf{u}}_{i-1} + F_1 F_2 \cdots F_D)| < |\mathcal{L}(\tilde{\mathbf{u}}_{i-1})|. \quad (9)$$

Moreover, the first mode is classically chosen such as it satisfies the Dirichlet and Neumann boundary conditions, and then next modes must only satisfy homogeneous boundary conditions.
The functions $F_J(\mathbf{x}_J)$ are calculated iteratively by a fixed-point algorithm that consists in alternatively considering, until convergence, all the functions known except one, e.g. $F_J(\mathbf{x}_J)$. Then, the trial function $\mathbf{u}^*$ is chosen as

$$ \mathbf{u}^*(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_D) = F_1(\mathbf{x}_1) F_2(\mathbf{x}_2) \cdots F_J^*(\mathbf{x}_J) \cdots F_D(\mathbf{x}_D). \quad (10) $$

So, the function $F_J(\mathbf{x}_J)$ becomes simply the solution of a differential problem in a $d_J$-dimensional space. In particular, we have an ordinary differential equation with initial (e.g. when $\mathbf{x}_J$ is the time coordinate) or boundary (e.g. when $\mathbf{x}_J$ is a space coordinate) conditions. The weak form of this problem is

$$ \int_{\Omega_1 \times \Omega_2 \times \cdots \times \Omega_D} F_1 F_2 \cdots F_J^* F_D \cdot \mathcal{L}(\tilde{\mathbf{u}}_{i-1} + F_1 F_2 \cdots F_D) d\mathbf{x}_1 d\mathbf{x}_2 \cdots d\mathbf{x}_D = 0, \quad (11) $$

which must be satisfied for all functions $F_J^*(\mathbf{x}_J)$. Noting $\Omega'$ the set of the known variables, i.e.

$$ \Omega' = \Omega_1 \times \Omega_2 \times \cdots \times \Omega_{J-1} \times \Omega_{J+1} \times \cdots \times \Omega_D, \quad (12) $$

and

$$ d\mathbf{x}' = d\mathbf{x}_1 d\mathbf{x}_2 \cdots d\mathbf{x}_{J-1} d\mathbf{x}_{J+1} \cdots d\mathbf{x}_D, \quad (13) $$

(11) can be written under the following form

$$ \int_{\Omega_J} F_J^* \cdot \left[ \int_{\Omega'} F_1 F_2 \cdots F_{J-1} F_{J+1} \cdots F_D \cdot \mathcal{L}(\tilde{\mathbf{u}}_{i-1} + F_1 F_2 \cdots F_D) d\mathbf{x}' \right] d\mathbf{x}_J, \quad (14) $$

and furthermore, the functions $(F_J)_{j=1,D \ j \neq J}$ being known,

$$ \int_{\Omega_J} F_J^* \cdot [\gamma_J + \mathcal{L}'(F_J)] d\mathbf{x}_J. \quad (15) $$

The new operators involved in this equation are

$$ \gamma_J = \int_{\Omega'} F_1 F_2 \cdots F_{J-1} F_{J+1} \cdots F_D \cdot \mathcal{L}(\tilde{\mathbf{u}}_{i-1}) d\mathbf{x}' \quad (16) $$

where $\mathcal{L}(\tilde{\mathbf{u}}_{i-1})$ is the residual of the solution with $(i-1)$ modes, and

$$ \mathcal{L}'(F_J) = \int_{\Omega'} F_1 F_2 \cdots F_{J-1} F_{J+1} \cdots F_D \cdot \mathcal{L}(F_1 F_2 \cdots F_D) d\mathbf{x}'. \quad (17) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1 The multilayer body: notations</figcaption></figure>

Remark that (15) can be seen as a particular<fnref n="1" /> strong form of the problem:

$$ \gamma_J + \mathcal{L}'(F_J) = 0. \quad (18) $$

### 2.1.2 Numerical Treatment

As mentioned above, the determination of a new mode consists in solving a series of ordinary differential equations (unknown $F_J$) in weak form (see (15)) with homogeneous initial or multi-point boundary conditions. In the fixed-point algorithm, each problem, i.e. the determination of each function $F_J(\mathbf{x}_J)$, is independent to each other, and can be solved by the more appropriate numerical method: finite difference method, finite element method, spectral method, ... As a classical example, if $\mathbf{x}_J$ is the time variable, the corresponding problem is classically solved by finite differences; if it is a space variable, finite elements are usually used.

## 2.2 A Simple Structural Problem

### 2.2.1 Definition

The problem consists in optimizing a multilayer body subjected to given temperatures on its external faces for different cost functions. Only the steady-state heat transfer problem is considered and the unknown temperature field is denoted $u$. As shown in Fig. 1, the body is composed of $n$ layers and the boundary temperatures are $u_0$ and $u_n$. The thicknesses of the layers are denoted $(\lambda_i)_{i=1,n}$. For sake of simplicity, we consider a one-dimensional problem in space, the space variable is $x$, the total space domain is denoted $\Omega^x$ and the space domains of the layers are $(\Omega_i^x)_{i=1,n}$ such as:

$$ \Omega^x = \bigcup_{i=1}^{n} \Omega_i^x = \left[ 0, \sum_{i=1}^{n} \lambda_i \right] \quad (19) $$

<footnote n="1">¹Particular because the strong form depends on the choice of the trial function (10).</footnote>
with

$$
\Omega_1^x = [0, \lambda_1] \quad \text{and} \quad \Omega_i^x = \left[ \sum_{j=1}^{i-1} \lambda_j, \sum_{j=1}^{i} \lambda_j \right]
$$
$$
\text{for } i=2, n. \tag{20}
$$

As proposed above, thicknesses of the layers are considered
as variables of the problem, then the temperature field can
be written as $u(x, (\lambda_i)_{i=1,n})$. The strong problem to solve is

$$
\frac{\partial}{\partial x} \left( D(x) \frac{\partial}{\partial x} u(x, (\lambda_i)_{i=1,n}) \right) + f(x) = 0, \quad (21)
$$

where $D(x)$ is the thermal conductivity and $f(x)$ is the heat
source, both of them can explicitly depends on $x$. Consider-
ing layers of different homogeneous materials, $D(x)$ is set
constant per layer ($D_i$)$_{i=1,n}$. So, referring to the notations
introduced in the previous section, the variables are:

$$
x_1 = x; x_2 = \lambda_1; \dots; x_i = \lambda_{i+1}; \dots; x_{n+1} = \lambda_n. \quad (22)
$$

Each variable $\lambda_i$ is defined in a finite interval $[\lambda_i^{\min}, \lambda_i^{\max}]$.
In fact, in the classical approach of optimization, these restrictions would have been expressed in terms of inequality constraints, e.g. $\lambda_i < \lambda_i^{\max}$.

If the variables $(\lambda_i)_{i=1,n}$ governing the physical dimensions of the problem were known, the solution would depend on $x$ only and the weak problem (after integration by parts) would read:

$$
\int_{\Omega^x} \frac{\partial}{\partial x} u^*(x) \left( D(x) \frac{\partial}{\partial x} u(x) \right) - u^*(x) f(x) dx = 0, \quad (23)
$$

where $u^*(x)$ is a trial function which satisfies the Dirich-
let boundary conditions. Since the physical domain depends
on $(\lambda_i)_{i=1,n}$, it is convenient to introduce a fixed parent do-
main $\Omega^s = [0, n]$ and a diffeomorphism $\phi : \Omega^s \to \Omega^x$ which
depends also on $(\lambda_i)_{i=1,n}$. More precisely, we take the fol-
lowing piecewise linear change of variable:

$$
\phi(s) = \begin{cases} s\lambda_1 & \text{if } 0 \le s \le 1, \\ (s-j+1)\lambda_j + \sum_{k=1}^{j-1} \lambda_k & \text{if } (j-1) \le s \le j. \end{cases} \tag{24}
$$

The parent subdomains $\Omega_i^s$, are defined as:

$$
\Omega_i^s = \phi^{-1}(\Omega_i^x) = [(i-1), i]. \tag{25}
$$

With this specific change of variable, the gradient of $\phi(s)$ is
constant over each subdomain:

$$
\frac{d\phi}{ds}(s) = \lambda_i \quad \text{for } s \in \Omega_i^s. \tag{26}
$$

We can now define the generalized weak problem which so-
lution is $u(s, (\lambda_i)_{i=1,n})$:

$$
\int_{\lambda_1^{\min}}^{\lambda_1^{\max}} \dots \int_{\lambda_n^{\min}}^{\lambda_n^{\max}} \int_0^n \left[ \left( \left( \frac{d\phi}{ds}(s) \right)^{-2} \frac{\partial u^*}{\partial s} D(s) \frac{\partial u}{\partial s} - u^* f \right) \frac{d\phi}{ds}(s) \right] ds \, d\lambda_1 \cdots d\lambda_n, \quad (27)
$$

where we have omitted the s- and $(\lambda_i)_{i=1,n}$-dependence of
both $u$ and $u^*$. This last equation can conveniently be rewritten as

$$
\int_{\lambda_1^{\min}}^{\lambda_1^{\max}} \dots \int_{\lambda_n^{\min}}^{\lambda_n^{\max}} \left[ \sum_{i=0}^{n} \int_{\Omega_i^s} \left( \frac{\partial u^*}{\partial s} \frac{D_i}{\lambda_i} \frac{\partial u}{\partial s} - \lambda_i u^* f \right) ds \right] d\lambda_1 \cdots d\lambda_n. \quad (28)
$$

Furthermore, in this specific example, we set $f(x) = 0$ in the second term of the integrand, although considering more complex sources $f(x)$ would not yield to major difficulties.

2.2.2 Derivation of the Separated Problem

From (28) it is straightforward using the PGD method to
seek a solution $u(s, (\lambda_i)_{i=1,n})$ of the form:

$$
u(s, (\lambda_i)_{i=1,n}) = \sum_{k=1}^{N} S^k(s) \Lambda_1^k(\lambda_1) \cdots \Lambda_n^k(\lambda_n). \quad (29)
$$

Note that the specific choice of mapping between $\Omega^s$ and $\Omega^x$
renders the definition of strong problems *per dimension* a bit
difficult because the derivative of the mapping is not contin-
uous everywhere. Nevertheless, it does not change the rest
of the method, because the iterative determination of $S^k(s)$
and $\Lambda_i^k(\lambda_i)_{i=1,n}$ is carried out by solving the weak form of
these problems. The solution procedure goes as follows:

1. Arbitrarily define the functions $S^1(s)$, $\Lambda_1^1(\lambda_1) \cdots \Lambda_n^1(\lambda_n)$ of the first mode, such that they satisfy the non-homogeneous boundary conditions of our problem. In our case we select $S^1(s) = u_0(1 - s/n) + u_n s/n$ and $(\Lambda_i^1(\lambda_i))_{i=1,n} = 1$.

2. Initialize the functions of the new mode $S^k(s)$, $\Lambda_1^k(\lambda_1)$, \dots, $\Lambda_n^k(\lambda_n)$ to some arbitrary values. The only constraint is on the boundary values of these functions. Indeed, in order to preserve the values at the boundary set by the first mode and therefore satisfy the boundary conditions of our problem we must ensure that $S^k(0) = S^k(n) = 0$.

3. Compute the value of the new mode using a fixed-point algorithm as described in the previous section. It is worth mentioning that our problem does not involve any partial derivatives with respect to $(\lambda_i)_{i=1,n}$ and therefore computing the value of $(\Lambda_i^k(\lambda_i))_{i=1,n}$ reduces to the resolution of an algebraic equation.
<table><caption>**Table 1** Discretization of the four variables</caption><thead><tr><th>dimension</th><th>range</th><th>number of grid points</th></tr></thead><tbody><tr><td>s</td><td>[0, 3]</td><td>301</td></tr><tr><td>λ<sub>1</sub></td><td>[0.1, 1.5]</td><td>25</td></tr><tr><td>λ<sub>2</sub></td><td>[0.1, 2.8]</td><td>25</td></tr><tr><td>λ<sub>3</sub></td><td>[0.1, 1.5]</td><td>25</td></tr></tbody></table>

4. Check for convergence and add a new mode if the norm of the residual is too large.

## 3 Results

In this section we present the numerical results obtained with the PGD method. All the results are computed for three physical subdomains of variable size ($n = 3$). The two extreme layers of thickness $\lambda_1$ and $\lambda_3$ are made of the same material and their thermal conductivity is set to 1 W/(m K). The middle layer of thickness $\lambda_2$ is supposed highly more insulating, its thermal conductivity is set to 0.05 W/(m K). Finally, the boundary conditions are of the Dirichlet type, and set to $u_0 = 0$ K and $u_3 = 1$ K.

### 3.1 Discretization and Convergence

As solving the four-dimensional ($s, \lambda_1, \lambda_2, \lambda_3$) problem only involves four one-dimensional problems along each dimension, simple one-dimensional discretization methods are considered. $s$-space as well as $(\lambda_i)_{i=1,3}$-spaces are discretized by piecewise linear finite elements and the corresponding weak problems are solved by the Galerkin method. Both range and number of grid points for each of the dimensions are given in Table 1. As no partial derivatives with respect to the thicknesses $(\lambda_i)_{i=1,3}$ are involved in the problem (see step 3 in Sect. 2.2.2), for each fixed-point iteration three algebraic systems (in terms of $(\lambda_i)_{i=1,3}$) and one differential problem (in terms of $s$) must be solved. Nevertheless, choosing a finite element discretization for the thicknesses leads to an easy interpolation for out-of-grid values.

34 modes are sufficient to satisfy the convergence criterion (6), as applied to the present problem

$$ \int_{0.1}^{1.5} \int_{0.1}^{2.8} \int_{0.1}^{1.5} \int_{0}^{3} |\mathcal{L}(\tilde{u})|^2 ds d\lambda_1 d\lambda_2 d\lambda_3 < 10^{-4}. \quad (30) $$

As mentioned above, the first mode is chosen to satisfy the non-homogeneous boundary conditions (see step 1 in Sec. 2.2.2); and for each of the 33 correction modes, between 9 and 81 fixed-point iterations are necessary to achieve convergence of the modes to machine precision.

Finally as we compute a variable-separated representation of the solution, the amount of memory required to store $\tilde{u}$ is negligible compared to what would have been required to store it on an equivalent four-dimensional grid.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Temperature field in three different multi-layer domains</figcaption></figure>

### 3.2 Solution of the Primary Problem: Obtaining $\tilde{u}(x, \lambda_1, \lambda_2, \lambda_3)$

As described in the Introduction, the use of the PGD method leads to the determination of an approximated solution of the problem $\tilde{u}$ for all combinations of the parameters $(s, \lambda_1, \lambda_2, \lambda_3)$ in the ranges given in Table 1. In order to illustrate the power of the method and the potential use of the computed solution beyond optimization problems, we present some particular results. For readability, they are restricted to cases in which only three variables are independent, although results with the four independent variables are available.

* As a first example, we present the evolution of the temperature $\tilde{u}$ as a function of the space $x$ for different thicknesses of the layers in Fig. 2. This figure shows that the introduction of the parent domain parametrized by the curvilinear coordinate $s$ is a relevant tool to introduce the change of geometry in the governing equations.

A similar conclusion can be drawn by considering Fig. 3. In this figure, the temperature field in the whole domain $\tilde{u}(x)$ is plotted as a function of the thickness of the middle layer. Here, the thicknesses of the extreme layers were set to 1, i.e. $\lambda_1 = \lambda_3 = 1$.

* The next two figures illustrate that the four variables are considered of the same nature in the method. To simplify the discussion, only results in which the whole thickness of the domain is set to 3, i.e. $\lambda_1 + \lambda_2 + \lambda_3 = 3$, are shown. Figures 4 and 5 present respectively the value of the temperature gradient in the origin, i.e. $\partial\tilde{u}/\partial x(x = 0)$, and the value of the temperature field in the middle of the whole domain, i.e. $\tilde{u}(x = 1.5)$, as functions of the thicknesses of the two extreme layers, i.e. $\lambda_1$ and $\lambda_3$. The power of the method is demonstrated by these two examples: the
<figure><img src="image_4.png" /><figcaption>Fig. 3 Temperature field in the domain as a function of the thickness of the middle layer $\lambda_2$ (with $\lambda_1 = \lambda_3 = 1$)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Temperature gradient at $x = 0$ (with $\lambda_1 + \lambda_2 + \lambda_3 = 3$)</figcaption></figure>

geometrical parameters $\lambda_1$ and $\lambda_3$ are considered as variables of the problem and then the solution of the thermal problem is available for all their possible values (in their respective range given in Table 1).

### 3.3 Solution of the Secondary Problem: Calculating the Cost Function

Once the primary problem solved, we can choose a cost function and minimize it with respect to any set of variables in $(x, \lambda_1, \lambda_2, \lambda_3)$ to determine the optimal set of parameters. In order to illustrate the versatility of our approach, we adopt successively two different cost functions.

- The first optimization problem consists in finding the set of thicknesses that minimize a complex cost function. We will not actually solve this problem but we demonstrate

<figure><img src="image_6.png" /><figcaption>Fig. 5 Temperature in the middle of the domain (with $\lambda_1 + \lambda_2 + \lambda_3 = 3$)</figcaption></figure>

here that the previously computed solution of the augmented structural problem can be used to compute even very complex cost functions:

$$C(\lambda_1, \lambda_2, \lambda_3) = \left( \sum_{i=1}^{3} c_i \lambda_i \right) \left( \max_i \gamma_i(\lambda_1, \lambda_2, \lambda_3) \right), \quad (31)$$

with

$$\gamma_i(\lambda_1, \lambda_2, \lambda_3) = d_i \int_{\Omega_i^x} \frac{\partial \tilde{u}}{\partial x} \tilde{u} dx, \quad (32)$$

under the constraint:

$$\lambda_1 + \lambda_2 + \lambda_3 = 3. \quad (33)$$

In the previous expressions, $c_i$ and $d_i$ are scalar parameters. In order to interpret this cost function, one should see $c_i$ as the cost of each layer per unit length and $\gamma_i$ as a reciprocal lifetime of each layer. The cost function can therefore be interpreted as the cost of the structure divided by the minimum lifetime of its parts. Figure 6 presents the value of the cost function as a function of the thicknesses of the extreme layers (the thickness of the middle layer is deduced from the previous constraint (33)). The following numerical values were used: $c_i = \{1, 50, 1\}$, $d_i = \{1, 0.05, 1\}$. As exhibited in the figure, the optimal set of parameters is easily determined by searching the minimum of the cost function in the two-dimensional space.

- The second cost function considered here is completely different than the previous one. We assume that $\lambda_1 = \lambda_2 = 1$; only $\lambda_3$ is considered as a parameter. The optimization problem consists in determining the position $x$ in the middle layer in which the temperature $\tilde{u}(x)$ is the most representative of the mean temperature in the whole domain.
<figure><img src="image_7.png" /><figcaption>Fig. 6 Cost function of (31)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Quadratic error between the whole temperature field and the temperature measured in each position of the middle layer</figcaption></figure>

So, the cost function can be defined as the quadratic error between the whole temperature field and the temperature in x

$$C(x, \lambda_3) = \frac{1}{\lambda_1 + \lambda_2 + \lambda_3} \int_0^{\lambda_1 + \lambda_2 + \lambda_3} (\tilde{u}(y) - \tilde{u}(x))^2 dy \quad (34)$$
for $x \in [\lambda_1, \lambda_1 + \lambda_2]$.

Such cost function can be used for example to determine the optimal position of a temperature sensor. Figure 7 shows the quadratic error as a function of the position in the middle layer and the thickness of the third layer. As in the previous example, the optimal set of parameters is easily obtained by directly minimizing the cost function.

## 4 Discussion

In the present paper, we have proposed a new way to handle structural optimization problems by using the Proper Generalized Decomposition method. Contrary to previous applications of the PGD method, we consider here geometrical, i.e. dimensions of the domain, as additional variables of an augmented structural problem. The ability of the PGD method to handle high dimensionality problems makes it possible to solve, with reasonable computing time, the structural problem for all possible values of these new variables in given ranges. Indeed, computing the solution of our problem using the PGD method requires to solve 937 differential problems. This number might be seen as huge but it is worth recalling that it leads to the solution of the structural problem for any value of the variables $(\lambda_i)_{i=1,3}$ in their respective range. A brute force approach, i.e. solving the problem for each combination of $(\lambda_i)_{i=1,3}$, would have required the solution of $25^3 = 15625$ differential problems. Additionally, the variable-separated representation for the solution maintains memory usage to very low levels. For the chosen discretization parameters, representing the full solution on a four dimensional grid would have required up to 35 megabytes of memory, while the 34 modes of the PGD solution fit in about 100 kilobytes of memory. With this approach, the structural problem is solved a priori and not in the optimization loop. In a sense it becomes the primary problem and then the secondary optimization step only consists in defining a cost function and finding its minimum values. Moreover, as the whole solution is already known, it is possible and numerically inexpensive to explore the problem before defining the cost function. Also, as evaluating the cost function is (almost) costless, potentially more efficient cost functions and optimization algorithms (including direct search) might become numerically tractable. This study demonstrates the relevance of the PGD method in the context of structural optimization. Nevertheless, it leaves some issues of fundamental importance unresolved. First, the ability of the method to handle mechanical problems, i.e. problems in which the displacement field is involved, has to be investigated. Second, the convergence of the PGD method on such problems must be thoroughly examined: here a fixed-point algorithm has been used to calculate the modes, but we also tested the minimization of the residual (see [5]) which failed to converge. Finally, the principal limitation of the approach is actually the necessity to define a parent domain that relates the geometry to the geometrical parameters; so it implies the existence of a diffeomorphism between the parent and the physical domains. Obviously, in other problems, the existence of such diffeomorphism can be questionable.
References

1. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some, classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. J Non-Newton Fluid Mech 139(3):153-176

2. Ammar A, Mokdad B, Chinesta F, Keunings R (2007) A new fam-
ily of solvers for some classes of multidimensional partial differ-
ential equations encountered in kinetic theory modelling of com-
plex fluids—part ii: Transient simulation using space-time sepa-
rated representations. J Non-Newton Fluid Mech 144(2-3):98-121

3. Arora JS, Wang Q (2005) Review of formulations for structural and mechanical system optimization. Struct Multidiscip Optim 30(4):251-272

4. Chinesta F, Ammar A, Joyot P (2008) The nanometric and micro-metric scales of the structure and mechanics of materials revisited:

an introduction to the challenges of fully deterministic numerical
descriptions. Int J Multiscale Comput Eng 6(3):191–213

5. Chinesta F, Ammar A, Cueto E (2009) Recent advances in the use of the proper generalized decomposition for solving multidimensional models. Arch Comput Methods Eng. doi:10.1007/s11831-010-9049-y

6. Fuchs M (1982) Explicit optimum design. Int J Solid Struct 18(1):13–22

7. Fuchs M (1983) Explicit optimum design technique for linear elastic trusses. Eng Optim 6(4):213–218

8. Rozvany GIN (2009) A critical review of established methods of structural topology optimization. Struct Multidiscip Optim 37(3):217–237

9. Schmit L, Fox R (1965) An integrated approach to structural synthesis and analysis. AIAA J 3(6):1104