Efficient limitation of resonant peaks by topology optimization including modal truncation augmentation

In many engineering applications, the dynamic frequency response of systems is of high importance. In this paper, we focus on limiting the extreme values in frequency response functions, which occur at the eigenfrequencies of the system, better known as resonant peaks. Within an optimization, merely sampling the frequency range and limiting the maximum values result in high computational effort. Additionally, the sensitivities of this method are not complete, since only information about the resonance peak amplitude is included. The design dependence with respect to the frequency of the extreme value is missed, thus hampering the convergence. To overcome these difficulties, we propose a constraint which can efficiently and accurately limit resonant peaks in a frequency response function. It has a close relation with eigenfrequency maximization; however, in that case, the amplitudes of the frequency response are unconstrained. In order to keep the computational time low, efficient implementation of this constraint is studied using reduced-order models based on modal truncation and modal truncation augmentation. Furthermore, approximated sensitivities are investigated, resulting in a large computational gain, while still yielding very accurate sensitivities and designs with almost equivalent performance compared with the non-approximated case. Conditions are established for the accuracy and computational efficiency of the proposed methods, depending on the number of peaks to be limited, numbers of inputs and outputs, and whether or not the system input and output are collocated.


Introduction
The dynamic behavior of structures is a key aspect of the design process for many engineering applications. A frequency response function (FRF) expresses the amplification of the system under harmonic dynamic excitation, which is a critical aspect of its functionality. In some applications, the maximum response must be limited, for example, in the case of sensitive equipment which needs protection (Tsai and Lin 1994), for limiting acoustic transmission (Fesina et al. 2017), or in systems for high positioning accuracy (van der Veen et al. 2017). Alternatively, a transmission ratio at a resonance frequency might be limited from below, for instance in sensor equipment which needs a minimum response (Xia et al. 2014).
Many engineering applications focus on maximizing eigenfrequencies in order to extend the bandwidth or the operating frequency range. Doing this manually is a timeconsuming and difficult iterative process. With structural optimization methods, this iterative process can be performed automatically. Specifically, topology optimization is a very powerful approach, since no initial concept needs to be given, and a wide range of resulting shapes and layouts is possible (Sigmund and Maute 2013). For the classical problem of optimizing eigenfrequencies, many approaches already have been proposed, see, e.g., Zargham et al. (2016) for an overview. Several methods exist to maximize specific eigenfrequencies (Díaz and Kikuchi 1992;Ma et al. 1995), to create a gap between two eigenfrequencies (Jensen and Pedersen 2006), or to obtain eigenfrequencies close to desired frequencies (Ma et al. 1994). However, resonance occurs when a structure is harmonically excited at the eigenfrequencies, causing extreme responses. This phenomenon has received no attention in all the aforementioned methods.
Besides eigenfrequencies and eigenmodes, the dynamic behavior of a system is also determined by the geometric location and direction of both the input load and the output at which displacement is observed. Some optimization formulations focus on minimizing the vibrational amplitude for steady-state periodic loading (i.e., dynamic compliance) at a specific operating frequency (Ma et al. 1993;Jog 2002;Du and Olhoff 2007a). Alternative objective functions have recently been studied more extensively to improve convergence of these problems, especially for operating frequencies above the first natural frequency, based on the input power (Niu et al. 2018;Silva et al. 2020). Instead of one specific working frequency, others focus on minimizing the frequency response amplitude over a range of frequencies (Ma et al. 1993;Yoon 2010;Shu et al. 2011;Liu et al. 2015). However, for some applications, only the peak amplitudes of a frequency response within a range of frequencies are of importance. The minimization of maximum frequency response amplitude in the entire frequency range (i.e., the H ∞ -norm) is shown for a sizing optimization of a beam model (Venini and Pingaro 2017), but not yet in a topology optimization setting.
The approach of Venini and Pingaro (2017) could be used to limit the maximum value of the FRF, schematically shown in Fig. 1. It relies on an iterative search to obtain the frequency corresponding to the maximum amplification (Bruinsma and Steinbuch 1990). However, this algorithm requires many evaluations of the FRF and the solution of additional eigenvalues, making it infeasible for practical use in large-scale continuum problems. Especially in cases with low damping, where very sharp resonance peaks are present, extra search iterations are needed to obtain the maximum value with sufficient accuracy, adding further to the computational cost. Furthermore, the optimization convergence with this method is slow, because sensitivity information regarding the resonant peak is incomplete. The sensitivities include information on the maximum peak amplitude, but not on the frequency at which the peak is located. Additionally, by only focusing on the Fig. 1 Resonant peaks of the frequency response function are not allowed to be above the indicated limit maximum peak frequency, the multi-modal nature of the peak amplitudes is not captured. This means that whenever the amplitude of one or more peaks are close to the maximum peak amplitude, a small design change could cause one of the other peaks to become the maximum peak, leading to a jump in peak frequency and a non-smooth behavior in maximum peak value. Only using the maximum peak value also limits practical uses such as individually constraining peaks at distinct levels, or upper limits which are a function of frequency (see Fig. 1). This motivates the present study of constraining a finite number of peak amplitudes, instead of considering the global maximum peak. Other implementations of a resonant peak constraint, which can overcome the above limitations, have not been studied, to the best of our knowledge.
For each point in the frequency response function, a different complex-valued linear system needs to be solved, which involves tremendous computation time for large-scale problems, even for only a few frequency points. Additionally, a second linear system needs solving in the case of non-self-adjoint problems, increasing the computational effort even more. Computation time could be saved by using a reduced-order model, requiring at each frequency point only a small system to be solved. Many different methods exist to create reduced-order models tailored for a wide variety of applications, see, e.g., Besselink et al. (2013). In the field of structural dynamics, the most common is modal truncation (MT), where eigenvectors are used to create the reduced-order model. An alternative is modal truncation augmentation (MTA), where a reduction basis consisting of eigenvectors is augmented with correction vectors, to compensate the errors introduced by the removal of higher frequency modes (for detailed information, see, e.g., Rayleigh (1945), Dickens et al. (1997, Craig andKurdila (2006), andBesselink et al. (2013)). Both these methods are suitable for approximating the low-frequency range with high accuracy.
When using a reduced-order model, calculating the frequency response function becomes inexpensive, but to calculate fully consistent sensitivities, the design dependence of the base vectors has to be taken into account (Hooijkamp and van Keulen 2018). These base vector sensitivities are expensive and usually involve the solution of a full linear system per base vector. A possible reduction in computational cost can be achieved by neglecting the sensitivities of the reduction basis. Using such approximate sensitivities, Han (2012) investigates frequency response sensitivities for a Krylov-based reduced-order model and concludes that the sensitivities are still usable although some degree of accuracy is lost. Furthermore, Yoon (2010) reports the approximate sensitivities hamper the optimization process due to their inaccuracy in non-self-adjoint problems. The direct reason for this lack of accuracy has not yet been clarified. The reduction method is usually chosen such that the response, in this case a resonance peak, is most accurate. However, in an optimization, the sensitivities are driving, thus in addition to an accurate response, the accuracy of the sensitivities is equally important.
We propose a constraint which can effectively limit extreme values in an FRF (Fig. 1), where our focus is on weakly damped structures. Using the eigenfrequencies of a mechanical system, which are related to the peak values of the FRF, a lower or upper limit can be set on the resonant peaks. Each eigenfrequency of interest is individually constrained to take care of the non-smoothness problem of the maximum amplitude described earlier, i.e., jumps in frequency corresponding to the maximum resonant peak are not possible. Additionally, having a constraint per resonant peak enables individual peaks to be limited by distinct limits, and the use of frequency-dependent limit functions becomes possible, thus providing a more flexible constraint than the approach of Venini and Pingaro (2017), who only use the maximum resonant peak value. Furthermore, by including the eigenvalue sensitivity information, our sensitivities become consistent with the resonance frequency and thus more accurate. To limit the time spent in calculating the resonance peak amplitudes and their sensitivities, we propose to use reduced-order models with approximated sensitivities. Using the reduction strategies MT and MTA, the accuracy and optimization convergence of the approximated sensitivities is investigated for both self-adjoint and non-self-adjoint problems. The implementation uses density-based topology optimization, but can be applied to other topology optimization approaches as well (Sigmund and Maute 2013). For clarity, we will limit ourselves to the single-input single-output (SISO) case, but the method is also extendable to multiple-input multiple-output (MIMO) cases.
The paper is organized as follows. In Section 2, the considered finite element model is introduced, followed by the definition of the optimization problem. The proposed constraint is explained in detail, by either solving full systems or by using model reduction techniques (MT and MTA). Additionally, for the reduction methods, both the consistent and approximate sensitivity calculation is described. Section 3 studies the different implementations using both self-adjoint and non-self-adjoint test cases. The performance of the approximated sensitivities are inspected and also their effect on the optimization is shown. Next to that, results of some variations in limit functions are given, to show potential in practical use.

Dynamic response modeling
Working towards a model suited for topology optimization, we first establish the considerations used regarding design parametrization and secondly, the numerical modeling of the dynamic response. Since density-based topology optimization is used, each element's elasticity and density is scaled continuously according to the design variables s between s min (void) and 1 (solid) (Bendsøe and Sigmund 2003). First, the design variables are filtered using a spatial density filter, resulting in filtered design variables ρ (Bruns and Tortorelli 2001). To force the optimizer to a clear solid/void design, intermediate density variables are penalized using scaling factors for the element matrices of stiffness κ and mass μ, respectively, as This penalization was investigated by Zhu et al. (2009) in order to prevent low-frequency eigenmodes of void areas, which often hamper topology optimizations using eigenfrequencies. For the scaling of stiffness, it uses a combination of a linear term (weighed by w) and a part with exponent p.
For the discretization, we use bi-linear quadrilateral finite elements, a 2 × 2 Gaussian quadrature, and assume a plane strain condition. The stiffness and mass matrices are assembled, respectively, by with the element matrices denoted K el and M el .
We introduce damping in the form of structural damping (i.e., hysteresis), often used in airplane vibrations and flutter analysis, which is proportional to displacement. Effectively, a damping factor of η is used to create a complex stiffness (Craig and Kurdila 2006). This kind of damping does not change the frequencies at which the peak amplitudes occur, which means the undamped eigenfrequencies can directly be used. A viscous damping such as Rayleigh or modal damping could also be used without adding computational effort, as the damped eigenfrequencies correspond to the peak amplitudes in that case, which can be calculated as a simple correction on the undamped eigenfrequencies.
Using a steady-state SISO system with harmonic inputs and outputs for the sake of simplicity, the discretized Ndimensional frequency-domain system of equations is where u denotes the state vector capturing the displacements and deformations of the entire structure. The input vector b describes the spatial distribution and direction of the unit input force, and the output vector c describes that of the observed unit displacement. The unit input vector is scaled with the input force q and the resulting output displacement is denoted y, both dependent on frequency ω. For a derivation, it is referred to any dynamics text book, such as Craig and Kurdila (2006). We can write this into a complex frequency-dependent transfer function G(ω), denoting the transmission between input force and output displacement, with Z(ω) the complex symmetric N × N frequency dependent dynamic stiffness matrix. The magnitude 1 of this function |G(ω)| is used to obtain the amplification of harmonic amplitudes from input to output, possibly scaled to decibel , denoted as |G(ω)| dB . The undamped eigenfrequencies Ω i and eigenvectors φ i of the system can be calculated by solving the generalized eigenvalue problem for which the eigenfrequencies are ordered as 0 ≤ Ω 1 ≤ Ω 2 ≤ . . . ≤ Ω n for the n eigenfrequencies of 1 The notation |•| means to take the complex norm or magnitude of the value • interest. The eigenvectors Φ = φ 1 , φ 2 , . . . , φ n are massorthonormalized according to Φ T MΦ = I. Since structural damping is used, the peak frequencies are equal to the undamped eigenfrequencies ω i = Ω i , at which the FRF amplitude |G(ω = ω i )| reaches its extreme values.

Optimization problem formulation
A general optimization problem involving resonance peak constraints can be formulated as The proposed constraints can be used as either an upper limit (7) or as a lower limit (8) for the response at any peak frequency, provided these frequencies are known from eigenvalue calculation. Any peak frequency (ω i ∀ i ≤ n) can be placed in subsets S upp or S low , or both. Additionally, the formulation is not limited to one single upper and lower limit function (g upp and g low ), e.g., the first resonance peak could be limited differently than the second. For robustness against mode switching, a mode tracking strategy (Kim and Kim 2000) is advisable to ensure continuity of the constraints.
Any reasonable choice of objective function f is possible, but in this paper, we will limit ourselves to an objective function in the form of a mean eigenvalue (Ma et al. 1995). To maximize n eigenfrequencies, the harmonic mean of those frequencies is taken as objective This objective function is relatively insensitive to mode switching, which otherwise could introduce discontinuities if not taken into account correctly (see, e.g., Kim and Kim (2000); Du and Olhoff (2007b)). In order not to obscure the scope of the paper, we choose to avoid using these techniques. Additionally, this objective helps preventing trivial solutions: in case n peaks corresponding to the n lowest frequencies are limited from above, a possible trivial solution would be to create n artificial rigid body modes (by means of disconnected parts or mechanisms), which have no effect on the point of interest and a very low transmission ratio, but do replace the lowest eigenfrequencies.
In order to prevent ill-conditioning of the system matrices, the lower bound on the design variables is set to s min . Secondly, a volume constraint is imposed to prevent trivial all-solid solutions. In total, this leads to the following optimization problem which is considered throughout this paper: For further use, we abbreviate the frequency response value at G(ω i ) as G i . In subsequent sections, three different methods are proposed to calculate the peak values |G i |. All methods require the eigenpairs (Ω i , φ i ) to be calculated beforehand.

Full method
The most straightforward method to calculate the FRF amplitudes at each required peak frequency, by solving the full linear system: The sensitivities of this function with respect to the filtered design variables are where u i is the state vector containing the solution of the harmonic response, and ξ i the adjoint vector at each peak frequency ω i . The last term in (17) adds the sensitivity information with respect to the peak frequency (and thus the eigenfrequency), which cannot be included in any method which iteratively finds the peak value in an FRF (such as Venini and Pingaro (2017)). Both the state and adjoint require a full complex-valued system to be solved, where the state is only dependent on the input vector b, and the adjoint depends on the output vector c. Hence, the importance of the output location on the sensitivities is explained by the adjoint system having the output vector as a right-hand side. Note that these equations could be solved using one matrix factorization. In case an iterative solver is used, the systems would have to be solved separately if b ∝ c.
To complete the sensitivity calculation, the derivatives of the interpolation functions (1) are The peak frequency sensitivities are equal to the undamped eigenfrequency sensitivities in case of structural damping, and are calculated by which does not require the solution of any extra linear systems (a derivation is found in, e.g., Haftka and Gürdal (1992)).
The peak responses of a damped dynamic system are complex values and so are their sensitivities. To obtain the magnitude of the frequency response, the complex norm is taken as 2 The sensitivity of the complex norm is calculated as resulting in a real-valued sensitivity. In case a transformation to decibel is used, the response and its sensitivity can be calculated as Finally, the sensitivities are treated with the density filter as described in the work of Bruns and Tortorelli (2001). These last five differentiation operations in (19)-(23) (material interpolation, eigenfrequency, complex norm derivatives, decibel transformation, and filter) are identical in the following methods using reduced-order models.

Modal truncation
In order to reduce the time spent in calculating all the required frequency response values, model reduction techniques can be used. Although it is very expensive to compute the eigenvectors, still the solution of the dynamic systems of equations will contribute significantly to the total computation time. The main reasons for this are twofold. First of all, the eigensolver only requires solution of real-valued matrices, while the resonance peaks involve a complex-valued matrix to be solved, which can be compared with a real-valued matrix of size 2N × 2N. Secondly, each peak requires an unique system of equations to be solved, while an eigensolver uses only one system of equations to iteratively converge towards the eigenpairs. When using a direct solver, this means that a factorization is required for each resonance peak in the optimization problem, while only one factorization is enough for the eigensolver. In case of an iterative solver, the same could be said about the preconditioner.
By using the eigenvectors which are already computed for the objective, modal truncation can be applied to obtain smaller (n << N) reduced system matricesK andM ∈ R n×n , and input-output vectorsb andc ∈ R n (see, e.g., Craig and Kurdila (2006)). Thus, the higher modes of the system are truncated, as is assumed that these do not greatly affect the lower frequency spectrum. The reduced system is obtained by a Galerkin projection of the full system matrices on all the known eigenvectors Φ, as Since the matrices are projected on the eigenvectors, the resulting system matrices are diagonal. The diagonal terms of the reduced stiffness matrix becomeZ kk (ω) = Ω 2 k (1 + iη) − ω 2 , which makes the response calculation very efficient using modal superposition. Again the frequency ω is chosen as peak frequency ω i , resulting in the reduced response as 3 In order to calculate the sensitivities of this function, two approaches are proposed. The first one is the consistent calculation of sensitivities, and the second method is an approximation of the sensitivities by ignoring the design dependence of the reduction basis.

Consistent sensitivities
The sensitivities for the consistent method become more involved than the full model sensitivities, since a reduction step has been added. The sensitivities now have to be calculated using in which the adjoint eigenvector sensitivity of dφ T k ds j needs to be calculated. The term involving the eigenvector derivatives is the sensitivity with respect to the reducing basis. By solving the adjoint saddlepoint problem the adjoints ν k and α k can be calculated. For a detailed explanation, the reader is referred to Lee (2007).
After solving for adjoints, the eigenvector sensitivities can be obtained as Note that for the calculation, a factorization of the matrix introduced in (27) is needed for each eigenvector in the base, or one iterative solution per eigenvector, for each peak to be observed. Instead of solving complex systems, now real-valued matrices can be used, which saves considerable computation time. Additionally, the number of linear systems to be solved is reduced by a factor of two in case the problem is not self-adjoint and an iterative solver is used.

Approximate sensitivities
From the exact sensitivities in (26), it can be seen that the first term is divided by the modal stiffness (Z kk (ω i ) = Ω 2 k (1 + iη) − ω 2 i ), and the second term is divided by the modal stiffness squared. The sensitivities are largest when the dynamic stiffness is very small (ω i ≈ Ω k ), causing the second term in (26) to become dominant, since it is squared. Therefore, in order to reduce computational effort, we propose to approximate the sensitivities by ignoring the first term containing eigenvector sensitivity terms with respect to input b and output c. Effectively, this means that the design dependency of the reduction basis in (24) is not considered and taken as constant, when taking the design sensitivities of (25). This results in where the reduced state and adjoint are now calculated using the reduced model by solving This method does not require any solution to the full linear system at all, but it may hamper convergence, because it is an approximation and information regarding the input and output is only contained via the eigenvector projection. This means that the sensitivities do not contain information anymore about the parts of b and c orthogonal to the eigenvectors in the basis, which was previously included in (28). In case the basis does not change with respect to the design (i.e., the eigenvectors do not change), these sensitivities are exact. Additionally, the above approximation makes implementation very easy, since the sensitivity of the approximated MT-based method (29) strongly resembles the sensitivity of the full method (17), with the following substitutions: which are simply the projections of the approximated problem. For comparison, Table 1 shows an overview of the number of full system solutions required for the different methods.

Modal truncation augmentation
MTA is an established concept of model reduction in the field of structural dynamics (Dickens et al. 1997). By augmenting the reduction basis with extra correction vectors, the reduced model becomes more accurate. The correction vectors add localized information, which was lost by removing the high-frequency content in MT. By extending the basis with specific local information, not only the response becomes more accurate but also the accuracy of sensitivities might be enhanced. Instead of only augmenting the solution to the input force b, we choose also to add a correction with respect to output vector c.
Since the adjoint vector ξ is determined by solving the system using the output vector (18), this should improve the approximated sensitivities. The reduction basis Φ is extended with undamped linear solutions v 1 and v 2 of the input and output vectors, at a shift frequency σ < Ω 1 , to obtain the augmented reducing basis Note that for a collocated system (b ∝ c), only one vector needs to be added as both vectors would be linearly dependent. In the case of MIMO systems, it is trivial to add more vectors for all distinct inputs and outputs (Dickens et al. 1997).
The augmentation vectors are orthonormalized with respect to the other vectors, by solving a small eigenvalue problem to diagonalize the matrix where Ψ = Φ v 1 v 2 , Q is an orthogonal matrix containing the eigenvectors, and Λ a diagonal matrix with the eigenvalues of the un-orthogonalized projected mass matrix. Using these, we can obtain a linear combination to get the orthonormal system The reducing basis thus becomes n · n p 0 0 0 Many in-and outputs, many peaks MTA n io + (n + n io )n p 0 n io 0 Few in-and outputs, many peaks In case direct solvers are used, the number of matrix factorizations is equal to n io = 1 and n p = 1, except for consistent MTA, which requires 1 + n factorizations. All methods require the solution of n eigenvalues and eigenvectors of the full system and the reduced model can be obtained as Note that the reduced stiffness matrixK now is not diagonal anymore. Therefore, we write the dynamic stiffness matrix asZ(ω) =K(1 + iη) − ω 2M . The peak values can be calculated using It is possible to calculate consistent sensitivities of this reduced model; however, it involves a lengthy derivation which is omitted here for the sake of compactness, and placed in the Appendix. In contrast, the approximated sensitivities, where we neglect the design dependence of all base vectors, are just as straightforward to derive as with the MT method. They can be calculated using which is merely a change of basis as compared with the approximated sensitivities of the MT method (29). The reduced linear system now involve a matrix of size (n+2)× (n + 2): This last method is not as inexpensive as the approximated sensitivity MT method, due to the additional solutions required to augment the basis, but also not as expensive as the full method or the consistent sensitivity MT method (an overview is given in Table 1). It only requires one extra matrix factorization, or one iterative linear solution per unique input and output. Since we only consider the SISO case, only one solution (collocated system) or two solutions (non-collocated system) are required. But in the case of MIMO, each unique in-and output vector of interest would have to be augmented, each requiring a linear solution (Dickens et al. 1997).

Results
Two test cases are used to study the optimization process and the influence of various model reductions. First of all is a cantilever problem (Fig. 2a), with a solid non-design region in the middle of the domain, and a collocated input and output. Secondly, a free-floating stage (Fig. 2b) is used, with non-collocated input and output, where the input is distributed along two square non-design domains which represent actuators. The output location is in the center of the top surface. Since the structure is free floating, it exhibits three rigid body modes. For the MTA approach, this means that the augmentation vectors cannot be static (σ = 0). Instead a quasi-static solution is obtained at σ = 50 rads −1 , which is well below the first eigenfrequency and will be used for both test cases. The objective involves n = 3 structural eigenvalues, which are also limited in amplitude (n p = 3). For simplicity, no mode tracking is applied in all examples. Further parameters used for the optimization are listed in Table 2.  The optimization problem is solved using the method of moving asymptotes (Svanberg 1987), of which the number of iterations is limited to a maximum of 200. All the results presented converged to a stationary solution, unless otherwise mentioned. Furthermore, the objective function is scaled to be − 100 in the first iteration and the volume constraint is scaled by a factor 10.
First, optimization results of the different peak constraint implementations with the consistent approach are given in Section 3.1. After that, in Section 3.2 a comparison is made between consistent and approximate sensitivities. Finally, optimizations using approximate sensitivity information are shown in Section 3.3.

Consistent optimization
To show the operation of the proposed peak constraint, the test cases are optimized using the full method and the two reduced-order models using consistent sensitivities. For reference, the results of an unconstrained optimization (i.e., only constrained in volume) are also shown.

Cantilever
The cantilever problem is constrained with an upper limit of g upp = − 1 dB for the first three peak amplitudes. This value is arbitrarily chosen here: it is physically achievable, and this limit will cause the constraint to be active. In practice, the designer would provide a limit based on operational targets. Using the method involving full system solutions for optimization (Section 2.3), we obtain the design shown in Fig. 3b. The reference design of an optimization without any peak constraints is shown in Fig. 3a.
An overview of the performance values is given in Table 3 and FRFs of the designs are shown in Fig. 4. From this, it can be seen that the peaks are indeed limited with g upp = − 1 dB. The lower peak values come at a cost, because the first two eigenfrequencies are significantly lower than in the design without peak constraints. This is reflected in the higher objective function, in which the lowest eigenfrequencies contribute most. When looking at the mode shapes in Fig. 5, it is evident that they are different from the reference case. At the left side of the constrained design, a mechanism can be recognized, which rotates the The number of inputs and outputs n io = 1 and n p = 3. Evaluated on the full model Fig. 4 FRFs of the different cantilever designs right part of the structure such that the tip displacements are reduced.
As for the results involving reduced-order models (Sections 2.4-2.5), the final designs of the MT and MTA method are respectively shown in Fig. 3c and d. These designs are hard to distinguish from the design obtained using the full method (Fig. 3b), and their performance is equivalent (Table 3). Looking at the convergence history of the cantilever optimization, as shown in Fig. 8a, it can clearly be seen that the three methods have similar convergence. In our Python implementation, using matrix factorizations whenever possible, solving the required eigenvectors takes about 13 s, one real-valued factorization 0.85 s, and a complex factorization 2.5 s. This means that factorizing three complex matrices, required for the full method, represent a significant portion (about a third) of the computational time. Using the reduced-order models, only 3 real factorizations for MT or 4 for MTA have to be made, hence the time saved per iteration (Table 3). In case of iterative solvers, this computational gain is more debatable, as the adjoints cannot be calculated with simple back-substitutions anymore and separate iterative solutions are needed. Although beyond the scope of this work, it should be noted that model order reduction with consistent sensitivities could become more viable using aggregation strategies. Aggregation already has been implemented successfully for instance in stress constraints (Yang and Chen 1996) to reduce computational time for the sensitivity calculation, and in eigenvalue optimization (Torii and de Faria 2017) to overcome differentiability issues. In our case, for the full solution strategy, each peak constraint is dependent on a different Z(ω j ), which requires an adjoint to be solved for each of the peaks. When a reducedorder model is used, the adjoints would not be expensive, as they are evaluated on the reduced modelZ(ω j ). The expensive adjoint solutions (in this case the model reduction basis sensitivities) are involving identical system matrices for each peak constraint, allowing the expensive sensitivities only to be calculated once when aggregating. This effectively results in only n linear solutions on the full system for MT and 2n io + n linear solutions for MTA, independent of the number of peaks considered, compared with the full method, still requiring the solution of 2n p full complex linear systems if the constraints would be aggregated. In Table 1, the resulting number of solutions on the full system for an aggregated constraint could be seen as n p = 1 for MT and MTA, but not for the full method.

Stage
The second example involves the optimization of a freefloating stage (Fig. 2b). In comparison with the cantilever case, this example has a non-collocated input and output vector. Additionally, there are rigid body modes present in this example, which means that the peaks of the 4th, 5th, and 6th eigenmode (the first three flexible modes) are constrained. We choose the constraint limit as g upp =-25dB.
The resulting design of the optimization without peak constraint (Fig. 6a) and the constrained designs (Fig. 6bd) again shows a trade-off between peak limitation and eigenfrequency values (Table 4). The peak limits are attained at cost of lower eigenfrequencies.
Between the designs resulting from the full method (Fig. 6b) and the consistent MTA (Fig. 6d), the difference in design is hardly recognizable. The design resulting from consistent MT (Fig. 6c) is different although its performance is equivalent to the other designs (Table 4 and Fig. 7), indicating a different local optimum. Also here, when looking at the convergence history in Fig. 8b, the use of reduced-order models with consistent sensitivities does not hamper optimization convergence and the convergence is very similar.

Comparison of consistent and approximate sensitivities
As already observed in previous section, the use of consistent sensitivities in optimization with reduced models yields comparable results to using a full model. However, using reduced-order models with consistent sensitivities in the previous examples results only marginally increased computational efficiency compared with the full method. Therefore, we investigate the effect of approximating the sensitivities of the reduced-order models, by ignoring the design dependency of the model reduction basis.
The effect of neglecting the reduction basis sensitivities, therefore not requiring any expensive solution of the adjoint problem, is visually demonstrated in Fig. 9. This figure shows the sensitivities of the first two peak constraints of the stage case in the first design iteration, thus consisting of a uniform density field. It is clear that the localized details Fig. 7 FRFs of the different stage designs around the output location (Fig. 9a, d) are not present in the approximate sensitivities (Fig. 9b, e). Identical observations could be made for the third peak constraint, not shown here. The local features are present again in the approximate sensitivities of the MTA method (Fig. 9c, f), where the additional vectors provide this information.
To quantify the error between approximate and consistent sensitivity fields, we introduce a sensitivity error norm as Again, the sensitivities are evaluated in the first design iteration for both the test cases. In Table 5, the error values are reported for both the cantilever and the stage case. Two observations can be made from these values. First of all, the sensitivity error is much lower for the MTA sensitivities than for the MT. Second, the sensitivity errors for the peak constraints go up for higher eigenfrequencies when using MTA. This might be related to the choice of shift frequency σ (chosen below the first eigenvalue) to evaluate the augmented response.

Optimization with approximate sensitivities
The effect on the optimization process when approximating the reduced-order model sensitivities is demonstrated in this section. Designs and performance are compared between the consistent and approximate approaches. Fig. 9 The exact sensitivities (a), and the approximate sensitivities using MT (b), of the first peak. Consistent and approximate sensitivities of the second peak, respectively (d) and (e). The approximate sensitivities of the first two peaks using MTA are shown in c and f

Cantilever
Starting again with the cantilever case, the results of optimization with approximate sensitivities are shown in Fig. 10b and c for MT and MTA, respectively. Especially the design optimized with approximate MTA is very similar to the consistent design (Fig. 10a). This can also be seen from the performance values in Table 6. The objective and eigenfrequencies of the consistent design and the approximated MTA design are very similar. However, the optimization with approximated MT did not even converge, as a feasible design was not reached (volume constraint violation). This can also be seen in the iteration history (Fig. 8a), which shows that the approximate MTA follows a similar path compared with the consistent methods, while the approximate MT completely different path as the optimization progresses. The introduced sensitivity approximations achieve further computational gain. The timings in Table 6 show that the approximate MT method requires virtually no extra time to calculate the peaks and their sensitivities, compared with the reference case without peak constraints (Table 3). The approximate MTA method saves about a quarter of total computation time, which means that the time required to calculate the peak values and their sensitivities is shortened almost an order of magnitude (8 s for the full method  versus 1.2 s for approximate MTA). When using iterative solvers, the computational gain might even be larger, since no iterative solutions are required for the sensitivities.

Stage
The optimization results of the stage case are shown in Fig. 11b and c for the MT and MTA methods using approximate sensitivities. For comparison, the design obtained from the consistent optimization (Fig. 11a) is also shown. Again, the design resulting from approximate MT is clearly distinct from the other designs. In contrast to the cantilever, the resulting MT stage design is feasible (Table 7), although it took 82 iterations instead of 15, which the other methods required. Again, the convergence of the approximate MTA is very similar to the other methods (Fig. 8b), while the approximate MT method causes a completely different convergence path and a worse optimum.

Case variations
In this section, variations of the limit function are explored to gain more insight into the behavior and possibilities of the proposed constraint. Both the full and approximate MTA methods are used for the optimizations.
To illustrate the individual control of peaks in the FRF, the first variation is to choose a lower limit instead of an upper limit on the cantilever case. For the first peak, we now use |G 1 | dB ≥ 10 dB. The other two peak constraints are  kept on their original upper limit of g upp = -1 dB. The resulting designs and FRF are shown in Fig. 12, from which can be seen that the requirements are fulfilled. Although the design optimized with approximate MTA is asymmetric, it performs a little better (higher eigenvalues) than the design optimized with the full method. The bulk of the material is distributed in the same manner for both designs, with the main difference the slender structure being removed at the top part for the approximate MTA design. Also non-constant limit functions can be used, which is demonstrated using the lower limit g low = −10 log(ω) + 20.
With a logarithmic frequency axis, it represents a sloped straight line in an FRF. Any other user-defined peak envelope function can also be used for the constraint. Note that the sensitivity of this function also needs to be taken Fig. 12 Results of choosing a lower limit for the first peak |G 1 | dB ≥ 10 dB. The design from the full method (a) and the approximate MTA method (b) and their FRF (c) into account as this is a function of frequency. The resulting designs are feasible and all three constraints are active, as is shown in Fig. 13. Both optimized designs are very similar to each other, with nearly identical dynamic behavior.
Next to a frequency-dependent constraint limit, another possibility is an adaptive limit. Instead of maximizing the eigenfrequencies, the maximum peak can be minimized. Practical implementation of this min-max problem can be done using a bound formulation, which results in the following optimization problem: The results of this optimization problem are shown in Fig. 14. The obtained designs feature appendages near the tip, that are weakly connected to the main structure. These appendages add low-frequency modes to the structure, that do not result in a large amplitude at the output point. This is advantageous in this example because the response of only the three lowest peaks is limited. It is clear that the optimizer exploits not having a penalization on low eigenfrequencies, by adding these artificial low-frequency modes in the process. This demonstrates the need for additional requirements on the optimization problem in order to obtain a meaningful design.
In practice, operational conditions might impose other requirements on the FRF. For example, in equipment operating at a constant frequency, a low response amplitude at exactly that frequency is desired. This can be achieved by extending the optimization problem in (43) with the constraint G(ω op , s) dB ≤-80, where ω op is the operational frequency, chosen as 300 rad s −1 , initially just in between the second and third natural frequency. Note that extra bounds on the FRF can be set without adding any computational expense when using approximate MTA.
The results of the optimization with this extra constraint can be seen in Fig. 15. The full design satisfies the operating frequency constraint with a value of − 80.01 dB at ω op , whereas the design optimized with approximate MTA is just infeasible with a value of − 79.87 dB evaluated on the full model, while it is feasible on the reduced model. This might be explained by the loss of accuracy when using reducedorder models. Especially when the response is close to zero, a small absolute error might introduce a large error in the decibel scale. For the full design, the optimum design is found with β = − 2.00 dB, while for the approximate MTA design β = − 2.40 dB. Fig. 15 Results a maximum peak minimization subject to an additional maximum amplitude constraint at 300 rad s −1 (a, b) and their corresponding FRF (c) In this work, we proposed a constraint to limit peak frequency response amplitudes. It is able to effectively and efficiently limit peak values of an FRF. By using the eigenvalues and eigenvectors of interest, versatile upper and lower limits can freely be selected per resonance peak.
Various ways to reduce the computational effort have been explored using three different implementations of this constraint, one based on solving the full system, and two involving reduced-order models (modal truncation (MT) and modal truncation augmentation (MTA)). It has been shown that by using consistent sensitivities, thus by including sensitivities of the reduction base vectors, the optimization converges to nearly identical designs with equal performance. In this case, no complex linear systems need solving anymore, but a larger number of real-valued linear solutions are required.
The sensitivities of the reduced models can be approximated by ignoring the design dependence with respect to the model reduction basis. Only having eigenmodes in the basis (MT) leads to inferior convergence during optimization and infeasible design. It is shown that the approximation in case of MT fails to reveal local details in the sensitivity field, since higher frequency eigenmodes are truncated. By augmenting the basis with correction vectors (MTA), important local details, which were previously truncated, are resolved in the sensitivities and results similar to the consistent optimization are obtained. Sensitivity approximation has two advantages, the first being implementation ease. The sensitivities directly use the projected reduced solution for evaluation, thus easy to implement in code which already uses sensitivities based on the full model. Secondly, fewer real-valued linear solutions are required as the number is not dependent on the number of peaks anymore, but dependent on the number of augmented vectors. In the examples shown, the time required to calculate the peaks and the associated sensitivities can be reduced by almost an order magnitude when using the proposed techniques, with a possibility of even further time reduction when iterative solvers are used.
Looking forward, an interesting research direction is the aggregation of peaks. Combining reduced models with aggregation of the resonant peaks, the number of full solutions can be reduced, which is not possible using the full method. Theoretically, the number of linear solutions could be independent of the number of peaks, while still using consistent sensitivities. This is especially interesting in cases where many peaks are constrained, or when many inputs and outputs are present (MIMO), making both directly evaluating the full system and the MTA method computationally expensive. In this case, MT consistent optimization might be very suitable. It needs to be noted that, depending on application demands, the damping model requires improvement for more accurate results. In order to improve accuracy, it is recommended to consider modal damping or even Voigt/Maxwell-type models, depending on the constituent material. Related to this, an interesting research direction is a multi-material design in a peak limitation context, where multiple materials with different damping properties can be placed (see, e.g., van der Kolk et al. (2017)). In this manner, the amplitude of a peak can be changed both by changing the eigenmodes (as in the current work) and by effective placement of material with different damping properties.

Replication of results
We are confident that the paper contains sufficient detail on methodology and implementation, such as in Table 2, that the results can be re-created. In case of questions or difficulties, readers are welcome to contact the authors. Additionally, raw data of the resulting designs and their FRFs can be found in the supplementary material. in which the operator ⊗ stands for an outer product A = u ⊗ v or in index notation A ij = u i v j . We know that the reduced dynamic stiffness matrix is written asZ i =K(1 + iη) − λ iM . Note that the previous calculated derivatives of G ĩ Z i are complex valued, but the derivatives with respect toK,M, and λ i should be real valued (van der Veen et al. 2015). Thus, we write

A.2 Orthogonalization
Without explicitly substituting dG i dV , the sensitivity with of the orthogonalization is calculated next. Observing that the basis is orthogonalized by solving a small eigenproblem as V = Ψ QΛ − 1 2 , or with index notation , the sensitivities with respect to the orthogonalizing eigenpairs (Q kj , Λ jj ) are determined as We know that Ψ T MΨ = QΛQ T , which is an eigenproblem. We know that the eigenvectors and eigenvalue sensitivities can be calculated by first solving the adjoint problem (from, e.g., Lee (2007)): Or for a multiplicit set of eigenvalues, also described by Lee (2007),

A.3 Response sensitivity
The final step is to calculate dG i dK and dG i dM , which has contributions from the full eigensolve Kφ j − λ j Mφ j = 0, the augumented vectors (K − σ 2 M)v 1 = b and (K − σ 2 M)v 2 = c, the reduced modelK = V T KV andM = V T MV. Additionally, the M matrix is also used in Ψ T MΨ . Just like the small eigenproblem, the eigenvectors and eigenvalue sensitivities can be calculated by first solving the adjoint (from, e.g., Lee (2007)) Or for a multiplicit set of eigenvalues subjected to KΦ = λMΦ and Φ T MΦ = I, where Φ = φ 1 , φ 2 , . . . , φ m , also described by Lee (2007), The contributions of the eigenpairs (λ j , φ j ) to the sensitivities are and The contributions of the augmented vectors can be calculated by solving the adjoints (K − σ 2 M)ξ 1 = − dG i dv 1 and (K − σ 2 M)ξ 2 = − dG i dv 2 , (A.20) in order to get and The contributions of the reduced-order model are simply and The final contribution to the mass matrix coming from the orthogonalization step is .25)