Study of interior point methods and multiple solutions in topology optimization

Study of interior point methods and multiple solutions in topology optimization Document: Report Author: Antonio Aristizábal Fernández Director/Co-director: Alex Ferrer Ferre/Jose Antonio Torres Lerma Degree: Bachelor in Aerospace Technology Engineering Examination session: Spring 2023 Acknowledgements I would like to express my sincere gratitude to the director of this thesis, Dr. Àlex Ferrer, for giving me the opportunity to work with him and introducing me to this thrilling world of topology optimization, a truly innovative and interesting methodology. Also, I would like to thank the co-director, Jose Antonio Torres, for all the hard work that he has done through the development of this project and how he has encouraged me to keep it going. The journey that I have experienced during these past years far from my home would not have been possible without the support and understandment of both of my parents, Clara and Luis, that they have always shown to me. I would also like to thank my two sisters, María José and Juanita, for being always there when I needed them and my two nephews for showing me that happiness can be found in the most little things. Someone that has been unconditionally by my side from the most difficult times to the moments of joy, is my significant other Eugenia. I want you to know that this project and my whole bachelors degree have a big part of you. Finally, I would like to dedicate this thesis to my grandmother Leticia for being without exception an inspiring and loved human being. You will always have a place inside my heart. Abstract The seek of lightweight, stiff structures has always been present in the aerospace industry, a sector where mass reduction pays off as efficiency in all areas. Topology Optimization is a methodology that was born in the early 1980s and started to fully develop during the late 1990s. It is a numerical methods technique that tries to find the optimal material distribution inside a sealed domain, given a certain criteria and subjected to a set of constraints. This thesis aims to implement an interior point methods constrained optimizer, able to solve structural Topology Optimization problems. To accomplish this objective, the project will start with the refactory of an open source interior point algorithm, and later on it will be introduced into the SwanLab GitHub repository to fit the structure and characteristics of the different algorithms that had been already implemented in this Topology Optimization tool. Thus, the implemented algorithm will be validated with simple academic optimization problems, before testing it with proper structural Topology Optimization problems. Therefore, the achieved optimal structures will be discussed and compared with optimal solutions reached by other optimizers. Resumen La búsqueda de estructuras ligeras y rígidas siempre ha estado presente en la industria aeroespacial, un sector dónde la reducción de masa se traduce en eficiencia en todas sus áreas. La optimización topológica es una metodología que nace a principios de los 80 y comienza a desarrollarse plenamente a finales de los 90. Es una técnica de métodos numéricos que intenta encontrar la distribución de material óptima dentro de un dominio cerrado, dado un criterio específico y un conjunto de resticciones. Esta tesis busca implementar un optimizador con restricciones basado en los métodos de punto interior, capaz de resolver problemas de optimización topológica, y aplicado en el ámbito estructural. Para conseguir este objectivo, el proyecto comenzaŕa con la refactorización de un algoritmo código abierto de métodos de punto interior y posteriormente será introducido en el repositorio de GitHub SwanLab para seguir la estructura y características de los distintos algoritmos que ya han sido implementados dentro de esta herramienta de optimización topológica. El algoritmo implementado será validado con problemas de optimización académicos simples antes de ser testeado frente a problemas de optimización topológica. Seguidamente, las soluciones óptimas obtenidas serán analizadas y comparadas con soluciones óptimas obtenidas por otros optimizadores. Contents 1 Introduction 1 1.1 Scope .......................................... 1 1.1.1 Introductorystage ............................... 1 1.1.2 Interior Point Methods algorithm development . . . . . . . . . . . . . . . 2 1.2 Requirements...................................... 3 1.3 Justification....................................... 3 2 State of the art 5 2.1 Applications....................................... 6 2.2 HowdoesTOwork? .................................. 8 2.3 Existence of solutions in TO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Density and Level-set Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.1 Density ..................................... 12 2.4.1.1 SIMPmethod ............................ 13 2.4.1.2 SIMP-ALL method . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.2 Level-set..................................... 15 2.4.2.1 Level-set theoretical background . . . . . . . . . . . . . . . . . . 15 2.4.2.2 Application to TO . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5 Reduced formulation of the problem . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Unconstrained optimization problem . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.6.1 Unconstrained optimizers . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Constrained optimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7.1 Thedualfunction ............................... 24 2.7.2 Constrained optimizers: x, λ ......................... 26 3 Interior Point Methods 28 3.1 Problemdefinition ................................... 28 3.2 Barrierfunction..................................... 29 3.3 Centralpath ...................................... 31 3.3.1 Dual points from central path . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.2 Interpretation by KKT conditions . . . . . . . . . . . . . . . . . . . . . . 33 i 3.3.3 Force field interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Barriermethod..................................... 34 3.4.1 Accuracy of centering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4.2 Choice of τ................................... 35 3.4.3 Choice of t(0) .................................. 35 3.5 Solving the Barrier Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.5.1 Problemdefinition ............................... 36 3.5.2 KKT Conditions for barrier problem . . . . . . . . . . . . . . . . . . . . . 37 3.5.3 Barrier problem resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 Code validation 40 4.1 Test1 .......................................... 40 4.2 Test2 .......................................... 42 4.3 Test3 .......................................... 44 4.4 Test4 .......................................... 45 5 Results 48 5.1 2DCantileverBeam .................................. 49 5.1.1 Problemdefinition ............................... 50 5.1.2 Casesofstudy ................................. 50 5.1.3 Case1...................................... 50 5.1.4 Case2...................................... 53 5.2 2DArch......................................... 56 5.2.1 Problemdefinition ............................... 57 5.2.2 Casesofstudy ................................. 58 5.2.3 Case1...................................... 58 5.2.4 Case2...................................... 62 5.3 2DBridge........................................ 65 5.3.1 Problemdefinition ............................... 66 5.3.2 Casesofstudy ................................. 67 5.3.3 Case1...................................... 67 5.3.4 Case2...................................... 69 5.3.5 Case3...................................... 71 5.4 Optimal structures comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4.1 2DCantileverBeam .............................. 73 5.4.2 2DArch..................................... 76 5.4.3 2DBridge.................................... 78 6 Conclusions and further research 79 ii A Algorithm development 81 A.1 Refactory ........................................ 81 A.1.1 Firstrefactory ................................. 81 A.1.2 Secondrefactory ................................ 82 A.1.2.1 Quasi-Newton methods . . . . . . . . . . . . . . . . . . . . . . . 83 B UML diagrams 86 B.1 Firstrefactory...................................... 87 B.2 Finalrefactory ..................................... 88 Bibliography 89 iii List of Figures 2.1 Continuous (left) vs discrete (right) structures in TO [9]. . . . . . . 5 2.2 Generalized shape design problem of finding the optimal material distribution, figure adapted from [10]. . . . . . . . . . . . . . . . . . 6 2.3 Material reduction with similar stiffness capabilities. . . . . . . . . . 7 2.4 Optimal structures regarding its application field and constraints [13]. 7 2.5 TO design variable representation in the material domain. . . . . . . 8 2.6 TO domain definition. . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.7 Multiple feasible solutions reducing compliance. . . . . . . . . . . . . 10 2.8 Filtering approach to reach uniqueness of solutions in TO. . . . . . . 11 2.9 Domain problem illustration. . . . . . . . . . . . . . . . . . . . . . . 12 2.10 Young’s modulus representation as a function of density. . . . . . . . 14 2.11 Interpretation of density values with different microstructures. . . . 14 2.12 A 2D boundary representation of the zero level-set of a 3D level-set function[25]................................ 16 2.13 Design variable χevolution with the level-set function ψ. . . . . . . 17 2.14 Variation of Ωθof a shape Ω⊂Rdaccording to a deformation field θ[27].................................... 18 2.15 Cost variation by boundary translation. . . . . . . . . . . . . . . . . 18 2.16 Cost variation of a shape Ωby a nucleated hole B(x, r). . . . . . . . 19 3.1 Indicator function representation. . . . . . . . . . . . . . . . . . . . . 30 3.2 Indication function and approximations via barrier functions, adapted from[4]................................... 30 4.1 Validation test n.1, IPM algorithm results comparison. . . . . . . . . 41 4.2 History curves for optimization test 1. . . . . . . . . . . . . . . . . . 41 4.3 Validation test n.2, IPM algorithm results comparison. . . . . . . . . 42 4.4 History curves for optimization test 2. . . . . . . . . . . . . . . . . . 43 4.5 Validation test n.3, IPM algorithm compared to other implemented optimization algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . 44 iv 4.6 History curves for optimization test 3. . . . . . . . . . . . . . . . . . 45 4.7 Validation test n.4, IPM algorithm compared to other implemented optimization algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.8 History curves for optimization test 4. . . . . . . . . . . . . . . . . . 46 5.1 Geometric setting for one load 2D cantilever beam test case. . . . . 49 5.2 2D Cantilever beam TO iterative evolution, case 1. . . . . . . . . . . 50 5.3 History curves of 2D Cantilever Case 1, iterative evolution of the monitoring parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.4 2D Cantilever beam TO iterative evolution, case 2. . . . . . . . . . . 53 5.5 History curves of 2D Cantilever Case 2, iterative evolution of the monitoring parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.6 Generic setting for one load 2D Arch test case. . . . . . . . . . . . . 57 5.7 2D Arch TO iterative evolution, case 1. . . . . . . . . . . . . . . . . 59 5.8 History curves of 2D Arch Case 1, iterative evolution of the monitoring parameters................................. 62 5.9 2D Arch TO iterative evolution, case 2. . . . . . . . . . . . . . . . . 63 5.10 History curves of 2D Arch Case 2, iterative evolution of the monitoring parameters................................. 65 5.11 Geometric setting for the 2D Bridge for the multiple load case test case, adapted from [13] . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.12 2D Bridge TO iterative evolution, case 1. . . . . . . . . . . . . . . . 68 5.13 History curves of 2D Bridge Case 1, iterative evolution of the monitoring parameters................................. 69 5.14 History curves of 2D Bridge Case 2, iterative evolution of the monitoring parameters................................. 71 5.15 History curves of 2D Bridge Case 3, iterative evolution of the monitoring parameters................................. 72 5.16 History curves of 2D Cantilever Case 1 optimal solutions comparison, iterative evolution of the monitoring parameters. . . . . . . . . . . . 74 5.17 History curves of 2D Cantilever Case 2 optimal solutions comparison, iterative evolution of the monitoring parameters. . . . . . . . . . . . 75 5.18 History curves of 2D Arch Case 1 optimal solutions comparison, iterative evolution of the monitoring parameters. . . . . . . . . . . . 77 5.19 History curves of 2D Arch Case 2 optimal solutions comparison, iterative evolution of the monitoring parameters. . . . . . . . . . . . 78 v Chapter 2 State of the art Topology optimization is a relatively new branch of numerical methods referring to the process of finding the optimal distribution of materials within a given design area [5]. TO is the structural optimization method least sensitive to the initial design, so has the greatest capacity for realizing improved design for new production techniques such as additive manufacture [6]. This structural optimization branch started with the establishment of its principles for discrete structures, firstly for truss-like and grilled structures, later on, bar systems were introduced to the methodology in [7]. Continuous structures appeared for the first time in TO studies with the homogenization technique where the determination of the optimum parameters of holes in unit cells that constituted such structure was consolidated [8]. For the sake of this thesis, only TO applied to continuous structures will be studied. Figure 2.1: Continuous (left) vs discrete (right) structures in TO [9]. Continuum topology design problems are defined on a fixed reference domain Ωcontained in R2 or R3. In this reference domain, the goal is set to obtain the optimal distribution of material within a certain volume of mass constraint, with this material distribution contained in a domain Ωmat, part of Ω, on which applied loads and boundary conditions are defined, see Fgiure 2.2. 5 2.1. Applications 6 The term “optimal” is defined by the choice of objective and constraint functions, and through the choice of design parametrization. These chosen functions involve some kind of physical modelling that provides a measure of efficiency within the framework of a given area of applications [10]. In this thesis TO will be applied in the structural mechanics area although the methodology can take part in multiple fields such as thermodynamics, fluid dynamics, etc. Reference domain =⇒No material Material domain ΩΩmat Ω⊇Ωmat Figure 2.2: Generalized shape design problem of finding the optimal material distribution, figure adapted from [10]. Structural TO can be considered as a procedure for optimizing the topological arrangement of material Ωmat into the reference domain Ω, eliminating the material volume that is not needed. It can also be seen as the problem of finding the structural layout that best transfers specific loading conditions to supports [11]. In another words, its purpose is to find Ωmat ⊆Ωwith a certain perimeter, volume or mass restriction and certain boundary conditions such that the stiffness of the structure is maximized. 2.1 Applications In the structural mechanics field, TO is mostly applied on the reduction of material loosing as little stiffness as possible. For this reason, as it will be seen throughout this thesis, the typical TO problem aims to minimize the compliance subject to a volumetric of mass constraint. The compliance of a structure can be seen as a measure of the structural flexibility, so minimizing the compliance will maximize the stiffness of the structure. Specifically, the compliance is defined as c=f·u, where fis the external forces vector and u the displacements vector. It yields a measure for the internal energy. TO works on removing the material that is not performing any structural work, this can also be defined as dead material. Performing different holes in the structure where the material can be considered as dead may solve the problem, Figure 2.3 illustrates the procedure. 2.1. Applications 7 FF =⇒ Figure 2.3: Material reduction with similar stiffness capabilities. Although TO was originally developed for mechanical design problems, over the years the approach has been applied to virtually every multiphysics problem. The essential requirement for any problem to be solved via TO is that its associated Partial Differential Equation (PDE) can be reliably discretized and modelled using finite element, finite difference, finite volume or other volume discretization schemes. Some of the recently reviewed areas can be fluidics, thermofluidics, wave propagation, electro-magnetism, etc [12]. This new perspective of the use of numerical methods can be employed in order to innovate. With a wide range of fields to be applied, TO algorithms could obtain new type of designs not only improving the desired functionalities of the structures, but also providing innovative conceptual designs at the early stages of the design process without assuming any prior structural configuration, see Figure 2.4. (a) Cantilever beam subjected to torsion and traction loads (b) Heat conductor with volume minimization Figure 2.4: Optimal structures regarding its application field and constraints [13]. 2.2. How does TO work? 8 2.2 How does TO work? TO mathematical formulation, comes by the definition of a simple minimization problem, with its design variable constrained within the interval {0; 1}, where zero denotes no material. For instance, a generic formulation of the problem can be defined as                          min χJ(u(χ), χ) =⇒cost s.t A(χ)·u=b(χ) =⇒PDE constraint g(χ) = 0 =⇒equality constraint h(χ)≤0 =⇒inequality constraint χ={0; 1}=⇒box constraints .(2.1) Previous Partial Differential Equation (PDE) constraint may appear from a physical problem differential equation, in discretized form, when a physical value u(χ)participates in the cost. Recall that in the material domain, χvalues can be represented as shown in Figure 2.5 χ= 0 χ= 1 Void (No material) Material Figure 2.5: TO design variable representation in the material domain. With this problem formulation and the design variable representation in the material domain, the general TO initial domain can be defined as depicted in Figure 2.6. Ω = {x∈D;x= 1} D Figure 2.6: TO domain definition. This problem can be solved using material distribution methods in order to find the optimum layout of linearly elastic isotropic material in a structural system [11]. Therefore, the question that arises is how to reduce material in the domain Ωin order to minimize the cost accomplishing the imposed constraints. Some of the typical structural TO cost and constraints parameters are: 2.3. Existence of solutions in TO 9 Cost J(χ)to minimize: – Compliance – Perimeter – Stress – Eigenvalues →ω dddddddddddConstraint C(χ)to achieve: – Volume – Perimeter – Minimum length scale All of these parameters can be called shape functionals being any abstract functional defined on some class of subsests Rd. These shape functionals depend directly on the geometrical domains, i.e. Ω, and implicitly by means of solutions, i.e. u=u(Ω). For the simplicity of its formulation, the minimization problem can be redefined to accomplish the requirements of the field of application. Closing this section, an optimizer will be needed in order to solve the TO minimization problem, which, for the sake of this thesis, is redefined as an structural problem. 2.3 Existence of solutions in TO In its beginnings, the main problem that TO faced was the no existence of unique solution, as proving convergence to a global optimum is generally not possible, as an infinite number of topologies may exist which exhibit similar performance [14]. As an example, it can be proposed a problem where a structure with some applied forces will be subjected to a reduction of volume, then several solutions are proposed accomplishing the volume constraint Example 2.1 Compliance minimization subject to a volume reduction constraint.            min χJ(χ) s.t V(χ) = 0.5V0 χ∈ {0,1} .(2.2) In this Example 2.1, the goal will be to minimize the compliance subjected to a volume reduction constraint of half of its initial volume. 2.3. Existence of solutions in TO 10 Reduction of 50% in volume JJJJ Less Compliance >> Figure 2.7: Multiple feasible solutions reducing compliance. As shown in Figure 2.7, each of this solutions will be better than the last one, in terms of compliance minimization, so the optimal solution will be found in the infinity of the number of bars, so it will be theoretically impossible to reach these optimal solutions. To ensure existence of solutions, there are different remedies for this particular problematic all of them involve combining the general formulation of the TO minimization problem with some kind of perimeter constraint, gradient constraint or filtering techniques. A further discussion of perimeter constraint and filtering techniques can be found below. Perimeter constraint Going back to the previous example, problem (2.2) can be redefined in order to accomplish uniqueness of solutions, for instance a perimeter constraint will be added to the problem formulation            min χJ(χ) + αPer(χ) s.tV(χ) = 0.5V0 χ∈ {0,1} ,(2.3) where Per(χ)is a functional computing the perimeter of holes boundaries and αits penalty factor. The general idea behind this perimeter control is to constraint the lengths or areas of all inner and outer boundaries of the structural topology, not allowing solutions with a high perimeter value. This restriction method to obtain the existence of solutions was proven in [15]. Filtering techniques Sensitivity and density filters make use of weighting values that incorporate the sensitivities or densities values of neighboring elements. The values of elements within the weighting matrix are dependent upon a user-defined filter radius as well as the distances between elements centroids in each finite element discretization. For centroids outside the filter radius the weighting value is zero, for centroids inside the filter size the averaged solution will appear gray, not anymore black nor white., see Figure 2.8. 2.4. Density and Level-set Approaches 11 Filter size Figure 2.8: Filtering approach to reach uniqueness of solutions in TO. The obtained weighting value can be computed from the following expression wij =     1 |Ni|j∈Ni 0j /∈Ni ,(2.4) where Niis the set of elements within the filtering radius for the ith element. Filtering allows existence of solutions as it stops the smaller scales of the problem and creates an average of the yet obtained parameters in order to solve these smaller scales. 2.4 Density and Level-set Approaches Before starting, recall that the domain problem is defined with the design variable χ, Besides, ¯x will represent position in the domain. Problem 2.1 Domain problem in TO:                    min χ, u c(χ, u) = f·u(χ) s.t K(χ)·u=f V(χ)−V∗= 0 χ∈ {0,1} .(2.5) The problem features the minimization of the structural compliance, c(χ, u) = F·u(χ), subjected to the equilibrium relation between forces and displacements and a volume constraint. In this initial formulation of the problem, χis defined as shown in Figure 2.9. 2.4. Density and Level-set Approaches 12 χ(¯x) = 1 χ(¯x) = 0 x1 x2 χ(¯x) =      1x∈Ω 0x /∈Ω Figure 2.9: Domain problem illustration. Nevertheless, the computation of this problem is extremely difficult, since χis not differentiable in the interfaces between voids and material. It is crucial, since the optimizers use gradient based methods. Thus, the problem will be relaxed by different approaches such as Density or Level-set. 2.4.1 Density Density-based TO operates on a fixed mesh of finite element methods penalizing the mechanical properties of the element. In this approach the design variable of the problem will be the density, ρ. The appearance of each element of the finite element mesh will be denoted by a density value, it will appear on black representing ρ= 1 while the void parts will appear white representing ρ= 0. The relaxation of the problem relies on ρ(x)∈[0,1] with the values comprehended between the lower and upper bounds represented as gray regions. The elements within these gray regions will be the ones penalized. Now, note that ρis a continuous, smooth and differentiable design variable. The goal of this approach is to design a structure as stiff as possible such that the displacement that the force is going to produce is very small. The mathematical formulation of the general Density-based TO problem, with compliance minimization and a volume reduction constraint, now comes by                    min ρ, u f·u(ρ) s.t K(ρ)·u=f V(ρ) = V∗ 0≤ρ≤1 .(2.6) Note that, these kind of TO problems, in which knowing one of the variables is enough to obtain the remaining unknown ones, are known as PDE-constraint optimization problems. 2.4. Density and Level-set Approaches 13 The main formulation of any physical problem is the PDE or strong form. In particular, the problem shall be rewritten into its weak form, as weak formulation tends to be easily discretized, in order to solve it numerically by, usually, an algebraic system of equations. Elasticity problem strong form (PDE): ∇ · ϕ:∇Su=f. (2.7) Elasticity problem weak form: Z∇Sv:ϕ:∇Su=Zfv, (2.8) being va test function [3]. The difference between the strong form and the weak form is that having in a structure unknown values, the strong form tends to obtain the exact solution of this points needing as much functions as unknowns, while the weak form averages the results by the surrounding points with known solutions. For instance, taking into account the previous weak formulation of the elasticity problem, Equation (2.6) will be rewritten as                    min ρ, u l(u) s.ta(ρ, u, v) = l(v) V(ρ) = V∗ 0≤ρ≤1 .(2.9) Between the different Density-based TO approaches, homogenization method and SIMP method are the most popular ones [16–18]. For the sake of this thesis, only SIMP method and its variant SIMP-ALL will be further studied. 2.4.1.1 SIMP method The Solid Isotropic Material with Penalization (SIMP) model, also called penalized, porportional “ficticious material” model, is the most popular approach in Density-based TO. This approach, makes use of the continuous variable ρresembling the density of material by the fact that the volume of the structure is evaluated as V(ρ) = ZΩ ρ(¯x)dΩ.(2.10) The relation between the weak formulation of the problem and its dependency on the density of material, ρ, comes by the constitutive tensor. In this constitutive tensor, one can find mechanical 2.4. Density and Level-set Approaches 14 properties of the material such as the Young’s Modulus and the Poisson ratio, that directly rely on ρas C(ρ) = E(ρ)       1−ν(ρ) 0 −ν(ρ) 1 0 0 0 1−ν(ρ) 2       .(2.11) In this case, as gray solutions will appear (with the relaxation of the problem there will be no longer only black or white regions) the Young’s Modulus of these gray areas or elements shall as well be computed. For instance, from Equation (2.11), let νbe constant, then a direct relation between Young’s Modulus and ρcan be defined, see Figure 2.10. E ρ E(1) = E E(ρ) = ρpE with p= 3 Figure 2.10: Young’s modulus representation as a function of density. With this Young’s Modulus interpolation, it is not immediately apparent that gray areas achieved with SIMP method can be interpreted in physical terms. However, any stiffness used in the SIMP model can be realized as the stiffness of a composite made of void and based material in portions that correspond to the relevant density. Figure 2.11, illustrates different homogenized microstructures with an intermediate density value. Note that infinite microstructures are found for a specific gray value. ρ= 0.5 E E E Figure 2.11: Interpretation of density values with different microstructures. SIMP method, as its name suggests, penalizes the intermediate values of ρ, the homogenized microstructures, which will contribute to the structural integrity with a very low percentage of 2.6. Unconstrained optimization problem 21 Then, ∇xJ=∂J ∂x + v z}|{ ∂J ∂y A−1(x) | {z } p   −∂A(x) ∂x A−1(x)b(x) | {z } y +∂b(x) ∂x   =∂J ∂x +p−∂A(x) ∂x y+∂b(x) ∂x , (2.31) with A(x)p=∂J ∂y being the system to compute the called adjoint variable p. In order to solve the adjoint problem based on the previously defined gradient computation, the following algorithm is proposed. Note that the only outcome of the adjoint problem will be to redefine the cost gradient, while reducing the amount of design variables and eliminating/simplifying the PDE constraint in any TO problem. Algorithm 1 Gradient method algorithm 1: Solve A(x)y=b(state, primal) 2: Solve A(x)p=∂J(x, y) ∂y (adjoint) 3: Compute ∇xJ=∂J ∂x +p∂A(x) ∂x y+∂b(x) ∂x  4: Update x;xk+1 =xk−α∇xJ, with αs.t. J(xk−α∇xJ)< J(xk) 2.6 Unconstrained optimization problem The mathematical formulation of a typical unconstrained optimization problem usually takes the following form      min xJ(x)x∈Rn Optim. condition ∇J(x∗)=0 ∇f(x∗)∈Rn .(2.32) When solving these type of problems by an iterative process, the updated design variable can be computed via the steepest descend method as xn+1 =xn−t· ∇J(xn),(2.33) where t=line search The line search parameter has a great importance in the optimization problem. Taking a really low value of the line search will assure J(xn)> J(xn+1)but will not accomplish great 2.6. Unconstrained optimization problem 22 computational speed or performance. On the other hand, large values of tcannot assure the desired descend condition of the variable. Hence, an optimal value of the line search must be computed in order to assure both proper descent behaviour and computational speed. Unconstrained optimization problems can as well have box constraints, so then, Equation (2.32) can be rewritten as follows            min xJ(x), x ∈Rn 0≤x≤1 ∇J(x∗)∈Rn Optim. condition         ∇J(x∗)>0x= 1 ∇J(x∗) = 0 0 ≤x≤1 ∇J(x∗)<0x= 0 .(2.34) Then, the steepest descend function can be complemented by the box projection as xn+1 =xn−t· ∇J(xn) xn+1 = max(0,min(1, xn+1)) .(2.35) 2.6.1 Unconstrained optimizers This type of optimizers solve the primal problem using the gradient of a defined merit function inside a steepest descent like algorithm. In this section, a density-based optimizer will be presented, as well as a level-set based optimizer, both of them already implemented inside the SwanLab environment [1]. Projected Gradient: x=ρ As the design variable of projected gradient suggests, this optimizer was designed for the resolution of the density optimization problem. The optimizer starts from an initial density ρ0∈[0,1]. Then, it produces a sequence of ρn∈ [0,1], n = 0,..., going from ρnto ρn+1 following the steepest descent direction supplied by the merit function gradient ∇Φ. Thus thresholding in such way that ρn+1 ∈[0,1] Projected Gradient can be formulated as ρn+1 =max(0,min(1, ρn−τn∇Φn)),(2.36) where τnis the step length. Its initial value can be computed as τ0=||∇Φ|| ||∇ρ|| . 2.6. Unconstrained optimization problem 23 In order to accomplish convergence criteria between updated values, this step length can be iteratively decreased until convergence is achieved. Also, seeking an increment in the computational speed, this step length can also be increased when a converged value is obtained. The solution obtained as a density function by the optimizer may have intermediate density values in certain regions, even after being penalized by the chosen interpolation profile. A popular method to eliminate this gray regions consists of projecting the latter onto the set of functions taking only values 0 and 1 in order to interpret the solution as a true “black and white” shape [29]. Slerp: x=ψ Slerp optimizers in TO estimate the variation of the cost function by means of the topological derivative, as seen at Section 2.4.2.2. As it was previously proven, the optimality condition of this methodology leads to the topological derivative DTJ(χ)being equivalent to the level-set function ψ. Thus, a fixed point algorithm can be introduced in order to fulfill the optimality condition as    DTJ(χ) = ψ ||ψ|| = 1 ,(2.37) with the updated value of the level-set function being computed as ψn+1 =DTJ(χ(ψn)) 7−→ ψn+1 =ψn+1 ||ψn+1||.(2.38) As this approach brings inestability to the algorithm resolution, it was mixed with the original Slerp application of quaternion interpolation in order to track 3D rotation. This combination resulted in the following expression to update the level-set function value ψn+1 =αn(K)ψn+βn(K)·DTJ(χ(ψn)),(2.39) where αn, βntaken such that ||ψn+1|| = 1, and θn=cos(ψn, DTJ(χ(ψn))) ||ψn|| · ||DTJ(χ(ψn))||. 2.7. Constrained optimization problem 24 2.7 Constrained optimization problem In this section, constrained optimization problems including equality constraints will be introduced and discussed. Optimization problems including inequality constraints are extendedly and fully discussed in the following Chapter 3. An equality constrained optimization problem can be mathematically expressed as      min xJ(x) s.t gi(x)=0, i = 1, . . . , p ,(2.40) with x∈Rn. It is assumed that a domain D=Tp i=1 dom Ciis nonempty, and denoting the optimal design variable as x∗. To abroad Equation (2.40), the idea will be to introduce a new problem, called the dual problem, in order to take into account the constraints by augmenting the objective function with a weighted sum of the constraint functions. This can be achieved by means of the Lagrangian function defined as follows L(x, λ) = J(x) + p X i=1 λigi(x),(2.41) with dom L=D × Rp. In Equation (2.41), λican be referenced as the Lagrange multiplier of the ith equality constraint. Thus, vector λis known as the Lagrange multiplier vector or dual variable while xis known as the primal variable of the problem. 2.7.1 The dual function Given any λ, the Lagrange function can be redefined in order to fulfill the minimization of the objective function J(x)over xas the dual function D(λ) = min xL(x, λ) = min xJ(x) + p X i=1 λ+ igi(x).(2.42) It can be seen that maximizing the dual variable λwill maximize the weighted sum of the constraint functions, maximizing as well its contribution to the objective function. So then, minimizing the dual problem with respect to xwill result on obtaining the optimum value x∗as both the objective function and the maximized weighted sum of the constraints are minimized. 2.7. Constrained optimization problem 25 This can be expressed by max λD(λ) = max λmin xL(x, λ) = max λmin xJ(x) + p X i=1 λ+ iCi(x).(2.43) The optimality conditions for both primal and dual variables are commonly known as the Karush- Kuhn-Tucker (KKT) conditions and take the following form                ∂L ∂x =∂J ∂x + p X i=1 λi ∂gi ∂x = 0 (Primal opt. condition) ∂L ∂λi =gi(x)=0 (Dual opt. condition) .(2.44) For a better understandment of the dual problem, the following example is proposed where the Lagrange function and the dual function are further discussed. Example 2.3 Dual problem example. Recalling the Lagrangian function of the equality constrained optimization problem, Equation (2.41), let ˜xbe feasible and g(˜x)=0, such that L(˜x, λ) = J(˜x) +  : 0 λg (˜x) = J(˜x).(2.45) Then, the dual function can be formulated as D(λ) = min xL(x, λ)≤ L (˜x, λ) = J(˜x),(2.46) since D(λ)≤J(x∗)≤J(˜x)∀˜xfeasible. It also holds for x∗that satisfies min xJ(˜x), this is J(x∗) = min xJ(˜x)≤J(˜x). Thus, min xL(x, λ) | {z } D(λ) ≤ L (x∗, λ) | {z } J(x∗) ≤ L (˜x, λ) | {z } J(˜x) (2.47) Now it is easier to understand the importance of maximizing the dual problem, finding between all feasible solutions the one that gives minimum compliance enhancing the effect of this constraints on the objective function by maximizing its Lagrange multiplier λ. 2.7. Constrained optimization problem 26 Concluding this section, the previously defined adjoint problem, see Section 2.5, will be solved approaching the problem from a Lagrange multiplier interpretation. The Lagrange function of the adjoint problem can be defined as max pmin x, y L(x, y, p) = max pmin x, y J(x, y) + pT[A(x)y−b(x)].(2.48) Now the KKT conditions for this Lagrange function will result in                            ∂TL(x, y, p) = A(x)y−b(x) = 0 =⇒A(x)p=∂J ∂y (fix point) ∂yL(x, y, p) = ∂J ∂y +pTA(x) = 0 =⇒A(x)y=b(x)(fix point) ∂xL(x, y, p) = ∂J ∂y +pT∂A(x) ∂x y−∂b(x) ∂x = 0 =⇒(gradient) .(2.49) The proposed algorithm for the resolution of this approach is the previously defined gradient method: Algorithm 2 Gradient method algorithm for the interpretation as Lagrange multiplier 1: Solve A(x)y=b(x) 2: Solve A(x)p=∂J ∂y 3: Compute ∂xL(x, y, p) = ∂J ∂y +pT∂A(x) ∂x y−∂b(x) ∂x  4: Update x;xn+1 =xn−α∂xLwith αs.t. J(xn−α∇xJ)< J(xn) 2.7.2 Constrained optimizers: x, λ Constrained optimizers will be crucial for the resolution of TO problems, as the minimization of the cost function will be subjected to certain physical constraints. There are a wide range of different constrained optimizers, the vast majority of them rely on the dual problem definition and try to reach the optimality conditions for both primal and dual variables. Between the constrained optimizers, there are some of them already implemented in the SwanLab environment such as: 2.7. Constrained optimization problem 27 – Dual Nested in Primal – Augmented Lagrangian – Null Space [13] ddddddddddd– Fmincon – MMA – IPOPT The optimizers on the left column were either revised or implemented in a bachelor final thesis also supervised by the director of this project [30]. It can be consulted for a further explanation of those constrained optimizers, regarding how they find xn+1 and λn+1 t each new step, along with convergence criteria. Chapter 3 Interior Point Methods 3.1 Problem definition In this chapter, interior point methods for solving constrained optimization problems with inequality constraints will be further discussed. Consistent the following generic nonlinear constrained optimization problem,          min xf0(x) s.tfi(x)≤0, i = 1, . . . , m Ax =b ,(3.1) where f0, . . . , fm:X→Rare convex and twice continuously differentiable functionals and x∈Xis a design variable from a functional space X: Ω →Rn. Assuming that the problem is solvable, an optimal x∗exists, then the optimal value can be denoted as f0(x∗) = p∗. Supposing that the problem is strictly feasible, there is an x∈Xthat satisfies a PDE Ax =b and fi(x)<0for i= 1, . . . , m meaning that the Slater’s constraint is achieved so there a dual optimal λ∗∈Rm,ν∗∈Rpexists. Slater’s condition is a specific example of a constraint qualification. In particular, if Slater’s condition holds for the primal problem, then the duality gap is 0, and if the dual value is finite then it is attained [31]. These last two parameters, λ∗and ν∗, combined with x∗satisfy the KKT conditions considering inequality constraints, as 28 3.2. Barrier function 29                    Ax∗=b, fi(x∗)≤0, i = 1, . . . , m ∇f0(x∗) + m X i=1 λ∗ i∇fi(x∗) + ATν∗= 0 λ∗⪰0 λ∗ ifi(x∗)=0, i = 1, . . . , m .(3.2) Interior point methods can solve both, general form of the problem and KKT conditions by applying Newton’s method to a sequence of equality constrained problems or a sequence of modified KKT conditions. For these problems, the KKT conditions are a set of linear equations which can be solved analytically. Newton’s method can be understood as a technique for solving a linear equality constrained optimization problem, with twice differentiable objective, by reducing it to a sequence of linear equality constrained quadratic problems. Therefore, interior point methods solve optimization problems with linear equality and inequality constraints by reducing it to a sequence of linear equality constrained problems [4]. 3.2 Barrier function The goal to reach is to approximately formulate the inequality constrained problem as an equality constrained problem that Newton’s method can solve. Equation (3.1) must be rewritten in order to make the inequality constraints implicit in the objective as      min xf0(x) + Pm i=1 I_(fi(x)) s.tAx =b ,(3.3) where I_: R→Ris the indicator function for nonpositive reals it consist of a translation of the constraint into the own minimization equation of the problem. The indicator function aims to “penalize” as much as possible the objective function of the problem when the inequality constraints are not accomplished. In the feasible region, this indicator function will not affect the objective function value. Then the idea of the indicator function can be represented as shown in Figure 3.1. 3.2. Barrier function 30 I_(u) =      ∞u > 0 0u≤0 u > 0 I_(u) u < 0 ∞ 0 Figure 3.1: Indicator function representation. Note that Equation (3.3) has no inequality constraints, but its objective function is not differentiable in general, so Newton’s method cannot be applied in first instance. The solution to this problematic is called the Barrier method which consists in approximating the indicator function by the barrier function, which is defined as ˆ I_(u) = −1 tlog(−u),dom ˆ I_=−R++, where parameter t > 0sets the accuracy of the approximation. The barrier function, as well as the indicator function are convex and nondecreasing. Now, unlike the indicator function, the barrier function is differentiable and closed. Figure 3.2 depicts the difference between these functions. Figure 3.2: Indication function and approximations via barrier functions, adapted from [4]. The approximation to the indicator function becomes more and more accurate as the value of t is increased. Substituting the barrier function for the indicator function in Equation (3.3) gives 3.5. Solving the Barrier Problem 37 As seen before, the barrier function can be introduced into this problem transforming the inequality constraints into an extra term of the objective function as        min xf(x)−µ m X i=1 log(gi) s.t ci(x)=0, i = 1, . . . , p .(3.23) Then, the dual problem can be defined, introducing the dual variable λas D(λ) = min xL(x, λ) = ∇f(x)− p X i=1 λi∇ci(x)− z z }| { µ m X i=1 ∇log(gi)=0.(3.24) 3.5.2 KKT Conditions for barrier problem The necessary KKT conditions that will determine optimality of the achieved solution can be expressed as                    ∇f(x∗)− p X i=1 λi∇ci(x∗)−z= 0 p X i=1 ci(x∗) = 0 XZe −µe = 0 ,(3.25) where e=       1 . . . 1       . With the problem and the optimality conditions already defined, the resolution process can be discussed. 3.5.3 Barrier problem resolution KKT conditions already shown in Equation (3.25) can be rearranged to fit the Newton-Raphson algorithm standards as follows       Wk∇c(xk)−I ∇c(xk)T0 0 Zk0Xk             dx k dλ k dz k       =−       ∇f(xk) + ∇c(xk)λk−zk c(xk) XkZke−µje       ,(3.26) 3.5. Solving the Barrier Problem 38 where Wk=∇2 xxL(xk, yk, zk) = ∇2 xx f(xk) + c(xk)Tλk−zk, Zk=       z10 0 0...0 0 0 zn       , and Xk=       x10 0 0...0 0 0 xn       . This last system of equations can be as well rearranged in order to obtain a linear system of equations    Wk+Pk∇c(xk) ∇c(xk)T0      dx k dλ k   =−   ∇f(xk) + ∇c(xk)λk−zk c(xk)   ,(3.27) being Pk=X−1 kZk. With the computed search directions for both primal and dual variable, one can compute the search direction for the zvariable dz kas dz k=µkX−1 ke−zk−Xkdx k.(3.28) Once the different directions are found, the line search process can begin in order to obtain the updated primal and dual variables as well as the updated barrier term z xk+1 =xk+τkdx k,(3.29) λk+1 =λk+τkdλ k,(3.30) zk+1 =zk+τkdz k,(3.31) where τkis the step length. This updated value must be checked in order to analyze if the step is acceptable in terms of the minimization process, there are two popular approaches: – Decrease in merit function merit(x) = objective z}|{ f(x) + ν p X i=1 constr. z }| { |ci(x)|.(3.32) – Filter methods. 3.5. Solving the Barrier Problem 39 In the developed algorithm, a decrease in the merit function will be taken as criteria to check the acceptance of the step. Then, the solution will converge when KKT conditions are satisfied with a previously defined tolerance shown below max|∇f(x) + ∇c(x)λ−z| ≤ εtol max|c(x)| ≤ εtol max|XZe −µe| ≤ εtol .(3.33) Finally, decrease the logarithmic barrier parameter µwith µ=κµ, (3.34) where κ < 1. Now, the previous IPM algorithm 3 can be modified to feature the equations defined throughout this section. Algorithm 4 APM’s IPM algorithm [32] 1: Initialize x0, λ0, z0 2: while not converged, E(x, λ, z)≤εtol do 3: Compute search direction with Equations (3.27) and (3.28) 4: while not acceptableStep do 5: Update x, λ and zwith Equations (3.29), (3.30) and (3.31) 6: if merit(xk+1)<merit(xk)then 7: acceptableStep = true 8: else 9: τk=τk 2 10: end if 11: end while 12: Check convergence with Equation (3.33) 13: Decrease µ, with Equation (3.34) 14: end while The development of the IPM algorithm with its initial refactory to fit OOP and clean code standards as well as the final refactory process to implement it and merge it into the SwanLab GitHub environment can be consulted at Appendix A. Chapter 4 Code validation In the course of this chapter, the refactored optimizer will accomplish a preliminary validation phase with analytical mathematical optimization problems. The first two tests were extracted directly from the open source algorithm, being two of the multiple example problems contained in it [32]. Hence, as the original algorithm included these two problems, the solution achieved with the refactored optimizer will be compared with the solution given by the original one. 4.1 Test 1 The first example features the following optimization problem              min (x1,x2)∈R2J(x1, x2) := x2 1+ (2x2)2 s.t h1(x1, x2) := 2x1+x2−9≤0 h2(x1, x2) := x1+ 2x2−10 = 0 .(4.1) This problem is composed by an elliptic paraboloid as objective or cost function and two linear constraints, being one of them inequality and the remaining one equality. The validation test n.1 is intended for proving the proper refactory of the algorithm through its different stages, accomplishing the purpose of any refactory process, restructuring and reorganizing the algorithm without changing its functionality and external behaviour. Figure 4.1 illustrates the obtained results. 40 4.1. Test 1 41 Figure 4.1: Validation test n.1, IPM algorithm results comparison. As can be seen, the final refactory could not achieve the same design variable path to the optimum point. After some investigation, it was concluded that the only factor that contributed in this discrepancy on the paths followed to reach the solution, was the Hessian function. Thus, forcing the refactored algorithm to use the correct Hessian function introducing it as an input to the algorithm, resulted in the exact same design variable paths to the solution. Being the Hessian of the Cost function a key variable on the calculation of the direction to follow when performing the line search, it is comprehensible that on the first iteration the refactored algorithm can not reach the same design variable value, as the initial Hessian function is chosen to be the identity matrix. With the following figures, the behaviour of the refactored algorithm can be further discussed. (a) Cost function Jevolution (b) Constraints hviolation evolution Figure 4.2: History curves for optimization test 1. 4.2. Test 2 42 Figure 4.2 helps to have a better understatement on the behaviour of the algorithm. On the one hand, as seen in Figure 4.1, both the refactored algorithm with the exact Hessian and the original algorithm coincide in cost and constraints evolution. On the other hand, the final algorithm with the Hessian approximation has a similar behaviour. However, on the first iteration, the offset between cost and constraints values makes it clear that the Hessian function is crucial in the optimization process. It is worth mentioning that the issue with the Hessian approximation will be only seen on the first iteration, as upcoming iterations will have the proper Quasi-Newton methods approximation. 4.2 Test 2 This second validation test aims to prove the proper functionality of the refactored algorithm with the approximation of the Hessian function. The optimization problem is formulated as              min (x1,x2)∈R2J(x1, x2) := x2 1−2x1·x2+ 4x2 2 s.t h1(x1, x2) := x2−0.1x1+ 1 ≤0 h2(x1, x2) := 10x1−2x2+ 1 ≤0 .(4.2) It can be seen that the cost function is more complex than in the previous case and both inequality constraints are linear. Figure 4.3: Validation test n.2, IPM algorithm results comparison. 4.2. Test 2 43 Figure 4.3 proves again that the refactory of the algorithm was successful with the exact Hessian. Nevertheless, the objective of this test is to further study the algorithm’s behaviour with the Hessian approximation. The design variable evolution between exact and approximated Hessian cases differs notably in comparison with Test 1. In spite of that, Figure 4.4 shows an interesting behaviour of the developed algorithm. (a) Cost function Jevolution (b) Constraints hviolation evolution Figure 4.4: History curves for optimization test 2. Previous figure shows that the tendency on the cost and constraints evolution is highly followed by the developed algorithm. In this case, the initial Hessian approximation did not have a notable effect on the cost and constraints value. From these tow initial tests, it can be concluded that the refactory was successful, taking into account that the Hessian approximation is the only factor that makes it differ from the original algorithm. Also, the three versions of the algorithm reach convergence in the same or similar number of iterations. The remaining two test have been extracted from Florian Feppon’s PhD thesis [13]. During these tests the behaviour of the IPM algorithm will be compared with other available algorithms at SwanLab. In this case, the comparison will be done with fmincon optimizer, in particular with SQP and IPOPT variants. It is also worth to mention that from now on, IPM algorithm will make use of the Hessian matrix approximation. 4.3. Test 3 44 4.3 Test 3 The third optimization case starts from an unfeasible point, x0= (1.5,2.25), and has the cost gradient not aligned with the constraints. The problem is as              min (x1,x2)∈R2J(x1, x2) := (x1−1)2+ (x2−2)2 s.t h1(x1, x2) := 1 x1 −x2≤0 h2(x1, x2) := x1−x2−3≤0 .(4.3) In this optimization test, the objective function is an elliptic paraboloid and the constraints form a parabolic feasible domain. The obtained results are illustrated in the following Figure 4.5. Figure 4.5: Validation test n.3, IPM algorithm compared to other implemented optimization algorithms. The path followed by the three optimizers of analysis is quite similar, as all of them tend to reach the feasible region directly. IPM algorithm seems to follow a more direct path to the solution but this can be discussed with the history curves of the results. 4.4. Test 4 45 (a) Cost function Jevolution (b) Constraints hviolation evolution Figure 4.6: History curves for optimization test 3. As can be seen in Figure 4.6, both fmincon variants have an outstanding peak in cost at its first iteration, while IPM tends to be more stable. Nonetheless, it is notable how quickly SQP reaches the solution in comparison with both interior point optimizers. From the constraints violation evolution, it can be remarked the stability of IPM algorithm while the remaining ones have an abrupter evolution. It can also be highlighted the number of iterations needed by IPM optimizer to reach the optimal solution, being comparable to the other optimizers analyzed. 4.4 Test 4 This final test is the more complex one, with the feasible region shaped from the outer region of a parabola and a half-space. The test has a saturated inequality constraint that becomes inactive along the optimization path. This will highlight the relevance of the dual problem which detects when a saturated inequality constraint ceases to be saturated. The definition of the problem comes by              min (x1,x2)∈R2J(x1, x2) := (x1+ 3)2+x2 2 s.t h1(x1, x2) := x2−x2 1≤0 h2(x1, x2) := −x1−x2−2≤0 .(4.4) 4.4. Test 4 46 Figure 4.7: Validation test n.4, IPM algorithm compared to other implemented optimization algorithms. In this case, the path followed to the solution by IMP algorithm is notably different from the one followed by both fmincon variants. IPM tends to reach rapidly a minimization of the cost function and then accomplishes the constraints violation, while fmincon cases always stay in the feasible region to reach the optimum value. (a) Cost function Jevolution (b) Constraints hviolation evolution Figure 4.8: History curves for optimization test 4. Figure 4.7 shows some interesting results regarding cost and constraints evolution as the behaviour of the IPM algorithm is remarkably different in both cases. Starting with Figure 4.8a, it can be seen that the path followed by the three optimizers is very similar and have the same tendency. 5.1. 2D Cantilever Beam 53 variations such as the value of µparameter. In this benchmark test, the three values of µstudied present a similar behaviour in terms of cost evolution and an almost exact constraint violation evolution. In the case of the dual variable evolution, the initial values of this parameter for µ= 10−8standout for being really high in comparison with the other two cases. However, with a few iterations it is able to reach a similar magnitude. It is worth to mention that, as these solutions were obtained by a Density-based unconstrained optimizer as it is Projected gradient, the results had to been post-processed in order to project the majority of the gray elements into black or white, making these solutions an accurate approximation of the real obtained structure. 5.1.4 Case 2 The only difference that this problem presents in front of the previous benchmark case is the final fraction of volume. Thus, Equation (5.3) takes the following form:      min χc(χ) = f·u(χ) s.t V(χ) = 0.4V0 (5.5) For this second benchmark test, Figure 5.4 shows a similar structural evolution as the previous case, see Figure 5.2. (a) Iteration 0 (b) Iteration 60 (c) Iteration 160 Figure 5.4: 2D Cantilever beam TO iterative evolution, case 2. In this case, as a lower final volume is required it can be seen how the upper and lower sides of the cantilever are turning thinner than in the previous case, also more lightgray areas appear by iteration 60, Figure 5.4b. From this last one to Figure 5.4c, the final material distribution is starting to show, with some gray areas to be replaced by white or black elements. 5.1. 2D Cantilever Beam 54 Table 5.2: Results for 2D Cantilever Beam Case 2, 80 ×40 mesh. 2D Cantilever Beam Results - Case 2 µOptimized structure 10−7 Volume 0.4V0 Compliance 2.2687 Iterations 313 10−8 Volume 0.4V0 Compliance 2.2621 Iterations 186 10−9 Volume 0.4V0 Compliance 2.2552 Iterations 181 5.1. 2D Cantilever Beam 55 10−10 Volume 0.4V0 Compliance 2.2630 Iterations 267 10−11 Volume 0.4V0 Compliance 2.2497 Iterations 616 In this second case of the 2D Cantilever Beam, a wider range of µvalues were able to reach a solution, in particular values from µ= 10−7up to µ= 10−11. Table 5.2, shows that every single µvalue tries to reach the same geometry accomplishing the imposed constraints. The only difference between the reached optimal solutions is the number of iterations and the final compliance value, although it can be seen that every solution has a similar compliance tending to minimize this value around 2.25. Figure 5.5 illustrates the iterative evolution of some of the monitoring parameters for µ={10−8,10−9and10−11}, being the values that achieved minimum compliance. 5.2. 2D Arch 56 Figure 5.5: History curves of 2D Cantilever Case 2, iterative evolution of the monitoring parameters. The three cases analyzed in this last figure present similar behaviour in both cost and constraints violation evolution, in the case of the dual variable evolution it is notable that this parameter takes the same evolutionary profile than in the previous benchmark test with an outstanding high value in the case of µ= 10−8, that rapidly follows the same tendence than the other two cases. Optimum compliance, and therefore the optimum solution, was obtained by µ= 10−11 in 616 iterations. However, the case of µ= 10−9obtained the second lowest compliance in only 181 iterations comprehending a really interesting optimal solution taking into account the computational speed and the efficiency of the minimization. 5.2 2D Arch The second benchmark test consists of a 2D arch clamped at its lower ends as can be seen in the Figure 5.6. This TO test has been extracted from [13] and can be originally found in [33]. 5.2. 2D Arch 57 ΓD ΓDΓ0 g Figure 5.6: Generic setting for one load 2D Arch test case. The 2D arch, Ωmat, is contained in the domain Ωmat ⊆Ω⊂R2with dimensions 2×1. As the previous benchmark case, the boundary ∂Ωis conformed by disjoints regions as ∂Ω=Γ∪ΓD∪Γ0(5.6) where •ΓDregions correspond to the Dirichlet boundary conditions. The structure Ωis clamped at this regions, made of its lower ends. The segments measure is about 10% of the structure’s length. •Γ0region is the Neumann boundary condition. This subset will be non-optimizable, it is the part of the structure Ωwhere a unit traction force gis applied in the downwards direction, uniformly distributed on the center bottom side of the arch as Γ0= 0.45xmin ≤x≤0.3xmax,(5.7) being xmin = 0 and xmax = 1. •Γis traction free and is the only region of the boundary ∂Ωsubject to optimization. 5.2.1 Problem definition For the 2D arch, two different optimization problems will be considered. In the first place, a simple compliance minimization with a volume constraint will be optimized by the algorithm, so then 5.2. 2D Arch 58 PArch1=     min χc(χ) = f·u(χ) s.t V(χ) = V∗ ,(5.8) where V∗is the final volume. This parameter will depend on each case of study. Then, the second problem will have a different approach, as, in this case, the volume will be the parameter subject to minimization with a compliance constraint, the problem can be defined as PArch2=     min χV(χ) s.t c(χ) = c∗ ,(5.9) where c∗is the final compliance. This parameter will be the result of Equation (5.8), as reference. 5.2.2 Cases of study Case 1. The first case of study will feature Equation (5.8) with a target volume of 20% of the structure initial volume. Case 2. This second case will feature Equation(5.9), for the proper formulation of the compliance constraint the solution of Case 1 will be needed. These two cases will also reach interesting results as they can be compared between each other. The behaviour of the algorithm can be further studied with the solutions as the formulation of the problem is the opposite between both cases. 5.2.3 Case 1 The mathematical formulation of this first arch-like benchmark case comes by      min χc(χ) = f·u(χ) s.t V(χ)=0.2V0 (5.10) Figure 5.7 illustrates the evolution of the material domain through different phases of the optimization process. 5.2. 2D Arch 59 (a) Iteration 0 (b) Iteration 80 (c) Iteration 140 Figure 5.7: 2D Arch TO iterative evolution, case 1. The iterative evolution of this arch-like benchmark test starts with Figure 5.7a representing the initial material domain among the Dirichlet and Neumann boundary conditions. From this initial state to Figure 5.7b the material distribution has been reduced until finding the basics of the optimal structure but mainly constituted by gray areas. Finally, the last iteration shown, iteration 140 in Figure 5.7c, has gotten rid of these gray regions leaving only a black and white structure. Table 5.3: Results of 2D Arch Case 1, 80 ×40 mesh. 2D Arch Results - Case 1 µOptimized structure 10−7 Volume 0.2V0 Compliance 4.6872 Iterations 244 10−8 Volume 0.2V0 Compliance 5.1763 Iterations 1000 5.2. 2D Arch 60 10−9 Volume 0.2V0 Compliance 4.8999 Iterations 1000 10−10 Volume 0.2V0 Compliance 4.6062 Iterations 884 10−11 Volume 0.2V0 Compliance 4.7454 Iterations 658 5.2. 2D Arch 61 10−12 Volume 0.2V0 Compliance 5.2585 Iterations 1000 This first case regarding the Arch configuration is the most extensive one in terms of µvalues to analyze, going from µ= 10−7to µ= 10−12. Nonetheless, not every µvalue was able to converge to an optimal solution and reached the maximum number of iterations imposed, 1000. In spite of that, it was considered that because of the variety of optimized geometries it was worth to include them all. Between the presented solutions, the ones that were not able to reach optimality are µ= 10−8, µ = 10−9and µ= 10−12. The solution corresponding to these values of µfeature a taller solution and their respective compliance highlights between the rest of the solutions reached. On the other hand, the monitoring parameters iterative evolution of the µvalues that were able to meet the optimality conditions are shown in Figure 5.8. 5.2. 2D Arch 62 Figure 5.8: History curves of 2D Arch Case 1, iterative evolution of the monitoring parameters. This case presents notable differences between µvalues in the cost and constraints iterative evolution in comparison with the previous cases of analysis. Different peaks can be seen in the cost evolution by µ={10−10,10−11}. Also, in the constraints violation evolution instabilities are notable for µ= 10−10, with the other two values of the µparameter following a smoother path. Finally, once again the lowest value of µpresents a initial values of the dual variable various magnitude orders higher than the rest. In this particular problem, the lowest compliance is achieved in 884 iterations with a value of µ= 10−10 presenting a compliance value that widely differs from the other optimal solutions. 5.2.4 Case 2 The formulation of the problem for this second case, see Equation (5.9), was dependent of the optimum compliance reached in the previous case. After obtaining and discussing the solutions, it was determined that the optimum compliance, and therefore, the compliance constraint to achieve in this problem will be of c∗= 4.60. This will be the first volume minimization subjected to compliance constraint problem that the optimizer will face. 5.3. 2D Bridge 69 Figure 5.13: History curves of 2D Bridge Case 1, iterative evolution of the monitoring parameters. As the complexity of the problem escalated in these bridge-like benchmark tests, an optimal solution could only be obtained with one value of µas can be seen in Table 5.5. The monitoring parameters iterative evolution shown in Figure 5.13, illustrates the complexity of the problem. The cost evolution starts with a huge minimization, this aggressive approach comes with a huge peak in the constraint violation. It can be seen that the cost starts to arise until reaching rapidly the lowest cost value and therefore the highest constraint violation. From there on, the constraint violation starts to descend until reaching a value of 0 along with the final optimal volume. 5.3.4 Case 2 The remaining two benchmark cases will not be presented with a structural iterative evolution because the applied loads occupy most of the structure’s surface, thus, it will be difficult or impossible to differentiate the states of the optimization process between iterations. The target compliance c∗for this particular problem was set to c∗= 1.45 after a compliance minimization subjected to a volume constraint of 50%was performed to the structure with the specific loads of this case of study. 5.3. 2D Bridge 70      min χV(χ) s.t c(χ)≤c∗(gi) = 1.45,for i= 3,4and 5 .(5.17) In this benchmark case, two values of µcould reach the optimality conditions. The obtained optimal structures can be seen at Table 5.6 and the history curves of some monitoring parameters can be consulted below in Figure 5.14. Table 5.6: Results of 2D Bridge Case 2, 200 ×10 mesh. 2D Bridge Results - Case 2 µOptimized structure 10−4Volume 0.4482V0 Compliance 1.45 Iterations 793 10−5Volume 0.4645V0 Compliance 1.45 Iterations 207 5.3. 2D Bridge 71 Figure 5.14: History curves of 2D Bridge Case 2, iterative evolution of the monitoring parameters. In the previous figure, one can appreciate how the behaviour of the optimizer with both parameters is similar. The cost evolution in both cases features a huge initial minimization being more protruding in the case of µ= 10−5. Analyzing the constraints violation evolution, both cases present a similar and smooth path. The optimal solution, in terms of volume minimization, corresponds to the case of µ= 10−4. Although, looking at the iterations needed to reach convergence µ= 10−5case highlights needing almost a third of iterations than the other case. 5.3.5 Case 3 In this case, the compliance to reach in the volume minimization will be set from a compliance minimization subjected to a volume reduction of 50%applying every single load represented in Figure 5.11. The achieved compliance value was c∗= 1.55.      min χV(χ) s.t c(χ)≤c∗(gi)=1.55,for i= 0,...,8 .(5.18) 5.3. 2D Bridge 72 Table 5.7: Results of 2D Bridge Case 3, 200 ×10 mesh. 2D Bridge Results - Case 3 µOptimized structure 10−5Volume 0.5014V0 Compliance 1.55 Iterations 216 Figure 5.15: History curves of 2D Bridge Case 3, iterative evolution of the monitoring parameters. Once again, only a value of µ= 10−5was able to reach an optimal solution. This optimal structure can be seen in Table 5.7. Also, the history curves of the monitoring parameters can be consulted in Figure 5.15, where a smooth and quick optimization is followed by the IPM optimizer as seen in cost, constraint violation and dual variable iterative evolution. 5.4. Optimal structures comparison 73 5.4 Optimal structures comparison In this section, the optimal solutions reached with the IPM algorithm will be compared with the optimal solutions obtained in [30]. The idea will be to select from each TO problem the structure that was able to reach the minimum cost value depending on the case of study. Regarding the implemented IPM algorithm, it will vary by means of the selected µ. On the other hand, as in [30] different optimization algorithms were studied, the optimal solution will vary in terms of the selected optimizer. Only solutions obtained with Projected Gradient unconstrained optimizer will be extracted from [30]. As the TO problems are the same, the optimizer that obtained the lowest cost will be employed to perform that specific benchmark test with the same mesh and load distribution to obtain a fair comparison between almost every Density-based optimizer in SwanLab. 5.4.1 2D Cantilever Beam Case 1. For the 2D Cantilever beam case, with a mass reduction of 40%the optimum result was obtained with an initial µvalue of 10−8by IPM algorithm. On the other hand, as detailed in [30], for this same benchmark TO problem the optimumm solution was reached with Bisection optimizer, also known as Dual Nested in Primal. Table 5.8: Optimal results of 2D Cantilever Beam case 1, 80 ×40 mesh. 2D Cantilever Beam Optimal Results - Case 1 IPM Bisection Volume 0.6V0 Compliance 1.4631 Iterations 414 Volume 0.6V0 Compliance 1.4658 Iterations 58 5.4. Optimal structures comparison 74 Figure 5.16: History curves of 2D Cantilever Case 1 optimal solutions comparison, iterative evolution of the monitoring parameters. Analyzing the results contained in Table 5.8 and in Figure 5.16, it can be seen that the IPM optimal solution can be considered better in terms of minimization of the cost, in this case the compliance of the structure. Also, IPM algorithm returns a more symmetrical solution, regarding manufacturing it can also be considered as a better solution as its production process is simplified. Nonetheless, Bisection optimizer presents an outstanding computational speed reaching the optimal solution in only 58 iterations, this can be crucial when solving more complex TO problems with a large number of elements. Case 2. This same TO problem was studied in [30]. There, it was found that the best minimization of the compliance was achieved by Augmented Lagrangian (AL) optimizer. Then, the comparison of optimal solutions can be performed between AL and IPM. Regarding the solutions obtained by IPM, the selected optimal structure did not present the lowest compliance value. However, due to the small difference in compliance and the huge gap between number of iterations the solution achieved by µ= 10−9was preferred in front of the µ= 10−11 solution. 5.4. Optimal structures comparison 75 Table 5.9: Optimal results of 2D Cantilever Beam case 2, 80 ×40 mesh. 2D Cantilever Beam Optimal Results - Case 2 IPM Augmented Lagrangian Volume 0.4V0 Compliance 2.2552 Iterations 181 Volume 0.4V0 Compliance 2.2368 Iterations 432 Figure 5.17: History curves of 2D Cantilever Case 2 optimal solutions comparison, iterative evolution of the monitoring parameters. Table 5.9 shows a comparison of optimal results of the 2D Cantilever Beam benchmark problem with a reduction of volume of 60%. Attending at the obtained results, it can be seen that the material distribution is almost the same for both IPM and AL optimizers. From Figure 5.17 it can be seen that AL algorithm presents more instabilities than IPM in terms of cost and constraints violation evolution, also this optimizer was capable to reach a lower compliance value 5.4. Optimal structures comparison 76 but, in spite of that, IPM needed less iterations to converge. In this case, it can be concluded that both optimizers have reached an optimal solution, IPM achieved an acceptable compliance minimization in less iterations while AL accomplished the lowest compliance. 5.4.2 2D Arch Case 1. The first TO benchmark problem of the arch case required a minimization of the compliance of the structure subjected to a volume reduction of 80%. As before, this benchmark test was also performed in [30], in that thesis, Bisection was the optimizer that reached the best solution in terms of minimizing the objective function of the problem. Table 5.10: Optimal results of 2D Arch case 1, 80 ×40 mesh. 2D Arch Optimal Results - Case 1 IPM Bisection Volume 0.2V0 Compliance 4.6062 Iterations 884 Volume 0.2V0 Compliance 4.6125 Iterations 86 Comparing the optimal structures, Table 5.10 shows both IPM and Bisection solutions. It can be seen that both optimizers tend to reach the same geometry, then comparing this shape with the ones shown in Table 5.3, it can be concluded that this is the optimal material layout for this particular problem. Also from Figure 5.18 a comparison between these optimizers can be made with cost, constraints violation and dual variable iterative evolution. IPM algorithm was able to reach a lower compliance value than Bisection, however it took IPM more than ten times the iterations that Bisection optimizer needed to reach the optimal solution. From these results and the ones collected in Table 5.8 it can be concluded that Bisection method is really efficient taking about computational speed. 5.4. Optimal structures comparison 77 Figure 5.18: History curves of 2D Arch Case 1 optimal solutions comparison, iterative evolution of the monitoring parameters. Case 2. This benchmark test had the best performance in terms of volume minimization and computational speed with the Method of Moving Asymptotes (MMA) optimizer, as detailed in [30]. From the analysis carried out throughout this thesis, the optimal solution was achieved by the IPM algorithm with a value of µ= 10−4. Both solutions can be consulted in Table 5.11. Table 5.11: Optimal results of 2D Arch case 2, 80 ×40 mesh. 2D Arch Optimal Results - Case 2 IPM MMA Volume 0.2127V0 Compliance 4.6 Iterations 292 Volume 0.2051V0 Compliance 4.6 Iterations 1500 5.4. Optimal structures comparison 78 Figure 5.19: History curves of 2D Arch Case 2 optimal solutions comparison, iterative evolution of the monitoring parameters. From the previous table and from Figure 5.19, it can be seen how the solution obtained with MMA algorithm could not reach convergence as the process was stopped by the maximum number of iterations. In spite of that, this optimizer was able to generate an interesting structure. Similarities between the obtained geometries are clear, being more unsymmetrical the solution provided by MMA. Also, it can be seen how to optimization process is smoother by this algorithm, as IPM minimizes the cost abruptly. 5.4.3 2D Bridge For the three different benchmark cases that conformed the bridge-like TO problem, no solution could be reached with other optimizers in order to perform a comparison between the IPM optimal solutions and other optimizer optimal solution. A further study of other optimizers can be performed to understand their behaviour in front of these kind of problems. As could be seen in the previous sections, when the volume is minimized subjected to a target compliance, the complexity of the TO problem increases and the behaviour of the optimizer completely changes in comparison to a case where the compliance is minimized subjected to a volume reduction. A.1. Refactory 85 For the first step, an approximate initial value B0will be needed, an initial estimation can be computed as B0=βI, (A.6) for instance this initial value will be assumed as B0=I, as there is no strategy for the selection of β[36]. The UML diagram of the second refactory can be consulted at Appendix B.2, where the difference on number of classes and simplicity of the algorithm is notable compared to the first refactory case. Appendix B UML diagrams 86 B.1. First refactory 87 B.1 First refactory B.2. Final refactory 88 B.2 Final refactory Bibliography [1] A. Ferrer Ferre, “Swan - Topology Optimization Laboratory,” Available: https://github. com/SwanLab/Swan, 2023. [2] T. Creus and J. A. Torres, “The student’s guide to clean code development,” https://github. com/SwanLab/Swan/wiki/The-student’s-guide-to-clean-code-development/_history, 2023. [3] A. Ferrer Ferre, “Lessons in Topology Optimization,” Available: https://sites.google.com/ view/alexferrer/teaching?authuser=0, 2022. [4] S. Boyd and L. Vandenberghe, Convex optimization, 7th ed. Cambridge University Press, 2004, pp. 555–777. [5] Z. Zhao, R. Yujian, L. Yongming, L. Zhibo, and T. Zhijian, “Topology Optimization of Continuum Structures Based on Binary Hunter-Prey Optimization Algorithm,” Symmetry, vol. 15, 2023. [6] A. Aremu, I. Ashcroft, R. Wildman, R. Hague, C. Tuck, and D. Brackett, “A hybrid algorithm for topology optimization of additive manufactured structures,” 22nd Annual International Solid Freeform Fabrication Symposium - An Additive Manufacturing Conference, SFF 2011, pp. 279–289, 2011. [7] G. Rozvany, “Aims, scope, methods, history and unified terminology of computeraided topology optimization in structural mechanics,” Structural and Multidisciplinary Optimization, vol. 21, pp. 90–108, 2001. [8] M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Computer Methods in Applied Mechanics and Engineering, vol. 71, pp. 197–224, 1988. [9] D. Budzyn, H. Zare-Behtash, A. Cowley, and A. Cammarano, “Topology optimization of compliant mechanisms as a design method to improve hardware performance in lunar dust environment,” in 19th European Space Mechanisms and Tribology Symposium, 2021. 89 Bibliography 90 [10] M. P. Bendsøe and O. Sigmund, “Material interpolation schemes in topology optimization,” Archive of Applied Mechanics, vol. 69, pp. 635–654, 1999. [11] G. Kazakis, I. Kanellopoulos, S. Sotiropoulos, and N. Lagaros, “Topology optimization aided structural design: Interpretation, computational aspects and 3d printing,” Heliyon, vol. 3, 2017. [12] O. Sigmund, “EML webinar overview: Topology Optimization - Status and Perspectives,” Extreme Mechanics Letters, vol. 39, 2020. [13] F. Feppon, “Shape and topology optimization of multiphysics systems,” Ph.D. dissertation, Thèse de doctorat de l’Universitè Paris-Saclay prèparèe à l’Ecole polytechnique, 2019. [14] J. A. Morelli, “Evaluation of Topology Optimization Filtering with Numeric Examples,” in Virginia Tech, Technical Reports, Aerospace and Ocean Engineering, 2019. [15] L. Ambrosio and G. Buttazzo, “An optimal design problem with perimeter penalization,” Calculus of Variations, vol. 1, pp. 55–69, 03 1993. [16] M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Computer Methods in Applied Mechanics and Engineering, vol. 71, no. 2, pp. 197–224, 1988. [Online]. Available: https: //www.sciencedirect.com/science/article/pii/0045782588900862 [17] M. Bendsoe and O. Sigmund, Topology Optimization - Theory, Methods, and Applications, 2nd ed. Berlin-Heidelberg: Springer-Verlag, 2004. [18] G. Rozvany, M. Zhou, and T. Birker, “Generalized shape optimization without homogenization,” Structural optimization, vol. 4, pp. 250–252, 1992. [19] A. Ferrer Ferre, “Simp-all: A generalized SIMP method based on the topological derivative concept,” International Journal for Numerical Methods in Engineering, vol. 120, no. 3, pp. 361–381, 2019. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme. 6140 [20] S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988. [21] J. A. Sethian, Level Set Methods and Fast Marching Methods, 2nd ed. Cambridge Monographs on Applied Computational Mathematics, 1999. [22] S. Osher and R. Fedkiw, The Level Set Methods and Dynamic Implicit Surfaces. Springer- Verlag, 05 2004, vol. 57, pp. xiv+273. Bibliography 91 [23] J. Sethian and A. Wiegmann, “Structural boundary design via level set and immersed interface methods,” Journal of Computational Physics, vol. 163, no. 2, pp. 489–528, 2000. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999100965811 [24] D. Herrero, “Level Set Method Applied to Topology Optimization,” Universidad Politécnica de Cartagena. Available: https://www.upct.es/goe/level-set-slides.pdf, 2012. [25] S. Chen, S. Gonella, W. Chen, and W. K. Liu, “A level set approach for optimal design of smart energy harvesters,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 37, pp. 2532–2543, 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045782510001258 [26] G. Allaire, C. Dapogny, and F. Jouve, “Chapter 1 - shape and topology optimization,” in Geometric Partial Differential Equations - Part II, ser. Handbook of Numerical Analysis, A. Bonito and R. H. Nochetto, Eds. Elsevier, 2021, vol. 22, pp. 1–132. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1570865920300181 [27] Éric Bonnetier and C. Dapogny, “An introduction to shape and topology optimization,” Université Grenoble-Alpes. Available: https://membres-ljk.imag.fr/Charles.Dapogny/ coursoptim/slides/PartIII_Hadamard_Method.pdf, 2020. [28] A. Novotny and J. Sokolowski, Topological Derivatives in Shape Optimization. Interaction of Mechanics and Mathematics, 01 2013. [29] S. Amstutz, C. Dapogny, and A. Ferrer Ferre, “A consistent relaxation of optimal design problems for coupling shape and topological derivatives,” Numerische Mathematik, vol. 140, 09 2018. [30] A. Pena, “Study of optimization algorithms for lightweight structures,” Bachelor Final Thesis, Universitat Politècnica de Catalunya. Available: https://upcommons.upc.edu/ handle/2117/374521, 2022. [31] J. J. Borwein and A. Lewis, Convex Analysis and Nonlinear Optimization: Theory and Examples. CMS Books in Mathematics (CMSBM), 2005. [32] APMonitor, “Interior point methods,” Available: http://apmonitor.com/me575/index.php/ Main/InteriorPointMethod, 2022. [33] G. Allaire, F. Jouve, and G. Michailidis, “Thickness control in structural optimization via a level set method,” Structural and Multidisciplinary Optimization, vol. 53, pp. 1349–1382, 2016. [34] Wikipedia, “Code refactoring,” Available: https://en.wikipedia.org/wiki/Code_refactoring, 2023. Bibliography 92 [35] L. Bottou, F. E. Curtis, and J. Nocedal, “Optimization Methods for Large-Scale Machine Learning,” 2018. [36] Wikipedia, “Quasi-Newton method,” Available: https://en.wikipedia.org/wiki/ Quasi-Newton_method, 2023.