Adjoint-based limit cycle oscillation instability sensitivity and suppression

Dynamical systems often exhibit limit cycle oscillations (LCOs), self-sustaining oscillations of limited amplitude. LCOs can be supercritical or subcritical. The supercritical response is benign, while the subcritical response can be bi-stable and exhibit a hysteretic response. Subcritical responses can be avoided in design optimization by enforcing LCO stability. However, many high-fidelity system models are computationally expensive to evaluate. Thus, there is a need for an efficient computational approach that can model instability and handle hundreds or thousands of design variables. To address this need, we propose a simple metric to determine the LCO stability using a fitted bifurcation diagram slope. We develop an adjoint-based formula to efficiently compute the stability derivative with respect to many design variables. To evaluate the stability derivative, we only need to compute the time-spectral adjoint equation three times, regardless of the number of design variables. The proposed adjoint method is verified with finite differences, achieving a five-digit agreement between the two approaches. We consider a stability-constrained LCO parameter optimization problem using an analytic model to demonstrate that the optimizer can suppress the instability. We also consider a more realistic LCO speed and stability-constrained airfoil problem that minimizes the normalized mass and stiffness. The proposed method could be extended to optimization problems with a partial differential equation (PDE)-based model, opening the door to other applications where high-fidelity models are needed.

In many cases, such as flutter, delaying LCO onset is desired [7]. This can be handled by design optimization (see the review by [4] and recent work by the authors [13]). For partial differential equation (PDE) constrained optimization problems with large numbers of design variables, the most promising approach is to use a gradient-based optimizer with gradients computed by the adjoint method [14]. This is because the adjoint method can compute the derivative at a cost that is independent of the number of design variables [15,Sec. 6.7]. There has been significant progress in applying the adjoint method to dynamical systems: fixed point [16][17][18][19], LCO [2,20], and chaotic dynamical system [21][22][23].
For derivatives related to LCO, most of the previous research focused on computing derivatives of LCO parameters, such as amplitude, phase, or period [2,20,24]. The stability of LCO and its derivative are also critical characteristics. The stability of LCO is determined by the behavior of its neighboring trajectories. If all the neighboring trajectories spiral toward the LCO, the LCO is stable. On the other hand, if all the neighboring trajectories spiral away from the LCO, it is unstable. Finally, if some trajectories spiral toward and others away from the LCO, it is semi-stable. In a supercritical bifurcation diagram, all of its neighboring trajectories are stable, while in a subcritical bifurcation diagram, a portion of its neighboring trajectories are unstable (see Fig. 1).
LCO stability is traditionally evaluated using Floquet theory [25][26][27], continuation methods [28,29], or central manifold-based methods [30][31][32][33][34]. Floquet theory states that the stability of a periodic dynamical system is determined by the eigenvalue distribution of the monodromy matrix. The monodromy matrix is obtained by evaluating the fundamental solution matrix of the periodic system at t = 0 and t = T , where t is the time, and T is one period of the system [25,Ch. 2.4]. In practice, the monodromy matrix is obtained by simulating the time-periodic dynamical system N times for one period, where N is the number of states of the dynamical system. This is prohibitively expensive for problems governed by PDEs, particularly when N is large, such as the high-fidelity aeroelastic systems considered by Thomas et al. [7] and He et al. [2].
Bauchau and Nikishkov [35] realized that it is possible to use the Arnoldi algorithm to compute the dominant eigenvalue instead of computing all eigenvalues to determine the stability of a periodic dynamical system. This approach is appealing because the matrix is never formed explicitly; only matrix-vector products are required by the Arnoldi algorithm, where the response vector can be evaluated using a time integration for one period for a given initial condition. However, an unsteady adjoint implementation is required to compute the derivative of this approach, which incurs significant development effort and computational time [36]. The Floquet analysis was differentiated in the forward mode by [37], but the computational cost of such a derivative computation method scales with the number of design variables.
The continuation method tracks the LCO response starting from the bifurcation solution [28,29]. Shukla and Patil [29] developed a method that continuously tracks the LCO commencement speed. They optimized the control variables to increase the LCO commencement speed to get a more benign response. The derivative was computed via a direct implicit analytic method [15,Sec. 6.7]. The computational cost of the direct method scales linearly with the number of design variables. Consequently, this approach can be computationally intractable if detailed high-fidelity models are considered in the design optimization.
The central manifold-based method states that the stability of a large-dimensional dynamical system solution near the bifurcation point can be captured by a lowdimensional manifold named the central manifold [30][31][32][33][34]. To construct the central manifold, higher-order (second-and third-order) derivatives are required ( [32], Ch. 5), making it challenging to implement in practice if the analytic form of the derivatives is not known a priori. Stanford and Beran [34] applied a multi-scale method to capture the bifurcation diagram near the bifurcation point and conducted optimization to suppress the subcritical response. However, as the magnitude increases, the true bifurcation deviates from the approximation. The finite differences used to compute the derivatives make this approach very expensive for problems governed by large-scale PDEs.
Besides the methods mentioned above, several other methods were proposed to simulate LCO. Shukla and Patil [38] developed a nonlinear normal mode method to approximate the LCO magnitude. The approximated LCO magnitude may cause numerical issues when conducting gradient-based optimization, impacting the robustness of such methods in design optimization. Riso et al. [39] proposed an approximate recovery rate as a metric for post-flutter analysis and gradient-based optimization that avoids the expensive transient simulation. The method is demonstrated in a simple optimization problem on a typical airfoil section, minimizing the mass ratio, subject to flutter and the proposed postflutter constraint. The proposed method successfully eliminates any subcritical characteristics and accounts for mode switches and hump modes. However, the method requires sampling multiple flight conditions and modal shapes to evaluate the post-flutter constraint, making it numerically expensive. Stanford et al. [40] proposed a cyclic-finite difference method to solve the periodic response of a nonlinear beam that can potentially capture LCOs.
To mitigate the computational cost and implementation challenges encountered when applying the Floquet theory, we propose a bifurcation diagram-based fitting method to capture the instability, including an adjoint method to evaluate the derivative of the instability metric efficiently with respect to the design variables. A third-order polynomial is fitted by sampling three postbifurcation points and enforcing the intersection angle with the parameter axis at the Hopf-point to be 90 • . In principle, the sampling points can be computed using various approaches, but in this work, they are evaluated by solving the so-called time-spectral LCO equation [7,41,42]. The slope is used as a metric to determine LCO stability. The derivative of the stability metric is then computed using the adjoint method. The fitted curve adjoint requires the solution of the time spectral LCO adjoint [2,20] for all three points. This is the first time the time spectral method has been applied for LCO optimization with the LCO stability constraint. The method can be used to solve optimization problems with many design variables because the derivative computational cost is independent of the number of design variables. However, the proposed method focuses on problems with bifurcation diagrams that can be approximated using a third-order polynomial. Problems with multiple branches and problems that switches to different bifurcation mechanisms in optimization are beyond the scope of the paper.
We consider two optimization problems to demonstrate the proposed method: (1) A simple analytic model and (2) a more practical airfoil section in quasisteady incompressible flow. For the analytic model, we maximize the LCO parameter while also ensuring the LCO stability. In the airfoil case, we minimize a function dependent on the airfoil's weight and softening effect. However, by doing so, the structure may potentially become more flexible and prone to LCO at a lower airspeed. Thus, to avoid this, we constrain the LCO parameter, in this case, airspeed, in addition to the LCO stability constraint.
To summarize, the new contributions in this work are (1) the bifurcation diagram-based fitting method to capture the instability; (2) an efficient adjoint approach to evaluate the derivative of the instability metric with respect to the design variables.
The paper is organized as follows. In Sect. 2, we gave an overview of the problem and optimization setup. In Sect. 3, we discuss the governing time spectral equations used for LCO computation. The proposed fitted curve-based stability criterion is presented in Sect. 4. Section 5 derives the adjoint equations used to compute the underlying LCO parameter derivative and the formula to compute the stability derivative. Then, in Sect. 6, we provide numerical results, including derivative verification and optimization results. Finally, we present our conclusions in Sect. 7.

Problem setup
This section provides an overview of the problem setup and optimization problem statement. A dynamic system can be described by the following general equation: where f : R N → R N is some nonlinear function, μ ∈ R is the LCO parameter, and w(x) ∈ R N is the vector of state variables. Also, when the motion magnitude approaches zero, μ is equal to the LCO parameter at bifurcation (see Fig. 1). The state variables depend on the design variable vector, x ∈ R n x , where n x is the number of design variables, and N is the dimension of the state space.
The optimization problem is formulated as follows: by varying x, subject to g stab (x) >ĝ, where g stab (x) is the stability constraint,ĝ is a nonnegative real number that can be tuned to make sure the solution of the dynamical system is stable, and x and x are the lower and upper bound vectors for the design variable, x, respectively. Here g stab (x) is a scalar since only one operation condition is considered. When performing a multipoint optimization, we can use a vector to store the stability of all points that need to be constrained.
It is also possible to choose a different objective function to be optimized and set both the LCO stability and the LCO parameter as constraints. In this case, the problem is formulated as follows, by varying x, where f (x) is the objective function and μ is a lower bound for the LCO parameter.

Time spectral LCO equation
In this paper, we solve the time spectral LCO equation first proposed by [7] to construct the bifurcation diagram. Alternatively, the bifurcation diagram can be constructed using other methods, e.g., a continuation method, as discussed in Sect. 1. Here, it comprises a time spectral dynamic equation, motion magnitude, and phase equation. The motion magnitude and phase equations are used to form constraints, which are necessary to ensure a unique solution.
The fundamental idea of the time spectral method is to approximate a periodic dynamical system using several equally-spaced time snapshots. The temporal residual is computed by first computing the Fourier coefficient of the periodic state variables, evaluating the time derivative in the frequency domain, and converting it back to the time domain snapshots. The spatial residual is computed in the time domain.
The dynamic system, Eq. (1), can be cast into the time-spectral form as follows. We write the LCO equation in residual form as . . .
where w n ∈ R n×N represents the state variables for all time instances, w i are the states for the i th time instance, and ω is the angular velocity defined as ω = 2π/T . The matrix P is a permutation matrix defined as The matrix D N is composed of N spectral differentiation matrices on the diagonal and is given by where i, j = 1, . . . , n. Here, we have N matrix blocks of D on the diagonal. In addition to the LCO equations, we constrain the magnitude and phase of the motion of one DOF. In residual form, this can be written as where |ŵ 1 | is the motion magnitude of the first DOF from w n , φ is the phase of the motion, andα,φ are the prescribed motion magnitude and phase, respectively.
Here,ŵ 1 is the state variables for all the first DOF. The two additional constraints are used to obtain unique LCO solutions. Without the magnitude constraint, any point from the bifurcation diagram satisfies the dynamical equation (see Fig 1). Similarly, any phase between [0, 2π] will be feasible without the phase constraint. The introduction of the two constraints eliminates the ambiguity, resulting in unique LCO solutions. He et al. [42] provide more details on how the motion magnitude and phase constraints are implemented. Combining Eqs. (4) and (7), the full system of equations becomes where q is the state vector defined as This system can be solved by segregated Newton-Raphson method [7], or by a coupled Newton-Krylov method [41,42]. In this work, we use the latter. Newtontype methods are appealing because of their quadratic convergence rate property. However, all Newton-type methods require the initial point to be within the basinof-attraction to converge, which often leads to issues with solution robustness. In the following sections, we discuss our proposed solution strategy.

Jacobian-free Newton-Krylov method
In this work, we apply a coupled Jacobian-free Newton-Krylov method to solve Eq. (8). The Newton increment q is computed by solving the following equation: where J is the Jacobian matrix of Eq. (8), and q (k) is the solution from the k th iteration. The updated state approximation is defined as where θ is determined by a line search.
When applying a Krylov subspace method to solve Eq. (10), only matrix-vector products are required. The product can be approximated by the directional finite difference thus never explicitly forming and storing the Jacobian. He et al. [42] provides more details on the solver and preconditioners used to obtain the Newton increment in Eq. (10). In this work, we use the Scipy [43] implementation of the Jacobian-free Newton-Krylov method [44].

Warm-start strategy based on Hopf bifurcation analysis
Newton-type methods require the initial point to be reasonably close to the final solution to converge robustly. However, determining or finding a suitable initial point can be challenging and problem-dependent. Globalization methods that address the issue of robustness exist, such as the approximate Newton-Krylov (ANK) [45], but this approach is not implemented here. Initialization with a prescribed motion is usually sufficient for a simple problem. However, this initialization strategy is insufficient for a more complex problem and requires a better strategy. Here we implement a warm-start strategy, proposed by [7], where the initial point is based on a Hopf bifurcation analysis. The idea is that we first linearize Eq. (1) around a steady solution for certain μ. We decompose the state variable into a steady component and a perturbed state where w is a steady solution and δw is a perturbed state. Inserting the decomposed state into Eq. (1), we obtain Linearizing about the steady component of the state, we obtain It follows that the steady solution, w, satisfies f(μ, w) = 0, (16) and the perturbed state variable δw satisfies The Jacobian matrix is defined as We can now identify the bifurcation point from the eigenvalues of the Jacobian matrix. The eigenmode of the solution δw takes the form where λ is the eigenvalue and δw is the eigenvector. Inserting Eq. (19) into Eq. (17) and ignoring the higherorder terms, we obtain The bifurcation happens when the maximum real part (the damping) of one conjugate pair of the eigenvalues is equal to zero, that is, We denote the steady solution of Eq. (16) at the bifurcation point as Similarly, we denote the eigenpair for the eigenvalue problem, Eq. (20), at the bifurcation point as where Re (λ bif ) = 0. There are multiple ways to solve for μ bif . One way is to define a search range and use bisection to find μ bif .
After the bifurcation problem is solved, we can use the solution to initialize the time spectral LCO solver. The process of converting an eigenvalue problem solution to an initial point for the time spectral LCO solver is detailed in Algorithm 1. In lines 2 and 3, we construct the initial guess of the LCO solution LCO parameter and the frequency using the values computed at the bifurcation point. Then, in lines 4-13, we use the bifurcation eigenvector to approximate the LCO motion by scaling the bifurcation eigenvector according to the prescribed motion magnitude and phase of the LCO. To be specific, in line 8, the phase difference between the jth DOF and the first DOF is measured and recorded. In line 9, the phase is computed taking into account the phase difference computed in the previous step. In line 10, the magnitude ratio is computed to compute the motion for the jth DOF using this ratio and the motion magnitude of the first DOF. Finally, the displacement is computed using the computed phase and magnitude information in line 11. The output q (0) is the initial solution found by this warm-start strategy.

Stability criterion via fitted bifurcation diagram
The stability of an LCO point can be determined by the slope of its bifurcation diagram. As a reminder,α is the prescribed motion magnitude of the LCO, and μ is the corresponding LCO parameter. For a bifurcation diagram with μ < μ bif the steady solution is stable, we have indicating that the LCO is stable; if otherwise, dα/dμ < 0, the solution of the system is unstable; and dα/dμ = 0 is a critical point. We assume that when μ < μ bif the steady solution is stable because, otherwise, by defining μ as −μ, the conclusion will be the opposite. This assumption can usually be met in many physical problems. For example, for aeroelastic LCO, where μ can be chosen as the airspeed, we know that the steady solution is stable when the airspeed is low. However, if a chosen parameter does not meet this condition, we can redefine a new LCO parameter as the additive inverse of the previous one, thus satisfying the condition. This can also be verified by checking the Jacobian of a steady solution with an LCO parameter slightly lower than the bifurcation point to ensure that all the eigenvalues are negative.
In this work, the bifurcation diagram is approximated using a third-order polynomial, where a, b, c, and d are coefficients to be determined by three sampling points and an intersection angle constraint. We chose the third-order polynomial because of its simplicity. Other interpolations, such as the non-uniform rational B-splines (NURBS), can also be applied. Three sampling points are considered, where (μ 0 ,α 0 ) denotes the parameters with the original LCO point, and (μ − ,α − ) and (μ + ,α + ) denotes parameters from two neighboring points. Theα 0 is chosen to simulate the original LCO, and the adjacent α − ,α + can be chosen close to the original magnitude to provide a better local estimation, or to be more uniformly distributed to provide a better global estimation. In practice, we found that settingα − ≈ 0.8α 0 and α + ≈ 1.2α 0 gave reasonable approximation. However, these values may differ for different applications. The intersection angle of the bifurcation diagram approximation curve and the x−axis (at the Hopfpoint) should be 90 • , as shown in Fig. 2. In terms of the third-order polynomial fitting, the slope dα/ dμ at α = 0 shall be zero. The slope constraint at the intersection can be satisfied by The remaining coefficients are found by solving the following equation The stability criterion from Eq. (24) can then be written as where g stab is the stability measure. As mentioned before,ĝ is a non-negative parameter that can be used to stabilize the system further. For problems with multiple unstable branches, we can enforce the constraint at several locations on the bifurcation diagram with different motion magnitudes. However, this is beyond the scope of the paper and remains to be investigated in the future. Furthermore, we assume that the bifurcation diagram changes smoothly with respect to the design variables. Problems with a change in the bifurcation mechanism are also out of the scope of this paper.

Derivative computation
For PDE-constrained gradient-based optimization, the derivatives must be computed efficiently, particularly when the function evaluations are expensive. The adjoint method is ideal for such a task because the computational cost is independent of the number of design variables [15]. In this section, we derive an adjointbased approach to compute the derivatives of the proposed constraint.

Computing dμ/ dx using adjoint method
To derive the total derivative of a parameter with respect to the design variables, we apply the adjoint method to Eq. (8). For more details about the adjoint method, we refer the reader to Sect. 6 of the textbook by [15]. The objective, f (μ, ω, w n ; x), is a functional of q (i.e., the states μ, ω, w n ), and the design variables x. For a given x, which is determined by the optimizer, q is found by solving the governing equation. Differentiating the function of interest and Eq. (8), applying the chain rule, we have where the second equality arises because the residual remains zero for any design variable x. Applying the adjoint method, we have where ψ ψ ψ is the adjoint variable. It is evident that by solving one linear equation, we obtain the derivative with respect to all the design variables.

Computing dg stab / dx
The derivative of the constraint dg stab /dx is computed using Eqs. (28), (29), and (31). We use the adjoint equation Eq. (31) by setting the entries of the μ μ μ as the function of interest f . It can be computed using the chain rule where The derivatives dg stab dc , dc dμ μ μ , (34) can be evaluated analytically. The most computationally demanding derivative is It can be computed by applying the above adjoint method three times, independent of the number of design variables. It is also possible to use the adjoint method to compute dα/dμ. However, this makes the computation of its derivative challenging because it requires a second-order adjoint method to compute design variable derivatives, which is expensive. The proposed fitting method circumvents this by formulating the constraint in terms of μ instead of its derivative.

Numerical results
In this section, we first verify the analysis and derivative computation capability of the proposed method in Sect. 6.1. Then, in Sect. 6.2, we consider a simple analytic dynamical system where the initial instability is suppressed, and a higher LCO parameter is sought. Finally, in Sect. 6.3, we demonstrate the proposed approach using an aeroelastic problem where we optimize the dimensionless weight and softening effect. The data of the problems are available online. 1 6.1 Verification

Stability analysis
In this section, we verify that the proposed fitting method can capture and represent the bifurcation diagram characteristics, and the proposed Eq. (32) can compute its slope accurately. Consider the following analytical dynamical systeṁ 1 https://github.com/SichengHe/LCO_stability_fit_data.  Fitted curve, g stab Actual curve, FD −0.74309383 −0.74319293 The dynamical system is constructed in such a way that the bifurcation point is at μ = 1, and the nonlinear term will make the system unstable. Here, we have N = 2, and the number of time instances used to solve the problem is set to n = 5. The center sampling points is arbitrarily chosen with the adjacent points selected following guidelines discussed in Sec. 4. For verification, we have (α − ,α 0 ,α + ) = (0.4, 0.5, 0.6). As a reminder, theα refers to the magnitude of the first DOF, |ŵ 1 |. The fitted and the actual curve for the bifurcation diagram are shown in Fig. 3. Overall, there is a good agreement between the two curves, particularly close to the sampling points, but further away at the extremê α = 1, there is a minimal visible difference. The stability measure, g stab , obtained from the fitted curve using analytic formula, Eq. (29), and by central finite differences at the center point, are compared in Table 1. We achieve a 4-digit match for this specific case, demonstrating that the proposed method can accurately capture the underlying physics.

Derivative verification
In the previous section, we verified the stability measure for the proposed method. In this section, we verify adjoint computation of the derivative of this stability measure, dg stab / dx. The dynamical system for this task is constructed for optimization purposes and now includes design variables and is defined as, The dynamical system is constructed in this way such that the bifurcation point is at μ = (x 1 + x 2 )/2. The nonlinear term can be either stabilizing or destabilizing depending on the design variables value determined by the optimizer. The prescribed motion magnitude is set toα 0 = 0.5, and the delta magnitude is the same as before (±0.1) for the adjacent sampling points. The numerical values are listed in Table 2. We obtain a 4-5 digit match for this case, demonstrating that the adjoint method can be used to compute the derivatives accurately.

Instability suppression using gradient-based optimization
Having verified the derivatives, we can now perform an optimization. The objective is to maximize the LCO parameter, subject to the proposed stability constraint that characterizes the LCO stability for a given design. The optimization problem statement is summarized by Eq. (2). The physical interpretation of this problem is that we want to postpone the onset of bifurcation while ensuring that the LCO branch is stable.
The dynamical system is defined in Eq. (37). Also, we need to specify a motion magnitude for the time spectral LCO solver. As before, we set the prescribed motion magnitude toα 0 = 0.5 and the number of time instances to n = 5. The design variable bounds we consider here are set to be 0 ≤ x ≤ 1. The lower bound of the constraint value from Eq. (29) is set tô g = 0.1, to make sure the solution of the dynamical system is indeed stable. The optimization is conducted using SNOPT [46] with the Python interface provided by pyOptSparse [47]. For more details about the problem definition, see Eq. (2).
It is possible to visualize this problem's design space because there are only two design variables. The contour plots for the constraint and objective functions are shown in Fig. 4. The major iterations, initial, and optimal points are shown in Fig. 4. The optimizer started from the baseline design (x 1 = 0.3, x 2 = 1) and found the optimal point to be at x * 1 = 1, x * 2 = 0.43. The baseline and optimized bifurcation diagrams are shown in Fig. 5. The two extreme bifurcation maps from left to right correspond to the baseline and optimized designs. The intermediate maps correspond to the analysis of the design variables that are equally spaced on a straight line between the baseline and optimized designs, demonstrating that bifurcation diagram characteristics change smoothly. The baseline design has a lower value of μ and turns out to be unstable; the optimized result has a higher value of μ and is stable. This demonstrates the capability of the algorithm for LCO stability optimization.

Problem setup
In this section, we consider an aeroelastic LCO optimization problem. The aeroelastic model used here is adopted from [39], consisting of a flat-plate typical section in a quasi-steady incompressible potential flow. The airfoil is restrained in translation and rotation by springs, and the elastic reaction by the rotational spring includes nonlinear terms. The dynamical system is given as, where the total function, f(μ, w, x), is composed of a linear coefficient matrix, A(μ, x), and a nonlinear vector, F nl (w, x). The state variable vector is defined as where h and α are dimensionless plunging and pitching variables, respectively. The design variable vector is defined as where m is a dimensionless mass per unit length defined as the ratio of the typical section mass per unit length and the mass of the fluid in a circle with the half chord as the radius, and κ (3) α is a dimensionless stiffness parameter defined as the ratio of the third-order and the firstorder spring rotational stiffness constants. In the optimization, we vary m by changing the mass per unit length where the mass of the fluid in a circle with the half chord as the radius is fixed.
The linear coefficient matrix is defined as The structural mass and stiffness matrices are defined as where x α , , and r α are constant parameters. The aerodynamic damping and stiffness matrices are defined as where μ is dimensionless LCO speed, e is a constant parameter. The nonlinear load is defined as where κ (5) α is another constant dimensionless stiffness parameter. The values of the constants are listed in Table 3.

Bifurcation diagram
Similar to the example shown in Sect. 6.1.1, we compare the fitted bifurcation diagram with the true bifurcation diagram for two sets of parameters corresponding to a subcritical and a supercritical behavior. For the time-spectral analysis, the number of time instances is n = 5. The prescribed motion center point is set to 5 • with a delta of ±4 • . Other motion magnitudes can be selected if the behavior of that point is of interest. Both cases use the non-dimensional parameters shown in Table 3. The parameters used for the subcritical and supercritical cases are shown in Table 4.
The results are shown in Figs. 6 and 7. The fitted curve agrees well with the actual curve for the supercritical case. However, for the subcritical case, there is some deviation. In the region between 1 • ≤ α ≤ 9 • , the fitted curve is close to the actual curve, but for the prescribed motion magnitude α > 9 • there is some deviation. Since our goal is to improve the stability property of one LCO point at α = 5 • , the current method is sufficient. Otherwise, if we want to improve the stability of the whole bifurcation diagram, like the current one, adding more sampling points and using a higher-order polynomial for the curve fitting is necessary. However, the agreement is sufficient for optimization as it captures the overall characteristics of the bifurcation diagram, particularly in the unstable branch. Furthermore, as the optimization progresses, the difference between the fitted and actual bifurcation diagrams should reduce as it approaches a supercritical response, which is well represented by the fitting method.

Optimization
With the bifurcation verified, we are ready to conduct an optimization. The detailed problem parameters are  Table 3 for the subcritical case. (Color figure online) Fig. 7 True (red line) and fitted (blue line) bifurcation diagrams with sampling points (red dots) using parameter from Table 3 for the supercritical case. (Color figure online) defined in Sect. 6.3.2. The optimization problem statement is given in Eq. (45) and is a special form of Eq. (3). The objective is to minimize a function of the dimensionless weight and softening effect with respect to the structural weight and cubic stiffness coefficient, subject to the LCO speed (not flutter speed) and LCO stability constraints. As mentioned before, the dimensionless mass is varied by changing the mass per unit length while keeping the mass of the fluid in a circle with the half chord as the radius fixed. We subtract the stiffness coefficient squared to penalize the nonlinearity of the design. It is because the nonlinear property, such as geometric nonlinearity, is more difficult to predict accurately than the linear property.
by varying m, κ (3) α , subject to g stab (m, κ (3) α ) ≥ĝ, Both constraints are bound, whereĝ is the LCO stability lower bound, and μ is the LCO speed lower bound. In this work, the lower bounds are as follows, The bounds for the design variables are defined as As before, the number of time instances is set to n = 5, corresponding to two temporal modes. This choice was shown to be sufficient by [42] for an airfoil. For more complex problems where the higher modes play a more significant role, the users can always add more modes at a price of computational time. As before, we conducted the optimization using SNOPT [46] with the Python interface provided by pyOptSparse [47].
The contour plot for the two constraints in the given design space is shown in Fig 8. The initial solution is set to be The initial design has a lower objective function value, However, this violates the stability constraint; the optimized solution has a higher objective function value of m ( * ) − κ (3) α ( * ) 2 = 9.11, but this design does not violate the stability constraint. This is evident in Fig. 9, which shows the initial and optimized bifurcation diagrams. Although it has a higher LCO speed, the baseline design is subcritical.

Conclusion
This paper proposes a fitting curve-based method as a metric for LCO stability and an adjoint-based derivative computation that can be used for design optimization. The advantage of this method is that the stability metric is computed with only three points, in this case, solved using a time-spectral method, making the approach computationally efficient for estimating the bifurcation diagram. Further, the computational cost of the adjoint approach is independent of the number of design variables. Three adjoint solutions are required to compute the total derivative. We demonstrated the proposed method in two problems: (1) An LCO stability constrained LCO parameter optimization problem where we found that the ini-tial instability was suppressed and the optimized solution was stable, (2) an aeroelastic optimization problem with LCO parameter and LCO stability constraints, where the baseline design subcritical behavior was successfully suppressed, resulting in a supercritical optimized design. This method can potentially be used to solve large-scale design optimization problems, such as a flutter speed optimization problem with a LCO stability constraint for an aircraft wing.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.