Comparison of heuristics and metaheuristics for topology optimisation in acoustic porous materials

When designing sound packages, often fully ﬁlling the available space with acoustic materials is not the most absorbing solution. Better solutions can be obtained by creating cavities of air pockets, but determining the most optimal shape and topology that maximises sound absorption is a computationally challenging task. Many recent topology optimisation applications in acoustics use heuristic methods such as solid-isotropic-material-with-penalisation (SIMP) to quickly ﬁnd near-optimal solutions. This study investigates seven heuristic and meta-heuristic optimisation approaches including SIMP applied to topology optimisation of acoustic porous materials for absorption maximisation. The approaches tested are hill climbing, constructive heuristics, SIMP, genetic algorithm, tabu search, covariance-matrix-adaptation evolution strategy (CMA-ES), and diﬀerential evolution. All the algorithms are tested on seven benchmark problems varying in material properties, target frequencies, and dimensions. The empirical results show that hill climbing, constructive heuristics, and a discrete variant of CMA-ES outperform the other algorithms in terms of the average quality of solutions over the diﬀerent problem instances. Though gradient-based SIMP algorithms converge to local optima in some problem instances, they are computationally more eﬃcient. One of the general lessons is that diﬀerent strategies explore diﬀerent regions of the search space producing unique sets of solutions.


I. INTRODUCTION
A. Background Historically, shape designs in engineering have been arrived at via a trial-and-error process, intuition, incremental improvements to old designs, human decisionmaking from numerical analyses, and recently, solely by computer analyses.Superior-to-human engineering designs have been achieved by computers using technologies such as structural topology optimisation.Topology optimisation involves finding the optimal topology (number of holes) and shape (size, dimensions) for a structure such that a given performance indicator is either maximised or minimised.Bendsøe and Kikuchi 1 introduced the concept of simultaneously optimising both shape and topology in the late 1980s.Since then, many theoretical developments have been made, and a community of researchers have actively been working in this field.One of the ways to formulate a topology optimisation problem is finding the optimal assignment of materials in each finite a) vivek.thaminniramamoorthy@nottingham.ac.uk element of a discretised structure.In principle, this for-  sation problems are SIMP 1,4-6 (solid-isotropic-material-with-penalisation), BESO [7][8][9] (bi-directional evolutionary structural optimisation), and the level-set method 10-12 .Among these, SIMP is the most commonly used and well-studied approach.In SIMP, the discrete material assignment problem is relaxed to the continuous space by allowing intermediate materials between solid and void.
A penalty-based material interpolation scheme is used to represent intermediate materials and gradient-based optimisation strategies such as optimality criteria 13 or method of moving asymptotes 14 is used to move across the design variable space to find a near-optimal design.
As SIMP is a derivative-based technique, it requires that a sensitivity analysis be carried out.BESO, not to be confused with evolutionary algorithms despite its name, is a type of constructive approach which iteratively adds material where stresses are high and removes material where stresses are low to arrive at a design.In the levelset method, a scalar field is associated with the design domain region and the isosurfaces of this scalar field are made the boundaries of the topology.This scalar field is then optimised to optimise the topology.

C. Metaheuristics
While heuristics are quick strategies to find nearoptimal solutions, it was realised by Glover 15 that many powerful heuristic approaches follow certain higher-level guidelines.These guidelines can be considered heuristics to design heuristic algorithms, and hence are termed metaheuristics.A popular example of a metaheuristic is genetic algorithms, wherein the guideline is to initiate a population of solutions, apply selection pressure to pick good individuals, recombine the selected individuals, mutate them and replace them into the population.Numerous metaheuristic techniques, such as genetic algorithms and CMA-ES, have also been studied on structural topology optimisation problems 16,17 .

D. Acoustic topology optimisation
Theoretical developments in structural topology optimisation have focused on the classical problem of compliance minimisation 18,19 .Nevertheless, the application of topology optimisation techniques to other problem domains is steadily on the rise 18,20,21 .These techniques have already been extended to acoustics, giving rise to a sub-field called acoustic topology optimisation.
At the time of writing this article, topology optimisation has been performed on a variety of acoustic applications, including horns, mufflers, rooms and sound barriers [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] .A majority of these applications use the gradient-based SIMP method or its variants, while a small fraction of them use BESO or level-set methods.These applications can be categorised into acoustic fluid-structure interaction problems and porous material problems.In acoustic fluid-structure interaction problems, the material choices are non-porous solid and fluid phases, and the wave propagation is modelled using mixed formulations 37,38 .Within acoustic fluid-structure interaction, problems other than topology 96 optimisation such as material parameter estimation 39 97 have also found the application of gradient-based meth-98 ods such as the method of moving asymptotes 40 .In 1. Finite element model of an impedance tube system with the design domain where the shape and topology of a poroelastic material is to be optimised.153

II. PROBLEM DESCRIPTION AND MODELLING
The optimisation formulation can be written as: where α(χ, f ) is the sound absorption coefficient in normal incidence for a given shape χ for frequency f , χ i are the decision variables represent the choice between air and porous material for the i th element, N is the number of elements in the design domain, and f 1 , f 2 , ..., f n are the target frequencies for which the mean absorption is to be maximised (where n is the number of frequencies considered).The symbol α is used to refer to the mean sound absorption coefficient (α) across the target frequencies.
In this paper, α may be referred to as simply absorption or fitness, which is to be maximised.Note that in the problem formulation1, a volume fraction constraint is not included, which is unlike in usual topology optimisation problems.One reason is that in porous material topology optimisation, often the optimal shapes need to be carved out from a large block of the base porous material.The removed material may not often constitute materialsaving, as the cost of recycling the carved out material could negate the material-saving benefit.Another reason is that more optimisation approaches to be tested as the formulation would resemble a conventional discrete optimisation problem.Without the volume constraint, since two choices are available (air or the base porous material) for each of the N elements in the design domain the search space size becomes 2 N .If a limit V f is imposed on the ratio of porous volume to the total volume in the design domain ( 1 N N i=1 χ i = V f ), the search space size would become N C (V f N ) .In both these cases, the number of feasible solutions grows quickly with an increase in N .Since discrete optimisation is considered difficult to solve, the problem is usually relaxed to a continuous problem allowing χ i to take values between 0 and 1, in other words allowing intermediate materials between air and porous material in the design domain.The problem is then solved using continuous optimisation approaches.Intermediate materials given by χ i ∈ (0, 1) are modelled using interpolating the material properties.One such interpolation scheme is the SIMP scheme (not to be confused with the SIMP approach).Using this scheme, the material property ψ for the intermediate material is given by equation 2.
Here, ψ is any material property from Young's modu- where (•) denotes the complex-valued nature of it's argument.The expressions for the state matrices K, M, H, coefficient α is then given by: The analytical gradient of absorption can be com-208 puted by using chain rule following a similar procedure 209 to that of Lee et al. 24 .From equation 7: Equation 9 computes the derivative of absolute of 217 The above step involves a large matrix inversion followed by sparse matrix and vector multiplications repeated for each element in the design domain.This step is performed efficiently by using the adjoint-based approach as detailed by Lee, Göransson and Kim 28 .Since only two elements in ∂ X ∂χi i.e. ∂χi need to be computed to compute the gradients, one can premultiply the equation 10 by the term ∂Px 1 ∂X , which is a vector of 0s except for one element with a value of 1 corresponding to the P x1 degree of freedom in equation 10.  along the horizontal and vertical directions respectively.
Within the unit cell, symmetry is assumed about the central horizontal line, and sliding boundaries (u x -free, u y = 0, P -free) are assumed at the top and bottom edges.
To save computational effort, only half of the system is modelled, and symmetry is imposed about the centerline (u x -free, u y = 0, P -free).It has been verified that this gives the same absorptions as obtained when modelling the full unit cell with sliding supports in the top and bottom edges.In all the problem instances, the mean sound   Often in practice, a particular topology optimisation strategy may be chosen, and one trial may be run to determine a near-optimal shape.In such cases, it is desirable to pick an algorithm that has the best median performance across trials.Hence, using the median absorption across trials, the algorithms are sorted best to worst from left to right in Figure 2

C. Performance across problem instances
For an overall comparison, the ranking is extended to other problem instances in Table IV.Such a comparison across many problem instances is essential as algorithms performing well on one problem instance need not necessarily perform well on other problem instances.The ranking scheme is such that if the median absorption values of two or more algorithms are the same correct to two decimal places, they are assigned the same rank.Although the progress of absorption in the initial stages of CH1 is slow compared to the other algorithms, the final absorption value makes CH1 one of the best algorithms.Notably, for many problem instances considered, the best absorption value from CH1 is higher than the absorption of the discretised solutions from both SIMPf0 and SIMPf2.CH1 seems to be better overall compared to CH2, indicating that constructing the solution from scratch may be better than removing material from a fully-filled solution.
Performance of CMA and DE were relatively poor in this benchmark.One reason could be that the number of design variables is large and these strategies do not exploit the correlation of the neighbouring-element design variables, a special attribute in topology optimisation problems.
Both CMAd and DEd seem to perform better than CMA and DE in general, indicating that rounding during the algorithm may be a better approach than rounding the solutions after the termination of continuous algorithms.While CMAd ranked among the top, the performance of DEd was similar to that of SIMP in terms of solutions quality.
Among the algorithms considered, GA performed the poorest.Though, scope for improvement exists in terms of using better mutation and crossover operators adapted to topology optimisation, focus may be diverted to other strategies which show better promise.

D. Best shapes obtained from algorithms
The best solutions from all the algorithms for all problem instances are plotted in In problem instance 3 with a high static airflow resistivity material, the shapes from all algorithms were seemingly random patterns but with sort of a cavity in the centre.SIMPf2 produces a result with a chunk of porous material suspended in the air.
For problem instance 4, the optimal solution seems to be a fully-filled design domain and most algorithms are able to find this except for GA.

20
mulation is discrete optimisation, and finding the exact 21 global optimum is computationally challenging.Exact 22 optimisation techniques that guarantee finding the global 23 optimum remain prohibitively expensive.Evaluating all 24 possible solutions becomes impractical due to the large 25 search space sizes and the expensive finite element eval-26 uations.A noteworthy effort towards topology optimisa-27 tion using an exact approach was by Stolpe and Bendsøe 2 28 on the Zhou and Rozvany problem instance 3 .But justi-29 fiably, the focus of previous work has mainly been on the 30 inexact or heuristic optimisation approaches.31 B. Heuristics 32 Heuristics are techniques that find solutions close 33 enough to the global optimum in a reasonable time.34 Though heuristics do not guarantee finding the optimal 35 solution, they are well-established and often the only 36 viable option to address hard problems, such as those 37 in NP-complete and NP-hard classes.The three most 38 popular heuristic approaches applied to topology optimi-

99
FIG.1.Finite element model of an impedance tube system with the design domain where the shape and topology of a poroelastic material is to be optimised.

147A.
Maximising sound absorption in normal incidence 148 Consider the following problem: Given a finite ele-149 ment model of an impedance tube as shown in Figure 1, 150 what is the best assignment of either air or a given poroe-151 lastic material to each element in the design domain that 152 maximises the sound absorption of an acoustic source.

183Q
and C are functions of the topological design/decision 184 variables χ .The construction of these matrices are ex-185 plained by Atalla et al. 41 and will not be detailed here.186 {u} and {P} denote the solid phase displacement and 187 fluid phase pressure degrees of freedom in the poroelastic 188 system respectively.The associated global stiffness ma-189 trix S(ω) and the load vector F are iteratively assembled 190 over each angular frequency ω = 2πf to yield a system of 191 linear equations.These equations are solved as given in 192 equation 5 to obtain the solution vector X(χ, ω) which 193 will contain the displacement and pressure fields of the 194 solid and fluid parts of the poroelastic material respec-195 tively.196 { X(χ, ω)} = [ S(χ, ω)] −1 { F} (5) For normal incidence, assuming plane waves, the sound 197 absorption coefficient can be computed using the two-198 microphone method.Considering two closely spaced 199 points x 1 and x 2 in the air region, the complex pres-200 sure amplitudes in frequency domain P x1 and P x2 can 201 be obtained from {P} in X.The plane wave reflection 202 coefficient Rc can then be computed from these pressures 203 as, 204 Rc (χ, ω) = P x1 (χ, ω)e (−ikx2) − P x2 (χ, ω)e (−ikx1) −P x1 (χ, ω)e (ikx2) + P x2 (χ, ω)e (ikx1) (6) Here, k is the wave number given by ω/c air with c air 205 being the speed of sound in air.The sound absorption 206

212
the complex-valued Rc , (•) is the real part and (•) is 213 the complex conjugate operator.The gradient ∂ Rc ∂χi is obtained from ∂Px 1 ∂χi and ∂Px 2 ∂χi , which in-turn are two ele-215 ments from the derivative vector ∂ X ∂χi .To find ∂ X ∂χi , equa-216 tion 5 is differentiated to get the following expression.

369
Multiple machines were used to run the optimisation 370 tests, and in order to remove the machine-dependence 371 on runtime in Figure 2(a), the best-so-far absorption val-372 ues were tracked against the number of function evalua-373 tions, and runtimes were then computed by using aver-374 age time-per-function-evaluation clocked on a reference 375 machine.The reference machine used features an In-376 tel(R) Core(TM) i7-3820 CPU 3.6 GHz processor, 32 377 GB RAM and a 64-bit Windows 10 operating system 378 running Matlab2019b 55 .Scales indicating the number of 379 function evaluations are also provided for benchmarking 380 purposes.For all non-deterministic algorithms, as mul-381 tiple trials were conducted, the absorption values shown 382 in Figure 2(a) are averaged across the 31 trials after each 383 generation of the algorithm.384 Firstly, note that initial absorption levels are differ-385 ent for the algorithms.While the discrete algorithms 386 HC, GA, TABU, CMAd and DEd are initiated from ran-387 dom discrete solutions with α around 0.71, the continu-388 ous algorithms CMA, DE, SIMPf0 and SIMPf2 are ini-389 tiated from random continuous solutions with α around 390 0.65.CH2 starts from fully porous design domain with 391 α around 0.84 and CH1 starts from an empty (air-filled) 392 design domain with no absorption.393 One of the first things to note is that the CH2 algo-394 rithm does not produce an improvement from the fully 395 porous-filled solution and hence the best-so-far absorp-396 tion value stays the same for this problem.For low CPU-397 time budgets, SIMPf0 and SIMPf2 produce higher qual-398 ity solutions than all the other algorithms except CH2.399 SIMPf0 and SIMPf2 converge to a higher absorption than 400 the porous-filled CH2 solution in under 5 minutes on 401 this problem instance highlighting that gradient-based 402 methods can be time-efficient.After about 20 minutes of 403 runtime, HC produces better solutions on average than 404 SIMP, but the difference is small.405 After the designated budget of 4096 function evalua-406 tions (1366 gradient-included function evaluations), HC, 407 SIMPf2, TABU, SIMPf0 and CH1 produce the top tier 408 solutions.CMAd follows closely by producing slightly 409 better-quality solutions compared to fully filled CH2 so-410 lution towards the end.Whereas for DEd and GA, the 411 runtime performance was considerably poor.412 It is important to appreciate that the solutions from 413 continuous algorithms (CMA, DE, SIMPf0 and SIMPf2) 414 consider intermediate materials between the porous ma-415 terial and air χ i ∈ (0, 1] whereas the discrete algorithms 416 consider only porous material or air solutions χ i ∈ {0, 1}.417 Since the solutions are from different search spaces, the 418 absorption levels cannot be directly compared between 419 the two.Although the final shapes from continuous algo-420 rithms are desired to be 0 or 1, they are often not.Hence, 421 they are forced to be discrete using a simple round-off 422 filter, and the absorption values are recomputed.Such 423 rounding leads to a drop or surge in the absorption val-424 ues at the end of all continuous algorithms as can be 425 observed noticeably in CMA and DE plotlines in Fig-426 ure 2(a).The rounded absorptions indicated by the end 427 markers are also trial-averaged.Rounding leads to no 428 significant changes in SIMPf0 and SIMPf2 solutions for 429 this problem instance.For CMA and DE, the rounded 430 solution absorption values were poorer than SIMP solu-

FIG. 2 .
FIG. 2. (color online) Optimisation trials on problem instance 6: (a) Progress of best absorption found vs runtime (trialaveraged).For continuous algorithms, the solutions are discretised in the end.(b) Distribution of final solution absorption across trials.(c) Distribution of solution quality vs volume fraction (d) Sound absorption vs frequency for final shapes from select algorithms.(e) Best shapes from different trials from top four algorithms and their absorption.
(b).HC and SIMPf2 turn out to be the top-performing algorithms for this problem instance followed by TABU, SIMPf0 and CH1 in the second tier.DE, CMA and CMAd follow with all trials producing better solutions than the fully-filled CH2 solution.DEd and GA performed the poorest with no trials producing better than the fully-filled solution.The shapes produced from 10 of the trials from the top four algorithms are displayed in Figure 2(e).Most shapes seem to have a thin layer of air near the rigid backing as this allows removing elastic resonance around 500 Hz as can be observed from the absorption curves in Figure 2(d).Without filtering, SIMPf0 produces intricate designs near this thin air layer compared to SIMPf2.

499
CMAd and CH1 ranked first in four problem in-500 stances.Although CMAd ranked 8 th in problem instance 501 6, its overall performance across the problem instances 502 puts the algorithm in second place.Notably, in problem 503 instance 3, which considers a high static airflow resistivity 504 material, CMAd performed the best.This problem in-505 stance likely has many local optima and the performance 506 of CMAd indicates its global topology optimisation po-507 tential.The poor performance of the SIMP algorithms in 508 this problem instance is likely due to the multi-modality 509 of the objective function and premature convergence to 510 local optima.

Figure 4 .
For nondeterministic algorithms, the solution with the highest absorption among the 31 trials is shown.It is recalled that manufacturability restrictions and morphological filters are not imposed in this study except for SIMPf2.Results show both SIMPf0 and SIMPf2 produce similar shapes for most problem instances.For problem instance 1, all algorithms except SIMPf2 result in irregular shapes.The best quality shapes from most algorithms are flat layers of air and porous material towards the rigid wall with a somewhat circular air cavity in the front.GA and DE produced checkerbox shapes.Moreover, shapes from GA for all problem instances are degenerate.For problem instance 2, HC, CH1 and SIMPf0 produce the best shape with an almost porous material filled design domain except for a layer of air next to the rigid wall.CH1, SIMPf2, CMA, CMAd, TABU produced similar shapes.CH2 resulted in a fully-filled shape with slightly less absorption.

FIG. 4 .•
FIG. 4. Optimised shapes obtained from all algorithms for each problem instance.The shapes are discretised by rounding for continuous algorithms.The values of mean absorption across frequencies (α) are printed at the top of each shape in bold font along with porous material volume fraction (V f ) in parentheses.White and black represent air and the porous, respectively, with the acoustic input on the left and rigid backing on the right.

654•
The optimal shapes produced by algorithms that 655 use stochastic components (GA, CMA, CMAd, DE, 656 DEd) tend to be irregular and unconnected, and 657 hence they might need additional filtering tech-658 niques.Although HC produced higher sound ab-659 sorption solutions in general, the optimal shapes produced were not smooth and crisp.On the other constraint by terminating the construction after the 668 desired volume fraction is reached.The material 669 removal heuristic (CH2) often returns a fully filled

TABLE II .
Acoustic and elastic properties of materials used in the benchmark problems in TableI.Here, φ is the open porosity, Λ is the thermal characteristic length, Λ is the viscous characteristic length, σ is the static airflow resistivity, α∞ is the tortuosity, k 0 is the thermal permeability, ρ is the bulk density, E is the solid elastic modulus, ν is the Poisson's ratio and η is the dissipation factor.
a length D is designated as the design domain.The design domain is followed by a fixed domain, which is just an air region in this case with a length L. The design domain is discretised into nelx and nely finite elements

TABLE III .
Optimisation approaches tested (pseudocodes are included in the supplementary material) 491 lem instance 1, a special material previously used by Lee, 272 Kim, Kim, and Kang 24 (LKKK material) is used on a 273 coarser 10 × 10 discretisation.Note that the LKKK ma-274 terial may not representative of a physical material due 275 to the high tortuosity value of 7.8.Problem instance 2 276 features a 45 mm long design domain representative of 277 a typical building application.Problem instance 3 uses 278 an artificial material with a high static airflow resistiv-279 ity.In problem instance 4, a thin design domain of 2 280 cm, representative of a foam layer in an automotive ab-281 sorber, is considered.In problem instance 5, a thin layer 282 is optimised for high-frequency absorption.Among the 283 problem instances, problem instance 6 has a relatively 284 fine mesh size with 50 × 20 elements featuring a thicker 285 design domain optimised on a broad frequency range.286Otherthan 1 and 3, all problem instances use Melamine 287 foam for control.In problem instance 7, a single target 288 frequency is considered.289III.OPTIMISATION APPROACHES 290 Several gradient-free heuristic and metaheuristic ap-291 proaches, including well known and novel, are evaluated 292 in this study alongside the state-of-the-art gradient-based 293 approach SIMP.Henceforth in this paper, all the heuris-294 tic and metaheuristic approaches will be referred to as 295 algorithms, and they are not to be confused with exact al-296 gorithms as used by some authors.The algorithms tested 297 and their settings are summarised in Table III.303flippedlike in a raster scan (row-by-row) until the func-304 tion evaluation budget is used up.CH1 is a constructive 305 heuristic that starts from an air-filled solution and pro-306 gressively adds porous material in elements of best im-307 provement in absorption.Similarly, CH2 starts from a 30849Start from a random continuous solution, follow the SIMP procedure49; use density filter ft=2.Use SIMP penalty p = 3; move update -optimality criteria; move limit m = 0.2; Volume fraction limit V f = 1; Filter radius r min =2.porous material-filled solution and progressively removes porous material from the elements where the decrease in absorption is the least.SIMPf0 and SIMPf2 are solidisotropic-material-with-penalisation approaches49which use gradients of absorption to modify the solution at each step.While SIMPf2 uses density filtering, SIMPf0 uses 314 no filtering techniques.315 Four popular metaheuristic approaches are tested, 316 including genetic algorithm (GA), tabu search (TABU), 317 covariance-matrix-adaptation evolution strategy (CMA) 338 CH2).The continuous algorithms, which allow interme-339 diate materials between air and porous materials in each 340 element, are initiated from solutions generated by assign-341 ing a random number uniformly distributed between 0 342 and 1 to the topological design variables.Such random 343 initialisation is done to ensure a fair comparison making 344 no apriori assumptions about the solution.345 Some of the newly proposed approaches, namely, hill 346 climbing, constructive heuristics, and the discrete vari-347 ants of CMA evolution strategy and differential evolu-348 tion, in the specific way used here are tested for the 349 first time in topology optimisation.The others are well-350 established algorithms, and resources including surveys, 351 tutorials and code implementations can be easily reached.352 More specific implementation details are included in the 353 supplementary material.It is noted that a thorough 354 knowledge of all the algorithms is not essential to under-355 stand the findings.These algorithms can be thought of 356 as black-boxes that optimise the shape design by search-357 ing for the optimal assignment of the decision variables 358 χ to maximise α(χ).

TABLE IV .
The algorithms are ranked based on median values of optimised shape absorption (α * ) across trails.Lesser the average rank, the better is the performance of the algorithm.Algorithms are sorted based on the average of the ranks across problem instances.This ranking scheme is provided for a quick lookup only and is not meant to be a precise indicator of the performance.The ranking could change if more problem instances and algorithms are considered.