Research

<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

Check for
updates

Article submitted to journal

**Subject Areas:**

mechanical engineering, applied
mathematics, computational
mechanics

**Keywords:**

harmonic balance method, bifurcation
analysis, stability, shape optimization,
nonlinear vibrations

**Author for correspondence:**

A. Mélot

e-mail: adrien.melot@inria.fr

# Multi-parametric optimization
for controlling bifurcation
structures

A. Mélot<fnref n="1" />, E. Denimal<fnref n="1" /> and L. Renson<fnref n="2" />

Bifurcations organise the dynamics of many natural and engineered systems. They induce qualitative and quantitative changes to a system's dynamics, which can have catastrophic consequences if ignored during design. In this paper, we propose a general computational method to control the local bifurcations of dynamical systems by optimizing design parameters. We define an objective functional that enforces the appearance of local bifurcation points at targeted locations or even encourages their disappearance. The methodology is an efficient alternative to bifurcation tracking techniques capable of handling many design parameters (> 10²). The method is demonstrated on a Duffing oscillator featuring a hardening cubic nonlinearity and an autonomous van der Pol-Duffing oscillator coupled to a nonlinear tuned vibration absorber. The finite element model of a clamped-free Euler-Bernoulli beam, coupled with a reduced-order modelling technique, is also used to show the extension to the shape optimization of more complicated structures. Results demonstrate that several local bifurcations of various types can be handled simultaneously by the bifurcation control framework, with both parameter and state target values.

<footnote n="1">¹Univ. Gustave Eiffel, Inria, COSYS-SII, I4S, Campus Beaulieu, 35042 Rennes, France</footnote>

<footnote n="2">²Dynamics Group, Department of Mechanical Engineering, Imperial College London, London, United Kingdom</footnote>
# 1. Introduction

Bifurcation phenomena are ubiquitous in natural and engineered systems. They play a key role in organising qualitative and quantitative changes to a system's dynamics, defining regions of multi-stability, or transitions from quiescent responses to oscillatory behaviours. For example, in biology, bifurcations underpin fundamental cellular processes like cell differentiation and spiking behaviours [1–3]; in chemical engineering, they serve as a powerful tool to identify the optimal operating conditions of chemical reactors [4]; in laser physics, they can be exploited to operate bi-stable lasers, which find applications as optical flip-flop memory apparatuses [5], or to produce pulsed beams via passive Q-switching [6,7]. In mechanics, the multi-stability induced by the presence of bifurcations has been exploited for mass detection [8,9], energy harvesting [10,11], shape morphing [12], and energy transfer [13]. The resonances of lightly damped nonlinear systems are also found near saddle-node bifurcations [14–16]. If unaccounted for, bifurcations can lead to undesirable transitions to quantitatively different regimes of motion, resulting in high stresses that can degrade components' life or even jeopardize the overall structure integrity [17–19].

In this paper, we propose a general computational methodology to control, through multi-parametric optimization, the local bifurcations of engineered systems, moving them at desired locations in parameter space or even removing them if necessary. The method is inspired by the work of Szep et al [20] where an optimization problem was formulated to estimate normal form model parameters of a gene model based on bifurcation diagrams. Here, we repurpose this work to the context of bifurcation control and generalise it to handle multiple types of bifurcations simultaneously. We develop a wider range of optimisation objectives that are more sophisticated than bifurcation parameter values at the bifurcation points and demonstrate the methodology can handle a relatively large number of design variables (> $10^2$).

Numerical continuation is a popular approach for conducting bifurcation analysis [21]. It is a predictor-corrector method that follows paths of equilibria or periodic solutions in, typically, a codimension-1 parameter space. Bifurcations along the solution curves can be detected, and their evolution can, in turn, be followed when a second parameter is permitted to vary (codimension-2) [22]. Bifurcation tracking methods have recently gained popularity within the mechanical vibration community as they enable the identification of stability boundaries and safe operating conditions [23–26]. Bifurcation tracking is arguably the simplest form of bifurcation optimisation in one design parameter. The first general methodology to optimise bifurcations was proposed by Doedel in [27,28] where the location of a fold bifurcation was controlled using a series of single-parameter continuations finding the extrema of the objective function with respect to one parameter at a time. This approach becomes tedious when a large number of parameters are considered and it necessitates repeating the procedure for each bifurcation point encountered on the solution branch. This introduces two notable challenges. Firstly, bifurcation tracking analyses are computationally more demanding compared to conventional nonlinear analyses, although some of the added computational burden could be alleviated by tracking bifurcation points in parallel. Secondly, several bifurcation curves may be disconnected from one another under variation of a design parameter [25]. This renders approaches based on bifurcation tracking less efficient. Overall, bifurcation tracking appears most suitable when one considers the control of a single bifurcation with respect to a single design variable. The methodology proposed here differs in that it only requires continuation with respect to the parameter(s) used to define the desired bifurcation location and does not require any continuation in all the other design parameters. Thus, it is well adapted for controlling all local bifurcations occurring in a specific bifurcation parameter interval with respect to a large number of design variables.
More recently, Boullé et al [29] formulated a shape optimization problem constrained by an augmented system of equations, similar to the one used for tracking fold bifurcations and branch points in [23,30], to control the location of a single branch point. This approach was extended to Hopf bifurcations [31] and successfully applied to a wide range of problems, including the complex Ginzburg-Landau and Navier-Stokes equations and a hyperelastic beam. The approach proposed here does not rely on bifurcation-defining constraints being verified at each iteration of the optimizer, giving more freedom to explore the design space.

Besides directly controlling bifurcations, numerous contributions investigated the optimization of nonlinear mechanical systems to obtain desired nonlinear behaviours. Pedersen [32] used topology optimization to mitigate the appearance of internal resonances in plates made of two materials. The methodology was able to handle several internal resonances simultaneously by expressing the objective function with a bound formulation, i.e. by maximizing the lowest of all combinations of frequencies leading to internal resonance. Dou et al [33] coupled the incremental harmonic balance method (HBM) and the method of moving asymptotes to optimize the non-uniform width of a clamped-clamped beam with constant thickness to minimize the amplitude of both primary and super-harmonic resonances. In [34], they successfully optimized coefficients defined with normal forms by employing the method of moving asymptotes by varying the element-wise thickness of both a clamped-clamped beam and a T-shaped structure. Later, they employed nonlinear normal modes computed with the HBM to optimize the hardening/softening behaviour of planar frame structures by varying the in-plane element-wise thickness of the constitutive elements [35]. The methodology developed in [34] was applied to the design of micro-electromechanical systems in [36]. More recently, Detroux et al [37] optimized the form of the nonlinearity to tailor the frequency-amplitude dependence of the resonances of nonlinear mechanical systems and suppress modal interactions and impose mode isochronicity with a 2-degree of freedom model. Grenat et al [38] implemented a recursive algorithm to compute the extrema of bifurcation curves in multi-parameter space for the design of a nonlinear tuned vibration absorber (NLTVA). Denimal et al [39,40] used topology optimization on under-platform dampers to reduce the amplitude at the resonance of a bladed disk. While all these methods are undoubtedly useful, they do not directly address bifurcations or only very specific types of bifurcations.

The paper is structured as follows. The optimization problem and the computational method for enforcing codimension-1 bifurcation points to appear at targeted locations are formulated in Section 2. The standard methods used to compute the solution curves and carry out bifurcation analysis are presented in Section 3. The method chosen in this paper are well-adapted for the mechanical applications considered here. However, the optimisation methodology is very general and other methods can be used. The capabilities of the proposed bifurcation control framework are demonstrated in Section 4 using three examples of increasing complexity, namely a forced Duffing oscillator, an autonomous van der Pol-Duffing (vdPD) oscillator coupled to a NLTVA [41], and the finite element model of a clamped-free Euler-Bernoulli beam with a nonlinear mechanism attached at its free end. The main conclusions are drawn in Section 5.

## 2. Formulation of the optimization problem

The aim of the proposed method is to determine a set of design parameters allowing a number of local bifurcation points to appear at targeted locations in the bifurcation diagram. The method relies on a series codimension-1 continuation analyses with respect to the bifurcation parameter in order to evaluate an objective functional whose purpose is to minimize the distance between the local bifurcations detected on the computed solution branch and a number of optimization targets. The objective functional is formulated in such a way that multiple local bifurcations of different types can be controlled simultaneously. This is achieved by systematically comparing all bifurcations of a specific type across all targets of the same type.
The optimization problem takes the general form

$$ 
\begin{array}{ll}
\underset{\boldsymbol{\theta}}{\text{minimize}} & f_o \\
\text{subject to} & b_i^l \le \theta_i \le b_i^u \quad \forall i \in [1, p],
\end{array}
\qquad (2.1) $$

where $f_o$ is the objective function to be minimized, $\theta \in \mathbb{R}^p$ is the vector of design (or decision) variables, e.g. the structural parameters of a mechanical system. The only constraints considered hereafter are box constraints defined with vectors $b^l$ and $b^u$ to ensure that the design parameters assume feasible values, e.g. positive mass. Any component $b_i^l$ or $b_i^u$ can be set to infinity if the parameter is only restricted to $\mathbb{R}^{+*}$ or $\mathbb{R}^{-*}$.

The proposed objective function is formulated as

$$ f_o(\mathcal{T}, \mathcal{P}) = \underbrace{|\mathcal{T} - \mathcal{P}|\Psi(\boldsymbol{\theta})}_{\text{Bifurcation measure}} + \underbrace{\frac{1}{|\mathcal{T}|} \sum_{\tau \in \mathcal{T}} \prod_{\pi(\boldsymbol{\theta}) \in \mathcal{P}} \left| \frac{\pi(\boldsymbol{\theta}) - \tau}{\tau} \right|^{1/|\mathcal{P}|}}_{\text{Error measure}}, \qquad (2.2) $$

where $\mathcal{T}$ represents the set of target bifurcations and $\mathcal{P}$ is the set of bifurcations predicted by the model at each step of the optimisation solver. The notation $|\mathcal{T}|$ simply represents the number of elements in the set $\mathcal{T}$.

The cost function comprises two types of terms that are hereafter referred to as the *bifurcation measure* and the *error measure*. They are further detailed in Sections 2(b) and 2(c), respectively. The computation of those terms for particular design parameters $\theta$ relies on a single one-parameter continuation presented in Section 2(a).

### (a) One-parameter continuation

At a fundamental level, numerical continuation is a method to trace out the solution manifold of a parameterised zero problem

$$ \mathbf{R}(\mathbf{q}, \boldsymbol{\mu}, \boldsymbol{\theta}) = \mathbf{0}, \qquad (2.3) $$

where $\mathbf{q} \in \mathbb{R}^n$ is the solution vector and $\boldsymbol{\mu}$ is the vector of bifurcation parameters. A wide range of responses can be captured by such a zero problem (equilibria, periodic responses, quasi-periodic responses, etc.). While the design parameters $\boldsymbol{\theta}$ could be considered as bifurcation parameters, they are here kept separate from $\boldsymbol{\mu}$ to emphasise the fact that no continuation is performed with respect to the design parameters. Furthermore, in this study, the bifurcation parameters are limited to a single parameter $\mu$, and the solution manifold is a one-dimensional curve. For the mechanical systems considered in Section 4, a typical choice of parameter is the applied excitation frequency, which gives the so-called forced (or frequency) response curves (FRC).

A bifurcation occurs at a specific point along the solution curve when the Jacobian of Eq. (2.3) becomes singular. Bifurcation points are usually detected by monitoring the sign of a scalar test function, hereafter denoted by $g$, whose zero indicates a bifurcation of a specific type. There are several ways to define such functions. In the present work, the approach based on minimally-extended systems is used (see Section 3). Although bifurcations are detected through changes in the parameter $\mu$, they also depend on the design variables $\boldsymbol{\theta}$.

### (b) Bifurcation measure

The bifurcation measure term incites the optimizer to explore regions of the optimization parameter space where the system exhibits multiple bifurcations of a given type. It circumvents the use of constraints that can be difficult to enforce during the optimization iterations, especially when the constraints are complicated functions, such as the Moore-Spence and Griewank and Riedden systems used in [29,31], that are computationally expensive to evaluate. With the present approach, the presence of bifurcations is not imposed but rather encouraged in a weak sense. This also facilitates the treatment of cases where bifurcations vanish during the optimization iterations.
The bifurcation measure is made of two terms. The first, $|T - P|$, serves as a penalization term when too few or too many bifurcations are detected and ensures the bifurcation measure disappears when the number of detected bifurcations exactly matches the number of targets. The second term, $\Psi(\theta)$, ought to tend to zero as the number of bifurcations increases and, for numerical reasons, be bounded between 0 and 1. Thus, it is desirable to formulate $\Psi(\theta)$ using the scalar test function $g$ defined in Section 3(b) as it changes sign every time a bifurcation occurs and is specific to a given bifurcation type. Therefore, the proposed, normalised bifurcation measure is:

$$G(\theta) = \frac{|g(\theta)|}{\max_{R=0} |g(\theta)|} \qquad (2.4)$$

In order to account for bifurcations occurring over a given bifurcation parameter interval, function $G(\theta)$ is integrated over the whole solution curve to define the global bifurcation measure. Since this integral depends on the path of the curve, defined by $R = 0$, it is also normalized:

$$\Psi(\theta) = \frac{\int_{R=0} G(\theta) ds}{\int_{R=0} ds} \qquad (2.5)$$

The bifurcation measure can be used independently of the error measure, i.e. without the need for imposing target locations. This can be used to penalise the appearance or presence of certain bifurcation types. For instance, in Section 4(b), it is used to completely suppress the presence of all local bifurcations of periodic orbits.

### (c) Error measure

Controlling the location of bifurcation points requires the evaluation of an error term to quantify how far detected bifurcations lie from the target locations. This constitutes a challenging task, especially when large variations of the structural parameters are allowed, as the number of bifurcations detected in the considered parameter range may change and even equal zero.

This issue is here circumvented by using a combination of geometric and arithmetic means as introduced in [20]. It allows one to consider all possible combinations of targets and detected bifurcations, with the added benefit of mitigating the risk of two bifurcations being matched to the same target as the term goes to zero only when all targets are matched with at least one bifurcation. The error measure in Eq. (2.2), can be broken down into two types of contributions as:

$$\text{error measure} = \frac{1}{|T|} \sum_{\tau \in T} \left( \prod_{\pi(\theta) \in P} \left| \frac{\pi_{\mu}(\theta) - \tau_{\mu}}{\tau_{\mu}} \right|^{1/|P|} + \prod_{\pi(\theta) \in P} \left| \frac{\pi_{\alpha}(\theta) - \tau_{\alpha}}{\tau_{\alpha}} \right|^{1/|P|} \right), \quad (2.6)$$

where the subscripts $\mu$ and $\alpha$ denote the terms associated with the bifurcation parameter and the response amplitude, respectively. Here, the term *response amplitude* refers to a measure of the system states which, depending on the application, can be considered as the infinity norm, the $L^2$ norm, etc.

### (d) Extension to multiple bifurcation types

The above objective function handles only a single type of bifurcations. However, it is common for a nonlinear forced response curve to exhibit several types of bifurcations in a given parameter range. Such a scenario can be handled with the present methodology by treating it as a regularized multi-objective optimization problem via the weighted sum method [42]. The optimization
problem for multiple bifurcation types can thus be formulated as:

$$ 
\begin{array}{ll}
\underset{\theta}{\operatorname{minimize}} & \alpha_{BP} f_o(T_{BP}, P_{BP}) + \alpha_{F} f_o(T_{F}, P_{F}) + \alpha_{PD} f_o(T_{PD}, P_{PD}) + \alpha_{NS} f_o(T_{NS}, P_{NS}) \\
\text { subject to } & b_i^l \le \theta_i \le b_i^u \quad \forall i \in [1, p] \\
&
\end{array}
\qquad (2.7)
$$

where $f_o$ is given by Eq. (2.2) and the subscripts BP, F, PD and NS denote branch points, fold, period doubling and Neimark-Sacker bifurcations, respectively. The extension of the cost function (2.7) to other types of bifurcations, such as bifurcations of equilibria (see Section 4(b)) and bifurcations of higher co-dimension (which would require to consider more bifurcation parameters), is straightforward.

### (e) Optimization algorithm

The proposed objective function is piece-wise continuous and non-convex. For this reason, we resort to gradient-free optimization algorithms. More specifically, this study was carried out using the Julia implementation of the NLopt nonlinear optimization package [43] since it offers numerous subroutines to choose from with support for bound constraints. Two gradient-free algorithms, namely COBYLA (Constrained Optimization BY Linear Approximation) [44,45] and Nelder-Mead [46–48] were employed. The latter was used for the computations in section 4(a) and was replaced by the COBYLA algorithm for the optimization of the beam model of section 4(c) as it was found to exhibit superior performances in higher dimensions.

## 3. Nonlinear analysis

### (a) Computation of the nonlinear forced response curves

A wide range of methods and software tools have been developed to compute nonlinear responses, study their stability and identify the presence of bifurcations [49–51]. For instance, time-domain methods such as the orthogonal collocation [49–52] or shooting methods [21] are popular methods to determine periodic responses. A powerful alternative considered here is the HBM [24,53–56]. Formulated in the frequency domain, the HBM is known to outperform time domain techniques when solving large-scale systems with strong nonlinearities [57] and has been employed in numerous fields of engineering and the applied sciences, including electrical engineering [58,59], aeroelasticity [60,61], fluid dynamics [62], chemical engineering [63] and biology [64].

The equation of motion of a wide range of n degrees-of-freedom mechanical models can be written as

$$ M\ddot{q}(t) + C\dot{q}(t) + Kq(t) + f_{nl}(q, \dot{q}, t) = f_{ex}(t), \qquad (3.1) $$

where M, C, K are respectively the mass, damping and stiffness matrices and $q(t)$ is the vector of unknown displacements. Nonlinear forces and external periodic forcing are described by vectors $f_{nl}$ and $f_{ex}$, respectively. We emphasize that the HBM can be applied to mathematical models expressed as first-order systems of equations, and the general optimization framework is independent of the method used to discretize periodic solutions. Other methods, such as orthogonal collocation, may be more suitable for problems stemming from other physics or with non-smooth nonlinearities. The only requirement for the continuation algorithm is that the set of equations be formulated as a zero problem (cf. Sect. 2(a)).

The HBM consists in expanding periodic solutions $q(t)$ of equation (3.1) as a truncated Fourier series of order H as

$$ q = \Re \left( \sum_{k=0}^{\infty} \tilde{q}_k e^{ik\Omega t} \right) \approx \Re \left( \sum_{k=0}^{H} \tilde{q}_k e^{ik\Omega t} \right), \qquad (3.2) $$

where $\tilde{q}$ contains the coefficients of the one-sided Fourier transform. In this paper, the response frequency, $\Omega$, coincides with the excitation frequency and corresponds to the bifurcation
<figure><img src="image_3.png" /><figcaption>Figure 1: Alternating frequency/time procedure using the fast Fourier transform and its inverse.</figcaption></figure>

parameter. The displacements, velocities and accelerations may be expressed as

$$ \mathbf{q} = \mathfrak{R}([\mathbf{T} \otimes \mathbf{I}_n] \tilde{\mathbf{q}}), \quad (3.3) $$

$$ \dot{\mathbf{q}} = \mathfrak{R} ([\mathbf{T}\Omega\nabla \otimes \mathbf{I}_n] \tilde{\mathbf{q}}), \quad (3.4) $$

$$ \ddot{\mathbf{q}} = \mathfrak{R} \left( \left[ \mathbf{T}\Omega^2 \nabla^2 \otimes \mathbf{I}_n \right] \tilde{\mathbf{q}} \right), \quad (3.5) $$

where $\otimes$ is the Kronecker product, $\mathbf{I}_n$ the identity matrix of size $n \times n$, $\mathbf{T}$ the vector of orthogonal exponential functions and $\nabla$ the complex frequency domain differential operator

$$ \nabla = \mathbf{diag}(0, i, 2i, \dots, Hi). \quad (3.6) $$

Similarly, the nonlinear forces and external forcing are written as

$$ \mathbf{f}_{nl} = \mathfrak{R} \left( \sum_{k=0}^{H} \tilde{\mathbf{f}}_{nl,k} e^{ik\Omega t} \right), \quad (3.7) $$

$$ \mathbf{f}_{ex} = \mathfrak{R} \left( \sum_{k=0}^{H} \tilde{\mathbf{f}}_{ex,k} e^{ik\Omega t} \right). \quad (3.8) $$

These expressions can be inserted into Eq. (3.1), and the time variable is eliminated by a Galerkin projection on the Fourier basis $\mathbf{T}$. This yields a residual $\mathbf{R}$ consisting of a set of $n(H + 1)$ nonlinear algebraic equations:

$$ \mathbf{R}(\tilde{\mathbf{q}}, \Omega) = \mathbf{Z}(\Omega)\tilde{\mathbf{q}} + \tilde{\mathbf{f}}_{nl}(\tilde{\mathbf{q}}) - \tilde{\mathbf{f}}_{ex} = \mathbf{0}, \quad (3.9) $$

where $\tilde{\mathbf{f}}_{nl}$ and $\tilde{\mathbf{f}}_{ex}$ contain the Fourier coefficients of the nonlinear forces and external forcing, respectively, and $\mathbf{Z}(\Omega)$ is the frequency domain dynamic stiffness matrix:

$$ \mathbf{Z}(\Omega) = \Omega^2 \nabla^2 \otimes \mathbf{M} + \Omega \nabla \otimes \mathbf{C} + \mathbf{I}_{H+1} \otimes \mathbf{K}. \quad (3.10) $$

To evaluate the nonlinear forces in the frequency domain, we resort to taking advantage of the existence of a closed-form expression of those forces in the time domain through the use of the alternating frequency/time procedure (AFT) as introduced in [65] and illustrated in Fig. 1. Note that the complex-valued residual is not necessarily holomorphic. It should therefore be converted to real and imaginary parts before calls to the solver.

The evolution of the solution with respect to a parameter (here, the excitation frequency) is followed by means of an arc-length continuation algorithm in conjunction with a tangent predictor. An extended residual $\mathbf{R}_{ex}$ is constructed by appending a constraint equation to account for the additional unknown and solved using a Newton-like solver. A thorough description of such a continuation method is omitted here for brevity. The interested reader may refer to [21] for additional details.

## (b) Bifurcation analysis
(i) Local stability analysis

The first step to detecting bifurcations consists in assessing the stability of the computed solutions during numerical continuation (Sect. 3(a)). This can be achieved with either the Floquet theory [66] or Hill's method [67–69]. Hill's method is employed here as it is particularly well suited to solutions computed with the HBM [67,70]. The main underlying idea is to introduce a perturbation of the form:

$$q_0(t) = q(t) + p(t)e^{At}. \qquad (3.11)$$

were $q(t)$ is a periodic solution of Eq. (3.1). This expression may be expanded as Fourier series with the Fourier basis $T$ used in the HBM formulation and introduced into the equations of motion. Eliminating the time variable with a Galerkin procedure and performing a first-order Taylor expansion on the nonlinear terms yields a quadratic eigenvalue problem:

$$[\Lambda^2 \tilde{M} + \Lambda \tilde{C} + \partial_q R] \tilde{p} = 0, \qquad (3.12)$$

with

$$\tilde{M} = I_{2H+1} \otimes M, \qquad (3.13)$$

$$\tilde{C} = \hat{\nabla} \otimes 2M + I_{2H+1} \otimes C, \qquad (3.14)$$

where $\partial_q R$ corresponds to the Jacobian of the HBM residual (3.9), which can be efficiently retrieved at each iteration of the continuation solver without additional computations, and $\hat{\nabla}$ is the real-valued differential operator in the frequency domain. The above equation can either be solved directly or linearized in order to alleviate the computational burden [71]:

$$\begin{bmatrix} \tilde{C} & \partial_{\tilde{q}} R \\ -I_{n(2H+1)} & 0 \end{bmatrix} + \lambda \begin{bmatrix} \tilde{M} & 0 \\ 0 & I_{n(2H+1)} \end{bmatrix} = 0. \qquad (3.15)$$

Solving this eigenvalue problem yields $2n(2H + 1)$ eigenvalues and eigenvectors. Their accuracy depends on the truncation order of the Fourier series — with higher orders being more accurate. Out of the complete set of solutions, $4nH$ are redundant due to the multiple harmonics used in the HBM approximation. An eigenvalue sorting algorithm is used to retain only the $2n$ eigenvalues with the smallest imaginary part. Monitoring the real part of all retained eigenvalues allows one to determine whether the solution is stable (all real parts are negative) or unstable (at least one real part is strictly positive).

(ii) Bifurcation detection

The bifurcation involved in the loss of stability of a periodic solution can be identified by studying the trajectory of the real parts crossing the imaginary axis. Figure 2 summarizes the three possible scenarios.

However, instead of monitoring the real part of the eigenvalues, it is common to define a scalar test function, $g$, that is equal to zero at specific bifurcation points. In this paper, this is achieved using a matrix bordering technique. The value of the test function can be evaluated by solving the following system:

$$\begin{bmatrix} \mathbf{A} & \mathbf{b} \\ \mathbf{d}^{\dagger} & 0 \end{bmatrix} \begin{pmatrix} \mathbf{w} \\ g \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix}, \qquad (3.16)$$

were $\mathbf{b}$ and $\mathbf{d}$ are chosen so that the system is nonsingular and $\dagger$ represents the hermitian transpose. The matrix $\mathbf{A}$ takes different values depending on the type of bifurcation to detect. For fold bifurcations, $\mathbf{A} = \partial_q \mathbf{R}$ is used. For branch points:

$$\mathbf{A} = \begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \mathbf{z}^T & \mathbf{0} \end{bmatrix}, \qquad (3.17)$$

where $\mathbf{z}$ is the tangent predictor used during the continuation analysis. Neimark-Sacker and period-doubling bifurcations are both detected using $\mathbf{A} = 2\mathbf{V} \odot \mathbf{I}$ where $\mathbf{V} = \text{Diag}(\lambda)$ is the
<figure><img src="image_4.png" /><figcaption>Figure 2: Bifurcation scenarios depending on the trajectory of the Floquet exponents. Fold bifurcation and branch point (a), period doubling bifurcation (b) and Neimark-Sacker bifurcation (c).</figcaption></figure>

diagonal matrix gathering the eigenvalues of Eq. (3.15) and $\odot$ is the bialternate matrix product [72]. Both bifurcation types can then be distinguished based on the imaginary part of the eigenvalues whose real part becomes positive.

## 4. Application of the proposed methodology

The capabilities of the proposed optimization framework are now illustrated in three examples. The first is a Duffing oscillator, which serves to show that different target types can be imposed and that bifurcations of different types can be handled simultaneously. The second example is a self-excited system comprised of a van der Pol oscillator and an NLTVA, and shows how the bifurcation measure (Section 2(b)) can be exploited to remove all local bifurcations of periodic orbits. Lastly, the methodology is applied to the finite element model of a cantilever beam with stiffness nonlinearity at its tip, which comprises > $10^2$ design parameters.

### (a) Duffing oscillator

The equation of motion of the Duffing oscillator is

$$m\ddot{q} + c\dot{q} + kq + k_{nl}q^3 = f_{ex} \quad (4.1)$$

where the initial parameter values for the mass, damping, linear and nonlinear stiffnesses are $m = 1$, $c = 0.1$, $k = 1$, and $k_{nl} = 0.5$, respectively. The external excitation is given by $f_{ex} = 3 \cos(\Omega t)$. Five harmonics and the mean value are considered in the analysis and the displacements and nonlinear forces are sampled with $N_s = 2^8$ points in the AFT procedure. Figure 3a shows that the system exhibits a hardening behaviour, as is expected since the nonlinear coefficient $k_{nl}$ is strictly positive. Eight bifurcations are present on the main solution branch — six are fold bifurcations (●) associated, from low to high frequencies, to the 5:1 and 3:1 superharmonic resonances and to the primary resonance. The other two bifurcations are branch points (●) associated with the 2:1 superharmonic resonance [73].

The various capabilities of the optimization framework are illustrated by setting targets on bifurcation frequencies, bifurcation amplitudes and both bifurcation frequencies and amplitudes simultaneously. The optimization parameters considered are $\theta = (m, c, k, k_{nl})^T$ and only $f_{ex}$ is kept fixed. Note that, for the Duffing oscillator, the number of independent model parameters can be reduced to two through rescaling. Here, the original four parameters are retained to show the good behaviour of the methodology in the presence of multiple optimum solutions arising from the over-parameterisation of the problem, which is a common occurrence in practice. The Nelder-Mead algorithm of the Julia implementation of the NLopt package is used. No stopping criterion other than the maximum number of iterations, set to 200, is defined. The optimization parameters
are allowed to take values within the interval $[b_i^l, \infty[$ where $b_i^l, i \in [1, 4]$ are the components of vector $\mathbf{b}^l = (0.01, 0.01, 0.01, 0.01)^T$.

Starting from the initial configuration, a first optimization run is carried out to enforce the appearance of two fold bifurcations at angular frequencies $\Omega = 1.59$ rad/s and $\Omega = 3$ rad/s, respectively. Figure 3a shows the location of the detected bifurcations with respect to the target frequencies for the initial parameter configuration. Figure 3b shows the nonlinear FRC of the oscillator with optimized parameters values $m = 1.85$, $c = 0.17$, $k = 1.77$ and $k_{nl} = 0.54$. The value of the objective function is $f_o = 2.0 \cdot 10^{-4}$. The computation of the dimensionless parameters $\epsilon = k_{nl}f_{ex}^2/k^3$ and $\delta = c/2\sqrt{k}$ reveals that the system with optimized parameters is significantly less nonlinear ($\epsilon = 0.9$) than the initial system ($\epsilon = 4.5$) and a little bit more damped ($\delta = 6.3 \cdot 10^{-2}$ compared to $\delta = 5.0 \cdot 10^{-2}$). One can see that the nonlinear FRC exhibits only two fold bifurcations, corresponding to the number of targets. Although superharmonic resonances are still visible, they no longer exhibit fold bifurcations which is consistent with the change in dimensionless parameters and suppress multi-stability and the risk of amplitude jumps. Furthermore, the angular frequencies at which the folds appear show excellent agreement with the imposed targets. This is associated with a shift of the resonance of the underlying linear system towards the lower frequencies and a reduction of approximately 13% of the resonance amplitude. Note that the amplitude reduction at resonance is expected when the interval between

<figure><img src="image_5.png" /><figcaption>Figure 3: Optimization of the structural parameters of the Duffing oscillator with respect to the frequencies (a, b) and amplitudes (c, d) at which fold bifurcations appear. Before optimization (a, c) and after optimization (b, d). Solid and dashed blue lines represent stable and unstable responses. Fold bifurcations and branch points are denoted by red (●) and orange (○) circle markers, respectively. Dashed black lines represent optimization targets.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the objective function with the number of iterations with the frequency of two folds as targets (a) leading to Fig. 3b and the amplitude of two folds as target (b) leading to Fig. 3d.</figcaption></figure>

two successive fold bifurcations decreases because of the frequency-amplitude dependence of the nonlinear modes. We have not identified any particular challenges when choosing targets resulting in an increase in oscillation amplitude.

A similar optimization may be carried out by imposing two fold bifurcations at targeted amplitudes equal to $\tilde{q}_1 = 1.5$ and $\tilde{q}_1 = 5$, respectively. Again, the FRC of the optimized system ($m = 1.29, c = 0.18, k = 1.60$ and $k_{nl} = 0.61$), shown in Fig. 3d, only exhibits two fold bifurcations appearing at the targeted amplitudes of the first harmonic, corresponding to a reduction of 26% of the resonance amplitude compared to the initial system (Fig. 3c). Besides, contrary to the optimization with targeted frequencies, the optimized parameters do not lead to a significant shift of the resonance frequency of the linear underlying system.

Figure 4 depicts the evolution of the objective function with the number of iterations of the optimization algorithm when the frequencies (Fig. 4a) and amplitudes (Fig. 4b) of two fold bifurcations are defined as targets. It appears that the objective function exhibits an overall rapid, although not monotonic, decrease with the number of iterations and bifurcations are matched quite accurately after a few iterations with the objective function reaching values smaller than 0.1. Note that the objective function assumes larger values when the bifurcations are matched quite accurately but their number exceeds the targets (cf. the FRCs shown in Fig. 4a). The oscillations cease in both cases after approximately 100 iterations of the optimization algorithm with the converged solutions showing only marginal improvement. The optimization with 4 variables required approximately 5 minutes of computational time on a standard laptop (i7-12700H @ 2.7 GHz, 16 Gb RAM) running Linux and Julia 1.8.

Using Eq. (2.7), the proposed optimization framework can also handle several types of bifurcation points simultaneously. Figure 5 shows the results of an optimization carried out with targets imposed on four fold bifurcations ($\Omega = 3.50$ rad/s, $\Omega = 2.00$ rad/s, $\Omega = 0.56$ rad/s, $\Omega = 0.54$ rad/s) and two branch points ($\Omega = 0.87$ rad/s, $\Omega = 0.80$ rad/s). One can see that the optimization leads to a smaller interval between two consecutive branch points while retaining a good match between predicted fold bifurcations and targets.

Figure 5 also illustrates the important observation that there is, in general, no complete freedom regarding the simultaneous optimization of multiple bifurcations. Indeed, the location of bifurcations is governed by physical phenomena, such as (super-/sub-harmonic) resonances or modal interactions, which occur at specific locations when particular ratios between the excitation frequency and the response's fundamental and higher-harmonic frequencies are realised. For instance, in the present example, it is not possible to obtain large variations of the frequencies at which the branch points occur while having active targets on other fold bifurcations since the
<figure><img src="image_7.png" /></figure>

Figure 5: Optimization of the structural parameters of the Duffing oscillator with respect to the frequencies at which fold bifurcations and branch points appear. Before optimization (a), after optimization (b). Solid and dashed blue lines represent the stable and unstable responses, respectively. Fold bifurcations and branch points are denoted by red (●) and orange (●) circle markers, respectively. Dashed coloured lines represent optimization targets.

<figure><img src="image_8.png" /></figure>

Figure 6: Optimization of the structural parameters of the Duffing oscillator with respect to the frequencies and amplitudes at which fold bifurcations appear. Before optimization (a), after optimization (b). Solid and dashed blue lines represent the stable and unstable responses, respectively. Fold bifurcations and branch points are denoted by red (●) and orange (●) circle markers, respectively. Dashed black lines represent optimization targets.

former are associated with even-order superharmonic resonances and the latter with odd-order superharmonic resonances. In this case, the optimisation algorithm finds a compromise between the different bifurcation targets.

Figure 6 shows that the proposed methodology can also simultaneously handle frequency and amplitude targets. An optimization was run with two pairs of frequency-amplitude targets to enforce fold bifurcations at $(\Omega, \tilde{q}_1) = (2, 1.5)$ and $(\Omega, \tilde{q}_1) = (3, 4.5)$, respectively. After optimization, the parameter values obtained are $m = 1.33$, $c = 0.22$, $k = 2.06$ and $k_{nl} = 0.63$. Once again, the optimized parameters lead to a significantly less nonlinear ($\epsilon = 0.65$) and more damped ($\delta = 7.7 \cdot 10^{-2}$) oscillator in order to reduce the number of bifurcations and the amplitude of the upper fold. The results show that only two bifurcations are visible on the bifurcation diagram of
the Duffing oscillator and that a very good agreement with the targets is achieved. The remaining error in the bottom fold re-emphasises the observation made in Fig. 5 that there is no complete freedom in setting the parameter and amplitude values for the two fold simultaneously.

(b) Autonomous van der Pol-Duffing oscillator coupled to a NLTVA

(i) Model description

The autonomous van der Pol oscillator is a well-known system exhibiting self-excited oscillations. It captures the dynamics of many systems in engineering and the applied sciences, such as aircraft wings [74], vehicle wheels [75], machine tools [76], and living cells [77,78]. A stable equilibrium point loses its stability through a Hopf bifurcation, leading to the appearance of limit cycle oscillations (LCO). The Hopf bifurcation can be either subcritical or supercritical, depending on the parameters. The former is dangerous, as stable limit cycles coexist with the stable equilibrium point. Various authors have attached nonlinear energy sinks, linear tuned vibration absorbers or NLTVA to delay the occurrence of the Hopf bifurcation, enforce the supercriticality of the Hopf bifurcation and reduce the amplitude of the LCOs. In the following, we use the mechanical system studied in [41], i.e. a van der Pol-Duffing oscillator coupled to a NLTVA with cubic nonlinearity. The equations of motion are

$$ 
\begin{align} 
m_1 \ddot{q}_1 + c_1 (q_1^2 - 1) \dot{q}_1 + k q_1 + k_{nl1} q_1^3 + c_2 (\dot{q}_1 - \dot{q}_2) + k_2 (q_1 - q_2) + k_{nl2} (q_1 - q_2)^3 = 0, \tag{4.2} \\ 
m_2 \ddot{q}_2 + c_2 (\dot{q}_2 - \dot{q}_1) + k_2 (q_2 - q_1) + k_{nl2} (q_2 - q_1)^3 = 0, 
\end{align} 
$$

where subscripts 1 and 2 denote quantities related to the vdPD oscillator and the NLTVA, respectively. For practical reasons, the mass ratio $\epsilon = m_2/m_1 = 0.05$ is kept fixed [41]. The mechanical properties of the vdPD oscillator are chosen as $m_1 = 1$, $k_1 = 1$ and $k_{nl1} = 0.03k_1$. The optimization parameters considered hereafter are the linear and nonlinear stiffness coefficients as well as the damping coefficient of the NLTVA, denoted respectively by $k_2$, $k_{nl2}$ and $c_2$ ($\theta = (k_2, k_{nl2}, c_2)^T$) and initially set at $k_2^0 = 4.625 \cdot 10^{-2}$, $k_{nl2}^0 = 6.8 \cdot 10^{-4}$ and $c_2^0 = 1.15 \cdot 10^{-2}$.

The continuation parameter is $c_1$ as in [41]. To apply continuation to this system, the procedure detailed in Section 3 must be modified to account for the autonomous nature of the system. To do so, a phase constraint $\dot{q}_1(t=0) = 0$ is added to the HBM residual equation (3.9). More details on the computation of LCOs with the HBM can be found in [79]. At each iteration of the optimizer, the Hopf bifurcation is located through linear stability analysis. The value of $c_1$ at which the Hopf bifurcation occurs, as well as the imaginary parts of the eigenvalues with positive real parts, are used as initial guess for the continuation algorithm. A small value of the Fourier coefficient corresponding to the first harmonic is used as the initial LCO amplitude. In the following, computations are carried out using 13 harmonics and $2^{10}$ samples in the AFT procedure.

(ii) Optimization of the NLTVA to suppress the fold-induced subcriticality

As illustrated in Fig. 7, the bifurcation diagram of the initial system described in Sect. 4(b)(i) exhibits a supercritical Hopf bifurcation at approximately $c_1 = 0.183$. The amplitude of the LCO increases with $c_1$ until the curve exhibits a first fold bifurcation (●), which causes the LCO to lose its stability and jump from an RMS amplitude equal to 0.27 to an amplitude equal to 1.0. More importantly, the second fold bifurcation occurs for a value of $c_1$ that is lower than that of the Hopf bifurcation ($c_1 = 0.175$). This results in the creation of a bi-stable region where a relatively high amplitude LCO coexists with the equilibrium point, which is often considered a dangerous situation in many engineering applications.

In order to try to suppress or mitigate this phenomenon, we optimize the NLTVA to remove any bifurcation appearing past the Hopf point on the bifurcation diagram. To this end, the targetted bifurcation measure is set to zero, i.e. $|T| = 0$, for all local bifurcation types of periodic orbits (folds, Neimark-Sacker, period doubling and branch points). No error measure
<figure><img src="image_9.png" /><figcaption>Figure 7: Bifurcation diagram of the van der Pol-Duffing oscillator coupled to the initial (purple) and optimized (blue) NLTVA. Solid and dashed lines indicate stable and unstable responses, respectively. Fold bifurcations are represented by red circle markers (●) and the Hopf bifurcation target is shown in dashed black line.</figcaption></figure>

is considered for these bifurcations. To preserve the location of the Hopf bifurcation with respect to the bifurcation parameter, a simplified error measure is used:

$$ \text{error measure} = \left| \frac{c_{1,\pi}^{\text{Hopf}} - c_{1,\tau}^{\text{Hopf}}}{c_{1,\tau}^{\text{Hopf}}} \right| \qquad (4.3) $$

where $c_{1,\pi}^{\text{Hopf}} = 0.183$ is the targeted Hopf location and $c_{1,\tau}^{\text{Hopf}}$ is the predicted location based on model parameters. No bifurcation measure is considered for the Hopf bifurcation since a unique Hopf bifurcation point is expected regardless of the chosen structural parameters.

Both a lower bound $\mathbf{b}^l = (0.5k_2^0, 0.1k_{nl2}^0, 0.8c_2^0)^T$ and upper bound $\mathbf{b}^u = (1.5k_2^0, 10k_{nl2}^0, 1.2c_2^0)^T$ are defined for the structural parameters. The initial guess for the optimization corresponds to the mechanical properties $\theta_0 = (k_2^0, k_{nl2}^0, c_2^0)^T$ given in Sect. 4(b)(i). The purpose of this example is to demonstrate that our methodology can be used to efficiently eliminate any local bifurcation of periodic orbit occurring within a specified parameter interval. Since no explicit target is defined for such bifurcations, the objective functional essentially behaves as a step function and is constant so long as the Hopf-parameter target is matched. Thus, the solution presented is but one within a family of possibilities, the size of which can be restrained by introducing additional design constraints.

The bifurcation diagram of the optimized system ($k_2 = 4.65 \cdot 10^{-2}$, $k_{nl2} = 6.8 \cdot 10^{-5}$ and $c_2 = 1.27 \cdot 10^{-2}$), plotted in blue in Fig. 7, shows that LCOs arise from a supercritical Hopf bifurcation at the same value of $c_1 = 0.183$ as the initial system. The fold-induced subcritical behaviour is effectively suppressed as no local bifurcation of periodic orbit appear on the diagram. Furthermore, although the LCO amplitude before the first fold bifurcation of the initial system is slightly higher, the amplitude of the stable branch is significantly smaller past that point.

### (c) Euler-Bernoulli beam with cubic nonlinearity

#### (i) Model description

The mechanical system considered in this section is a planar cantilever beam featuring a nonlinear attachment consisting of a cubic spring located at its free end (Fig. 8). The geometrical and mechanical properties (Table 1) are taken from the benchmark system studied in [15,80,81]. The initial model has a length $L = 0.7$ m, a height and a width of 0.014 m. The beam is discretized with
$N_e = 59$ Euler-Bernoulli beam elements with two degrees of freedom per node, i.e. one translation and one rotation. The cubic spring is attached to the translational degree of freedom. An external periodic forcing $F = 90 \sin(\Omega t)$ is applied at 0.21 m from the clamped end.

To reduce the computational effort, a reduced-order model is constructed using the Craig-Bampton method. The translational degree of freedom connecting with the cubic spring at the beam tip is kept as master degrees of freedom and 5 fixed-interface modes are kept in the reduction basis. A constant modal damping matrix is defined such that the damping ratio of the first linear mode is equal to 1% and that of all subsequent modes is equal to 5%.

<figure><img src="image_10.png" /><figcaption>Figure 8: Schematic representation of the Euler-Bernoulli beam.</figcaption></figure>

<table><caption>Table 1: Reference geometrical and mechanical properties of the considered beam [80].</caption><thead><tr><th>Length (m)</th><th>Width (m)</th><th>Height (m)</th></tr></thead><tbody><tr><td>0.7</td><td>0.014</td><td>0.014</td></tr><tr><td>Young's modulus (N/m<sup>2</sup>)</td><td>Density (kg/m<sup>3</sup>)</td><td>Nonlinear coeff. (N/m<sup>3</sup>)</td></tr><tr><td>2.05 · 10<sup>11</sup></td><td>7800</td><td>6.0 · 10<sup>9</sup></td></tr></tbody></table>

## (ii) Shape optimization

The structural parameters considered for the optimization are $\theta = (L, \mathbf{w}_e^T, \mathbf{h}_e^T, k_{nl})^T$ where $\mathbf{w}_e$ and $\mathbf{h}_e$ denote the vectors containing the width and height, $w_e^{(i)}$ and $h_e^{(i)}$, $i \in [1, N_e]$, of all the elements. The total number of optimization variables is $2N_e + 2 = 120$. We assume that any modification of the total length $L$ of the beam is evenly distributed over all elements. New FE and reduced-order models are generated at each iteration of the optimization algorithm and the error on the eigenfrequencies of the reduced and full-order models is monitored to ensure the convergence and accuracy of the results. Lower and upper parameter bounds are defined so that the length of the beam is allowed to vary between 0.4 m and 1.0 m, the element-wise width and height can exhibit variations as large as 40% of the reference values and the nonlinear coefficient is bounded between $6.0 \cdot 10^8$ and $6.0 \cdot 10^{10}$ N/m$^3$.

A preliminary reference computation is carried out with the parameters of Table 1. The nonlinear forced response (see Fig. 9a) is computed around the primary resonance corresponding to the first linear mode whose frequency is $f_1 = 23.3$ Hz. To account for changes in the eigenfrequencies of the underlying linear system due to the updated geometry and ensure that the whole resonance peak is computed, the continuation interval is expressed as a function of the first eigenfrequency $\omega_1$ of the underlying linear system, i.e. $[\omega_1, 4\omega_1]$. Note that the very same procedure could be applied to other modes. Besides, the continuation procedure is initiated at $\Omega = 4\omega_1$, i.e. far from resonance where the response amplitude is small to avoid potential convergence issues of the initial guess.

It is well known that this beam system possesses an isolated branch of solutions that merges with the main branch for sufficiently high excitation amplitudes [15]. Such a scenario, illustrated in Figure 9a, can be extremely detrimental to the structural integrity of the system as the
<figure><img src="image_11.png" /><figcaption>Figure 9: Forced response curve of the beam with initial (a) and optimized (b) geometries for $F = 90 \sin(\Omega t)$. Solid and dashed lines indicate stable and unstable responses, respectively. Fold and Neimark-Sacker bifurcations are denoted by red circle (●) and green triangle markers (▽), respectively. Dashed black lines represent the optimization targets. The FRC of the initial geometry with unmerged isola is shown in light grey (a) for $F = 80 \sin(\Omega t)$.</figcaption></figure>

maximum resonance frequency and amplitude can be significantly underestimated (here, by almost 50%). In the following, we run the optimization procedure to mitigate the frequency increase induced by the merged isola by defining two frequency targets for fold bifurcations at $\Omega = 47.7$ Hz and $\Omega = 71.6$ Hz. As can be seen in Figure 9b, the FRC of the optimized beam exhibits four folds (●), two of which are at the targeted locations. This number of folds exceeds the number of targets, which could be due either to the use of a local optimization algorithm or to inherent limitations in the modifications of the structure of the resonance. Interestingly, the shape of the resonance is altered (see the unstable region near the bottom fold). This indicates that the optimized design induced a shift of the forcing level at which the isola merges with the main solution branch.

The FRC in Fig. 9b also exhibits four NS bifurcations (▽) — two more than in Fig. 9a. The two new NS bifurcations occur outside the resonance peak at approximately $\Omega = 75$ Hz and $\Omega = 83$ Hz. Their appearance was made possible because NS bifurcations were not accounted for in the objective function. To eliminate these bifurcations, two additional target parameter values for NS bifurcations are set to $\Omega = 48.3$ Hz and $\Omega = 54.5$ Hz, i.e. values corresponding to the NS bifurcations located on the upper part of the FRC after the first optimization run. Since the number of detected NS bifurcations exceeds the number of targets, the bifurcation measure of the objective function pushes the optimization algorithm towards designs that only exhibit two NS bifurcations on the FRC. Although it is not possible to explicitly choose which bifurcation will be eliminated, as all prediction/target pairs of a given bifurcation type are considered, the elimination of the two NS bifurcations located on the lower part of the FRC is more likely for they are further from the targets. Figure 10 effectively shows that the FRC of the optimized system exhibits only two folds and two NS bifurcations, as intended.

A small error ($\approx 2$ Hz) on the location of the lower fold bifurcation (near $\Omega = 45$ Hz) can be observed. Similarly, the NS bifurcations do not perfectly match their target values. The reason for this discrepancy is similar to that for the branch points on the Duffing oscillator. Here, the two NS bifurcations cannot be arbitrarily moved in parameter space because they are associated with the 3:1 modal interaction that exists between the first and second bending modes of the beam [15,82,83]. It is therefore quite difficult to obtain significant changes in their locations when considering only the first resonance and constraining its corresponding fold bifurcations.
<figure><img src="image_12.png" /></figure>

Figure 10: Forced response curve of the beam with initial geometry (a) and geometry optimized for Neimark-Sacker bifurcations (b). Solid and dashed lines indicate stable and unstable responses, respectively. Fold and Neimark-Sacker bifurcations are denoted by red circle (●) and green triangle (▽) markers, respectively. Dashed red and green lines represent target locations for fold and Neimark-Sacker bifurcations, respectively.

<figure><img src="image_13.png" /></figure>

Figure 11: Geometry of the optimized beam. Side view (a) and top view (b). The outline of the reference geometry is displayed in red.
Following this example, it appears desirable to systematically consider all potential types of bifurcations in the objective function and exploit the bifurcation measure to avoid introducing additional unwanted bifurcations (as in Section 4(b)). However, this usually comes at an increased computational cost. Indeed, in the present beam example, accurate detection of the NS bifurcations requires the use of at least $H = 15$ harmonics in the stability analysis whereas only $H = 5$ harmonics are necessary to reliably detect the fold bifurcations. On this example, the computational time for 120 optimization variables with 5 harmonics for 200 iterations was approximately 2 h. The optimization of NS bifurcations required significantly more harmonics (23) to reach convergence on all continuation runs which led to a computational time of approximately 20 h.

The optimal design associated with this optimization is shown in Fig. 11. One can see that the optimized beam is longer than the reference configuration, which is coherent with the shift of the resonance towards lower frequencies (Fig. 10). The height of all elements is slightly reduced except for that of the 11th, 46th, 47th, 58th and 59th which is slightly increased. The top view of the beam reveals a more complex geometry. All the elements of the optimized beam are slightly wider than those of the reference configuration, with a few (elements 9, 35, 47, 50 53 58 and 59) showing a noticeable increase in width. Overall, most elements have the same width and length. Although this is most likely associated to the fact that the initial guess has uniform width and height, there is no clear explanation as to why this occurs as the local increases in width/height cannot be easily linked to a specific mechanical property that would offer leverage on bifurcation location. The obtained geometry is also likely to be a local optimum and other, potentially qualitatively different, designs could be obtained if a more extensive search of the design space was conducted with, for instance, a global optimization algorithm. The manufacturing of such a design can be achieved using, for instance, additive manufacturing methods. A similar approach was taken to make NLTVAs [84,85].

## 5. Conclusion

Bifurcations lead to significant changes in a system's dynamics, which can be exploited (switching, energy harvesting, etc.) to improve a system's performance, or result in degraded life and even catastrophic failures if unaccounted for. In this paper, a computational method to optimize the presence and location of local bifurcations in the response of a dynamical system was proposed and successfully demonstrated on three mechanical systems, including a finite element model with more than one hundred design parameters. The particularity of the objective function is to include a bifurcation error term, which promotes the right number and the right types of bifurcations. The objective function also includes an error measure to minimise the distance between the bifurcations and their desired location in terms of the state and bifurcation parameter spaces. The method was shown to simultaneously tackle multiple bifurcations of different types, including fold, branch points and Neimark-Sacker.

The proposed approach is also extremely general and a priori applicable to a wide range of problems in engineering and the applied sciences. Given the mechanical applications considered here, the harmonic balance method was used to discretise periodic solutions, the arclength continuation method was exploited to follow their evolution in the codimension-one parameter space, and Hill's stability analysis method with bordering techniques was exploited to detect and classify bifurcations along the solution curve. However, the proposed method does not rely on the use of these methods, and others could be considered if deemed more suitable for the application at hand.

**Acknowledgements.** This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.
References

1. Marco E, Karp RL, Guo G, Robson P, Hart AH, Trippa L, Yuan GC. 2014 Bifurcation analysis of single-cell gene expression data reveals epigenetic landscape. *Proceedings of the National Academy of Sciences* **111**, E5643–E5650. (10.1073/pnas.1408993111)

2. Gerstner W, Kistler WM. 2002 *Spiking Neuron Models: Single Neurons, Populations, Plasticity*. Cambridge University Press. (10.1017/CBO9780511815706)

3. Hay E, Hill S, Schürmann F, Markram H, Segev I. 2011 Models of Neocortical Layer 5b Pyramidal Cells Capturing a Wide Range of Dendritic and Perisomatic Active Properties. *PLOS Computational Biology* **7**, 1–18. (10.1371/journal.pcbi.1002107)

4. Guckenheimer J. 1986 Multiple bifurcation problems for chemical reactors. *Physica D: Nonlinear Phenomena* **20**, 1–20. (https://doi.org/10.1016/0167-2789(86)90093-X)

5. Liu L, Kumar R, Huybrechts K, Spuesens T, Roelkens G, Geluk EJ, de Vries T, Regreny P, Van Thourhout D, Baets R, Morthier G. 2010 An ultra-small, low-power, all-optical flip-flop memory on a silicon chip. *Nature Photonics* **4**, 182–187. (10.1038/nphoton.2009.268)

6. Arimondo E, Casagrande F, Lugiato LA, Glorieux P. 1983 Repetitive passive Q-switching and bistability in lasers with saturable absorbers. *Applied Physics B* **30**, 57–77. (10.1007/BF00688479)

7. Doedel EJ, Oldeman BE, Pando L. CL. 2011 Bifurcation structures in a model of a CO₂ laser with a fast saturable absorber. *International Journal of Bifurcation and Chaos* **21**, 305–322. (10.1142/S021812741102843X)

8. Kumar V, Boley JW, Yang Y, Ekowaluyo H, Miller JK, Chiu GTC, Rhoads JF. 2011 Bifurcation-based mass sensing using piezoelectrically-actuated microcantilevers. *Applied Physics Letters* **98**, 153510. (10.1063/1.3574920)

9. Nguyen VN, Baguet S, Lamarque CH, Dufour R. 2015 Bifurcation-based micro/nanoelectromechanical mass detection. *Nonlinear Dynamics* **79**, 647–662. (10.1007/s11071-014-1692-7)

10. Arrieta AF, Hagedorn P, Erturk A, Inman DJ. 2010 A piezoelectric bistable plate for nonlinear broadband energy harvesting. *Applied Physics Letters* **97**, 104102. (10.1063/1.3487780)

11. Abdelkefi A, Nayfeh AH, Hajj MR. 2012 Design of piezoaeroelastic energy harvesters. *Nonlinear Dynamics* **68**, 519–530. (10.1007/s11071-011-0233-x)

12. Kuder IK, Arrieta AF, Ermanni P. 2015 Design space of embeddable variable stiffness bi-stable elements for morphing applications. *Composite Structures* **122**, 445–455. (https://doi.org/10.1016/j.compstruct.2014.11.061)

13. Hurel G, Baguet S, Lamarque CH. 2022 Neimark Sacker bifurcations and non-linear energy exchange in chains of non-linear oscillators. *International Journal of Non-Linear Mechanics* **144**, 104057. (https://doi.org/10.1016/j.ijnonlinmec.2022.104057)

14. Hill T, Cammarano A, Neild S, Wagg D. 2015 Interpreting the forced responses of a two-degree-of-freedom nonlinear oscillator using backbone curves. *Journal of Sound and Vibration* **349**, 276–288. (https://doi.org/10.1016/j.jsv.2015.03.030)

15. Kuether R, Renson L, Detroux T, Grappasonni C, Kerschen G, Allen M. 2015 Nonlinear normal modes, modal interactions and isolated resonance curves. *Journal of Sound and Vibration* **351**, 299–310. (https://doi.org/10.1016/j.jsv.2015.04.035)

16. Cenedese M, Haller G. 2020 How do conservative backbone curves perturb into forced responses? A Melnikov function analysis. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* **476**, 20190494. (10.1098/rspa.2019.0494)

17. Rose J, Jinu G. 2016 Influence of aeroelastic flow-induced oscillations on fatigue life of an airplane wing structure. *International Journal of Modelling, Identification and Control* **25**, 199–216.

18. Nilsson KF, Thesken J, Sindelar P, Giannakopoulos A, Stoåkers B. 1993 A theoretical and experimental investigation of buckling induced delamination growth. *Journal of the Mechanics and Physics of Solids* **41**, 749–782.

19. Renson L, Noël J, Kerschen G. 2015 Complex dynamics of a nonlinear aerospace structure: numerical continuation and normal modes. *Nonlinear Dyn* **79**, 1293–1309.

20. Szep G, Dalchau N, Csikasz-Nagy A. 2021 Parameter Inference with Bifurcation Diagrams. Preprint. (10.48550/ARXIV.2106.04243)

21. Seydel R. 2010 *Practical Bifurcation and Stability Analysis*. Springer-Verlag New York. (https://doi.org/10.1007/978-1-4419-1740-9)

22. Y. A. K. 2023 *Elements of Applied Bifurcation Theory*. Springer. (https://doi.org/10.1007/978-3-031-22007-4)
23. Xie L, Baguet S, Prabel B, Dufour R. 2017 Bifurcation tracking by Harmonic Balance Method for performance tuning of nonlinear dynamical systems. *Mechanical Systems and Signal Processing* **88**, 445-461. (https://doi.org/10.1016/j.ymssp.2016.09.037)

24. Detroux T, Renson L, Masset L, Kerschen G. 2015 The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. *Computer Methods in Applied Mechanics and Engineering* **296**, 18-38. (https://doi.org/10.1016/j.cma.2015.07.017)

25. Mélot A, Rigaud E, Perret-Liaudet J. 2022 Bifurcation tracking of geared systems with parameter-dependent internal excitation. *Nonlinear Dynamics* **107**, 413-431. (10.1007/s11071-021-07018-6)

26. Mélot A, Rigaud E, Perret-Liaudet J. 2023 Robust design of vibro-impacting geared systems with uncertain tooth profile modifications via bifurcation tracking. *International Journal of Non-Linear Mechanics* **149**, 104336. (https://doi.org/10.1016/j.ijnonlinmec.2022.104336)

27. Doedel E, Keller HB, Kernevez JP. 1991a Numerical analysis and control of bifurcation problems (I): bifurcation in finite dimensions. *International Journal of Bifurcation and Chaos* **01**, 493-520. (10.1142/S0218127491000397)

28. Doedel E, Keller HB, Kernevez JP. 1991b Numerical analysis and control of bifurcation problems (II): bifurcation in infinite dimensions. *International Journal of Bifurcation and Chaos* **01**, 745-772. (10.1142/S0218127491000555)

29. Boullé N, Farrell PE, Paganini A. 2022 Control of Bifurcation Structures using Shape Optimization. *SIAM Journal on Scientific Computing* **44**, A57–A76. (10.1137/21M1418708)

30. Moore G, Spence A. 1980 The Calculation of Turning Points of Nonlinear Equations. *SIAM Journal on Numerical Analysis* **17**, 567-576.

31. Boullé N, Farrell PE, Rognes ME. 2023 Optimization of Hopf Bifurcation Points. *SIAM Journal on Scientific Computing* **45**, B390-B411. (10.1137/22M1474448)

32. Pedersen N. 2005 Designing plates for minimum internal resonances. *Structural and Multidisciplinary Optimization* **30**, 297-307. (10.1007/s00158-005-0529-x)

33. Dou S, Jensen JS. 2015 Optimization of nonlinear structural resonance using the incremental harmonic balance method. *Journal of Sound and Vibration* **334**, 239-254. (https://doi.org/10.1016/j.jsv.2014.08.023)

34. Dou S, Strachan BS, Shaw SW, Jensen JS. 2015 Structural optimization for nonlinear dynamic response. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **373**, 20140408. (10.1098/rsta.2014.0408)

35. Dou S, Jensen JS. 2016 Optimization of hardening/softening behavior of plane frame structures using nonlinear normal modes. *Computers & Structures* **164**, 63-74. (https://doi.org/10.1016/j.compstruc.2015.11.001)

36. Li LL, Polunin PM, Dou S, Shoshani O, Scott Strachan B, Jensen JS, Shaw SW, Turner KL. 2017 Tailoring the nonlinear response of MEMS resonators using shape optimization. *Applied Physics Letters* **110**, 081902. (10.1063/1.4976749)

37. Detroux T, Noël JP, Kerschen G. 2021 Tailoring the resonances of nonlinear mechanical systems. *Nonlinear Dynamics* **103**, 3611-3624. (10.1007/s11071-020-06002-w)

38. Grenat C, Baguet S, Lamarque CH, Dufour R. 2019 A multi-parametric recursive continuation method for nonlinear dynamical systems. *Mechanical Systems and Signal Processing* **127**, 276-289. (https://doi.org/10.1016/j.ymssp.2019.03.011)

39. Denimal E, El Haddad F, Wong C, Salles L. 2021 Topological Optimization of Under-Platform Dampers With Moving Morphable Components and Global Optimization Algorithm for Nonlinear Frequency Response. *Journal of Engineering for Gas Turbines and Power* **143**, 021021 (10.1115/1.4049666)

40. Denimal E, Renson L, Wong C, Salles L. 2022 Topology optimisation of friction under-platform dampers using moving morphable components and the efficient global optimization algorithm. *Structural and Multidisciplinary Optimization* **65**, 56. (10.1007/s00158-021-03158-w)

41. Habib G, Kerschen G. 2015 Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* **471**, (10.1098/rspa.2014.0976)

42. Marler R, Arora J. 2004 Survey of multi-objective optimization methods for engineering. *Structural and Multidisciplinary Optimization* **26**, 369-395. (10.1007/s00158-003-0368-6)

43. Steven G. Johnson The NLopt nonlinear-optimization package. http://github.com/stevengj/nlopt.

44. Powell MJD. 1994 pp. 51-67. In *A Direct Search Optimization Method That Models the*
*Objective and Constraint Functions by Linear Interpolation*, pp. 51–67. Springer Netherlands.
(https://doi.org/10.1007/978-94-015-8330-5_4)

45. Powell MJD. 1998 Direct search algorithms for optimization calculations. *Acta Numerica* 7, 287–336. (10.1017/S0962492900002841)

46. Nelder JA, Mead R. 1965 A Simplex Method for Function Minimization. *The Computer Journal* 7, 308–313. (10.1093/comjnl/7.4.308)

47. Box MJ. 1965 A New Method of Constrained Optimization and a Comparison With Other Methods. *The Computer Journal* 8, 42–52. (10.1093/comjnl/8.1.42)

48. Richardson JA, Kuester JL. 1973 Algorithm 454: The Complex Method for Constrained Optimization [E4]. *Commun. ACM* 16, 487–489. (10.1145/355609.362324)

49. Doedel EJ, Champneys AR, Fairgrieve TF, Kuznetsov YA, Sandstede B, Wang X. 1997 AUTO 97: Continuation And Bifurcation Software For Ordinary Differential Equations (with HomCont).

50. Govaerts W, Kuznetsov Y, Dhooge A. 2005 Numerical continuation of bifurcations of limit cycles in MATLAB. *SIAM Journal ON Scientific Computing* 27, 231–252.

51. Dankowicz H, Schilder F. 2010 An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints. *Journal of Computational and Nonlinear Dynamics* 6. (https://doi.org/10.1115/1.4002684)

52. Bernd Krauskopf, Hinke M. Osinga JGV. 2007 *Numerical Continuation Methods for Dynamical Systems*. Springer Dordrecht. (https://doi.org/10.1007/978-1-4020-6356-5)

53. Urabe M. 1965 Galerkin's procedure for nonlinear periodic systems. *Archive for Rational Mechanics and Analysis* 20, 120–152. (10.1007/BF00284614)

54. Karkar S, Cochelin B, Vergez C. 2013 A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities. *Journal of Sound and Vibration* 332, 968–977. (https://doi.org/10.1016/j.jsv.2012.09.033)

55. Guillot L, Cochelin B, Vergez C. 2019 A Taylor series-based continuation method for solutions of dynamical systems. *Nonlinear Dynamics* 98, 2827–2845. (https://doi.org/10.1007/s11071-019-04989-5)

56. Krack M, Gross J. 2019 *Harmonic Balance for Nonlinear Vibration Problems*. (https://doi.org/10.1007/978-3-030-14023-6)

57. Karkar S, Cochelin B, Vergez C. 2014 A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems. *Journal of Sound and Vibration* 333, 2554–2567. (https://doi.org/10.1016/j.jsv.2014.01.019)

58. Gilmore RJ, Steer MB. 1991a Nonlinear circuit analysis using the method of harmonic balance—A review of the art. Part I. Introductory concepts. *International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering* 1, 22–37. (https://doi.org/10.1002/mmce.4570010104)

59. Gilmore RJ, Steer MB. 1991b Nonlinear circuit analysis using the method of harmonic balance—a review of the art. II. Advanced concepts. *International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering* 1, 159–180. (https://doi.org/10.1002/mmce.4570010205)

60. Shahrzad P, Mahzoon M. 2002 Limit cycle flutter of airfoils in steady and unsteady flows. *Journal of Sound and Vibration* 256, 213–225. (https://doi.org/10.1006/jsvi.2001.4113)

61. Ekici K, Kielb RE, Hall KC. 2013 The effect of aerodynamic asymmetries on turbomachinery flutter. *Journal of Fluids and Structures* 36, 1–17. (https://doi.org/10.1016/j.jfluidstructs.2012.08.009)

62. Hall KC, Thomas JP, Clark WS. 2002 Computation of Unsteady Nonlinear Flows in Cascades Using a Harmonic Balance Technique. *AIAA Journal* 40, 879–886. _eprint_:
https://doi.org/10.2514/2.1754 (10.2514/2.1754)

63. Alvarez J, Meraz M, Valdes-Parada FJ, Alvarez-Ramirez J. 2012 First-harmonic balance analysis for fast evaluation of periodic operation of chemical processes. *Chemical Engineering Science* 74, 256–265. (https://doi.org/10.1016/j.ces.2012.02.048)

64. Linkens D. 1979 Modulation analysis of forced non-linear oscillators for biological modelling. *Journal of Theoretical Biology* 77, 235–251. (https://doi.org/10.1016/0022-5193(79)90356-4)

65. Cameron TM, Griffin J. 1989 An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems. *Journal of Applied Mechanics*. (https://doi.org/10.1115/1.3176036)

66. Peletan L, Baguet S, Torkhani M, Jacquet-Richardet G. 2013 A comparison of stability
computational methods for periodic solution of nonlinear problems with application to
rotordynamics. *Nonlinear Dynamics* **72**, 671–682. (https://doi.org/10.1007/s11071-012-0744-0)

67. Groll GV, Ewins D. 2001 The harmonic balance method with arc-length continuation in rotor/stator contact problems. *Journal of Sound and Vibration* **241**, 223 – 233. (https://doi.org/10.1006/jsvi.2000.3298)

68. Mélot A, Benaïcha Y, Rigaud E, Perret-Liaudet J, Thouverez F. 2021 Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train. *Journal of Sound and Vibration* p. 116495. (https://doi.org/10.1016/j.jsv.2021.116495)

69. Bayer F, Leine RI. 2023 Sorting-free Hill-based stability analysis of periodic solutions through Koopman analysis. *Nonlinear Dynamics* **111**, 8439–8466. (10.1007/s11071-023-08247-7)

70. Guillot L, Lazarus A, Thomas O, Vergez C, Cochelin B. 2020 A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems. *Journal of Computational Physics* **416**, 109477. (https://doi.org/10.1016/j.jcp.2020.109477)

71. Tisseur F, Meerbergen K. 2001 The Quadratic Eigenvalue Problem. *SIAM Review* **43**, 235–286.

72. Govaerts W, Sijnave B. 1999 Matrix manifolds and the Jordan structure of the bialternate matrix product. *Linear Algebra and its Applications* **292**, 245 – 266. (https://doi.org/10.1016/S0024-3795(99)00039-7)

73. Volvert M, Kerschen G. 2021 Phase resonance nonlinear modes of mechanical systems. *Journal of Sound and Vibration* **511**, 116355. (https://doi.org/10.1016/j.jsv.2021.116355)

74. Dimitriadis G. 2017 *Introduction to Nonlinear Aeroelasticity*. John Wiley & Sons.

75. Beregi S, Takacs D, Stepan G. 2019 Bifurcation analysis of wheel shimmy with non-smooth effects and time delay in the tyre-ground contact. *Nonlinear Dynamics* **98**, 841–858.

76. Kalmár-Nagy T, Stépán G, Moon FC. 2001 Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations. *Nonlinear Dynamics* **26**, 121–142.

77. Adimy M, Crauste F, Ruan S. 2005 Stability and Hopf bifurcation in a mathematical model of pluripotent stem cell dynamics. *Nonlinear Analysis: Real World Applications* **6**, 651–670.

78. Guo S, Huang L. 2003 Hopf bifurcating periodic orbits in a ring of neurons with delays. *Physica D: Nonlinear Phenomena* **183**, 19–44.

79. McGurk M, Yuan J. 2022 Computation of limit cycle oscillations and their stabilities in nonlinear aeroelastic systems using harmonic balance methods. In *International Forum on Aeroelasticity and Structural Dynamics (IFASD)* 2022.

80. Thouverez F. 2003 Presentation of the ECL benchmark. *Mechanical Systems and Signal Processing* **17**, 195–202. (https://doi.org/10.1006/mssp.2002.1560)

81. Peeters M, Viguié R, Sérandour G, Kerschen G, Golinval JC. 2009 Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques. *Mechanical Systems and Signal Processing* **23**, 195–216. Special Issue: Non-linear Structural Dynamics (https://doi.org/10.1016/j.ymssp.2008.04.003)

82. Shaw A, Hill T, Neild S, Friswell M. 2016 Periodic responses of a structure with 3:1 internal resonance. *Mechanical Systems and Signal Processing* **81**, 19–34. (https://doi.org/10.1016/j.ymssp.2016.03.008)

83. Renson L, Shaw A, Barton D, Neild S. 2019 Application of control-based continuation to a nonlinear structure with harmonically coupled modes. *Mechanical Systems and Signal Processing* **120**, 449–464. (https://doi.org/10.1016/j.ymssp.2018.10.008)

84. Grappasonni C, Habib G, Detroux T, Wang F, Kerschen G, Jensen J. 2014 Practical design of a nonlinear tuned vibration absorber. In *Proceedings of the International conference on Noise and Vibration Engineering (ISMA)*.

85. Grappasonni C, Habib G, Detroux T, Kerschen G. 2016 Experimental Demonstration of a 3D-Printed Nonlinear Tuned Vibration Absorber. In *Nonlinear Dynamics, Volume 1* pp. 173–183. Springer International Publishing.