EDP Sciences logo
Submit your paper
Search
Menu
Home All issues Volume 8 (2017) Int. J. Simul. Multidisci. Des. Optim., 8 (2017) A9 Full HTML
Open Access
Issue
Int. J. Simul. Multidisci. Des. Optim.
Volume 8, 2017
Article Number A9
Number of page(s) 11
DOI https://doi.org/10.1051/smdo/2017002
Published online 01 March 2017

© B. A. El Majd et al., Published by EDP Sciences, 2017

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

We focus on improving the computational efficiency on optimum shape algorithms for application in which the cost functional relies on the prior solution of a complex set of partial-differential equations (PDEs), such as those governing compressible aerodynamics (e.g. the Euler equations), or related coupled disciplines [1] such as structural mechanics (e.g. elasticity), or electromagnetics [2] (e.g. the Maxwell equations). Hence, each evaluation functional is computationally costly. In several papers [3–8], we introduced a conceptual parallel from grid to geometrical parametrization by constructing a full and adaptive multi-level optimum-shape algorithm (FAMOSA), analogous to the classical full multigrid method (FMG) known to have optimum linear complexity with respect to the number of degrees of freedom). We demonstrated the efficiency of this construction for problems of drag reduction for transonic flight, external noise reduction for supersonic flight, and aerodynamic and structural optimization of a business-jet wingshape. Figure 1 illustrates an unstructured mesh, composed of 31,124 nodes and 173,445 elements, generated around the wing, including a refined area in the vicinity of the shock. Figure 2 shows a comparison of the convergence for the three strategies; FMOSA is the most efficient, yielding a shape of better fitness using a smaller computational effort. Using such setting, the wing geometry of an aircraft in transonic flow has been optimized by the minimization of a functional of the form:(1) CD0 and CL0 are respectively the drag and lift coefficients corresponding to the initial shape (NACA 0012 section).

thumbnail Figure 1.

Initial wing shape and mesh in the symmetry plane.

thumbnail Figure 2.

Comparison of the convergence history for the three strategies.

Until now, free-optimization techniques, such as the classical simplex search method, the genetic algorithms (GAs) and the particular swarm optimizers (PSOs), has been employed to reach the optimum design. Despite their greater robustness, they require a great number of function evaluations and the cost increases with the number of design variable. Another alternative, according to the literature, consists to utilize a derivative-based optimization algorithms that are able to locate the optimum solutions in the expense of only few function evaluations (i.e. with reasonable CPU time). They are often supported by tools that compute the derivatives of the objective function; such as the adjoint techniques [9, 10], the Automatic Diffrentiation (AD) [11, 12], the surrogate models (Radial Basis Function(RBF), Artificial Neural Networks (ANN), etc.) [13].

We note that prior to us, we find many applications employing the derivative-based methods with Bézier parametrization. In particular, Marino et al. [14] has used triangular Bézier surface and the optimization is numerically solved by a constrained gradient-based minimization algorithm. In this work, Marino et al. find the gradient-based algorithm more suitable than the evolutionary algorithm. In reference [15], we find a combination of NURBS (Non-Uniform Rational B-Spline) and a gradient-based optimization algorithm to conduct shape optimization of stratospheric airships.

In this paper, we propose to solve the shape optimisation problem, within multilevel and adaptive strategies, by the BFGS (Broyden Fletcher Goldfarb and Shanno) Quasi Newton method. One of the key properties of BFGS method is its superlinear convergence, which is guarantee under reasonable assumptions [16]. In addition, it works well in practice, even when the initial iterate is far from the solution. Despite these advantages, the initialization of BFGS algorithm by identity may affects the performance of our hierarchical and adaption strategies; this later depend on the quality of transfer operators in order to rigorously construct embedded and nested search space. Thus, our proposed approach consists to well re-initialize the Hessian and the gradient at each change of the search parametric space.

This paper is organized as follow. First, we review the construction of multilevel algorithms, in the context of parametric shape optimization. Second, we present the proposed approach for multilevel and adaptive parametric shape optimization. We then present a numerical results on a model problem from caculus of variations. Finally, we conclude and mention some perspectives.

2 Nested Bézier parameterizations for multilevel shape representation

We begin with the simplest situation of a two-dimensional geometry for which we employ a Bézier shape representation:(2)in which the parameter t varies from 0 to 1, n is the degree of the parameterization,(3)is a Bernstein polynomial, , and(4)is the generic control point. The coordinates of these control points are split into two vectors(5)and we refer to the vector X as the support of the parameterization, and the vector Y as the design vector. Typically, we optimize the design vector for fixed support according to some physical criterion, such as drag reduction in aerodynamics. The somewhat unsymmetrical roles dispensed to the vectors X and Y are chosen to reduce (to n essentially) the dimension of the search space in the optimization phase, which is the most numerically costly and subject to numerical stiffness.

We also use the notation:(6)in which the vector . In all this article, only supports for which the sequence {xk} is monotone increasing are said to be admissible and considered throughout. Thus, the function x(t) is monotone-increasing and defines a one-to-one mapping of, say, [0,1] onto itself. Recall also the simple formula for the derivative:(7)in which Δ denotes the forward-difference operator (Δxk = xk+1 − xk) as well as the associated n × (n + 1) matrix.

In the prototypical case of an airfoil, we use such a parametric representation for both the upper and lower surfaces separately. The vertical slope at the leading edge is enforced by the conditions:(8)for both surfaces which assures a smooth match; at the trailing edge, we simply have:(9)for a continuous match.

Our geometrical construction employs the degree-elevation process, well-known in the Computer-Aided Design literature (see, e.g. [17]). This process permits to cast (Eq. (6)) into the following equivalent Bézier parameterization of degree n + 1:(10)in which the new control points are obtained from the former by convex combinations:(11)obtained by multiplying equation (6) by (1 − t) + t and grouping together the monomials in tk(1 − t)n+1−k, for each k.

From a theoretical viewpoint, our construction guarantees rigorously nested search spaces, and exact upward transfer operators (from low to high-degree parameterization). Note that, apart from the specified endpoints, the abscissas of the initial support X are not a subset of the abscissas of any support of a higher degree parameterization. Nevertheless, any Bézier curve given on the initial support can be expressed exactly on any other support of higher degree provided it results from the degree elevation process. The parameterizations are nested, or embedded in one another in this sense precisely.

Based on this concept of nested Bézier parameterizations, the following theorem introduces the upward and downward transfer operators.

Theorem 1. Let and X = {xk} two nested (by degree elevation process) parametrization supports of degree n′ = n + 1 and n respectively, ti(i = 0, …, n′) and qi(i = 0, …, n) two arbitrary partitions of the interval [0, 1]. Then, the upward relation (elevation) between X′ and X is given by the relation: (12) where the generalized Vandermonde matrix Bn′ and the rectangular matrix Bn,n′ are: (13)

(14) and the downward relation (reduction) between X′ and X is given by the relation: (15) where the generalized Vandermonde matrix and the rectangular matrix are: (16)

(17)

Proof. Let X = {xk} and X′ ={xk′} two nested Bézier parametrization supports of degree n′ and n respectively. It means that:(18)

Let {ti} (i = 0, …, n′) arbitrary partition of [0,1] and {ei} (i = 0, …, n′) a canonic base of . We first muliply the equation (18) by ek (k = 0, …, n′), then:(19)and by adding these equations (19) with respect to k, we find the following formula:(20)where(21)and(22)

Since the generalized Vandermonde matrix Bn′ is invertible, it follows that:(23)by the same way, we prove the second statement (15) by interchanging n and n′.

In what follows, we note by the upward operator from the lower level n to N (N > n) and by the downward operator. By construction, we have this following property:(24)

To test these two transfer operators, we use the following model problem:(25)in which x(t) is given, smooth and monotone increasing,(26)and(27) p and are, for specified ω(t) > 0 and α > 1, the pseudo-length of the arc, and the pseudo-area below the arc. This model problem has been studied in [18] to which we refer for a full description of the numerical test-case which corresponds to α ≈ 2.03, for the functional is known to be convex, and a certain ω(t) for which the minimizing shape is the half-thickness distribution of the RAE2822 airfoil.

Figure 3a (resp. 3b) illustrates the upward process using (resp. the downward process ). The curves corresponding to the both parametrizations are superimposed.

thumbnail Figure 3.

Illustration of and by using two embedded control polygons of degree n = 6 and N = 8. (a) Upward process, (b) downward process.

3 Hessian and gradient transfer

In this section, we present our approach within a multilevel algorithms and adaption procedure. The subject is to good initialize the optimization process when the parametric search space is changed (Sect. 3.1) and the control polygon is regularized (Sect. 3.2).

3.1 First case

As explained in the previous section, the transfer operators guarantee the construction of nested control polygon with higher degree (N > n). To achieve a good transfer of informations, we need to better update the Hessian and the gradient in order to conserve the efficiency of our multilevel strategies.

For fixed support X0, let (X0, Y0) the solution of the minimization of the problem (25) and γ0 the corresponding shape. By using the upward operator, we construct a new control polygon (X, Y) of degree N corresponding to the same shape γ0. It follows that:(28)and(29)where jN is the cost function of the problem (25) corresponding to the Bézier parametrization of degree N.

To express the Hessian and the gradient on the coarse level (of degree n), we differentiate two times the equation (29) with respect to Y0:(30)

(31)where .

By differentiating the equation (28) with respect to Y, we find:(32)it follows that(33)and(34)where . We note that (33) and (34) can also be obtained using the property (24).

Proposition 1. The Hessian matrix given by (34) is positive definite.

Proof. Let(35)

The matrix Bn is the BFGS approximation corresponding to the first optimization phase on the coarse level; thus, it is positive definite.

Let , then(36)where .

Since Wn ≠ 0 and Bn is positive definite, then xTBNx > 0. It follows that BN is positive definite.

As consequence, the following expression holds:(37)

Finally, we can update the Hessian by the matrices BN or , it depends on the BFGS method used, instead of identity when the parametrization degree is augmented or reduced. The second phase of the optimization in the fine level of degree N is initialized by an exact Hessian matrix.

3.2 Second case

Suppose that we are placed on the search space of degree n. the parametrisation adaption consists to adapt the support X by alternating two complementary phases:

  1. Optimization: optimize the design vector Y for fixed support X = X0 according to some criterion; let Y0 be the result of this phase.

  2. Regularization: given the parametrization (X0, Y0) of an approximate optimum shape γ0, the new support X1 is taken to be the better support for which the total variation (TV) in the components of the corresponding vector Y1 is minimal, such that the correspondent shape γ1 to (X1, Y1) approximates γ0 in the sense of least squares; substitute X1 to X0.

This adaption procedure has been studied extensively in [3] for two-dimentional parametrization. The extension to three-dimentional optimum design in aerpdynamics within the framework of the so-called free-form deformation has been proposed in reference [19]. It improves noticeably the convergence rate and reduces the numerical stifness. The main idea consists to redefine the geometrical representation by using regularization techniques. Thus, a new formulation of the Hessian and the gradient is required in order to re-initialize the adaption process by a good regularized solution with equivalent performance.

As the both control polygons (X0,Y0) and (X1,Y1) correspond to the same performance, it follows that:(38)i.e.(39)

This intrinsic formulation can be transformed into a parametric optimization by:(40)where and τ = τ(t, X) is related to the change of support X0 → X, and defined by this condition:(41)

To satisfy the unicity of the solution, we impose that the components of X0 are monotone increasing.

The derivative of the cost function is given by:(42)where(43)and(44)

The matrix A(X) is real-symmetric positive definite and depends linearly upon the vector X, thus(45)where  = A′(X) is a tensor of order 3, independent of X and ⊗ stands for the contracted product.

The minimization of with respect to Y is equivalent to . It follows that:(46)

In particular, (X0, Y0) and (X1,Y1) verify this linear system.

For X = X0, the first and second derivative of (40) with respect to Y0 are given by:(47)and(48)

By differentiating the linear system (46) with respect to Y0, we obtain:(49)where(50)

As A(X1) and do not depend on Y0, then:(51)

Proposition 2. If all the components of the support X are monotone increasing, then the matrix A(X) given by the formula (43) is positive definite.

Proof. Let (Y ≠ 0). We have(52)where Z = Bn(t)TY.

As the components of X are monotone increasing, then ΔX > 0 (Δ is the forward difference operator). So, its enough to prove that ZTZ ≠ 0. We proceed by contradiction and assume that ZTZ = 0, so Z = 0. As form a basis of , then Y = 0, which contradicts the first assumption.

Since X1 satisfy the condition of the Proposition 2, then the matrix A(X1) is invertible. Thus:(53)

As consequence,(54)and(55)

The formulas express the relation between the Hessian matrices corresponding to the initial control polygon (X0, Y0) and the regularized one (X1, Y1). It permits to re-initialize the optimization process after each parametrization adaption.

4 Numerical experiments

The proposed approach is tested using the model problem, described in Section 3.1, for the first case and the shape-reconstruction problem, described in Section 3.2, for the second one. The numerical optimum is determined by the Quasi-Newton algorithm [20] based on BFGS (Broyden-Fletcher-Goldfarb-Shanno) updated formula. This algorithm is a line search family method, one of the most powerful methods to solve unconstrained optimization problem. It exposes superlinear convergence; resource-intensivity is estimated as O(n2) per iteration for n-component argument vector. The search direction is given by sk = −Hk Gk, where Gk and Hk are respectively the gradient and the Hessian of the objective function. Hk a symmetric definite matrix, constructed at each iteration, which approximates the Hessian matrix and/or its inverse. In the literature, there are several formulas for the calculation of Hk+1 from Hk. We use here, the BFGS update formula [16],(56)where yk = Gk+1 − Gk and dk = xk+1 − xk.

As a starting point, H0 can be set to any symmetric positive definite matrix and very often the identity matrix. This initialization has a bad impact on the efficiency of our multilevel and adaption strategies.

In the following numerical experiments, the gradient is calculated analytically and we focus only on the impact of the Hessian transfer. For the model problem, the derivative of the cost function in the search space of degree n is given by:withand

We use here the formula of the derivative of y′(t) given by , in which Δ is the n × (n + 1) matrix associated with the forward-difference operator (ΔY = yk+1 − yk).

For the shape-reconstruction problem, the gradient corresponds to the linear system (46):

To Calculate the gradient on the search space of degree N > n (resp. N < n), we apply the upward operator (resp. the downward operator ) to X and Y and we substitute.

Figure 4 provides the convergence history of the basic algorithm for three optimization strategies: at full convergence, with Hessian transfer after 20 iterations and without Hessian transfer (the Hessian is initialized automatically by identity). And Figure 4 gives the corresponding optimal shapes. Clearly, all the optimization achieve the convergence and the main point is the fact that the strategy based on Hessian transfer has no effect on the convergence rate; so, this experiment validate the formulas of the Hessian and gradient.

thumbnail Figure 4.

Basic algorithm for N = 12; Comparison between three strategies: full convergence, with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.

In this experiment, the multilevel optimization is considered here for n = 4 and N = 12. Figures 5a and 5b give respectively the iterative performance and the correspondingoptimal shapes for the multilevel algorithm with and without Hessian transfer. Clearly, the algorithm based on Hessian transfer is faster in term of the convergence rate. This experiment confirms the impact on the performance of the multilevel algorithm when the Hessian matrix is initialized by identity.

thumbnail Figure 5.

Multilevel algorithm for n = 4 and N = 12; Comparison between two strategies: with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.

Lastly, we measure the effect of the hessian transfer on the accuracy of the optimization when using the adaption procedure. Figure 6a shows a comparison of the basic algorithm for n = 8, with the adaption procedure with and without Hessian transfer. The Figure 6b gives the corresponding optimal shape. As expected, the Hessian transfer guarantees a best convergence rate of the optimization.

thumbnail Figure 6.

Adaption algorithm for n = 8; Comparison between two strategies: with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.

5 Conclusion

The Hessian transfer is essential to preserve and enhance the quality of our multilevel strategies and adaption procedure. This approach permits to better initialize the BFGS Quasi-Newton method. The numerical experiments demonstrate clearly its efficiency in the case of two-dimentional model problem. This study, encourage developing a new multilevel geometrical structure in which, for example in the case of two-level ideal V-cycle [22, 23], the problem is solved to complete convergence on the coarse level by robust optimization algorithms (Genetic algorithm [GA], Particle swarm optimization [PSO], simplex method, etc.) and few iterations are enough (acts as preconditionners) in the fine level to alleviate the numerical stifness; BFGS Quasi-Newton can be used in this case. This variant is illustrated by the following schematic. Others strategies can be elaborated mimicking the multigrid-type strategies (saw-touth, V and W cycle, FMG, etc.) [21].

The proposed approach will be generalized for 3D optimum design in aerodynamics, structural mechanics, or coupled discplines within the framework of the so-called Free-Form deformation (FFD) [24]. Testing a great number of promising algorithmic strategies is also a part of our ongoing work.

References

  1. Abou El Majd B, Desideri J-A, Habbal A. 2010. Aerodynamic and structural optimization of a business-jet wingshape by a Nash game and an adapted split of variables. Mécanique & Industries, 11(3–4), 209–214. [CrossRef] [EDP Sciences] [Google Scholar]
  2. Ouchetto O, Abou El Majd B, Ouchetto H, Essakhi B, Zouhdi S. 2016. Homogenization of periodic structured materials with chiral properties. IEEE Transactions on Antennas and Propagation, 64(5), 1751–1758. [CrossRef] [Google Scholar]
  3. Abou El Majd B, Désidéri J-A, Janka A. 2004. Nested and self-adaptive Bézier parameterizations for shape optimization, International Conference on Control, Partial Differential Equations and Scientific Computing (dedicated to late Professor J.-L. Lions), Beijing, China, 13–16 September, 2004. [Google Scholar]
  4. Abou El Majd B, Désidéri J-A, Duvigneau R. 2008. Multilevel strategies for parametric shape optimization in aerodynamics. European Journal of Computational Mechanics, 17(1–2). [Google Scholar]
  5. Désidéri J-A, Duvigneau R, Abou El Majd B, Tang Z. 2007. Algorithms for efficient shape optimization in aerodynamics and coupled disciplines. 42nd AAAF Congress on Applied Aerodynamics, Sophia-Antipolis, France. [Google Scholar]
  6. Abou El Majd B, Désidéri JA, Do TT, Fourment L. 2005. Multilevel strategies and hybrid methods for shape optimization and application to aerodynamics and metal forming. Evolutionary and Deterministic Methods for Design, Optimization and Control with Applications to Industrial and Societal Problems Conference (EUROGEN 2005), September, 2005, pp. 12–14. [Google Scholar]
  7. Duvigneau R, Abou El Majd B, Désidéri J-A. 2008. Towards a self-adaptive parameterization for aerodynamic shape optimization. ESAIM: Proceedings. Vol. 22, pp. 169–174, EDP Sciences, 2008. [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  8. Duvigneau R, Abou El Majd B, Désidéri J-A. 2008. Aerodynamic design using hierarchical shape parameterizations for descent and Particle Swarm Optimization Methods. Numerical Analysis and Scientific Computing for Partial Differential Equations and Their Challenging Applications. CIMNE. [Google Scholar]
  9. Giannakoglou KC, Papadimitriou DI. 2008. Adjoint Methods for Shape Optimization. Archives of Computational Methods in Engineering (State of the Art Reviews), 15(4), 447–488. [CrossRef] [MathSciNet] [Google Scholar]
  10. Jameson A, Hu R, Wang Q. 2012. Adjoint-based aerodynamic optimization of supersonic biplane airfoils. Journal of Aircraft, 49(3), 802–814. [CrossRef] [Google Scholar]
  11. Courty F, Dervieux A, Dervieux B, Hascoet L. 2003. Hascoet, Reverse automatic differentiation for optimum design: from adjoint state assembly to gradient computation. Optimization Methods and Software, 18(5), 615–627. [CrossRef] [MathSciNet] [Google Scholar]
  12. Hascoet L. 2004. TAPENADE: A Tool for Automatic Differentiation of Programs. Proceedings of the ECCOMAS Conference, Jyvaskyla, Finland. [Google Scholar]
  13. Duvigneau R, Chandrashekar P. 2012. Kriging-based optimization applied to flow control. Int. J. for Numerical Methods in Fluids, 69(11), 1701–1714. [CrossRef] [Google Scholar]
  14. Marino E, Salvatori L, Orlando M, Borri C. 2016. Two shape parametrizations for structural optimization of triangular shells. Computers & Structures, 166, 1–10. [CrossRef] [Google Scholar]
  15. Marino E, Salvatori L, Orlando M, Borri C. 2016. Two shape parametrizations for structural optimization of triangular shells. Comput. Struct., 166, 1–10. [Google Scholar]
  16. Nocedal J, Wright SJ. 2006. Numerical Optimization, 2nd edn. Springer-Verlag: Berlin, Germany. [Google Scholar]
  17. Farin G. 1990. Curves and surfaces for computer-aided geometric design – A practical guide, Rheinboldt W, Siewiorek D, Editors. Academic Press, Boston. [Google Scholar]
  18. Désidéri J-A, Zolésio J-P. 2005. Inverse shape optimization problems and application to airfoils. Control and Cybernetics, 34(1), 165. [MathSciNet] [Google Scholar]
  19. Abou El Majd B. 2014. Parameterization adaption for 3D shape optimization in aerodynamics. International Journal of Science and Engineering, 6(1), 61–69. [Google Scholar]
  20. Nielsen HB. 2000. UCMINF – An Algorithm For Unconstrained, Nonlinear Optimization, Informatics and Mathematical Modelling (IMM). Technical University of Denmark. [Google Scholar]
  21. Wesseling P. 1992. An introduction to Multigrid Methods. John Wiley & Sons: Chichester. [Google Scholar]
  22. Zhao J, Abou El Majd B, Desideri JA. 2015. Two level correction algorithm for parametric shape inverse optimization. International Journal of Engineering and Mathematical Modelling, 2(1), 17–30. [Google Scholar]
  23. Zhao J, Abou El Majd B, Desideri JA. 2007. Two level correction algorithms for model problems. INRIA. [Google Scholar]
  24. Sederberg T, Parry S. 1986. Free-from deformation of solid geometric models. ACM SIGGRAPH Computer Graphics, 20(4), 151–160. [CrossRef] [Google Scholar]

Cite this article as: El Majd BA, Ouchetto O, Désidéri JA & Habbal A: Hessian transfer for multilevel and adaptive shape optimization. Int. J. Simul. Multisci. Des. Optim., 2017, 8, A9.

All Figures

thumbnail Figure 1.

Initial wing shape and mesh in the symmetry plane.
In the text
thumbnail Figure 2.

Comparison of the convergence history for the three strategies.
In the text
thumbnail Figure 3.

Illustration of and by using two embedded control polygons of degree n = 6 and N = 8. (a) Upward process, (b) downward process.
In the text
thumbnail Figure 4.

Basic algorithm for N = 12; Comparison between three strategies: full convergence, with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.
In the text
thumbnail Figure 5.

Multilevel algorithm for n = 4 and N = 12; Comparison between two strategies: with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.
In the text
thumbnail Figure 6.

Adaption algorithm for n = 8; Comparison between two strategies: with Hessian transfer, without Hessian transfer. (a) Iterative performance, (b) optimum shape.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.