Florian Feppon
email: florian.feppon@kuleuven.be

Density based Topology Optimization with the Null Space Optimizer: a tutorial and a comparison

Keywords: Nonlinear constrained optimization, density based Topology Optimization, null space gradient flows, Python

The Null Space Optimizer is a constrained optimization solver that has been developed in the context of level-set based Topology Optimization. One of its appealing aspects comes from its relative independence to the need for tuning unintuitive algorithm parameters. The first contribution of this paper is to introduce an upgrade of the Null Space Optimizer that enables to solve optimization problems featuring a large number of constraints with sparse Jacobian matrix. This allows to include in particular bound constraints, making it possible to use the Null Space Optimizer for solving density based Topology Optimization problems. The second contribution of the paper is to present three tutorials giving an educational view on how to use the open source Python implementation of the Null Space Optimizer for solving Topology Optimization problems in structural mechanics and conductive heat transfer, on structured and unstructured meshes. Elegant Python programming features are used for automating the implementation of density filters and the assembly of sparse Jacobian matrices. Numerical results are presented on three design problems and compared to those obtained with the Method of Moving Asymptotes (MMA), the Interior Point Optimizer IPOPT and the Optimality Criteria (OC) method. We found that on the situations considered, (i) the Null Space optimizer is able to compute optimized designs with performances comparable to its competitors with very little parameter tuning, (ii) the OC method or IPOPT with default parameters sometimes converge to nonoptimal designs, (iii) MMA sometimes converge to slightly better design with a faster decay than the Null Space Optimizer during the first iterations, but may also require case dependent fixes to converge to satisfactory solutions.

Introduction

The Null Space optimization algorithm (Feppon et al (2020a)) was introduced in the context of level-set based Topology Optimization (TO) [START_REF] Allaire | Casting constraints in structural optimization via a level-set method[END_REF]) as a convenient alternative to the Augmented Lagrangian (AL) [START_REF] Nocedal | Numerical optimization, 2nd edn[END_REF]; [START_REF] Wang | The augmented Lagrangian method in structural shape and topology optimization with RBF based level set method[END_REF]; [START_REF] Xia | Topology optimization of thermoelastic structures using level set method[END_REF]; [START_REF] Allaire | Casting constraints in structural optimization via a level-set method[END_REF]; [START_REF] Deng | Multi-constrained 3D topology optimization via augmented topological level-set[END_REF] ;[START_REF] Emmendoerfer | Stress-constrained level set topology optimization for compliant mechanisms[END_REF]), and to the Sequential Linear Programming (SLP) methods (Dunning and Kim (2015); Dunning et al (2015); [START_REF] Liu | Topological design of microstructures using periodic material-field series-expansion and gradient-free optimization algorithm[END_REF][START_REF] Liu | Multi-material topology optimization considering interface behavior via XFEM and level set method[END_REF]). One of the motivations for developing the Null Space algorithm stemmed from the need for efficient constrained optimization solvers that work with little to no tuning of case dependent algorithm parameters. Such tuning is indeed associated with a number of inconveniences: (i) multiple runs of the optimization are needed in order to find the right parameters leading to a satisfactory optimized design, which is unaffordable for large-scale applications; (ii) the tuning may be quite complex to achieve in situations involving a large number of design constraints; (iii) it limits TO to a set of predefined applications for its use by non-experts in industrial contexts.

The Null Space algorithm is a first-order constrained optimization method that evolves the design variable by solving an Ordinary Differential Equation (ODE). The trajectories of this ODE minimize the objective function while smoothly correcting the violation of the constraints. As a consequence, the sole true algorithm parameter is the time step, which only needs to be taken sufficiently small for ensuring the minimization property along the integrated trajectories. The Null Space algorithm has proved to be quite effective for solving a number of large-scale multiphysics design problems involving multiple physical constraints with level-set based methods (see e.g. [START_REF] Feppon | Shape and topology optimization of multiphysics systems[END_REF]; Feppon et al (2020b[START_REF] Feppon | Body-fitted topology optimization of 2D and 3D fluid-to-fluid heat exchangers[END_REF]; Salazar De Troya et al ( 2021)). It allows to deal with both inequality and equality constraints and corrects unfeasible initializations. Its implementation has been released as a freely available open-source Python package called "Null Space Optimizer" 1 .

When it was introduced in Feppon et al (2020a), the Null Space Optimizer did not allow to solve optimization problems featuring a very large number of constraints: indeed the original implementation required to explicitly compute the inverse of a matrix of size the number of active constraints. In the Appendix C, we detail a revised implementation that takes advantage of sparsity to avoid this explicit inversion. This allows to extend the Null Space algorithm to constrained optimization problems with a large number of constraints provided these have a sparse Jacobian matrix.

As a straightforward by-product of this upgrade, the Null Space Optimizer now handles bound-constraints, which makes it applicable for solving density based Topology Optimization problems. In the density-based TO framework [START_REF] Bendsoe | Topology Optimization: Theory, Methods, and Applications[END_REF]), the design is 1 https://gitlab.com/florian.feppon/null-space-optimizer/ parameterized with density variables (ρ i ) 1≤i≤n which need to satisfy 0 ≤ ρ i ≤ 1 for every 1 ≤ i ≤ n, where n is the number of elements discretizing the computational domain. In contrast with the level-set framework where these bound constraints do not exist, this entails the need to deal with twice as many additional inequality constraints as the number n of design variables.

With this extension at hand, this paper provides a tutorial view on how to use the Python implementation of the Null Space Optimizer for conveniently solving Topology Optimization problems in the density framework. The Null Space Optimizer offers an elegant paradigm for implementing Topology Optimization problems. It automates the use of filters for the design variable, the memoization of the physical solver and the assembly of the Jacobian matrices of the bound constraints. Moreover, it is interfaced with Python implementations of the popular Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF]; [START_REF] Deetman | GCMMA-MMA-Python[END_REF]), of the Optimality Criteria (OC) method [START_REF] Bendsoe | Topology Optimization: Theory, Methods, and Applications[END_REF]), as well as with the primal-dual interior point method optimizer IPOPT [START_REF] Curtis | A note on the implementation of an interior-point algorithm for nonlinear optimization with inexact step computations[END_REF]), which allows us to conveniently provide comparisons with these algorithms.

So far, MMA has been considered as the "gold standard" algorithm for solving density based TO problems with arbitrary inequality constraints [START_REF] Gersborg-Hansen | Topology optimization of channel flow problems[END_REF]; [START_REF] Wang | On projection methods, convergence and robust formulations in topology optimization[END_REF]; [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]; [START_REF] Jensen | Topology optimization for nano-photonics[END_REF]; [START_REF] Liang | Nonlinear topology optimization of parallel-grasping microgripper[END_REF]; [START_REF] Zheng | Topology Optimization of Self-supporting Porous Structures Based on Triply Periodic Minimal Surfaces[END_REF]; [START_REF] Liang | Macro-microscale topological design for compliant mechanisms with special mechanical properties[END_REF]). Alternatively, IPOPT has also been used in a number of applications [START_REF] Swartz | Topology optimization of 3D photonic crystals with complete bandgaps[END_REF]; [START_REF] Alonso | Topology optimization applied to the design of 2D swirl flow devices[END_REF]; [START_REF] Alonso | Topology optimization for blood flow considering a hemolysis model[END_REF]; [START_REF] Sá | Topology optimization of turbulent rotating flows using Spalart-Allmaras model[END_REF]). Although these methods have proved effective in many situations, both MMA and IPOPT remain sensitive to algorithm parameter tuning, requiring sometimes case-dependent fixes to obtain satisfactory optimized designs. For instance, MMA requires the objective and the constraint functions to be rescaled in order to map their values to a bounded interval [START_REF] Svanberg | MMA and GCMMA -two methods for nonlinear optimization[END_REF]), which can be challenging when these functions encounter large variations. In some situations, it may additionally require structural changes of the moving asymptote rules to lead to satisfactory convergence. IPOPT, on the other hand, may require expert tuning of algorithm parameters (such as the strategy for updating the barrier parameter) to obtain good performances [START_REF] Rojas-Labanda | Benchmarking optimization solvers for structural topology optimization[END_REF]). The OC method, is often favored in situations featuring a single equality constraint [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF]; [START_REF] Aage | Giga-voxel computational morphogenesis for structural design[END_REF]), although some generalizations to multiple constraints have been proposed [START_REF] Yin | Optimality criteria method for topology optimization under multiple constraints[END_REF]; [START_REF] Kim | Generalized optimality criteria method for topology optimization[END_REF]).

In the present paper, we consider four benchmark cases in which, the Null Space Optimizer appears to be competitive with standard methods with little to no parameter tuning. Depending on the situations, we found that OC and IPOPT converge to a less performant design, while MMA may require some case dependent tuning to achieve performances similar or slightly better to that of the Null Space algorithm. Our numerical results suggest that (i) the Null Space Optimizer can be a suitable alternative to MMA and IPOPT for density based Topology Optimization, and (ii) it allows to compute designs that have comparable performance to its competitors without the need for any particular structural change to the algorithm, rescaling of the constraints or the tuning of nonphysical parameters.

Please note that this paper does not aim to compare more extensively various optimization solvers: we do not claim the superiority of the Null Space Optimizer in every scenario since we have considered only a limited range test cases and optimization solvers, and we have not invested significant effort in fine-tuning the parameters. Due to these limitations, we have chosen not to provide runtime comparisons, since they can be heavily influenced by parameter tuning and other factors, such as the use of precompiled routines. Actually, our numerical results of Sections 4 and 5 suggest a superiority of MMA to find more efficient initial design updates, but this advantage seems to be lost after 70 to 80 iterations. The reader interested in further comparisons between various optimizers on benchmark test cases in density based TO may refer to [START_REF] Fanni | A Comparison Between Different Topology Optimization Methods[END_REF]; [START_REF] Rojas-Labanda | Benchmarking optimization solvers for structural topology optimization[END_REF].

An overview on the principles of the Null Space optimization algorithm is given in Section 2. The necessary technical updates to account for optimization constraints with sparse Jacobian, including bound constraints are included in the Appendix C. This appendix also presents a first benchmark numerical example challenging the enforcement of bound constraints by the Null Space Optimizer and compares the results with MMA and IPOPT.

The next three sections provide tutorials for implementing three 2D Topology Optimization test cases with the Null Space Optimizer as well as comparisons with MMA, IPOPT and the OC method.

Section 3 details how to implement the classical MBB beam compliance minimization problem with a single volume constraint from the popular 88 lines Matlab code [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF]; [START_REF] Aage | A 200 line topology optimization code[END_REF]). The complete code for solving this example with the Null Space Optimizer is provided in Appendix A. Then, we test the different optimizers on a more challenging version of this problem obtained by imposing a very small volume fraction. This example is interesting in that the MMA algorithm needs a major structural fix to obtain a satisfactory convergence, while the IPOPT optimizer struggles to find a performant optimum. We find that the Null Space Optimizer only requires a sufficiently small time step to find a satisfactory solution.

Section 4 considers the design optimization of a bridge featuring multiple loads as more complex test case that precludes the use the OC method. This test case is implemented by adding very minor modifications to the code of Appendix A. On this example, we find that without extensive tuning, MMA converges to a solution that is slightly better than the one computed by the Null Space Optimizer, while the one computed by IPOPT is slightly worse.

Section 5 provides an additional tutorial describing the implementation of the Topology Optimization of a heat sink on unstructured meshes, involving thus a different physical model and a different finite element discretization. The complete source code for this example is reproduced in Appendix B; it is constituted of a main Python file interfaced with two scripts written in the FreeFEM language [START_REF] Hecht | New development in FreeFem++[END_REF]), which is especially convenient for dealing with the finite element method on unstructured meshes. For this particular example that involves a single volume constraint, we find that the Null Space algorithm, MMA and IPOPT converge to designs with similar performances, while the popular Optimality Criteria method converges to a local minimum with worse performance.

Throughout the paper, we relied on the following implementations of MMA, IPOPT and OC for the purpose of our numerical comparisons:

(i) we use the Python implementation of the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The Method of Moving Asymptotes (MMA) with Some Extensions[END_REF]) provided by Arjen Deetman2 that closely matches the original Matlab implementation distributed by K. Svanberg [START_REF] Svanberg | MMA and GCMMA code[END_REF][START_REF] Svanberg | MMA and GCMMA -two methods for nonlinear optimization[END_REF]). In order to improve MMA's performance, we systematically rescale the objective function so that the initial value is equal to 10. Since MMA does not handle by default equality constraints, a volume equality constraint is treated as an inequality constraint when using this optimizer in Sections 3 and 4. (ii) We rely on the Python interface cyipopt3 in order to use the primal-dual Interior Point Optimizer IPOPT [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF]; [START_REF] Curtis | A note on the implementation of an interior-point algorithm for nonlinear optimization with inexact step computations[END_REF]). Since the optimizer involves many parameters that are not straightforward to tune, we systematically run the optimizer with the default options. (iii) The Optimality Criteria (OC) method [START_REF] Bendsoe | Topology Optimization: Theory, Methods, and Applications[END_REF]) is implemented following the source code of [START_REF] Aage | A 200 line topology optimization code[END_REF].

Presentation of the Null Space Optimizer

The Null Space algorithm is a numerical method for solving generic nonlinear constrained optimization problems of the form

min x∈X J(x) s.t. g i (x) = 0 for 1 ≤ i ≤ p, h j (x) ≤ 0 for 1 ≤ j ≤ q.
(1)

Here, X denotes the design set, x is the design variable and the functions J, g := (g i ) 1≤i≤p and

Null Space gradient flows

The basis of the Null Space algorithm for solving the optimization problem ( 1) is to numerically integrate the following Ordinary Differential Equation (ODE), or so-called "null space gradient flow":

ẋ(t) = -α J ξ J (x(t)) -α C ξ C (x(t)).
(2)

The trajectories x(t) are determined by two orthogonal vector fields ξ J and ξ C referred to as the "null space" and "range space" directions, following denominations of null space optimization methods for equality constrained problems, see e.g. [START_REF] Nocedal | Numerical optimization, 2nd edn[END_REF]; Nie (2004); [START_REF] Liu | A null-space primal-dual interior-point algorithm for nonlinear optimization with nice convergence properties[END_REF]. The null space step -ξ J is the best locally feasible ascent direction. The range space step -ξ C is a Gauss-Newton direction that gradually corrects the violation of the constraints.

The precise definitions of ξ J and ξ C are reviewed in Sections 2.2 and 2.3 below. All in all, (2) is a generalization of the standard gradient flow ẋ = -∇J(x) minimizing J without constraints. Similarly to the unconstrained case, the trajectories x(t) of (2) tend to converge to feasible designs x * that are locally optimal in the sense that they satisfy the first order Karush-Kuhn-Tucker (KKT) conditions, namely x * is feasible and

∃λ ∈ R p , ∃µ ∈ R q + , ∇J(x * ) + Dg(x * ) T λ + Dh(x * ) T µ = 0, µ j > 0 ⇒ h j (x * ) = 0 for all 1 ≤ j ≤ q, (3)
where R q + denotes the q-dimensional nonnegative orthant and ∇J(x * ), Dg(x * ) and Dh(x * ) are respectively the gradient of J and the Jacobian matrices of g and h at x * .

The ODE (2) is parameterized by two positive constants α J > 0 and α C > 0. The parameter α J tunes the rate at which the objective function values decrease while not worsening the constraints.

The parameter α C tunes the pace at which the violation of the constraints decreases (the violation decreases precisely at the exponential rate e -α C t ). In principle, the success of the method to find a local minimizer does not depend on the choice of α J and α C : it depends only on the selection of a sufficiently small time step. However, setting the ratio α J /α C to a custom value sometimes helps to find a better optimization path because it tunes the trade-off between guiding quickly the design to a region with small objective function values and driving it quickly to the feasible set. In any case, we emphasize that these parameters have a clear intuitive meaning and in many situations (including those considered in this paper), it is possible to set α J /α C = 1. For a more flexible usage, the Python implementation of the Null Space Optimizer implements automatic rescaling rules on the parameters α J and α C , see Section 2.5.

Remark 1. The Null Space algorithm is not restricted to optimization on finite-dimensional spaces but also allows to consider much more general sets equipped with a kind of manifold structure, see [START_REF] Feppon | Shape and topology optimization of multiphysics systems[END_REF], chapter 6. This feature is essential for level-set based Topology Optimization, where the design set is the set of shapes X = {Ω ⊂ D} for a given computational domain D.

Remark 2. In (3) and throughout the paper, we denote by DJ(x) and Dg(x) the Fréchet derivatives of a function J and a vector g, where DJ(x) := ∂ x1 J(x) . . . ∂ xn J(x) is understood as a row vector, and where Dg(x) := (∂ j g i (x)) is the Jacobian matrix of g. We use the gradient notation to refer to the transpose ∇J(x) := DJ(x) T . This distinction between Fréchet derivative and gradient is needed in the Hilbertian framework of level-set based Topology Optimization; see Feppon et al (2020a), Definition 2.1.

Definition of the null space direction ξ J

In what follows, || • || denotes the standard Euclidean norm. For a given design x ∈ X , we consider the set I(x) of inequality constraints that are active or violated, and we denote respectively by q(x) and h I(x) their number and the associated vectors of constraints:

I(x) = {i |h i (x) ≥ 0}, (4)
q(x) := Card( I(x)), (5)

h I(x) := (h i (x)) i∈ I(x) . (6)
The null space direction ξ J (x(t)) is defined by the formula

ξ J (x) := ∇J(x)+Dg(x) T λ * (x)+Dh I(x) (x) T µ * (x), (7) where (λ * (x), µ * (x)) ∈ R p × R q(x) +
are optimal Lagrange multipliers obtained by solving the "dual" problem min

(λ,µ)∈R p ×R q(x) + ||∇J(x)+Dg(x) T λ+Dh I(x) (x) T µ||, (8)
In practice, ( 8) is a nonnegative least-squares problem which can be solved by using any standard quadratic programming solver; the current Python implementation of the Null Space Optimizer allows to choose between the CVXOPT [START_REF] Vandenberghe | The CVXOPT linear and quadratic cone program solvers[END_REF]) and the OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]) solvers.

Remark 3. The problem (8) has a clear and intuitive interpretation. Either its minimum value is zero and the KKT optimality conditions (3) are satisfied; or it is not zero and the vector ξ J (x) obtained from (7) with the dual variables (λ * (x), µ * (x)) is the "best" locally feasible descent direction. More precisely, ξ * (x) := -ξ J (x)/||ξ J (x)|| minimizes the linearized objective function while not worsening the first order variation of the constraints among all possible unit directions:

- ξ J (x) ||ξ J (x)|| = arg min ξ∈R n DJ(x) • ξ s.t.      Dg(x) • ξ = 0, Dh I(x) (x) • ξ ≤ 0, ||ξ J (x)|| ≤ 1. ( 9
)
This alternative is the generalization to the well-known fact that in the unconstrained case, ∇J(x) is either zero at a local minimum or -∇J(x)/||∇J(x)|| is the best unit descent direction. The problem (8) is called "dual" because it is the dual of (9) [START_REF] Nocedal | Numerical optimization, 2nd edn[END_REF]).

Definition of the range space direction ξ C

The range space direction ξ C (x(t)) is a Gauss-Newton direction for reducing the violation of the constraints, thus gradually correcting unfeasible initializations. It is mathematically defined as

ξ C (x) := DC T I(x) (DC I(x) DC T I(x) ) -1 C I(x) (10)
where

C I(x) := g(x) h I(x)
is the column vector gathering the equality and active inequality constraints. The definition (10) ensures that DC I(x) ξ C (x) = C I(x) . By using the properties Dg(x)•ξ J (x) = 0 and Dh I(x) •ξ J (x) ≤ 0, it is found that the violation of the constraints decreases exponentially as

g(x(t)) = e -α C t g(x(0)), max(h(x(t)), 0) ≤ max(h(x(0)), 0)e -α C t ,
see Feppon et al (2020a).

Remark 4. Formula (10) assumes DC I(x) to be full rank for the matrix DC I(x) DC T I(x) to be invertible. In order to address the issue of rank degeneracy, the Python implementation of the Null Space Optimizer uses the SciPy routine lsqr to compute the right pseudo inverse of DC I(x) in place of DC T I(x) (DC I(x) DC T I(x) ) -1 .

Treatment of constraints with sparse Jacobian and bound constraints

A naive implementation of the quadratic norm in (8) requires to compute the large matrix DC I(x) DC T I(x) . This matrix needs also to be inverted in the formula (10). When the total number p + q(x) of active or violated inequality constraints is large, performing these operations explicitly may be computationally intractable. However, it is possible to reformulate (8) and ( 10) in terms of matrix-vector products that lend themselves to iterative methods when the constraints have a sparse Jacobian matrix. The necessary adaptations are detailed in Appendix C.

Bound constraints are a special instance of constraints with sparse Jacobian matrix. With the upgrade of Appendix C, the Null Space Optimizer treats them as any other constraints, by simply including their Jacobian matrix in the assembly of DC I(x) . For instance, bound constraints of the form l i ≤ x i ≤ u i with 1 ≤ i ≤ n can be reformulated into 2n inequality constraints l i -x i ≤ 0 and x i -u i ≤ 0. The associated Jacobian is the block matrix

DC = I -I ,
where I is the n × n identity matrix. For the ease of the user, the Nullspace Optimizer Python package provides a decorator @bound constraints optimizable for automating this assembly which is described in more details in Section 3.2.4. Additionally, the decorator ensures that the design variable is projected to the admissible set at every iteration by using the min/max operator

x → min(u i , max(l i , x)). (11)
This operation ensures a more conservative enforcement of the bound constraints.

Implementation aspects and variations of the algorithm

The gradient flow (2) is discretized by using a Forward Euler scheme with time step ∆t. The Python implementation of the Null Space Optimizer introduces further adaptations for improved numerical stability and usage flexibility.

Automatic rescaling of α J and α C

The parameters α J ≡ α n J and α C ≡ α n C are automatically scaled at every iteration n in order to ensure

||α n J ξ J (x n )|| ∞ ∆t ≤ A J ∆t, ( 12
)
||α n C ξ C (x n )|| ∞ ∆t ≤ A C ∆t, (13)
where x n is the iterate after n time steps, A J and A C are two fixed constants set by the user and || • || ∞ is the supremum or L ∞ -norm (thus not the Euclidean norm || • || of Section 2.2). Furthermore, the parameter α n J is chosen such that (12) holds as an equality during for the first N it iterations:

||α n J ξ J (x n )|| ∞ ∆t = A J ∆t, for all 0 ≤ n ≤ N it . (14)
The number N it is fixed by the user (see the correspondence table at the end of this subsection). The equality is not enforced afterwards in order to let ξ J (x n ) → 0 as x n gets closer to a local minimum. These rescalings are quite convenient since A J and A C can be understood as fractions of the time step. They have therefore a clear intuitive meaning for the user who does not need to worry about the magnitudes of ξ J (x n ) and ξ C (x n ).

Stabilization near active inequality constraints

Eq. ( 2) is an ODE with discontinuous right-hand side, which may introduce numerical instabilities near the boundaries of inequality constraints. When a novel constraint becomes active in the set I(x), the gradient of the objective function is suddenly projected on a smaller set when crossing the constraint barrier, resulting in a discontinuous change in the orientation vector of the trajectory. These instabilities can be addressed by enlarging the set I(x) to include inequality constraints {0 ≤ j ≤ q | h j (x) ≤ ϵ j } that are not yet but likely to become active when computing the null space direction ξ C (x) with ( 7) and ( 8). The parameter ϵ j are small threshold values that are automatically computed to detect constraints located at a distance less than K∆t, where K is a parameter set by default to 0.1 (which is working well in practice independently of the optimization problem). Then, the set I(x) considered in (10) for computing the range step ξ C (x) is also enlarged with constraints that are likely to remain active based on the information given by the optimal multiplier µ * . We refer the reader to Feppon et al (2020a) for the detailed implementation of these rules.

Weights for the violated constraints

The direction ξ C (x) in ( 10) is designed to ensure that the violation of the constraints decays as e -α C t . In practice, the user may desire different decay rates for each of the constraints for a more flexible control of the optimization path. This is possible by defining ξ C (x) as [START_REF]Fh0 dkappa = p * rho ^( p -1)*( kappaf -kappas )[END_REF] where (α i ) 1≤i≤p+q are p + q weights ensuring that the violation of the corresponding equality or inequality constraint i decreases faster than e -αiα C t .

ξ C (x) := DC T I(x) (DC I(x) DC T I(x) ) -1 × diag(α 1 , . . . , α p , α p+1 , . . . , α p+q )C I(x)

Correspondence table for the algorithm parameters

The Null Space Optimizer accepts various algorithm parameters in the form of a Python dictionary provided by the optimization routine (see e.g. Section 3.3 below). The following table provides the correspondence between the mathematical notation of the parameters mentioned in this section and their name (as dictionary keys) in the Python implementation.

Parameter

Implementation Default ∆t dt 0.1

A J alphaJ 1 A C alphaC 1 N it itnormalisation 1 K K 0.1 (α i ) 1≤i≤p+q alphas (1) 1≤i≤p+q
Note that these algorithm parameters are not tuning parameters, they are rather used to provide an increase level of flexibility allowing for faster convergence towards an optimum. They have a clear physical meaning and many practical situations do not require a custom value. In the numerical examples considered in this paper, we leave the parameters to the default values unless otherwise specified.

Topology Optimization of the 88 lines Matlab MBB beam example with the Null Space Optimizer

This section and the next to follows provide three tutorials giving an in-depth view on the use of the Null Space Optimizer for density based Topology Optimization. In this section, we revisit the classical MBB beam example from the popular 99 lines Matlab Topology Optimization code originally introduced in Sigmund ( 2001), which was refactored in a more efficient 88 lines version in [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF].

We provide a step-by-step tutorial for writing an elegant implementation of this classical example with the Null Space optimizer, without bringing substantial modifications to the state solver and the sensitivities. The proposed implementation is claimed "elegant" because it uses elegant paradigms of the Python language such as object-oriented programming and decorators to isolate the various building blocks of the Topology Optimization code. As the following examples demonstrate, we obtain as such succinct and flexible implementations that closely matches the mathematical formulation (1) of the optimization problem, sets a clear separation between the implementation of the state solver, the filtering of the design variable and the optimization solver. This programming paradigm is easily reproducible to more general situations. This section is organized as follows. The physical setting of the MBB beam test case and the formulation of the associated design optimization problem is given in Section 3.1. Section 3.2 details how to implement this test case in Python and to solve it with the Null Space Optimizer. Numerical results and comparisons with the Optimality Criteria, MMA and IPOPT solvers are presented in Section 3.4 for two different benchmark situations: one with a rather large target volume fraction corresponding to the setting of the 88 line Matlab code paper of Andreassen et al ( 2011), and a second example imposing a tiny target volume fraction that is more challenging.

The classical MBB beam

problem of the 88 lines Matlab Topology Optimization code. Throughout the paper, the SIMP exponent is set to p := 3, Emax := 1 is the Young modulus of the solid material and Emin := 10 -9 is the Young modulus of an ersatz material representing void. The left-side boundary is subject to no horizontal displacement and the bottom right corner to no vertical displacement. A unit downward force g is applied at the top left corner.

The structural design optimization problem for minimizing the compliance of the structure subject to a single volume equality constraint reads:

min x=(xe) 1≤e≤n J(x) = U T K(x)U = n e=1 E e (x e )u T e KEu e s.t.            V (x)/V 0 := 1 n n e=1 x e = volfrac, K(x)U = F , 0 ≤ x e ≤ 1, for all 1 ≤ e ≤ n, (16)
where K(x) is the global stiffness matrix, U and F are the global displacement and force vectors, u e is the displacement vector in the element e and KE is the element stiffness matrix for an element with unit Young's modulus. V (x) refers to the volume of the optimized design, V 0 is the total volume of the design domain and volfrac is the target volume fraction.

A classical issue in density based TO when solving the minimization problem ( 16) is the arising of checkerboard patterns in the design variable [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, meshdependencies and local minima[END_REF]). In order to eliminate these nonphysical solutions from the design space, the original 88 lines Matlab Topology Optimization code implements a density filter on the optimization variable. More precisely, the final optimization problem considered reads min

x=(xe) 1≤e≤n J( x) s.t.      V ( x)/V 0 = volfrac, K( x)U = F , 0 ≤ x e ≤ 1 for all 1 ≤ e ≤ n, (17)
where the filtered density x is a local averaging of elements that are neighbors by a distance smaller than a prescribed filter radius r min :

x e := 1 n i=1 H ei n i=1 H ei x i for all 1 ≤ e ≤ n, (18)
with H ei = max(0, r min -∆(e, i)), where ∆(e, i) is the distance between the centers of elements e and i.

Python implementation of the MBB beam problem with the Null Space Optimizer

In order to closely match the original Matlab implementation of this test case, we reshape a Python translation of the 88 lines code of [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF] that is provided by [START_REF] Aage | A 200 line topology optimization code[END_REF]. In the subsequent subsections, we provide a comprehensive description of the implementation using the Null Space Optimizer. The full source code is provided in Appendix A and we refer to the associated line numbering throughout the discussion.

The code can be comprehended in the following parts:

• lines 1-9: Python modules imports;

• lines 11-109: initializations of the element stiffness and density filter matrices; • lines 111-142: declaration of the state solver;

• lines 144-152: declaration of the filtering functions; • lines 154-191: definition of the optimization problem • lines 199-203: setting of optimization parameters and running the optimization.

Declaration of the optimization problem

The overall code is best understood by looking first at the declaration of the class TO problem in lines 154-191. The class TO problem inherits the EuclideanOptimizable class from the Null Space Optimizer Python package, which indicates to the optimization routine that the design domain is a finite dimensional design space R n (represented by NumPy arrays). The declaration of this class closely follows the mathematical definition of the Topology Optimization problem ( 17) by implementing the following methods (required by the Null Space Optimizer):

• x0 returns the initial design;

• J returns the objective function computed on the design variable x; • G returns the array of equality constraints, here the single volume constraint V ( x)/V 0volfrac; • dJ returns the sensitivity of the objective function (a row vector); • dG returns the Jacobian matrix of the equality constraints (here a row vector); Arbitrary inequality constraints could be added similarly by implementing methods H and dH. Additionally, the class implements a function accept for initializing and updating plots of the design. The function accept is called at the beginning of every iteration by the Null Space Optimizer, it takes as inputs a dictionary params containing the optimization algorithm parameters and a dictionary results storing the optimization histories. The optimized design from the current iteration is accessed as results

['x'][-1].
The implementation of the Topology Optimization problem relies then on three key ingredients:

• solve state (lines 111-142): a function solving the linear elasticity problem and returning the compliance obj and its sensitivity dc, • @filtered optimizable (line 156): a decorator which implements the filtering of the design variable, • @bound constraints optimizable (line 155): a decorator which supplements the optimization problem with the bound constraints 0 ≤ x e ≤ 1 for all 1 ≤ e ≤ n. These three ingredients are described in the next Sections 3.2.2 to 3.2.4.

Computations of the state solution and the sensitivities

The computation of the compliance and its sensitivity is achieved by the function solve state at lines 111-142. The function depends on a preliminary part of the code (lines 11-109) which precomputes the element stiffness matrix KE, the index pointers iK, jK as well as the element stiffness matrix edofMat describing the degrees of freedom attached to every element. The array dofs references the finite element degrees of freedom attached to every vertex. The total number of degrees of freedom is stored in the variable ndof, the degrees of freedom associated to the Dirichlet and the free boundary are stored in the arrays named fixed and free, and the downward unit load is stored in the finite element vector f (which implements the right-hand side F of the linear elasticity system).

Then, the function solve state takes as an input an array x representing the density variables x = (x e ) 1≤e≤n . It internally assembles the finite element matrix and solves the linear system K(x)U = F by applying a sparse Cholesky factorization provided by the CVXOPT package. This part follows almost textually the Python translation of the 88-lines Matlab code of [START_REF] Aage | A 200 line topology optimization code[END_REF] and the reader is referred to [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF] for more details.

However, a few minor useful changes have been implemented. First, the function solve state allows the force vector f to be a matrix with N loads columns corresponding to different loads. This parameter is useful for the multiple load case of Section 4 below. In the present situation, f is a single column vector. The function returns an array obj for multiple loads N loads > 1, and the real number obj[0] when N loads = 1 (the present case).

Second, the NumPy function np.einsum is used at line 131 for efficient and vectorized computation of the terms u T e KEu e in [START_REF][END_REF]. Third, the definition of the function solve state is preceded by the decorator @memoize() provided by the Null Space Optimizer package. This decorator caches the function's results, effectively returning a stored value if the computationally intensive solve state function is invoked multiple times with the same argument. This feature is particularly advantageous in the declaration of the optimization problem which requires to separation of the objective function and the sensitivities into distinct functions, resulting in separate calls at lines 162 and 166: this syntax ensures that the linear elasticity system is solved only one time per iteration. The decorator can be understood as equivalent to appending the following lines of codes after the definition of solve state:

old_solve_state = solve_state solve_state = \ memoize ()( old_solve_state )
The function memoize() relies on a hashing algorithm that determines whether the function solve state has already been called on the argument x, and in that case, returns a stored value.

Filtering of the density variable

In the proposed implementation, the density filter x → x of ( 18) is conveniently implemented at line 156 by the decorator @ f i l t e r e d _ o p t i m i z a b l e ( filter , diff_filter )

Its arguments filter and diff filter respectively implement the filtering function ( 18) and its sensitivity.

The decorator @filtered optimizable converts the problem ( 16) into the filtered problem (17). It automatically creates a copy of the class TO problem, substituting the methods J(self,x) and G(self,x) with J(self,filter(x)) and G(self,filter(x)).

The filtering function filter is implemented at lines 144-152 with the following lines of code: @memoize () def filter ( x ):

return Hs @ ( H @ x )

where H is the matrix (H ei ) of ( 18) and Hs is the sum of its elements; H and Hs are precomputed at lines 80-105 reproducing the implementation of Aage and Johansen (2016); [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF]. Memoization is used once again for efficiency, which allows for the proposed implementation paradigm without the need to worry about redundant calls to the filter function.

To obtain fully consistent sensitivities, the derivatives need to be updated according to the chain rule

d dx [J( x)] • dx = dJ( x) • d x dx • dx.
The decorator @filtered optimizable automatically applies this transformation to the methods dJ and dG by exploiting the second argument diff filter, which is required to be a Python function implementing the left vectormatrix product:

(x, v) → v • d x dx , ( 19
)
where v is an input row vector (or a p×n matrix if G is a function returning p constraints) that is to be replaced with dJ( x) or dG( x) by the optimization routine. In the present case, the filter is linear which makes the implementation of diff filter straightforward (lines 143-146):

@memoize () def diff_filter (x , v ): return ( H @ ( Hs @ v . T )). T

Bound constraints

The 2n bound constraints 0 ≤ x e ≤ 1 for 1 ≤ e ≤ n are implemented with the decorator @bound constraints optimizable(l=0,u=1) at line 155. This decorator returns a copy of the filtered class TO problem, by: (i) appending a method H(self,x) implementing the 2n bound constraints -x e ≤ 0,

x e ≤ 1, for 1 ≤ e ≤ n.
This operation is equivalent to adding a method H(self,x) in the definition of the class TO problem4 :

def H ( self , x ): # Returns the vector # of constraints # [ x_0 -1 , x_1 -1 ,... , x_n -1 , # -x_0 , -x_1 ,... , -x_n ] return np . hstack (( x -1 , -x ))
The use of the decorator allows to simplify this process for the user's convenience; (ii) appending a method dH(self,x) returning the constraints sensitivities, that is the block diagonal matrix

dH(x) = I -I
where I is the n × n identity matrix, This function returns this Jacobian matrix in a sparse format. This operation is equivalent to adding the following method in the class TO problem:

def dH ( self , x ): I = sp . eye (n , format = " csc " ) return sp . vstack (( I , -I ))

(iii) adding a method retract(self,x,dx) implementing the projection rule (11) that is applied after every update by the Null Space Optimizer. We refer to [START_REF] Absil | Projection-like Retractions on Matrix Manifolds[END_REF] and to the documentation of the package regarding the definition of retractions and their use in optimization. Notice that unlike MMA or IPOPT, the Null Space optimizer treats bound constraints almost as any other inequality constraints. The only difference in this treatment with other constraints is the use of the projection rule (11) that enables a stricter enforcement than the range space step ξ C .

Running the optimizer

The optimization problem is solved with the Null Space Optimizer at lines 193-203. Line 200 initializes finite element matrices by calling the function init(). The class TO problem is instantiated at line 195. Then, the optimization problem is solved at line 201 by the routine nlspace solve by discretizing the null space gradient flow (2). The arguments of the routine are an instance of the optimization problem and the dictionary of algorithm parameters optimization params. Here, we specify a time step dt = 0.3. The parameter itnormalisation=50 ensures that the null space step α J ξ J (x) is rescaled to unity during the first 50 iterations, before being allowed to smoothly decrease (see Section 2.5). The parameter maxit stops the optimization after 150 iterations. The parameters save only N iterations and save only Q constraints are used for memory efficiency: by default, the Null Space Optimizer stores the values of all constraints and designs of all iterations in the output dictionary results. In the present case, we require that only the last iteration and the first 5 constraints (which include the volume constraints and the first 4 bound constraints) are stored.

Numerical results and comparisons with OC, MMA and IPOPT

We solve the MBB Topology Optimization problem with the Null Space optimizer and we compare the convergence to that of the Optimality Criteria method, MMA and IPOPT. We consider two situations: first, the exact configuration of the MBB problem considered in the 88 lines Matlab code (Section 3.4.1), then a more challenging version where the volume fraction target is set to a very small value (Section 3.4.2).

MBB beam test case with medium volume fraction

In this first situation corresponding to the setting of [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF], the volume fraction target is set to volfrac:=0.4 and the filter radius to r min := 5.4. The number of elements in the x and y directions are set to nelx := 180 and nely := 60. The Null Space optimizer time step dt, the move limit m OC of the Optimality Criteria method and the move limit m MMA of MMA are set to the same value dt = m OC = m MMA = 0.3. The number of normalisation iterations for the null space direction ξ J is set to itnormalisation=50. We keep the default settings for IPOPT. The optimization is stopped for every algorithm after 150 iterations.

The final objective and constraint values are reported in Table 1. Convergence histories and intermediate designs at identical iterations are displayed respectively on Figs. 2 and4. A close-up on the convergence histories for the first 30 iterations is shown on Fig. 3.

On this basic example, all the algorithms obtain the same final design with very similar performances. The Null Space Optimizer (abbreviated NLSPACE) performs similarly well as its competitors, the Optimality Criteria method, MMA and IPOPT. With the chosen settings, NLSPACE manages to enforce the volume constraint at machine precision. The performance of the design found by MMA is slightly better, while the performance of IPOPT and OC are almost identical to that of the Null Space Optimizer.

IPOPT takes slightly more iterations to converge than the other algorithms and exhibits a rather surprising behavior around iteration 20, it still manages to drive the design to the optimum after a few more iterations. This analysis, however, must be tempered by the fact that we didn't tune IPOPT parameters to find the most efficient convergence setting. It can also be observed that MMA takes a rather large deviation to the volume constraint during the first 20 iterations, while OC maintains it close to zero throughout the iterations. The Null Space Optimizer as well as IPOPT only slightly violate the constraint during the first iterations before correcting it after approximately 20 iterations.

Looking at the first iterations on the convergence histories of Fig. 3, we observe that the Optimality Criteria method, MMA and the Null Space Optimizer have similar initial decays, with a small advantage for OC and MMA. This effect is also visible on the other numerical examples considered in the next sections. This small advantage may be surprising because we used identical move and time step dt parameters while the definition (9) of the null space step ξ J should ensure that it produces the best locally feasible descent direction, hence theoretically leading to best update among all possible directions with identical norm. This apparent paradox is explained by the fact that the norm || • || involved in the definition (9) of ξ J (x) is a quadratic norm, while the implementation normalizes the update steps with the supremum (L ∞ ) norm ||•|| ∞ (see Section 2.5). Here, MMA and OC use descent directions that are more optimal than ξ J in the L ∞ -norm, but we checked that these descent directions are of strictly larger L 2 -norm than ξ J . Still, this advantage remains negligible on this test case as NLSPACE catches up MMA and OC after about 7 iterations.

Algorithm J G = V (x)/V 0 -

MBB beam test case with a tiny volume fraction

We now investigate a more challenging example where the volume fraction is set to the small value volfrac := 0.01. For this test case, the resolution parameters of the computational mesh are set to nelx:=150 and nely:=50, and the filter radius to r min := 1.9. As the numerical results demonstrate below, the test case is somehow unrealistic since a completely black and white design cannot be achieved for such a low amount of material given the resolution. Still, the test case is an interesting challenging benchmark for testing the performance of optimization methods.

For such a small value of the volume fraction target, the classical implementation of MMA based on the Matlab code released by Svanberg (2014) fails both at finding an optimum value and ensuring the correct volume fraction, see Fig. 5 and Table 2. Rescaling the objective function is not sufficient for ensuring convergence towards a satisfactory design. The issue comes from what seems to us a limitation of the heuristic update rules of the moving asymptotes. Indeed, a suitable convergence can be obtained by the following for k = 0, 1. This setting avoids that the lower asymptotes to become lower than 0, and it "pushes" the upper asymptotes to 10 making the initial approximation of the objective and constraint functions flatter; • discarding the moving asymptote rule, which reads (eq. (3.12) in Svanberg ( 2014)): We note that the use of the line-search provided by the Globally Convergent MMA (GCMMA) algorithm of [START_REF] Svanberg | MMA and GCMMA -two methods for nonlinear optimization[END_REF] would also be an appropriate fix for this problem. Still, even GCMMA would require the smart initialization of the asymptotes to lead to satisfactory designs.

l (k) e = x (k) e -γ (k) e (x (k-1) e -l (k-1)
We run the different optimizers with dt=0.02 for the Null Space Optimizer, move=0.2 for MMA and move=0.05 for the OC solver. These move parameters yield satisfactory results for MMA and OC. The NLSPACE time step is the sole parameter to tune, it needs to be taken sufficiently small to accommodate the small volume fraction and the initial design (a uniform gray design x 0 = (x 0 e ) with x 0 e = 0.01/(nelx × nely) for 1 ≤ e ≤ n) close to 0. IPOPT is executed with default options. The numerical values of the objective function and volume constraint of the final designs are listed in Table 2. Convergence histories are shown on Fig. 6 as well as a close-up on the first 30 iterations on Fig. 7. Figure 8 displays a sequence of intermediate designs for each optimizer.

As it is visible on Fig. 8, the target volume fraction could not be achieved with completely black and white designs, leading "gray" solid bars having a density that do not reach the maximum value x e = 1. The numerically computed optimized designs suggest that this benchmark example seems to admit several locally minimizing designs with different topologies achieving comparable performances. The values of Table 2 indicate that the fixed version of MMA happens to find a design with the best performance among the tested algorithms. The Optimality Criteria solver does not need any fix and is able to find a slightly different design with similar performances. However, this solver is less accurate regarding the enforcement of the volume constraint as illustrated on Fig. 6b. The Null Space algorithm converges to a different design with a performance very similar to the one found by the fixed version of MMA, albeit slightly worse and slightly better than OC. Noticeably, it enforces the volume constraint at machine precision at convergence. On this example, the IPOPT solver with default parameters converges to a design that satisfies the constraint, but that is much less performant than the solutions provided by the other successful algorithms.

Algorithm J G = V (x)/V 0 -

Implementing Robust Topology Optimization for multiple load cases

In this section, we describe the implementation of robust design optimization of a bridge subject to multiple loads using a min/max formulation. We demonstrate the suitability of the Null Space Optimizer for solving this problem featuring multiple inequality constraints, and we compare the numerical results to those of MMA and IPOPT (for this example, it is not possible to use the standard Optimality Criteria method which is adapted to optimization problems featuring only a single equality constraint). The Null Space Optimizer has already proven to be effective for multiple load case Topology Optimization with the levelset method: a similar bridge test case has been considered in Feppon et al (2020a). We find that the performance of the method is not altered by adding the bound constraints that are peculiar to density-based Topology Optimization.

The physical setting and the formulation of the optimization problem are outlined in Section 4.1. Section 4.2 describes the Python implementation of the multiple load case with the Null Space Optimizer, that is obtained after only a few modifications of the MBB beam test case of Section 3. A decorator @minmax optimizable is provided by the Python package to conveniently implement the robust min/max formulation of the optimization problem. Finally, numerical results comparing the Null Space Optimizer, MMA and IPOPT are given in Section 4.3.

Setting of the bridge Topology Optimization problem with multiple loads

The bridge test case that we consider is identical to that of the MBB problem illustrated on Fig. 1, up to the following changes:

• the design domain features now nelx × nely elements with nelx=180 and nely=90; • the bottom right corner is now fixed both in the horizontal vertical direction; g 1 g 2 g 3 g 4 g 5 g 6 g 7 g 8

Fig. 9: Setting of the bridge design problem of Section 4.

• that eight different unit loads (g i ) 1≤i≤8 are independently applied on the top boundary. Here, the optimized domain actually represents half of the final structure, assuming symmetry with respect to the y-axis; the whole structure is thus subject to 8 loads consisting of two vertical forces distributed symmetrically. Each load spreads continuously over one eighth of the top boundary with a magnitude set to -1/nelx. The setting of the multiple load case is illustrated on Fig. 9.

The goal is to minimize the maximum of the compliances of the structure subject to the 8 possible loads and a prescribed volume fraction of material. This robust design formulation corresponds mathematically to solving the min/max optimization problem min

x=(xe) 1≤e≤n max 1≤i≤8 J i ( x) := U T i K(x)U i s.t.      V ( x)/V 0 = volfrac, K( x)U i = F i , for all 1 ≤ i ≤ 8, 0 ≤ x e ≤ 1 for all 1 ≤ e ≤ n, ( 20
)
where x is the filtered density ( 18), F i is a column vector corresponding to the finite element representation of the i-th load g i , and U i is the finite element representation of the associated displacement vector.

In order to treat (20) without the need for differentiating the maximum function, we solve the classical equivalent formulation that introduces a slack variable m ∈ R:

min x=(xe) 1≤e≤n , m∈R m s.t.          J i ( x) -m ≤ 0, for all 1 ≤ i ≤ 8, V ( x)/V 0 = volfrac, K( x)U i = F i , for all 1 ≤ i ≤ 8, 0 ≤ x e ≤ 1 for all 1 ≤ e ≤ n. (21)

Python implementation with the Null Space Optimizer

Thanks to a few automations provided by the Null Space optimizer Python package, the implementation of the problem ( 21) requires only minor updates of the code of the MBB beam test case of Appendix A.

First, the variable nely can now be assigned by default to 90 by updating line 50: nely = kwargs . get ( ' nely ' , 90)

Then, the vector fixed containing the finite element degrees of freedom corresponding to zero Dirichlet boundary conditions needs to be updated to match the situation of Fig. 9. This is done by updating line 77 as follows:

fixed = np . union1d ( dofs [0 ,: ,0] , dofs [ -1 , -1 ,:])
Second, we assemble the matrix f= F 1 . . . F 8 whose columns contain the finite element representation of the 8 different loads. This can be achieved by substituting lines 108-109 with the following code which determines the degrees of freedom at the top boundary corresponding to each load g i and assigns them the force -1/nelx:

nforces = 8 n = np . ceil ( nelx / nforces ) positions = np . array (\ [ k * n for k in range ( nforces )]\ +[ nelx ] , dtype = int ) f = np . zeros (( ndof , nforces )) for i in range ( nforces ): start = positions [ i ] end = positions [ i +1] f [ dofs [ start : end ,0 ,1] , i ] \ = -1.0/ nelx
Since the function solve state of lines 111-142 is designed to allow for such matrix argument f, the (unfiltered) method J(self, x) at line 161 returns without any further change the vector (J 1 (x), J 2 (x), . . . , J 8 (x)), and the (unfiltered) method dJ(self, x) at line 165 returns the 8

× n matrix      dJ 1 (x) dJ 2 (x) . . . dJ 8 (x)      , ( 22
)
where dJ i (x) is the sensitivity of the i-th compliance.

Third, the multi-objective optimization problem needs to be converted into the min/max formulation ( 21). This can be conveniently achieved using the decorator @minmax optimizable provided by the package, which needs to be imported at the beginning of the code as follows:

from nu l l s pa c e _o p t im i z er import mi nm ax _o pt im iz ab le

Then, the conversion of the multi-objective problem to the min/max formulation is implemented with a single line of code by inserting the decorator before the definition of the class TO problem just before line 155: The decorator @minmax optimizable provided by the Null Space Optimizer Python allows thus the user to implement the robust formulation (20) in a multi-objective fashion with the multivalued objective function J(self,x). This paradigm could be easily extended to further multi-objective formulations, for instance it would not be difficult to implement a decorator @weighted optimizable converting TO problem into a class corresponding to the minimization of a weighted average of the objective functions.

In the present case, the decorator @minmax optimizable returns a copy of the class TO problem with the following changes:

• the method x0(self) returns (x 0 , max i (J i (x 0 ))) where x 0 is the output of the original method x0 of TO problem;

• the method J(self,x) returns m ≡ x[1]
since the optimization variable is now a tuple x ≡ (x, m); • a method H(self,x) is added to return the inequality constraints (J i (x) -m) 1≤i≤8 ; • methods dJ(self,x) and dH(self,x) are added for automatically returning sensitivity Jacobian matrices updated with respect to the variable m. dJ(self,x) simply returns the n + 1 dimensional row vector 0 . . . 0 1 . The method dH(self,x) returns the 8

× (n + 1) matrix      dJ 1 (x) -1 dJ 2 (x) -1 . . . . . . dJ 7 (x) -1      . ( 23
)
• finally, a method retract(self,x,dx) is added, which implements the projection rule

m ← max 1≤i≤8 J i (x n+1 ).
This projection is executed by the Null Space Optimizer after each update of the design variable. Although not fully necessary, this ensures a more conservative enforcement of the inequality constraints J i (x) ≤ m than the gradual correction of the range step ξ C .

Numerical results and comparisons with MMA and IPOPT

We run the Null Space Optimizer on the bridge test case with the time step dt=0.1 and the setting itnormalisation=50. We also solve the design test case with the standard implementation of MMA (with move=0.1) and IPOPT (with default parameters). The final values of the compliance for all the load cases and of the volume constraint are reported in Table 3. Convergence histories of the objective and volume constraint are shown on Fig. 11; the whole bridge is represented by symmetrizing the numerically computed design for a better visualization.

G NLSPACE MMA IPOPT (b) Volume constraint G(x) = V (x)/(V 0 volfrac) -1.
On this example, the three optimizers tested converge to three different designs with similar performances. Noticeably, MMA takes a significantly distinct optimization path leading to a very different geometry while IPOPT and NLSPACE find designs that are similar up to some small differences at the middle of the structure. Here, the optimum found by MMA is slightly better than the one found by the Null Space Optimizer, which is itself slightly better than the one computed by IPOPT. It seems that there is a bifurcation between two possible local minimizers and that for some unclear reason, MMA selects the path leading to a better optimum. We observed that MMA

Heat sink Topology Optimization on unstructured meshes

This section presents an additional tutorial for implementing Topology Optimization of 2D heat sinks. This classical test case was initially considered by [START_REF] Bejan | Constructal-theory network of conducting paths for cooling a heat generating volume[END_REF], and then solved with density based Topology Optimization in [START_REF] Gao | Topology optimization of heat conduction problem involving design-dependent heat load effect[END_REF]; [START_REF] Zhang | Design of conducting paths based on topology optimization[END_REF]; [START_REF] Marck | On some shape and topology optimization problems in conductive and convective heat transfers[END_REF]; [START_REF] Yan | On the nonoptimality of tree structures for heat conduction[END_REF].

The test case substantially deviates from the setting of the 88 lines Matlab TO paper of [START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF] considered in the previous parts, because we consider a different physics, unstructured mesh discretization and the use of continuous sensitivities rather than discrete ones.

The physical setting of the heat sink design problem is described in Section 5.1. A step by step description of the implementation of this test case within the paradigm of the Null Space Optimizer Python package is outlined in Section 5.2. We describe how to implement the physical state solver with FreeFEM [START_REF] Hecht | New development in FreeFem++[END_REF]) and a Laplacian smoothing density filter is implemented to avoid the checkerboard effect on triangular meshes. The interfacing of the optimizer with the Finite Element software FreeFEM is especially advantageous as it could be used to deal with more challenging physics or for implementing density based TO on computational domains with arbitrary geometries. Finally, Section 5.3 presents numerical results of the proposed implementation including a comparison with the OC, MMA and IPOPT optimizers.

The content of this section partly follows the work of [START_REF] Krasniqi | Density-based topology optimization with the Null Space optimiser[END_REF] where the same test case was solved using a simpler extension of the Null Space algorithm for dealing with bound constraints. (namely, which did not consider an arbitrary inner product A as described in Appendix C).

Setting of the heat Topology Optimization problem

We consider a square domain D = (0, 1) × (0, 1) constituted of a conductive material with temperature field T subjected to a uniform heat source temperature on D, with a prescribed maximum volume fraction volfrac for the most conductive material κ f . We use an "optimize-then-discretize" density based Topology Optimization approach, whereby the design variable is a density function ρ : D → (0, 1) such that ρ(x) corresponds to the local volume fraction of material of conductivity κ f in an infinitesimal cell centered around x. The conductivity κ(ρ) inside the whole domain D is determined by the SIMP interpolation law with penalization exponent p = 3:

Q. A prescribed temperature T = 0 is applied on a part Γ D of
κ(ρ) = ρ p (κ f -κ s ) + κ s . (24)
Of course, the purpose of the exponent p is to penalize intermediate densities so that numerically optimized designs ideally satisfy ρ(x) = 0 or ρ(x) = 1 in the whole domain.

The mathematical formulation of the heat sink design problem reads

min ρ J(ρ) = D T dx s.t.    V (ρ)/V 0 := 1 |D| D ρ dx ≤ volfrac 0 ≤ ρ(x) ≤ 1, for all x ∈ D; (25)
where T ≡ T (ρ) is the solution to the heat equation

       -div(κ(ρ)∇T ) = Q in D, T = 0 on Γ D , ∂T ∂n = 0 on Γ N . (26)
Note that in this test case, we treat the volume constraint as an inequality constraint.

In the following numerical applications, the numerical values of the physical parameters are set to

κ f = 401, κ s = 1, Q = 10 4 , volfrac = 0.1.
Equation ( 26) is solved with the finite element method by assembling the weak formulation [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]).

find T ∈ V such that ∀S ∈ V, D κ(ρ)∇T • ∇S dx = D QS dx, (27) where V = {v ∈ H 1 (D) | v |Γ D = 0} is the subspace of H 1 functions vanishing on Γ D

Following

the "optimize-then-discretize" approach, we discretize the Fréchet derivatives of J(ρ) and V (ρ) rather than computing the derivatives of the discretized finite element functionals. For this classical self-adjoint problem, the Fréchet derivatives (Allaire, 2007, Definition 10.1.1) of the objective function and volume constraint respectively read:

DJ(ρ) • δρ = D DT (ρ) • δρ dx = 1 Q D κ(ρ)∇T • ∇(DT (ρ) • δρ) dx = - 1 Q D [κ ′ (ρ)∇T • ∇T ]δρ dx, (28) D[V (ρ)/V 0 ] • δρ = 1 |D| D δρ dx, (29)
where ( 28) is obtained by differentiating the variational equation ( 27) with respect to ρ.

Finally, in order to avoid checkerboard effects on the unstructured discretization mesh, we actually solve the following filtered version of (25):

min ρ J( ρ) = D T ( ρ) dx s.t. V ( ρ)/V 0 ≤ volfrac, 0 ≤ ρ(x) ≤ 1, for all x ∈ D, (30)
where ρ is a regularization of the design variable ρ obtained by solving the variational problem

Find ρ ∈ H 1 (D) such that ∀v ∈ H 1 (D), D (γ 2 ∇ ρ • ∇v + ρv) dx = D ρv dx, (31)
with γ a smoothing parameter set to the minimum edge length. The variational problem [START_REF][END_REF] is the weak formulation of the equation

-γ 2 ∆ ρ + ρ = ρ
with Neumann boundary conditions on ∂D. For this reason, the mapping ρ → ρ is sometimes called "Helmholtz type filter" despite ( 31) is not a wave equation. The density filter ( 31) is often preferred over (18) due to its ease of implementation on unstructured meshes. We refer the reader to [START_REF] Lazarov | Filters in topology optimization based on Helmholtz-type differential equations[END_REF] for an extensive comparison of filters in density based Topology Optimization.

Python implementation with the Null Space Optimizer and FreeFEM

The complete implementation of the heat Topology Optimization test case is provided in Appendix B. The code is divided into a main Python file (Appendix B.1) and two FreeFEM scripts filter matrices.edp and solve state.edp (respectively Appendix B.2 and Appendix B.3). The first one performs finite element operations for the filtering of the design variable, and the second one solves the state equations before computing the objective and constraint functions and their sensitivities. These two scripts are written in the FreeFEM language [START_REF] Hecht | New development in FreeFem++[END_REF]) for the convenience of the implementation, and they are interfaced with the main Python script through the library PyFreeFEM 5 importing the class FreeFemRunner at line 2 of the main script. In what follows, we systematically refer to the line numbering of the main script of Appendix B.1 unless otherwise specified.

Declaration of the optimization problem

Similarly to the MBB beam test case, the code is best understood by looking first at the declaration of the class Heat TO at lines 61-95. The implementation carefully follows the mathematical definition of the design optimization problem (25): the class implements methods x0, J, H returning respectively the initial design, the values of the objective function and of the volume inequality constraint.

Similarly to the implementation of the MBB test case of Section 3, the method x0 returns a constant array of length the number of mesh triangles corresponding to the initial design density ρ(x) = volfrac. In the present implementation, the design variable is called rho in the Python script of Appendix B.1; it is a P0 discretization of the continuous variable ρ, that is constant on every triangle of the discretization mesh.

The implementation of the optimization class Heat TO at lines 61-95 requires several ingredients which are described in the next subsections:

• a computational mesh Th with Th.nt triangles (Section 5.2.2); • a state solver solve state that returns a dictionary with the objective and constraint values as well as their sensitivities (Section 5.2.3); • the declaration of an appropriate filter function filter as well as its derivative diff filter on the triangular mesh Th (Section 5.2.4).

Construction of design domain mesh

The design domain D is discretized with a triangular finite element mesh Th that is generated at lines 15-19. The mesh build in FreeFEM with (line 17)

5 https://gitlab.com/florian.feppon/pyfreefem where $N is a prescribed number of edges on each side and flags=1 allows to obtain a symmetric mesh. The FreeFEM code is executed at line ( 19) which returns a Python variable Th of the Mesh type defined the library PyMedit6 (line 20). The Mesh type complies with the INRIA mesh format and could be easily exported to or imported from a .mesh or .meshb file. The variable Th is then manually processed to change the boundary labels, so that Γ D is assigned to the label 5 and Γ N to the label 1 (lines 21-30).

The mesh obtained for N=10 is represented on Fig. 13. It would be straightforward to run the test case on any other triangular mesh by updating the variable Th.

Implementation of the physical state solver

The The method import variables of the library PyFreeFEM automatically stores the Python array x and the mesh Th in a temporary directory, which are then made accessible in the FreeFEM script solve state.edp of Appendix B.3 through the instructions importMesh("Th") and importArray("rho") (lines 3 and 9 of this script). The script solve the heat equation (lines [START_REF]solve heat (T , S )= 22 int2d ( Th )( kappa * grad ( T ) '* grad ( S )) 23 -int2d ( Th )( Q * S ) on (5 , T =0)[END_REF][22][23][24] and compute the objective function (line 26), the volume constraint (line 29) and their sensitivities (lines 32-33 and 35-36). As these lines demonstrate, the use of FreeFEM is especially convenient for implementing ( 25) and ( 27) to [START_REF]* volfrac ))[END_REF] with a language very close to the continuous mathematical formulation.

The objective function, volume constraints and their sensitivities are then exported towards the Python script at lines 38-42 with the instructions exportVar and exportArray. Then, the execution of the FreeFEM script occurs at line 54 of the main script of Appendix B.1 with the method execute which returns a Python dictionary exports gathering the exported data structures: for instance exports['J'] and exports ['dJ'] are respectively a floating number and a NumPy array containing the objective function value and its sensitivity.

Implementation of a smoothing density filter

The implementation of the density filter ρ → ρ defined by ( 31) is automated by the decorator @ f i l t e r e d _ o p t i m i z a b l e ( filter , diff_filter )

at line 63 of the main Python script. Here, the implementation of the filtering procedure requires a little caution because ρ is implemented as a P0 function rho ≡ ρ P0 while the variational problem [START_REF][END_REF] solved with P1 finite elements returns a P1 function ρ P1 . To address this issue, the proposed implementation converts this P1 function ρ P1 to a P0 function ρ P0 by solving the variational problem: find a P0 function ρ P0 such that for any P0 function v P0 ,

D ρ P0 v P0 dx = D ρ P1 v P0 dx. ( 32
)
In a shorthand matrix notation, the numerical filtering procedure mapping P0 to P0 densities reads

ρ P0 := M -1 P0P0 M T P1P0 A -1 M P1P0 ρ P0 . ( 33
)
where we have denoted by M P1P0 and M P0P0 the mass matrices discretizing the bilinear form (u, v) → D uv dx with v a P0 function and u respectively a P1 or a P0 function.

The assembly of the finite element matrices A, M P 0P 0 and M P 1P 0 is performed by the FreeFEM script filter matrices.edp of Appendix B.2. This script is executed between lines 32-37 of the main Python script of Appendix B.1. The matrices are exported to sparse Scipy matrices and a factorization of A and M P 0P 0 is computed at lines 37-38 in order to compute efficiently the mapping (33).

Finally, the filtering map filter and its Fréchet derivative diff filter that are required by the decorator @filtered optimizable are implemented at lines 41-47 of the main script by copying the formula (33).

Running the optimizer

The test case is instantiated at line 100 with the function init(100) which constructs a computational mesh Th with N=100 edges on each boundary. The optimization parameters are set at lines 100-104, setting a time step dt=0.3, 50 normalization iterations and stopping the optimization after 150 iterations. The design optimization problem is solved with the Null Space Optimizer at line 105 by calling the function nlspace solve. The execution of the code generates plots of the intermediate densities thanks to the instructions given in the method accept at lines 87-95, relying on the plotting features of PyMedit (line 92).

Numerical results and comparisons with OC, MMA and IPOPT

We solve the heat Topology Optimization problem (30) with the Null Space Optimizer, MMA, IPOPT and the Optimality Criteria method. The move parameters of the OC method and MMA are set to 0.1. IPOPT runs with the default settings.

The numerical values of the optimized designs found by the different optimizers are reported in For this example, the OC method converges to a design which performs more poorly (by a factor of approximately 30%) than its competitors. MMA, NLSPACE and IPOPT converge to different designs with very similar performance. Notably, the four optimizers drive the design variable on significantly different optimization paths.

Algorithm J H = V (x)/V 0 -
Let us conclude by noticing that none of the optimizers converge to needle like structures that are known to be optimal for this problem [START_REF] Yan | On the nonoptimality of tree structures for heat conduction[END_REF]. Appendix A Source code of the MBB beam example of Section 3 . For many practical situations including density based Topology Optimization, the reader may assume that A is the identity matrix. However, we need to consider an arbitrary inner product A to account for more general contexts such as the Hilbertian setting of level-set based Topology Optimization where such matrix is used to regularize shape derivatives [START_REF] De Gournay | Velocity Extension for the Level-set Method and Multiple Eigenvalues in Shape Optimization[END_REF]).

[ k [0] , k [1] , k [2] , k [3] , k [4] , k [5] , k [6] , k [7]] , 19 [ k [1] , k [0] , k [7] , k [6] , k [5] , k [4] , k [3] , k [2]] , 20 [ k [2] , k [7] , k [0] , k [5] , k [6] , k [3] , k [4] , k [1]] , 21 [ k [3] , k [6] , k [5] , k [0] , k [7] , k [2] , k [1] , k [4]] , 22 [ k [4] , k [5] , k [6] , k [7] , k [0] , k [1] , k [2] , k [3]] , 23 [ k [5] , k [4] , k [3] , k [2] , k [1] , k [0] , k [7] , k [6]] , 24 [ k [6] , k [3] , k [4] , k [1] , k [2] , k [7] , k [0] , k [5]] , 25 [ k [7] , k [2] , k [1] , k [4] , k [3] , k [6] , k [5] , k [0]]]) 26
Using an arbitrary inner product has for consequence that the transpose operator • T must be preceded by A -1 in the definitions (7), ( 8) and ( 15) of respectively ξ J (x), Λ * (x) := (λ * (x), µ * (x)) and ξ C (x). The null space direction ξ J (x) now reads

ξ J (x) = A -1 DJ(x) T + A -1 DC I(x) (x) T Λ * (x), ( C1
)
where the multiplier Λ * (x) := (λ * (x), µ * (x)) is the minimizer of the problem min

Λ=(λ,µ)∈R p ×R q(x) + ||A -1 DJ(x) T +A -1 DC I(x) (x) T Λ|| A .
(C2) The range space direction is updated as follows:

ξ C (x) = A -1 DC T I(x) (DC I(x) A -1 DC T I(x) ) -1 C I(x) .
(C3) The reader may verify that these formulas coincide with ( 7), ( 8) and (10) when substituting A with the n × n identity matrix.

We now show that it is possible to rewrite (C1) to (C3) only in terms of DC I(x) and A, avoiding altogether the need for assembling DC I(x) A -1 DC I(x) or computing its explicit inverse. Proposition 1. (i) The null space direction ξ J (x) defined in (C1) is equivalently given by ) where X * (x) is obtained from the solution (X * (x), Λ * (x)) to the quadratic minimization program min

ξ J (x) = A -1 DJ(x) T + X * (x), ( C4
(X,Λ)∈R n ×R p+ q(x) 1 2 X T AX + DJ(x)X s.t.
AX -DC T I(x) Λ = 0, Λ = (λ, µ) ∈ R p × R q(x) with µ ≥ 0.

(C5) (ii) The minimizer Λ * (x) of (C5) is also the optimal Lagrange multiplier Λ * (x) = (λ * (x), µ * (x)) as defined by (C2). (iii) The range space direction ξ C (x) defined in (C3) can be obtained from the solution to the symmetric indefinite linear system

A DC T I(x) DC I(x) 0 ξ C (x) Λ = 0 C I(x)
. (C6)

Proof. The square of the right-hand side of (C2) reads

||A -1 DJ(x) T + A -1 DC T I(x) Λ|| 2 A = DJ(x)A -1 DJ(x) T + 2DJ(x)A -1 DC T I(x) Λ + Λ T DC I(x) A -1 DC T I(x) Λ = X T AX + DJ(x) T X, where we have introduced X := A -1 DC T I(x) Λ. Since X is equivalently characterized by AX -DC T I(x) Λ = 0, this proves the point (ii) and the point (i) follows from (C1). The point (iii) is obtained by observing that DC I(x) ξ C (x) = C I(x) and Aξ C (x) + DC T I(x) Λ = 0 where Λ := -(DC I(x) A -1 DC T I(x) ) -1 C I(x) . When both the matrices A and DC I(x) are sparse, Proposition 1 provides a definition of ξ J (x) and ξ C (x) that becomes accessible to iterative solvers. A number of strategies can be devised for efficiently solving the symmetric indefinite system (C6), also called KKT system, see e.g. [START_REF] Freund | A new Krylovsubspace method for symmetric indefinite linear systems[END_REF]; [START_REF] Ashcraft | Accurate symmetric indefinite linear equation solvers[END_REF]; [START_REF] Toh | Block preconditioners for symmetric indefinite linear systems[END_REF]; [START_REF] Bai | Constraint preconditioners for symmetric indefinite matrices[END_REF]; [START_REF] Nocedal | Numerical optimization, 2nd edn[END_REF]. The present implementation of the Null Space Optimizer solves (C6) with Scipy lsqr method with machine precision. The quadratic minimization problem (C5) with the nonnegativity constraint µ ≥ 0 is solved with OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]) or CVXOPT [START_REF] Vandenberghe | The CVXOPT linear and quadratic cone program solvers[END_REF]), which both internally rely on sparse linear solvers. Finally, let us note that the proposed implementation also addresses the issue of degenerate constraints, which are either dealt with the internal iterative solvers of OSQP or CVXOPT, or with Scipy least-squares lsqr solver. All the examples of Sections 3 to 5 considered in this paper relied on the OSQP solver as it seems to provide in general better runtime performances than CVXOPT (see Appendix C.2).

Remark 5. The quadratic program (C5) is rather large if the matrix DC I(x) features a large number of constraints, resulting in a potentially higher cost per optimization iteration than solving the MMA subproblem [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF]). However, for the test cases considered in this article, we found the OSQP solver quite efficient at solving the quadratic problem (C5) in reasonable time, which is therefore not a bottleneck of the method.

C.2 A benchmark example: generating a checkerboard

In order to illustrate the numerical treatment of bound constraints and the influence of the quadratic programming solver, we consider the benchmark minimization problem of generating a checkerboard on a n × n grid: min

(xij ) 1≤i,j≤n J((x ij ) 1≤i,j≤n ) s.t. 0 ≤ x ij ≤ 1 for all 1 ≤ i ≤ n. (C7)
where the function J((x ij ) 1≤i,j≤n ) measures the opposite of the discrepancies of the variables (x ij ) 1≤i,j≤n between adjacent cells on the n × n grid:

J((x ij ) 1≤i,j≤N ) := - 1 2 n-1 i=1 n j=1 (x i+1,j -x i,j ) 2 + n i=1 n-1 j=1 (x i,j+1 -x i,j ) 2 (C8)
The global minimum is clearly attained for a checkerboard, namely for a variable (x ij ) 1≤i,j≤n satisfying |x i+1,j -x ij | = |x i,j+1 -x ij | = 1 for all 1 ≤ i, j ≤ n. However, due to the quadratic exponents, the problem (C7) admits many local minima. The minimization is thus very sensitive to numerical noise and to the chosen optimization path, which makes it an interesting benchmark for comparing several optimizers.

We run the Null Space Optimizer with the time step dt = 0.3, using either OSQP or CVX-OPT as the quadratic programming solver for (3). The number of normalization iterations for the null space step is set to itnormalisation=10. The tolerance parameters of both the OSQP and CVXOPT solver are set to 10 -15 . The number of scaling iterations as well as the maximum number of iterations allowed for the OSQP solver is set to 400. The maximum number of iterations for the CVXOPT solver is set to 30 (these are observed to be computationally more costly).

The initialization is defined to be the density x 0 ij := ((i-1)/(n-1)-0.5)((j -1)/(n-1)-0.5), 1 ≤ i, j ≤ n, (C9)

which is symmetric with respect to the grid medial axes. We also solve (C7) with the optimizers IPOPT and MMA for comparison (we do not use the OC method which is not relevant in this context devoid of equality constraints). The move parameter of MMA is set to 0.3 and IPOPT is run with the default parameters.

The numerical values of the computed optima are listed in Table C1 and the convergence histories are shown on Fig. C1. Intermediate design variables are shown on Fig. C2 for the various optimization strategies. On this test case, IPOPT finds the best optimum. The version of MMA provided by [START_REF] Deetman | GCMMA-MMA-Python[END_REF] converges to a very poor solution, probably due to the fact that the default moving asymptotes rules is not suited on this particular example that only features bound constraints.

The CVXOPT variant of the Null Space Optimizer converges to a rather suboptimal local minimum, although it provides a slightly faster decay than its OSQP variant at the first iterations: CVXOPT seems unable to update the design variables in the middle black region (Fig. C2b). The OSQP variant should in theory find the same numerical result than the with the CVXOPT one since both solvers solve the same dual problem (3). The difference between the two can be attributed to the fact that the OSQP solver seems less conservative on the nonnegativity constraint µ ≥ 0 of (C5) than CVXOPT, accepting negligible violations. This seems to be advantageous to escape the local minimizer and to help finding a design that has a performance rather comparable to the one found by IPOPT in a similar number of iterations.

The significant discrepancy between the numerically computed designs illustrates the strong dependency on the optimization path chosen by the optimizer and on numerical round off errors. This also suggests some limitations to the standard implementation of the MMA algorithm, which has difficulty to converge on a case devoid of constraints other than the bound constraints on the design variable. It would be tempting to conclude that IPOPT performs better than the Null Space Optimizer, and that the Method of Moving Asymptotes is inappropriate, however, we have seen in the tutorial sections that this is of course very situation specific. Once again, we emphasize that the performance of an algorithm to find an optimum as close as possible to the global optimum is situation dependent and cannot constitute a single selection criterion between several optimizers for their use on more general situations.

Algorithm

J

Iteration number at convergence NLSPACE (OSQP) -6.27600 • 10 3 65 NLSPACE (CVXOPT) -3.70000 • 10 3 65 MMA -1.07864 • 10 0 30 IPOPT -6.88800 • 10 3 58 Table C1: Final objective values for the checkerboard example of Appendix C.2. Comparison between the Null Space Optimizer using either CVXOPT or OSQP for solving the quadratic subproblem (3), MMA and IPOPT. asymptotes. This solution came from discussions between N. Aage, Michael Stingl and Fabian Wein. I also thank very much R. Krasniqi for his Master thesis work [START_REF] Krasniqi | Density-based topology optimization with the Null Space optimiser[END_REF]) upon which part this work is built.

Fig. 1 :

 1 Fig. 1: Setting of the MBB beam Topology Optimization problem of Section 3.

  Volume constraint G(x) = V (x)/V 0volfrac.

Fig. 2 :

 2 Fig. 2: Convergence history of the different optimizers on the MBB beam example of Section 3.4.1.

  Volume constraint G(x) = V (x)/V 0volfrac.

Fig. 3 :

 3 Fig. 3: Convergence history of the different optimizers on the MBB beam example of Section 3.4.1 (zoom on the first 30 iterations).

Fig. 4 :

 4 Fig. 4: Comparison of the designs obtained by the different optimizers for the MBB beam test case of Section 3.4.1 at iterations 0, 10, 15, 20, 30, 80 and 150.

.

  In our "fixed" implementation of MMA, we manually disable this rule and we leave the remainder of the code provided by Deetman (2020) unchanged.• using tighter moves between iterations. Our default implementation of MMA sets x min e = 0 and x max e = 1 for all 1 ≤ e ≤ n. The fixed implementation forbids updates larger than the parameter move by settingx min e = max(0, x (k)move),x max e = min(1, x (k) + move). at every iteration.

Fig. 5 :

 5 Fig. 5: Convergence history of the standard implementation of MMA on the MBB beam example of Section 3.4.2 with volume fraction target volfrac=0.01 (iterations 0, 20, 30, 80, 100, 150).

  constraint G(x) = V (x)/V 0volfrac.

Fig. 6 :

 6 Fig. 6: Convergence history of the different optimizers on the MBB beam example of Section 3.4.2 with tiny volume fraction target.

  constraint G(x) = V (x)/V 0volfrac.

Fig. 7 :

 7 Fig. 7: Convergence history of the different optimizers on the MBB beam example of Section 3.4.2 with tiny volume fraction target (close-up on the first 30 iterations).

Fig. 8 :

 8 Fig. 8: Intermediate designs obtained with the different optimizers at iterations 0, 10, 15, 20, 40 and 80 for the MBB beam example of Section 3.4.2 with tiny volume fraction target.

@

  mi n m ax _ o pt i m i za b l e @ b o u n d _ c o n s t r a i n t s _ o p t i m i z a b l e ( l =0 , u =1) @ f i l t e r e d _ o p t i m i z a b l e ( filter , diff_filter ) class TO_problem ( E u c l i d e a n O p t i m i z a b l e ): # ... # identical to lines 152 -185

Fig. 10 .

 10 Intermediate designs obtained during the iterations of the three algorithms are displayed on

Fig. 10 :

 10 Fig. 10: Convergence history of the various optimizers for the multiple load bridge test case of Section 4.

  • 10 -1 9.65891 • 10 -1 1.00830 • 10 0 G = V (x)/V 0volfrac -3.10 • 10 -6 -1.60 • 10 -6 -9.80 • 10 -7 Table3: Final objective and constraint values obtained for the multiple load case of Section 4 with the Null Space Optimizer, MMA and IPOPT.

Fig. 11 :

 11 Fig. 11: Intermediate designs obtained by the different optimizers at iterations 0, 10, 40, 80, and 150 for the bridge test case of Section 4.

DFig. 12 :

 12 Fig.12: Geometry for the heat sink test case of Section 5.

Fig. 13 :

 13 Fig. 13: Finite element mesh for the heat Topology Optimization problem of Section 5 with boundary labels corresponding to Γ D and Γ N .

  )/(V 0 volfrac) -1.

Fig. 14 :

 14 Fig. 14: Convergence history of the different optimizers tested on the heat Topology Optimization test case of Section 5.

Fig. 15 :

 15 Fig. 15: Intermediate designs of the heat Topology Optimization test case of Section 5 obtained with the different optimizers at iterations 0, 10, 15, 20, 40 and 80.

  Fig. C1: Convergence history for the checkerboard example of Appendix C.2 for the three optimizers tested.

  (a) Null Space (OSQP), iterations 0, 22, 45 and 68. (b) Null Space (CVXOPT), iterations 0, 24, 49 and 74. (c) MMA, iterations 0, 23, 46 and 70. (d) IPOPT, iterations 0, 20, 41 and 62.

Fig. C2 :

 C2 Fig. C2: Intermediate optimized designs for the various optimizers tested in the checkerboard example of Appendix C.2.

Table 1 :

 1 Final objective and constraint values for the MBB example for the different optimization algorithms.

			volfrac
	NLSPACE 288.702	0
	OC	288.830	-1.30 • 10 -5
	MMA	286.383	-2.20 • 10 -6
	IPOPT	288.801	-9.30 • 10 -8

Table 2 :

 2 Final objective and constraint values for the MBB example of Section 3.4.2 with tiny volume fraction target for the different optimization algorithms.

Table 4 :

 4 Final objective and constraint values for the heat Topology Optimization test case of Section 5 obtained for the tested optimization algorithms.

			volfrac
	NLSPACE 2.45701 • 10 2	9.60 • 10 -11
	OC	3.60279 • 10 2	-7.90 • 10 -6
	MMA	2.46293 • 10 2	-1.90 • 10 -5
	IPOPT	2.45544 • 10 2	-1.90 • 10 -6

Table 4 .

 4 Convergence histories of the objective function and volume constraint are reported on Fig.14. Figure15depicts a sequence of intermediate designs obtained for each of the optimizers.

https://github.com/arjendeetman/GCMMA-MMA-Python

3 https://github.com/mechmotum/cyipopt h := (h j ) 1≤j≤q are respectively the objective function, p equality constraints and q inequality constraints. For density based Topology Optimization, the design set is the finite-dimensional space X = R n , where n is the number of cells or mesh elements discretizing the design domain.

Note that H(self,x) is a method of TO problem which is not to be confused with the filtering matrix also denoted by H (a global variable). We use H for the filtering matrix to match the notation of[START_REF] Andreassen | Efficient topology optimization in MAT-LAB using 88 lines of code[END_REF].

https://gitlab.com/florian.feppon/pymedit

B.2 File filter matrices.edp 1 IMPORT " io . edp " 2 mesh Th = importMesh ( " Th " ); 3 4 fespace Fh0 ( Th , P0 ); 5 fespace Fh1 ( Th , P1 ); 6 7 varf mass (u , v )= int2d ( Th )( u * v ); 8 9 matrix MP1P0 = mass ( Fh0 , Fh1 ); 10 matrix MP0P0 = mass ( Fh0 , Fh0 ); 11 12 macro grad ( u ) [ dx ( u ) , dy ( u )] // 13 14 real gamma = Th . hmin ; 15 varf helmholtz (u , v ) = int2d ( Th )( gamma ^2* grad ( u ) '* grad ( v )+ u * v ); 16 matrix A = helmholtz ( Fh1 , Fh1 ); 17 18 exportMatrix ( MP1P0 ); 19 exportMatrix ( MP0P0 ); 20 exportMatrix ( A ); B.3 File solve state.edp 1 IMPORT " io . edp " 2 3 mesh Th = importMesh ( " Th " ); 4 5 fespace Fh0 ( Th , P0 ); 6 fespace Fh1 ( Th , P1 );

Appendix C Treatment of constraints with sparse Jacobian matrix

In this appendix, we describe the modifications of the Null Space Optimizer that enable to solve optimization programs with constraints with sparse Jacobian matrix, deviating from the original description of this algorithm in Feppon et al (2020a). The theoretical changes for the computation of the null space and range space directions are detailed in Appendix C.1. Then, in Appendix C.2, we challenge the Null Space Optimizer on a benchmark numerical example featuring many bound constraints and we provide qualitative comparisons with the standard IPOPT and MMA optimizers.

Appendix D Statements and declarations

• The author has no relevant financial or nonfinancial interests to disclose.

Appendix E Reproducibility of numerical results

The integral source code corresponding to the examples presented in Sections 3 to 5 are included in the examples/topopt examples directory of the GitLab repository of the Null Space Optimizer. Since the source code published on this repository is subject to future updates, a copy of the version used for the present manuscript is available as a supplementary material to this article.

Appendix F Acknowledgements

I sincerely thank Niels Aage for interesting discussions about the MMA algorithm, for suggesting the test case of Section 3.4.2 and sharing the solution of using the fixed initialization for the