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

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

#### 2. Method Formulation

##### 2.1. Definitions

Let us consider a *N*-degrees of freedom (DOFs) undamped vibrating system with denoting the displacement vector ( is the acceleration vector). The system is modelled through its symmetric positive definite mass matrix **M** = **M**^{T} > 0 , and its positive semidefinite and symmetric stiffness matrix **K** = **K**^{T} ≥ 0. The external forces are denoted by . The equation of motion is therefore

The model parameters collected in **M** 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 (*ω*):and the receptance matrix is introduced, such that . The *r, c* entry of , denoted by *h*_{rc}(*ω*), 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]:where and 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 the *N* natural frequencies (*i* = 1, … , *N*) together with their related mode shapes (*i* = 1, …, *N*) by solving the following eigenproblem [28]:

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]where *N*_{p} denotes the number of eigenpairs to be assigned, i.e., the pairs made by the natural frequency and its related mode shape . 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]:where are weighing coefficients and 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 (*n*_{x} is their cardinality), to be computed for achieving the desired eigenstructure. The feasible domain is introduced too, which includes upper and lower bounds (denoted by *x*^{L} and *x*^{U}, 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 : . 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:

A new feasible domain is defined as the union of the constraints on and , .

In writing equation (7) and in defining , it has been implicitly assumed that all the entries of are unassigned, to obtain a simpler formulation; the extension to the partial assignment of just some entries of 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 the* n*_{p} parameters affecting **M** and **K** into vector . The sensitivity the of *i*-th natural frequency with respect to the *j*-th parameter , denoted by for brevity, is computed as follows [50]:

*Proof. *Let us consider the natural frequency eigenproblem in equation (4) and let us assume that **M** and **K** depend on an uncertain parameter , i.e., and . Developing the partial derivatives leads to the following relation:By premultiplying equation (9) by , it yieldsSince **M** and **K** are symmetric, then . 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 and (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 . 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 prescribed *N*_{p} eigenpairs together with the sensitivities of the natural frequencies with respect to some critical parameters (whose cardinality is denoted by *N*_{s}). Such sensitivities , for the indexes *i* and *j* of interest, are required to assume the prescribed values making the natural frequency variations admissible as the parameter changes within an expected range of variation. Indeed, represents the local, first-order approximation of the effect of each parameter on the natural frequencies. The smaller the is, the more robust the assigned is with respect to changes of .

In the case of the system modified through **ΔM** and **ΔK**, the sensitivity is inferred from equation (8) by taking into account the structural modification matrices:

The inclusion of requirements on 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

collects the specification on the sensitivities of with respect to the *N*_{s} parameters, written in nonfractional form:where

By including it within the inverse eigenvalue problem, the following problem is obtained: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 Sections 2.2 and 2.3 can be extended to perform robust assignment of antiresonances, by considering that antiresonance frequencies together with their related right eigenvectors can be computed by solving the eigenproblem of the adjunct system:

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):where is the *i*-th prescribed antiresonance frequency and 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 could be used, in principle, to define set interval specifications on the antiresonant response [41, 51]. However, imposing some entries of usually makes the achievement of the desired antiresonance very challenging. As a matter of fact, just one paper in the literature assigns through DSM [52]. In contrast, it is always neglected (see, e.g., [22, 53]). The difficulties in assigning 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 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, , are computed as follows [51]:where and 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 , 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 sensitivity is defined as follows:

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 (Section 3.2)(2)Definition of the homotopy map (Section 3.3)(3)Introduction of bilinear and trilinear McCormick’s constraints (Section 3.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 form (for proper indexes and ) and variables may include entries of as well as of . Bilinear terms are replaced in the convex approximation of the cost function by auxiliary variables denoted by , i.e.,

Secondly, the sensitivity constraints introduce bilinear terms of the form (for proper indexes *o* and ) that are replaced in the convex approximation of by the auxiliary variables :

Lastly, trilinear terms also appear in the sensitivity constraint as (for proper indexes *o*, , and ), and without loss of generality, may include entries of as well as entries of . These terms are replaced in the convex approximation of (both and ) by the auxiliary variables :

The proposed variable lifting approach leads to the convex approximations of the nonconvex cost function (denoted by *f*_{nc} and convexified through *f*_{c}) and of the sensitivity constraint **S** (convexified through **S**_{c}) in equation (15), by lifting the optimization problem to a higher-dimensional one, which is linear in the auxiliary variables , , and (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: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 *f*_{c} and **S**c is solved; in the last one (i.e., *λ* = 1), the original nonconvex problem defined by *f*_{nc} and **S** is solved. In the intermediate steps, the problems defined by *f*_{h}(*λ*) and **S**_{h}(*λ*) are solved.

The initial guess of the nonconvex minimization of in the *i*-th iteration is the optimal solution of the previous step, i.e., . Finally, the values of the optimal design variables are computed as and . In this way, a path of solutions converging to the optimal one is obtained.

In this work, *f*_{h} and are obtained by replacing each bilinear or trilinear terms with an affine combination of themselves and their related auxiliary variable, i.e.,

A flowchart of the method is proposed in Figure 1.

##### 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]):where the upper and lower bounds for are denoted by and , respectively, while the for , the upper and lower bounds are and . For brevity, the trivial definition of bilinear McCormick’s constraints for the auxiliary variable that lifts 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 *n*_{k}, such that , is introduced; then, the following inequalities are defined:

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 *n*_{k}, 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 of *N*_{p} natural frequencies and *N*_{z} antiresonances, and let us assume that *N*_{s} 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 terms ** b** 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

**(that are related to products between a design variable and an eigenvector entry) depends on the structure of matrices and . Its worst-case upper bound can be estimated by assuming that each design variable in**

*b***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 . Usually, this is an overestimate of the actual cardinality, indeed each design variable just affects few entries of and .

Similarly, by analysing the sensitivity constraint in equations (11) and (19), an upper limit of the number of the bilinear variables (i.e., those related to products between two eigenvector entries) can be provided since it depends on the structure of and and on the design variables considered for the sensitivity. Hence, their maximum number is .

Finally, the maximum number of trilinear terms ** t**, that depends on the structure of the modification matrices as well, is . 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.

Conversely, 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 *f*_{h} is equal to . Similarly, the number of sensitivity constraints *S*_{h} in equation (23) is equal to .

The computational aspects of each test case proposed in this work are summarized in Table 1. 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.

The assignment task prescribes to assign the first natural frequency at 25 Hz and to decrease its sensitivity with respect to *m*_{1}, *m*_{2}, and *k _{g}*

_{1}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 to ; 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 Table 2: it is interesting noticing that all masses and stiffness have been increased, except for *m*_{3} 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.

To assess the effectiveness of the proposed method, the percentage variation of *ω*_{p,1} is investigated in the presence of simultaneous perturbations of *m*_{1}, *m*_{2}, and *m*_{3} 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.

A different graphical representation is provided in Figure 4, 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.

**(a)**

**(b)**

**(c)**

A Montecarlo statistical analysis is proposed in Figure 5: 10^{6} random simultaneous perturbations of *m*_{1}, *m*_{2}, and 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.

##### 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 receptance *h*_{12}(*jω*): such a frequency should be shifted to 40 Hz and its sensitivity with respect to parameters *m*_{5} and 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: . The same bounds assumed in Section 4.2 are exploited, as reported in Table 6.

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

The effectiveness of the method is further evidenced by Figure 6 and by Table 9. The less steep curve is the one achieved through the robust design.

The optimal modified parameters, summarized in Table 6, 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 10^{6} random perturbations within the ±25% variation ranges admitted for *m*_{5} and 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 Figure 7. 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.

The sample RPEA problem here formulated requires maintaining unchanged the third natural frequency and to double the robustness with respect to *k*_{1}, *k*_{2}, and *k*_{3} (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.

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

Figure 8 proves the effectiveness of the proposed technique by showing the variations of the natural frequency if the three parameters *k*_{1}, *k*_{2}, and *k*_{3} 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.

**(a)**

**(b)**

The 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. Figure 9 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.

**(a)**

**(b)**

**(c)**

The Montecarlo analysis with 10^{6} 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 Figure 10. 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 m^{2}.

Let us assume that a 100 Hz harmonic disturbance is applied at the tip-end of the cantilever structure. The goal is to isolate *y*_{26} by imposing an antiresonance frequency for receptance *h*_{26, 26}(*jω*). Additionally, to handle variations of the mass placed to node K, it is required to reduce the sensitivity 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 Table 14. 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.

The 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 Figure 11, 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 Δ*m*_{K} = 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 *y*_{26} are smaller than the one of the original one.

The statistical Montecarlo analysis is performed for 10^{6} random samples by varying mass *m*_{K} within the range [−1.5; +1.5] kg. The results are reported in Figure 12.

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

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by Università degli Studi di Padova, Program BIRD 2020, project code: BIRD205841/20, “Active and Passive Approaches for the Optimization of the Dynamic Response of Flexible Systems.”