A Hilbertian projection method for constrained level set-based topology optimisation

We present an extension of the projection method proposed by Challis et al. (Int J Solids Struct 45(14–15):4130–4146, 2008) for constrained level set-based topology optimisation that harnesses the Hilbertian velocity extension-regularisation framework. Our Hilbertian projection method chooses a normal velocity for the level set function as a linear combination of (1) an orthogonal projection operator applied to the extended optimisation objective shape sensitivity and (2) a weighted sum of orthogonal basis functions for the extended constraint shape sensitivities. This combination aims for the best possible first-order improvement of the optimisation objective in addition to first-order improvement of the constraints. Our formulation utilising basis orthogonalisation naturally handles linearly dependent constraint shape sensitivities. Furthermore, use of the Hilbertian extension-regularisation framework ensures that the resulting normal velocity is extended away from the boundary and enriched with additional regularity. Our approach is generally applicable to any topology optimisation problem to be solved in the level set framework. We consider several benchmark constrained microstructure optimisation problems and demonstrate that our method is effective with little-to-no parameter tuning. We also find that our method performs well when compared to a Hilbertian sequential linear programming method.


Introduction
The field of topology optimisation has enjoyed rapid growth owing to improved computing power, new optimisation techniques, and application to a wide range of design problems (Bendsøe and Sigmund, 2004;Deaton and Grandhi, 2013;Sigmund and Maute, 2013).Classical computational methods for topology optimisation include densitybased methods (Bendsøe, 1989;Rozvany et al, 1992) in which the design variables are material densities of elements/nodes in a mesh, and level set-based methods (Wang et al, 2003;Allaire et al, 2004) in which the boundary of the shape is implicitly tracked as the zero level set of a higher dimensional level set function.Conventional level set methods rely on the Hamilton-Jacobi evolution equation to update the design according to a normal velocity field defined on the boundary (e.g., Wang et al, 2003;Allaire et al, 2004).
1 An important aspect of these methods is extending the normal velocity away from the boundary.To this end, the Hilbertian velocity extensionregularisation framework, which is well known in the context of level set methods (see discussion by Allaire et al (2021b)), can be used to generate a velocity field that guarantees a descent direction and has additional regularity (smoothness) over the whole computational domain.
Topology optimisation problems often include multiple constraints.In density-based topology optimisation (Bendsøe, 1989;Rozvany et al, 1992) the application of constraints is usually straightforward and handled by the optimisation algorithm (e.g., Method of Moving Asymptotes (Svanberg, 1987)).In the context of conventional level set-based methods, applying constraints is more complicated.The augmented Lagrangian method (Nocedal and Wright, 2006;Birgin and Martínez, 2009) is a classical approach for constrained optimisation problems.It converts a constrained optimisation problem into a sequence of unconstrained problems that are a combination of the classical Lagrangian method and quadratic penalty method.In a level set framework, applying the method is straightforward: the shape sensitivity of the augmented Lagrangian is used to inform the normal velocity for the Hamilton-Jacobi equation (e.g., Guo et al, 2014;Allaire et al, 2016;Cao et al, 2021).However, the difficulty associated with tuning the accompanying parameters is problem dependent and scales with the number of constraints (Allaire et al, 2021b).The levelset sequential linear programming method (SLP) (Dunning et al, 2015;Dunning and Kim, 2015) involves linearising the optimisation problem into a number of sub-problems that are then solved using a linear programming method (e.g., the simplex method (Kambampati et al, 2020)).For level set-based topology optimisation, applying SLP is fairly straightforward except for implementing appropriate trust region constraints (Dunning and Kim, 2015).The projection method (Wang and Wang, 2004;Challis et al, 2008) projects the objective shape sensitivity onto a space that will leave the constraints unchanged and combines this with constraint shape sensitivities.This approach has been used to successfully design material microstructures subject to isotropy constraints (Challis et al, 2008) but has not been widely adopted in the literature.However, similar methods have more recently been proposed in the literature by Barbarosie et al (2020) and Feppon et al (2020) for level set topology optimisation.The projection method and these two recent works are examples of general null-space gradient methods (e.g., Nocedal and Wright, 2006).
It is natural to consider methods of constrained optimisation that take advantage of the Hilbertian extension-regularisation framework.For example, Allaire et al (2021b) recently presented an SLP method in the Hilbertian framework.In this paper we revisit the projection method from Challis et al (2008) and combine it with the Hilbertian extension-regularisation procedure.Our method constructs an orthogonal basis that spans the set of extended constraint shape sensitivities and an orthogonal projection operator that projects onto the set perpendicular to the extended constraint shape sensitivities.We then define the normal velocity for the level set function as a linear combination of the orthogonal projection operator applied to the extended objective function shape sensitivity and a weighted sum of basis functions for the extended constraint shape sensitivities.This normal velocity is naturally extended onto the bounding domain and endowed with additional regularity due to the Hilbertian extensionregularisation.While our method is similar to other recently proposed approaches (Barbarosie et al, 2020;Feppon et al, 2020), our formulation utilising an orthogonal basis provides significant benefits.
To demonstrate our presented Hilbertian projection method we consider several linear elastic microstructure optimisation (i.e., inverse homogenisation) problems with multiple constraints.The constraints naturally arise under the enforcement of symmetries for the effective material properties, such as isotropy.Irrespective of optimisation method, microstructure optimisation has been used successfully for a range of design problems including linear elastic materials with extremal properties (e.g., Gibiansky and Sigmund, 2000;Andreassen et al, 2014), multifunctional composites (e.g., Challis et al, 2008), auxetic materials (e.g., Vogiatzis et al, 2017), piezoelectric materials (e.g., Silva et al, 1998;Wegert et al, 2022), and multi-material composites (e.g., Zhuang et al, 2010;Faure et al, 2017).In this work we consider maximising the bulk modulus with and without isotropy constraints and the design of auxetic and multi-phase materials.We compare our Hilbertian projection method with a Hilbertian sequential linear programming (SLP) method (Sec. 5.3.2, Allaire et al, 2021b) and show that the Hilbertian projection method is able to successfully handle these optimisation problems with little-to-no parameter tuning.
The remainder of the paper is as follows.In Section 2 we discuss the mathematical background for the level set method, Hilbertian extension-regularisation procedure, and linear elastic microstructure optimisation.In Section 3 we formulate the Hilbertian projection method and compare our formulation to the null space method presented by Feppon et al (2020).In Section 4 we discuss our numerical implementation.In Section 5 we present and discuss the example optimisation problems and results.Finally, in Section 6 we present our concluding remarks.

Mathematical background
In this section we give a brief introduction to the level set method for topology optimisation and shape derivatives.We then discuss the Hilbertian extension-regularisation framework.We conclude by describing linear elastic microstructure optimisation for single-and multi-phase materials.

The level set method
Level set methods track the boundary of a domain Ω inside a bounding domain D ⊂ R d implicitly via the zero level set of a function ϕ : D → R (Sethian, 1996;Osher and Fedkiw, 2006).For a domain Ω inside a bounding domain D, the level set function ϕ is typically defined as (1) Using this definition and assuming that the interface may evolve in time, a material derivative of ϕ on ∂Ω gives where v is the normal velocity of the interface.In practice, the above is solved over the whole bounding domain D instead of only on the interface ∂Ω by extending the velocity v away from the boundary.Assuming that the time interval (0, T ) is small so that the velocity does not vary in time gives the Hamilton-Jacobi evolution equation (Sethian, 1996;Osher and Fedkiw, 2006;Allaire et al, 2021b): where ϕ 0 (x) is the initial condition for ϕ at t = 0.It is often useful to reinitialise the level set function as the signed distance function d Ω .This ensures that the level set function is neither too steep nor too flat near the boundary of Ω (Osher and Fedkiw, 2006).The signed distance function may be defined as (Allaire et al, 2021b): where d(x, ∂Ω) := min p∈∂Ω |x−p| is the minimum Euclidean distance from x to the boundary ∂Ω.
Several methods are available for constructing the signed distance function and the reader is referred to Osher and Fedkiw (2006) and Allaire et al (2021b) and the references therein for a detailed discussion.In this work we use the following reinitialisation equation (Peng et al, 1999;Osher and Fedkiw, 2006) to reinitialise a pre-existing level set function ϕ 0 (x) as the signed distance function: Here S is the sign function and Equation 5 is solved until close to steady state.Similar numerical schemes may be used to solve for Hamilton-Jacobi evolution and reinitialisation (Eqs.3 and 5).

Shape derivatives
To find a normal velocity v that reduces some functional J(Ω) via solution of the Hamilton-Jacobi equation (Eq. 3) we use the notion of shape derivatives.We recall the following useful results from Allaire et al (2004Allaire et al ( , 2021b)).Suppose that we consider smooth variations of the domain Ω of the form Ω θ = (I + θ)(Ω), where θ ∈ W 1,∞ (R d , R d ).Then the following definition and lemma follow: Definition 1 (Allaire et al ( 2004)) The shape derivative of J(Ω) at Ω is defined as the Fréchet derivative with lim θ→0 |o(θ)| ∥θ∥ = 0, where the shape derivative Lemma 1 (Allaire et al (2004)) Let Ω be a smooth bounded open set and f ∈ W 1,1 (R d ).Define Then J is differentiable at Ω and Céa's formal method (Céa, 1986) can be applied to find the shape derivative of a functional J that depends on fields that satisfy specified state equations (e.g., Allaire et al, 2004Allaire et al, , 2021b)).The method relies on defining a Lagrangian functional L that satisfies the two following properties: 1.The state equations are generated by stationarity of L under variations of the fields.2. L is equal to the functional of interest J at the solution to the state equations.
Once these properties are satisfied the shape derivative of the functional of interest can be found using Lemma 1 (Allaire et al, 2004).

Hilbertian extensionregularisation
To infer a descent direction from J ′ (Ω) we utilise the Hilbertian extension-regularisation method as discussed by Allaire et al (2021b).This involves solving an identification problem over a Hilbert space H on D with inner product ⟨ For an unconstrained optimisation problem the resulting field g Ω is the extended shape sensitivity that is used to evolve the interface with θ = τ g Ω n where τ > 0 is sufficiently small.The Hilbertian extension-regularisation method provides two important benefits: it naturally extends the shape sensitivity from ∂Ω onto the bounding domain D, and ensures a descent direction for J(Ω) with additional regularity (i.e., H as opposed to L 2 (∂Ω)) (Allaire et al, 2021b).As discussed by Allaire et al (2021b), this may be viewed as an analog to the sensitivity filtering used in density-based topology optimisation algorithms.
A common choice for the Hilbert space H is H 1 (D) with the inner product where β is the so-called regularisation length scale (e.g., Allaire et al, 2016;Feppon et al, 2019;Allaire et al, 2021a).For microstructure optimisation we use the periodic Sobolev space H = H 1 per (D) and use the inner product defined in Equation 10.

Linear elastic microstructure optimisation
In this section we briefly discuss computational homogenisation and topology optimisation in the context of periodic microstructure design.
The state equations for linear elastic homogenisation over a domain Ω contained in a representative volume element (RVE) D ⊂ R d under an applied strain field εij are (e.g., Yvonnet, 2019) −σ ij,i = 0 in Ω, (11) where σ ij is the stress tensor, ε ij = ε ij (u) is the D-periodic strain field with displacement u, and C ijkl is the spatially dependent elasticity tensor.Note that in the above we use summation notation for indices and comma notation for derivatives.
To compute the homogenised stiffness tensor Cijkl of a periodic material, the above state equations are solved over Ω for three (d = 2) or six (d = 3) different combinations of macroscopic strain fields.These macroscopic strain fields are applied by decomposing the strain into the constant macroscopic strain field and fluctuation strain field as ε ij = εij + εij .The macroscopic strain fields are then given by the unique components of ε(kl in k and l.For example, in two dimensions the unique macroscopic strains are: The notation ε(kl) ij is used to denote the strain field fluctuation arising from the applied strain field ε(kl) ij .In practice, Equations 11-15 are solved using a finite element method and the weak formulation given here: Weak form 1 For each unique constant macroscopic strain field ε(kl

Single-phase problems
For single-phase problems (one solid and a void phase), once the solution ũ(ij) to Weak Form 1 has been found for each unique macroscopic strain ε(ij) pq , the resulting homogenised stiffness tensor may be computed via (Yvonnet, 2019) To evaluate Weak Form 1 and the homogenised stiffness tensor above, we utilise the ersatz material approximation.This method, which is classical in the literature (e.g., Allaire et al, 2004), fills the void phase with a soft material so that the state equations can be resolved without a body-fitted mesh.To this end, for small ε void we take and relax integration to be over D. We can provide a smooth approximation to Equation 18 using a smoothed Heaviside function where η is half the length of the small transition region of H η (ϕ) between 0 and 1. Equation 18can then be replaced with (20) It is important to note that the ersatz material approximation is consistent (Allaire et al, 2021b).That is, as ε void → 0, the approximation becomes exact.
We conclude this section by stating the shape derivative of the homogenised stiffness tensor: Lemma 2 The shape derivative of Equation 17 is given by Proof See Appendix A. □

Multi-phase problems
For multi-phase problems, we utilise the colour level set method in which up to 2 M phases can be represented by the sign of M level set functions (Wang and Wang, 2004;Allaire et al, 2014a).For example, in the case of four phases Ω 1 , Ω 2 , Ω 3 and Ω 4 with two level set functions ϕ 1 and ϕ 2 , we have We further denote the domains associated with each level set function ϕ 1 and ϕ 2 as D 1 and D 2 , respectively.Figure 1 shows an illustration of this case.
In a similar way to Equation 20we can interpolate the value of the stiffness tensor between these domains via where C ijkl,α is the elasticity tensor for the phase occupying Ω α .In this multi-phase case we have replaced the level set functions ϕ 1 and ϕ 2 with d D1 and d D2 denoting their respective signed distance functions.This change facilitates shape differentiation.Unlike the situation discussed by Allaire et al (2016), periodicity of D ensures that this replacement is valid provided the level set functions are reinitialised often.Additional care should then be taken regarding the calculation of certain quantities and their shape derivatives.Namely, the homogenised stiffness tensor (Eq.17) becomes (25) Similar expressions are used for Ω 2 , Ω 3 and Ω 4 .
Integration over the whole cell D in Equations 24 and 25 with dependence on the signed distance functions differs from the single-phase case where integration is over the domain Ω occupied by the solid phase (see Eq. 17).Shape differentiability of the signed distance function and a coarea formula can be used in this multi-phase case to derive the shape derivative.We utilise the "approximate" formula discussed by Allaire et al (2014a).This assumes that the Heaviside smoothing parameter η is small and that the principal curvatures of ∂Ω vanish.
Lemma 3 The approximate shape derivatives of Equations 24 and 25 under variation of the domain where g = Hη(d D1 ), and Analogous expressions follow for θ 2 and Ω 2 , Ω 3 and Ω 4 .
Proof See Appendix B. □ We note that comparisons between the "true" formula, "Jacobian-free" formula (zero principal curvatures), and "approximate" formula have been discussed for compliance elsewhere in the literature (Allaire et al, 2014a(Allaire et al, , 2016)).It suffices to mention that Allaire et al (2016) found that the "approximate" formula does not capture the distortion that arises due to the ray integration and approximation of the principal curvatures in the "true" formula.

Hilbertian projection method
The Hilbertian framework yields a descent direction θ = τ g Ω n for unconstrained optimisation problems.However, for constrained optimisation problems such as min the choice of θ is more difficult.
In the literature a variety of optimisation methods deal with this problem but few of these take advantage of the Hilbertian framework.Allaire et al (2021b) recently presented a sequential linear programming (SLP) method in the Hilbertian framework.The projection method uses orthogonal projections to evolve the design in a direction that aims for best possible improvement of the objective functional while improving the constraint functionals (Challis et al, 2008).In the following we present a Hilbertian extension of the projection method for constrained topology optimisation.

Preliminaries
We proceed by first solving the following set of scalar Hilbertian extension-regularisation problems over H for an objective functional J(Ω) and constraint functionals C i (Ω): Next we use Gram-Schmidt orthogonalisation to remove linearly dependent constraints from the set {µ Ωi } N i=1 to obtain the set {µ Ωp } N p=1 , where N ≤ N .We use {μ Ωp } to denote the corresponding orthogonal basis that spans the set C ⊂ H of extended constraint shape sensitivities.The basis {μ Ωp } can be used to construct an orthogonal projection operator P C ⊥ that projects the shape sensitivity g Ω onto the set C ⊥ perpendicular to the set of extended constraint shape sensitivities.We define this operator as Then, by construction, evolving the level set function using normal velocity P C ⊥ g Ω would to first order improve the objective functional J(Ω) while leaving the constraint functionals C i (Ω) unchanged.On the other hand, the set of basis functions {μ Ωp } describes directions that to first order improve the constraint functionals C i (Ω).

Formulation
The Hilbertian projection method can then be formulated as follows: For some rate parameter λ ∈ R, suppose we choose v Ω ∈ H in the deformation field θ = τ v Ω n so that J(Ω) and C i (Ω) decrease via (Challis et al, 2008) It is important to note that we purposefully pose the former requirement as "min possible" so that the objective functional may increase when required to improve constraints.Furthermore, linearly dependent constraints in the optimisation problem need to be consistent to ensure that the second line of Equation 32 is well posed.Specifying constraints that have linearly dependent shape sensitivities but for which the directions of improvement are in contradiction would violate this requirement.
In the Hilbertian framework, Equation 32 may be rewritten as Note that the change in sign comes from the application of Equations 29 and 30.We choose the following linear combination for v Ω : where α p ∈ R are determined using ⟨µ Ωp , v Ω ⟩ H = λC p for p = 1, . . ., N .This generates a lowertriangular linear system of the form which can easily be solved via forward substitution (Challis et al, 2008).
To first order the first term of Equation 34 improves the objective while leaving the constraints unchanged due to orthogonality, while the second term improves the constraints with extended shape sensitivities that contribute to the basis {μ Ωp }.In the numerical examples we have observed satisfaction of all constraints at convergence of the optimisation algorithm, including those that have linearly dependent shape sensitivities.The square root term of Equation 34 is included to facilitate a balance between improving the objective and constraints.

Parameters
As discussed by Challis et al (2008), the rate parameter λ should be chosen to ensure 1 The new parameter α min then controls the balance between improving the objective or constraints.For example, α min = 1 ignores the objective function in Equation 28 and instead solves a constraint satisfaction problem.As a result, the method only has a single parameter α min while λ is dictated by the inequalities above.In general we find that α min does not require fine tuning and unless otherwise stated we choose α 2 min = 0.1.

Comparison to null space methods
Our formulation is similar to null space methods recently developed by Barbarosie et al (2020) and Feppon et al (2020), both of which present a similar formulation.In the following we discuss some differences between our Hilbertian projection method and the null space method proposed by Feppon et al (2020).Most notably, our formulation makes use of an orthogonal basis for the set of extended constraint shape sensitivities.This avoids a possibly expensive matrix inversion that appears in the algorithm presented by (Feppon et al, 2020).Our use of the orthogonal basis also avoids reliance on the 'Linear Independence Constraint Qualification' (LICQ) condition (Sec.2.1, Feppon et al, 2020).Our method therefore naturally handles linearly dependent constraint shape sensitivities.Such dependencies often appear in microstructure optimisation problems when symmetries are imposed on the effective material properties (e.g., Secs.5.2 and 5.4 below).Multi-phase level setbased topology optimisation via the colour level set method (Wang and Wang, 2004;Allaire et al, 2014a) can also give rise to such linear dependency (e.g., Sec.5.4 below).The ability of the Hilbertian projection method to naturally handle these situations gives the user more freedom in how topology optimisation problems are posed and avoids additional special treatment of linearly dependent constraint sensitivities.
For equality constrained problems when LICQ is satisfied, both the null space method (Feppon et al, 2020) and our Hilbertian project method give equivalent directions for improvement of the objective and violated constraints, up to coefficients α p and constant λ.However, the method of attaining this improvement is quite different.In this work, the second term of Equation 34 is a linear combination of the orthogonal basis of the set of extended constraint shape sensitivities.The coefficients α p are chosen to improve constraints exponentially, as per the second requirement in Equation 33.The null space instead uses a Gauss-Newton direction for ensuring exponential decay of violated constraints (Lemma 2.5 & Prop.2.6, Feppon et al, 2020).This again relies on LICQ and possibly expensive matrix inversion discussed above.
Finally, unlike the null space methods proposed by Barbarosie et al (2020) and Feppon et al (2020), the Hilbertian projection method as formulated above is unable to handle inequality constraints.Such an extension could be considered in future using slack variables or by adopting the procedure used by either Barbarosie et al (2020) or Feppon et al (2020).Interestingly, in terms of the total length covered to reach the optimum, Feppon et al (2020) found that their dual quadratic programming method for handling inequality constraints yielded equivalent performance when compared to the method of slack variables.However, using slack variables introduces additional computational cost and possibly further parameter tuning (Feppon et al, 2020).For our implementation this additional cost should be small owing to the use of orthogonalisation.As such the slack variable method would be an appropriate first recourse for implementing inequality constraints within the Hilbertian projection method.This would be similar to the approach taken by Schropp and Singer (2000).

Numerical implementation
In the following we describe the numerical implementation of our topology optimisation algorithm.We first discuss the resolution of state equations and Hamilton-Jacobi type equations followed by an overview of the optimisation algorithm.We finish with a brief discussion of the Hilbertian SLP method which is compared to our presented Hilbertian projection method.

Resolving state and Hamilton-Jacobi-type equations
To solve the state equations and the Hilbertian extension-regularisation problems we use the finite element package Gridap (Badia and Verdugo, 2020;Verdugo and Badia, 2022) in the programming language Julia.In particular, we discretise the periodic domain D ⊂ R d into n d linear quadrilateral (d=2) or hexahedral (d=3) elements with element width ∆x and discretise the level set function at the nodes of the triangulation.
To reduce computational cost when solving the state equations, we remove any elements that are completely void phase and leave a strip of ersatz material near the phase interface.The resulting linear systems for the state equations and Hilbertian extension-regularisation problems are then solved using a direct method in 2D or a GPUbased Jacobi pre-conditioned conjugate gradient method in 3D.
For the Hamilton-Jacobi evolution equation and signed distance reinitialisation equation (Eqs.3 and 5) we use standard first order Godunov upwind finite difference schemes (Peng et al, 1999;Allaire et al, 2004;Osher and Fedkiw, 2006) that have been implemented on GPUs using CUDA.jl(Besard et al, 2018).For the Hamilton-Jacobi evolution equation we use ⌊n/10⌋ or ⌊n/3⌋ number of time steps in two dimensions or three dimensions, respectively.We are more conservative in three dimensions because we have less elements along each axial direction.For the reinitialisation equation we iterate until reaching a stopping condition where q is the iteration number.In addition, for the sign function S we use the common approximation that applies on a Cartesian grid with square elements of side length ∆x (Osher and Fedkiw, 2006).

Algorithm overview
In Algorithm 1 we present our optimisation algorithm that is based on the theory discussed in Sections 2 and 3. Algorithm 1 is similar to Algorithm 5 presented by Allaire et al (2021b), with the addition of a line search method for determining the Courant-Friedrichs-Lewy (CFL) coefficient γ for solving the Hamilton-Jacobi evolution equation (e.g., Allaire et al, 2021a) with time step (Osher and Fedkiw, 2006) Note that we omit the indices on γ that appear in Algorithm 1 for sake of clarity.In general, the line search method helps to remove oscillations in the optimisation history and improve convergence.For the stopping criterion we require that the current objective value compared to the previous five is stationary and that the constraints are satisfied within specified tolerances.The Hilbertian projection method is implemented using the package DoubleFloats.jlthat gives a machine epsilon of roughly 5 × 10 −32 .This prevents accumulation of round-off error when generating the projection operator that can affect the optimisation history.All other computations are completed in standard double precision.
Table 1 gives the parameter values used for all the optimisation examples unless otherwise stated in Section 5.

Algorithm 1 Optimisation Algorithm
Initialisation: Initialise Ω 0 inside a computational domain D with mesh T and a level set function ϕ 0 .1: Find the initial solution u (ij) to the homogenisation problem for each unique ε(ij) .2: for q = 1, . . ., qmax do 3: Calculate the shape sensitivity of the objective J(Ω q−1 ) and constraints C i (Ω q−1 ).

4:
Solve scalar Hilbertian extension-regularisation problems for the objective and constraints with length scale β.

5:
Apply the Hilbertian projection method to find v Ω .

9:
Solve state equations for linear elastic homogenisation and calculate the new objective Jnew.

16:
Reject the new iteration and continue loop.

18:
end for return Ω q and ϕ q . 21: end if 22: end for

Comparison to sequential linear programming (SLP)
We compare the Hilbertian projection method to the Hilbertian SLP method presented by Allaire et al (2021b) (Sec.5.3.2).We make two adjustments to the method.Firstly, to match our formulation we replace the inequality constraints with equality constraints.In addition, we change the trust region constraints for the constraint functionals to be for i = 1, . . ., N .For several two-dimensional problems we find that this choice promotes convergence and better optimisation results.However, as we will discuss later, choosing trust region constraints is not straightforward for our example optimisation problems.We implement Hilbertian SLP by adjusting line 5 of Algorithm 1 accordingly.To solve the resulting linearised optimisation problem we use the Julia packages JuMP.jl (Dunning et al, 2017) and Ipopt.jl(Wächter and Biegler, 2006).

Example problems
In the following we give the optimisation results for several example problems that have been solved with both the Hilbertian projection method and Hilbertian SLP method.

Example 1: Maximum bulk modulus
In this example we consider a bounding domain D = [0, 1] d that contains a solid phase and void phase.The solid phase is constructed from an isotropic medium with E = 1 and ν = 0.3.Subject to a volume constraint Vol(Ω) = 1/2, we maximise the effective bulk modulus κ(Ω) of the material.The bulk modulus is a measure of stiffness to volumetric strain given by in three dimensions.In other words, we seek to solve the optimisation problem: The last line represents the satisfaction of the state equations.
In two dimensions we use a periodic starting structure with four equally spaced holes.For three dimensions the initial boundary between void and solid material is given by a Schwarz P minimal surface.It is well-known that in two dimensions hole nucleation is not possible under Hamilton-Jacobi evolution (e.g., Allaire et al, 2004).For this reason we initialise the two-dimensional optimisation problems with more holes than required.Topological derivatives could be incorporated to rectify this, but this is outside of the scope of the current paper.We also note that different starting structures could be used provided they are periodic and have non-zero stiffness.
Figures 2 and 3 show the starting structures and optimisation results for two and three dimensions, respectively for both the Hilbertian projection method and SLP.In addition, we compare the objective value with the Hashin-Shtrikman (HS) upper bound (Hashin and Shtrikman, 1963).Table 2 shows a summary of the results.In two dimensions both the Hilbertian projection method (Fig. 2b/2d) and Hilbertian SLP (Fig. 2c/2e) perform well, converging in 32 and 12 iterations respectively at 99.71% and 99.51% of the HS bound.In addition, the resulting structures are geometrically similar and match classical results in the literature (Sec.2.10.3,Bendsøe and Sigmund, 2004).It is worth noting that no topological changes occur in this example.
In three dimensions the Hilbertian projection method converges to 99.86% of the HS bound while Hilbertian SLP converges to 71.07% of the bound.Over the course of the iteration history, topological changes occur for the Hilbertian projection method, while SLP fails to evolve from the initial topology.The resulting structure from the  Hilbertian projection method matches other classical results (Sec.2.10.3,Bendsøe and Sigmund, 2004).

Example 2: Maximum bulk modulus with isotropy
In Example 2 we consider the same problem setup as Example 1 with the addition of macroscopic isotropy constraints.This ensures that the resulting homogenised stiffness tensor is invariant under rotation.
In two dimensions the optimisation problem is given by min The constraints C i (Ω) are given by where μ is the isotropic shear modulus given by Analogous expressions appear in three dimensions and the number of isotropy constraints increases from 6 to 21 (Challis, 2009).It should be noted that the term 4κ 2 + 8μ 2 appears as a normalisation coefficient and is considered constant for the purpose of the shape differentiation.Furthermore, owing to the symmetry of the isotropy constraints the extended constraint shape sensitivities have a nullity of two.Our method handles these with no special treatment.
To visualise the behaviour of the isotropy constraints in the iteration history, we define the effective anisotropy Ā to be the sum of squares of the violation of these constraints.That is, For the two-dimensional problem we implement both the full set of isotropy constraints as well as the single constraint Ā = 0. We find that the SLP method struggles with the full set of constraints.For this reason we instead use Ā = 0 as the isotropy constraint for the SLP method.For the Hilbertian projection method we find that the full set of constraints is more effective.This matches previous literature (e.g., Challis et al, 2008;Challis, 2009).
Figure 4 shows the two-dimensional results for the Hilbertian projection method with the full set of constraints and single anisotropy measure constraint, and SLP with the anisotropy measure constraint.We omit the starting structure as this is the same as previously (Fig. 2a). Figure 5 shows the three-dimensional results for the Hilbertian projection method and SLP.Table 3 shows a summary of the results for this example.
In two dimensions the Hilbertian projection method performs well for the full set of constraints and the single anisotropy measure constraint.The number of iterations is markedly lower when using the full set of constraints (78 vs 220).This is unsurprising because including the individual symmetry constraints enables the projection method to improve each constraint separately (Challis et al, 2008;Challis, 2009).The final optimisation values with the full set or single constraint are very similar, being 99.68% and 99.69% of the HS bound, while the final structures are geometrically identical under a periodic shift of half the cell edge length along each coordinate direction.The resulting structures match other classical results (Sec.2.10.3,Bendsøe and Sigmund, 2004).In contrast, the SLP method does not manage to reduce the measure of anisotropy to zero and reaches the maximum number of iterations.However, the final structure obtained with SLP (Fig. 4c) is very close to those obtained with the Hilbertian projection method (Fig. 4a and 4b).We suspect that SLP fails to converge due to the difficulty associated with choosing the trust region constraints.
For the three-dimensional results, we find that the Hilbertian projection method with the full    Bendsøe and Sigmund, 2004).For SLP, the optimisation algorithm fails to converge in 1000 iterations and the final structure is clearly not optimal.
These results demonstrate that the Hilbertian projection method is able to handle a large number of constraints (22 in 3-dimensions) without requiring any parameter tuning of α min or λ.In addition, the optimisation histories for the Hilbertian projection method with the full set of isotropy constraints (Fig. 4d and 5c) are smooth and converge fairly quickly.

Example 3: Auxetic materials
In this example we consider two-dimensional minimum volume auxetic materials with a Poisson's ratio of −0.5.For the problem set up, we consider a bounding domain D = [0, 1] 2 that contains a solid phase and void phase.As previously, the solid phase is constructed from an isotropic material with E = 1 and ν = 0.3.
It may be noted that this is similar to the approach taken by Vogiatzis et al (2017).However, they instead minimise a weighted sum of Cijkl and prescribed value subject to a volume constraint.
For the purpose of this example we choose C1111 = 0.1, which results in C1122 = −0.05.The resulting optimisation problem is then min For this example we use γ max = 0.05 and α 2 min = 0.5.This means that the optimiser favors improvement of the constraints rather than the objective and does not move too quickly to avoid disconnecting in the first few iterations.
Figure 6 shows the starting structure (Fig. 6a) and optimisation results (Fig. 6b and 6c) for this problem using the Hilbertian projection method.We use a periodic starting structure with sixteen equally spaced holes.The method converges in 61 iterations with a final volume of 0.3159 and Poisson's ratio of -0.4998.In contrast, SLP fails because the optimiser prioritises the objective leading to a disconnected solid phase and the algorithm is unable to recover.

Example 4: Multi-phase materials
For our final example we consider two-dimensional multi-phase maximum bulk modulus problems with and without isotropy constraints.We consider a bounding domain D = [0, 1] 2 that contains two solid phases and void.The solid phase contained in Ω 2 is an isotropic material with E = 1 and ν = 0.3, while Ω 3 contains isotropic material with E = 1/2 and ν = 0.3.The void phase is contained in Ω 1 and Ω 4 .As previously, most void phase is removed from the mesh (see Sec. 4) while any material close to the interface is specified as weak material with E = 10 −3 and ν = 0.3.We consider two optimisation problems.The first is maximum bulk modulus subject to volume constraints on Ω 2 and Ω 3 given by min The second is maximum bulk modulus subject to macroscopic isotropy constraints and volume constraints on Ω 2 and Ω 3 given by min where the constraints C i are as in Example 2. For SLP we use the single anisotropy constraint and for the Hilbertian projection method we use the full set of isotropy constraints.For these examples we use γ max = 0.05 so that the optimiser does not evolve the designs too quickly.It should be noted that for the case of only volume constraints the extended constraint shape sensitivities have a nullity of one (in θ 1 ) and zero (in θ 2 ) owing to the structure of the shape derivatives for the volume constraints.For the case of volume and isotropy constraints the extended constraint shape sensitivities have a nullity of three (in θ 1 ) and two (in θ 2 ) due to the underlying symmetry of the shape derivatives for the volume and isotropy constraints.Our method handles these with no special treatment.
We initialise with two overlapping level set functions that give a starting structure completely comprised of the less stiff material (Ω 3 ) and the void phase.Regions of stiffer material (Ω 2 ) are readily generated during the optimisation via independent evolution of the two level set functions.We use a starting structure of four equally spaced holes for the optimisation problem without isotropy constraints and nine equally spaced holes for the problem with isotropy constraints.
Figures 7 and 8 show the optimisation results and history without and with isotropy constraints, respectively.We denote the stiff and less stiff material phase by blue and green, respectively, while the smooth interface is given by the dark green overlap.We include the Hashin-Shtrikman-Walpole (HSW) upper bound (Walpole, 1966) for this problem as the black dashed line in these figures.Table 4 shows a summary of the results.
For the case of no isotropy constraints, both the Hilbertian projection method and SLP converge to roughly 96% of the HSW bound while the resulting structures (Fig. 7b and 7c) are geometrically similar apart from the thin interface that presents in the results for the Hilbertian projection method.The difference between our results and the theoretical upper bound is likely due to the use of the approximate formula for    et al (2014a) showed that the approximate formula yields slightly less optimal results than the true or Jacobian-free counterpart.The iteration history for the Hilbertian projection method (Fig. 7d) is fairly smooth while the history for SLP (Fig. 7e) moves rapidly at the beginning of the optimisation to satisfy the volume constraints and increase the objective.As previously noted, this is likely due to the trust region constraints that, in this case, need to be chosen to be more conservative.
With the addition of isotropy constraints, we again find that the Hilbertian projection method works well without significant parameter tuning.In 148 iterations the optimisation algorithm is able to satisfy all constraints and obtain a final objective value that is 93.44% of the HS bound.In contrast, SLP does not converge within 1000 iterations.

Conclusion
In this paper we have presented a Hilbertian extension of the projection method for constrained level set-based topology optimisation.At its core, the method relies on the Hilbertian extensionregularisation method in which a set of identification problems are solved over a Hilbert space H on D with inner product ⟨•, •⟩ H .This procedure naturally extends shape sensitivities onto the bounding domain D and enriches them with the regularity of H.For a constrained optimisation problem the projection method framework aims for the best first-order improvement of the objective in addition to first-order improvement of the constraints.These requirements for the projection method may then be reposed in the Hilbertian framework in terms of H and ⟨•, •⟩ H .We satisfy these reposed requirements by defining the normal velocity of the level set function as a linear combination of the orthogonal projection operator applied to the extended objective shape sensitivity and basis functions for the extended constraint shape sensitivities.Owing to the Hilbertian extension-regularisation of shape sensitivities, the chosen normal velocity is already extended onto the bounding domain D and endowed with the regularity of H.
To demonstrate the Hilbertian projection method for constrained level set-based topology optimisation we have solved several example microstructure optimisation problems with multiple constraints.We showed that the Hilbertian projection method successfully handled all of these optimisation problems with little-to-no tuning of the parameter α min that controls the balance between improving the objective and constraints.The Hilbertian projection method also naturally handles linearly dependent constraint shape sensitivities.Such linear dependencies often appear in microstructure optimisation and multi-phase optimisation problems.
We found that our method performs well when compared to a Hilbertian sequential linear programming (SLP) method (Allaire et al, 2021b).For problems only involving volume constraint(s), both methods converged to appropriate optimised microstructures.However, SLP did not successfully solve some of the more complex example optimisation problems.These results demonstrate the capacity of the Hilbertian projection method and likely other projection/null space methods (e.g., Barbarosie et al, 2020;Feppon et al, 2020) for solving constrained level set-based topology optimisation problems.
Applying multiple constraints is challenging in level set-based topology optimisation owing to the reliance on implicitly defined domains.Alongside other recent work (Barbarosie et al, 2020;Feppon et al, 2020;Allaire et al, 2021b), our proposed Hilbertian projection method makes significant progress towards improving the capacity of conventional level set-based methods for constrained topology optimisation.Furthermore, due to its generality, the Hilbertian projection method is not confined to microstructure optimisation.It may be applied to any topology optimisation problem to be solved in a level-set framework, including macroscopic or multi-physics problems.Inequality constraints could likely be incorporated into the method using slack variables.In addition, the method is not confined to Eulerian level set methods and can be readily applied to Lagrangian or body-fitted level set methods (e.g., Allaire et al, 2014b).These extensions could be considered in future work.∂C pqrs ∂g × (ε pq (ũ (ij) ) + ε(ij) pq )(ε rs (ũ (kl) ) + ε(kl) rs ) dz dΓ (B7) where g(x) = H η (x).
Finally, the support of H ′ η (x) is |x| < 2η, so the integral over ray ∂D1 (x) ∩ D is restricted to a tubular region about ∂D 1 (Allaire et al, 2014a).Therefore, for small η, we may assume that (B14) which concludes the proof.

Fig. 1 :
Fig. 1: An illustration of colour level sets with two level set functions and four phases.(a) and (b) show the domain D 1 and D 2 represented via the level set function ϕ 1 and ϕ 2 respectively.(c) shows the colour representation of Ω 1 , Ω 2 , Ω 3 , and Ω 4 for different signs of ϕ 1 and ϕ 2 .

Fig. 2 :
Fig. 2: Two-dimensional optimisation results for Example 1: maximum bulk modulus.For the starting structure (a), (b) and (c) show the final structures for the Hilbertian projection method and SLP respectively, while (d) and (e) show the respective iteration histories.The Hashin-Shtrikman upper bound for the bulk modulus is given by the dashed black line.

Fig. 3 :
Fig. 3: Three-dimensional optimisation results for Example 1: maximum bulk modulus.For the starting structure (a), (b) and (c) show the final structures for the Hilbertian projection method and SLP respectively, while (d) and (e) show the respective iteration histories.The Hashin-Shtrikman upper bound for the bulk modulus is given by the dashed black line.

Fig. 4 :
Fig. 4: Two-dimensional optimisation results for Example 2: maximum bulk modulus with isotropy.(a), (b), and (c) show the final structures for the Hilbertian projection method with the full set of constraints, single anisotropy constraint, and SLP with the single anisotropy constraint respectively, while (d), (e), and (f) show the respective iteration histories.The Hashin-Shtrikman upper bound for the bulk modulus is given by the dashed black line.Due to very little change between iteration 600 and 1000, the upper bound on the x-axis in (f) has been reduced to 700.

Fig. 5 :
Fig. 5: Three-dimensional optimisation results for Example 2: maximum bulk modulus with isotropy.(a) and (b) show the final structures for the Hilbertian projection method and SLP respectively, while (c) and (d) show the respective iteration histories.The Hashin-Shtrikman upper bound for the bulk modulus is given by the dashed black line.

Fig. 6 :
Fig. 6: Optimisation results for Example 3: auxetic materials.For the starting structure (a), (b) shows the final structure while (c) shows the iteration history for the Hilbertian projection method.The desired Poisson's ratio of -0.5 is given by the dashed black line.

Fig. 7 :
Fig. 7: Optimisation results for Example 4: maximum bulk modulus multi-phase materials without isotropy constraints.For the starting structure (a), (b) and (c) show the final structures for the Hilbertian projection method and SLP respectively, while (d) and (e) show the respective iteration histories.The Hashin-Shtrikman-Walpole upper bound for the bulk modulus is given by the dashed black line.Note that the volume constraint for each phase is 0.25.

Fig. 8 :
Fig. 8: Optimisation results for Example 4: maximum bulk modulus multi-phase materials with isotropy constraints.For the starting structure (a), (b) and (c) show the final structures for the Hilbertian projection method and SLP respectively, while (d) and (e) show the respective iteration histories.The Hashin-Shtrikman-Walpole upper bound for the bulk modulus is given by the dashed black line.Note that the volume constraint for each phase is 0.25.

Table 1 :
Parameter values used for optimisation examples.Note: Any variation in the parameters from these values will be specifically stated in the relevant example problem section.

Table 2 :
Summary of optimisation results for Example 1: maximum bulk modulus.

Table 3 :
Summary of optimisation results for Example 2: maximum bulk modulus with isotropy.The asterisk denotes failure to converge.

Table 4 :
Summary of optimisation results for Example 4: two-dimensional maximum bulk modulus multi-phase materials.The asterisk denotes failure to converge.