$\mathcal{L}_{\mathrm{UV}}\supset\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{-(4\kappa_{1}\dot{\xi}+\kappa_{2}\dot{\sigma})(\partial\xi)^{2}-2(\kappa_{2}\dot{\xi}+\kappa_{4}\dot{\sigma})(\partial^{\mu}\xi\partial_{\mu}\sigma)-(2\kappa_{3}\dot{\xi}+\kappa_{5}\dot{\sigma})(\partial\sigma)^{2}}.$
(98)

#####  $\expectationvalue{\xi\xi\xi}$ correlator.

The inflaton three-point function receives a contribution from the
$\kappa_{1}$ term in Eq. 98 and is given by

$\begin{split}&\expectationvalue{\xi(\mathbf{k}_{1})\xi(\mathbf{k}_{2})\xi(\mathbf{k}_{3})}\\\
&\qquad=-\frac{2\cdot
4i\kappa_{1}}{8k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\int_{-\infty}^{0}\frac{\differential{\eta}}{\eta^{4}}(\eta^{2}k_{1}^{2})\bqty{k_{2}^{2}k_{3}^{2}\eta^{4}+\eta^{2}(\mathbf{k}_{2}\cdot\mathbf{k}_{3})(1-ik_{23}\eta-
k_{2}k_{3}\eta^{2})}e^{ik_{123}\eta}\\\
&\qquad\qquad+\text{c.c.}+\text{perms.}\end{split}$ (99)

As in the main text, we have denoted $k_{ij\dotsm
m}=k_{i}+k_{j}+\dotsb+k_{m}$. To compute the above integral, we can use the
relation,

$\int_{-\infty}^{0}\differential{\eta}\eta^{\alpha}e^{ik\eta}=-\frac{i^{1+\alpha}}{k^{1+\alpha}}\,\Gamma(1+\alpha).$
(100)

After further simplification, this leads to

$\begin{split}&\expectationvalue{\xi(\mathbf{k}_{1})\xi(\mathbf{k}_{2})\xi(\mathbf{k}_{3})}\\\
&\qquad=\frac{2\kappa_{1}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{2}\cdot\mathbf{k}_{3})k_{1}^{2}}{k_{123}}\pqty{1+\frac{k_{23}}{k_{123}}+\frac{2k_{2}k_{3}}{k_{123}^{2}}}}+\text{2
perms.}\end{split}$ (101)

#####  $\langle\sigma\xi\xi\rangle$ correlator.

This correlator has the same structure as above since both $\sigma$ and $\xi$
are massless fields, except we need to account for the two terms proportional
to $\kappa_{2}$ in Eq. 98. The result is still given by

$\begin{split}&\expectationvalue{\sigma(\mathbf{k}_{1})\xi(\mathbf{k}_{2})\xi(\mathbf{k}_{3})}\\\
&\qquad=\frac{\kappa_{2}}{2k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{2}\cdot\mathbf{k}_{3})k_{1}^{2}}{k_{123}}\pqty{1+\frac{k_{23}}{k_{123}}+\frac{2k_{2}k_{3}}{k_{123}^{2}}}}\\\
&\qquad\qquad+\Bqty{\frac{\kappa_{2}}{2k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{2}\cdot\mathbf{k}_{1})k_{3}^{2}}{k_{123}}\pqty{1+\frac{k_{12}}{k_{123}}+\frac{2k_{1}k_{2}}{k_{123}^{2}}}}+(\mathbf{k}_{2}\leftrightarrow\mathbf{k}_{3})}.\end{split}$
(102)

#####  $\langle\sigma\sigma\xi\rangle$ correlator.

The structure remains same as above, except the individual contributions are
controlled by $\kappa_{3}$ and $\kappa_{4}$:

$\begin{split}&\expectationvalue{\sigma(\mathbf{k}_{1})\sigma(\mathbf{k}_{2})\xi(\mathbf{k}_{3})}\\\
&\qquad=\frac{\kappa_{3}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{1}\cdot\mathbf{k}_{2})k_{3}^{2}}{k_{123}}\pqty{1+\frac{k_{12}}{k_{123}}+\frac{2k_{1}k_{2}}{k_{123}^{2}}}}\\\
&\qquad\qquad+\Bqty{\frac{\kappa_{4}}{2k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{1}\cdot\mathbf{k}_{3})k_{2}^{2}}{k_{123}}\pqty{1+\frac{k_{13}}{k_{123}}+\frac{2k_{1}k_{3}}{k_{123}^{2}}}}+(\mathbf{k}_{1}\leftrightarrow\mathbf{k}_{2})}.\end{split}$
(103)

#####  $\langle\sigma\sigma\sigma\rangle$ correlator.

This case is identical to the $\langle\xi\xi\xi\rangle$ correlator, except
that $4\kappa_{1}$ is replaced by $\kappa_{5}$:

$\begin{split}&\expectationvalue{\sigma(\mathbf{k}_{1})\sigma(\mathbf{k}_{2})\sigma(\mathbf{k}_{3})}\\\
&\qquad=\frac{\kappa_{5}}{2k_{1}^{3}k_{2}^{3}k_{3}^{3}}\frac{\dot{\phi}_{0}}{\Lambda^{4}}\bqty{\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{2}\cdot\mathbf{k}_{3})k_{1}^{2}}{k_{123}}\pqty{1+\frac{k_{23}}{k_{123}}+\frac{2k_{2}k_{3}}{k_{123}^{2}}}}+\text{2
perms.}\end{split}$ (104)

## References

  * (1) E. Fermi, Tentativo di una Teoria dei Raggi $\beta$, Nuovo Cim. 11 (1934) 1.
  * (2) W. Heisenberg and H. Euler, Consequences of Dirac’s theory of positrons, Z. Phys. 98 (1936) 714, [physics/0605038].
  * (3) S. Weinberg, Phenomenological Lagrangians, Physica A 96 (1979) 327.
  * (4) W. E. Caswell and G. P. Lepage, Effective Lagrangians for Bound State Problems in QED, QCD, and Other Field Theories, Phys. Lett. B 167 (1986) 437.
  * (5) A. V. Manohar and M. B. Wise, Heavy Quark Physics.  Cambridge University Press, 2000.
  * (6) J. Polchinski, Effective field theory and the Fermi surface, in Theoretical Advanced Study Institute (TASI 92): From Black Holes and Strings to Particles, p. 235, 6, 1992.  hep-th/9210046.
  * (7) C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, An effective field theory for collinear and soft gluons: Heavy to light decays, Phys. Rev. D 63 (2001) 114020, [hep-ph/0011336].
  * (8) W. D. Goldberger and I. Z. Rothstein, Effective field theory of gravity for extended objects, Phys. Rev. D 73 (2006) 104029, [hep-th/0409156].
  * (9) A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis, and R. Rattazzi, Causality, analyticity and an IR obstruction to UV completion, JHEP 10 (2006) 014, [hep-th/0602178].
  * (10) T. N. Pham and T. N. Truong, Evaluation of the Derivative Quartic Terms of the Meson Chiral Lagrangian From Forward Dispersion Relation, Phys. Rev. D 31 (1985) 3027.
  * (11) B. Ananthanarayan, D. Toublan, and G. Wanders, Consistency of the chiral pion-pion scattering amplitudes with axiomatic constraints, Phys. Rev. D 51 (1995) 1093, [hep-ph/9410302].
  * (12) M. R. Pennington and J. Portoles, The chiral lagrangian parameters, $\ell_{1}$, $\ell_{2}$, are determined by the $\rho$-resonance, Phys. Lett. B 344 (1995) 399, [hep-ph/9409426].
  * (13) C. de Rham, S. Kundu, M. Reece, A. J. Tolley, and S.-Y. Zhou, Snowmass White Paper: UV Constraints on IR Physics, in 2022 Snowmass Summer Study, 3, 2022.  arXiv:2203.06805.
  * (14) C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, The Effective Field Theory of Inflation, JHEP 03 (2008) 014, [arXiv:0709.0293].
  * (15) L. Senatore and M. Zaldarriaga, The Effective Field Theory of Multifield Inflation, JHEP 04 (2012) 024, [arXiv:1009.2093].
  * (16) T. Grall and S. Melville, Positivity Bounds without Boosts, arXiv:2102.05683.
  * (17) P. Creminelli, O. Janssen, and L. Senatore, Positivity bounds on effective field theories with spontaneously broken Lorentz invariance, JHEP 09 (2022) 201, [arXiv:2207.14224].
  * (18) Planck Collaboration, Y. Akrami et al., Planck 2018 results. IX. Constraints on primordial non-Gaussianity, Astron. Astrophys. 641 (2020) A9, [arXiv:1905.05697].
  * (19) N. N. Bogoliubov, D. V. Shirkov, and S. Chomet, Introduction to the theory of quantized fields, vol. 59.  Interscience New York, 1959.
  * (20) H. J. Bremermann, R. Oehme, and J. G. Taylor, Proof of Dispersion Relations in Quantized Field Theories, Phys. Rev. 109 (1958) 2178\.
  * (21) K. Hepp, On the analyticity properties of the scattering amplitude in relativistic quantum field theory, Helv. phys. acta 37 (1964) 639\.
  * (22) C. de Rham, S. Melville, A. J. Tolley, and S.-Y. Zhou, UV complete me: Positivity Bounds for Particles with Spin, JHEP 03 (2018) 011, [arXiv:1706.02712].
  * (23) D. Baumann, D. Green, H. Lee, and R. A. Porto, Signs of Analyticity in Single-Field Inflation, Phys. Rev. D 93 (2016) 023523, [arXiv:1502.07304].
  * (24) S. Melville and J. Noller, Positivity bounds from multiple vacua and their cosmological consequences, arXiv:2202.01222.
  * (25) V. Chandrasekaran, G. N. Remmen, and A. Shahbazi-Moghaddam, Higher-Point Positivity, JHEP 11 (2018) 015, [arXiv:1804.03153].
  * (26) C. Cheung and G. N. Remmen, Positive Signs in Massive Gravity, JHEP 04 (2016) 002, [arXiv:1601.04068].
  * (27) G. N. Remmen and N. L. Rodd, Consistency of the Standard Model Effective Field Theory, JHEP 12 (2019) 032, [arXiv:1908.09845].
  * (28) G. N. Remmen and N. L. Rodd, Flavor Constraints from Unitarity and Analyticity, Phys. Rev. Lett. 125 (2020) 081601, [arXiv:2004.02885]. [Erratum: Phys. Rev. Lett. 127 149901 (2021)].
  * (29) G. N. Remmen and N. L. Rodd, Signs, spin, SMEFT: Sum rules at dimension six, Phys. Rev. D 105 (2022) 036006, [arXiv:2010.04723].
  * (30) G. N. Remmen and N. L. Rodd, Spinning sum rules for the dimension-six SMEFT, JHEP 09 (2022) 030, [arXiv:2206.13524].
  * (31) C. Zhang and S.-Y. Zhou, Convex Geometry Perspective on the (Standard Model) Effective Field Theory Space, Phys. Rev. Lett. 125 (2020) 201601, [arXiv:2005.03047].
  * (32) T. Trott, Causality, unitarity and symmetry in effective field theory, JHEP 07 (2021) 143, [arXiv:2011.10058].
  * (33) N. Arkani-Hamed, Y.-t. Huang, J.-Y. Liu, and G. N. Remmen, Causality, unitarity, and the weak gravity conjecture, JHEP 03 (2022) 083, [arXiv:2109.13937].
  * (34) X. Li, H. Xu, C. Yang, C. Zhang, and S.-Y. Zhou, Positivity in Multifield Effective Field Theories, Phys. Rev. Lett. 127 (2021) 121601, [arXiv:2101.01191].
  * (35) N. Arkani-Hamed, T.-C. Huang, and Y.-T. Huang, The EFT-Hedron, JHEP 05 (2021) 259, [arXiv:2012.15849].
  * (36) A. Nicolis, R. Rattazzi, and E. Trincherini, Energy’s and amplitudes’ positivity, JHEP 05 (2010) 095, [arXiv:0912.4258]. [Erratum: JHEP 11, 128 (2011)].
  * (37) H. P. Stapp, Crossing, hermitian analyticity, and the connection between spin and statistics, J. Math. Phys. 9 (1968) 1548.
  * (38) J. L. Miramontes, Hermitian analyticity versus real analyticity in two-dimensional factorized $S$-matrix theories, Phys. Lett. B 455 (1999) 231, [hep-th/9901145].
  * (39) N. Arkani-Hamed, P. Creminelli, S. Mukohyama, and M. Zaldarriaga, Ghost inflation, JCAP 04 (2004) 001, [hep-th/0312100].
  * (40) B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev. Mod. Phys. 78 (2006) 537, [astro-ph/0507632].
  * (41) K. A. Malik and D. Wands, Cosmological perturbations, Phys. Rept. 475 (2009) 1, [arXiv:0809.4944].
  * (42) A. D. Linde and V. F. Mukhanov, Nongaussian isocurvature perturbations from inflation, Phys. Rev. D 56 (1997) R535–R539, [astro-ph/9610219].
  * (43) K. Enqvist and M. S. Sloth, Adiabatic CMB perturbations in pre-Big-Bang string cosmology, Nucl. Phys. B 626 (2002) 395, [hep-ph/0109214].
  * (44) D. H. Lyth and D. Wands, Generating the curvature perturbation without an inflaton, Phys. Lett. B 524 (2002) 5, [hep-ph/0110002].
  * (45) T. Moroi and T. Takahashi, Effects of cosmological moduli fields on cosmic microwave background, Phys. Lett. B 522 (2001) 215, [hep-ph/0110096]. [Erratum: Phys.Lett.B 539, 303 (2002)].
  * (46) K. M. Smith, L. Senatore, and M. Zaldarriaga, Optimal analysis of the CMB trispectrum, arXiv:1502.00635.
  * (47) D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 (2016) 1, [arXiv:1510.07633].
  * (48) Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [arXiv:1807.06211].
  * (49) C. de Rham, S. Melville, A. J. Tolley, and S.-Y. Zhou, Massive Galileon Positivity Bounds, JHEP 09 (2017) 072, [arXiv:1702.08577].
  * (50) B. Bellazzini, F. Riva, J. Serra, and F. Sgarlata, Beyond Positivity Bounds and the Fate of Massive Gravity, Phys. Rev. Lett. 120 (2018), no. 16 161101, [arXiv:1710.02539].
  * (51) J. M. Maldacena, Non-Gaussian features of primordial fluctuations in single field inflationary models, JHEP 05 (2003) 013, [astro-ph/0210603].
  * (52) D. Langlois and B. van Tent, Isocurvature modes in the CMB bispectrum, JCAP 07 (2012) 040, [arXiv:1204.5042].
  * (53) S. Kumar and R. Sundrum, Heavy-Lifting of Gauge Theories By Cosmic Inflation, JHEP 05 (2018) 011, [arXiv:1711.03988].
  * (54) X. Chen, Primordial Non-Gaussianities from Inflation Models, Adv. Astron. 2010 (2010) 638979, [arXiv:1002.1416].
  * (55) S. Weinberg, Quantum contributions to cosmological correlations, Phys. Rev. D 72 (2005) 043514, [hep-th/0506236].


[1]Peter Dunning

1]School of Engineering, University of Aberdeen, King’s College, Aberdeen,
AB24 3UE, UK

# Achieving binary topology optimization solutions via automatic projection
parameter increase

<EMAIL_ADDRESS>[

###### Abstract

A method is created to automatically increase the threshold projection
parameter in three-field density-based topology optimization to achieve a near
binary design. The parameter increase each iteration is based on an
exponential growth function, where the growth rate is dynamically changed
during optimization by linking it to the change in objective function. This
results in a method that does not need to be tuned for specific problems, or
optimizers, and the same set of hyper-parameters can be used for a wide range
of problems. The effectiveness of the method is demonstrated on several 2D
benchmark problems, including linear buckling and geometrically nonlinear
problems.

###### keywords:

Topology optimization, Threshold projection, Linear buckling, Nonlinear
geometry

##  1 Introduction

Topology optimization aims to find the ideal layout, shape and size of a
structure for given objectives and constraints. These problems are typically
framed by deciding where material should exist within the design domain,
leading to a binary value problem, where a value of 1 indicates material, and
0 is empty space (or void). However, solving binary optimization problems is
challenging, and it is therefore popular to relax the problem using continuous
design variables, as introduced by [1]. The most popular topology optimization
approach using this principle is the Simple Isotropic Material with
Penalization (SIMP) method, where a pseudo-density value is assigned to each
element, which can vary between 0 (void) and 1 (solid) [2, 3].

The pseudo-density method is usually coupled with filtering techniques to
avoid numerical issues such as checkerboard patterns and mesh dependency [4].
A popular approach is to use a density filter, as this results in a consistent
problem formulation, where the derivatives of the optimization objectives and
constraints can be explicitly computed [5, 6]. However, using a density filter
can lead to solutions with a significant amount intermediate pseudo-density
values (between 0 and 1), far from the ideal binary design [7]. Intermediate,
or ”gray” elements can also appear in problems where the optimizer may exploit
the properties of these elements to improve the objective, or meet a
constraint. For example, intermediate material may appear to create a low-
density core to improve buckling performance, similar to a sandwich structure.

To force a design to be binary, a further operation can be used to threshold
pseudo-density values such that any value above the threshold becomes 1
(solid) and any below becomes 0 (void) [8]. However, to use a threshold during
gradient-based optimization, the threshold projection function is also relaxed
to take a continuous form [8, 9, 7]. These relaxed formulations use a
projection parameter that defines the sharpness of the projection, where a
small value results in little, or no change to pseudo-density values, and a
value of infinity results in a binary design (albeit non-differentiable). This
introduces an optimization hyper-parameter that must be chosen by the user.
From experience, many researchers use a continuation scheme for the projection
parameter, where it is kept at a small value (e.g. 1) during early
optimization iterations to allow a design to emerge. It is then incrementally
increased to push the design towards a binary solution, for example see [10,
11]. These continuation schemes are usually ad-hoc and heuristic, needing to
be tuned for each new problem. It is also not clear if these heuristic schemes
achieve the best balance between computational cost and performance of the
design (including how discrete it is).

An alternative approach was introduced by [12] where the threshold projection
parameter is included as an extra element-wise design variable. The lower
limit for each variable is dynamically updated based on the current grayness
measure (a measure on non-discreetness) for the element. In addition, material
penalisation parameters are treated in the same way. The method is an
improvement on heuristic continuation schemes, and works well for the problem
studied (fibre reinforced structures with stress constraints), but still has
several tuning parameters and it is unclear if the method would work on a
wider class of problems.

Another idea was recently proposed by [13] where a surrogate optimization
method is used as an alternative to manually determining topology optimization
hyper-parameters, including continuation parameters. This involves solving the
original topology optimization problem as an inner problem in a nested loop,
where hyper-parameters are optimized in an outer loop aided by machine
learning techniques. This can reduce the effort required to manually tune
hyper-parameters, but may not be feasible for computationally expensive
problems (e.g. nonlinear and multi-disciplinary problems).

Note that several topology optimization methods have been developed that use
binary design variables, such as Bi-directional Evolutionary Structural
Optimization (BESO) [14], and more recently Topology Optimization of Binary
Structures (TOBS) [15], amongst others. However, these methods have their own
heuristic hyper-parameters that may be problem dependent. There are also some
methods that allow intermediate densities in the solution by linking their
properties to a porous micro-structure, as originally explored in the seminal
paper by [1]. This idea has been recently applied to develop multi-scale
methods, see for example [16, 17].

In this paper, a simple alternative method for obtaining near binary topology
optimization designs is proposed. Starting with the three-field density-based
method (where a density filter is combined with threshold projection) a method
is developed to automatically increase the threshold projection parameter
during optimization by linking it to the progress of the objective function.
The three-field density-based method is reviewed in Section 2.1 and the new
automatic continuation scheme is introduced in Sections 2.2 and 2.3. The
method is demonstrated using several benchmark problems in Section 3,
including comparison with heuristic, ad-hoc continuation schemes from the
literature.

##  2 Methodology

###  2.1 Three-field topology optimization

The three-field density-based topology optimization method starts with the
design variable field, where a scalar value $x_{e}\in[0,1]$ is assigned to
each element, $e$, in the mesh. A density-based filter is then applied to
avoid well-known issues, such as checkerboard patterns and mesh dependent
solutions [5, 6].

$\tilde{x}_{e}=\frac{\sum_{i\in N_{i,e}}{H_{i,e}x_{i}}}{\sum_{i\in
N_{i,e}}{H_{i,e}}}$ (1)

where $N_{i,e}$ is the set of elements $i$ where the distance between the
center of element $i$ and center of element $e$, $d(e,i)$, is less than the
filter radius $r_{\min}$. $H_{i,e}$ is a weight factor, defined as:

$H_{i,e}=\max(0,r_{\min}-d(e,i))$ (2)

The filtered field $\tilde{\boldsymbol{x}}$ creates a band of intermediate
density values around the boundary of the structure. Thus, a threshold
projection function is used to obtain a more crisp (or binary) design
description [9, 7]:

$\bar{x}_{e}=\frac{\tanh{(\beta\eta)}+\tanh{(\beta(\tilde{x}_{e}-\eta))}}{\tanh{(\beta\eta)}+\tanh{(\beta(1-\eta))}}$
(3)

where $\bar{x}$ is the physical density field, $\beta$ is the projection
parameter that controls the sharpness of the projection, and $\eta$ is the
threshold value, typically set as $\eta=0.5$.

The value of $\bar{x}_{e}$ is used to compute the physical properties of the
element, where the volume is directly proportional to $\bar{x}_{e}$ and the
modified SIMP method is used to compute the Young’s Modulus:

$E_{e}(\bar{x}_{i})=E_{\min}+\bar{x}_{e}^{p}(E_{0}-E_{\min})$ (4)

where $E_{0}$ is the Young’s modulus of the material, $E_{\min}$ is a small
value to prevent the stiffness matrix becoming singular and $p$ is the
penalization factor, where the typical value of $p=3$ is used in this paper.
Note that the physical pseudo-density values $\bar{x}_{e}$ are used to plot
solutions.

###  2.2 Binary solutions via threshold projection

The solution to a topology optimization problem should ideally be binary,
indicating where material exists ($\bar{x}_{e}=1$), and where is does not
($\bar{x}_{e}=0$). Therefore, we aim to avoid intermediate density values in
the solution, although they are necessary during optimization to make the
problem differential such that it can be solved by efficient gradient-based
optimizers.

The gray level indicator function [18] can be used to measure how close a
design is to being binary:

$G(\bar{x})=\frac{4}{N}\sum_{e=1}^{N}{\bar{x}_{e}(1-\bar{x}_{e})}$ (5)

where $N$ is the number of elements. Smaller values of $G(\bar{x})$ indicate
the design is closer to being completely binary, with $G(\bar{x})=0$
characterizing a true binary design.

The threshold projection function, Eq. (3), can be used to reduce the gray
level of a design by increasing the projection parameter $\beta$, thus making
it closer to being binary, regardless of other parameters, such as the
penalization power $p$ in Eq. (4). When $\beta$ is increased, then if
$\eta<\bar{x}_{e}<1$ it increases to be closer to 1, and if
$0<\bar{x}_{e}<\eta$ it decreases to become closer to 0. In the limit of
$\beta=\infty$, only three physical density values are possible:
$\bar{x}_{e}\in\\{0,0.5,1\\}$, where $\bar{x}_{e}=0.5$ is only possible when
$\tilde{x}_{e}=\eta$, thus the design is practically binary.

However, starting the optimization with a large value of $\beta$ can cause
convergence problems and limits design exploration. Therefore, researchers
often start with a small value of $\beta$ and employ continuation schemes to
gradually increase it during optimization. This allows a design to emerge in
the early iterations, whilst being projected close to a binary design at the
end. Current continuation approaches for the projection parameter are
dependent on several hyper-parameters, which are often tuned for each new
problem because they can significantly influence the performance of the
solution and number of iterations. This issue is addressed in the following
section.

###  2.3 Automatic projection parameter increase

Some continuation schemes increase $\beta$ after a number of iteration steps.
However, this can cause a discontinuity in the optimization process. There are
also questions on how frequent $\beta$ should be increased, by how much should
it be increased, and when to stop increasing. The method proposed here
addresses the question of how frequent, by checking for $\beta$ increase every
iteration. The amount of increase, $\Delta\beta$, is linked to the rate of
decrease in the objective function, as this indicates whether a design is
still emerging (large relative change in objective), or is settling (small
relative change in objective). Assuming the objective is to minimize a
function, $f$, the increase in $\beta$ at iteration $k$ is calculated by:

$\Delta\beta^{k}=\max\left(-\frac{\gamma}{2}~{}\frac{f^{k}+f^{k-1}}{f^{k}-f^{k-1}},0\right)$
(6)

where $\gamma$ is a tuning parameter that sets how fast $\beta$ increases. A
smaller value leads to a slower increase, allowing the optimizer more time to
refine the design, before being pushed towards a binary solution, but
potentially taking more iterations. Whereas a bigger value pushes the design
towards a binary solution quicker, at the possible sacrifice in objective
function value. The value of $\beta$ in the next iteration is then:

$\beta^{k+1}=\beta^{k}+\min\left(\Delta\beta^{k},\Delta\beta_{max}\right)$ (7)

where $\Delta\beta_{max}$ is the maximum possible increase. Note that a
maximum value on $\beta$ should also be set to avoid numerical issues with
extremely large numbers, and that $\beta=1.0$ in the first iteration.

Although the proposed automatic increase method is not free from tuning
parameters, testing has shown that setting $\gamma=10^{-4}$ and
$\Delta\beta_{max}=0.2\beta^{k}$ works well for a variety of problems, as
shown in Section 3 below. Therefore, the need to tune these parameters for new
problems is reduced, compared to existing methods, and in most cases tuning is
not required.

The final part of the proposed method is to use the current value of the gray
level indicator function, Eq. (5), to determine if $\beta$ should be increased
and if the optimization process can be stopped. Thus, $\beta$ is only
increased if, $G(\bar{x})>\epsilon$, and the optimization process should
continue until, $G(\bar{x})\leq\epsilon$, and all other stopping, or
convergence criteria are met. The value of $\epsilon$ can be set by the user
to trade-off computational time (bigger $\epsilon$ leads to fewer iterations)
with how close the solution is to a binary design (smaller $\epsilon$ means
less gray in the solution). For all examples in this work $\epsilon=0.01$,
providing a reasonable balance between computational cost and grayness.

##  3 Examples

Several benchmark problems are used to demonstrate the method described in
Section 2.3. The first is a compressed column, Figure 1a, and the second is a
cantilever with aspect ratio 4, Figure 1b.

The compressed column is discretized using $120\times 240$ four-node plane
stress elements with a thickness of 1. The applied load is $10^{-3}$, which is
spread over 8 elements, $E_{0}=1$, $E_{min}=10^{-6}$ and Poisson’s ratio
$\nu=0.3$. In addition, a small region around the applied load ($8\times 4$
elements) is fixed to remain part of the structure by setting $\bar{x}_{e}=1$.
Two different filter radii in Eq. (2) (4, or 8 times element edge length, $h$)
are used to test the ability of the proposed method to achieve near binary
designs with different length scales, noting that a larger filter radius
creates more gray elements around the boundary.

The cantilever is discretized using three mesh sizes: $160\times 40$,
$320\times 80$, and $640\times 160$ to test the ability of the proposed method
to achieve a near binary design whilst refining the mesh. This problem is
solved considering geometrically nonlinear behaviour using the corotational
element proposed by [19]. Two different applied loads are used: $2\times
10^{5}$ and $3\times 10^{5}$. The element thickness is 0.1, $E_{0}=3\times
10^{9}$, $E_{min}=3$ and $\nu=0.4$. The filter radius in Eq. (2) is set to
0.075.

Note that the same stopping criteria are used for all problems, where the gray
level indicator function, Eq. (5), must be below 0.01, the absolute value of
relative change in objective function between two iterations must be below
$10^{-5}$, and all constraints are satisfied.

Figure 1: Example problems, a) compressed column, b) cantilever.

###  3.1 Compressed column results

The compressed column example is solved for two different problems. The first
is to maximize the buckling load factor (or minimize the inverse) with a 35%
volume constraint, using a uniform density starting design equal to the volume
constraint. The second is to minimize volume with a buckling load factor
constraint (minimum value of 15) and compliance constraint - set at 2 times
the compliance for the initial fully solid design. Both problems are solved
using the Matlab code developed by [20], which includes an inbuilt
continuation method for incrementally increasing $\beta$ during optimization.
For comparison, the problems are also solved using this in-built continuation
method. The default settings are that $\beta=1$ for the first 400 iterations,
it is then increased by 2 every 25 iterations until a maximum value of 25. A
further comparison is made to a “hand tuned” modified version of this approach
following analysis of the optimization history using the default settings,
where initially $\beta=1$ for the first 200 iterations, it is then increased
by 2 every 25 iterations up to a possible maximum value of 500. The problems
are solved using the inbuilt optimally criteria method in the Matlab code
[20]. To avoid issues with mode switching, 30 modes are considered for the
maximization of buckling load factor problem, and 20 for the buckling
constraint in the volume minimization problem. In both cases, the maximum
buckling load factor is approximated using a K-S function with an aggregation
parameter of 160.

The solutions for the maximum buckling load factor problem are shown in Figure
2. Note that the solutions for the default continuation approach are stopped
after 2000 iterations, as the maximum value of $\beta=25$ could not meet the
gray level stopping criteria - i.e. the default maximum value of $\beta$ is
too small to achieve a near binary design. This is most noticeable in the
solution using a filter radius of $8h$ in Figure 2d, where large areas of
intermediate density material exist, demonstrating the challenge in achieving
near binary designs for this problem. The modified scheme allows for a bigger
$\beta$ values and a solution for the filter radius $8h$ problem that meets
the gray level criterion is obtained after 2751 iterations, Figure 2e. The
proposed automatic scheme also produces a design that meets the gray level
criterion, Figure 2f, but it only needs 1029 iterations. The designs are
slightly different, but the maximum buckling load factor for the design
obtained using the automatic scheme is slightly better (16.26 versus 16.22 -
see Table 1). Similar results can be seen when the filter radius is $4h$,
although the reduction in number of iterations comparing the automatic and
modified schemes is smaller, as the smaller filter radius makes it easier to
meet the gray level stopping criterion. This suggests that the modified scheme
could be further improved for the $8h$ filter radius problem, but would
require further optimization runs and hand tuning, which are not required for
the proposed automatic scheme.

Figure 2: Compressed column results for maximization of buckling load factor.
Solutions with filter radius = $4h$: a) Default continuation scheme, b)
modified scheme, c) automatic scheme. Solutions with filter radius = $8h$: d)
Default continuation scheme, e) modified scheme, f) automatic scheme.

Optimization histories for all problems, including the progress of the the
gray level function and $\beta$ value, are shown in Figure 3. They demonstrate
how the automatic scheme modifies $\beta$ in relation to the progress of the
objective function, and is able to push the design towards a near binary
solution (as shown by the gray level indicator) quickly, without compromising
the performance of the solution, in comparison with the other schemes.

Figure 3: Compressed column convergence history for maximization of buckling load factor. Convergence of volume objective for: a) filter radius = $4h$, c) filter radius = $8h$. Convergence of gray level indicator function, $G(x)$ , and threshold projection parameter, $\beta$, for b) filter radius = $4h$, d) filter radius = $8h$. Filter radius = $4h$ |  |  |  |
---|---|---|---|---
Continuation method | Buckling load factor | Iterations | $\beta$ | Gray value
Default | 16.72 | 2000 | 25 | $3.68^{-2}$
Modified | 16.54 | 941 | 61 | $1.00^{-2}$
Automatic | 16.72 | 728 | 229.7 | $9.10^{-3}$
Filter radius = $8h$ |  |  |  |
Continuation method | Buckling load factor | Iterations | $\beta$ | Gray value
Default | 15.75 | 2000 | 25 | $1.27^{-1}$
Modified | 16.22 | 2751 | 207 | $9.87^{-3}$
Automatic | 16.26 | 1029 | 385.9 | $9.93^{-3}$
Table 1: Results for compressed column for maximization of buckling load
factor.

The solutions of the volume minimization problem are shown in Figure 4. Again,
the default scheme is unable to produce designs that meet the gray level
stopping criteria, and are stopped after 2000 iterations. The modified and
automatic schemes do produce designs that meet the criterion, although the
automatic scheme achieves this in significantly fewer iterations
(approximately a third - see Table 2). However, for this examples, the slower,
modified scheme obtains slightly better solutions with smaller volume
fractions (by approximately 1-2%). This is evident from the simpler designs
obtained by the modified scheme (e.g. compare Figure 4b with c, or Figure 4e
with f). However, this small increase in performance, comes at significantly
more computational cost, as seen by the required number of iterations in Table
2. The optimization histories, Figure 5, show the same as the previous
problem, where the automatic scheme is able to adapt to the optimization
progress and push the design towards a near binary solution quicker than the
other two schemes, which rely on a predetermined $\beta$ increase scheme.

Figure 4: Compressed column results for minimization of volume. Solutions with filter radius = $4h$: a) Default continuation scheme, b) modified scheme, c) automatic scheme. Solutions with filter radius = $8h$: d) Default continuation scheme, e) modified scheme, f) automatic scheme. Figure 5: Compressed column convergence history for minimization of volume. Convergence of volume objective for: a) filter radius = $4h$, c) filter radius = $8h$. Convergence of gray level indicator function, $G(x)$, and threshold projection parameter, $\beta$, for b) filter radius = $4h$, d) filter radius = $8h$. Filter radius = $4h$ |  |  |  |
---|---|---|---|---
Continuation method | Volume fraction | Iterations | $\beta$ | Gray value
Default | 0.384 | 2000 | 25 | $1.36^{-2}$
Modified | 0.384 | 606 | 35 | $9.45^{-3}$
Automatic | 0.388 | 232 | 54.5 | $7.29^{-3}$
Filter radius = $8h$ |  |  |  |
Continuation method | Volume fraction | Iterations | $\beta$ | Gray value
Default | 0.393 | 2000 | 25 | $3.39^{-2}$
Modified | 0.387 | 1076 | 73 | $9.71^{-3}$
Automatic | 0.394 | 334 | 114.8 | $8.59^{-3}$
Table 2: Results for compressed column for minimization of volume.

###  3.2 Cantilever results

The cantilever problem, Figure 1b, is solved to minimize nonlinear end
compliance, with a 40% volume constraint, and a constraint on stability, where
the minimum critical load factor is 2. The stability constraint is formulated
using an eigenvalue analysis at a load factor of 1, which provides a good
approximation of the true critical load factor, whilst being robust and
efficient [11]. The eigenvalue analysis also allows for multiple modes to be
included via a K-S function to help avoid issues with mode switching. For all
cantilever problems, 6 modes are used with as K-S aggregation parameter of 50.
The problem is solved using the Method of Moving Asymptotes [21], with a
conservative move limit of 0.05.

The solutions for a load of $2\times 10^{5}$ and the three different mesh
sizes are shown in Figure 6. The designs are very similar, due to the mesh
independent filter, but with denser meshes producing smoother solutions. The
automatic continuation scheme for $\beta$ is able to meet the gray level
indicator stopping criterion, by choosing a larger $\beta$ value for a denser
mesh - see Figure 7d and Table 3. It is also noticed that the value of $\beta$
follows a similar increase for about the first 90 iterations for all meshes.
After this, $\beta$ is increased quicker and to a larger value for denser
meshes to meet the gray level stopping criterion. Convergence histories for
the three meshes are shown in Figure 7, which show smooth convergence, except
for some oscillations in the stability constraint. This is mainly due to some
mode switching, that could be improved by including more eigenvalues in the
constraint formulation. However, feasible designs are obtained in a reasonable
number of iterations. Finally, it is observed that the objective (end
compliance) reduces as the mesh is refined, primarily due to the better
resolution of topological features, but this comes at more computational cost,
for both the finite elements analysis and number of iterations, see Table 3.

Figure 6: Cantilever results with load = $2\times 10^{5}$ using mesh size: a) $160\times 40$, b) $320\times 80$, c) $640\times 160$. Figure 7: Cantilever with load = $2\times 10^{5}$. Convergence history for mesh size: a) $160\times 40$, b) $320\times 80$, c) $640\times 160$. d) Convergence of gray level indicator function, $G(x)$, and threshold projection parameter, $\beta$. Mesh size | End compliance | Iterations | $\beta$ | Gray value
---|---|---|---|---
$160\times 40$ | 57,068 | 149 | 38.7 | $8.77\times 10^{-3}$
$320\times 80$ | 56,155 | 154 | 63.7 | $8.67\times 10^{-3}$
$640\times 160$ | 55,840 | 202 | 93.2 | $9.59\times 10^{-3}$
Table 3: Results for cantilever with load = $2\times 10^{5}$.

Cantilever solutions for a load of $3\times 10^{5}$ and the three different
mesh sizes are shown in Figure 8. Again, the designs are very similar,
although the finest mesh ($640\times 160$) has some additional thin members,
which cannot be represented accurately using the coarser meshes. Compared with
the results for a load of $2\times 10^{5}$, similar observation are made about
the final value of $\beta$ being larger for a finer mesh to meet the gray
level stopping criterion, and that more iterations are needed for the finer
mesh, but the objective function is lower - see Table 4. Convergence histories
are also shown in Figure 9, which again show some oscillations in the
stability constraints. This is most noticeable for the finest mesh ($640\times
160$) and is caused by mode switching due to the longest thin member becoming
even thinner. This can be eliminated by increasing the number of modes in the
constraint formulation to 12, but the result for 6 modes is presented for
consistency, and a feasible design is obtained anyway.

Figure 8: Cantilever results with load = $3\times 10^{5}$ using mesh size: a) $160\times 40$, b) $320\times 80$, c) $640\times 160$. Figure 9: Cantilever with load = $3\times 10^{5}$. Convergence history for mesh size: a) $160\times 40$, b) $320\times 80$, c) $640\times 160$. d) Convergence of gray level indicator function, $G(x)$, and threshold projection parameter, $\beta$. Mesh size | End compliance | Iterations | $\beta$ | Gray value
---|---|---|---|---
$160\times 40$ | 134,415 | 140 | 35.0 | $7.73\times 10^{-3}$
$320\times 80$ | 131,305 | 159 | 65.1 | $8.56\times 10^{-3}$
$640\times 160$ | 129,422 | 221 | 101.5 | $9.67\times 10^{-3}$
Table 4: Results for cantilever with load = $3\times 10^{5}$.

##  4 Conclusions

A simple method is proposed to to achieve near binary designs for three-field
density-based topology optimization, by automatically increasing the value of
the threshold projection parameter during optimization. The method has two key
parts. Firstly the increase in the threshold projection parameter is linked to
the progress of the objective function, and secondly the gray level indicator
function is used as a stopping criterion.

The method is tested on some 2D benchmark problems for linear buckling and a
geometrically nonlinear case. The results show that the proposed method can
achieve near binary solutions (as measured by the gray level indicator
function). When compared with a conventional continuation approach, where the
increase in $\beta$ is predetermined and not linked to optimization progress,
improved results in terms of reduced grayness, and/or reduced number of
iterations is observed. Further testing is required to determine the utility
of the proposed method for a wider range of problems, including different
physics, and when using different optimizers.

## References

  *   * Bendsøe and Kikuchi [1988] Bendsøe, M.P., Kikuchi, N.: Generating optimal topologies in structural design using a homogenization method. Computer methods in applied mechanics and engineering 71(2), 197–224 (1988)
  * Bendsøe [1989] Bendsøe, M.P.: Optimal shape design as a material distribution problem. Structural optimization 1, 193–202 (1989)
  * Rozvany et al. [1992] Rozvany, G.I., Zhou, M., Birker, T.: Generalized shape optimization without homogenization. Structural optimization 4, 250–252 (1992)
  * Sigmund and Petersson [1998] Sigmund, O., Petersson, J.: Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural optimization 16, 68–75 (1998)
  * Bruns and Tortorelli [2001] Bruns, T.E., Tortorelli, D.A.: Topology optimization of non-linear elastic structures and compliant mechanisms. Computer methods in applied mechanics and engineering 190(26-27), 3443–3459 (2001)
  * Bourdin [2001] Bourdin, B.: Filters in topology optimization. International journal for numerical methods in engineering 50(9), 2143–2158 (2001)
  * Wang et al. [2011] Wang, F., Lazarov, B.S., Sigmund, O.: On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization 43(6), 767–784 (2011)
  * Guest et al. [2004] Guest, J.K., Prévost, J.H., Belytschko, T.: Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International journal for numerical methods in engineering 61(2), 238–254 (2004)
  * Xu et al. [2010] Xu, S., Cai, Y., Cheng, G.: Volume preserving nonlinear density filter based on heaviside functions. Structural and Multidisciplinary Optimization 41(4), 495–505 (2010)
  * Ferrari and Sigmund [2020] Ferrari, F., Sigmund, O.: A new generation 99 line matlab code for compliance topology optimization and its extension to 3d. Structural and Multidisciplinary Optimization 62, 2211–2228 (2020)
  * Dunning [2023] Dunning, P.: Stability constraints for geometrically nonlinear topology optimization. Structural and multidisciplinary optimization (2023)
  * da Silva et al. [2023] Silva, A.L.F., Salas, R.A., Silva, E.C.N.: Topology optimization of fiber reinforced structures considering stress constraint and optimized penalization. Composite Structures 316, 117006 (2023)
  * Ha and Carstensen [2024] Ha, D., Carstensen, J.: Automatic hyperparameter tuning of topology optimization algorithms using surrogate optimization (2024) https://doi.org/10.21203/rs.3.rs-3944293/v1
  * Huang and Xie [2007] Huang, X., Xie, Y.: Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. Finite elements in analysis and design 43(14), 1039–1049 (2007)
  * Sivapuram and Picelli [2018] Sivapuram, R., Picelli, R.: Topology optimization of binary structures using integer linear programming. Finite Elements in Analysis and Design 139, 49–61 (2018)
  * Christensen et al. [2023] Christensen, C.F., Wang, F., Sigmund, O.: Topology optimization of multiscale structures considering local and global buckling response. Computer Methods in Applied Mechanics and Engineering 408, 115969 (2023)
  * Hübner et al. [2023] Hübner, D., Wein, F., Stingl, M.: Two-scale optimization of graded lattice structures respecting buckling on micro-and macroscale. Structural and Multidisciplinary Optimization 66(7), 163 (2023)
  * Sigmund [2007] Sigmund, O.: Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33(4-5), 401–424 (2007)
  * Crisfield and Moita [1996] Crisfield, M., Moita, G.: A co-rotational formulation for 2-d continua including incompatible modes. International Journal for Numerical Methods in Engineering 39(15), 2619–2633 (1996)
  * Ferrari et al. [2021] Ferrari, F., Sigmund, O., Guest, J.K.: Topology optimization with linearized buckling criteria in 250 lines of matlab. Structural and Multidisciplinary Optimization 63(6), 3045–3066 (2021)
  * Svanberg [2002] Svanberg, K.: A class of globally convergent optimization methods based on conservative convex separable approximations. SIAM journal on optimization 12(2), 555–573 (2002)