# Robust Assignment of Natural Frequencies and Antiresonances in Vibrating Systems through Dynamic Structural Modification

Authors: Roberto Caracciolo; Dario Richiedei; Iacopo Tamellin
Journal: Shock and Vibration
 (2021)

---

## Abstract

This paper proposes a novel method for the robust partial assignment of natural frequencies and antiresonances, together with the partial assignment of the related eigenvectors, in lightly damped linear vibrating systems. Dynamic structural modification is exploited to assign the eigenvalues, either of the system or of the adjoint system, together with their sensitivity with respect to some parameters of interest. To handle with constraints on the feasible modifications, the inverse eigenvalue problem is cast as a minimization problem and a solution method is proposed through homotopy optimization. Variables lifting for bilinear and trilinear terms, together with bilinear and double-McCormick’s constraints, are exploited to provide a convexification of the problem and to boost the attainment of the global optimum. The effectiveness of the proposed method is assessed through four numerical examples.


---

## Body

## 1. Introduction

### 1.1. State of the Art and Motivations

Eigenstructure assignment (EA) has been a core research area in structural dynamics. Several direct and inverse dynamic structural modification (IDSM) methods have been proposed up to now in the literature. In particular, inverse approaches are the most fascinating and challenging for researchers since those require to compute the structural modifications that enable to obtain a prescribed eigenstructure for the modified system [1].EA has been tackled through passive, active, and hybrid approaches (see, e.g., [2–4]). Active approaches exploit the energy introduced in the system by external actuators based on the information of some sensors to perform the EA. The state of the art in active control spans from model-based techniques (see, e.g., [5–9]) to receptance-based methods (see, e.g., [10–13]).Passive approaches are widely investigated since they are stable due to the symmetric and definite positive nature of the structural modifications. Several works have been proposed in the literature to address different problems, such as the assignment of frequencies [14–18], vibration nodes [19], mode shapes and natural frequencies [20, 21], antiresonance frequencies [22–24], or the optimal design of dampers [25]. These approaches are useful both in the design of new systems and in the modification of existing ones, to meet the desired dynamic performances. For example, assigning natural frequencies is useful in resonators [26] to ensure large displacements with small excitation forces, or to place resonances far from excitation sources whenever damping them is not possible or very difficult [27]. Assigning antiresonances, i.e., frequencies where the system experiences negligible steady-state force oscillations at some particular coordinates when excited by harmonic excitations, is useful to reduce unwanted vibrations [22]. Assigning the mode shapes is useful to vibration confinement [3] or to shape the forced response in resonators [26].Both analytical and numerical methods have been proposed. The latter ones are usually adopted to solve more complex and general problems that cannot be handled analytically due to their highly nonlinear nature or to the absence of exact solutions. Hence, approximate solutions should be found through minimization-based formulations. For example, the partial assignment of mode shapes and natural frequencies has been performed by the authors and their co-workers in [20, 26, 28] by solving a constrained nonlinear nonconvex optimization through the homotopy optimization technique.A limitation of most of the proposed techniques to EA is that robustness is not addressed, and hence, it is neglected how the assigned natural frequencies are affected by perturbations of the system parameters about the nominal ones. On the other hand, making the performances minimally sensitive to parameter uncertainty is always of paramount concern in optimal design of machines and structures, leading to the idea of “robust design” [29]. Robust design is widely adopted for the optimal design of the strength of materials and components (such as static problems, buckling, or other structural applications) where the goal is trading between oversizing and ensuring reliability in the range of possible parameter variability in “everyday service conditions” [30], or in the worst-case scenario. In contrast, in the case of EA, operating with worst-case scenario is not reasonable, because it would lead to mistuned eigenstructures in the more likely set of system parameters.In the frame of EA, robustness is often interpreted as finding the optimal eigenstructure, or just the optimal natural frequencies, ensuring adequate performances even under perturbed conditions, and some ad hoc applications have been proposed. For example, in [31, 32], the robust design of brakes with respect to self-excited vibrations has been oriented to split the eigenfrequencies as much as possible by maximizing the difference between consecutive eigenvalues. Robust mechanical design to prevent self-excited vibrations due to frictional forces in screw jacks has been exploited in [33] through the analytical analysis of the eigenvalues of a two degrees of freedom (DOFs) model. Robustness, with respect to uncertainties in the system parameters as well as in the disturbance input force, in the design of tuned mass dampers for vibration absorption is often handled by ensuring a large range of frequencies about the imposed antiresonance where the receptance from the disturbance to the response has a low gain (see, e.g., [34, 35] and the references therein).Robust EA has been sometimes discussed in the case of active approaches that are out of the scope of this paper. Indeed, the control force exerted by the actuators can be changed based on the sensed (or estimated) plant modifications. This possibility, at the cost of a more complicate system (that requires actuators, sensors, and a controller) makes usually active EA more robust than passive EA [13, 36].A different goal, that is the one pursued in this work, is ensuring that the imposed natural frequencies and antiresonance frequencies are robust with respect to parameter changes, e.g., due to changing loads, or to uncertainty on parameters that are not exactly known at the design stage. This goal is of practical interest in structures and machines in the presence of either wanted or unwanted vibrations. For example, ensuring insensitive natural (or resonance) frequencies is of great importance in the design of resonators (such as feeders [26, 37] or sonotrodes [20, 38]), where the resonance frequency should match the excitation frequency regardless of the changes in some parameters, to ensure large motion with small actuation effort. In the case of vibration absorption of harmonic excited systems, assigning antiresonances at the excitation frequency is an effective approach ([39–41]), that should be ensured even in the case of changes of some system parameters. Generally speaking, any application where passive EA is performed can benefit of this approach, since increasing robustness allows to get rid of small parameter changes, compared to the nominal model, of uncertainty in the nominal parameters.An effective nonprobabilistic way to incorporate such a robustness requirement into design is using sensitivity, that is useful whenever probabilistic distributions of uncertain parameters are not available and when the optimal design should be, first of all, well-tuned about a nominal condition (i.e., the more likely system configuration).The use of sensitivity as a measure of robustness in EA has been proved to be effective in the case EA through active control, in particular when higher-rank control allows for more combinations of gains leading to the desired eigenstructure. In [42], the eigen-sensitivities with respect to measurement errors in receptances are adopted to increase robustness in state-feedback robust pole placement, by proposing a receptance-based formulation. In [43], robust pole assignment through state-feedback of friction-induced asymmetric systems is made by assigning both eigenvalues and their eigen-sensitivities with respect to the friction coefficient. In [44], eigen-sensitivity to parameter variations in second-order linear systems under state-feedback and state-derivative feedback is discussed and assessed. In [45], the problem of robust eigenstructure assignment with respect to uncertain parameters in the fitted receptances is solved, by minimizing sensitivity with respect to the fitting parameters. In [46], it is developed an approach to increase the robustness of the closed-loop eigenvalues based on the notion of spectrum sensitivities, for both state-feedback and state-derivative feedback control.In the case of EA through IDSM, a few works consider sensitivities: indeed, incorporating robustness with respect to parameters changes is not straightforward due to the complicated mathematical problem. In contrast, sensitivities are largely adopted in the “static” structural optimization. For example, in [47], it is proved that when the response is linear in the uncertain variables, reduction of sensitivity is effective in reducing the probability of malfunctioning and failure; the static deformation of a three-bar truss is assumed as the test case. A sensitivity analysis scheme is used in [30] for improving robustness of the static displacement field in a twenty-five-bar space truss.


### 1.2. Contributions of the Paper

In the light of the limitations of the state of the art in the EA through IDSM, this paper tackles the problem of robustness by proposing a new method to perform robust partial EA (RPEA) for lightly damped linear vibrating systems. Assignment is formulated through the system mass and stiffness matrices and the inverse eigenvalue problem: the assignment of both some natural frequencies, together with the partial definition of the related mode shapes, and some antiresonances can be tackled. Robustness specifications are embedded in the problem through constraints on the sensitivities of the assigned frequencies (natural frequencies of antiresonance frequencies) with respect to the system parameters of interest. Indeed, as can be inferred from the discussion on the literature, the use of sensitivity fits well with the goal of precisely assigning in the nominal conditions and reducing the perturbation due to changes in the system parameters, i.e., of some masses and stiffnesses. Additionally, it does not require the knowledge of the probabilistic distributions of uncertain parameters, which is often difficult to obtain accurately.The system eigenvectors (mode shapes or antiresonant response vectors) could be partially assigned (although this requirement imposes increasing the number of modifications), and the admissible structural modification design variables are constrained to belong to a feasible domain. The goal of the proposed method is to be general, and not targeted to one specific case, leading to a nonlinear, nonseparable, nonconvex constrained least-square problem. Numerical solution is proposed through an algorithm that leverages on homotopy maps to boost the attainment of a global minimum. Variables lifting and McCormick’s constraints are employed to provide a convexification of the bilinear and trilinear terms appearing in the formulation.Due to the adopted optimization-based formulation, the proposed method provides a general mathematical framework suitable for the design of vibrating systems as well as for modification of existing ones. Hence, it is a general-purpose design tool, that can be applied to the robust EA for several systems (e.g., either lumped or finite element models), also with a medium/large number of model coordinates.Four numerical test cases are adopted to assess the effectiveness of the proposed method, and a comparison with a state-of-the-art method for PEA that does not embed robustness is here proposed. Statistical analysis is adopted to assess robustness under random parameter changes.


## 1.1. State of the Art and Motivations

Eigenstructure assignment (EA) has been a core research area in structural dynamics. Several direct and inverse dynamic structural modification (IDSM) methods have been proposed up to now in the literature. In particular, inverse approaches are the most fascinating and challenging for researchers since those require to compute the structural modifications that enable to obtain a prescribed eigenstructure for the modified system [1].EA has been tackled through passive, active, and hybrid approaches (see, e.g., [2–4]). Active approaches exploit the energy introduced in the system by external actuators based on the information of some sensors to perform the EA. The state of the art in active control spans from model-based techniques (see, e.g., [5–9]) to receptance-based methods (see, e.g., [10–13]).Passive approaches are widely investigated since they are stable due to the symmetric and definite positive nature of the structural modifications. Several works have been proposed in the literature to address different problems, such as the assignment of frequencies [14–18], vibration nodes [19], mode shapes and natural frequencies [20, 21], antiresonance frequencies [22–24], or the optimal design of dampers [25]. These approaches are useful both in the design of new systems and in the modification of existing ones, to meet the desired dynamic performances. For example, assigning natural frequencies is useful in resonators [26] to ensure large displacements with small excitation forces, or to place resonances far from excitation sources whenever damping them is not possible or very difficult [27]. Assigning antiresonances, i.e., frequencies where the system experiences negligible steady-state force oscillations at some particular coordinates when excited by harmonic excitations, is useful to reduce unwanted vibrations [22]. Assigning the mode shapes is useful to vibration confinement [3] or to shape the forced response in resonators [26].Both analytical and numerical methods have been proposed. The latter ones are usually adopted to solve more complex and general problems that cannot be handled analytically due to their highly nonlinear nature or to the absence of exact solutions. Hence, approximate solutions should be found through minimization-based formulations. For example, the partial assignment of mode shapes and natural frequencies has been performed by the authors and their co-workers in [20, 26, 28] by solving a constrained nonlinear nonconvex optimization through the homotopy optimization technique.A limitation of most of the proposed techniques to EA is that robustness is not addressed, and hence, it is neglected how the assigned natural frequencies are affected by perturbations of the system parameters about the nominal ones. On the other hand, making the performances minimally sensitive to parameter uncertainty is always of paramount concern in optimal design of machines and structures, leading to the idea of “robust design” [29]. Robust design is widely adopted for the optimal design of the strength of materials and components (such as static problems, buckling, or other structural applications) where the goal is trading between oversizing and ensuring reliability in the range of possible parameter variability in “everyday service conditions” [30], or in the worst-case scenario. In contrast, in the case of EA, operating with worst-case scenario is not reasonable, because it would lead to mistuned eigenstructures in the more likely set of system parameters.In the frame of EA, robustness is often interpreted as finding the optimal eigenstructure, or just the optimal natural frequencies, ensuring adequate performances even under perturbed conditions, and some ad hoc applications have been proposed. For example, in [31, 32], the robust design of brakes with respect to self-excited vibrations has been oriented to split the eigenfrequencies as much as possible by maximizing the difference between consecutive eigenvalues. Robust mechanical design to prevent self-excited vibrations due to frictional forces in screw jacks has been exploited in [33] through the analytical analysis of the eigenvalues of a two degrees of freedom (DOFs) model. Robustness, with respect to uncertainties in the system parameters as well as in the disturbance input force, in the design of tuned mass dampers for vibration absorption is often handled by ensuring a large range of frequencies about the imposed antiresonance where the receptance from the disturbance to the response has a low gain (see, e.g., [34, 35] and the references therein).Robust EA has been sometimes discussed in the case of active approaches that are out of the scope of this paper. Indeed, the control force exerted by the actuators can be changed based on the sensed (or estimated) plant modifications. This possibility, at the cost of a more complicate system (that requires actuators, sensors, and a controller) makes usually active EA more robust than passive EA [13, 36].A different goal, that is the one pursued in this work, is ensuring that the imposed natural frequencies and antiresonance frequencies are robust with respect to parameter changes, e.g., due to changing loads, or to uncertainty on parameters that are not exactly known at the design stage. This goal is of practical interest in structures and machines in the presence of either wanted or unwanted vibrations. For example, ensuring insensitive natural (or resonance) frequencies is of great importance in the design of resonators (such as feeders [26, 37] or sonotrodes [20, 38]), where the resonance frequency should match the excitation frequency regardless of the changes in some parameters, to ensure large motion with small actuation effort. In the case of vibration absorption of harmonic excited systems, assigning antiresonances at the excitation frequency is an effective approach ([39–41]), that should be ensured even in the case of changes of some system parameters. Generally speaking, any application where passive EA is performed can benefit of this approach, since increasing robustness allows to get rid of small parameter changes, compared to the nominal model, of uncertainty in the nominal parameters.An effective nonprobabilistic way to incorporate such a robustness requirement into design is using sensitivity, that is useful whenever probabilistic distributions of uncertain parameters are not available and when the optimal design should be, first of all, well-tuned about a nominal condition (i.e., the more likely system configuration).The use of sensitivity as a measure of robustness in EA has been proved to be effective in the case EA through active control, in particular when higher-rank control allows for more combinations of gains leading to the desired eigenstructure. In [42], the eigen-sensitivities with respect to measurement errors in receptances are adopted to increase robustness in state-feedback robust pole placement, by proposing a receptance-based formulation. In [43], robust pole assignment through state-feedback of friction-induced asymmetric systems is made by assigning both eigenvalues and their eigen-sensitivities with respect to the friction coefficient. In [44], eigen-sensitivity to parameter variations in second-order linear systems under state-feedback and state-derivative feedback is discussed and assessed. In [45], the problem of robust eigenstructure assignment with respect to uncertain parameters in the fitted receptances is solved, by minimizing sensitivity with respect to the fitting parameters. In [46], it is developed an approach to increase the robustness of the closed-loop eigenvalues based on the notion of spectrum sensitivities, for both state-feedback and state-derivative feedback control.In the case of EA through IDSM, a few works consider sensitivities: indeed, incorporating robustness with respect to parameters changes is not straightforward due to the complicated mathematical problem. In contrast, sensitivities are largely adopted in the “static” structural optimization. For example, in [47], it is proved that when the response is linear in the uncertain variables, reduction of sensitivity is effective in reducing the probability of malfunctioning and failure; the static deformation of a three-bar truss is assumed as the test case. A sensitivity analysis scheme is used in [30] for improving robustness of the static displacement field in a twenty-five-bar space truss.


## 1.2. Contributions of the Paper

In the light of the limitations of the state of the art in the EA through IDSM, this paper tackles the problem of robustness by proposing a new method to perform robust partial EA (RPEA) for lightly damped linear vibrating systems. Assignment is formulated through the system mass and stiffness matrices and the inverse eigenvalue problem: the assignment of both some natural frequencies, together with the partial definition of the related mode shapes, and some antiresonances can be tackled. Robustness specifications are embedded in the problem through constraints on the sensitivities of the assigned frequencies (natural frequencies of antiresonance frequencies) with respect to the system parameters of interest. Indeed, as can be inferred from the discussion on the literature, the use of sensitivity fits well with the goal of precisely assigning in the nominal conditions and reducing the perturbation due to changes in the system parameters, i.e., of some masses and stiffnesses. Additionally, it does not require the knowledge of the probabilistic distributions of uncertain parameters, which is often difficult to obtain accurately.The system eigenvectors (mode shapes or antiresonant response vectors) could be partially assigned (although this requirement imposes increasing the number of modifications), and the admissible structural modification design variables are constrained to belong to a feasible domain. The goal of the proposed method is to be general, and not targeted to one specific case, leading to a nonlinear, nonseparable, nonconvex constrained least-square problem. Numerical solution is proposed through an algorithm that leverages on homotopy maps to boost the attainment of a global minimum. Variables lifting and McCormick’s constraints are employed to provide a convexification of the bilinear and trilinear terms appearing in the formulation.Due to the adopted optimization-based formulation, the proposed method provides a general mathematical framework suitable for the design of vibrating systems as well as for modification of existing ones. Hence, it is a general-purpose design tool, that can be applied to the robust EA for several systems (e.g., either lumped or finite element models), also with a medium/large number of model coordinates.Four numerical test cases are adopted to assess the effectiveness of the proposed method, and a comparison with a state-of-the-art method for PEA that does not embed robustness is here proposed. Statistical analysis is adopted to assess robustness under random parameter changes.


## 2. Method Formulation

### 2.1. Definitions

Let us consider aN-degrees of freedom (DOFs) undamped vibrating system with qt denoting the displacement vector (q¨t is the acceleration vector). The system is modelled through its symmetric positive definite mass matrix M = MT > 0 , and its positive semidefinite and symmetric stiffness matrix K = KT ≥ 0. The external forces are denoted by ft. The equation of motion is therefore(1)Mq¨t+Kqt=ft.The model parameters collected inM and K are assumed as the nominal ones, i.e., the most likely ones. For example, in the case of existing systems to be modified, they can be obtained through robust model updating (see, e.g., [48]) to ensure a reliable representation of the true system. In the design of new systems, they can be estimated through finite element models. It should be stressed that the availability of a robust design method is effective in handling the uncertainty on such parameters.Equation (1) can be also transformed in the frequency domain (ω):(2)−ω2M+Kqω=fω,and the receptance matrix Hω=−ω2M+K−1 is introduced, such that qω=Hωfω. The r, c entry of Hω, denoted by hrc(ω), is the transfer function from the external force applied at the c-th coordinate to the displacement of the r-th coordinate, which can be computed through Cramer’s rule of linear system as follows [22, 23, 40]:(3)hrcω=qrωfcω=−1r+cdet−ω2Mc¯r¯+Kc¯r¯det−ω2M+K,where Mc¯r¯ and Kc¯r¯ are the mass and stiffness matrices of the, so-called, adjunct (or adjoint) system and are obtained by removing the c-th row and the r-th column from matrices M and K, respectively. The roots of the numerator in equation (3) are the antiresonances [49]. The zeros of the denominator in equation (3) are the natural frequencies of the system that coincide with its resonance frequencies in the case of undamped systems [49].


### 2.2. Structural Modification for the Assignment of Natural Frequencies and Mode Shapes

Besides solving the characteristic polynomials, it is possible to compute theN natural frequencies ωp,i (i = 1, … , N) together with their related mode shapes up,i∈ℝN (i = 1, …, N) by solving the following eigenproblem [28]:(4)−ωp,i2M+Kup,i=0.Starting from such an equation, the structural modification problem for the assignment of natural frequencies and mode shapes is to find the additive structural modification matricesΔM and ΔK such that [20](5)−ω¯p,i2M+ΔM+K+ΔKu¯p,i=0i=1,...,Np,where Np denotes the number of eigenpairs to be assigned, i.e., the pairs made by the natural frequency and its related mode shape ω¯p,i,u¯p,i. The assumption of additivity for the structural modifications is widely employed in the literature (see, e.g., [18, 49]) and will be exploited also in this paper without lack of generality. Indeed, it represents most of the modifications, such as variations of lumped and distributed mass and stiffness parameters.The solvability of equation (5) is not ensured for all the choices of the structural modification matrices, because of several reasons. First, the rank-matching condition must be satisfied [2], meaning that one unbounded, independent design variable must be introduced to assign one eigenvalue, with no specification on the eigenvector. If the related mode shape should be assigned too, the number of design variables must be increased. Secondly, the presence of constraints on the design variables due to technical and economical limitations on the feasible values drastically limits the attainable eigenpairs and often imposes more design variables to assign just one eigenvalue. Hence, the inverse eigenvalue problem can be conveniently cast as a constrained least-square minimization [20, 49]:(6)minx∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,x∈Γx,where αp,i are weighing coefficients and ⋅2 denotes the Euclidean norm (2-norm) of a vector. Equation (6) also introduces the dependence of the modification matrices on the design variables, collected in vector x=x1,...,xnxT∈ℝNx (nx is their cardinality), to be computed for achieving the desired eigenstructure. The feasible domain Γx is introduced too, which includes upper and lower bounds (denoted by xL and xU, respectively) as well as other convex polyhedra relating different and inhomogeneous variables. It should be noted that formulating the assignment problem as a minimization is also convenient under a numerical point of view, since several algorithms for optimization have been developing in the literature and in commercial software, thus making easier its solution.If mode shapes are not assigned, or just partially assigned for some coordinates, i.e., a PEA problem is to be solved, the unassigned entries of the eigenvector are treated as additional unknowns, that are constrained to belong to the domainΓu¯: u¯pL≤u¯p≤u¯pU. Such a constrain can be used to define set interval specifications on the desired mode shapes, as well as to allow for the use of McCormick’s constraints for an effective numerical solution (see Section 3), even if no specifications on the mode shapes are desired. The least-square problem in equation (6) is further recast into the following form:(7)minx,u¯p∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,xu¯p∈Γ.A new feasible domain is defined as the union of the constraints onx and u¯p, Γ=Γx∪Γu¯.In writing equation (7) and in defining Γu¯, it has been implicitly assumed that all the entries of u¯p,i are unassigned, to obtain a simpler formulation; the extension to the partial assignment of just some entries of u¯p,i is obvious under a mathematical point of view. However, imposing some mode shape entries usually makes achieving the desired natural frequencies more challenging, if a few design variables are assumed and constraints on the physical modifications are tight.An effective approach to solve the nonlinear and nonconvex equation (7) is adopted in Section 3, by exploiting the methods proposed by the authors and their co-workers, in [28, 48].


### 2.3. Inclusion of Robustness Constraints

The problem introduced in equation (7) can be improved by including robustness specifications, by exploiting the concept of sensitivity of natural frequencies with respect to some critical parameters that are uncertain or might vary. Introducing sensitivities is useful to handle, for example, the presence of uncertainty on the parameters of the nominal model in equation (1), whenever robust model updating is not performed if it does not deliver reliable results. On the other hand, some parameters might change during operations (e.g., due to the presence of loads).Let us collect thenp parameters affecting M and K into vector θ=θ1,...,θnpT. The sensitivity the of i-th natural frequency ωp,i with respect to the j-th parameter θj, denoted by S¯p,ij=∂ω¯p,i2/∂θj for brevity, is computed as follows [50]:(8)∂ωp,i2∂θj=up,iT−ωp,i2∂M/∂θj+∂K/∂θjup,iup,iTMup,i.Proof.
Let us consider the natural frequency eigenproblem in equation (4) and let us assume that M and K depend on an uncertain parameter θj, i.e., Mθj and Kθj. Developing the partial derivatives leads to the following relation:(9)−ωp,i2∂M∂θj+∂K∂θj−∂ωp,i2∂θjMup,i+−ωp,i2M+K∂up,i∂θj=0.
By premultiplying equation (9) by up,iT, it yields(10)up,iT−ωp,i2∂M∂θj+∂K∂θj−∂ωp,i2∂θjMup,i+up,iT−ωp,i2M+K∂up,i∂θj=0.
SinceM and K are symmetric, then up,iT−ωp,i2M+K=0. Hence, the second term on the left-hand side of equation (10) is null, leading to equation (8) after some simple algebraic manipulations.Equation (8) reveals that sensitivities depend on the Jacobian matrices JMθj=∂Mθ/∂θj and JKθj=∂Kθ/∂θj (henceforth, J is adopted to represent the Jacobian of the matrix written as the subscript with respect to the parameter written as the superscript), and on the mode shapes up,i. It should be noted that arbitrary normalization of the mode shapes is assumed in equation (8) to allow for a more general representation of sensitivity within the frame of the IDSM problem.The goal of the RPEA problem solved in this section is to assign the prescribedNp eigenpairs ω¯p,i,u¯p,i together with the sensitivities of the natural frequencies with respect to some critical parameters (whose cardinality is denoted by Ns). Such sensitivities S¯p,ij, for the indexes i and j of interest, are required to assume the prescribed values making the natural frequency variations admissible as the parameter θj changes within an expected range of variation. Indeed, S¯p,ij represents the local, first-order approximation of the effect of each parameter on the natural frequencies. The smaller the S¯p,ij is, the more robust the assigned ω¯p,i is with respect to changes of θj.In the case of the system modified throughΔM and ΔK, the sensitivity S¯p,ij is inferred from equation (8) by taking into account the structural modification matrices:(11)S¯p,ij=u¯p,iT−ω¯p,i2JM+ΔMθj+JK+ΔKθju¯p,iu¯p,iTM+ΔMu¯p,i.The inclusion of requirements onS¯p,ij within the IDSM problem is not trivial, due to its fractional form, which would lead to a fractional function of the problem unknown x, hence exacerbating its nonlinear, nonseparable, nonconvex nature. This issue can be overcome by translating it as a set of constraints, having the form(12)Sp,1x,u¯pT⋯Sp,Ndx,u¯pT︸Spx,u¯p=0.Sp,ix,u¯p∈ℝNs collects the specification on the sensitivities of ω¯p,i2 with respect to the Ns parameters, written in nonfractional form:(13)Sp,ix,u¯p=S¯p,i1x,u¯pDix,u¯p−Np,i1u¯p⋮S¯p,iNSx,u¯pDix,u¯p−Np,iNSu¯p,where(14)Niju¯p=u¯p,iT−ω¯p,i2JMθj+JKθju¯p,i,Dix,u¯p=u¯p,iTMθ+ΔMθ,xu¯p,i.By including it within the inverse eigenvalue problem, the following problem is obtained:(15)minx,u¯p∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,xu¯p∈Γ,subjSpx,u¯p=0,where the prescribed sensitivities are introduced as an additional set of nonlinear and nonconvex constraints. To tackle the nonconvexity of the constraint, an effective solution approach is proposed in Section 3.


### 2.4. Robust Assignment of Antiresonance Frequencies

The method proposed in Sections2.2 and 2.3 can be extended to perform robust assignment of antiresonances, by considering that antiresonance frequencies ωz,i together with their related right eigenvectors uz,i∈ℝN−1 can be computed by solving the eigenproblem of the adjunct system:(16)−ωz,i2Mc¯r¯+Kc¯r¯uz,i=0.Hence, as proposed by the authors in [22], antiresonance assignment can be cast as inverse eigenvalue problem of the adjunct system, which recalls the one in equation (7):(17)minx,u¯z∑i=1Nzαz,i−ω¯z,i2Mc¯r¯+ΔMc¯r¯x+Kc¯r¯+ΔKc¯r¯xu¯z,i22,xu¯z∈Γ,where ω¯z,i is the i-th prescribed antiresonance frequency and αz,i are the weighing coefficients. As in equation (7), both the parameter modifications and the right eigenvectors are constrained through the convex domain Γ. The box constraint on u¯z could be used, in principle, to define set interval specifications on the antiresonant response [41, 51]. However, imposing some entries of u¯z usually makes the achievement of the desired antiresonance very challenging. As a matter of fact, just one paper in the literature assigns u¯z through DSM [52]. In contrast, it is always neglected (see, e.g., [22, 53]). The difficulties in assigning u¯z are exacerbated if it is wanted to assign antiresonance sensitivity, unless a lot of modifiable parameters are assumed, and large constraints are set. Therefore, in practice, box constraint on u¯z is just set to allow for the use of McCormick’s constraints for an effective numerical solution (see Section 3), even if no specifications on the antiresonant response are wanted.Antiresonance sensitivities,Sz,ij=∂ωz,i2/∂θj, are computed as follows [51]:(18)Sz,ij=lz,iT−ωz,i2∂Mc¯r¯/∂θj+∂Kc¯r¯/∂θjuz,iIz,iTMc¯r¯uz,i,where lz,i and uz,i are, respectively, the left and right eigenvectors of the adjunct system, which is asymmetric in the case of cross receptances, i.e., if r ≠ c. By considering the eigenproblem of the adjunct system in equation (16) and that lz,iT−ωz,i2Mc¯r¯+Kc¯r¯=0, the proof of equation (18) is similar to the one of equation (8) and is here omitted for brevity.By considering the system with the mass and stiffness modifications, the sensitivityS¯z,ij=∂ω¯z,i2/∂θj is defined as follows:(19)S¯z,ij=l¯z,iT−ω¯z,i2JMc¯r¯+ΔMc¯r¯θj+JKc¯r¯+ΔKc¯r¯θju¯z,il¯z,iTMc¯r¯+ΔMc¯r¯u¯z,i.By interpreting antiresonance assignment as a RPEA problem in the adjunct system (as proposed in [22]), the extension of the method proposed in equations from (12) to (15) for the robust assignment of the antiresonance frequencies is obvious and here omitted for brevity and leads to a problem similar to equation (15).


## 2.1. Definitions

Let us consider aN-degrees of freedom (DOFs) undamped vibrating system with qt denoting the displacement vector (q¨t is the acceleration vector). The system is modelled through its symmetric positive definite mass matrix M = MT > 0 , and its positive semidefinite and symmetric stiffness matrix K = KT ≥ 0. The external forces are denoted by ft. The equation of motion is therefore(1)Mq¨t+Kqt=ft.The model parameters collected inM and K are assumed as the nominal ones, i.e., the most likely ones. For example, in the case of existing systems to be modified, they can be obtained through robust model updating (see, e.g., [48]) to ensure a reliable representation of the true system. In the design of new systems, they can be estimated through finite element models. It should be stressed that the availability of a robust design method is effective in handling the uncertainty on such parameters.Equation (1) can be also transformed in the frequency domain (ω):(2)−ω2M+Kqω=fω,and the receptance matrix Hω=−ω2M+K−1 is introduced, such that qω=Hωfω. The r, c entry of Hω, denoted by hrc(ω), is the transfer function from the external force applied at the c-th coordinate to the displacement of the r-th coordinate, which can be computed through Cramer’s rule of linear system as follows [22, 23, 40]:(3)hrcω=qrωfcω=−1r+cdet−ω2Mc¯r¯+Kc¯r¯det−ω2M+K,where Mc¯r¯ and Kc¯r¯ are the mass and stiffness matrices of the, so-called, adjunct (or adjoint) system and are obtained by removing the c-th row and the r-th column from matrices M and K, respectively. The roots of the numerator in equation (3) are the antiresonances [49]. The zeros of the denominator in equation (3) are the natural frequencies of the system that coincide with its resonance frequencies in the case of undamped systems [49].


## 2.2. Structural Modification for the Assignment of Natural Frequencies and Mode Shapes

Besides solving the characteristic polynomials, it is possible to compute theN natural frequencies ωp,i (i = 1, … , N) together with their related mode shapes up,i∈ℝN (i = 1, …, N) by solving the following eigenproblem [28]:(4)−ωp,i2M+Kup,i=0.Starting from such an equation, the structural modification problem for the assignment of natural frequencies and mode shapes is to find the additive structural modification matricesΔM and ΔK such that [20](5)−ω¯p,i2M+ΔM+K+ΔKu¯p,i=0i=1,...,Np,where Np denotes the number of eigenpairs to be assigned, i.e., the pairs made by the natural frequency and its related mode shape ω¯p,i,u¯p,i. The assumption of additivity for the structural modifications is widely employed in the literature (see, e.g., [18, 49]) and will be exploited also in this paper without lack of generality. Indeed, it represents most of the modifications, such as variations of lumped and distributed mass and stiffness parameters.The solvability of equation (5) is not ensured for all the choices of the structural modification matrices, because of several reasons. First, the rank-matching condition must be satisfied [2], meaning that one unbounded, independent design variable must be introduced to assign one eigenvalue, with no specification on the eigenvector. If the related mode shape should be assigned too, the number of design variables must be increased. Secondly, the presence of constraints on the design variables due to technical and economical limitations on the feasible values drastically limits the attainable eigenpairs and often imposes more design variables to assign just one eigenvalue. Hence, the inverse eigenvalue problem can be conveniently cast as a constrained least-square minimization [20, 49]:(6)minx∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,x∈Γx,where αp,i are weighing coefficients and ⋅2 denotes the Euclidean norm (2-norm) of a vector. Equation (6) also introduces the dependence of the modification matrices on the design variables, collected in vector x=x1,...,xnxT∈ℝNx (nx is their cardinality), to be computed for achieving the desired eigenstructure. The feasible domain Γx is introduced too, which includes upper and lower bounds (denoted by xL and xU, respectively) as well as other convex polyhedra relating different and inhomogeneous variables. It should be noted that formulating the assignment problem as a minimization is also convenient under a numerical point of view, since several algorithms for optimization have been developing in the literature and in commercial software, thus making easier its solution.If mode shapes are not assigned, or just partially assigned for some coordinates, i.e., a PEA problem is to be solved, the unassigned entries of the eigenvector are treated as additional unknowns, that are constrained to belong to the domainΓu¯: u¯pL≤u¯p≤u¯pU. Such a constrain can be used to define set interval specifications on the desired mode shapes, as well as to allow for the use of McCormick’s constraints for an effective numerical solution (see Section 3), even if no specifications on the mode shapes are desired. The least-square problem in equation (6) is further recast into the following form:(7)minx,u¯p∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,xu¯p∈Γ.A new feasible domain is defined as the union of the constraints onx and u¯p, Γ=Γx∪Γu¯.In writing equation (7) and in defining Γu¯, it has been implicitly assumed that all the entries of u¯p,i are unassigned, to obtain a simpler formulation; the extension to the partial assignment of just some entries of u¯p,i is obvious under a mathematical point of view. However, imposing some mode shape entries usually makes achieving the desired natural frequencies more challenging, if a few design variables are assumed and constraints on the physical modifications are tight.An effective approach to solve the nonlinear and nonconvex equation (7) is adopted in Section 3, by exploiting the methods proposed by the authors and their co-workers, in [28, 48].


## 2.3. Inclusion of Robustness Constraints

The problem introduced in equation (7) can be improved by including robustness specifications, by exploiting the concept of sensitivity of natural frequencies with respect to some critical parameters that are uncertain or might vary. Introducing sensitivities is useful to handle, for example, the presence of uncertainty on the parameters of the nominal model in equation (1), whenever robust model updating is not performed if it does not deliver reliable results. On the other hand, some parameters might change during operations (e.g., due to the presence of loads).Let us collect thenp parameters affecting M and K into vector θ=θ1,...,θnpT. The sensitivity the of i-th natural frequency ωp,i with respect to the j-th parameter θj, denoted by S¯p,ij=∂ω¯p,i2/∂θj for brevity, is computed as follows [50]:(8)∂ωp,i2∂θj=up,iT−ωp,i2∂M/∂θj+∂K/∂θjup,iup,iTMup,i.Proof.
Let us consider the natural frequency eigenproblem in equation (4) and let us assume that M and K depend on an uncertain parameter θj, i.e., Mθj and Kθj. Developing the partial derivatives leads to the following relation:(9)−ωp,i2∂M∂θj+∂K∂θj−∂ωp,i2∂θjMup,i+−ωp,i2M+K∂up,i∂θj=0.
By premultiplying equation (9) by up,iT, it yields(10)up,iT−ωp,i2∂M∂θj+∂K∂θj−∂ωp,i2∂θjMup,i+up,iT−ωp,i2M+K∂up,i∂θj=0.
SinceM and K are symmetric, then up,iT−ωp,i2M+K=0. Hence, the second term on the left-hand side of equation (10) is null, leading to equation (8) after some simple algebraic manipulations.Equation (8) reveals that sensitivities depend on the Jacobian matrices JMθj=∂Mθ/∂θj and JKθj=∂Kθ/∂θj (henceforth, J is adopted to represent the Jacobian of the matrix written as the subscript with respect to the parameter written as the superscript), and on the mode shapes up,i. It should be noted that arbitrary normalization of the mode shapes is assumed in equation (8) to allow for a more general representation of sensitivity within the frame of the IDSM problem.The goal of the RPEA problem solved in this section is to assign the prescribedNp eigenpairs ω¯p,i,u¯p,i together with the sensitivities of the natural frequencies with respect to some critical parameters (whose cardinality is denoted by Ns). Such sensitivities S¯p,ij, for the indexes i and j of interest, are required to assume the prescribed values making the natural frequency variations admissible as the parameter θj changes within an expected range of variation. Indeed, S¯p,ij represents the local, first-order approximation of the effect of each parameter on the natural frequencies. The smaller the S¯p,ij is, the more robust the assigned ω¯p,i is with respect to changes of θj.In the case of the system modified throughΔM and ΔK, the sensitivity S¯p,ij is inferred from equation (8) by taking into account the structural modification matrices:(11)S¯p,ij=u¯p,iT−ω¯p,i2JM+ΔMθj+JK+ΔKθju¯p,iu¯p,iTM+ΔMu¯p,i.The inclusion of requirements onS¯p,ij within the IDSM problem is not trivial, due to its fractional form, which would lead to a fractional function of the problem unknown x, hence exacerbating its nonlinear, nonseparable, nonconvex nature. This issue can be overcome by translating it as a set of constraints, having the form(12)Sp,1x,u¯pT⋯Sp,Ndx,u¯pT︸Spx,u¯p=0.Sp,ix,u¯p∈ℝNs collects the specification on the sensitivities of ω¯p,i2 with respect to the Ns parameters, written in nonfractional form:(13)Sp,ix,u¯p=S¯p,i1x,u¯pDix,u¯p−Np,i1u¯p⋮S¯p,iNSx,u¯pDix,u¯p−Np,iNSu¯p,where(14)Niju¯p=u¯p,iT−ω¯p,i2JMθj+JKθju¯p,i,Dix,u¯p=u¯p,iTMθ+ΔMθ,xu¯p,i.By including it within the inverse eigenvalue problem, the following problem is obtained:(15)minx,u¯p∑i=1Npαp,i−ω¯p,i2M+ΔMx+K+ΔKxu¯p,i22,xu¯p∈Γ,subjSpx,u¯p=0,where the prescribed sensitivities are introduced as an additional set of nonlinear and nonconvex constraints. To tackle the nonconvexity of the constraint, an effective solution approach is proposed in Section 3.


## 2.4. Robust Assignment of Antiresonance Frequencies

The method proposed in Sections2.2 and 2.3 can be extended to perform robust assignment of antiresonances, by considering that antiresonance frequencies ωz,i together with their related right eigenvectors uz,i∈ℝN−1 can be computed by solving the eigenproblem of the adjunct system:(16)−ωz,i2Mc¯r¯+Kc¯r¯uz,i=0.Hence, as proposed by the authors in [22], antiresonance assignment can be cast as inverse eigenvalue problem of the adjunct system, which recalls the one in equation (7):(17)minx,u¯z∑i=1Nzαz,i−ω¯z,i2Mc¯r¯+ΔMc¯r¯x+Kc¯r¯+ΔKc¯r¯xu¯z,i22,xu¯z∈Γ,where ω¯z,i is the i-th prescribed antiresonance frequency and αz,i are the weighing coefficients. As in equation (7), both the parameter modifications and the right eigenvectors are constrained through the convex domain Γ. The box constraint on u¯z could be used, in principle, to define set interval specifications on the antiresonant response [41, 51]. However, imposing some entries of u¯z usually makes the achievement of the desired antiresonance very challenging. As a matter of fact, just one paper in the literature assigns u¯z through DSM [52]. In contrast, it is always neglected (see, e.g., [22, 53]). The difficulties in assigning u¯z are exacerbated if it is wanted to assign antiresonance sensitivity, unless a lot of modifiable parameters are assumed, and large constraints are set. Therefore, in practice, box constraint on u¯z is just set to allow for the use of McCormick’s constraints for an effective numerical solution (see Section 3), even if no specifications on the antiresonant response are wanted.Antiresonance sensitivities,Sz,ij=∂ωz,i2/∂θj, are computed as follows [51]:(18)Sz,ij=lz,iT−ωz,i2∂Mc¯r¯/∂θj+∂Kc¯r¯/∂θjuz,iIz,iTMc¯r¯uz,i,where lz,i and uz,i are, respectively, the left and right eigenvectors of the adjunct system, which is asymmetric in the case of cross receptances, i.e., if r ≠ c. By considering the eigenproblem of the adjunct system in equation (16) and that lz,iT−ωz,i2Mc¯r¯+Kc¯r¯=0, the proof of equation (18) is similar to the one of equation (8) and is here omitted for brevity.By considering the system with the mass and stiffness modifications, the sensitivityS¯z,ij=∂ω¯z,i2/∂θj is defined as follows:(19)S¯z,ij=l¯z,iT−ω¯z,i2JMc¯r¯+ΔMc¯r¯θj+JKc¯r¯+ΔKc¯r¯θju¯z,il¯z,iTMc¯r¯+ΔMc¯r¯u¯z,i.By interpreting antiresonance assignment as a RPEA problem in the adjunct system (as proposed in [22]), the extension of the method proposed in equations from (12) to (15) for the robust assignment of the antiresonance frequencies is obvious and here omitted for brevity and leads to a problem similar to equation (15).


## 3. Problem Solution

### 3.1. Outline of the Proposed Method

The proposed method relies on the idea of homotopy optimization that has been adopted by the authors and their co-workers, in solving PEA problems similar to equation (7) [28, 52], antiresonance assignment [22, 52], as well as to the response optimization of underactuated multibody systems [37]. In this paper, the technique is extended to solve equation (15), i.e., the PEA problem embedding a robustness condition. The problem is not straightforward due to the robustness constraints, which imposes an ad hoc approach to handle it.Homotopy is an optimization technique that effectively solves nonconvex optimizations. It has been adopted for global optimization since it has been widely proved to be effective in boosting the convergence to the global optimal solution [54, 55]. The idea is to approximate the original nonconvex minimization problem (including both the function and the constraints) through a convex minimization one. Then, by means of a finite number of optimization step, the minimization is morphed back to the original nonconvex problem through the so-called homotopy map, by leading to a path of solutions.In this work, three basic steps are exploited to solve the RPEA problems in equation (15) and the related extension to antiresonance assignment (see Section 2.4):(1)
Variable lifting (Section3.2)(2)
Definition of the homotopy map (Section3.3)(3)
Introduction of bilinear and trilinear McCormick’s constraints (Section3.4)


### 3.2. Variables Lifting

Variable lifting is the introduction of some additional variables that make the problem convex with respect to such variables. It is often adopted in optimization of nonconvex problems to transform them into convex ones with respect of other additional variables [22, 56]. Nonconvexity in equations (12), (13), and (15) arises for three reasons.First, the cost function has some bilinear terms of the formxwu¯g (for proper indexes w and g) and variables u¯g may include entries of u¯p,i as well as of u¯z,i. Bilinear terms are replaced in the convex approximation of the cost function by auxiliary variables denoted by bl, i.e.,(20)bl⟶xwu¯g.Secondly, the sensitivity constraints introduce bilinear terms of the formu¯ou¯g (for proper indexes o and g) that are replaced in the convex approximation of Sx,u¯ by the auxiliary variables vm:(21)vm⟶u¯ou¯g.Lastly, trilinear terms also appear in the sensitivity constraint asu¯oxwu¯g (for proper indexes o, w, and g), and without loss of generality, u¯o may include entries of u¯p,i as well as entries of l¯z,i. These terms are replaced in the convex approximation of S (both Spx,u¯p and Szx,u¯z) by the auxiliary variables tk:(22)tk⟶u¯oxwu¯g.The proposed variable lifting approach leads to the convex approximations of the nonconvex cost function (denoted byfnc and convexified through fc) and of the sensitivity constraint S (convexified through Sc) in equation (15), by lifting the optimization problem to a higher-dimensional one, which is linear in the auxiliary variables bl, vm, and tk (for proper indexes l, m, and k).


### 3.3. Homotopy Maps

The cost function and the constraints are morphed back to the original nonconvex one by means of the following homotopy maps:(23)fh=1−λfc+λfnc,Sh=1−λSc+λS,where λ is a parameter that discretely spans from 0 to 1. In practice, a sequence of Nλ optimizations is performed: in the first one (i.e., λ = 0), the convex problem defined by fc and Sc is solved; in the last one (i.e., λ = 1), the original nonconvex problem defined by fnc and S is solved. In the intermediate steps, the problems defined by fh(λ) and Sh(λ) are solved.The initial guess of the nonconvex minimization offh in the i-th iteration is the optimal solution of the previous step, i.e., x0,i=xopt,i−1. Finally, the values of the optimal design variables are computed as ΔMxλ=1 and ΔKxλ=1. In this way, a path of solutions converging to the optimal one is obtained.In this work,fh and Sh are obtained by replacing each bilinear or trilinear terms with an affine combination of themselves and their related auxiliary variable, i.e.,(24)xwu¯g⟶1−λbl+λxwu¯gu¯ou¯g⟶1−λvm+λu¯ou¯gu¯oxwu¯z⟶1−λtk+λu¯oxwu¯g.A flowchart of the method is proposed in Figure1.Figure 1
Flowchart of the solution method.


### 3.4. Bilinear and Trilinear McCormick’s Constraints

Additional convex constraints should be defined to make the auxiliary variables correctly approximate the nonlinear terms. The tightest convexification of the constraints for bilinear terms, i.e., the convex hull, is obtained through bilinear McCormick’s constraints (see, e.g., [57]):(25)bl≥xwLu¯g+u¯gLxw−xwLu¯gLbl≥xwUu¯g+u¯gUxw−xwUu¯gUbl≤xwUu¯g+u¯gLxw−xwUu¯gLbl≤xwLu¯g+u¯gUxw−xwLu¯gU,where the upper and lower bounds for xw are denoted by xwU and xwL, respectively, while the for u¯g, the upper and lower bounds are u¯gU and u¯gL. For brevity, the trivial definition of bilinear McCormick’s constraints for the auxiliary variable vm that lifts u¯ou¯g is omitted.In this paper, a convexification of the constraints for the trilinear terms is also provided based on the idea proposed in [58]. The convexification of the constraints for trilinear terms is obtained by performing a double-McCormick convexification: first, an additional auxiliary variable nk, such that tk=u¯oxwu¯g=nkxw, is introduced; then, the following inequalities are defined:(26)nk≥u¯oLu¯g+u¯gLu¯o−u¯oLu¯gLnk≥u¯oUu¯g+u¯gUu¯o−u¯oUu¯gUnk≤u¯oUu¯g+u¯gLu¯o−u¯oUu¯gLnk≤u¯oLu¯g+u¯gUu¯o−u¯oLu¯gU,(27)tk≥xwLnk+u¯oLu¯gLxw−xwLu¯oLu¯gLtk≥xwUnk+u¯oUu¯gUxw−xwUu¯oUu¯gUtk≤xwUnk+u¯oLu¯gLxw−xwUu¯oLu¯gLtk≤xwLnk+u¯oUu¯gUxw−xwLu¯oUu¯gU.Inequalities in (26) and (27) can be clustered together through the Fourier–Motzkin elimination, which is a procedure that enables to eliminate the auxiliary variable nk, leading to 14 nonredundant linear constraints. The formulation of the inequalities obtained after the Fourier–Motzkin elimination is here omitted for brevity; for further details, refer to [58].


## 3.1. Outline of the Proposed Method

The proposed method relies on the idea of homotopy optimization that has been adopted by the authors and their co-workers, in solving PEA problems similar to equation (7) [28, 52], antiresonance assignment [22, 52], as well as to the response optimization of underactuated multibody systems [37]. In this paper, the technique is extended to solve equation (15), i.e., the PEA problem embedding a robustness condition. The problem is not straightforward due to the robustness constraints, which imposes an ad hoc approach to handle it.Homotopy is an optimization technique that effectively solves nonconvex optimizations. It has been adopted for global optimization since it has been widely proved to be effective in boosting the convergence to the global optimal solution [54, 55]. The idea is to approximate the original nonconvex minimization problem (including both the function and the constraints) through a convex minimization one. Then, by means of a finite number of optimization step, the minimization is morphed back to the original nonconvex problem through the so-called homotopy map, by leading to a path of solutions.In this work, three basic steps are exploited to solve the RPEA problems in equation (15) and the related extension to antiresonance assignment (see Section 2.4):(1)
Variable lifting (Section3.2)(2)
Definition of the homotopy map (Section3.3)(3)
Introduction of bilinear and trilinear McCormick’s constraints (Section3.4)


## 3.2. Variables Lifting

Variable lifting is the introduction of some additional variables that make the problem convex with respect to such variables. It is often adopted in optimization of nonconvex problems to transform them into convex ones with respect of other additional variables [22, 56]. Nonconvexity in equations (12), (13), and (15) arises for three reasons.First, the cost function has some bilinear terms of the formxwu¯g (for proper indexes w and g) and variables u¯g may include entries of u¯p,i as well as of u¯z,i. Bilinear terms are replaced in the convex approximation of the cost function by auxiliary variables denoted by bl, i.e.,(20)bl⟶xwu¯g.Secondly, the sensitivity constraints introduce bilinear terms of the formu¯ou¯g (for proper indexes o and g) that are replaced in the convex approximation of Sx,u¯ by the auxiliary variables vm:(21)vm⟶u¯ou¯g.Lastly, trilinear terms also appear in the sensitivity constraint asu¯oxwu¯g (for proper indexes o, w, and g), and without loss of generality, u¯o may include entries of u¯p,i as well as entries of l¯z,i. These terms are replaced in the convex approximation of S (both Spx,u¯p and Szx,u¯z) by the auxiliary variables tk:(22)tk⟶u¯oxwu¯g.The proposed variable lifting approach leads to the convex approximations of the nonconvex cost function (denoted byfnc and convexified through fc) and of the sensitivity constraint S (convexified through Sc) in equation (15), by lifting the optimization problem to a higher-dimensional one, which is linear in the auxiliary variables bl, vm, and tk (for proper indexes l, m, and k).


## 3.3. Homotopy Maps

The cost function and the constraints are morphed back to the original nonconvex one by means of the following homotopy maps:(23)fh=1−λfc+λfnc,Sh=1−λSc+λS,where λ is a parameter that discretely spans from 0 to 1. In practice, a sequence of Nλ optimizations is performed: in the first one (i.e., λ = 0), the convex problem defined by fc and Sc is solved; in the last one (i.e., λ = 1), the original nonconvex problem defined by fnc and S is solved. In the intermediate steps, the problems defined by fh(λ) and Sh(λ) are solved.The initial guess of the nonconvex minimization offh in the i-th iteration is the optimal solution of the previous step, i.e., x0,i=xopt,i−1. Finally, the values of the optimal design variables are computed as ΔMxλ=1 and ΔKxλ=1. In this way, a path of solutions converging to the optimal one is obtained.In this work,fh and Sh are obtained by replacing each bilinear or trilinear terms with an affine combination of themselves and their related auxiliary variable, i.e.,(24)xwu¯g⟶1−λbl+λxwu¯gu¯ou¯g⟶1−λvm+λu¯ou¯gu¯oxwu¯z⟶1−λtk+λu¯oxwu¯g.A flowchart of the method is proposed in Figure1.Figure 1
Flowchart of the solution method.


## 3.4. Bilinear and Trilinear McCormick’s Constraints

Additional convex constraints should be defined to make the auxiliary variables correctly approximate the nonlinear terms. The tightest convexification of the constraints for bilinear terms, i.e., the convex hull, is obtained through bilinear McCormick’s constraints (see, e.g., [57]):(25)bl≥xwLu¯g+u¯gLxw−xwLu¯gLbl≥xwUu¯g+u¯gUxw−xwUu¯gUbl≤xwUu¯g+u¯gLxw−xwUu¯gLbl≤xwLu¯g+u¯gUxw−xwLu¯gU,where the upper and lower bounds for xw are denoted by xwU and xwL, respectively, while the for u¯g, the upper and lower bounds are u¯gU and u¯gL. For brevity, the trivial definition of bilinear McCormick’s constraints for the auxiliary variable vm that lifts u¯ou¯g is omitted.In this paper, a convexification of the constraints for the trilinear terms is also provided based on the idea proposed in [58]. The convexification of the constraints for trilinear terms is obtained by performing a double-McCormick convexification: first, an additional auxiliary variable nk, such that tk=u¯oxwu¯g=nkxw, is introduced; then, the following inequalities are defined:(26)nk≥u¯oLu¯g+u¯gLu¯o−u¯oLu¯gLnk≥u¯oUu¯g+u¯gUu¯o−u¯oUu¯gUnk≤u¯oUu¯g+u¯gLu¯o−u¯oUu¯gLnk≤u¯oLu¯g+u¯gUu¯o−u¯oLu¯gU,(27)tk≥xwLnk+u¯oLu¯gLxw−xwLu¯oLu¯gLtk≥xwUnk+u¯oUu¯gUxw−xwUu¯oUu¯gUtk≤xwUnk+u¯oLu¯gLxw−xwUu¯oLu¯gLtk≤xwLnk+u¯oUu¯gUxw−xwLu¯oUu¯gU.Inequalities in (26) and (27) can be clustered together through the Fourier–Motzkin elimination, which is a procedure that enables to eliminate the auxiliary variable nk, leading to 14 nonredundant linear constraints. The formulation of the inequalities obtained after the Fourier–Motzkin elimination is here omitted for brevity; for further details, refer to [58].


## 4. Numerical Examples

### 4.1. Implementation Issues

Four numerical test cases are proposed in this paper. The assignment problems are solved on a 64-bit laptop equipped with a 16 GB RAM and a i7 9750-H, 2.60 GHz, 6 core processor. No parallel computing is exploited, since optimizing the CPU time is not the goal of the proposed work; hence, just one core is adopted in practice.To provide an estimation of the computational effort, let us consider the general scenario involving the assignment ofNp natural frequencies and Nz antiresonances, and let us assume that Ns sensitivities are prescribed for each assigned frequency. Additionally, let us consider the more widespread case where the eigenvectors are let unassigned. Since eigenvectors are defined regardless of their normalization, one component can be arbitrarily assigned (e.g., imposed to 1), as discussed in [22, 41].Bilinear and trilinear terms arise both in the cost function and in the sensitivity constraint. Let us first consider the bilinear termsb due to the inverse eigenvalue problem to perform natural frequency and/or antiresonance assignment. In accordance with equations (7) and (17), the exact computation of the number of the bilinear variables b (that are related to products between a design variable and an eigenvector entry) depends on the structure of matrices ΔMx and ΔKx. Its worst-case upper bound can be estimated by assuming that each design variable in x belongs to all the entries of the modification matrices, and it follows that the upper bound on the cardinality of the b variables is NxNpN−1+NzN−2. Usually, this is an overestimate of the actual cardinality, indeed each design variable just affects few entries of ΔMx and ΔKx.Similarly, by analysing the sensitivity constraint in equations (11) and (19), an upper limit of the number of the bilinear variables v (i.e., those related to products between two eigenvector entries) can be provided since it depends on the structure of ΔMx and ΔKx and on the design variables considered for the sensitivity. Hence, their maximum number is NsNpN−12+NzN−22.Finally, the maximum number of trilinear termst, that depends on the structure of the modification matrices as well, is NxNsNpN−12+NzN−22. Again, this evaluation just provides a worst-case estimation and in practice the number of those terms is significantly smaller, as it can be inferred from Table 1.Table 1
Summary of the computational details of the test cases.
Test-case 1Test-case 2Test-case 3Test-case 4CPU runtime (s)52.120.7137.0180.4Number of DOFs:N551326Number of assigned natural frequencies:Np1010Number of assigned antiresonances:Nz0101Number of assigned sensitivities:Ns3231Number of design variablesNx106825Number of unassigned right eigenvector entriesu431224Number of unassigned right eigenvector entriesl-3-24Number of bilinear variablesb, v12656181Number of trilinear variablest4130135Conversely, the size of the cost function in equation (23) is exactly determined by the number of DOFs of the system and by the number of assigned natural frequencies and/or antiresonances. In the general case, it follows that the number of fh is equal to NpN+NzN−1. Similarly, the number of sensitivity constraints Sh in equation (23) is equal to NsNp+Nz.The computational aspects of each test case proposed in this work are summarized in Table1. The CPU runtime grows if the number of DOFs and the assignment specifications increases. However, the CPU runtime is reasonably small, since IDSM is an offline design technique and does not require real-time computing.The algorithm is implemented in MATLAB through the Yalmip toolbox [59]. Bilinear and trilinear terms are automatically detected, and the variable lifting and McCormick’s convexification are handled through a MATLAB automatic routine. The convex optimization in the first step of the homotopy is solved through the solver Mosek, while in the following iterations, fmincon is adopted. Nonetheless, any numerical solver for linear and nonlinear optimization can be adopted depending on the homotopy stage considered. The Jacobian matrices in equations (11) and (19) have been computed in a fast and accurate way through the complex-step derivatives [48, 60].


### 4.2. Test-Case 1: Robust Assignment of a Natural Frequency in a Five-Mass System

The first test case is a five-mass system widely adopted in the literature as a benchmark for IDSM methods (see, e.g., [16, 23, 61]). The system is sketched in Figure 2, and its parameters are summarized in Table 2.Figure 2
Test-case 1: sketch of the five-mass system.Table 2
Test-case 1: original and modified system design variables and constraints.
ParameterValueModification boundsModified system without robustness (method [28])Modified system with robustness (proposed method)m1 (kg)1.727[0; 1.080]2.3132.751m2 (kg)5.123[0; 1.080]5.7626.202m3 (kg)8.214[0; 1.080]8.8138.214m4 (kg)2.609[0; 1.080]3.1573.671m5 (kg)1.339[0; 1.080]1.8802.417kg1 (Nm−1)94.26e3[0; 144.9e3]140.82e3129.52e3kg2 (Nm−1)94.26e3[0; 144.9e3]106.74e3231.74e3kg3 (Nm−1)94.26e3[0; 144.9e3]173.15e3135.58e3kg4 (Nm−1)94.26e3[0; 144.9e3]126.78e3101.51e3kg5 (Nm−1)94.26e3[0; 144.9e3]115.95e397.429e3k12 (Nm−1)75.14e3—75.14e375.14e3k23 (Nm−1)67.74e3—67.74e367.74e3k34 (Nm−1)75.47e3—75.47e375.47e3k45 (Nm−1)83.40e3—83.40e383.40e3The assignment task prescribes to assign the first natural frequency at 25 Hz and to decrease its sensitivity with respect tom1, m2, and kg1 by 60%, 60%, and 70%, respectively, compared to the original system. All the mode shape entries related to the first natural frequency are free to assume any value. Just to perform the McCormick relaxation, they have been bounded to range between −1 and +1 for numerical purposes; indeed, it should be noted that mode shapes are defined regardless their normalization.Structural modifications are admitted for all the masses as well as all the grounding stiffnesses, leading tox=Δm1,Δm2,Δm3,Δm4,Δm5,Δkg1,Δkg2,Δkg3,Δkg4,Δkg5T; such modifications are bounded in the ranges stated in Table 2. The proposed method is compared with the assignment performed through the PEA method without robustness constraints ([28]), which exploits homotopy optimization for the numerical solution. To ensure a fair comparison, the same assignment specifications, design variables, and constraints, as well as the same solver parameters, are adopted.The optimal modifications of the design variables are reported in Table2: it is interesting noticing that all masses and stiffness have been increased, except for m3 in the robust design. The natural frequencies of the original and modified systems are reported in Table 3. The sensitivities are instead reported in Table 4. The results of Tables 3 and 4 clearly corroborate the effectiveness of the proposed method in solving the RPEA problem. Indeed, exact assignment of the first natural frequency is obtained with the achievement of the required increased robustness. In contrast, the application of the method without robustness requirement leads to a modified system featuring the prescribed natural frequency with small robustness, which is even smaller than the original system.Table 3
Test-case 1: original, desired, and modified system natural frequencies.
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)ωp,1 (Hz)21.8725.0025.0025.00ωp,2 (Hz)32.13—32.6434.33ωp,3 (Hz)42.48—41.2836.84ωp,4 (Hz)52.44—50.5647.05ωp,5 (Hz)64.36—58.3150.91Table 4
Test-case 1: original, desired, and modified system sensitivities of interest.
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)∂ωp,12/∂m1 (rad2s−2 kg−1)−230.9−92.4−322.3−92.4∂ωp,12/∂m2 (rad2s−2 kg−1)−765.2−306.1−1441.3−306.1∂ωp,12/∂kg1 (rad2s−2N−1 m)0.01220.00370.01310.0037To assess the effectiveness of the proposed method, the percentage variation ofωp,1 is investigated in the presence of simultaneous perturbations of m1, m2, and m3 within the range ±25% with respect to their nominal values. The results are shown in Figure 3 where the original system and the two modified ones are compared. Table 5 summarizes the result by revealing that the maximum positive and negative percentage variations are significantly reduced by the robust design, as expected because of the successful assignment.Figure 3
Test-case 1: percentage variation of the first natural frequency in the original and modified systems without and with robustness constraint. Simultaneous perturbations ofm1, m2, and kg1.Table 5
Test-case 1: the original and modified systems without and with robustness: percentage variation of the first natural frequency.
Percentual variations (%) [Δm1; Δm2; Δ kg1]Original system Δωp,1 (%)Mod. nonrobust (method [28]) Δωp,1 (%)Mod. robust (proposed method) Δωp,1 (%)[−25; −25; −25]1.883.140.72[+25; −25; −25]1.192.100.35[−25; +25; −25]−3.60−5.54−1.34[+25; +25; −25]−4.55−6.91−2.02[−25; −25; +25]2.934.231.01[+25; −25; +25]2.643.890.89[−25; +25; +25]−1.87−3.48−0.78[+25; +25; +25]−2.31−4.05−1.00max (Δωp,1(%))2.934.230.89min (Δωp,1(%))−4.55−6.91−2.02A different graphical representation is provided in Figure4, to corroborate the results through a three-dimensional plot. In such a figure, the percentage variations of ωp,1 are analysed by fixing the value of one of the three parameters at a time, while varying the remaining two parameters in the range ±25% with respect to their nominal values. It is evident that the obtained surface is less steep in the case of robust assignment. In contrast, the nonrobust assignment increases the sensitivity with respect to the original system. This result corroborates the effectiveness of the proposed approach.Figure 4
Test-case 1: percentage variation of the first natural frequency in the original and modified systems. (a) With fixedm2, (b) with fixed m1, and (c) with fixed kg1.
(a)(b)(c)A Montecarlo statistical analysis is proposed in Figure5: 106 random simultaneous perturbations of m1, m2, and kg1 in the range ±25% about the nominal values are simulated and the probability density functions of the percentage variations of the assigned natural frequency, Δωp,1, are evaluated. The mean value of Δωp,1 [%], shown through a dashed red line in Figure 5, is equal to −0.16%, −0.29%, and −0.08% for the original system, the nonrobust design, and the robust one, respectively. The mean of the absolute value of Δωp,1 is equal to 1.35%, 2.14%, and 0.51%, respectively. The effect of the robustness improvement can be further appreciated through the evaluation of the interval containing the 95% of values: [−3.0%; +2.6%] in the original system, [−4.8%; +4.4%] in the modified system without robustness, and [−1.2%; +1.0%] in the case of robust design. The intervals containing the 100% of values are instead [−4.5%; +3.0%], [−6.8%; 4.4%], and [−1.9%; 1.0%], respectively. All these results corroborate the effectiveness of the proposed method in reducing the sensitivity of the assigned natural frequency to random simultaneous variations of the parameters of interest.Figure 5
Test-case 1: probability density functions of Δωp,1 obtained through Montecarlo simulations.


### 4.3. Test-Case 2: Robust Assignment of an Antiresonance Frequency in a Five-Mass System

The five-mass system adopted in the previous test is here adopted to assess the method for robust assignment of an antiresonance. The sample test case chosen consists of the robust assignment of the second antiresonance frequency of receptanceh12(jω): such a frequency should be shifted to 40 Hz and its sensitivity with respect to parameters m5 and kg5 should be decreased by 50% and 40%, respectively. The example might represent a practical application where forced vibration of coordinates 2 should be absorbed in the presence of a 40 Hz harmonic excitation of coordinate 1. The following design parameters are assumed, to make a severe test case: x=Δm3,Δm4,Δkg3,Δkg4,Δk34,Δk45T. The same bounds assumed in Section 4.2 are exploited, as reported in Table 6.Table 6
Test-case 2: original and modified system design variables and constraints.
ParameterValueModification boundsModified system without robustness (method [22])Modified system with robustness (proposed method)m3 (kg)8.214[0; 1.080]8.7618.744m4 (kg)2.609[0; 1.080]3.2133.582kg3 (Nm−1)94.26e3[0; 144.9e3]98.44e3207.51e3kg4 (Nm−1)94.26e3[0; 144.9e3]99.85e3105.25e3k34 (Nm−1)75.47e3[0; 144.9e3]75.90e379.91e3k45 (Nm−1)83.40e3[0; 144.9e3]85.27e389.18e3The robust assignment is compared with the nonrobust assignment performed through the method proposed by the authors in [22]: the obtained antiresonances and related sensitivities are, respectively, listed in Tables 7 and 8: both the designs lead to the desired antiresonance frequency. Additionally, the required sensitivities are achieved for the robust assignment together with the prescribed antiresonance frequency. It should be noted that the application of the nonrobust method leads to reduced sensitivities with respect to the original system; however, the increase of robustness is not as large as the one provided the proposed method. The achievement of small sensitivities is obtained in the robust design by means of larger modifications of the two modifiable ground springs. Similar modifications of masses and springs connecting masses are instead obtained by both the designs. These different modifications are difficult to predict and to trace back to intuitive physical explanations: this confirms the need of a numerical approach like the one here proposed.Table 7
Test-case 2: original, desired, and modified antiresonance frequencies ofh12(jω).
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)ωz,1 (Hz)24.68—24.0828.99ωz,2 (Hz)42.4240.0040.0040.00ωz,3 (Hz)64.36—63.0463.23Table 8
Test-case 2: original, desired, and modified systems ofh12(jω).
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)∂ωz,22/∂m5(rad2s−2 kg−1)–16606–8303–10583–8415∂ωz,22/∂kg5(rad2s−2N−1 m)0.23380.14030.16760.1332The effectiveness of the method is further evidenced by Figure6 and by Table 9. The less steep curve is the one achieved through the robust design.Figure 6
Test-case 2: percentage variation of the second antiresonance frequency in the original and modified systems.Table 9
Test-case 2: original and modified systems without and with robustness: percentage variation of the second antiresonance frequency ofh12(jω).
Percentual variations (%) [Δm5; Δ kg5]Original system Δωz,2 (%)Mod. nonrobust (method [22]) Δωz,2 (%)Mod. robust (proposed method) Δωz,2 (%)[−25; −25]0.04−0.35−0.27[+25; −25]−9.24−7.35−5.76[−25; +25]6.214.683.77[+25; +25]−0.030.300.24max (Δωz,2(%))6.214.683.77min (Δωz,2(%))−9.24−7.35−5.76The optimal modified parameters, summarized in Table6, show that a fundamental role in increasing the second antiresonance frequency robustness with respect to the parameters of interest is related to the magnitude of the third grounding spring whose value is more than doubled to cope with the prescribed task. Lastly, also the ratio of the other masses and springs is slightly changed to correctly match the assignment specifications.The Montecarlo analysis with 106 random perturbations within the ±25% variation ranges admitted for m5 and kg5 has been performed. The mean value and the interval containing the 95% of the values of Δѡz,2 [%] are −0.26% and [−6.4%; +6.0%] for the original system, −0.22% and [−5.0%; +4.6%] for the modified system without robustness requirements, and −0.17% and [−4.0%; +3.6%] for the modified system with robustness constraints. The intervals containing the 100% of values are instead [−9.2%; +6.2%], [−7.4%; 4.8%], and [−5.8%; 3.8%], respectively. Finally, the mean absolute variations of the assigned antiresonance frequency are equal to 2.59%, 2.00%, and 1.58%, respectively.


### 4.4. Test-Case 3: Robustification of a Natural Frequency in a System with Continuum and Lumped Flexible Elements

In the third test case, it is required to keep the same natural frequency of the original system, while increasing its robustness with respect to some parameters. This problem could be of interest in resonators (such as feeders [37] or sonotrodes [38]) that are excited by an external tuned harmonic force matching one resonance frequency. Since the system parameters might slightly change during operation, a robust design of these devices is of great importance to ensure that the resonance frequency is robust to parameter variations that could mistune the system and hence downgrade the system performances.Let us consider the system composed of continuum and lumped flexible elements sketched in Figure7. The system parameters are listed in Table 10. EI and ρA denote the linear mass density and the flexural stiffness of the beam, respectively, which has been modelled through four Euler–Bernoulli finite elements.Figure 7
Test-case 3: sketch of the system.Table 10
Test-case 3: original and modified system design variables and constraints.
ParameterValueModification boundsModified system with robustness (proposed method)m1 (kg)25.000[−12.50; 25.00]27.500m2 (kg)25.000[−12.50; 25.00]39.866m3 (kg)25.000[−12.50; 25.00]27.531ρA (kgm−1)15.800[−7.90; 7.90]21.813kl (Nm−1)1.00e3—1.00e3kr (Nm−1)1.00e3—1.00e3k1 (Nm−1)8.70e5—8.70e5k2 (Nm−1)8.70e5—8.70e5k3 (Nm−1)8.70e5—8.70e5kg1 (Nm−1)1.00e5[−0.5e5; 7e5]1.22e5kg2 (Nm−1)1.00e5[−0.5e5; 7e5]7.76e5kg3 (Nm−1)1.00e5[−0.5e5; 7e5]1.24e5EI (Nm2)1.50e5[−0.75e5; 0.75e5]0.77e5The sample RPEA problem here formulated requires maintaining unchanged the third natural frequency and to double the robustness with respect tok1, k2, and k3 (i.e., sensitivities must be reduced by the 50% with respect to the original system). The requirements are summarized in Tables 11 and 12. The admissible modifications together with the lower and upper bounds of the design variables are reported in Table 10.Table 11
Test-case 3: original, desired, and modified system lower natural frequencies.
Original systemDesired valueModified system with robustness (proposed method)ωp,1 (Hz)7.480–7.329ωp,2 (Hz)8.426–11.812ωp,3 (Hz)30.14330.14330.143ωp,4 (Hz)43.716–38.608ωp,5 (Hz)52.034–42.355Table 12
Test-case 3: original, desired, and modified system sensitivities of interest.
Original systemDesired valueModified system with robustness (proposed method)∂ωp,32/∂k1 (rad2s−2N−1 m)0.00500.00250.0025∂ωp,32/∂k2 (rad2s−2N−1 m)0.02320.01160.0116∂ωp,32/∂k3 (rad2s−2N−1 m)0.00500.00250.0025Tables11 and 12 list the obtained natural frequencies and the sensitivities of interest and corroborate the fulfilment of the assignment task, i.e., robustness with respect to the parameters of interest of the third natural frequency (which is the one of particular interests for the system operations) is exactly doubled with respect to the original system design.Figure8 proves the effectiveness of the proposed technique by showing the variations of the natural frequency if the three parameters k1, k2, and k3 are perturbed simultaneously in the intervals ±25% about their nominal values. As expected, the use of the proposed approach makes the system more robust than the original one. Table 13 summarizes the results of the perturbation analysis of the system and corroborates the plots. Indeed, the interval collecting all the values of Δωp,3 [%] is reduced from [−10.9%, 9.4%] in the original system to the smaller one [−5.6% 4.5%] in the system modified through the robust method proposed.Figure 8
Test-case 3: percentage variation of the third natural frequency in the original and modified systems. Simultaneous perturbations ofk1, k2, and k3.
(a)(b)Table 13
Test-case 3: original and modified systems: percentage variation of the third natural frequency.
Percentual variations (%) [Δk1; Δk2; Δk3]Original system Δωp,3 (%)Mod. robust (proposed method) Δωp,3 (%)[−25; −25; −25]−10.92−5.62[+25; −25; −25]−9.11−4.69[−25; +25; −25]0.42−0.15[+25; +25; −25]3.701.44[−25; −25; +25]−9.11−4.69[+25; −25; +25]−6.78−3.35[−25; +25; +25]3.701.44[+25; +25; +25]9.364.54max (Δωp,3(%))9.364.54min (Δωp,3(%))−10.92−5.62The advantages of the proposed method are even more evident in the representation if the parameters of interest are perturbed and the percentage variation of the third natural frequency with respect to its nominal value is considered. Figure9 shows the increase of robustness obtained for the third natural frequency in the modified system in the case of simultaneous perturbations of two parameters, while the third parameter is kept unchanged.Figure 9
Test-case 3: percentage variation of the third natural frequency in the original and modified systems. (a) With fixedk3, (b) with fixed k2, and (c) with fixed k1.
(a)(b)(c)The Montecarlo analysis with 106 random perturbations within the ±25% variation ranges has been performed. The statistical analysis of the random variations of the system parameters leads to the following values of the mean, the absolute mean, and the interval containing the 95% of the values of Δωp,3 [%]: −0.77%, 3.62%, and [−8.7%; +6.9%] in the original system; −0.45%, 1.82%, and [−4.4%; +3.6%] in the modified system with robustness constraints. The intervals containing the 100% of values are finally [−10.8%; +9.3%] and [−5.6%; 4.6%], respectively.


### 4.5. Test-Case 4: Robustification of an Antiresonance Frequency in a Cantilever Truss Structure

In this test case, the proposed method is applied to a more complex system: a cantilever truss structure composed of 27 elements, leading to 26 DOFs. The goal is to discuss the effectiveness of the proposed method whenever it is applied to structures typically employed in civil engineering, that are modelled through a high number of coordinates. Just the main results are reported for brevity.The system is sketched in Figure10. The length of the structure is 5 m, equally divided by 7 members; each member in turn is composed of equal bars made of steel (E = 190 GPa), with uniform mass density (ρ = 7800 kgm−3). The cross-section area of the original system is assumed equal for all the bars, A = 0.0013 m2.Figure 10
Test-case 4: sketch of the cantilever truss structure.Let us assume that a 100 Hz harmonic disturbance is applied at the tip-end of the cantilever structure. The goal is to isolatey26 by imposing an antiresonance frequency for receptance h26, 26(jω). Additionally, to handle variations of the mass placed to node K, it is required to reduce the sensitivity ∂ωz,12/∂mK of −80% with respect to the original system.The design variables include the cross-section area of the upper horizontal and diagonal bars (1, …, 14) and eleven lumped masses that could be attached in nodes A, …, J and L. The constraints on the design variables are listed in Table14. The application of the proposed method leads to the optimal system parameters listed in Table 14, which allow exactly obtaining the desired antiresonance and sensitivity.Table 14
Test-case 4: original and modified system design variables and constraints.
ParameterValueModification boundsModified system with robustness (proposed method)mA (kg)0[0; 10.000]2.429mB (kg)0[0.000; 10.000]5.452mC (kg)0[0.000; 8.000]7.363mD (kg)0[0.000; 8.000]5.167mE (kg)0[0.000; 6.000]1.773mF (kg)0[0.000; 6.000]2.143mG (kg)0[0.000; 5.000]3.987mH (kg)0[0.000; 5.000]3.625mI (kg)0[0.000; 4.000]1.073mJ (kg)0.000[0.000; 4.000]3.595mL (kg)0.000[0.000; 2.000]0.948A1 (m2)0.0013[−0.0009; 0.0009]0.0005A2 (m2)0.0013[−0.0009; 0.0009]0.0014A3 (m2)0.0013[−0.0009; 0.0009]0.0016A3 (m2)0.0013[−0.0009; 0.0009]0.0014A5 (m2)0.0013[−0.0009; 0.0009]0.0011A6 (m2)0.0013[−0.0009; 0.0009]0.0009A7 (m2)0.0013[−0.0009; 0.0009]0.0012A8 (m2)0.0013[−0.0009; 0.0009]0.0008A9 (m2)0.0013[−0.0009; 0.0009]0.0018A10 (m2)0.0013[−0.0009; 0.0009]0.0005A11 (m2)0.0013[−0.0009; 0.0009]0.0015A12 (m2)0.0013[−0.0009; 0.0009]0.0015A13 (m2)0.0013[−0.0009; 0.0009]0.0009A14 (m2)0.0013[−0.0009; 0.0009]0.0014The effectiveness of the proposed method is finally corroborated through the simulation of the steady-state response of the structure excited by the disturbance force, as shown in Figure11, by assuming that the amplitude is 1000 N: the amplitude decreases from 0.140 mm to 0.000 mm in the case of the unperturbed system. The simulation is also provided for the original and modified systems perturbed through a mass perturbation ΔmK = 1.5 kg. Clearly, the effectiveness of the vibration absorption decreases and a small oscillation appears in the modified system; however the vibrations experienced by coordinate y26 are smaller than the one of the original one.Figure 11
Test-case 4: steady-state response of the original and modified systems without and with massmK perturbation.The statistical Montecarlo analysis is performed for 106 random samples by varying mass mK within the range [−1.5; +1.5] kg. The results are reported in Figure 12.Figure 12
Test-case 4: variation of the antiresonance frequency of interest through Montecarlo simulation for the original and modified systems vs massmK perturbation.


## 4.1. Implementation Issues

Four numerical test cases are proposed in this paper. The assignment problems are solved on a 64-bit laptop equipped with a 16 GB RAM and a i7 9750-H, 2.60 GHz, 6 core processor. No parallel computing is exploited, since optimizing the CPU time is not the goal of the proposed work; hence, just one core is adopted in practice.To provide an estimation of the computational effort, let us consider the general scenario involving the assignment ofNp natural frequencies and Nz antiresonances, and let us assume that Ns sensitivities are prescribed for each assigned frequency. Additionally, let us consider the more widespread case where the eigenvectors are let unassigned. Since eigenvectors are defined regardless of their normalization, one component can be arbitrarily assigned (e.g., imposed to 1), as discussed in [22, 41].Bilinear and trilinear terms arise both in the cost function and in the sensitivity constraint. Let us first consider the bilinear termsb due to the inverse eigenvalue problem to perform natural frequency and/or antiresonance assignment. In accordance with equations (7) and (17), the exact computation of the number of the bilinear variables b (that are related to products between a design variable and an eigenvector entry) depends on the structure of matrices ΔMx and ΔKx. Its worst-case upper bound can be estimated by assuming that each design variable in x belongs to all the entries of the modification matrices, and it follows that the upper bound on the cardinality of the b variables is NxNpN−1+NzN−2. Usually, this is an overestimate of the actual cardinality, indeed each design variable just affects few entries of ΔMx and ΔKx.Similarly, by analysing the sensitivity constraint in equations (11) and (19), an upper limit of the number of the bilinear variables v (i.e., those related to products between two eigenvector entries) can be provided since it depends on the structure of ΔMx and ΔKx and on the design variables considered for the sensitivity. Hence, their maximum number is NsNpN−12+NzN−22.Finally, the maximum number of trilinear termst, that depends on the structure of the modification matrices as well, is NxNsNpN−12+NzN−22. Again, this evaluation just provides a worst-case estimation and in practice the number of those terms is significantly smaller, as it can be inferred from Table 1.Table 1
Summary of the computational details of the test cases.
Test-case 1Test-case 2Test-case 3Test-case 4CPU runtime (s)52.120.7137.0180.4Number of DOFs:N551326Number of assigned natural frequencies:Np1010Number of assigned antiresonances:Nz0101Number of assigned sensitivities:Ns3231Number of design variablesNx106825Number of unassigned right eigenvector entriesu431224Number of unassigned right eigenvector entriesl-3-24Number of bilinear variablesb, v12656181Number of trilinear variablest4130135Conversely, the size of the cost function in equation (23) is exactly determined by the number of DOFs of the system and by the number of assigned natural frequencies and/or antiresonances. In the general case, it follows that the number of fh is equal to NpN+NzN−1. Similarly, the number of sensitivity constraints Sh in equation (23) is equal to NsNp+Nz.The computational aspects of each test case proposed in this work are summarized in Table1. The CPU runtime grows if the number of DOFs and the assignment specifications increases. However, the CPU runtime is reasonably small, since IDSM is an offline design technique and does not require real-time computing.The algorithm is implemented in MATLAB through the Yalmip toolbox [59]. Bilinear and trilinear terms are automatically detected, and the variable lifting and McCormick’s convexification are handled through a MATLAB automatic routine. The convex optimization in the first step of the homotopy is solved through the solver Mosek, while in the following iterations, fmincon is adopted. Nonetheless, any numerical solver for linear and nonlinear optimization can be adopted depending on the homotopy stage considered. The Jacobian matrices in equations (11) and (19) have been computed in a fast and accurate way through the complex-step derivatives [48, 60].


## 4.2. Test-Case 1: Robust Assignment of a Natural Frequency in a Five-Mass System

The first test case is a five-mass system widely adopted in the literature as a benchmark for IDSM methods (see, e.g., [16, 23, 61]). The system is sketched in Figure 2, and its parameters are summarized in Table 2.Figure 2
Test-case 1: sketch of the five-mass system.Table 2
Test-case 1: original and modified system design variables and constraints.
ParameterValueModification boundsModified system without robustness (method [28])Modified system with robustness (proposed method)m1 (kg)1.727[0; 1.080]2.3132.751m2 (kg)5.123[0; 1.080]5.7626.202m3 (kg)8.214[0; 1.080]8.8138.214m4 (kg)2.609[0; 1.080]3.1573.671m5 (kg)1.339[0; 1.080]1.8802.417kg1 (Nm−1)94.26e3[0; 144.9e3]140.82e3129.52e3kg2 (Nm−1)94.26e3[0; 144.9e3]106.74e3231.74e3kg3 (Nm−1)94.26e3[0; 144.9e3]173.15e3135.58e3kg4 (Nm−1)94.26e3[0; 144.9e3]126.78e3101.51e3kg5 (Nm−1)94.26e3[0; 144.9e3]115.95e397.429e3k12 (Nm−1)75.14e3—75.14e375.14e3k23 (Nm−1)67.74e3—67.74e367.74e3k34 (Nm−1)75.47e3—75.47e375.47e3k45 (Nm−1)83.40e3—83.40e383.40e3The assignment task prescribes to assign the first natural frequency at 25 Hz and to decrease its sensitivity with respect tom1, m2, and kg1 by 60%, 60%, and 70%, respectively, compared to the original system. All the mode shape entries related to the first natural frequency are free to assume any value. Just to perform the McCormick relaxation, they have been bounded to range between −1 and +1 for numerical purposes; indeed, it should be noted that mode shapes are defined regardless their normalization.Structural modifications are admitted for all the masses as well as all the grounding stiffnesses, leading tox=Δm1,Δm2,Δm3,Δm4,Δm5,Δkg1,Δkg2,Δkg3,Δkg4,Δkg5T; such modifications are bounded in the ranges stated in Table 2. The proposed method is compared with the assignment performed through the PEA method without robustness constraints ([28]), which exploits homotopy optimization for the numerical solution. To ensure a fair comparison, the same assignment specifications, design variables, and constraints, as well as the same solver parameters, are adopted.The optimal modifications of the design variables are reported in Table2: it is interesting noticing that all masses and stiffness have been increased, except for m3 in the robust design. The natural frequencies of the original and modified systems are reported in Table 3. The sensitivities are instead reported in Table 4. The results of Tables 3 and 4 clearly corroborate the effectiveness of the proposed method in solving the RPEA problem. Indeed, exact assignment of the first natural frequency is obtained with the achievement of the required increased robustness. In contrast, the application of the method without robustness requirement leads to a modified system featuring the prescribed natural frequency with small robustness, which is even smaller than the original system.Table 3
Test-case 1: original, desired, and modified system natural frequencies.
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)ωp,1 (Hz)21.8725.0025.0025.00ωp,2 (Hz)32.13—32.6434.33ωp,3 (Hz)42.48—41.2836.84ωp,4 (Hz)52.44—50.5647.05ωp,5 (Hz)64.36—58.3150.91Table 4
Test-case 1: original, desired, and modified system sensitivities of interest.
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)∂ωp,12/∂m1 (rad2s−2 kg−1)−230.9−92.4−322.3−92.4∂ωp,12/∂m2 (rad2s−2 kg−1)−765.2−306.1−1441.3−306.1∂ωp,12/∂kg1 (rad2s−2N−1 m)0.01220.00370.01310.0037To assess the effectiveness of the proposed method, the percentage variation ofωp,1 is investigated in the presence of simultaneous perturbations of m1, m2, and m3 within the range ±25% with respect to their nominal values. The results are shown in Figure 3 where the original system and the two modified ones are compared. Table 5 summarizes the result by revealing that the maximum positive and negative percentage variations are significantly reduced by the robust design, as expected because of the successful assignment.Figure 3
Test-case 1: percentage variation of the first natural frequency in the original and modified systems without and with robustness constraint. Simultaneous perturbations ofm1, m2, and kg1.Table 5
Test-case 1: the original and modified systems without and with robustness: percentage variation of the first natural frequency.
Percentual variations (%) [Δm1; Δm2; Δ kg1]Original system Δωp,1 (%)Mod. nonrobust (method [28]) Δωp,1 (%)Mod. robust (proposed method) Δωp,1 (%)[−25; −25; −25]1.883.140.72[+25; −25; −25]1.192.100.35[−25; +25; −25]−3.60−5.54−1.34[+25; +25; −25]−4.55−6.91−2.02[−25; −25; +25]2.934.231.01[+25; −25; +25]2.643.890.89[−25; +25; +25]−1.87−3.48−0.78[+25; +25; +25]−2.31−4.05−1.00max (Δωp,1(%))2.934.230.89min (Δωp,1(%))−4.55−6.91−2.02A different graphical representation is provided in Figure4, to corroborate the results through a three-dimensional plot. In such a figure, the percentage variations of ωp,1 are analysed by fixing the value of one of the three parameters at a time, while varying the remaining two parameters in the range ±25% with respect to their nominal values. It is evident that the obtained surface is less steep in the case of robust assignment. In contrast, the nonrobust assignment increases the sensitivity with respect to the original system. This result corroborates the effectiveness of the proposed approach.Figure 4
Test-case 1: percentage variation of the first natural frequency in the original and modified systems. (a) With fixedm2, (b) with fixed m1, and (c) with fixed kg1.
(a)(b)(c)A Montecarlo statistical analysis is proposed in Figure5: 106 random simultaneous perturbations of m1, m2, and kg1 in the range ±25% about the nominal values are simulated and the probability density functions of the percentage variations of the assigned natural frequency, Δωp,1, are evaluated. The mean value of Δωp,1 [%], shown through a dashed red line in Figure 5, is equal to −0.16%, −0.29%, and −0.08% for the original system, the nonrobust design, and the robust one, respectively. The mean of the absolute value of Δωp,1 is equal to 1.35%, 2.14%, and 0.51%, respectively. The effect of the robustness improvement can be further appreciated through the evaluation of the interval containing the 95% of values: [−3.0%; +2.6%] in the original system, [−4.8%; +4.4%] in the modified system without robustness, and [−1.2%; +1.0%] in the case of robust design. The intervals containing the 100% of values are instead [−4.5%; +3.0%], [−6.8%; 4.4%], and [−1.9%; 1.0%], respectively. All these results corroborate the effectiveness of the proposed method in reducing the sensitivity of the assigned natural frequency to random simultaneous variations of the parameters of interest.Figure 5
Test-case 1: probability density functions of Δωp,1 obtained through Montecarlo simulations.


## 4.3. Test-Case 2: Robust Assignment of an Antiresonance Frequency in a Five-Mass System

The five-mass system adopted in the previous test is here adopted to assess the method for robust assignment of an antiresonance. The sample test case chosen consists of the robust assignment of the second antiresonance frequency of receptanceh12(jω): such a frequency should be shifted to 40 Hz and its sensitivity with respect to parameters m5 and kg5 should be decreased by 50% and 40%, respectively. The example might represent a practical application where forced vibration of coordinates 2 should be absorbed in the presence of a 40 Hz harmonic excitation of coordinate 1. The following design parameters are assumed, to make a severe test case: x=Δm3,Δm4,Δkg3,Δkg4,Δk34,Δk45T. The same bounds assumed in Section 4.2 are exploited, as reported in Table 6.Table 6
Test-case 2: original and modified system design variables and constraints.
ParameterValueModification boundsModified system without robustness (method [22])Modified system with robustness (proposed method)m3 (kg)8.214[0; 1.080]8.7618.744m4 (kg)2.609[0; 1.080]3.2133.582kg3 (Nm−1)94.26e3[0; 144.9e3]98.44e3207.51e3kg4 (Nm−1)94.26e3[0; 144.9e3]99.85e3105.25e3k34 (Nm−1)75.47e3[0; 144.9e3]75.90e379.91e3k45 (Nm−1)83.40e3[0; 144.9e3]85.27e389.18e3The robust assignment is compared with the nonrobust assignment performed through the method proposed by the authors in [22]: the obtained antiresonances and related sensitivities are, respectively, listed in Tables 7 and 8: both the designs lead to the desired antiresonance frequency. Additionally, the required sensitivities are achieved for the robust assignment together with the prescribed antiresonance frequency. It should be noted that the application of the nonrobust method leads to reduced sensitivities with respect to the original system; however, the increase of robustness is not as large as the one provided the proposed method. The achievement of small sensitivities is obtained in the robust design by means of larger modifications of the two modifiable ground springs. Similar modifications of masses and springs connecting masses are instead obtained by both the designs. These different modifications are difficult to predict and to trace back to intuitive physical explanations: this confirms the need of a numerical approach like the one here proposed.Table 7
Test-case 2: original, desired, and modified antiresonance frequencies ofh12(jω).
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)ωz,1 (Hz)24.68—24.0828.99ωz,2 (Hz)42.4240.0040.0040.00ωz,3 (Hz)64.36—63.0463.23Table 8
Test-case 2: original, desired, and modified systems ofh12(jω).
Original systemDesired valueModified system without robustness (method [28])Modified system with robustness (proposed method)∂ωz,22/∂m5(rad2s−2 kg−1)–16606–8303–10583–8415∂ωz,22/∂kg5(rad2s−2N−1 m)0.23380.14030.16760.1332The effectiveness of the method is further evidenced by Figure6 and by Table 9. The less steep curve is the one achieved through the robust design.Figure 6
Test-case 2: percentage variation of the second antiresonance frequency in the original and modified systems.Table 9
Test-case 2: original and modified systems without and with robustness: percentage variation of the second antiresonance frequency ofh12(jω).
Percentual variations (%) [Δm5; Δ kg5]Original system Δωz,2 (%)Mod. nonrobust (method [22]) Δωz,2 (%)Mod. robust (proposed method) Δωz,2 (%)[−25; −25]0.04−0.35−0.27[+25; −25]−9.24−7.35−5.76[−25; +25]6.214.683.77[+25; +25]−0.030.300.24max (Δωz,2(%))6.214.683.77min (Δωz,2(%))−9.24−7.35−5.76The optimal modified parameters, summarized in Table6, show that a fundamental role in increasing the second antiresonance frequency robustness with respect to the parameters of interest is related to the magnitude of the third grounding spring whose value is more than doubled to cope with the prescribed task. Lastly, also the ratio of the other masses and springs is slightly changed to correctly match the assignment specifications.The Montecarlo analysis with 106 random perturbations within the ±25% variation ranges admitted for m5 and kg5 has been performed. The mean value and the interval containing the 95% of the values of Δѡz,2 [%] are −0.26% and [−6.4%; +6.0%] for the original system, −0.22% and [−5.0%; +4.6%] for the modified system without robustness requirements, and −0.17% and [−4.0%; +3.6%] for the modified system with robustness constraints. The intervals containing the 100% of values are instead [−9.2%; +6.2%], [−7.4%; 4.8%], and [−5.8%; 3.8%], respectively. Finally, the mean absolute variations of the assigned antiresonance frequency are equal to 2.59%, 2.00%, and 1.58%, respectively.


## 4.4. Test-Case 3: Robustification of a Natural Frequency in a System with Continuum and Lumped Flexible Elements

In the third test case, it is required to keep the same natural frequency of the original system, while increasing its robustness with respect to some parameters. This problem could be of interest in resonators (such as feeders [37] or sonotrodes [38]) that are excited by an external tuned harmonic force matching one resonance frequency. Since the system parameters might slightly change during operation, a robust design of these devices is of great importance to ensure that the resonance frequency is robust to parameter variations that could mistune the system and hence downgrade the system performances.Let us consider the system composed of continuum and lumped flexible elements sketched in Figure7. The system parameters are listed in Table 10. EI and ρA denote the linear mass density and the flexural stiffness of the beam, respectively, which has been modelled through four Euler–Bernoulli finite elements.Figure 7
Test-case 3: sketch of the system.Table 10
Test-case 3: original and modified system design variables and constraints.
ParameterValueModification boundsModified system with robustness (proposed method)m1 (kg)25.000[−12.50; 25.00]27.500m2 (kg)25.000[−12.50; 25.00]39.866m3 (kg)25.000[−12.50; 25.00]27.531ρA (kgm−1)15.800[−7.90; 7.90]21.813kl (Nm−1)1.00e3—1.00e3kr (Nm−1)1.00e3—1.00e3k1 (Nm−1)8.70e5—8.70e5k2 (Nm−1)8.70e5—8.70e5k3 (Nm−1)8.70e5—8.70e5kg1 (Nm−1)1.00e5[−0.5e5; 7e5]1.22e5kg2 (Nm−1)1.00e5[−0.5e5; 7e5]7.76e5kg3 (Nm−1)1.00e5[−0.5e5; 7e5]1.24e5EI (Nm2)1.50e5[−0.75e5; 0.75e5]0.77e5The sample RPEA problem here formulated requires maintaining unchanged the third natural frequency and to double the robustness with respect tok1, k2, and k3 (i.e., sensitivities must be reduced by the 50% with respect to the original system). The requirements are summarized in Tables 11 and 12. The admissible modifications together with the lower and upper bounds of the design variables are reported in Table 10.Table 11
Test-case 3: original, desired, and modified system lower natural frequencies.
Original systemDesired valueModified system with robustness (proposed method)ωp,1 (Hz)7.480–7.329ωp,2 (Hz)8.426–11.812ωp,3 (Hz)30.14330.14330.143ωp,4 (Hz)43.716–38.608ωp,5 (Hz)52.034–42.355Table 12
Test-case 3: original, desired, and modified system sensitivities of interest.
Original systemDesired valueModified system with robustness (proposed method)∂ωp,32/∂k1 (rad2s−2N−1 m)0.00500.00250.0025∂ωp,32/∂k2 (rad2s−2N−1 m)0.02320.01160.0116∂ωp,32/∂k3 (rad2s−2N−1 m)0.00500.00250.0025Tables11 and 12 list the obtained natural frequencies and the sensitivities of interest and corroborate the fulfilment of the assignment task, i.e., robustness with respect to the parameters of interest of the third natural frequency (which is the one of particular interests for the system operations) is exactly doubled with respect to the original system design.Figure8 proves the effectiveness of the proposed technique by showing the variations of the natural frequency if the three parameters k1, k2, and k3 are perturbed simultaneously in the intervals ±25% about their nominal values. As expected, the use of the proposed approach makes the system more robust than the original one. Table 13 summarizes the results of the perturbation analysis of the system and corroborates the plots. Indeed, the interval collecting all the values of Δωp,3 [%] is reduced from [−10.9%, 9.4%] in the original system to the smaller one [−5.6% 4.5%] in the system modified through the robust method proposed.Figure 8
Test-case 3: percentage variation of the third natural frequency in the original and modified systems. Simultaneous perturbations ofk1, k2, and k3.
(a)(b)Table 13
Test-case 3: original and modified systems: percentage variation of the third natural frequency.
Percentual variations (%) [Δk1; Δk2; Δk3]Original system Δωp,3 (%)Mod. robust (proposed method) Δωp,3 (%)[−25; −25; −25]−10.92−5.62[+25; −25; −25]−9.11−4.69[−25; +25; −25]0.42−0.15[+25; +25; −25]3.701.44[−25; −25; +25]−9.11−4.69[+25; −25; +25]−6.78−3.35[−25; +25; +25]3.701.44[+25; +25; +25]9.364.54max (Δωp,3(%))9.364.54min (Δωp,3(%))−10.92−5.62The advantages of the proposed method are even more evident in the representation if the parameters of interest are perturbed and the percentage variation of the third natural frequency with respect to its nominal value is considered. Figure9 shows the increase of robustness obtained for the third natural frequency in the modified system in the case of simultaneous perturbations of two parameters, while the third parameter is kept unchanged.Figure 9
Test-case 3: percentage variation of the third natural frequency in the original and modified systems. (a) With fixedk3, (b) with fixed k2, and (c) with fixed k1.
(a)(b)(c)The Montecarlo analysis with 106 random perturbations within the ±25% variation ranges has been performed. The statistical analysis of the random variations of the system parameters leads to the following values of the mean, the absolute mean, and the interval containing the 95% of the values of Δωp,3 [%]: −0.77%, 3.62%, and [−8.7%; +6.9%] in the original system; −0.45%, 1.82%, and [−4.4%; +3.6%] in the modified system with robustness constraints. The intervals containing the 100% of values are finally [−10.8%; +9.3%] and [−5.6%; 4.6%], respectively.


## 4.5. Test-Case 4: Robustification of an Antiresonance Frequency in a Cantilever Truss Structure

In this test case, the proposed method is applied to a more complex system: a cantilever truss structure composed of 27 elements, leading to 26 DOFs. The goal is to discuss the effectiveness of the proposed method whenever it is applied to structures typically employed in civil engineering, that are modelled through a high number of coordinates. Just the main results are reported for brevity.The system is sketched in Figure10. The length of the structure is 5 m, equally divided by 7 members; each member in turn is composed of equal bars made of steel (E = 190 GPa), with uniform mass density (ρ = 7800 kgm−3). The cross-section area of the original system is assumed equal for all the bars, A = 0.0013 m2.Figure 10
Test-case 4: sketch of the cantilever truss structure.Let us assume that a 100 Hz harmonic disturbance is applied at the tip-end of the cantilever structure. The goal is to isolatey26 by imposing an antiresonance frequency for receptance h26, 26(jω). Additionally, to handle variations of the mass placed to node K, it is required to reduce the sensitivity ∂ωz,12/∂mK of −80% with respect to the original system.The design variables include the cross-section area of the upper horizontal and diagonal bars (1, …, 14) and eleven lumped masses that could be attached in nodes A, …, J and L. The constraints on the design variables are listed in Table14. The application of the proposed method leads to the optimal system parameters listed in Table 14, which allow exactly obtaining the desired antiresonance and sensitivity.Table 14
Test-case 4: original and modified system design variables and constraints.
ParameterValueModification boundsModified system with robustness (proposed method)mA (kg)0[0; 10.000]2.429mB (kg)0[0.000; 10.000]5.452mC (kg)0[0.000; 8.000]7.363mD (kg)0[0.000; 8.000]5.167mE (kg)0[0.000; 6.000]1.773mF (kg)0[0.000; 6.000]2.143mG (kg)0[0.000; 5.000]3.987mH (kg)0[0.000; 5.000]3.625mI (kg)0[0.000; 4.000]1.073mJ (kg)0.000[0.000; 4.000]3.595mL (kg)0.000[0.000; 2.000]0.948A1 (m2)0.0013[−0.0009; 0.0009]0.0005A2 (m2)0.0013[−0.0009; 0.0009]0.0014A3 (m2)0.0013[−0.0009; 0.0009]0.0016A3 (m2)0.0013[−0.0009; 0.0009]0.0014A5 (m2)0.0013[−0.0009; 0.0009]0.0011A6 (m2)0.0013[−0.0009; 0.0009]0.0009A7 (m2)0.0013[−0.0009; 0.0009]0.0012A8 (m2)0.0013[−0.0009; 0.0009]0.0008A9 (m2)0.0013[−0.0009; 0.0009]0.0018A10 (m2)0.0013[−0.0009; 0.0009]0.0005A11 (m2)0.0013[−0.0009; 0.0009]0.0015A12 (m2)0.0013[−0.0009; 0.0009]0.0015A13 (m2)0.0013[−0.0009; 0.0009]0.0009A14 (m2)0.0013[−0.0009; 0.0009]0.0014The effectiveness of the proposed method is finally corroborated through the simulation of the steady-state response of the structure excited by the disturbance force, as shown in Figure11, by assuming that the amplitude is 1000 N: the amplitude decreases from 0.140 mm to 0.000 mm in the case of the unperturbed system. The simulation is also provided for the original and modified systems perturbed through a mass perturbation ΔmK = 1.5 kg. Clearly, the effectiveness of the vibration absorption decreases and a small oscillation appears in the modified system; however the vibrations experienced by coordinate y26 are smaller than the one of the original one.Figure 11
Test-case 4: steady-state response of the original and modified systems without and with massmK perturbation.The statistical Montecarlo analysis is performed for 106 random samples by varying mass mK within the range [−1.5; +1.5] kg. The results are reported in Figure 12.Figure 12
Test-case 4: variation of the antiresonance frequency of interest through Montecarlo simulation for the original and modified systems vs massmK perturbation.


## 5. Conclusions

This paper proposes a novel method to perform the robust assignment of natural frequencies and antiresonance frequencies in lightly damped linear vibrating systems. Robustness is obtained by decreasing the sensitivity of the imposed eigenfrequency with respect to variations of some system parameters, i.e., masses and stiffnesses. This problem statement allows tuning correctly the system in the case of nominal parameters, which are assumed as the more likely ones, and reduces the frequency variations under finite perturbations of the considered critical parameters. Besides assigning the eigenfrequency, the related eigenvectors can be partially assigned too.The inverse eigenvalue problem is cast as a constrained least-square problem based on the system matrices and embeds a sensitivity constraint. The sensitivities are formulated by exploiting the system matrices and their Jacobians that are computed through the complex-step derivatives, ensuring quick and precise implementation. A nonconvex, nonlinear optimization problem is obtained, and a reliable numerical solution method is proposed. Such an algorithm is based on the homotopy optimization to boost the converge to a global optimum by exploiting a tight convex approximation of the original cost function by means of a variable lifting of the bilinear and trilinear terms. Constraints are convexified through single- and double-McCormick’s constraints.The effectiveness of the proposed method is assessed through four numerical test cases. A five-mass system where its first natural frequency is assigned to a prescribed value and the increase of robustness is compared with respect to the original system and to the system modified with one of the state-of-the-art methods. In the second test case, the second antiresonance frequency of a prescribed receptance for the same five-mass system is robustly assigned and the results are compared with a method from the literature. In the third test case, the third natural frequency of a system, composed of a flexible beam and some lumped elements, is kept unchanged while its robustness is increased avoiding modification of the parameters involved in the sensitivity analysis. Finally, in the fourth test case, an antiresonance of a prescribed receptance is robustly assigned for a cantilever truss structure with 26 degrees of freedom, thus showing the method applicability to large-dimensional models.The examples show the effectiveness of the proposed method in achieving precise assignment and increasing of robustness, hence overcoming some limitations of the existing literature.



---

*Source: 5593473-2021-06-16.xml*