A Mélot

E Denimal

L Renson

Multi-parametric optimization for controlling bifurcation structures

Keywords: mechanical engineering, applied mathematics, computational mechanics harmonic balance method, bifurcation analysis, stability, shape optimization, nonlinear vibrations

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 2 ). 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.

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 [START_REF] Marco | Bifurcation analysis of single-cell gene expression data reveals epigenetic landscape[END_REF][START_REF] Gerstner | Spiking Neuron Models: Single Neurons, Populations, Plasticity[END_REF][START_REF] Hay | Models of Neocortical Layer 5b Pyramidal Cells Capturing a Wide Range of Dendritic and Perisomatic Active Properties[END_REF]; in chemical engineering, they serve as a powerful tool to identify the optimal operating conditions of chemical reactors [START_REF] Guckenheimer | Multiple bifurcation problems for chemical reactors[END_REF]; in laser physics, they can be exploited to operate bi-stable lasers, which find applications as optical flip-flop memory apparatuses [START_REF] Liu | An ultra-small, low-power, all-optical flip-flop memory on a silicon chip[END_REF], or to produce pulsed beams via passive Q-switching [START_REF] Arimondo | Repetitive passive Qswitching and bistability in lasers with saturable absorbers[END_REF][START_REF] Doedel | Bifurcation structures in a model of a CO 2 laser with a fast saturable absorber[END_REF]. In mechanics, the multi-stability induced by the presence of bifurcations has been exploited for mass detection [START_REF] Kumar | Bifurcationbased mass sensing using piezoelectrically-actuated microcantilevers[END_REF][START_REF] Nguyen | Bifurcation-based micro-/nanoelectromechanical mass detection[END_REF], energy harvesting [START_REF] Arrieta | A piezoelectric bistable plate for nonlinear broadband energy harvesting[END_REF][START_REF] Abdelkefi | Design of piezoaeroelastic energy harvesters[END_REF], shape morphing [START_REF] Kuder | Design space of embeddable variable stiffness bi-stable elements for morphing applications[END_REF], and energy transfer [START_REF] Hurel | Neimark Sacker bifurcations and non-linear energy exchange in chains of non-linear oscillators[END_REF]. The resonances of lightly damped nonlinear systems are also found near saddle-node bifurcations [START_REF] Hill | Interpreting the forced responses of a twodegree-of-freedom nonlinear oscillator using backbone curves[END_REF][START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF][START_REF] Cenedese | How do conservative backbone curves perturb into forced responses? A Melnikov function analysis[END_REF]. 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 [START_REF] Rose | Influence of aeroelastic flow-induced oscillations on fatigue life of an airplane wing structure[END_REF][START_REF] Nilsson | A theoretical and experimental investigation of buckling induced delamination growth[END_REF][START_REF] Renson | Complex dynamics of a nonlinear aerospace structure: numerical continuation and normal modes[END_REF].

In this paper, we propose a general computational methodology to control, through multiparametric 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 [START_REF] Szep | Parameter Inference with Bifurcation Diagrams[END_REF] 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 [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. 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) [START_REF]Elements of Applied Bifurcation Theory[END_REF]. Bifurcation tracking methods have recently gained popularity within the mechanical vibration community as they enable the identification of stability boundaries and safe operating conditions [START_REF] Xie | Bifurcation tracking by Harmonic Balance Method for performance tuning of nonlinear dynamical systems[END_REF][START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF][START_REF] Mélot | Bifurcation tracking of geared systems with parameter-dependent internal excitation[END_REF][START_REF] Mélot | Robust design of vibro-impacting geared systems with uncertain tooth profile modifications via bifurcation tracking[END_REF]. 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 [START_REF] Doedel | a Numerical analysis and control of bifurcation problems (I): bifurcation in finite dimensions[END_REF][START_REF] Doedel | b Numerical analysis and control of bifurcation problems (II): bifurcation in infinite dimensions[END_REF] 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 [START_REF] Mélot | Bifurcation tracking of geared systems with parameter-dependent internal excitation[END_REF]. 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 [START_REF] Boullé | Control of Bifurcation Structures using Shape Optimization[END_REF] 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 [START_REF] Xie | Bifurcation tracking by Harmonic Balance Method for performance tuning of nonlinear dynamical systems[END_REF][START_REF] Moore | The Calculation of Turning Points of Nonlinear Equations[END_REF], to control the location of a single branch point. This approach was extended to Hopf bifurcations [START_REF] Boullé | Optimization of Hopf Bifurcation Points[END_REF] 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 [START_REF] Pedersen | Designing plates for minimum internal resonances[END_REF] 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 [START_REF] Dou | Optimization of nonlinear structural resonance using the incremental harmonic balance method[END_REF] 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 [START_REF] Dou | Structural optimization for nonlinear dynamic response[END_REF], 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 [START_REF] Dou | Optimization of hardening/softening behavior of plane frame structures using nonlinear normal modes[END_REF]. The methodology developed in [START_REF] Dou | Structural optimization for nonlinear dynamic response[END_REF] was applied to the design of micro-electromechanical systems in [START_REF] Li | Tailoring the nonlinear response of MEMS resonators using shape optimization[END_REF]. More recently, Detroux et al [START_REF] Detroux | Tailoring the resonances of nonlinear mechanical systems[END_REF] 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 [START_REF] Grenat | A multi-parametric recursive continuation method for nonlinear dynamical systems[END_REF] 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 [START_REF] Denimal | Topological Optimization of Under-Platform Dampers With Moving Morphable Components and Global Optimization Algorithm for Nonlinear Frequency Response[END_REF][START_REF] Denimal | Topology optimisation of friction underplatform dampers using moving morphable components and the efficient global optimization algorithm[END_REF] 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 [START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF], 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.

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.

l i ≤ θ i ≤ b u i ∀i ∈ 1, p , (2.1)
where fo is the objective function to be minimized, θ θ θ ∈ 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 l i or b u i can be set to infinity if the parameter is only restricted to R + * or R - * .

The proposed objective function is formulated as

fo(T , P) = |T -P|Ψ (θ θ θ) Bifurcation measure + 1 |T | τ ∈T π(θ θ θ)∈P π(θ θ θ) -τ τ 1/|P| Error measure , (2.2)
where T represents the set of target bifurcations and P is the set of bifurcations predicted by the model at each step of the optimisation solver. The notation |T | simply represents the number of elements in the set 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 θ θ θ 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 R(q, µ µ µ, θ θ θ) = 0, (2.3) where q ∈ R n is the solution vector and µ µ µ 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 θ θ θ could be considered as bifurcation parameters, they are here kept separate from µ µ µ 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 µ, 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 minimallyextended systems is used (see Section 3). Although bifurcations are detected through changes in the parameter µ, they also depend on the design variables θ θ θ.

(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 [START_REF] Boullé | Control of Bifurcation Structures using Shape Optimization[END_REF][START_REF] Boullé | Optimization of Hopf Bifurcation Points[END_REF], 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, Ψ (θ θ θ), 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 Ψ (θ θ θ) 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(θ θ θ) = |g(θ θ θ)| max R=0 |g(θ θ θ)| . (2.4)
In order to account for bifurcations occurring over a given bifurcation parameter interval, function G(θ θ θ) 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:

Ψ (θ θ θ) = R=0 G(θ θ θ)ds R=0 ds . (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 [START_REF] Szep | Parameter Inference with Bifurcation Diagrams[END_REF]. 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:

error measure = 1 |T | τ ∈T   π(θ θ θ)∈P πµ(θ θ θ) -τµ τµ 1/|P| + π(θ θ θ)∈P πα(θ θ θ) -τα τα 1/|P|   , (2.6)
where the subscripts µ and α 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 [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]. The optimization

P D ) + α N S fo(T N S , P N S ) subject to b l i ≤ θ i ≤ b u i ∀i ∈ 1, p (2.7
) where fo is given by Eq. (2.2) and the subscripts BP , F , P D and N S 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 [START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF] since it offers numerous subroutines to choose from with support for bound constraints. Two gradient-free algorithms, namely COBYLA (Constrained Optimization BY Linear Approximation) [START_REF] Powell | A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation[END_REF][START_REF] Powell | Direct search algorithms for optimization calculations[END_REF] and Nelder-Mead [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF][START_REF] Box | A New Method of Constrained Optimization and a Comparison With Other Methods[END_REF][START_REF] Richardson | Algorithm 454: The Complex Method for Constrained Optimization [E4[END_REF] 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.

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 [START_REF] Doedel | AUTO 97: Continuation And Bifurcation Software For Ordinary Differential Equations[END_REF][START_REF] Govaerts | Numerical continuation of bifurcations of limit cycles in MATLAB[END_REF][START_REF] Dankowicz | An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints[END_REF]. For instance, time-domain methods such as the orthogonal collocation [START_REF] Doedel | AUTO 97: Continuation And Bifurcation Software For Ordinary Differential Equations[END_REF][START_REF] Govaerts | Numerical continuation of bifurcations of limit cycles in MATLAB[END_REF][START_REF] Dankowicz | An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints[END_REF][START_REF] Krauskopf | Numerical Continuation Methods for Dynamical Systems[END_REF] or shooting methods [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] are popular methods to determine periodic responses. A powerful alternative considered here is the HBM [START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF][START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF][START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF][START_REF] Guillot | A Taylor series-based continuation method for solutions of dynamical systems[END_REF][START_REF] Krack | Harmonic Balance for Nonlinear Vibration Problems[END_REF]. Formulated in the frequency domain, the HBM is known to outperform time domain techniques when solving large-scale systems with strong nonlinearities [START_REF] Karkar | A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems[END_REF] and has been employed in numerous fields of engineering and the applied sciences, including electrical engineering [START_REF] Gilmore | a Nonlinear circuit analysis using the method of harmonic balance-A review of the art. Part I. Introductory concepts[END_REF][START_REF] Gilmore | Nonlinear circuit analysis using the method of harmonic balance-a review of the art. II. Advanced concepts[END_REF], aeroelasticity [START_REF] Shahrzad | Limit cycle flutter of airfoils in steady and unsteady flows[END_REF][START_REF] Ekici | The effect of aerodynamic asymmetries on turbomachinery flutter[END_REF], fluid dynamics [START_REF] Hall | Computation of Unsteady Nonlinear Flows in Cascades Using a Harmonic Balance Technique[END_REF], chemical engineering [START_REF] Alvarez | First-harmonic balance analysis for fast evaluation of periodic operation of chemical processes[END_REF] and biology [START_REF] Linkens | Modulation analysis of forced non-linear oscillators for biological modelling[END_REF].

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

Mq(t) + C q(t) + Kq(t) + f nl (q, q, t) = fex(t), (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 fex, 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 = R ∞ k=0 qk e ikΩt ≈ R H k=0 qk e ikΩt , (3.2)
where q contains the coefficients of the one-sided Fourier transform. In this paper, the response frequency, Ω, coincides with the excitation frequency and corresponds to the bifurcation parameter. The displacements, velocities and accelerations may be expressed as

q = R ([T ⊗ In] q) , (3.3)
q = R ([TΩ∇ ⊗ In] q) , (3.4)
q = R TΩ 2 ∇ 2 ⊗ In q , (3.5)
where ⊗ is the Kronecker product, In the identity matrix of size n × n, T the vector of orthogonal exponential functions and ∇ the complex frequency domain differential operator

∇ = diag(0, i, 2i, . . . , Hi). (3.6)
Similarly, the nonlinear forces and external forcing are written as

f nl = R H k=0 fnl,k e ikΩt , (3.7
)

fex = R H k=0 fex,k e ikΩt . (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 T. This yields a residual R consisting of a set of n(H + 1) nonlinear algebraic equations: R(q, Ω) = Z(Ω)q + fnl (q) -fex = 0, (

where fnl and fex contain the Fourier coefficients of the nonlinear forces and external forcing, respectively, and Z(Ω) is the frequency domain dynamic stiffness matrix:

Z(Ω) = Ω 2 ∇ 2 ⊗ M + Ω∇ ⊗ C + I H+1 ⊗ K. (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 [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] 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 Rex 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 [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] 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 [START_REF] Peletan | computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF] or Hill's method [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF][START_REF] Mélot | Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train[END_REF][START_REF] Bayer | Sorting-free Hill-based stability analysis of periodic solutions through Koopman analysis[END_REF]. Hill's method is employed here as it is particularly well suited to solutions computed with the HBM [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF][START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF]. The main underlying idea is to introduce a perturbation of the form: q 0 (t) = q(t) + p(t)e Λt .

(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:

Λ 2 M + Λ C + ∂ qR p = 0, (3.12)
with

M = I 2H+1 ⊗ M, (3.13)
C = ∇ ⊗ 2M + I 2H+1 ⊗ C, (3.14)
where ∂ qR 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 ∇

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 [START_REF] Tisseur | The Quadratic Eigenvalue Problem[END_REF]:

C ∂ qR -I n(2H+1) 0 + λ M 0 0 I n(2H+1) = 0. (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 scenarii. 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:

A b d † 0 w g = 0 1 , (3.16)
were b and d are chosen so that the system is nonsingular and † represents the hermitian transpose. The matrix A takes different values depending on the type of bifurcation to detect. For fold bifurcations, A = ∂ qR is used. For branch points: diagonal matrix gathering the eigenvalues of Eq. (3.15) and ⊙ is the bialternate matrix product [START_REF] Govaerts | Matrix manifolds and the Jordan structure of the bialternate matrix product[END_REF]. Both bifurcation types can then be distinguished based on the imaginary part of the eigenvalues whose real part becomes positive.

A = ∂ qR ∂ Ω R z T , (3.17

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

mq + c q + kq + k nl q 3 = fex (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 fex = 3 cos(Ωt).

Five harmonics and the mean value are considered in the analysis and the displacements and nonlinear forces are sampled with Ns = 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 [START_REF] Volvert | Phase resonance nonlinear modes of mechanical systems[END_REF].

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 θ θ θ = (m, c, k, k nl ) T and only fex 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 Starting from the initial configuration, a first optimization run is carried out to enforce the appearance of two fold bifurcations at angular frequencies Ω = 1.59 rad/s and Ω = 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 fo = 2.0 • 10 -4 . The computation of the dimensionless parameters ε = k nl f 2 ex /k 3 and δ = c/2 √ k reveals that the system with optimized parameters is significantly less nonlinear (ε = 0.9) than the initial system (ε = 4.5) and a little bit more damped (δ = 6.3 • 10 -2 compared to δ = 5.0 • 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 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 q1 = 1.5 and q1 = 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 (Ω = 3.50 rad/s, Ω = 2.00 rad/s, Ω = 0.56 rad/s, Ω = 0.54 rad/s) and two branch points (Ω = 0.87 rad/s, Ω = 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 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 (Ω, 1 q1 ) = ( 2 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 [START_REF] Dimitriadis | Introduction to Nonlinear Aeroelasticity[END_REF], vehicle wheels [START_REF] Beregi | Bifurcation analysis of wheel shimmy with non-smooth effects and time delay in the tyre-ground contact[END_REF], machine tools [START_REF] Kalmár-Nagy | Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations[END_REF], and living cells [START_REF] Adimy | Stability and Hopf bifurcation in a mathematical model of pluripotent stem cell dynamics[END_REF][START_REF] Guo | Hopf bifurcating periodic orbits in a ring of neurons with delays[END_REF]. 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 [START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF], i.e. a van der Pol-Duffing oscillator coupled to a NLTVA with cubic nonlinearity. The equations of motion are

m 1 q1 + c 1 q 2 1 -1 q1 + kq 1 + k nl1 q 3 1 + c 2 ( q1 -q2 ) + k 2 (q 1 -q 2 ) + k nl2 (q 1 -q 2 ) 3 = 0, m 2 q2 + c 2 ( q2 -q1 ) + k 2 (q 2 -q 1 ) + k nl2 (q 2 -q 1 ) 3 = 0, (4.2)
where subscripts 1 and 2 denote quantities related to the vdPD oscillator and the NLTVA, respectively. For practical reasons, the mass ratio ϵ = m 2 /m 1 = 0.05 is kept fixed [START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF]. 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 (θ θ θ = (k 2 , k nl2 , c 2 ) T ) and initially set at k 0 2 = 4.625 • 10 -2 , k 0 nl2 = 6.8 • 10 -4 and c 0 2 = 1.15 • 10 -2 . The continuation parameter is c 1 as in [START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF]. 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 q1 (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 [START_REF] Mcgurk | Computation of limit cycle oscillations and their stabilities in nonlinear aeroelastic systems using harmonic balance methods[END_REF]. 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 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:

error measure = c Hopf 1,π -c Hopf 1,τ c Hopf 1,τ , (4.3)
where c Both a lower bound b l = (0.5k 0 2 , 0.1k 0 nl2 , 0.8c 0 2 ) T and upper bound b u = (1.5k 0 2 , 10k 0 nl2 , 1.2c 0 2 ) T are defined for the structural parameters. The initial guess for the optimization corresponds to the mechanical properties θ θ θ 0 = (k 0 2 , k 0 nl2 , c 0 2 ) 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 • 10 -2 , k nl2 = 6.8 • 10 -5 and c 2 = 1.27 • 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 [START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF][START_REF] Thouverez | Presentation of the ECL benchmark[END_REF][START_REF] Peeters | Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques[END_REF] 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%. (ii) Shape optimization

The structural parameters considered for the optimization are θ θ θ = (L, w T e , h T e , k nl ) T where we and he denote the vectors containing the width and height, w (i) e and h (i) e , i ∈ 1, Ne , of all the elements. The total number of optimization variables is 2Ne + 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 • 10 8 and 6.0 • 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 ω 1 of the underlying linear system, i.e. [ω 1 , 4ω 1 ]. Note that the very same procedure could be applied to other modes. Besides, the continuation procedure is initiated at Ω = 4ω 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 [START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF]. Such a scenario, illustrated in Figure 9a, can be extremely detrimental to the structural integrity of the system as the 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 Ω = 47.7 Hz and Ω = 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 Ω = 75 Hz and Ω = 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 Ω = 48.3 Hz and Ω = 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 (≈ 2 Hz) on the location of the lower fold bifurcation (near Ω = 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 [START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF][START_REF] Shaw | Periodic responses of a structure with 3:1 internal resonance[END_REF][START_REF] Renson | Application of control-based continuation to a nonlinear structure with harmonically coupled modes[END_REF]. It is therefore quite difficult to obtain significant changes in their locations when considering only the first resonance and constraining its corresponding fold bifurcations. 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 [START_REF] Xie | Bifurcation tracking by Harmonic Balance Method for performance tuning of nonlinear dynamical systems[END_REF] 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 [START_REF] Grappasonni | Practical design of a nonlinear tuned vibration absorber[END_REF][START_REF] Grappasonni | Experimental Demonstration of a 3D-Printed Nonlinear Tuned Vibration Absorber[END_REF].

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.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The optimization problem takes the general form minimize θ θ θ fo subject to b

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Alternating frequency/time procedure using the fast Fourier transform and its inverse.

Figure 2 :

 2 Figure 2: Bifurcation scenarii depending on the trajectory of the Floquet exponents. Fold bifurcation and branch point (a), period doubling bifurcation (b) and Neimark-Sacker bifurcation (c).

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . are allowed to take values within the interval [b l i , ∞[ where b l i , i ∈ 1, 4 are the components of vector b l = (0.01, 0.01, 0.01, 0.01) T .

Figure 3 :

 3 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.

Figure 4 :

 4 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.

Figure 5 :Figure 6 :

 56 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.

  , 1.5) and (Ω, q1 ) = (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 (ε = 0.65) and more damped (δ = 7.7 • 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 royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 7 :

 7 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.

Hopf 1 ,

 1 τ = 0.183 is the targeted Hopf location and c Hopf 1,π 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.

  . The initial model has a length L = 0.7 m, a height and a width of 0.014 m. The beam is discretized with royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ne = 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(Ωt) is applied at 0.21 m from the clamped end.

FFigure 8 :

 8 Figure 8: Schematic representation of the Euler-Bernoulli beam.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 9 :

 9 Figure 9: Forced response curve of the beam with initial (a) and optimized (b) geometries for F = 90 sin(Ω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(Ωt).

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 10 :Figure 11 :

 1011 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.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	𝐪	IFFT	𝐪 t , ሶ 𝐪(t)
	ሚ 𝐟 nl	FFT	𝐟 nl (𝐪, ሶ 𝐪, t)
	Frequency domain		Time domain

Table 1 :

 1 Reference geometrical and mechanical properties of the considered beam[START_REF] Thouverez | Presentation of the ECL benchmark[END_REF].

	Length (m)	Width (m)	Height (m)
	0.7	0.014	0.014
	Young's modulus (N/m 2 ) Density (kg/m 3 ) Nonlinear coeff. (N/m 3 )
	2.05 • 10 11	7800	6.0 • 10 9

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