Pole Assignment for Active Vibration Control of Linear Vibrating Systems through Linear Matrix Inequalities

This paper proposes a novel method for pole placement in linear vibrating systems through state feedback and rank-one control. Rather than assigning all the poles to the desired locations of the complex plane, the proposed method exactly assigns just the dominant poles, while the remaining ones are free to assume arbitrary positions within a pre-specified region in the complex plane. Therefore, the method can be referred to as “regional pole placement”. A two-stage approach is proposed to accomplish both the tasks. In the first stage, the subset of dominant poles is assigned to exact locations by exploiting the receptance method, formulated for either symmetric or asymmetric systems. Then, in the second stage, a first-order model formulated with a reduced state, together with the theory of Linear Matrix Inequalities, are exploited to cluster the subset of the unassigned poles into some stable regions of the complex plane while keeping unchanged the poles assigned in the first stage. The additional degrees of freedom in the choice of the gains, i.e., the non-uniqueness of the solution, is exploited through a semidefinite programming problem to reduce the control gains. The method is validated by means of four meaningful and challenging test-cases, also borrowed from the literature. The results are also compared with those of classic partial pole placement, to show the benefits and the effectiveness of the proposed approach.


State of the Art
The assignment of the dynamic response of vibrating mechanical systems, such as structures, mechanisms, or multibody systems, is often performed by properly assigning the poles of the controlled systems. Both active feedback control [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] and passive approaches (i.e., parameter modifications, see e.g., [19][20][21]) or are exploited to accomplish this important task. Indeed, the system poles, which are often denoted as the eigenvalues, define the system stability as well as the properties of the transient response such as the damping ratio, the rise time, and the settling time. Therefore, pole placement (also denoted as "pole assignment" or "eigenvalue assignment") is a broadly adopted vibration control technique, due to the simple definition of the aforementioned specifications.
The goal of pole placement through active control is to find the gains, usually exploiting state feedback, to ensure that the closed-loop poles are exactly the desired ones, to feature the desired dynamic performances. The literature proposes several computational methods to compute the gains,

Motivations and Contributions of This Paper
In the present work, a different and new method for RPP through state-feedback is proposed and validated, to solve the issues related to the features of vibrating systems, as previously discussed. The goal is to compute the control gains that exactly assign the set of dominant poles to the prescribed locations of the complex plane, while the remaining non-dominant poles are clustered into LMI regions to feature some dynamic properties. The method relies on two-stages. In the first stage, the receptance method is adopted to exactly place the dominant poles and to define the non-spillover condition on such poles. In the second stage, by using a first-order model on a reduced basis, the gains are modified to cluster the unassigned poles into the desired LMIs without spillover on the previously assigned dominant poles. The degrees of freedom in the assignment of the poles within the LMI is, indeed, useful to reduce the control effort. Hence, a norm minimization problem is formulated to solve the assignment problem in the second stage. The proposed method, therefore, is aimed at improving performances achievable by complete pole placement and by PPP too, by proposing a new mathematical framework that effectively solves a wide range of pole assignment tasks and exploits the presence of fewer specifications on the pole assignment. Four meaningful and challenging numerical test-cases, also borrowed from the recent literature, are exploited to validate the method that handles also asymmetric and unstable or marginally stable systems.

Definitions
Let us consider a linear time-invariant vibrating system, modelled through ordinary differential equations: M ..
where M, C, K∈R N × N are the mass, damping, stiffness matrices, respectively, and b∈R N × Nb is the input matrix of the control forces (often named the control force distribution matrix), with N being the number of degrees of freedom and Nb the number of independent control forces. q is the vector of the displacement of the N independent coordinates. It is here assumed that the system is controllable, and that full state feedback is available. The issue related to state estimation are here neglected, and it is assumed that a wise tuning of the observer does not cause spillover on the assigned poles [17].

Control Specifications
The method is developed for the case of rank-one control, i.e., Nb = 1. Control force is defined through state feedback as: where f, g ∈ R N are the control gain vectors. The reference is omitted since it creates a feedthrough term that does not alter the closed-loop poles. The problem to be solved in this paper can be summarized as follows. Problem: Given the self-conjugate set of desired poles, Σ d = µ 1 , µ 2 , . . . , µ p , with p ≤ 2 N, compute the controller gains f and g such that: 1.
the controlled system has exactly the p desired poles in µ 1 , µ 2 , . . . , µ p , 2. the remaining 2N − p poles must belong to some prescribed regions of the complex plane.
The novel method proposed in this paper is performed through two steps. In the first one, the gains that solve requirement 1 are computed, and the receptance method can be exploited. In the second one, the gains are modified to meet requirement 2.

Formulation for Symmetric Systems
The method is first introduced with reference to a system with symmetric matrices. To understand the genesis of the method, let us consider the case of complete pole placement, i.e., p = 2N as proposed in [1]. The gains ensuring exact placement of all the poles are those that solves the following linear system where it is defined that with H being the receptance matrix of the open-loop system, Appl. Sci. 2020, 10, 5494

of 21
If p = 2N, the coefficient matrix of the linear system in Equation (3), henceforth denoted as T is a square matrix. Whenever the system is controllable, i.e., rank λ 2 i M + λ i C + K b = N for all the open-loop poles λ i (i = 1, . . . , 2N) and µ 1 , . . . , µ 2N ∈ C is self-conjugate, matrix T is invertible, hence the solution of Equation (3) exists, it is real and unique. A linear system can be formulated in the case of higher-rank control too (see e.g., [18]).
In contrast, if p < 2N then matrix T is rectangular and Equation (3) is an underdetermined linear system. Hence, infinite solutions can be found. The parametrization of the set of solutions is helpful to perform RPP, as required in the second stage, while accomplishing a secondary task, such as reducing the controller gains. Without loss of generality, it can be assumed that rank(T) = p. Hence, the solutions can be written as where k 0 = f 0 g 0 ∈ R 2N is a particular solution of Equation (7), V is a 2N × (2N − p) matrix whose columns span ker(T) and k r ∈ R 2N−p is a vector of parameters that are to be properly determined.
The particular solution k 0 can be calculated either through a least-squares problem, to compute the minimum norm solution, or by imposing arbitrary values to 2N − p entries of the gain vector (such as 0, for example). Additionally, 2N − p rows can be wisely added to T. If a least-squares solution is sought, the analytic solution is although numerical solutions of the minimum-norm problem can be exploited too.
Regardless of how f 0 g 0 is computed, and for arbitrary choices of k r ∈ R 2N−p , the gains f, g solving Equation (7) lead to a controlled system that has the poles µ 1 , . . . , µ p .

Extension to Asymmetric Systems
Even though the linear system in Equation (3) has been introduced in [1] for symmetric systems, the receptance method has been extended to some asymmetric systems with non-symmetric damping and stiffness matrices in [2] and [4]. In these cases, the linear system in Equation (3) should be slightly adjusted. First of all, matrices C and K should be decomposed into the sum of the symmetric parts C s , K s ∈ R N×N and asymmetric ones C a , K a ∈ R N×N , namely: Even if such a decomposition is not unique, the result is not affected by such a choice. Then, it is introduced H S (λ) = λ 2 M + λC s + K s This formulation of H(λ) can be exploited to obtain vector t i as in Equation (4) and hence the linear system as in Equation (3). More details can be found in [4,6,19].

Overview
The aim of this section is to formulate the condition ensuring that the remaining 2N − p poles belong to certain regions of the complex plane in a way that is computationally reliable, and then to explain how to properly choose k r that fulfils such a requirement. First, the first-order representation of the controlled system is employed: where Then, the system is transformed into a reduced one representing the 2N − p poles whose location is not exactly imposed. Such a reduction introduces the state of the reduced system,x ∈ R 2N−p , by exploiting the following transformation where V has been defined in Equation (7). Let U ∈ R (2N−p)×2N be the pseudoinverse of V. Since the columns of V are linearly independent and U is its left inverse, it follows that UV = I 2N−p . Hence the reduced subsystem is defined as follows: Since it has been defined k = k 0 + Vk r (see Equation (7) and (13)) can be further manipulated as follows: By definingÂ = U(A − Bk 0 )V,B = UB andk = V T Vk r , for brevity of representation, the reduced system can be finally written as

Linear Matrix Inequalities
Once the reduced system is written in a state-space form, LMIs can be exploited. A brief explanation of the basic theory of LMI is provided in the present section. More details can be found in the literature, such as in book [22].
Let D be a subset of the complex plane. A first-order system . x = Ax is D-stable if all its poles lie in D. There is a class of regions that have a straightforward characterization of D-stability: the so-called LMI regions. D ⊆ C is denoted as a LMI region if there exist a symmetric matrix R ∈ R d×d and a matrix S ∈ R d×d (for a suitable size d depending on the shape of the LMI region) such that Appl. Sci. 2020, 10, 5494 6 of 21 where f D (z) = R + zS + zS T , and ≺ 0 means negative definite. LMI regions are convex and are symmetric respect to the real axis. For example, conic sectors, vertical half-planes, horizontal strips, ellipses, parabolas and hyperbolic sectors are LMI regions, as well as any intersection of the above. For instance, the left half-plane is an LMI region obtained by setting R = 0 and S = 1. Some samples of LMI regions, and their relations with the properties of the system dynamic response, are shown in Figure 1.
lie in  . There is a class of regions that have a straightforward characterization of  -stability: the : . Some samples of LMI regions, and their relations with the properties of the system dynamic response, are shown in Figure 1. The following theorem, whose proof is provided in [13], characterizes  -stability for LMI regions.
The symbol ⊗ is the Kronecker product and , R S are the matrices used to define the LMI.
Theorem 1 generalizes the Lyapunov stability, that is a special case of Equation (18) (by setting 0 = R and 1 = S ).

Control Gain Synthesis
In accordance with the LMI theory proposed, k ensures that the closed-loop poles of the system in Equation (11) It should be noted that such a matrix inequality becomes non-linear. Indeed, both k and X are unknown and hence Equation (19) has bilinear terms of the unknowns. Linearity is restored, defining the auxiliary variables ˆT = X l k , thus leading to the following LMI: The following theorem, whose proof is provided in [13], characterizes D-stability for LMI regions.
Theorem 1: Let D be an LMI region. The system The symbol ⊗ is the Kronecker product and R, S are the matrices used to define the LMI.
Theorem 1 generalizes the Lyapunov stability, that is a special case of Equation (18) (by setting R = 0 and S = 1).

Control Gain Synthesis
In accordance with the LMI theory proposed,k ensures that the closed-loop poles of the system in Equation (11) lie in the prescribed LMI region D if Equation (18) holds forÂ CL =Â −Bk T : It should be noted that such a matrix inequality becomes non-linear. Indeed, bothk and X are unknown and hence Equation (19) has bilinear terms of the unknowns. Linearity is restored, defining the auxiliary variables l =k T X, thus leading to the following LMI: The LMI can be exploited as a constraint in the control gain synthesis, by defining a suitable cost function to be accounted for in a semidefinite programming (SDP) problem. It is here proposed to solve the following problem, since it is convex and hence global optimal solution can be effectively found: The LMI can be exploited as a constraint in the control gain synthesis, by defining a suitable cost function to be accounted for in a semidefinite programming (SDP) problem. It is here proposed to solve the following problem, since it is convex and hence global optimal solution can be effectively found: Once l and X are obtained, k is computed from the relation l . Finally, f and g are computed from Equation (7). Under mild hypotheses, which will be discussed in the Section 4.4, the closed-loop system with such gains (i.e., , and the other ones belong to  , i.e., Once l and X are obtained,k is computed from the relationk T = lX −1 . Finally, f and g are computed from Equation (7). Under mild hypotheses, which will be discussed in the Section 4.4, the closed-loop system with such gains (i.e., A CL = A − Bk) features a set of poles Σ which includes Σ d = µ 1 , . . . , µ p , and the other ones belong to D, i.e., Σ\Σ d ⊂ D. Additionally, it should be noticed that and, since X is a positive definite matrix, minimizing the norm of l guarantees that the controller effort is kept small.

Insights on the Reduced Model and on the Inclusion Principle
In order to justify the procedure employed for gain computation, a discussion on the Inclusion Principle [23] is proposed for the case under investigation, i.e., on the fact that a "part" of the overall system behavior is reproduced by its contraction having smaller dimension. It should be noted that the formulation of the first-order model through a reduced state is beneficial to reduce the size of the matrix adopted to define the LMI. Hence, the numerical computation of the gains benefits of this system scale reduction. This is another important and novel feature of the proposed method.
The following notation is introduced: let us denote V a matrix whose columns span ker(U), and U the unique left inverse of V such that ker U = im(V). Let us also define W = V V , which is clearly invertible. We want to prove that: First, it must be noted that VU + V U = I (the proof is tedious but straightforward, so it will be omitted). Hence, The claim is a consequence of the fact that UA CL V =Â CL . The result of Equation (23) can be further simplified, by using the well-known fact that the null space of the transpose of a matrix is equal to the null space of the pseudoinverse. In our case, ker V T = ker(U), and then V T V = 0. Therefore, Provided that: it is demonstrated that the set of closed-loop poles is constituted by the poles ofÂ CL (which are constrained into the LMI region D) and the poles of U A − Bk T 0 V, which are µ 1 , . . . , µ p , as can be seen from the numerical results obtained in this work. Although a proof that Equation (26) always holds is not provided, it has been experienced in all the test cases studied (also including dozens of examples not shown in the paper, with random requirements and parameters) that such a condition is satisfied.

Results
In this section, the proposed method is applied and validated by means of four test-cases with increasing complexity, as summarized in Table 1. In the first stage of the method k 0 is solved by means of Equation (8), i.e., computing the minimum Euclidean norm solution. In the second stage, the SDP programming is solved exploiting the Yalmip toolbox [24] for Matlab and the MOSEK solver. The negative definiteness of the LMI and the positive definiteness conditions for X have been approximated as follows to improve the numerical solvability of the SDP problem: is satisfied.

Results
In this section, the proposed method is applied and validated by means of four test-cases with increasing complexity, as summarized in Table 1. In the first stage of the method 0 k is solved by means of Equation (8), i.e., computing the minimum Euclidean norm solution. In the second stage, the SDP programming is solved exploiting the Yalmip toolbox [24] for Matlab and the MOSEK solver. The negative definiteness of the LMI and the positive definiteness conditions for X have been approximated as follows to improve the numerical solvability of the SDP problem: where I denotes identity matrices of suitable dimensions and i ε is a scalar approaching zero. All the tests have been solved by setting

Test-Case 1: Three-Mass System
Let us consider a simple three-mass system modeled as follows: This simple system is often used in the literature to provide simplified models of more complex structures, to represent their low-frequency dominant dynamics (see e.g., [25]).
The assignment task consists of assigning two dominant pair of complex conjugate poles at μ = − ± 1,2 0.001 1.500i and μ = − ± 3,4 0.001 3.000i , while the remaining closed-loop poles i μ must (27) where I denotes identity matrices of suitable dimensions and ε i is a scalar approaching zero. All the tests have been solved by setting ε i = 1 · 10 −8 .

Test-Case 1: Three-Mass System
Let us consider a simple three-mass system modeled as follows: This simple system is often used in the literature to provide simplified models of more complex structures, to represent their low-frequency dominant dynamics (see e.g., [25]).
The assignment task consists of assigning two dominant pair of complex conjugate poles at µ 1,2 = −0.001 ± 1.500i and µ 3,4 = −0.001 ± 3.000i, while the remaining closed-loop poles µ i must belong to the intersection of the LMI regions D 1 , D 2 and D 3 ensuring that Re(µ i ) ≤ −0.10, ξ(µ i ) ≥ 0.02 and max ω np,i ≤ 5, where ω np,i denotes the i-th pole natural frequency and ξ(µ i ) its damping ratio. Figure 2 shows the resulting open-loop and closed-loop poles in the complex plane. The pair of closed-loop poles computed by means of k 0 is exactly assigned, while the unassigned ones are unstable and do not feature the desired dynamic properties. Conversely, the modified control gain vector k stabilizes the system and clusters the unassigned poles within the prescribed LMI region. The open-loop and closed-loop poles are also summarized in Table 2, while the control gains are reported in Table 3.
closed-loop poles computed by means of 0 k is exactly assigned, while the unassigned ones are unstable and do not feature the desired dynamic properties. Conversely, the modified control gain vector k stabilizes the system and clusters the unassigned poles within the prescribed LMI region. The open-loop and closed-loop poles are also summarized in Table 2, while the control gains are reported in Table 3.

Test-Case 2: Slider-Belt System
The second test-case exploits the asymmetric model of a slider-belt system with friction, proposed in [2,6,19,26]. The system is sketched in Figure 3. This system has marginally stable, or even Im( )

Test-Case 2: Slider-Belt System
The second test-case exploits the asymmetric model of a slider-belt system with friction, proposed in [2,6,19,26]. The system is sketched in Figure 3. This system has marginally stable, or even unstable, poles depending on the friction coefficient, here denoted as η. For this reason, this test-case is here considered to show that the proposed method can handle both the instability and the asymmetry issues.
The model independent coordinates are the translations of the three masses: q = [x 1 , y 3 , x 2 , y 2 ] T . The asymmetric stiffness matrix K is the sum of two parts: the symmetric part K s and the asymmetric part K a . The latter is caused by the pre-compression normal force and the friction force that acts at the slider-belt interface. The resulting system matrices are reported in Equation (29): The model parameters have been assumed as in [19]: The value of the friction coefficient is the one at the flutter instability, i.e., η = 0.3868. Indeed, a pair of complex conjugate poles with real part appears: 0.00 ± 8.73i [19].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 22 unstable, poles depending on the friction coefficient, here denoted as η. For this reason, this test-case is here considered to show that the proposed method can handle both the instability and the asymmetry issues. The asymmetric stiffness matrix K is the sum of two parts: the symmetric part s K and the asymmetric part a K . The latter is caused by the pre-compression normal force and the friction force that acts at the slider-belt interface. The resulting system matrices are reported in Equation (29): The model parameters have been assumed as in [19]: Indeed, a pair of complex conjugate poles with real part appears: 0.00 ± 8.73i [19]. The control is designed to exactly place two pairs of dominant poles into prescribed locations of the complex plane, to ensure stability and the desired dynamic responses, while the remaining two pairs of dominated poles are constrained into an LMI.
The results are compared with the standard PPP methods used in vibration control, such as for example those proposed in [6] and [7]. Traditional PPP techniques differ from the RPP technique proposed in this paper, since they assign p poles while the remaining 2N p − poles are kept The control is designed to exactly place two pairs of dominant poles into prescribed locations of the complex plane, to ensure stability and the desired dynamic responses, while the remaining two pairs of dominated poles are constrained into an LMI.
The results are compared with the standard PPP methods used in vibration control, such as for example those proposed in [6,7]. Traditional PPP techniques differ from the RPP technique proposed in this paper, since they assign p poles while the remaining 2N − p poles are kept unchanged, i.e., 2N − p open-loop poles will coincide with the closed-loop ones. To this purpose, the two desired complex conjugate pole pairs in Table 4 are assigned at −1 ± 9i and −1 ± 13.5i, while the remaining two complex conjugate pole pairs (−0.51 ± 16.75i and −0.19 ± 19.86i) are kept unchanged. The control gains computed by this benchmark, denoted as k B , are reported in Table 5.  To ensure a fair comparison, the proposed method performs RPP within the LMI Re(µ i ) ≤ −0.19 that includes the two open-loop pairs of poles at −0.51 ± 16.75i and −0.19 ± 19.86i. Clearly, k B is a possible solution of the assignment problem, although it is expected that the proposed method is able to improve it. The gains after the first and second stage of the proposed method are summarized in Table 5 and compared with k B . The desired and the obtained poles are shown in Table 4 and in Figure 4, where they are compared with those provided by the traditional partial pole placement as well.  possible solution of the assignment problem, although it is expected that the proposed method is able to improve it. The gains after the first and second stage of the proposed method are summarized in Table 5 and compared with B k . The desired and the obtained poles are shown in Table 4 and in Figure 4, where they are compared with those provided by the traditional partial pole placement as well.
The comparison of the control gains in Table 5 reveals that RPP leads to a smaller norm compared to PPP:  To better compare the control effort, the actuator energy in a free response of the controlled system tracking the zero-vibration condition can be evaluated ( [19,27]), by assuming a meaningful set of initial conditions of the state and by simulating a time interval T covering several oscillation periods: The comparison of the control gains in Table 5 reveals that RPP leads to a smaller norm compared to PPP: k B 2 = 177.9, while k 0 2 = 26.7 and k 2 = 67.6.
To better compare the control effort, the actuator energy in a free response of the controlled system tracking the zero-vibration condition can be evaluated ( [19,27]), by assuming a meaningful set of initial conditions of the state and by simulating a time interval T covering several oscillation periods: The following velocity and displacement initial conditions have been set: . q(0) = 1 · 10 −3 1 1 1 1 T and q(0) = 1 · 10 −3 1 1 1 1 T . The system response has been computed by numerically solving the system ODEs with the ODE45 solver, as shown in Figures 5 and 6. The resulting control effort is drastically lower, since it decreases from 0.0123 J to 0.0053 J, as shown in Figure 7. . The system response has been computed by numerically solving the system ODEs with the ODE45 solver, as shown in Figures 5 and 6. The resulting control effort is drastically lower, since it decreases from 0.0123 J to 0.0053 J, as shown in Figure 7.    The benefits of gain reduction can be also appreciated through the theory of robust control. Indeed, although the proposed method is not explicitly addressed to robustness, the gain decreases due to the RPP in lieu of the exact pole placement can be very effective to robustness. The robustness of the two controllers is compared by means of the H-infinity norm of the transfer function H ZW from the model uncertainty to the nominal state .

q(t) q(t)
, i.e., through H ZW ∞ . By assuming, without loss of generality, the structure of a state-multiplicative perturbation, H ZW is defined through the following state space realization [4]: In accordance with the Small Gain Theorem, smaller values of H ZW ∞ denote higher robustness with respect to the uncertainty. The application of such a theory to the example under investigation shows that the proposed method leads to higher robustness compared to the benchmark since H ZW (k) ∞ = 20.5 while H ZW (k B ) ∞ = 190.2. Figure 8 corroborates the results, by showing the maximum singular value for both the controllers: σ(H ZW (k)) is smaller than σ(H ZW (k B )) in the whole frequency range, leading to a system with increased robustness.
showing the maximum singular value for both the controllers: σ ZW H k in the whole frequency range, leading to a system with increased robustness. All these results are remarkable and proves the advantages and the novelty of the method proposed in this paper over traditional PPP. Indeed, relaxing the constraint on the less relevant poles allows optimizing a secondary task without worsening the system response, or even by improving it.

Test-Case 3: Aircraft Wing System
The third test-case is another common benchmark in the literature on pole assignment. It consists of an aircraft wing under an air stream, as proposed in [15,28,29]. The system matrices are: The system is unstable, and the stiffness and damping matrices are both asymmetric, and are hence decomposed into a symmetric and an asymmetric part: K=K s +K a and C=C s +C a .
The control design task prescribed to assign the dominant complex conjugate pole pair µ 1,2 = −1.500 ± 3.000i, while the remaining closed-loop poles µ i must belong to the LMI region proposed in [15], which is a disk of radius equal to 1 and center at coordinates (−3,0).
The gain k 0 , computed by means of the receptance method for asymmetric systems (see Section 3.2), is unable to stabilize the unstable open-loop system, leading to an unstable right half complex conjugate pole pair. Finally, the gain k computed through the second stage can stabilize the system and cluster the unassigned poles into the prescribed LMI region. The gains are reported in Table 6 while the open-loop and closed-loop poles are reported in Table 7 and shown in Figure 9. It should be pointed out that, unlike [15], the method proposed here is able to modify the open-loop pole in −0.918 ± 1.7606i, by moving it, for example, to −1.500 ± 3.000i. This modification is effective in speeding up the closed-loop system response due to the increase of the pole natural frequency compared to the open-loop poles. This is a relevant feature and novelty of the method proposed in this paper, compared to the methods already available in the literature.

Test-Case 4: Four-Bar Flexible-Link Multibody System
The method has been applied to the four-bar flexible-link multibody system proposed in [4,17]. The dynamic of these systems is nonlinear and disturbed by gravity forces. Even if the dynamic of such a system is non-linear, linearized or piecewise-linear models have been proved to be effective in the synthesis of control schemes [4,17,[30][31][32].
To apply the proposed method, the system has been modeled through ODE by exploiting an Equivalent Rigid Link System (ERLS) from which elastic displacements are defined. Euler-Bernoulli beam finite elements are adopted to model the link flexibility. The resulting non-linear model, formulated through the Principle of Virtual Work, has then been linearized about an equilibrium configuration of interest. The model accounts for the mutual coupling between the "rigid-body motion" of the ERLS and the small elastic displacements with respect of the ERLS itself. Additionally, it accounts for the effect of gravity on both the "rigid-body" and the elastic dynamics. More details are provided in [4,17,28]. The finite element model of the mechanism is shown in Figure 10 while the system parameters are reported in Table 8. All the links, except the frame, are flexible.

Test-Case 4: Four-Bar Flexible-Link Multibody System
The method has been applied to the four-bar flexible-link multibody system proposed in [4,17]. The dynamic of these systems is nonlinear and disturbed by gravity forces. Even if the dynamic of such a system is non-linear, linearized or piecewise-linear models have been proved to be effective in the synthesis of control schemes [4,17,[30][31][32].
To apply the proposed method, the system has been modeled through ODE by exploiting an Equivalent Rigid Link System (ERLS) from which elastic displacements are defined. Euler-Bernoulli beam finite elements are adopted to model the link flexibility. The resulting non-linear model, formulated through the Principle of Virtual Work, has then been linearized about an equilibrium configuration of interest. The model accounts for the mutual coupling between the "rigid-body motion" of the ERLS and the small elastic displacements with respect of the ERLS itself. Additionally, it accounts for the effect of gravity on both the "rigid-body" and the elastic dynamics. More details are provided in [4,17,28]. The finite element model of the mechanism is shown in Figure 10 while the system parameters are reported in Table 8. All the links, except the frame, are flexible.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 22 Figure 10. Four-bar linkage: system model (test-case 4) (reproduced from [16]). Sixteen DOFs are adopted to model the system: one coordinate is the rotation of the ERLS, to represent the "rigid-body" motion. The remaining DOFs are the displacement of the 15 elastic coordinates. The mechanism is moved by a motor driving the crank rotation. Such an actuator ensures controllability of all the modes, although higher frequency modes are less controllable [17]. The open-loop system is unstable due to the presence of gravity, as it can be inferred from Table 9. The system stiffness and damping matrices are both asymmetric and can be decomposed into a symmetric and an asymmetric part, as done in the second test-case: s a = + K K K and s a = + C C C .   Sixteen DOFs are adopted to model the system: one coordinate is the rotation of the ERLS, to represent the "rigid-body" motion. The remaining DOFs are the displacement of the 15 elastic coordinates. The mechanism is moved by a motor driving the crank rotation. Such an actuator ensures controllability of all the modes, although higher frequency modes are less controllable [17]. The open-loop system is unstable due to the presence of gravity, as it can be inferred from Table 9. The system stiffness and damping matrices are both asymmetric and can be decomposed into a symmetric and an asymmetric part, as done in the second test-case: K = K s + K a and C = C s + C a . To tackle the risk of ill numerical conditioning of the SDP programming, due to the different magnitudes of the entries of the mass, damping and stiffness matrices, the system matrices are scaled by defining the scaling parameters γ and β: which leads to the scaled receptance matrix (adopted in the first stage): The issue of numerical conditioning is exacerbated in this kind of problem due to the medium-large number of DOFs and to the simultaneous presence of low and high frequency vibrational modes. Therefore, as shown in [4], some numerical methods for pole placement often might fail.
Accordingly to the scaled receptance matrix, scaled gain vectors f 0 and g 0 are computed in the first stage and then in the unscaled model: In the second stage, the state space system matrices in Equation (12) are scaled accordingly to the definitions provided in Equation (33) The LMI conditions have been scaled too. The following scaling parameters have been here adopted: γ = 1 · 10 5 and β = 1 · 10 −4 .
The control specifications prescribe to assign the two real poles, that in practice represent the "rigid-body" motion, and the three lowest-frequency complex conjugate pole pairs, as reported in Table 9. The remaining closed-loop poles must belong to the LMI ensuring that ξ(µ i ) ≥ 0.05.
The gain k 0 and k are reported in Table 10, while the open-loop and the closed-loop poles are summarized in Table 9 and shown in the complex plane in Figure 11. The application of k 0 alone leads to an unstable pole pair at 28.6 ± 244.7i. The second stage of the method ensures the stabilization of the unstable pole pair, which is placed at −74.0 ± 215.7i. At the same time, it damps all the remaining unassigned poles accordingly with the requirements, as shown in the low frequency detail of the pole map in Figure 12

Conclusions
In this work a novel, two-stage, pole placement method is proposed and an embedded a-priori stabilizing condition is introduced. In the first stage the receptance method for symmetric or asymmetric systems is applied to compute the gains that exactly assign a subset of the system poles. These poles are usually the dominant ones, such as those with the lowest frequency or with higher controllability and observability. In the second stage, by leverage of a first-order model formulated with a reduced state and by exploiting the Inclusion Principle, the unassigned poles are clustered into some region of the complex plane exploiting the LMI theory. An SDP programming is solved to reduce the control effort, by exploiting the additional degrees of freedom introduced by the RPP, in lieu of the "exact" pole placement. A novel formulation has been proposed by exploiting the aforementioned mathematical tools and formulations and some methods to improve the numerical reliability and solvability of the problem.
Four meaningful and challenging test-cases, also taken from common benchmark in the literature on pole placement in vibrating systems, are proposed. In the first test-case, a simple three-mass system and some features of the method are introduced. In the second test-case, it is provided that the method handles unstable (or marginally stable) and asymmetric systems as well. Moreover, the comparison with traditional partial pole placement techniques, under comparable requirements, shows the benefits of relaxing the requirements on the less relevant poles. Not only can the control effort be reduced and the robustness increased, but the dynamic response of the system can also be improved in term of reducing the settling time. The third test case tackles the control of an aircraft wing under an air stream by stabilizing its unstable open-loop behavior and by improving its closed-loop dynamic response. The fourth test case is the challenging linearized model of a nonlinear flexible multibody system perturbed by gravity, which makes it unstable. Due to the medium-large number of degrees of freedom, to the presence of an unstable pole, to the presence of real poles and to a wide spectrum that includes high-frequency poles, this test case is numerically ill-conditioned and hence very challenging.
Additionally, if PPP is not performed correctly, spillover on the unassigned poled might lead to unstable poles. Again, the proposed method has been shown to be effective.
The results in the four test-cases, and the comparison with the existing literature too, confirm the benefits of the proposed method.