An existence condition for the inverse dynamics solution of a slewing Euler–Bernoulli beam

Mechanism and Machine Theory 46 (2011) 845–860

Contents lists available at ScienceDirect

Mechanism and Machine Theory j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / m e c h m t

An existence condition for the inverse dynamics solution of a slewing Euler–Bernoulli beam Qiao Sun ⁎ Department of Mechanical and Manufacturing Engineering, University of Calgary, Calgary, Alberta, Canada

a r t i c l e

i n f o

Article history: Received 28 July 2010 Received in revised form 1 March 2011 Accepted 7 March 2011 Available online 2 April 2011 Keywords: Euler–Bernoulli beam Causal inverse dynamics solution Stability Existence condition

a b s t r a c t Causal inverse dynamics problem of a slewing flexible beam has been identified as being illposed because it violates the stability of solution. Such problem has been studied intensively for the past three decades due to its application to the tip tracking control of flexible-link manipulators. A well known remedy has been to modify the location of the boundary point, or the output function in tracking control. This paper re-examines the problem by analyzing a closed form solution. Such a solution is made possible by truncating the assumed modes beam deflection model after the first mode. An existence condition for the solution is established, which links to the deflection mode shape, payload inertia, and the location of boundary point. This condition indicates that for certain system parameters, causal solutions may exist. Existing remedies as suggested by many published results can be explained by the condition. Extension to a multi-mode model is discussed and examined through numerical simulations. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction We consider, in the horizontal plane, a slewing flexible beam of Euler–Bernoulli type that is fixed to a rigid hub at one end and carries a payload at the other end, as shown in Fig. 1. An external moment τ is applied to the hub to directly control the angular displacement θ(t) at time t (t ≥ 0). Due to beam elasticity, the displacement of beam's free end (tip), measured by β(t) is different from θ(t). Finding angular displacement β in response to an applied moment function τ(t) is known as the forward dynamics problem. A typical solution procedure involves approximating the continuous beam deflection by a finite number of spatial degrees of freedom to obtain a set of algebraic and ordinary differential equations. Boundary conditions, both geometric and natural, are satisfied if the deflection shape functions are those of normal modes. Increasing the number of spatial degrees of freedom increases solution accuracy with convergence to the exact solution ensured. Conversely, for a given C2 function βd(t) defined over t ∈ [0, tf] to be imposed on a boundary point (tip here), finding the control moment τ is known as the inverse dynamics problem. However, it turns out to be significantly more difficult. Following the same procedure as that used in the forward dynamics but switching input with output, solving the initial value problem of the same ODE's results in unbounded control torque τ. In the past, a large amount of effort has been made to remedy the problem although little has been understood regarding the nature and cause of the problem. One family of remedies is to seek non-causal solutions first proposed by Bayo et. al. [1–3]. The proponents of this approach used time delay between the input and output to justify non-causality. However, the beam equation of an Euler–Bernoulli type, being non-hyperbolic, permits infinite speed of disturbance propagation [4]. Time delay is not captured by the governing equations to justify the use of time delay as the cause of solution behavior. Another family of remedies is to

⁎ Tel.: + 1 403 220 4141; fax: + 1 403 282 8406. E-mail address: [email protected]. 0094-114X/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmachtheory.2011.03.001

846

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

Fig. 1. A slewing flexible beam.

redefine the output function [5–7]. In particular, boundary point is either chosen at an interior point rather than the tip of the beam or as a linear combination of joint angles and elastic coordinates. Proponents of this approach identify the reason of unbounded boundary control solution as due to the non-collocation between input and output [8]. Again, spatial separation of the boundary points is considered as the root cause of the solution instability. In special applications where actuation redundancy is available, an exact solution can be obtained with control actions applied at both boundary points [9]. In tip tracking control applications, solutions are further complicated by the instability of the closed-loop system when tip tracking error is to be zeroed. Most of the published results of remedies are guided by numerical modeling which makes the identification of an existence condition and the cause of solution instability difficult to understand. In this study, we seek to understand the nature of the problem by examining the governing equations explicitly expressed in system parameters. This enables the determination of an existence condition in physical parameters. From the existence condition, we can identify possible remedies to ensure a bounded inverse problem solution. One of them points to the existing method of boundary point re-definition. Strictly speaking, this is not a solution as it modifies the definition of the “tip”. Another remedy is by having certain amount of payload rotary inertia [5]. Without a unified existence condition, it would be hard to see a logical connection between the two proposed remedies. In fact, both remedies are to compensate for the effect of deflection shape properties that causes the dissatisfaction of the existence condition. The parametric analysis is conducted based on a single-mode model. Numerical simulations are carried out to verify conclusions led by the theoretical analysis and to illustrate the utility of the single-mode model in generating solutions for a more general model. The rest of the paper is organized as follows. In Section 2, model equations are introduced. We first extend the clamped-free beam deflection model [10] by including payload mass and rotary inertia. Mode shape and frequency equations are parameterized in beam and payload inertia properties. Dynamics and kinematics equations are then formulated using the assumed-modes deflection model. Such a model, while truncated after the first mode, makes it possible for the inverse dynamics problem to be dealt with parametrically. Subsequently in Section 3, solution existence condition is established and closed form solutions are formulated. In addition, solutions are examined for their tracking performance through forward dynamics. We show the exact regeneration of the tip trajectory with a sinusoidal acceleration profile with a particular choice of initial conditions in the inverse dynamics solution. In Section 4, simulations results are given to verify conclusions derived from the analytical solutions. Through numerical simulations, we also show the effectiveness of the single-mode model in providing inverse dynamics solutions to a multi-mode system. Concluding remarks are given in Section 5. 2. Model equations Refer to Fig. 1, the slewing elastic beam has a length L, constant mass density ρ, uniform cross sectional area A and bending stiffness EI. Let M be the total mass of the beam: M = ρAL. The beam is clamped to a rigid hub with a moment of inertia Jh much higher than that of the beam Jb = 13 ML2 . A payload with a total mass Me and rotary inertia Je is located at the free end of the beam. Frame Oxy is a moving frame fixed to the hub and therefore θ(t) measures the rigid-body motion with respect to the inertial frame OXY. The system lies in a horizontal plane, permitting the convenient omission of gravity. 2.1. Fundamental frequency and mode shapes Let v(x, t) be the relative displacement of a cross section located at a distance of x from the origin of the moving frame Oxy. Assuming small beam deflection, v(x, t) represents the transverse beam deflection. Neglecting rotary inertia and shear

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

847

deformation as hinted by the choice of Euler–Bernoulli beam, a clamped-free beam deflection is governed by the following equation [10] 4

EI

2

∂ v ∂ v ðx; t Þ + ρA 2 ðx; t Þ = 0: ∂x 4 ∂t

ð1Þ

Geometric boundary conditions specified at the clamped end and natural boundary conditions specified at the payload end are: vð0; t Þ = 0;

ð2aÞ

∂v ð0; t Þ = 0; ∂x

ð2bÞ

2

3

EI

∂ v ∂ v ðL; t Þ = −Je ; and ∂x2 ∂x∂t 2

ð2cÞ

EI

∂3 v ∂2 v ðL; t Þ = Me 2 ðL; t Þ: 3 ∂x ∂t

ð2dÞ

Here, we have adopted the constrained modal analysis where the influence of rigid body motion in beam deflection is ignored owing to the assumption of large inertia ratio between hub and beam [11]. Using the standard method of separation of variables, v(x, t) = ϕ(x)q(t), Eq. (1) can be split into two ordinary differential equations, one in space and the other in time: ð4Þ

EI ϕ q̈ 2 =− =ω : ρA ϕ q

ð3Þ

In the above, superposed dots are used to denote total derivatives with respective to time variable t and primes are used to denote total derivatives with respect to space variable x. The above equation suggests that q ̈ = −ω2 qðt Þ and therefore, ∂2 v = −ω2 vðx; t Þ. Thus, boundary conditions (2) can be simplified and imposed on the spatial function ϕ(x) only. These are ∂t 2 ϕð0Þ = 0;

ð4aÞ

ϕ′ ð0Þ = 0;

ð4bÞ

2 EIϕ″ ðLÞ = Je ω ϕ′ ðLÞ;

ð4cÞ

2

EIϕ‴ ðLÞ = −Me ω ϕðLÞ:

ð4dÞ

The general expression for the spatial function that satisfies Eq. (3) is ϕðxÞ = C1 ðcos γx + cosh γxÞ + C2 ðcos γx − cosh γxÞ + C3 ðsin γx + sinh γxÞ + C4 ðsin γx− sinh γxÞ

ð5Þ

ρA 2 ω . with γ defined according to γ = EI Applying boundary conditions (4a) to (4d) in an attempt to solve for mode shape coefficients in Eq. (5) results in the following fundamental or modal frequency equation through parameter γ: 4


 4 2 4 2 1−ðγLÞ ðR eÞ + 1 + ðγLÞ ðR eÞ




  1 4 2 4 2 − γLR 1 + ðγLÞ e tan γL− 1−ðγLÞ e tanh γL = 0 : ð6Þ cos γL cosh γL rffiffiffiffiffiffiffi Je = eL to denote the In the above, we use R = Me/M to denote the mass ratio of payload with respect to the beam and re = Me radius of gyration of the payload in ratio of the beam length. Both payload mass and rotary inertia modify the beam's modal frequency. Table 1 lists the first two roots of Eq. (6), γ1L and γ2L, for different values of R and e. They determine the first two EI ðγLÞ4 . The general effect is that as R and e increase, fundamental frequencies fundamental frequencies of the beam by ω2 = ρAL4 decrease. This can be observed in the table by moving from left to right, or top to bottom. The first column represents a beam with no payload and therefore γ1L and γ2L stay constant moving down the column. The boundary condition (4) also allows us to define the deflection shape function (5) as:   sin γL− sinh γL + γLRðcos γL− cosh γLÞ ðsin γx− sinh γxÞ : ϕðxÞ = C2 cos γx− cosh γx + cos γL + cosh γL− γLRðsin γL− sinh γLÞ

ð7Þ

848

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

Table 1 First two fundamental frequencies influenced by R and e. γ1L

R

γ2L 0 0.001 e

0.01 0.1 1

0

0.01

0.1

1

10

1.8751 4.6941 1.8751 4.6941 1.8751 4.6941 1.8751 4.6941 1.8751 4.6941

1.8568 4.6497 1.8568 4.6497 1.8568 4.6496 1.8564 4.6394 1.8227 3.7745

1.7227 4.3995 1.7227 4.3995 1.7227 4.3987 1.7203 4.3215 1.5078 2.6064

1.2479 4.0311 1.2479 4.0311 1.2479 4.0271 1.2424 3.6386 0.9316 1.8414

0.73578 3.9385 0.73578 3.9381 0.73574 3.9053 0.73178 2.4612 0.5294 1.1015

Here C2 is arbitrary and we choose C2 = 1/ϕ(L) so that all shape functions are scaled to result in ϕ(L) = 1. Eq. (7) defines the deflection shape function corresponding to a particular modal frequency. Referring to Appendix A, the above mode shapes are orthogonal with one another, as one would expect so. Since Eq. (6) permits infinite number of solutions for modal frequency, a general solution for the deflection function is the sum of all fundamental solutions,
 ∞ vðx; t Þ = ∑ϕj ðxÞ Aj cos ωj t + Bj sin ωj t

ð8Þ

j=1

where constants Aj and Bj can be determined from the initial conditions of the beam. Practically, only a finite number of modes are used to approximate beam deflection. Eq. (8) can be interpreted as discretizing a continuous beam deflection function with a finite number of fundamental modes. The contribution of each mode to the entire beam deflection at any time t is indicated by a time dependent magnitude qj, that is, n

vðx; t Þ = ∑ϕj ðxÞqj ðt Þ:

ð9Þ

j=1

Here n is the number of modes. Comparing Eq. (9) with Eq. (8), it should be understood that qj oscillates in its corresponding modal frequency ωj. 2.2. System equations To derive the dynamics equations for the system in Fig. 1, we start with the total kinetic energy: !2    2 1 1 L ∂v 2 1 ∂v 1 ∂2 v 2 ˙ ˙ ˙ ˙ dx + Me Lθ + T = JH θ + ∫ ρA xθ + ðL; t Þ + Je θ + ðL; t Þ : 2 2 0 2 2 ∂t ∂t ∂x∂t

ð10Þ

Notice that the material's internal damping has been omitted for simplicity. This is justified in this study since the existence of damping is not essential to the existence of a solution but rather to ensure a decayed oscillation should the solution exist. Applying Eq. (9) to the beam deflection v in Eq. (10) results in n L 1 2 2 T = ρA∫ x θ˙ + 2xθ˙ ∑ ϕj ðxÞ q˙ j + 0 2 j=1 n 1 2 + Je θ˙ + 2θ˙ ∑ ϕj′ðLÞ q˙ j + 2 j=1

!2 !

n

∑ ϕj ðxÞ q˙ j

j=1

n

! !

∑ ϕj′ðLÞq˙ j

2

j=1

n 1 1 2 2 2 dx + JH θ˙ + Me L θ˙ + 2Lθ˙ ∑ q˙ j + 2 2 j=1

n

!2 !

∑ q˙ j

j=1

ð11Þ

:

Note that ϕj(L) = 1 is implied in the above. The potential energy is due to the bending deformation in the beam. Assuming small deflections, it is

V=

1 L ∫ EI 2 0

2

∂ v ∂x2

!2 dx =

n L 1 EI∫ ∑ ϕj″ðxÞqj 0 2 j=1

!2 dx:

ð12Þ

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

849

The system has n + 1 generalized coordinates, θ and qj, j = 1, … n. Using the Euler–Lagrange formulation and ignoring higher order deflection terms, the equations of motion are obtained as:

Mθθ θ̈ + Mθq q̈ = τ

ð13aÞ

Mqθ θ̈ + Mqq q̈ + Kq q = 0;

ð13bÞ

where 2

Mθθ = Jh + Jb + Me L + Je n L o Mθq = row ∫ ρAxϕj ðxÞdx + Me L + Je ϕj′ðLÞ 0

Mqθ =

MTθq ;

n L o n o 2 Mqq = diag ∫ ρAϕj ðxÞdx + Me + Je ϕi′ðLÞϕj′ðLÞ; ∀i≠j ði; jÞ  0
2  L Kq = diag ∫ EI ϕj″ðxÞ dx : 0

In the above, rowf⋅g∈R1×n is a row vector; diagf⋅g∈Rn×n is a diagonal matrix; and f⋅gði;jÞ ∈Rn×n gives the element (i, j) in an n × n matrix. Given an external driving moment τ and system's initial conditions θ(0) and q(0), one can determine the generalized coordinates θ and q at time t (t ≥ 0) using Eq. (13) by integrating forward in time. Subsequently, one can find the motion of the beam's free end using the following kinematics relation: β=θ+

vðL; t Þ = θ + Ψq q L

ð14Þ

1 ½1 1 ⋯ 1 since ϕj(L) = 1, j = 1, …, n. Eq. (13) governs the motion of the system. Depending on the input and L output, it can be used to solve different problems. One particularly of interest is the inverse dynamics problem which is the focus of the following sections.

where Ψq =

3. Inverse dynamics solutions Given a C2 function βd(t), 0 ≤ t ≤ tf, and initial values of the state variables θ(0), q(0), θ˙ ð0Þ, and q˙ ð0Þ, find a moment function τ(t), 0 ≤ t ≤ tf, so that the tip of the elastic beam displaces according to the desirable trajectory: β = βd(t), 0 ≤ t ≤ tf. This is defined as d the inverse dynamics problem. Since βd(t) is a C2 function, β̈ ðt Þ is continuous. An equivalent problem is to find τ(t), 0 ≤ t ≤ tf so d that β̈ = β̈ ðt Þ for 0 ≤ t ≤ tf with initial values of β(0) and β˙ ð0Þ consistent with the initial states θ(0), q(0),θ˙ ð0Þ, and q˙ ð0Þ of the beam. 3.1. Solution formulation through a single-mode model In order to gain insight into the inverse dynamics problem, it is convenient if a parametric closed-form solution can be formulated. To this end, we truncate the assumed-modes model after the first mode. Qualitatively, such a model preserves the fundamental characteristics of the problem, i.e., solution existence factors, although coupling effects among multiple modes cannot be explored through this model. For notation simplicity, we also drop subscripts that are used to index the elastic coordinates. Thus, beam deflection is expressed in the moving frame Oxy as: vðx; t Þ = ϕðxÞqðt Þ;

ð15Þ

where ϕ(x) and q(t) are the shape function and elastic coordinate of the first normal mode, that is, n = 1. For convenience, we x normalize the space variable by beam length, η = . Thus 0 ≤ η ≤ 1 represents the space domain and η = 1 indicates the boundary L point at the tip. In addition, we introduce the following notations:

2 1 1 2 1 I1 = ∫ ηϕðηÞdη; I2 = ∫ ϕ ðηÞdη; I3 = ∫ ϕ″ ðηÞ dη: 0

0

0

ð16Þ

850

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

These integrals have no units. They can be pre-calculated if the shape function is known. Derivatives of the shape function with respect to x and η can be converted according to: ϕ′(x) = ϕ′(η)/L and ϕ″(x) = ϕ″(η)/L2. Owing to the simple deflection model (15), coefficients in the equations of motion (13) are scalars. They can be expressed using notions from Eq. (16) as: Mθq = MLI1 + Me L + Je ϕ′ ð1Þ = L; Mqq = MI2 + Me + Je ðϕ′ ð1ÞÞ2 = L2 ;

ð17Þ

EI Kq = 3 I3 : L Remark 3.1. The natural clamped-free beam vibration can be recovered from Eq. (13b) by setting θ̈ = 0. This corresponds to a sffiffiffiffiffiffiffiffiffi Kq beam with no rigid body rotation. The natural frequency is indicated by the dynamics Eq. (13b) and is thus ωq = . From Mqq Eq. (A5) in Appendix A with i = j = 1 and using Eqs. (16) and (17), one confirms that ωq = ω1. To solve the inverse dynamics problem, we start with the kinematics Eq. (14). After twice differentiating with respect to time, it d

is used to eliminate θ̈ from Eq. (13b). Letting g ðt Þ = β̈ ðt Þ, we obtain the following ODE of a single variable that governs the elastic degree of freedom if the tip is to follow a prescribed displacement function βd(t), Mi q ̈ + Kq q = −Mqθ g ðt Þ

ð18aÞ

where Mi = Mqq −Mqθ Ψq = M ðI2 −I1 Þ + MRe

2






2 ϕ′ ð1Þ −ϕ′ ð1Þ :

ð18bÞ

Notice that Ψq = 1/L. Crucial to the property of Eq. (18a) is the sign of Mi. The existence of a bounded solution for a bounded input function g(t) requires the existence of the homogeneous solution. Since Kq N 0, a necessary and sufficient condition is that Mi N 0. In detail, the condition translates to

2 I2 −I1 + R e ϕ′ð1Þ ϕ′ð1Þ−1 N 0:

ð19Þ

This condition puts a necessary requirement on beam deflection mode shape. The effect of payload mass and rotary inertia is not only through parameters R and e, they also modify the shape function and therefore values of I1, I2, and ϕ ' (1). Fig. 2 shows several mode shapes with different values of R and e. Recall that all shape functions are scaled to yield ϕ(1) = 1. Except for large payload inertia, e.g. R = 5, and e = 1 in Fig. 2, the shape function curve situates below the straight line that connects the two end points of the beam, ϕðηÞ ≤ η;∀0 ≤ η ≤ 1:

ð20Þ

1

R=0, e=0 R=5, e=0.1 R=5, e=1

φ(η)

0.8

φ(η)=η

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5

η

0.6

Fig. 2. Mode shape function.

0.7

0.8

0.9

1

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860 1

1

1

0

0

0

In this case, ∫ ϕ2 dη ≤ ∫ ηϕdη ≤ ∫ η2 dη = I2 ≤ I1 ≤

851

1 , which means the following is true: 3

1 : 3

ð21Þ

We can now explain why in general, the inverse dynamics problem has no solution. If payload rotary inertia is zero, e = 0, condition (19) simply becomes I2 −I1 N 0:

ð22Þ

In the absence of a payload, Eq. (20) is true and Eq. (21) applies. Condition (22) is immediately violated. This conclusion can also be verified by taking γ1L = 1.8751 from Table 1 and substituting into Eq. (7) for the mode shape function and carrying out integrations according to Eq. (16). One obtains I1 = 0.2844 and I2 = 0.2500, confirming the violation of the existence condition (19). It is worth noting that although beam physical parameters affect beam natural frequency and mode shape function, they do not change the fact that the first mode shape curve stays under the connecting line of the end points. Even with the inclusion of a large payload mass, as long as the payload rotary inertia is zero, the shape curve remains under the line. Large values of payload rotary inertia can change the mode shape curve to cross the connecting line which makes I2 − I1 term less negative. In addition, only when the payload rotary inertia is not zero are extra terms added to I2 − I1 in Eq. (19) which gives more possibilities for condition (19) to be satisfied. Using Eq. (18b), it is straightforward to explain the choosing of an imaginary boundary point in Reference [7]. There, a “reflective” tip point is used meaning Ψq = − 1/L. The effect is to turn the negative term in Mi of Eq. (18b) into a positive term. The existence condition becomes I2 + I1 N 0 for R = 0, which is always guaranteed. We could also explain the effect of moving the boundary to an interior point of distance d from the hub. In such case, Ψq = ϕ(d)/d. The existence condition (19) with R = 0 becomes: I2 −I1 ϕðdÞðL = dÞN 0:

ð23Þ

With the same reason that the shape curve is below the connecting line of end points, a simple geometric relation shows that ϕ (d) b (ϕ(L)/L)d. Therefore, ϕ(d)(L/d) b 1 which reduces the magnitude of the negative term when comparing Eq. (23) with Eq. (22). In fact, one can easily solve for a critical value of d by finding the roots of I2 − I1ϕ(d)(L/d) = 0. For any desired trajectory specified for a point left of d, stability condition (23) can be satisfied. If d = 0, Ψq = 0, Mi = I2 N 0 explains the choice of a joint-based output function [6] that guarantees a stable solution. However, modifying the boundary point modifies the problem. Exact regeneration of the prescribed tip trajectory is not expected. Numerical solutions can be sought to determine the values of R and e that uphold condition (19). In Fig. 3, the solid curve separates the R-e space into acceptable and unacceptable regions. Any point of (R, e) above the solid curve results in the satisfaction of condition (19) and therefore a bounded solution of Eq. (18a). It is also important to realize, from Fig. 3, that a minimum radius of gyration of nearly 10% of the beam length is required in any case. In practical applications such as robotic manipulators, it is likely that a gripper is located at the free end of the beam with non-zero rotary inertia. Eq. (19) may become a design consideration for a flexible manipulator.

1 0.9

1 mode 2 modes

0.8 0.7

re/L

0.6 0.5 0.4 0.3 0.2 0.1 0

0

1

2

3

4

5

6

7

Me/M Fig. 3. Boundaries of bounded elastic motion.

8

9

10

852

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

Due to condition (19), Eq. (18a) can be re-written as: 2

q ̈ + ωi q = Ri g ðt Þ where ω2i =

Kq Mi ,

Ri = −

ð24Þ Mqθ Mi

:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Remark 3.2. The natural frequency of Eq. (24) as given by ωi is higher than the beam modal frequency ω1 = Kq = Mqq since Mi b Mqq. As it will become clear later, ωi is also different from the slewing beam natural frequency. Therefore, solutions to the inverse dynamics problem entail a forced oscillation at a frequency different from beam's natural frequency. This is re-assuring since any external excitation in beam's natural frequency could potentially cause unbounded response. If the system permits a bounded solution to the inverse dynamics problem, that is, with the existence condition (19) satisfied, a closed-form solution of Eq. (18a) is

qðt Þ =

    q˙ 0 R t R t + i ∫ g ðιÞ cos ωi ι dι sin ωi t + q0 − i ∫ gðιÞ sin ωi ι dι cos ωi t; ωi ωi 0 ωi 0

ð25Þ

where q0 = qð0Þand q˙ 0 = q˙ ð0Þ are initial conditions of the elastic degree of freedom. Subsequently, one can proceed with finding the corresponding angular displacement of the hub, θ(t), from the kinematics Eq. (14) with Ψq = 1/L, t

ξ

θðt Þ = ∫ ∫ g ðιÞdιdξ−Ψq 0 0

    q˙ 0 R t R t + i ∫ gðιÞ cos ωi ι dι sin ωi t−Ψq q0 − i ∫ gðιÞ sin ωi ι dι cos ωi t : ωi ωi 0 ωi 0

ð26Þ

Finally, once both degrees of freedom are determined, Eq. (13a) can be used to find the solution for a control torque. Simple algebraic manipulations conclude that the torque is a linear combination of the given tip acceleration function and the resulting elastic motion from Eq. (25): τðt Þ = Jq g ðt Þ + Cq qðt Þ

ð27Þ

where
 2 Mi−1 Jq = Mθθ Mqq −Mθq
 Cq = − Mθq −Mθθ Ψq Kq Mi−1 : Due to Eq. (25), it can be confirmed from the above equation that the control torque applied to the hub does not contain oscillatory elements that would excite the beam natural mode. 3.2. Boundary tracking performance From the inverse dynamics solution given by Eqs. (25) and (27), one realizes that in order to achieve a prescribed motion at the beam's boundary point (tip), the control torque needs to force a beam oscillation at a frequency different from its natural frequency so that any unwanted oscillation at the tip can be zeroed. We now examine the tracking performance of the inverse dynamics solution by substituting Eq. (27) back into the beam's governing Eq. (13a) and finding the response motion of the boundary point. To do this, we again first eliminate the rigid degree of freedom from the governing equations, this time by using Eq. (13a) and bringing it to Eq. (13b). To avoid confusion of variables with those used in the inverse problem solutions, we use ξ in places of q for the elastic coordinate. The following equation describes the motion of the elastic degree of freedom in response to the applied torque: 2 ξ̈ + ωf ξ = Rf τðt Þ

where 2

−1

ωf = Mf Kq −1 Rf = −Mf−1 Mqθ Mθθ ; and −1 Mf = Mqq −Mqθ Mθθ Mθq :

ð28Þ

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

853

Similar to the inverse dynamics problem, the existence of a solution to the above equation requires that Mf N 0. Using Eqs. (16) and (17), this condition can be detailed as
2 I1 + R + R e2 ϕ′ ð1Þ

2 N 0: I2 + R + R e ϕ′ ð1Þ − b 2 1 3 + 3 + R + Re 2

ð29Þ

Jh is the moment of inertia ratio between the hub and the beam. In a simple setting when there is no payload Jb involved, Eq. (29) simplifies to Here, b =

I2 −

3 2 I N 0: b+1 1

ð30Þ

1 Since both I1 and I2 are valued less than according to Eqs. (20), (21) indicates I12 b I32 . This proves that condition (30) is 3 guaranteed, meaning forward dynamics solution is stable in the absence of a payload. Without a further rigorous proof, we state that, in general, condition (19) is a lot easier to be satisfied than condition (19) and hence, Eq. (28) is most likely stable. In addition, a large the hub inertia improves the chance of satisfying condition (29). Remark 3.3. The natural frequency of Eq. (28) ωf represents the natural frequency of the slewing beam. It differs from the beam's 1 modal frequency ω1 due to the influence of rigid body motion. Comparing ωf with ωq, since Mf N 0, MqθM− θθ Mθq N 0, thus, Mf b Mqq, indicating ωf N ωq = ω1. Therefore, rigid body rotation stiffens the elastic beam. It is also worth noting that ωf ≠ ωi ≠ ω1. Given a C0 function of τ applied at the hub and beam initial conditions ξ(0) = ξ0 and ξ˙ ð0Þ = ξ˙ 0 , the solution of Eq. (28) is: ξðt Þ =

! Rf t ξ˙ 0 + ∫0 τðιÞ cos ωf ι dι sin ωf t + ωf ωf

ξ0 −

! Rf t ∫0 τðιÞ sin ωf ι dι cos ωf t: ωf

ð31Þ

To obtain the tip response, we again differentiate the kinematics Eq. (14) twice. Equation (13a) is used to solve for θ̈ with ξ̈ given by twice differentiating Eq. (31). The tip response is: −1 β̈ ðt Þ = Jq τðt Þ + Dq ξðt Þ;

ð32Þ

where

  −1 −1 −1 = 1− Ψq −Mθθ Mθq Mf Mqθ Mθθ ;
 −1 −1 Dq = − Ψq −Mθθ Mθq Mf Kq : −1

Jq

1 Here we show that the right hand side of J− is indeed the inverse of Jq. From Eq. (27), q


 2 −1 −1 Jq = Mθθ Mqq −Mθq Mi = Mi Mf Mθθ :

ð33Þ

From Eq. (32),
 −1 −1 Mθq Mqθ Mf−1 Mθθ Jq−1 = Mf −Ψq Mqθ + Mθθ −1 = Mi Mf−1 Mθθ
−1 = Jq :

ð34Þ

To examine the tracking quality of the inverse dynamics solution, we apply the torque found in Eqs. (27) to (32). Note that the resulting elastic motion is given by Eq. (31). Therefore, oscillations of both frequencies ωi and ωf are potentially present in the response to a torque aimed to achieve a prescribed tip motion: β̈ ðt Þ = Aβ g ðt Þ + Bβ sin ωi t + Cβ cos ωi t + Dβ sin ωf t + Eβ cos ωf t:

ð35Þ

For exact regeneration of the prescribed tip trajectory, one desires that Aβ = 1 and all other coefficients are zeros. To show the details of these coefficients, we take a sinusoid function as an example for the prescribed tip acceleration, g ðt Þ = h sin αt:

ð36Þ

854

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

Here h is a constant representing the maximum angular acceleration (rad/s2) of the tip with respect to the fixed frame and α is angular frequency. With Eq. (36), one can carry out integrations in Eq. (25), which defines the beam elastic motion in the inverse dynamics solution. It is qðt Þ = Ai g ðt Þ + Bi sin ωi t + Ci cos ωi t

ð37Þ

with Ri ; ω2i −α2   q˙ α Ri h Bi = 0 − ; ωi ω2i −α2 ωi

Ai =

C i = q0 : Obviously, α ≠ ωi is a requirement on the given function g(t). Consequently, substituting the above into Eq. (27) gives the control torque
 τðt Þ = Jq + Cq Ai g ðt Þ + Cq Bi sin ωi t + Cq Ci cos ωi t:

ð38Þ

Taking Eq. (38) as the input and carrying out integrations in Eq. (31), the response of the elastic coordinate should contain oscillations of three distinct frequencies, ξðt Þ = Af h sin αt + Bf sin ωi t + Cf cos ωi t + Df sin ωf t + Ef cos ωf t;

ð39Þ

where

Af = Bf = Cf =


 Rf Jq + Cq Ai ω2f −α2 Rf Cq Bi ω2f −ω2i Rf Cq Ci ω2f −ω2i

; ;

α ξ˙ Df = 0 − ωf ωf Ef = ξ0 −

;

!


 Rf h Jq + Cq Ai ω2f −α2

Rf Cq Ci ω2f −ω2i

−

ωi ωf

!

Rf Cq Bi ω2f −ω2i

;

:

It is also apparent that α ≠ ωf must be demanded of the input g(t) in addition to α ≠ ωi concluded earlier. Finally, coefficients in Eq. (35) for the tip motion are Aβ =

−1

Jq

Bβ = Cq

−1 Jq

Jq +

ω2f −α2 +

−1

Cβ = Cq q0 Jq

Dβ = Dq

!

Dq Rf

+

+

D q Rf

;

Rf Cq ω2f −ω2i

Rf Cq q0 ω2f −ω2i

;

ð40aÞ

ð40bÞ

!

ω2f −ω2i !

ω2i −α2

!

! α Ri h q˙ 0 − ; ωi ωi ω2i −α2

ω2f −ω2i

ξ˙ 0 q˙ − 0 ωf ωf

Eβ = Dq ξ0 −

!

Dq Rf

Cq Ri

ð40cÞ ! ;

ð40dÞ

! :

ð40eÞ

It is not difficult to show that the right hand side in Eq. (40a) can indeed be reduced to 1 and therefore, Aβ = 1 (see eq. (B3) in Dq Rf = 0 from Eq. (B5), which results in Bβ = Cβ = 0. Therefore, the fact Appendix B). In Appendix B, we also show that Jq−1 + 2 ωf −ω2i that the inverse dynamics solution results in beam oscillation of frequency ωi is not a problem in the resulting tip trajectory.

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

855

Finally, in order for the last two terms in the tip response to be zero, we can choose the initial conditions of the inverse dynamics solution according to 2

2

q0 =

ωf −ωi ξ ; and Rf Cq 0

ð41aÞ

q˙ 0 =

ω2f −ω2i ˙ ξ0 : Rf C q

ð41bÞ

It is interesting to note that the initial conditions for the inverse dynamics solution are required to be proportional to the respective actual initial conditions rather than being identical in order to produce exact tip tracking. Consequently, we show that the inverse dynamics solution results in a tip motion response β̈ ðt Þ = g ðt Þ, providing exact tracking. Finally, we mention that although a sinusoidal acceleration function (36) is used in the above illustration, using a step function results in the same conclusion, that is, β̈ ðt Þ = g ðt Þ with initial conditions chosen by Eq. (41). 4. Numerical simulations Numerical simulations of the inverse dynamics problem are carried out to verify the theoretical solutions. We use a common ODE solver “ode45” in Matlab to solve Eqs. (13) and (14). For the inverse dynamics problem, the inputs are prescribed tip motion and system initial conditions. Solutions for τ are saved as a data file and later read as the input for simulations on tracking performance. Table 2 lists system parameters used in the simulation. The desired tip motion is to rotate 90 degrees about the hub center in 1 second with its acceleration following a one-cycle sinusoidal function in 1 s and maintaining zeros afterwards. Using Eq. (36), the desired trajectory parameters are h = π2 and α = 2π. This trajectory is used for all simulation cases discussed below. First, we consider a payload with inertial properties of R = 3 and e = 0.2, meaning Me = 0.1651 Kg and re = 86 mm. Referring to Fig. 3, such a combination of R and e ensures the existence of a bounded torque in inverse dynamics. Besides qualitative assurance, frequency characteristics serve as a good quantitative evaluation. According to the theoretical solutions derived above, with R = 3 and e = 0.2, the smallest root of the frequency Eq. (6) is γ1L = 0.9615, translating to the first modal frequency at ω1 = 5.67 rad/s (or 0.903 Hz). The inverse dynamics torque should contain oscillations in ωi = 44.3 rad/s (or 7.05 Hz). The natural frequency of the driven beam is at ωf = 8.03 rad/s (or 1.28 Hz), slightly above ω1. Fig. (4a) shows the torque solution. The frequency element of ωi at 7 Hz can be identified over the entire time period. The frequency element of α at 1 Hz is also visible for the first second when g(t) is non-zero. Observations are consistent with what the analytical solution (35) suggests. Subsequently, the torque solution is used as the input to the forward dynamics simulation to evaluate tracking performance. Angular displacements of the tip (β) and the hub (θ) are compared in Fig. 4b. From the kinematics relation (14), hub rotation oscillates in similar fashions as the tip deflection (q) so that the tip can follow the prescribed trajectory exactly. Comparing Fig. 4b with c, we see that the hub angle contains the same frequency contents as the tip deflection. Persistent oscillations can be seen in all responses due to the absence of damping in the system. All three frequencies, namely α, ωi, and ωf, should be present in both θ and q responses. However, because of the dominance of the forced response over the natural response, ωf element is not visually obvious in either θ or q. Nevertheless, it is much easier to identify ωf in the tip tracking error plot from Fig. 4d. The first half of the response is in α while the second half is in ωf, both are close to 1 Hz. Tracking error is on the order of 10− 4 degrees limited by the numerical tolerance, indicating an exact tracking. Solution procedure as described in the previous section applies to a general case where n modes (n N 1) are included in the model. Using Eqs (13a) and (14), the elastic self-motion in the inverse dynamics problem is governed by n equations Mi q̈ + Kq q = −Mqθ g ðt Þ

ð42Þ

where Mi = Mqq −Mqθ Ψq . To obtain bounded solutions of elastic coordinates q, a necessary and sufficient condition is Mi N 0. This is a more strict condition than condition (19). In Fig. 3, the dashed line boundary curve applies to the assumed-modes model with n = 2. Any combinations of (R,e) above the dashed boundary curve lead to a bounded inverse dynamics solution. Obviously, the size of the

Table 2 System parameters. Properties

Notation

Value (Unit)

Bending stiffness Mass density Cross section area Length Hub moment of inertia

EI ρ A L Jh

0.1647 (N-m2) 8000 (Kg/m3) 16 (mm2) 430 (mm) 0.0339 (Kg-m2)

856

100

4

90

3

80 2 60

deg

torque (N-m)

70 1 0

50 40

-1

30 20 -3 -4

10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0

2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

1.8

2

time (s)

time (s)

a 1.2

0.2 0.15

b

x 10-4

1

0.8

error (deg)

tip deflection (m)

0.1 0.05 0 -0.05

0.6

0.4

-0.1 0.2 -0.15 -0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8

1

time (s)

time (s)

c

d

Fig. 4. Simulation results with R = 3, e = 0.2. (a) — Actuator torque; (b) — Angular displacement; (c) — Tip deflection; and (d) — Tracking error.

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

-2

100

1.5

90 1

80

θ β

70 60

deg

torque (N-m)

0.5

0

50 40

-0.5

30

10 -1.5

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0

2

0

0.2

0.4

0.6

0.8

time (s)

a

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

1.8

2

b

0.02

4.5

0.015

4

0.01

3.5

x 10-6

3

0.005

error (deg)

elastic cooridantes (m)

1

time (s)

0 -0.005

q2x103

-0.01

2 1.5 1

q1

-0.015

2.5

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

20

-1

0.5

tip deflection

-0.02

0

-0.025 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-0.5

0

0.2

0.4

0.6

0.8

1

time (s)

time (s)

c

d 857

Fig. 5. Simulation results with n = 2, R = 5, e = 0.9. (a) — Actuator torque; (b) — Angular displacement; (c) — Tip deflection; and (d) — Tracking error.

858

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

2.5

x 10-3

5

x 10-3

2

difference in torque (N-m)

4 1.5 3

1

pseudo-rigid

exact (x103)

0.5

2

0 1

-0.5 -1

0

-1.5 -1 -2 -2.5 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-2 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time (s)

a

b Fig. 6. Comparative results n = 1, R = 5, e = 0.9.

region is much smaller than that of the single-mode model. With increased number of elastic modes included in the model, the size of the solution region gets smaller, requiring larger amount of payload inertia. Although such demand may not be practical in some cases, systems with large payloads and long-slender arms are common in space manipulators. The second case of simulation is conducted on a system with n = 2. We consider R = 5 and e = 0.9. According to Fig. 3, a bounded inverse dynamics solution should be obtained. Fig. 5a shows the control torque. Comparing Figs. 5a with 4a, it seems counterintuitive that a less amount of torque is required to maneuver a heavier payload that is to follow the same trajectory. However, the reason is because a smaller level of vibration is involved in the latter case. Both rigid angle θ (Fig. 5b) and tip deflection (q1+q2) (Fig. 5c) show smaller magnitudes of oscillation compared with Fig. 4b and c, respectively. As a result, rigid angle deviates less from the tip angular displacement as indicated in Fig. 5b. Fig. 5c also shows elastic coordinates q1 and q2. Tip deflection overlaps with q1 because of the absolute dominance of q1 in the overall beam deflection. This can be noted in Fig. 5c by the label indicating that q2 is at least three orders of magnitude smaller than q1. Modal frequencies f1 = 1.69 Hz and f2 = 22.23 Hz are noticeable in Fig. 5c although the exact frequencies of elastic coordinates would have been modified by the rigid-body motion to be slightly higher than the corresponding modal frequencies. Finally, tracking error as shown in Fig. 5d indicates an exact tracking once again. Finally, we wish to show that a single-mode model is a preferred model when solving the inverse dynamics problem. To check how well a single-mode model performs, we use the single-mode model to determine the inverse dynamics torque. The torque is then applied to the assumed modes model with n = 2 (consider exact here) to evaluate tracking performance. We use the same values of R and e as used in the previous case (Fig. 5) for comparison purpose. Fig. 6 shows such comparison on torque (Fig. 6a) and tracking error (Fig. 6b). One can see that the difference in torque as the inverse dynamics solutions is not significant, on the order of one thousandth of the applied torque. From the high frequency contents in the torque difference in Fig. 6a, we conclude that the difference is due to the second mode vibration, although small in magnitude. Although the single-mode solution results in a tracking error that is roughly three orders of magnitude larger, it's still within the tolerance of practical sense. After all, it is the first mode of vibration that dominates the beam deflection. The advantage of adopting the single-mode model is that it gives not only a simpler and faster solution, but also the ease to avoid potential numerical stiffness when both low and high frequency dynamics are involved in the numerical integration. This is of particular interest to control applications.

5. Conclusions In this paper, the inverse dynamics problem of a slewing flexible beam is examined. In particular, a closed form causal solution is derived through the use of a single-mode beam deflection model. It is realized that the property of the deflection shape function determines whether or not the inverse dynamics problem has a solution. A necessary condition is followed which reveals the crucial role the payload rotary inertia can play in ensuring the existence of a bounded solution. Physically, the payload rotary inertia serves to constrain the beam rotation at the boundary point which alters the deflection shapes. This finding disagrees with the explanation that time delay between input and output is the cause of the non-existence of a causal solution. Through the analysis of inverse dynamics solutions and their tracking performance, it is realized that the inverse dynamics solution oscillates in a fundamental frequency different from the beam's natural frequency. They are carefully paired to provide a zero oscillation at the boundary point. Numerical simulations show a good tracking performance of the inverse dynamics solution

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

859

by the single-mode model. Such solution can be used in command shaping or inversion techniques for the tip tracking control of flexible manipulators. Acknowledgment This research is partially funded by the Natural Sciences and Engineering Research Council of Canada. Appendix A Due to the periodic solution of the beam deflection, beam equation can be written as

EI

d2 ″ 2 ϕi ðxÞ−ρAωi ϕi ðxÞ = 0: dx2

ðA1Þ

Multiplying the above equation with ϕj and integrate over the entire length gives L d2 ″ 2 ϕ ðxÞdx = ωi ρA∫ ϕj ðxÞϕi ðxÞdx: 2 i 0 dx Integrating by parts on the left hand side yields: L

EI∫ ϕj ðxÞ

ðA2Þ

0


 L L 2 EI ϕj ðLÞϕ‴i ðLÞ−ϕj′ðLÞϕi″ðLÞ + ∫ ϕj″ðxÞϕi″ðxÞdx = ωi ρA∫ ϕj ϕi dx: 0

0

ðA3Þ

Apply the boundary conditions as given by Eq. (4a) and note that ϕ(L) = 1. We obtain
 L L 2 EI∫ ϕj″ðxÞϕi″ðxÞdx = ωi ρA∫ ϕj ðxÞϕi ðxÞdx + Me + Je ϕj′ðLÞϕi′ðLÞ : 0

0

ðA4Þ

Switching the order of subscripts i with j in the above derivation, Eq. (A4) can also read:
 L L 2 EI∫ ϕi″ðxÞϕj″ðxÞdx = ωj ρA∫ ϕi ðxÞϕj ðxÞdx + Me + Je ϕi′ðLÞϕj′ðLÞ : 0

0

ðA5Þ

Subtracting Eq. (A5) from Eq. (A4) leads to
 L 2 2 ωi −ωj ρA∫ ϕi ðxÞϕj ðxÞdx = 0: 0

ðA6Þ

Since ωi and ωj are distinct, we obtain the orthogonality relation L

∫ ϕi ðxÞϕj ðxÞdx = 0 0

∀i≠j:

ðA7Þ

Appendix B We first show that Jq +

Cq Ri ω2i

ðB1Þ

= Mθθ

by expressing out the left hand side using Eqs. (25), (28), and (35):

Jq +

Cq Ri ω2i

=

−1 Mi Mf Mθθ

+


 Mi−1 Mθq −Mθθ ψq Kq ⋅Mi−1 Mθq

Mi−1 Kq

  −1 = Mi Mf Mθθ + Mθq −Mθθ ψq Mθq :

Using Eq. (29) to replace Mf, the above expression reduces to Mθθ, which proves Eq. (B1).

860

Q. Sun / Mechanism and Machine Theory 46 (2011) 845–860

Similarly, by using Eqs. (29), (33), and (34), we can show that −1

Jq

+

Dq Rf

−1

= Mθθ :

ω2f

ðB2Þ

Assuming α2 ≠ ω2f and α2 ≠ ω2i , using Eqs. (B1) and (B2), Jq +

Cq Ri

! Jq−1

ω2i −α2

+

Dq Rf

!

2

=

ω2f −α2

=

ω2i −α2

−1

2

Mθθ ωi −Jq α

⋅

2

−1

2

Mθθ ωf −Jq α ω2f −α2

−1 ω2i ω2f −Mθθ Jq−1 ω2i α2 −Mθθ Jq ω2f α2 + α4
 : 2

ωi −α2 ω2f −α2

1 2 2 −1 2 2 Using Eqs. (34) and (35), we can show that Mθθ J− q ωi = ωf , and ,Mθθ Jqωf = ωi the right hand side of the above expression reduces to one. Therefore, we prove that

Jq +

Cq Ri

! −1

Jq

ω2i −α2

+

Dq Rf ω2f −α2

! = 1:

ðB3Þ

1 2 1 1 −1 2 −2 = ω2f ω− . From Eq. (34), MiM− = J− and By definition from Eqs. (25) and (29), MiM− f i f q Mθθ. Therefore, Jq Mθθ = ωf ωi equivalently,

−1

2

−1

2

Jq ωi = Mθθ ωf :

ðB4Þ

1 2 −1 2 Substituting Eq. (B2) into the right hand side leads to J− q ωi = Jq ωf + DqRf. Therefore,

−1

Jq

+

Dq Rf ω2f −ω2i

= 0:

ðB5Þ

References [1] E. Bayo, P. Papadopoulos, J. Stubbe, M. Serna, Inverse dynamics and kinematics of multi-link elastic robots: an iterative frequency domain approach, International Journal of Robotics Research 8 (6) (1989) 49–62. [2] S. Devasia, E. Bayo, Inverse dynamics of articulated flexible structures: simultaneous trajectory tracking and vibration reduction, Journal of Dynamics and Control 4 (1994) 299–309. [3] D. Kwon, W. Book, A time-domain inverse dynamic tracking control of a single-link flexible manipulator, ASME Journal of Dynamic Systems, Measurement, and Control 116 (2) (1994) 193–200. [4] S. Taylor, S. Yau, Boundary control of a rotating Timoshenko beam, ANZIAM Journal 44 (E) (2003) E143–E184. [5] A. De Luca, P. Lucibello, G. Ulivi, Inversion techniques for trajectory control of flexible robot arms, Journal of Robotic Systems 6 (1989) 325–344. [6] C. Damaren, Approximate inverse dynamics and passive feedback for flexible manipulator with large payloads, IEEE Transactions on Robotics and Automation 12 (1) (1996) 131–138. [7] D. Wang, M. Vidyasagar, Transfer functions for a single flexible link, International Journal of Robotics Research 10 (5) (1991) 540–549. [8] J. Park, H. Asada, Dynamic analysis of noncollocated flexible arms and design of torque transmission mechanisms, ASME Journal of Dynamic Systems, Measurement, and Control 116 (2) (1994) 201–211. [9] Q. Sun, Control of flexible-link multiple manipulators, ASME Journal of Dynamic Systems, Measurement, and Control 124 (1) (2002) 64–75. [10] S.S. Rao, Vibration of Continuous Systems, John Wiley & Sons Inc, New Jersey, 2007, pp. 325–341. [11] E. Barbieri, U. Uzguner, Unconstrained and constrained mode expansions for a flexible slewing link, ASME Journal of Dynamic Systems, Measurement, and Control 110 (2) (1988) 416–421.