A First Step Towards the Use of Proper General Decomposition Method for Structural Optimization

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.


Introduction
A classical problem of structural optimization consists in determining the optimal set of n parameters α opt = (α opt i ) i=1,n such that where C is a cost function, and α = (α i ) i=1,n are parameters of different types: geometrical and material; the cost function can be written as In this equation, F is a functional of the space (x)-time (t) field u which can be a scalar, a vectorial or a tensorial field, e.g. temperature, displacement, stress. u is usually the solution of a structural problem, i.e. a set of partial differential equations that can be non-linear; the subscript ·| α signifies that the parameters (α 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 u as an implicit function of the design variables α. 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 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ũ (x, t, α), (ii) defining a cost function and minimizing it to obtain α 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.

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

General Statements
Consider a differential problem with both Dirichlet and Neumann boundary conditions. u is the unknown field and depends on a large number of variables ,D and then the dimensionality of the problem is D j =1 d j . Note that the variables (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 The product F i is called the i th mode of the solution, and the number of modes N is chosen large enough to ensure convergence at a given precision ε, i.e.
Practically the solution (5) is built incrementally, i.e. mode by mode. We consider the weak form of the problem (4) where u * satisfies the Dirichlet boundary conditions. We assume that the first (i − 1) modes have been determined and we notẽ the previous solution. Then, determining the i th mode consists in calculating functions 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 (x j ) are calculated iteratively by a fixedpoint algorithm that consists in alternatively considering, until convergence, all the functions known except one, e.g. F J (x J ). Then, the trial function u * is chosen as So, the function F J (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 x J is the time coordinate) or boundary (e.g. when x J is a space coordinate) conditions. The weak form of this problem is which must be satisfied for all functions F * J (x J ). Noting the set of the known variables, i.e. and (11) can be written under the following form and furthermore, the functions (F j ) j =1,D j =J being known, The new operators involved in this equation are where L(ũ i−1 ) is the residual of the solution with (i − 1) modes, and Remark that (15) can be seen as a particular 1 strong form of the problem:

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 (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 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.

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 (λ 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 x and the space domains of the layers are ( x i ) i=1,n such as: As proposed above, thicknesses of the layers are considered as variables of the problem, then the temperature field can be written as u(x, (λ i ) i=1,n ). The strong problem to solve is where D(x) is the thermal conductivity and f (x) is the heat source, both of them can explicitly depends on x. Considering 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: Each variable λ i is defined in a finite interval [λ min i , λ max i ]. In fact, in the classical approach of optimization, these restrictions would have been expressed in terms of inequality constraints, e.g. λ i < λ max i . If the variables (λ 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: where u * (x) is a trial function which satisfies the Dirichlet boundary conditions. Since the physical domain depends on (λ i ) i=1,n , it is convenient to introduce a fixed parent domain s = [0, n] and a diffeomorphism φ : s → x which depends also on (λ i ) i=1,n . More precisely, we take the following piecewise linear change of variable: The parent subdomains s i , are defined as: With this specific change of variable, the gradient of φ(s) is constant over each subdomain: We can now define the generalized weak problem which solution is u(s, (λ i ) i=1,n ): where we have omitted the s-and (λ i ) i=1,n -dependence of both u and u * . This last equation can conveniently be rewritten as 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.

Derivation of the Separated Problem
From (28) it is straightforward using the PGD method to seek a solution u(s, (λ i ) i=1,n ) of the form: Note that the specific choice of mapping between s and x renders the definition of strong problems per dimension a bit difficult because the derivative of the mapping is not continuous everywhere. Nevertheless, it does not change the rest of the method, because the iterative determination of S k (s) and k i (λ 1 ) 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), 1 1 (λ 1 ) · · · 1 n (λ 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 ( 1 i (λ i )) i=1,n = 1.

Initialize the functions of the new mode S k (s), k
1 (λ 1 ), . . . , k n (λ 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 (λ i ) i=1,n and therefore computing the value of ( k i (λ i )) i=1,n reduces to the resolution of an algebraic equation.

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 λ 1 and λ 3 are made of the same material and their thermal conductivity is set to 1 W/(m K). The middle layer of thickness λ 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.

Discretization and Convergence
As solving the four-dimensional (s, λ 1 , λ 2 , λ 3 ) problem only involves four one-dimensional problems along each dimension, simple one-dimensional discretization methods are considered. s-space as well as (λ 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 (λ 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 (λ 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 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 storẽ u is negligible compared to what would have been required to store it on an equivalent four-dimensional grid. As described in the Introduction, the use of the PGD method leads to the determination of an approximated solution of the problemũ for all combinations of the parameters (s, λ 1 , λ 2 , λ 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ũ 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ũ(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. λ 1 = λ 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. λ 1 + λ 2 + λ 3 = 3, are shown. geometrical parameters λ 1 and λ 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).

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, λ 1 , λ 2 , λ 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 here that the previously computed solution of the augmented structural problem can be used to compute even very complex cost functions: under the constraint: 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 γ 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 λ 1 = λ 2 = 1; only λ 3 is considered as a parameter. The optimization problem consists in determining the position x in the middle layer in which the temperatureũ(x) is the most representative of the mean temperature in the whole domain. So, the cost function can be defined as the quadratic error between the whole temperature field and the temperature in x 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.

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 (λ i ) i=1,3 in their respective range. A brute force approach, i.e. solving the problem for each combination of (λ 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.