Integrated Inverse Dynamics and Optimized Mechanical Design in Underactuated Linear Vibratory Feeders Under Periodic Excitation

This paper proposes an integrated method for optimizing the response of underactuated linear vibratory feeders operating in open-loop control, under generic periodic excitations. The goal is ensuring a uniform motion of the tray, despite the presence of less actuators than degrees of freedom and of several specifications of the desired motion. To cope with the underactuated nature of these systems and with their non-minimum phase behavior, dynamic structural modification and the inverse dynamics approach are properly integrated by exploiting a common definition of the system internal dynamics. In the inverse dynamics problem, the inverse dynamics is stabilized through output redefinition and the resulting ordinary differential equations are integrated to compute causal actuation forces, ensuring almost-exact tracking for as many coordinates as the number of actuators. The tracking of the remaining coordinates of interest is improved through a proper design of the mechanical parameters, based on the modification of the internal dynamics. The effectiveness of the proposed method is assessed through numerical simulations performed on the challenging case of a 14-degrees of freedom underactuated non-minimum phase linear vibratory feeder adopted in manufacturing plants to convey products. The results evidence the benefits obtained by integrating structural modification together with inverse dynamics. Inverse dynamics is effective, since the tracking error on the imposed coordinates is negligible. On the other hand, the benefits introduced by dynamic structural modification are proved as well, by the reduction of the tracking error also for the non-imposed coordinates.


Motivations and State of the Art
Vibratory feeders are widely used in the industry to convey products in manufacturing plants by exploiting mechanical vibrations [1].Linear vibratory feeders convey the products along a linear path [2] and the flexibility of these devices is exploited to achieve a large amplitude of vibration through a low actuation effort [3].In general, the number of actuators is smaller than the number of degrees of freedom (DOFs) arising in the system due to its flexibility, leading to an underactuated system.It is well known in the literature that the control of underactuated systems is not trivial, in particular when they are operated in open loop, as it is usually done for linear vibratory feeders.Hence, it is fundamental to carefully plan a-priori the actuation forces to achieve high 1 3 performances [4].Open loop control is usually adopted in feeders (as well as other resonators or vibration generators), because these systems, with the goal of reducing their cost, are typically not equipped with sensors, with microprocessors implementing feedback control, or with actuators whose forces can be modified with a high rate and in real-time control as commanded by a controller.
Vibratory feeding is often performed by means of a sinusoidal vibration of the surface employed to convey the products (see, e.g., [5,6]); however, harmonic excitation does not always allow achieving high conveying velocity [7].It has been proved that non-sinusoidal vibrations might provide high conveying speed as well as other beneficial features such as reliable sorting and orientating motion [8].
The development of feedforward and motion planning algorithms for underactuated multibody systems is not trivial due to the presence of a rectangular force distribution matrix that exacerbates the difficulties in model inversion.In the case of non-minimum phase systems, model inversion is not trivial; indeed, zeros with positive real part arise, and once the model is inverted, those become right-half complex plane poles leading to unstable internal dynamics that does not enable to compute the actuation forces, since they diverge.
Trajectory tracking is tackled in [9] through an inverse dynamics approach together with the Byrnes-Isidori regulator, a non-causal solution is obtained, i.e., pre-actuation is required, to track the desired trajectory.More recently, the servo-constraint method has been introduced in [10] and applications to meaningful test cases have been proposed in [11,12].Such technique represents the trajectory to be tracked as an algebraic constraint that transforms the system model into a set of differential algebraic equations.Optimization-based approaches have been proposed in [4] where the model inversion approach through the non-linear input-output normal form and the servo-constraints are analyzed.Other formulations based on the solution of twopoint boundary value optimization problems are proposed to tackle both the inverse dynamics and motion planning; see, e.g., [13][14][15] and the references therein.Recently, an inverse dynamics method for non-linear systems that exploits the non-linear output redefinition has been proposed by the authors in [16] to compute the actuation forces, and the method has been adopted in [17] to perform the motion planning.
In the case of highly underactuated systems, the number of non-actuated coordinates is much larger with respect to the number of independent actuators.In this scenario, achieving satisfying results is not trivial.This is the typical scenario of long and flexible linear vibratory feeders employed in manufacturing plants where several coordinates are adopted to model the flexible tray over which products flow, while a limited number of actuators are installed.In this case, dynamic structural modification (DSM) can be employed to improve the system performances.DSM methods compute the optimal structural modifications (mass, damping and stiffness system parameter modifications) that enable to achieve the desired specifications.Traditionally, structural modification has been widely adopted over the decades to assign natural frequencies [18][19][20], mode shapes [21,22] and antiresonance frequencies [23][24][25].Recently, DSM has been adopted also to increase the robustness against the system uncertainties in motion planning [26] and to improve the subspace of the allowable motion in linear vibratory feeders excited through harmonic excitation [27].In addition, a mechanical design approach to convert non-minimum phase systems into minimum phase has been proposed in [4].

Contributions of the Paper
This paper proposes the general mathematical formulation for the inverse dynamics of non-minimum phase underactuated mechanical systems, where QR decomposition is exploited to transform the system into the actuated-unactuated model.Then, the internal dynamics is characterized exploiting the definition of the output reference in terms of the actuated and unactuated coordinates.Further, since the system is non-minimum phase, a strategy to stabilize the internal dynamics (relating the desired output trajectory and the unactuated coordinates) is proposed through the output redefinition method.In particular, the paper extends the inverse dynamics method for linear vibratory feeders proposed by the authors in [27] and [28] in the simpler case of a sinusoidal trajectory reference, to the more general case of a generic periodic reference trajectory.Conversely to the above-mentioned papers where an algebraic problem was solved to compute the actuation forces, here, a differential-algebraic problem is adopted because of the presence of the internal dynamics that is governed by ordinary differential equations.
The manuscript also provides a novel structural modification strategy compared to the usual DSM approaches proposed in the literature, that is therefore another novelty introduced by this paper.Indeed, the proposed approach directly exploits the internal dynamics (formulated for the undamped system, and in the presence of harmonic excitations), thus considering the forced response in the presence of the reference trajectories, to compute the optimal mechanical design for the system that enables to increase the tracking error performances of the whole system.
The proposed method is applied to the challenging case of a non-minimum phase underactuated linear vibratory feeder employed in manufacturing plants to convey products and the results confirm the effectiveness of the inverse dynamics 1 3 algorithm, of the DSM approach, and of the overall approach proposed in this paper.
The proposed method is of interest under a theoretical point of view, because of the challenges in controlling underactuated systems that are here exacerbated by the presence of more unactuated variables than the actuated ones and because of the unstable internal dynamics.Additionally, it is of interest for several practical applications: the overall approach is an effective tool for the optimal concurrent mechanical design of the feeders (through the DSM) and of the actuation forces (through the inverse dynamics) to improve the feeding capabilities of these widely adopted devices, as well as for other vibration generators sharing similar dynamic models.

Dynamic Model
Let us consider the linear vibratory feeder, like the one sketched in Fig. 1, employed in manufacturing plants to convey products (see, e.g., [3]).The tray of the feeder is modeled as a flexible beam through N b Euler-Bernoulli planar beam elements, as proposed in [28].The vertical displacements and the nodal rotations of the tray are denoted as y i and i , respectively (i = 1,…,N b + 1).The horizontal displacement of the tray is x t , and the beam is assumed as axially stiff, thus just one coordinate is needed in this direction.Translations that are orthogonal to the xy plane of motion are not considered.The tray is actuated through N A independent actuators modeled as suspended masses (as usually done in linear feeder [28]), whose absolute translations are x a i = x t + s i cos f and y a i = y i + s i sin f (i = 1,…,N A ). y i is the tray vertical translation in the point of attachment of the actuators and s i is the tray-actuator relative displacement along the direction defined by the angle f = 2 − f , with f being the throw angle, that is a key parameter in the feeder design to ensure the desired flow.The model of the actuators, whose mass is m a i and k a i is the stiffness of the leaf spring connecting it to the tray, is here omitted for brevity (for more details about the system model, please refer to [3,27,28]).
The equations of motion of the tray and of the actuators can be merged together and a second-order model of the N-DOFs system is obtained, with N = 2(N b + 1) + N A + 1 where , , ∈ ℝ N×N are the mass, damping and stiff- ness matrices of the overall system, respectively.The displacements, velocities and accelerations are collected in vectors z(t) ∈ ℝ N , ̇z(t) ∈ ℝ N and z(t) ∈ ℝ N , with collects the N A independent forces exerted by the actuators that are distributed along the system through the actuator influence matrix ∈ ℝ N×N A .As it will be shown in the section "System Description and Performance Specifications", B is a dense matrix, i.e., it has more than N A not-null entries.The partition of the coordinates into actuated and unactuated ones, to perform inverse dynamics, is not straightforward.In this work, coordinate transformation based on QR decomposition of B is performed.This approach leads to an orthonormal matrix (satisfying T = T = , where I is the identity matrix) and an upper triangular nonsingular matrix , with rank = rank( ) = N A , such that: = .By defining a new coordinate vector q = q A q U , and then by performing the pre multiplication of the lefthand side and of the right-hand side of Eq. (1) by T , then a (1) (2) z = q A q U , Fig. 1 Simplified scheme of the linear vibratory feeder assumed as the test case partitioned form of the equations of motion is obtained (for more details, refer to [27]) The coordinate vector is partitioned into actuated q A ∈ ℝ N A and unactuated q U ∈ ℝ N−N A coordinates, and the system submatrices are consistently defined: The dynamics of the actuated sub-system is defined through the first N A equations in Eq. ( 3), by showing how it is directly affected by the actuation forces f through the input matrix B A .In contrast, the unactuated coordinates are just indirectly controlled by the evolution of the actuated coordinates through the remaining N -N A equations Hereafter, the dependence with respect to time will be avoided, if not necessary, for brevity.Equation (4) defines constraints on the allowable state trajectories, and depends on both positions, speeds and accelerations.Thus, it is a set of second-order nonholonomic constraints [29], and any commanded trajectory should satisfy such constraints at any time.For this reason, the number of independent outputs to be imposed for the inverse dynamics is assumed to be equal to the number of independent actuation forces f.

The Algebraic-Differential Scheme of Model Inversion for Generic Periodic References
The inverse dynamics problem is aimed at computing the N A actuation forces collected in f , so that N A output coor- dinates, collected in the output vector , track the prescribed periodic desired displacement des (t) .In particular, the cho- sen coordinates should be defined to impose the desired tray motion, and therefore are defined as N A absolute tray displacements.A common requirement is to get uniform displacement along the tray, with a prescribed throw angle [27], thus imposing a relationship between the horizontal and the vertical displacements of the tray.
While the problem of dynamic inversion for undamped feeders tracking harmonic des has been solved by the authors in [27] through a purely algebraic scheme, the use of generic periodic (and hence non-harmonic) references and in the presence of damping is not straightforward, as briefly outlined in [16], and requires an algebraic-differential scheme.
Additionally, the possibility to impose just N A coordinates might severely affect the uniformity of the tray displacement, thus reducing the overall feeding capability of the device [3].
If z is adopted in the dynamic model, the definition of such output coordinates would be straightforward, since they are x t and two vertical tray displacements that belong to the set of coordinates y i (i = 1,…,N b + 1).Because of the coordinate partitioning introduced in Eq. ( 2), becomes a linear combination of the actuated-unactuated coordinates where ∈ ℝ N A ×N A and ∈ ℝ N A ×(N−N A) are two full-rank, constant matrices related to the choice of the three physical coordinates and to Q.It is possible to exploit Eq. ( 5) together with the desired output trajectory des to express the position q A , speed ̇qA and acceleration qA of the actuated coordinates as follows: Equation ( 6) allows the second-order nonholonomic constraint in Eq. ( 4) to be written as a function of q U , des and their derivatives, leading to the so-called "internal dynamics" It should be noted that the internal dynamics depends on the chosen output to be controlled, and this choice affects its poles because of the term −1 .

Stabilization of the Internal Dynamics
In the feeder under investigation, as it will be shown in the section "Application Example" through the numerical example, and in similar cases as well, the system in Eq. ( 7) is unstable because of poles with positive real parts.Hence, (5) the numerical integration of the internal dynamics leads to divergent solutions even when forced by bounded input ( des , ̇ des , ̈ des ).Stabilization should be therefore performed.An effective approach is the output redefinition, that is approximating y through a proper fictitious output ̃ so that Eq. ( 7) is approximated through a stabilized internal dynamics.̃ is, again, written as a linear combination of the actu- ated and unactuated coordinates, through proper matrices Provided that a proper choice of stabilizing ̃ and ̃ is done, the following stabilized internal dynamics is finally obtained: The choice of ̃ and ̃ is not trivial and should pursue a twofold goal: on the one hand, ̃ should ensure negative real parts of the poles of the approximated internal dynamics; on the other hand, it should lead to a negligible perturbation of the stable poles of the internal dynamics to limit the tracking error.
Forward numerical integration of the ODEs in Eq. ( 9) is adopted to compute qU , ̇qU and q U .The numerical integra- tion scheme employed, and the choice of the discretization time step, rely on the usual practices adopted in the numerical simulations of multibody systems (see, e.g., [30]).

Computation of the Actuation Forces
Once the unactuated coordinates accelerations, speeds and displacements are obtained, and will be henceforth denoted as q des U , ̇qdes U and qdes U , it is possible to compute the reference trajectory for the actuated coordinates q des A , ̇qdes A , qdes A by inverting the actual maps in Eq. ( 6), i.e., using the actual values of and ; indeed, the method just relies on algebraic computations in the subsequent steps.
Finally, the actuation forces collected in f are obtained through the solution of the first N A algebraic equations of the system in Eq. ( 3) that represent the invertible inverse dynamics of the actuated sub-system for imposed values of the motion of all the coordinates (8)

Aims of Dynamic Structural Modification
Although the presence of the second-order nonholonomic constraints in Eq. ( 4) imposes to assign the motion of just N A outputs, the proper functioning of long-tray feeders needs the correct assignment of more coordinates of interest.Indeed, a uniform motion of the conveyed parts is obtained if uniform oscillations are achieved all along the tray, not just in the (N A − 1) positions whose vertical displacements belong to .
It should be noted that also the use of feedback control is not effective in imposing the correct shape of vibration because of the limitation in assigning eigenvectors, and hence the vibrational mode shape, in underactuated systems [31,32].
To overcome these difficulties and to obtain a correct feeding, this paper exploits the concepts of dynamic structural modification (DSM), i.e., the modification of elastic and inertial physical parameters, to achieve the desired displacements for all the N d coordinates of interest (N A < N d < N).DSM is usually adopted for assigning poles, mode shapes and antiresonances in vibrating systems, with particular regard to undamped ones.In this work, it is exploited in a slightly different manner, by considering the system forced response in the presence of the desired output des .The application of DSM for feeders is significantly simplified, and improved in term of solvability and solution reliability, by neglecting damping and by considering harmonic motion [27].Under these circumstances, both the second-order nonholonomic constraint and the internal dynamics become algebraic equations (that define holonomic constraints relating the oscillation amplitudes of the unactuated coordinates), and the DSM problem can be solved as a purely algebraic optimization.
To apply this idea to the case under investigation, where the output reference des is an arbitrary periodic function, the Fourier series is exploited together with the system linearity, that allows splitting the DSM problem into more subproblems, each one addressing one of the N H dominant harmonic components of des .

3
Let des (t) be periodic with period T F = 2 ∕ F , then its jth entry, j = 1,…,N A , can be approximated as the following truncated series ( j,k and A j,k denotes the phase and the amplitude, respectively, of each harmonic) The phase of each harmonic is not of interest for the DSM, and therefore can be discarded.The continuous component can be discarded as well.
The selection of N H should provide an accurate signal representation, also considering the actuators bandwidths, while trading off with the need of reducing the size of the DSM problem.It should be further stressed that such an approximation of des is just used for the computation of the physical parameter modifications, to aid the inverse dynamics method proposed in the section "Inverse Dynamics" to get a quasiuniform displacement along the tray, despite the high value of the degree of underactuation (i.e., N − N A ). Once the system is modified through proper mass and stiffness modification, such a method will provide the almost-exact solution of the inverse dynamic problem; the sole approximation is the one required in stabilizing the internal dynamics.

DSM Problem Formulation for a Single Harmonic
In this section, the DSM problem is formulated and discussed for an arbitrary harmonic, indexed through k = 1,…, N H .The overall problem, considering all the dominant harmonics, is formulated in the section "DSM Problem Formulation for NH Harmonics".
The internal dynamics formulated in Eq. ( 7), under the hypotheses of neglecting damping and in the presence of harmonic reference, is written as follows: where des 0,k collects the amplitudes A j,k j = 1,…,N A , for the kth harmonic.Hereafter, subscript 0 denotes the amplitude of the signals.The inverse of the receptance matr ices , are introduced and defined as (11) des Hence, Eq. ( 12) is therefore recast in the following more compact form: Equation ( 14) can be exploited to modify the system internal dynamics, when forced by the desired output, such that all the N d coordinates of interest track the desired reference.The following structural modification matrices are defined to represent the model changes due to the N p design variables collected in vector p: ΔK AU , ΔM AU , ΔK UU , ΔM UU are the mass and stiffness modification matrices obtained as The mass and stiffness modification matrices of the model with physical coordinates z, are denoted by ΔK, ΔM.The topology of the modification matrices is a-priori assumed and it represents the admissible choices in the system design, as usually done in the state of the art (see, e.g., [3,27]).
Once structural modifications are considered, Eq. ( 14) is recast, in terms of the modified system, as follows: with:

DSM Problem Formulation for N H Harmonics
Since the system is linear, the superposition principle can be applied.Therefore, for all the N H harmonics adopted to model the output reference signal in the frequency domain, the structural modification problem introduced in Eq. ( 17) becomes ( 14) Let us consider the scenario, introduced in the section "Aims of Dynamic Structural Modification", where a set z des s 0,k of N d coordinates of the system are of interest (k = 1,…,N H ). The remaining N − N d are unassigned coordinates, whose amplitudes are collected in vector z f0,k , k = 1,…,N H , as additional unknowns in the IDSM problem, together with the design parameters p.The desired displacement amplitude vector in the transformed coordinate q is obtained through the transformation matrix Q, thus making q des A0,k and q des U0,k functions of the unknown where superscript des denotes the desired displacements.
The N p design variables are constrained to belong to a feasible domain defined through lower and upper bounds as follows p L ≤ p ≤ p U (with element-wise inequalities), and p L , p U ∈ ℝ N p are related to technological and economical constraints on the admissible system modifications.Constraints on the unassigned displacements are introduced as well to define their admissible values through the following element-wise inequalities, i.e., Such constraints provide several benefits in the development of an effective solution method for finding the optimal values of p, as explained in the section "Solution of the Minimization Problem".All the constraints on the unknowns are merged to define the feasible domain .
The IDSM problem is defined by recasting Eq. ( 19) into the following minimization: where ∈ ℝ N H ×N H is a diagonal matrix that collects N H sca- lar coefficients that are introduced to reflect different levels of concern on each harmonic adopted to approximate the output reference signal.The cost function also includes a Tikhonov regularization term ∈ ℝ N p ×N p , i.e., a positive- definite diagonal matrix that collects N p scalar coefficients, introduced to improve numerical conditioning, and it is adopted to properly weigh the structural modification parameters [33].
Equation ( 21) features bilinear terms due to the products between the two unknowns, and hence, it yields to a non-linear and non-convex minimization problem.Handling non-convex optimization is cumbersome, and (20) , the main threat is related to the risk of attaining a local minimum instead of the desired global optimal solution.Hence, it requires a proper solving strategy, such as the one explained in the section "Solution of the Minimization Problem", to boost the attainment of the global optimum.Once the optimal values p opt for the design variables are found, it is possible to compute the modified system matrices as follows: Then, the actuation forces are computed through the procedure proposed in the section "Inverse Dynamics", by adopting the modified system matrices obtained by solving the DSM problem in Eq. (21).

Solution of the Minimization Problem
The minimization problem in Eq. ( 21) is solved through homotopy optimization.This technique has been widely adopted to solve structural modification problems [26,27,34].Homotopy relies on the idea of solving a finite set of optimization problems f i , starting from a convex relaxation f c of the non-convex function to be solved and moving to the non-convex function f nc to be solved by means of a path, called the homotopy map defined as where ∈ [0;1] is a morphing parameter that varies dis- cretely from 0 to 1. Since is monotonically increasing, f i varies from convexity to non-convexity.
In the last optimization step ( = 1 ), the func- tion f nc is solved and the solutions obtained at the last step (i.e., opt = p opt = p =1 and opt = p opt = p =1 ) are the optimal structural modification matrices.It is worth to notice that homotopy procedures prescribe to exploit the solution of the previous optimization step as the initial guess of the following minimization.This feature is fundamental to boost the attainment of a global minima instead of a local one.
More details regarding homotopy optimization and its applications to DSM problems are provided in [27,34] and here omitted for brevity.

Outline of the Method
The method proposed in this paper is outlined in Fig. 2, where the required equations are referenced.The method requires, as the input, the system dynamic model through matrices M, C, K, B, together with the required motion profile for the output coordinates des (t) ; the method output are the actuation forces f (t).
Two decisions may be needed: first, once the desired outputs are defined, if the dynamic system features unstable internal dynamics, it is necessary to stabilize it as shown in the flowchart in Fig. 2; otherwise, stabilization is not required.For these systems, it is quite common experiencing a non-minimum phase behavior, and hence, stabilization is usually likely to be required.Second, if the tracking errors obtained in the coordinates not belonging to des (and hence not imposed in the desired output) are larger than desired, it is possible to take advantage of the DSM modification algorithm to improve trajectory tracking in more than N A coordinates.

System Description and Performance Specifications
Let us consider the linear vibratory feeder sketched in Fig. 1.It recalls a typical linear feeder adopted in packaging lines.The system properties and the performance requirements of the proposed test case are based on typical specifications set for proper operating these devices in industrial plants [3,35].
The number of actuators is N A = 3 and the force distribution matrix is  1.
To ensure a uniform flow of the conveyed parts, it is required that the tray behaves as a rigid body with uniform displacements ensuring a uniform throw angle α f = 20°.This requirement would impose assigning the 5 vertical displacements y 1 , y 2 , y 3 , y 4, y 5 and the horizontal one x t , and the nodal elastic rotations as well.All the vertical displacements should track the same references, while the elastic nodal rotations of the tray should be equal to zero, to ensure the rigid-like motion.These specifications yield to N d = 11 coordinates of interest.
However, since N A = 3, just three outputs can be assigned.To boost the achievement of the desired goal, the output vector is defined as = y 2 y 4 x t T : two vertical displacements in the attachment points of the first and third actuators, and the horizontal displacement to set the desired throw angle.des (t) is defined as a periodic asymmetric reference profile, whose period is T F = 28.60 ms, i.e., ω F = 2π 35 rads −1 , as often adopted in feeding of small parts [3,27].
The peak oscillation amplitude in the throw angle direction is s d max = 5 mm , and hence, y d max = s d max sin − f and The required motion profiles for y 2 , y 4 , x t are shown in Fig. 3a.

Application of the Inverse Dynamics on the Original System
The inverse dynamic approach proposed in the section "Inverse Dynamics" is applied to the system with the original parameters, to meet the motion requirements listed in the section "System Description and Performance Specifications".The linear vibratory feeder, for the chosen output coordinates , has unstable internal dynamics, as shown in Fig. 4a, where a right-half complex plane pole appears.A  proper selection of the tuning matrices ̃ and ̃ enables to stabilize its internal dynamics as evidenced in Fig. 4a; further, the original internal dynamics poles are only slightly perturbed in the approximate internal dynamics.The actuation forces are computed through the method outlined in the section "Outline of the Method" without applying structural modification, and the time history of the obtained actuation forces is shown in Fig. 5.The method implementation is not demanding under a computational point of view: the average computational time employed to compute the actuation forces for an interval of 0.3 s, by means of MATLAB/Simulink running on an i7 16 GB RAM laptop, is equal to 1.7 s.It is worth to notice, however, that open-loop control can be computed offline.
The multibody system is simulated with T s = 0.01 ms, adopting the Runge-Kutta ODE4 solver for the numerical integration, and the time histories of the desired output coordinates y 2 , y 4 , x t are reported in Fig. 6, where it can be noticed that the reference signals and the simulated signals are almost overlapped.This aspect is confirmed by the tracking errors shown in Fig. 6, where e i = y des i − y i is the tracking error related to coordinate y i .Similarly, e x is the tracking error related to x t .Figure 6 highlights the tracking errors  Figure 10a shows the overall RMS error for the imposed coordinates, denoted as e RMS des = RMS des − .Its time his- tory evidences that the inverse dynamic method proposed in this paper is effective and, indeed, it enables to achieve low tracking error values: the maximum value of e RMS des is equal to 0.16 mm.
This paper pursues a further goal that is to achieve a rigid-like motion of the entire tray, as already pointed out in the previous section.This aspect is analyzed by considering the remaining vertical displacements (y 1 , y 3 , y 5 ) and comparing them with the reference signals for the desired vertical output displacements.The simulated time histories, together with the tracking errors, for these signals are reported in Fig. 7. Higher tracking errors with respect to the imposed output coordinates are obtained.This yields to an unsatisfying behavior for the feeder tray.These results are summarized by the following performance indexes: Additionally, the nodal rotations at each node of the tray are reported in Fig. 9, together with the desired nodal rotations that should be null to enforce uniform displacements all along the tray of the feeder.In this case The behavior of the non-imposed coordinates is evaluated through two additional error vectors adopted as synthetizing performance indexes.The RMS error of the remaining vertical displacements: e RMS T is also considered, which contains all the nodal rotations, allowing for a concise evaluation of the rotations at each node of the tray.Its RMS value is reported in Fig. 10c and the maximum value of e RMS is equal to 0.13°.An overall performance index is adopted to evaluate the tracking error between the desired and the obtained vertical displacements of the tray: it is computed as The time history of the RMS value of e y t is shown in Fig. 10d and the maximum value of e RMS y t is equal to 1.6 mm.The results here discussed highlight that the error in the trajectory tracking of the imposed coordinates approaches zero.Conversely, for the remaining non-imposed coordinates, the tracking errors are larger and these should be reduced to achieve a more uniform and rigid-like behavior of the tray.In this light and following the outline of the proposed method provided in Fig. 2, it is possible to exploit structural modification to improve the performances of the system.

Application of the Dynamic Structural Modification
The required motion profile is decomposed through Fourier Series and its dominant harmonics are: 35, 70 and 105 Hz, as evidenced in Fig. 3b.These harmonics are adopted to define the DSM problem in Eq. ( 21) with = diag 10 5 ;10 3 ;1 , which reflects the relative importance of the three harmonics.
The following system parameters are assumed to be modifiable: the actuator masses and the leaf spring stiffnesses, respectively, m ai and k ai , i = 1,…3, the tray support stiffnesses in the vertical (k l and k r ) and horizontal directions (k x ).In addition, 5 nodal lumped masses m i , i = 1,…,5, to be placed to the tray, are introduced in the DSM problem.Modifications of the tray flexural stiffness EI and linear mass density ρA are not allowed, although the method can handle their modifications too.Indeed, for an existing feeder, it is very difficult to modify its tray.Conversely, including these parameters would allow for better results.The Tikhonov regularization matrix adopted is = diag 1;1;1;1;1;10 −3 ;10 −3 ;10 −3 ;5 ⋅ 10 −3 ;5 ⋅ 10 −3 ;5 ⋅ 10 −3 ;10 3 ;10 3 ;10 3 , and these weights are chosen to penalize modifications of the lumped masses and of the leaf springs of the actuators.The optimal values of the modified variables together with the structural modification constraints adopted are listed in Table 1.The average computational time needed to compute the optimal structural modifications in MATLAB running on a laptop with an i7 16 GB RAM is equal to 5.2 s.

Inverse Dynamics on the Modified System
The modified parameters reported in Table 1 are adopted to update the mass and stiffness matrices accordingly to Eq. (22).The same motion specifications listed in the section "System Description and Performance Specifications" are required.For the chosen output coordinates , the internal dynamics is unstable once again, as shown through the pole map in Fig. 4b, where it is evident the presence of a righthalf complex plane pole.Once ̃ and ̃ are chosen, the internal dynamics is stable, as shown in Fig. 4b, indeed the right-half complex plane pole migrates beyond the imaginary axis and its real part becomes negative; further, only a negligible spillover is obtained on the remaining stable poles.
The inverse dynamic approach proposed in the section "Inverse Dynamics" is exploited once again for the modified system, to compute the actuation forces through the procedure outlined in the section "Outline of the Method".The time history of the obtained actuation forces is shown in Fig. 5 and compared with the actuation forces computed for the original system.It shows that the adoption of DSM does not increase the average RMS value of the actuation forces, that is 1.15e5 N for the original system, and 1.14e5 N for the one modified with the optimal structural modification.While f 1 and f 3 decrease (both in term of RMS and maximum values), f 2 increases.Indeed, the actuation forces of the original system are Those of the modified system are It should be clearly stated that there is no a-priori control on the changes in the actuation forces, that can either increase or decrease, since DSM is focused on just reducing the tracking error.
The same sample time and numerical integrator adopted for the original system (see the section "Application of the Inverse Dynamics on the Original System") are used for the modified system to provide a fair comparison.The application of the actuation forces to the dynamic model of the modified vibratory feeder yields to the response of the imposed coordinates y 2 , y 4 , x t shown in Fig. 6.Once again, the reference signals and the simulated ones are almost overlapped, as corroborated by the small values of the tracking errors reported in Fig. 6.Hence, the proposed inverse dynamics approach is effective.Additionally, the RMS and maximum values of the tracking errors for the imposed coordinates are the following: These results evidence that the adoption of DSM does not alter the trajectory tracking performances for the desired output coordinates.Indeed, the tracking error is smaller for y 2 and x t, in the modified system and only an increase of the RMS and maximum values of the tracking error is obtained for the output coordinate y 4 .Table 2 summarizes the obtained errors together with the percentage variation of these parameters with respect to those obtained for the original system.
The time history of the overall RMS error for the desired output coordinates is reported in Fig. 10a.It confirms the effectiveness of the inverse dynamics method proposed in this paper and it shows, once again, that the mechanical redesign of the feeder is almost negligible for The increase of performances provided by the adoption of structural modifications is shown in Fig. 7, where the trajectory tracking responses for the non-imposed vertical displacements (y 1 , y 3 , y 5 ) are reported together with their tracking errors.An improvement of the motion of the entire tray is obtained with the goal of achieving the desired rigid-like motion.This result is confirmed by the following values for the RMS and maximum tracking errors: • e RMS The percentage variations with respect to the original system are summarized in Table 2 and confirm the benefits obtained through the adoption of the DSM.
The Cartesian displacements of each tray node for the original and modified system are reported in Fig. 8, together with the desired reference signals.It can be noticed that the inverse dynamic algorithm is effective in obtaining the desired throw angle for the imposed coordinates, i.e., y 2 , y 4 and x t .Further, structural modification allows to significantly improve the throw angle for the non-imposed nodes, i.e., those denoted through y 1 , y 3 and y 5 .The effectiveness of the proposed method is assessed by the elastic rotations at each node of the tray, which are shown in Fig. 9 and take the following values for the modified system: The time history of e RMS y ni is reported in Fig. 10b: it clearly shows that the tracking error for the non-imposed coordinates is remarkably reduced for the modified system.Indeed, it is characterized by a maximum value equal to 1.00 mm which yields to a percentage reduction of the 50% with respect to the original system.The tracking error related to the nodal elastic rotations, whose reference is imposed to be null to enforce the uniformity requirement for the tray displacements, is shown in Fig. 10c.The maximum value of e RMS is equal to 0.07°, i.e., it is 45.2% smaller if compared with the original system.The time history of the overall performance index e RMS y t is reported in Fig. 10d and it shows that the tracking error for the vertical coordinates of the tray remarkably decreases; indeed, its maximum value is equal to 0.81 mm leading to a percentage reduction equal to 49.1% with respect to the original system.
These results, together with those summarized in Table 2, evidence that the concurrent usage of the proposed inverse dynamics approach together with the optimized mechanical design obtained through the DSM algorithm is effective.Indeed, both the imposed output coordinates and the remaining non-imposed coordinates of interest lead to a uniform rigid-like behavior of the tray.

Conclusions
This paper proposes an integrated method for ensuring uniform motion of the tray of linear vibratory feeders, operated in open-loop control with periodic references, by means of the integrated use of dynamics and dynamic structural modification.The goal is challenging, since these systems are underactuated and non-minimum phase.Additionally, although inverse dynamics enables to impose just a number of desired outputs equal to the number of independent actuation forces, specifications on the motion of the whole tray is usually set to ensure proper motion.This problem cannot be therefore solved if internal dynamics is used alone.To overcome this limitation, this paper exploits DSM.
The inverse dynamics technique exploits the QR decomposition to transform the model of the system into a representation featuring actuated and unactuated coordinates.Then, the output reference coordinates are defined, and the system internal dynamics is analyzed.Since unstable internal dynamics occurs for non-minimum phase systems, the output redefinition method is exploited to stabilize it.The stabilized internal dynamics ordinary differential equations are then integrated to obtain the reference values for the unactuated coordinates and, finally through algebraic calculations, the actuation forces.The proposed approach leads to a causal, almost-exact solution of the inverse dynamics problem, thanks to the use of the actual output map in the computation of the forces.
The formulation of the internal dynamics with the actual output map is then exploited to perform structural modification, by considering the undamped system in the presence of harmonic excitations, thus leading to a straightforward problem formulation.Hence, the periodic output reference is decomposed through the Fourier Series truncated to a subset of harmonics of interest.Then, a DSM problem for the internal dynamics is formulated to improve the system performances by properly shaping such a constraint on the allowable motion.The solvability of the DSM problem for all the choices of the design parameters is ensured by recasting it into a non-linear non-convex optimization solved through the homotopy optimization technique.
To summarize, the main novel contributions of this integrated method for optimal control and design of feeders are: • Generic periodic references are handled both in the inverse dynamics and in the DSM; • The presence of viscous damping is considered in the inverse dynamics problem; • The second-order nonholonomic constraints, representing the motion of the unactuated coordinates leads to unstable internal dynamics, that is coped with in this paper through an output redefinition strategy; • DSM is performed by considering the internal dynamics, to better represent the effect of the desired reference on the desired output.
The effectiveness of the proposed method is assessed by its application to the challenging test case of a 14-DOFs industrial linear vibratory feeder employed in packaging plants.The system is excited through 3 independent actuators and its tray, over which products flow, is modeled through 11 coordinates.To ensure a uniform flow of the transported parts over the tray, it is required that the tray behaves as a rigid-like beam even if it is flexible.Underactuation exacerbates the problem, since only 3 tray coordinates can be imposed through the inverse dynamics algorithm, while an effective achievement of the remaining 8 1 3 coordinates of the tray relies on the DSM algorithm.The results proposed in the paper highlight the effectiveness of the proposed method.Indeed, the proposed inverse dynamics technique enables to obtain low tracking error for the imposed coordinates.Further, the optimized mechanical design obtained through structural modification enables to drastically reduce the tracking error also for the non-imposed coordinates of interest.

Fig. 2
Fig. 2 Outline of the proposed method

Fig. 3 Fig. 4
Fig. 3 Output reference for y 2 , y 4 , x t : time history (a) and their frequency spectrum (b)

Fig. 5 6
Fig. 5 Actuation forces f 1 (a), f 2 (b) and f 3 (c) for the original and the modified systems

=
RMS y des ni − y ni , where y des ni collects the desired motions of the not imposed vertical displacements, which are equal to the displacements required for y 2 and y 4 , and y ni = y 1 y 3 y 5 T .The time history of e RMS y ni is shown in Fig. 10b and its maximum value is equal to 2 .0 m m .A d d i t i o n a l l y, t h e e r r o r v e c t o r e =

Fig. 8
Fig. 8 Cartesian displacements of the tray nodes for the original and the modified systems

Fig.Fig. 10
Fig. Elastic rotations for the tray nodes for the original and the modified systems

Funding
Open access funding provided by Università degli Studi di Padova within the CRUI-CARE Agreement.This research was funded by the Italian Ministry of University and Research by the research Grant "SISTEMA-Dipartimenti di Eccellenza".

Table 1
Original and modified system parameters

Table 2
Comparison of the tracking errors for the original and the modified systems des is equal to 1.94 × 10 −4 m.