Response optimization of underactuated vibration generators through dynamic structural modification and shaping of the excitation forces

Resonant vibration generators, such as vibratory feeders or ultrasonic sonotrodes, are often employed in manufacturing to generate harmonic vibrations with suitable amplitude, spatial shape, and frequency, in order to meet the process requirements. These underactuated systems are usually excited in open loop by few actuators, and therefore, it is not ensured that the desired response is correctly achieved, since the feasible motions should belong to the subset of the allowable motions. To achieve the closest approximation of the desired vibrations, some new solutions are here proposed. The first strategy is the optimal shaping of the harmonic forces exerted by the actuators, by solving an inverse dynamic problem through a coordinate transformation and the projection of the desired response onto the subspace of the allowable motion. By exploiting the formulation of such a subspace, a second approach that involves concurrently both the force shaping and the modification of the inertial and elastic system parameters is proposed. The idea of this approach is to exploit the modification of the elastic and inertial parameters to properly shape the allowable subspace in such a way that it spans the desired response. A solution method is developed, and analytical sensitivity analysis is proposed to choose the design variables. Validation is proposed through a linear vibratory feeder with a long flexible tray, taken from the literature. The results show the effectiveness of the proposed strategies that lead to a very precise approximation of the desired response.


Nomenclature b n
Approximation of the nth bilinear term B Force distribution matrix B A Upper-triangular matrix obtained from the QR decomposition of B D Auxiliary matrix for the sensitivity analysis EJ Flexural stiffness of the tray f, f 0 Generalized force vector and force amplitude vector f 0 i , f 0 M ith excitation force amplitude and maximum excitation force amplitude f c Convex approximation of the nonconvex function to optimize F ext Generalized forces for the Lagrange equation of motion G AA , G AU , G UU Inverse of the receptance matrix partitions in the transformed system G AA , G AU , G UU Inverse of the receptance matrix partitions in the transformed system, evaluated for the optimal design parameters I Identity matrix K Stiffness matrix K AA , K AU , K UU Partitions of the stiffness matrix in the transformed system k a i Spring stiffness of the ith actuator K a i Stiffness matrix of the ith actuator K B Stiffness matrix of the tray k l , k r , k Number of design parameters p, p i , p 0 Generic design parameter, ith design parameter and its nominal value p, p opt , p λ = 1 Vector of the design parameters, optimal values and values at the last homotopy iteration P Subspace of the feasible displacements in the transformed system in the case of PFS Desired generalized displacement of interest (set) and desired amplitude vectors x t , x d t 0 , x t 0 Horizontal displacement coordinate of the tray, its desired amplitude and its amplitude y d 0 ; y 0 i Desired amplitude and amplitude of the vertical tray displacement for the ith node of the tray y, y d 0 Desired displacement and amplitudes vectors in the transformed system y a i Vertical coordinate of the ith actuator y A , y A 0 , y d A 0 , e y d A 0 Actuated coordinates displacement, amplitude, desired amplitude and modified desired amplitude vectors y i Vertical coordinate of the ith node of the tray y U , y U 0 , y d U 0 Unactuated coordinates displacements, amplitude and desired amplitude vectors z d 0 Desired amplitude of the tray displacement in the throw angle direction α f , α f i Nominal throw angle and its value at the ith node of the tray Γ Feasible displacement for the constrained PFS ΔG AA, ΔG AU, ΔG UU Partitions of the inverse of the structural modification receptance matrix in the transformed system ΔM Mass structural modification matrix ΔM AA, ΔM AU, ΔM UU Partitions of mass structural modification matrix in the transformed system ΔK Stiffness structural modification matrix ΔK AA, ΔK AU , ΔK UU Partitions of stiffness structural modification matrix in the transformed system θ f Angle defining the direction of the relative actuator displacement  [5,6]. On the other hand, the shape of the vibration of the trough should be consistent with the requirements on the conveying velocity. For example, a uniform flow of the product can be obtained by ensuring uniform vibrations of the trough. Actually, bending of the trough at the frequency of excitation might cause large elastic deformation that severely compromise the part flow, if the design and the control of the feeder do not account for this issue [1]. Other examples are ultrasonic sonotrodes (horns) adopted for welding, atomizing, deagglomerating, machining, cleaning, or cutting [7][8][9][10]. These devices should be designed and excited to vibrate at the prescribed ultrasonic frequency with the proper shape and amplitude [10]. For example, plastic welding and solid-state bonding usually need high uniformity of the tip vibration of the sonotrode to ensure an effective process due to uniform heating and to prevent damage of the sonotrode itself and the products too [8]. These examples show that both a proper excitation force and a proper design of the mechanical part should be adopted to achieve the desired behavior. The difficulties in achieving the desired vibration are exacerbated in the presence of large and flexible devices, such as the long tray of linear feeders, which are excited by just few actuators. This problem has been up to now tackled in the literature through the optimal design of the elastic and inertial properties of the system, in the frame of the "dynamic structural modification" (DSM), by considering the free response of the system, i.e., without considering the excitation forces. A rich literature on DSM has been recently proposed, and several methods to assign the desired modal properties are available (see [11][12][13][14]). In contrast, the optimal shaping of the excitation forces in this kind of systems, as well as its relationship with the mechanical design, has not been investigated.
The optimal force shaping is an open-loop feedforward control. Feedforward control has been usually exploited in the literature in different fields of motion control, especially to overcome bandwidth limitations. The idea of feedforward is to invert the input-output map to compute the input profile leading to the desired output trajectory, by exploiting the system dynamic model. Feedforward is usually developed to track the desired trajectory of the end effector of flexible multibody systems (see [15][16][17][18]), which is defined through just few coordinates. The difficulties in the inversion of such a model are significantly exacerbated in the case of underactuated systems, i.e., when the number of independent actuators is smaller than the number of degrees of freedom (DOFs) of the system, and when the number of prescribed outputs is greater than the number of the actuators. These two situations are typical occurrences in vibration generators, where the response is set to many DOFs that represent different positions of the flexible system, while few actuators are employed.

Contributions of this paper
The analysis of the state-of-the-art highlights that developing a method for shaping the force to be exerted by the actuators is crucial, since the methods in the literature of feedforward control are not appropriate for the case under investigation. On the other hand, the presence of more output specifications than the number of independent control forces restricts the achievable time trajectories of the coordinates of interest, thus often preventing the achievement of the desired system response. Therefore, it is necessary to evaluate the set of allowable responses and to define the best feasible approximation of such a requirement. As for the use of DSM, it is usually addressed by neglecting the forces and are focused on the eigenstructure.
In the light of these needs and of the limitations of the existing literature, this paper proposes four novel approaches for improving the harmonic response of underactuated vibration generators under harmonic open-loop excitations. The goal is to approximate, as much as possible, the desired forced response for all or some coordinates (whose number is greater than the number of actuators), to meet the requirements of the manufacturing process. The first strategy investigated is the proper shaping of the harmonic excitation provided by the actuators. This technique, denoted "force shaping" (FS), computes the best achievable response and the suitable force by solving an inverse dynamic problem through a coordinate transformation and projection of the desired response onto the subspace of the allowable motion. The formulation of such a subspace reveals its relationship with the physical parameters of the system. Hence, a new concurrent approach involving the FS and the mechanical design is proposed. The idea of this approach is to exploit the modification of the elastic and inertial parameters to properly shape the allowable subspace in such a way that it spans the desired response. The benefit of concurrent approaches that handle the design of both the mechanical part and of the feedback control has been already shown to be very effective in the field of eigenstructure assignment together with state feedback control [19][20][21]. This idea is often denoted as "hybrid approach," since the design problem is tackled with a multi-domain perspective and can be interpreted as a model-based mechatronic design [22]. In this work, a novel formulation is proposed, to solve a new problem that is optimizing the forced response through the structural modification and the feedforward control.
All the approaches proposed in this paper, both involving just FS and FS with DSM, are original. The effectiveness and the comparison of the results they provide is proposed through a meaningful test case taken from the literature, that is, the optimization of the response of an industrial vibratory feeder employed in manufacturing plant proposed in [1].

Definitions and statements of the problems
Let us model the vibration generator as an N DOF undamped linear system whose equations of motion are: M ∈ ℝ N × N and K ∈ ℝ N × N are the mass and stiffness matrices, respectively. x ∈ ℝ N is the vector of the displacement of the independent coordinates, and B∈ℝ NÂN F is the force distribution matrix, that is, assumed to be full column rank. The independent external forces are collected in vectorf t ð Þ∈ℝ N F , with N F ≤ N. Damping can be neglected, since vibration generators are usually built with low damping to ensure large vibration with reduced actuation effort.
Let us assume that the system is excited by in-phase harmonic forces f(t) = f 0 cos(ω f t) (f 0 is the amplitude vector, ω f the radian excitation frequency) to obtain the harmonic vibratory motion required by the process. Hence, the steady-state response is x(t) = x 0 cos(ω f t), (x 0 is the amplitude vector). The relation between the amplitude vectors is written in the frequency domain as follows: As briefly mentioned in Section 1.2, two different strategies are adopted in this paper to optimize the system response, in such a way that it feats the desired one: 1. the optimal shaping of the amplitude of the excitation force f 0 to ensure the attainment of the desired displacements for all the DOFs or just a reduced subset of meaningful coordinates of interest. This strategy is denoted "FS." 2. the simultaneous shaping of the amplitude of the excitation force f 0 and modification of the inertial and elastic properties of the system (i.e., of matrices M and K) to ensure the attainment of the desired displacements for all the DOFs or just a reduced subset of meaningful coordinates of interest.
Each strategy, in turn, can be devoted to two different goals: a. the achievement of the desired response for all the coordinates b. the achievement of the desired response for a reduced subset of coordinates.
Combining strategies "1" and "2" together with the goals "a" and "b" leads to four problems solved in the paper: 3.1 FS with fully assigned displacements: the "full force shaping"

Problem statement
The problem consists in shaping the vector of the force amplitudes, f 0 , to ensure the attainment of the desired displacements , which is defined for all the N DOFs of the system. The problem is therefore finding f 0 , ensuring that: All along the paper, superscript d denotes the desired values, while subscript 0 denotes the vector of amplitudes.

Subspace of the allowable motions
If the system is fully actuated (i.e., N F = N), the solution is straightforward since B ∈ ℝ N × N is an invertible matrix: An exact solution exists for any x d 0 . Hence, any displacement vector can be achieved, provided that f 0 belongs to the characteristic curves of the actuators and the actuator bandwidths are greater than ω f .
If the system is underactuated (i.e., N F < N), two critical issues should be tackled. First, B is rectangular, and therefore, it is not possible to calculate its inverse. On the other hand, using the pseudoinverse of B (denoted as B † ) in lieu of B −1 , i.e., calculating the force vector as , just leads to a rough approximation of x d 0 . Second, it should be observed that the motion of the coordinates is constrained to lie on a N F -dimensional subspace, the so-called allowable subspace. These issues impose to develop a proper way to compute f 0 that provides the best approximation of x d 0 among the allowable ones.
In the general case, where B is a sparse matrix with more than N F not-null entries, let us perform the QR decomposition of B as follows: where Q ∈ ℝ N × N is an orthonormal matrix (Q T Q = QQ T = I, I is the identity matrix) and R∈ℝ N ÂN F is upper triangular. R, in turn, can be written by highlighting its upper-triangular, invertible sub-matrix B A ∈ℝ N F ÂN F which collects all its nonnull entries, rank (B A ) = N F . By substituting Eq. (5) into Eq. (2), pre-multiplying it by Q T , and finally introducing the new coordinate vector y = Q T x, Eq. (2) is transformed into the following one: Although y has no physical meaning, it simplifies the problem solution. Indeed, y can be partitioned into N F actuated coordinates, y A ∈ℝ N F , and N − N F unactuated coordinates, y U The model in Eq. (6) is written accordingly, by defining the following partitions of M and Equation (7) can be written in term of the inverse of the receptance matrices It follows that: The second matrix equation in Eq. (9) reveals that the motion of the unactuated DOFs is fully imposed by the motion of the actuated ones: The following transmission matrix can be therefore defined: Equation (11) shows that the feasible displacement vectors are constrained to belong to a subspace spanned by the N F The allowable subspace depends on the inertial and elastic properties of the underactuated subsystem and of the frequency of the excitation force. Additionally, it implicitly depends on matrix B. Indeed, not only the rank of the force distribution matrix imposes the dimension of the allowable subspace but also B itself sets the transformation from the physical coordinates x to the auxiliary ones y, and hence the partitioning into y A and y U . Clearly, the possibility to achieve the desired response is strongly related to the number of independent external forces compared with the number of DOFs (i.e., the so-called degree of underactuation).
By taking advantage of Eq. (10), the first matrix equation in Eq. (9) can be written as a fully actuated N F -DOF subsystem: Since the input matrix B A computed through the QR decomposition is square and full rank, Eq. (13) can be easily inverted, and therefore, any arbitrary desired motion of the actuated DOFs, y d A 0 , can be obtained through the following forces:

Solving strategy
Let us assume that all the N entries of x 0 are required to assume the prescribed values x d 0 . The desired displacement vector is transformed from x to y, leading to the desired amplitudes of the actuated and unactuated DOFs: Imposing y d A 0 in Eq. (14) does not ensure that y U 0 equals y d U 0 too. Therefore, should be projected onto the subspace of the allowable motion spanned by the columns of L to compute the modified desired displacements of the actuated DOFs, denoted e y d A 0 , that should be imposed in Eq. (14) in lieu of y d A 0 to provide the best approximation of x d 0 : Weighed projections can be also performed to reflect differing levels of concern about each entry, i.e., to boost the achievement of the desired displacements for some meaningful coordinates. Additionally, numerical methods for improving the calculation of (L T L) −1 L T can be exploited.
The approach proposed in this section will be denoted for brevity in the following of the paper as full force shaping (FFS).
3.2 FS with partial assignment of the displacements: the "partial force shaping"

Problem statement
In large-scale systems, it is common that the desired displacements are imposed on just some coordinates [10,11]. A different solution of the FS problem should be therefore developed, compared with the one in Section 3.1, leading to the partial force shaping (PFS).
Let us assume that only n d < N harmonic in-phase DOF displacements, x d s ∈ℝ n d (x d s 0 is the amplitude vector), are "set" to assume desired values. In contrast, the remaining N − n d coordinate displacements are "free" to assume arbitrary values and are collected into The goal is to properly define e y d A 0 that should be imposed in Eq. (14) in lieu of y d A 0 to provide the best approximation of x d s 0 .

Subspace of the allowable motions for underactuated systems
By exploiting the coordinate transformation, such as the one in Eq. (15) and Eq. (11), it is possible to write: Equation (17) shows that the allowable x d s 0 are constrained to belong to a reduced subspace spanned by the column of the full-rank matrix This proves that reducing the number of imposed coordinates enlarges the dimension of the achievable subspace and hence boosts the achievement of the desired response for the n d DOFs of interest.

Solving strategy
In accordance with Eq. (18), x d s 0 should be projected onto the subspace of the allowable motion spanned by the columns of P. Hence, the following e y d A 0 should be imposed: Again, the optimal excitation forces f 0 are obtained exploiting Eq. (14), introducing e y d A 0 provided by Eq. (19) in lieu of y d A 0 . Constraints on the admissible values of the unassigned displacements can be also set. In this case, no analytical solution for e y d A 0 can be found. In contrast, a numerical least-square solution of the following problem should be computed (Γ denotes the set of constraints on the feasible displacements for the unassigned DOFs): A 0 If Γ is a convex set, Eq. (20) is a convex minimization problem.

Definitions and problem statement
Both the FS approaches proposed in Section 3 look for the closest approximation of the desired displacement vector by projecting it onto the allowable subspace. This approximation might be inadequate. This limitation can be overcome by modifying the allowable subspace in such a way that it leads to a closer approximation of the desired system response, or, in the best case, it spans such a vector. Whenever the number and the locations of the external forces cannot be modified (i.e., matrix B is not modified), as it often happens in practice, the modification of the elastic and inertial parameters should be exploited.
The DSM problem can be formulated by taking advantage of the relation stated in the second equation in Eq. (9) that shows how the motion of the unactuated DOFs is imposed by the motion of the actuated ones. By defining the following structural modification matrices: then, such a relation can be written for the modified system as follows: The DSM problem exploits Eq. (22). Two different formulations should be adopted for the cases of full or partial assignment of the response.
The mass and stiffness modification matrices, ΔK AU , ΔM AU , ΔK UU , ΔM UU , are obtained by transforming and partitioning the mass and stiffness modification matrices of the model with physical coordinates x, denoted as ΔM and ΔK, through matrix Q: It is assumed, as usually done in the literature on DSM, that the topology of the modification matrices is a-priori assumed, to represent design choices on the admissible modifications. The values of these modifications are, instead, unknowns and are the design variables.

Full FS and dynamic structural modification
If both y d A 0 and y d U 0 are prescribed, the DSM problem is formulated starting from the following equation: Vector p∈ℝ n p collects the n p design variables, i.e., the structural modification parameters, which are constrained to belong to a feasible domain Ψ p .
Since the exact solvability of Eq. (24) is not always ensured, it is convenient to formulate the DSM problem in Eq. (24) as a constrained least-square minimization: Equation (25) is a convex optimization problem whenever Ψ p is convex. Hence, the global optimal solution can be found regardless of the initial guess.
Once the optimal values of the design variables, p opt , are found, it is possible to apply the method outlined in Section 3.1 to compute the optimal excitation forces for the modified system. First, e y d A 0 is find out for the modified system through Eq. (16). Then, f 0 is computed as follows: where: Hereafter, the concurrent FFS enhanced through DSM will be denoted for brevity as FFS-DSM.

Partial FS and dynamic structural modification
Let us consider the case of Section 3.2, and therefore, let us regard the N − n d unassigned coordinates, whose amplitudes are collected in vector x f 0 , as additional unknowns in the DSM problem, together with the design parameters p. This concurrent use of PFS together with DSM is denoted henceforth as PFS-DSM. The desired displacement amplitude vector in the transformed coordinate y is obtained through matrix Q, thus making y d A 0 and y d U 0 functions of the unknown: The unassigned displacement amplitudes are bounded by means of upper and lower element-wise inequalities, i.e., Equation (29) is non-linear and non-convex, due to the presence of some products of two unknowns. Therefore, a proper solving strategy, such as the one proposed in Section 4.3.1, should be adopted to boost the convergence to a global optimal solution. Once Eq. (29) has been solved and the optimal values of the parameter modifications, p opt , and of the unassigned displacements x f 0 are found, then the solving strategy proposed in Section 3.2.3 is exploited. First, e y d A 0 is computed through Eq. (19) with matrix P (and hence L) of the modified system (i.e., adopting the modified system matrices introduced in Eq. (27)). Then, the optimal excitation forces is computed through Eq. (26), by adopting once again the modified system matrices introduced in Eq. (27).

Homotopy optimization
The proposed solution method for Eq. (29) exploits homotopy optimization and is here briefly described. The proposed technique takes advantage of two methods developed by the Authors of this paper for assigning eigenfrequencies and partial mode shapes [11] or antiresonances [23][24][25][26] in vibrating systems. Homotopy is an optimization technique for nonconvex functions that boosts the convergence to a global minimum by solving a sequence of optimization problems that starts with a convex relaxation of the original problem and then morphs it back to the original non-convex one through a path of solutions.
Variable lifting is here exploited to obtain the convex approximation of the bilinear terms in the form p i x f 0 ; j (for some indexes i and j). The problem can be, however, written as a convex function of the auxiliary variables b n (for some indexes n), b n ¼ p i x f 0 ; j . This approach is denoted as variable lifting.
The homotopy map that morphs back to the original problem of Eq. (29) is defined by replacing each bilinear term p i x f 0 ; j with an affine combination of itself and its related auxiliary variable b n , i.e., λp i x f 0 ; j þ 1−λ ð Þb n . The scalar λ ∈ [0; 1] that grows in a discrete way from 0 to 1 is the morphing parameter. In this way, the first problem (λ = 0) consists in solving the convex approximation f c , while the final one (λ = 1) is Eq. (29). In practice, homotopy is adopted to set effective initial guesses for the intermediate non-convex problems. Indeed, at any step of the homotopy transformation, the initial guess is set as the optimal solution of the previous step. The solution obtained at the last step (i.e., ΔG AU (p λ = 1 ) = ΔG AUopt and ΔG UU (p λ = 1 ) = ΔG UUopt ) is the optimal modification of the system parameters.
To provide an effective approximation of the bilinear terms p i x f 0 ; j in all the intermediate steps, the McCormick envelopes are adopted too [27], by imposing upper and lower bounds to p i and x f 0 ; j .

Sensitivity analysis on the system displacements
The numerical efficiency and reliability of DSM algorithms can be increased by removing the design parameters whose variations have a weak effect on the system response (see [28,29]). Sensitivity analysis is often formulated in DSM as the sensitivity of eigenfrequencies or mode shapes with respect to the design parameters. In contrast, this paper proposes the sensitivity analysis of the response of the DOFs of interest with respect to the design parameters. By denoting as p an arbitrary parameter belonging to p (either a mass or a stiffness), sensitivity S p is: where p 0 is the nominal value. By exploiting Eq. (10), the following relation is obtained: The demonstration of Eq. (31) is proposed in Appendix 1, with the analytical expression of D∈ℝ N −N F ð Þ Â N F .

Model of the system
This section applies the proposed methods to a meaningful test case: an industrial linear vibratory feeder. The main part concerning the operations of a linear vibratory feeders is the oscillating tray, where the conveyed products flow. Uniform flow along the tray is obtained by imposing uniform vibration of the tray, despite its flexibility. Vibratory feeders are usually driven by electromagnetic actuators [1,30,31], linear actuators [32,33], and rotating motor with eccentric masses [34], although some new promising techniques are being investigated in the literature, such as actuators based on dielectric elastomers [35]. Regardless of the technology adopted for actuation, it is here assumed that the force amplitude of each actuator can be set to the desired one.
The linear vibratory feeder assumed for validation is sketched in Fig. 1, and its parameters are taken from [1] and summarized in Table 1. Small vibratory feeders are often modeled in the literature as single-DOF [35,36] or twomass [7] systems. In the case under investigation, due to the presence of more actuators and of a long and flexible tray, such a model is not adequate [1]. Therefore, a multi-DOF model with finite elements to model the continuum masses and elasticity should be adopted.
The tray is a flexible beam supported by two vertical grounding springs (whose stiffnesses are k l , k r respectively) and an horizontal grounding spring (k x ). The feeder is excited by three actuators that are represented in the model as three lumped masses (with the same value m a 1 ¼ m a 2 ¼ m a 3 ) connected to the tray by means of springs whose stiffnesses are k a 1 ¼ k a 2 ¼ k a 3 . These springs are adopted and tuned to create a resonance frequency where the tray should vibrate in a quasi-rigid way in the direction of the desired throw angle, i.e., the angle between the speed vector and the tray longitudinal direction. For the same reason, the actuators are linked to the tray with an inclination equal to α f , which is the desired throw angle. In practice, tray flexibility does not allow for this rigid-body ideal behavior, and hence, both the vibration amplitude and the throw angle change along the tray. This is the main cause of non-uniform flow. The problem is exacerbated for long trays. To account for this issue, the tray is modeled through four Euler-Bernoulli finite elements with equal length. In contrast, the beam is assumed to be infinitely rigid in the axial direction; hence, it behaves as a lumped mass m B = ρAl whose horizontal displacement is defined through coordinate x t . The resulting mass and stiffness matrices of the subsystem made by the tray are denoted M B , K B ∈ ℝ 11 × 11 . Three independent coordinates are adopted to model the relative motion of the actuators and the tray: s 1 , s 2 , and s 3 along a direction tilted as the throw angle. These assumptions lead to a 14-dimensional model which includes the 5 nodal translational coordinates of the beam (y), 5 nodal elastic rotations (φ), the horizontal translation of the tray (x t ), and the 3 relative displacements of the actuators (s): x = {y 1 , φ 1 , y 2 , φ 2 , y 3 , φ 3 , y 4 , φ 4 , y 5 , φ 5 , x t , s 1 , s 2 , s 3 } T .
Three independent control forces, f 0 i , are applied. The mass and stiffness matrices (M a i ∈ℝ 3Â3 and K a i ∈ℝ 3Â3 , respectively) of the ith actuator (i = 1, 2, 3), developed in Appendix 2, are:

Specifications
In accordance with the typical requirement set for feeders, the beam should have uniform displacement in the direction of the desired throw angle α f , when excited with a 35-Hz harmonic force.
By assuming an oscillation with prescribed amplitude z d 0 in the throw angle direction, i.e., parallel to the actuators in the rest configuration, equal displacements are imposed for all the nodes: Uniformity is enforced by imposing a rigid motion of the beam with null elastic nodal rotations, i.e., φ d 0 ¼ 0. In the case of FFS, the actuator oscillation amplitude should be assigned too. A reasonable choice is imposing the relative motion s d 0 to be scaled with respect to z d 0 in accordance with a simplified rigid-body model (i.e., a four-mass system). It follows that the amplitudes of the desired displacements at the frequency of operations are collected in x d 0 : It should be noted that changing z d 0 leads to the same scaling of f 0 , since the same scaling of x d 0 is obtained and the system is linear time-invariant. In this test case, let us assume z d 0 ¼ 5 mm and s d 0 ¼ 13:70 mm, which lead to x d 0 stated in the first column of Table 2. A sketch of the system undergoing such desired displacements is reported in Fig. 2 (different scaling are adopted to represent the X and Y axis).
In the case of PFS, the desired coordinates are just those of the tray, while arbitrary motion is allowed for the actuators:

Modal analysis
The modal analysis reveals that there are some vibrational modes close to the excitation frequency ω f , in which the beam undergoes large elastic deformations. The mode shapes u i ∈ ℝ N and the natural frequencies ω i of the 9 lowest-frequency ones are shown in Fig. 3 (by assuming mass normalization). These modes are excited by the actuator, by perturbing the desired shape of the tray vibration. Attention should be paid to the 6th and the 7th vibrational.
In particular, if the tray is rigid, the 6th mode shape is the one that best describes the motion of interest, i.e., the relative opposite motion between the tray and the actuators in the throw angle direction.
In the light of the results of the modal analysis, the basic idea of the proposed force shaping can be therefore interpreted as finding the optimal force vector that optimally compensates for these unwanted elastic behaviors, leading to the best achievable response. As for the hybrid approaches involving FS and DSM, they can be interpreted in finding the optimal feasible modifications that allow the force to properly compensate for these elastic behaviors.

Application of FFS and PFS
FFS and PFS are here applied to the test case and compared with the results provided by the application of three equal forces, as often done by practitioners in commanding vibratory feeders whenever a uniform response along the tray is wanted [1]. It should be noted that equal forces are also the ones computed by solving Eq. (3) through the pseudoinverse f MþK x d 0 , and are optimal in the case of rigid beams. This widespread benchmark approach will be named "PsI method" (PsI) in tables and in figures.
Performances are evaluated through these parameters, listed in Table 2 and Table 3: & the cosine between the obtained displacements and the desired ones that should approach 1. Such a value is computed for just the tray displacements, i.e., cos x d s 0 ; x s 0 ; & the maximum difference among the vertical displacements of the tray, max y 0 i À Á −min y 0 i À Á , and the maximum elastic rotations max φ 0 i À Á that provide two measures of the elastic deformations; & the throw angles for each node of the tray α f i ¼ tan −1 y 0 i = x t 0 j j À Á that should be 20°, as well as the maximum difference among the throw angles of the tray, max α f i À Á −min α f i À Á The steady-state forced responses obtained by using the methods discussed in this section are reported in Table 2 and depicted in Fig. 4. It is evident, by comparing all the evaluation parameters, that PFS closely approximates the desired response, which is instead very roughly approximated if the benchmark is adopted. The displacement of the tray is a quasi-rigid motion, with negligible elastic rotations and high uniformity. A minor improvement is obtained through FFS, which exploits the simplified rigid-body model as a heuristic rule to set s d 0 . The comparison of the throw angles in Table 3 confirms the results, since just small variations about the target are obtained through PFS.
The force vectors leading to such results are listed in Table 4, by also normalizing the harmonic excitation amplitudes with respect to their maximum value, i.e., f 0 M ¼ max f 0 1 À ; f 0 2 ; f 0 3 Þ. While the PsI method employs three equal Fig. 2 Sketch of the desired displacement amplitudes in the case of FFS forces, the two proposed methods lead to different amplitudes. In the case of PFS, forces with 180-degree relative phases are obtained too, to compensate for the beam flexibility. Clearly, it causes a higher control effort, as corroborated by ‖f 0 ‖ 2 . The results are analyzed through the modal participation factor (MPF) of each vibrational mode in the forced response. Indeed, when a linear system is excited by harmonic forces at the frequency ω f and whose amplitudes is f 0 ∈ℝ N F , its displacements are a linear combination of the mode shapes u i . The coefficients of this linear combination are the MPFs: The MPFs of the three methods are reported in Fig. 5, where each value has been normalized through the maximum  participation factor for ease of comparison: the 6th and the 7th modes are those that mainly contributes in the case of PsI and FFS. These mode shapes involve large elastic deformations of the tray (see Fig. 3), thus causing the non-uniform displacement.
In contrast, PFS excites remarkably other vibrational modes, where the tray experiences a quasi-rigid motion, such as the 4th and the1st modes. The higher contributions of these modes partially compensate for the large displacement at the central point of the beam due to the 6th mode. This new mix of mode shapes ensures therefore the desired quasi-rigid motion of the tray.

FS with structural modification
Fourteen design variables are initially chosen (p ∈ ℝ 14 ), including the actuator masses that can be modified by adding or removing counterweight masses [1]; 5 additional nodal lumped masses m i , i = 1, …, 5, to be placed to the tray; and the modification of all the lumped springs: In contrast, it is assumed that the tray flexural stiffness and linear mass density cannot be modified, although the method can handle their modifications, too. Clearly, in the case of an existing tray, it is very difficult modifying such parameters. On the other hand, including these parameters would allow for better results. A sketch of the feeder highlighting the chosen design variables is reported in Fig. 6. Constraints on the feasible modifications are imposed, as stated in the first column of Table 6. Additionally, a constraint on the overall mass increase to be lower than 15 kg is set.

Sensitivity analysis
The sensitivity analysis is here applied to discard the design variables with the lowest sensitivity. In the case of the nodal masses m i , i = 1, …, 5 whose initial values are 0, the value adopted for normalization p 0 is assumed as the upper bound of the feasible modification. The sensitivity analysis for PFS is reported in Table 5.
To show the correctness and the usefulness of the sensitivity analysis, PFS-DSM will be evaluated considering both the complete set of available structural modifications reported in Eq. (38) and a reduced set not including m a 1 , m a 3 , and k x , which are the masses and stiffness with the smallest sensitivities. The test with a lower number of design parameters will be denoted as PFS-DSM 2 ( Table 6).

Application of FFS-DSM and PFS-DSM
The methods have been applied with the same requirements given for the FS alone. Additionally, as required by the use of the McCormick envelopes, the actuator displacements s 0 i are constrained as follows: −20 mm ð Þ≤ s 0 i ≤ 20 mm ð Þ. Columns from the 3rd to the 5th ones report the optimal modifications computed through FFS-DSM, PFS-DSM, and PFS-DSM 2. The displacements obtained are listed in Table 2.
In the case of FFS, the concurrent use of DSM leads to a considerable improvement, as corroborated by the increase of cos x d s 0 ; x s 0 from 0.6526 to 0.9873, the decrease of max y 0 i À Á −min y 0 i À Á from 3.29 to 0.62 (mm), and the decrease of max φ 0 i À Á from 0.0023 to 0.0005 (rad). In the case of PFS, the concurrent use of DSM leads to a smaller improvement compared with FFS, since PFS alone provides very good performances. Nonetheless, the improvement is evident by the cosine that approaches 1 and by comparing max y 0 i À Á −min y 0 i À Á that decreases from 0.38 to 0.05 (mm) and max φ 0 i À Á that decreases from 0.0003 to 0.0001 (rad). In practice, the tray behaves as a rigid beam under the excitation of the computed optimal force, despite its flexibility.
In the case of PFS-DSM 2, i.e., with reduced design variables, just a negligible reduction of the performances  is obtained, as expected: the partial cosine decreases from 0.9998 to 0.9990. The elastic rotations are almost unchanged, while the maximum variation of the vertical displacements along the tray increases for just 5%. On the other hand, a simpler set of modification has been employed. Figure 7 corroborates the results showing the system displacements obtained in the cases of concurrent DSM and FS.
The analysis of the throw angles along the tray, proposed in Table 3, confirms these results. High uniformity is obtained in both the strategies of partial FS with DSM, with negligible differences between the maximum and minimum angles (0.6°and 0.7°for PFS-DSM and PFS-DSM2 respectively).
The analysis of the optimal forces in Table 4 reveals that FFS-DSM and PFS-DSM employ three forces having equal phases and different amplitudes. In contrast, PFS-DSM 2 leads to forces with different amplitudes and 180degree relative phases. FFS-DSM leads to higher actuation forces, as corroborated by ‖f 0 ‖ 2 , while PFS-DSM 2 enables the lowest actuator effort among the hybrid strategies.
Finally, it should be noted that all the constraints on the design variables are satisfied. The possibility to handle different constraints is another strength of the proposed methods.

Eigenstructure and modal participation factors in the modified systems
The modal participation factors of the mode shapes of the modified systems, excited with the optimal forces, are shown in Fig. 8. It highlights that DSM notably alters the participation of each mode. It is interesting to notice that elimination of some design variables in PFS-DSM 2 leads to a significantly different mix of mode shapes involved in the forced response.
Finally, the presence of several vibrational modes that significantly contribute to the forced response confirms that the formulation of DSM problem through the subspace of the allowable motions of the forced response, in lieu of the traditional assignment of the vibrational modes, is more straightforward.

Response of the modified system to three equal forces
Finally, it should be noted that these structural modifications are not enough to obtain the desired response when used alone. Indeed, if three equal forces are applied to the modified systems, as those of the PsI method, a bad uniformity is obtained, as

Conclusions
This paper proposes four novel strategies for improving the response of underactuated vibration generators operating under harmonic, open-loop excitation. The goal is to assign the desired steady-state response to all the coordinates (the "full assignment") or to a subset of meaningful ones of interest (the "partial assignment") to accomplish the requirements of the manufacturing process. Two strategies are exploited, leading to 4 methods: the optimal shaping of the excitation force (FS) and the simultaneous use of FS and the modification of the inertial and elastic properties of the system (DSM). Underactuated systems are investigated, where the number of specifications is greater than the number of independent control forces.
Since the input matrix can be sparse, its QR decomposition is employed to transform the model in a convenient set of coordinates. Although this representation loses physical meaning, it defines the subspace of the allowable motions. The dimension of the basis of such a space depends on the number of actuators and of coordinates whose motion is not imposed. Consequently, partial assignment solved through the proposed PFS enables a better approximation for the coordinate of interest, compared with the FFS where the motion of all the coordinates is imposed.
An even closer approximation of the desired response is obtained by modifying the subspace of the allowable motions, in such a way that it could span the desired response or, at least, provide a closer approximation of it. Two "hybrid" approaches are proposed: FFS-DSM and PFS-DSM. A strategy to solve the concurrent design of the control forces and of the modifications of the inertial and elastic parameters is proposed for both the full and the partial assignments. Constraints on the design variables can be introduced. Analytical formulation of the sensitivity analysis is proposed too, to discard the parameters that weakly affect performances.
The method validation is proposed through a linear vibratory feeder, taken from the literature, whose 14 DOFs are controlled through 3 independent forces, and where specifications are set on the motion of the flexible tray. The application of the 4 techniques that are compared with a common benchmark used by practitioners shows the benefits of the proposed methods. In particular, PFS-DSM provides the best approximation of the desired response and makes the system behave as a "quasirigid" system, despite the presence of relevant flexible dynamics of the tray.
Funding Open access funding provided by Università degli Studi di Padova within the CRUI-CARE Agreement. The third author acknowledges the financial support of the Cariparo foundation ("Fondazione Cassa di Risparmio di Padova e Rovigo") through a Ph.D. scholarship.
In order to provide By exploiting once again the equation for the derivative of the inverse of a matrix, already adopted in Eq. (45) and Eq. (46), it is possible to obtain that: where, in turn, the derivative ∂L ∂p is: The second block of the matrix on the right-hand side of Eq. (50) has been already computed in the previous equations. Hence, the substitution of Eq. (49) and Eq. (50) in Eq. (48) leads to: In the case of PFS, by exploiting the derivative of Eq. (19), the following equation holds: where the derivative ∂P ∂p is: The derivative ∂L ∂p is the same already provided in Eq. (50). To sum up, the derivative components of the generalized forces applied to the system are collected in F ext .
To compute the energy of the ith actuator, it is convenient to define the absolute displacements of the actuator through the generalized coordinates: Hence, the kinetic energy and potential energy of the actuator are respectively: Finally, by computing the derivatives, reported in the Lagrange equation of motion, it is possible to obtain the equation of motion for each actuator: Recasting in a matrix form the equation of motions leads to the mass and stiffness matrices of Eq. (32).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.