- Email: [email protected]
force control of flexible manipulators
Mechatronics Vol. 2, No. 2, pp. 129-148, 1992 Printed in Great Britain
0957-4158/92 $5.00+0.00 © 1992 Pergamon Press plc
COMPREHENSIVE DYNAMIC MODELING AND MOTION/FORCE CONTROL OF FLEXIBLE MANIPULATORS YUEH-JAw LIN and TIAN-SOoN LEE Department of Mechanical Enginccring, The University of Akron, Akron, OH 44325, U.S.A.
(Received 16 September 1991; accepted 30 October 1991) Abstract--This paper presents a comprehensive dynamic modeling of a flexible link manipulator by using Hamilton's principle. The model is formulated by considering the flexible arm as a Timoshenko beam model. Hence, the shear, bending as well as rotational inertia effect of the manipulator are all taken into account in the dynamic formulation. In order to obtain a comprehensive dynamic model, the gravitational effect, which is missing in most of the flexible manipulator dynamic models, is also included in the formulation. Then, an efficient motion and force controller design method utilizing root contour analysis is proposed. With the proposed method the multi controller gains of the position/force controller for the flexible manipulator can be chosen analytically, as opposed to the method of trial and error used conventionally for selecting manipulator controller gains. The tuning process of these controllers' gains proves the simplicity of gains selection. And the motion/force simulation results verify the effectiveness of the controllers. The simulation results also show that the derived control scheme, which is based on a linear control system analysis, can drive the highly non-linear flexible arm to achieve the desired position and force satisfactorily.
INTRODUCTION T h e main objective of this w o r k is to advance the state-of-the-art in the area of m o t i o n and force control of flexible a r m manipulators. T h e study is based on a flexible single-link m a n i p u l a t o r o p e r a t e d in a vertical plane as shown in Fig. 1. T h e driving t o r q u e is p r o v i d e d by an a c t u a t o r m o u n t e d on the m a n i p u l a t o r base joint. Position, velocity and force sensors are used to m e a s u r e joint angular position, angular velocity as well as wrist position, velocity and contact force. The study in this w o r k is aimed at developing an effective and systematic m e t h o d o l o g y to accurately control the m o t i o n o f flexible manipulators subjected to external contact forces. C u r r e n t research on flexible link manipulators has the following problems, namely: (1) most flexible m a n i p u l a t o r models u n d e r study are either considering the manipulator o p e r a t e d on a horizontal plane or the shear and gravitational effects are excluded in the d y n a m i c model, which m a k e s it difficult to control the m a n i p u l a t o r m o t i o n accurately; (2) the m o t i o n of the flexible manipulators being considered is usually a contact-free operation, i.e. without any contact forces, which m a k e s it impractical when p e r f o r m i n g industrial assembly tasks, and (3) the selection of controller gains is usually by empirical m e t h o d s or trial-and-error approaches, which m a k e s the controller design very time consuming. In view o f the a f o r e m e n t i o n e d p r o b l e m s the flexible link manipulator, u n d e r study in this paper, is constructed in such a way that it is o p e r a t e d vertically and in 129
130
Y.-J. LIN and T.-S. LEE (F.)
Y0
~+
T X0
Fig. 1. Configuration of a single link flexible arm manipulator.
addition, subjected to a contact force at the tip. The dynamic model is comprehensively formulated by considering shear, bending, rotational inertial and gravitational effects. The inclusion of gravitational effect is considered to be necessary since most industrial robots are operated in three-dimensional workspace. In the derivation of control algorithm of the flexible link manipulator, this paper also deals with external contact forces which cannot be avoided when performing assembly tasks. However, the research in this area is missing in the literature of flexible link manipulators as far as the authors are aware. Moreover, in this paper we propose a quick and systematic way of selecting good controller gains utilizing the root contour method. This method greatly reduces the laborious trial-and-error approach conventionally used in gains selection. Although the research in force control of flexible manipulators is rare, many research articles studying the force control problems of rigid link manipulators, which may be valuable to this study, can be found in the literature [1-9]. Whitney [1] in his early work proposed a feedback force control algorithm called accommodation control which modeled the end-point impedance as a dashpot. The controller is not considered accurate and is vulnerable to stability problem because it does not compensate for the manipulator dynamics. Hogan et al. [2, 3] developed a force control scheme which they called impedance control. This control method is evolved from observation of behavior of admittance and impedance subjected to an input. According to their definition, when a manipulator is interacting with an environment, the manipulator should assume the behavior of an impedance while the environment is an admittance. Another prominent force control method is the active stiffness control [4]. It started by defining a linear dynamic relationship between the contact
Comprehensive dynamic modeling of flexible manipulators
131
force and the deviation of the hand position from its desired position. This relationship in turn is transformed into joint space for implementation. In order to control position and force simultaneously, the hybrid position/force control technique [5] is developed. Since the controller of this method only consists of error driven servo, it is expected that the controller will not perform its tasks accurately. Recognizing the limitation of the hybrid position/force algorithm, Shin and Lee [6] have derived an extended resolved acceleration control scheme [7] which takes the manipulator dynamics into account. However, this method is devised for position control rather than for force control. Hence, to control position and force simultaneously, the extended resolved control method is proposed [8, 9]. The difference between the resolved acceleration and extended resolved acceleration control method is that in the latter, in addition to the position error, the force error term is included in the error driven servo. Since the controller implements its tasks in joint space, while the desired position and force are specified with respect to an external frame, considerable computational time will be consumed in calculating inverse Jacobian and homogeneous transformation matrices. As mentioned previously, the main objective of this paper is to develop an efficient controller which simultaneously deals with the flexible link manipulator's position and contact force and takes the manipulator dynamics into consideration to assure accuracy. A quick and systematic approach for tuning controller gains using the root contour method is proposed [10, 11]. Basically, the controller is developed based on the so-called partitional control law architecture [12]. However, it is extended in this paper to include both force and position errors driven servos. The artificial and natural constraint concept is used to specify and identify force and position control directions. This control scheme is different from others [1-9], which are used to control position and force of rigid link manipulators, in the sense that the Cartesian space position and force are indirectly controlled by joint space position and force. In addition, the dynamics of gravitational effect, flexible effect, and contact force are included in this controller design to reduce their influences on the accuracy of the motion of the flexible arm manipulator. In the following sections the detailed derivation of the control scheme based on the proposed methodology is described. Since the shear force depends on the shape of the cross-section of the arm, the numerical factors of three typical manipulator cross-sectional shapes, namely, hollow round, rectangular, and square, are derived in the next section. Following that, a complete set of dynamic equations of the flexible single-link manipulator with an endmass are formulated using the Hamilton's principle. The arm is considered to be a Timoshenko beam model when deriving the dynamic equations. The inclusion of shear force effect is required only in the cases of flexible robots with medium slenderness ratios subjected to very high speeds and with high precision requirements, as was addressed in [13]. For other cases, the shear force effect can usually be neglected for simplicity. The derived dynamic equations are solved by using the Galerkin's method. In Section 4, the control law that is used to control the motion of the flexible link manipulator subjected to external contact forces is developed. In addition, an efficient and systematic approach for tuning the controller gains is proposed. Finally, the motion simulation results of the flexible single link manipulator and the most
132
Y.-J. LIN and T.-S. LEE
important findings with the simulation are reported to verify the effectiveness of the proposed methodology.
N U M E R I C A L FACTORS As mentioned in the previous section, the upon the shape of the cross-section of the called numerical factors which represent the be needed in the dynamic formulation. Referring to [14], the numerical factor of a
shear force of the flexible arm depends arm. Therefore, the physical quantities characteristics of arm cross-sections, will flexible b e a m is defined as
Ib AQ
k' --
(1)
where I is defined as the second area m o m e n t of inertia of the cross-sectional shape of the b e a m computed with respect to its neutral axis. Q denotes the first m o m e n t about the neutral axis of the area contained between an edge of the cross-section of the b e a m parallel to the neutral axis and the surface at which the shear stress is to be computed. A and b represent the cross-sectional area and the width of the cross-sectional area at which the shear stress is required, respectively. Due to the fact that the numerical values of k' for most flexible beams are less than unity, there is a need to define the so-called correction factors k, which are inversely proportional to their corresponding numerical factors k' ]14]. Hence, k -
1
k'
-
AQ Ib
(2)
In this section we will derive the correction factors k for three typical robot arm cross-sectional shapes, namely, hollow round, hollow rectangular, and hollow square. Let us start with considering a b e a m with a hollow round cross-sectional shape. The inner and outer radius of this shape are given as r I and r 2, respectively. Thus, the approximate k for this shape is
4 (r~ + r2r 1 + r~) k
3
(r~ + r~)
(3)
for the thick type. However, for the thin type of this shape, the approximate k is obtained to be k = 2. Similarly, for a b e a m with a hollow rectangular cross-sectional shape, the height and width of the inner rectangle are given by h 1 and bl, respectively. The height of the outer rectangle is denoted by h2 and the width is represented by b2. The approximate ks for the thick and thin type of this shape are obtained as:
3 h l ( b 2 h ~ - b~h~) k - - ~ ( b 2 h ~ - blh 3) '
(4)
and k = 1, respectively. Notice that in obtaining these two ks the web area is used instead of the cross-sectional area. The hollow square cross-sectional b e a m means that the height and the width of the cross-section of the b e a m are equal. Thus, by setting
Comprehensive dynamic modeling of flexible manipulators
133
bl = h i and b 2 = h2, the correction factor k of the hollow square is therefore given by 3 hl(h~ + hlh2 + h~) k - - ~ (h2 + hl)(h~ + h12) "
(5)
If the web and the flange of the hollow square are thin, then it implies that ht ~- h 2 . Thus, k = 1.125. The values of the correction factors k of the hollow round, rectangular and square shapes are tabulated in Table 1.
COMPREHENSIVE
DYNAMIC
MODELING
Figure 1 shows a flexible arm with a concentrated mass at the tip of the arm. It also shows that the arm is subjected to an external force F e. On the one hand, the concentrated mass can be considered as the mass of a tool which is carried by the arm to p e r f o r m a constrained task. On the other hand, the external force can be considered as the contact force generated by interacting the arm with its environment. In the two-dimensional case, the contact force only consists of x and y components. The other end of the arm is clamped into the m o t o r shaft. The figure also shows that the length of the arm is L, the mass per unit length of the arm is p ( x ) , the mass m o m e n t of inertia of the actuator is Jh and the mass m o m e n t of inertia per unit length of the arm is Ja" The X o Y o Z o coordinates are chosen to denote the inertia frame of the system where the directions of X0-, Y0- and Z0-axis are the X0-axis points to the right, Y0-axis is normal to the X0-axis and Z0-axis is determined by completion of the right-hand rule. A rotational frame X ~ Y I Z 1 is attached and rotated simultaneously with the flexible arm. This rotational frame is used to describe the deflection of the arm. D ( x , t) denotes the transverse deflection of the arm from its equilibrium state at distance x along the Xl-axis and at time t. The assumptions of the dynamic derivation are as follows. (1) (2) (3) (4) (5)
Axial elongation, coriolis and torsional effects are assumed to be negligible. The frictional force and backlash of the m o t o r are ignored in modeling. The motion of the arm is in a vertical plane. The arm is considered to be a T i m o s h e n k o beam. The arm is considered to be uniform and prismatic. Table 1. Correction factor k for hollow round, rectangular and square Type
Thick
Thin
Hollow round
4 (r~ + r2rl + r~) 3 (r~ + r~)
2
Hollow rectangular
3 h~(b2h~- b,h~) 2 (b2h~- blh~)
1
Hollow square
3 h~(h~ + hlh2 + h~) (h2 + h,)(h~ + h~)
1.125
134
Y.-J. LIN and T.-S. LEE
Hamilton's method is used to derive dynamic equations of motion of the flexible manipulator. For a non-conservative system, Hamilton's principle is stated as
t2 f,
i
((~(T e - ge) -t-- a W e ) d t = 0,
(6)
where Te = kinetic energy of the system, V e = potential energy of the system, We = work done by the non-conservative forces to the system. Figure 2 illustrates the displacement of the arm from its equilibrium state at location x and the time t after the joint has rotated 0 degrees. In view of this figure, the kinematic relationships of any point on the flexible link described in terms of the inertial frame is given by, Y =xsin0+
Dcos0,
(7)
X = xcos0-
Dsin0.
(8)
The kinetic energy of the flexible arm system is the combination of the kinetic energy of the flexible arm, endmass and actuator. In view of these, the total kinetic energy of the manipulator T e expressed in terms of the inertial frame is therefore given by,
1 "2 + re= 1 f~ p(?2+22)dx+ m(X
?e) + 21 ft~
Ja~2dx + 1 fL
aO2dx +
] "2
(9) where the first term is the translational kinetic energy of the flexible arm. The second term is the translational kinetic energy of the endmass. The third term is the rotational kinetic energy due to the angular velocity caused by the bending moment. The fourth term is the rotational kinetic energy due to joint angular velocity. The last term is the kinetic energy of the joint itself. For simplicity, in Eqn (9) we neglect shear force effect and the Coriolis effect. The potential energy of the arm is caused by the internal bending moment and shear force as well as the gravity force. Therefore, the total potential energy of the system is
Ve - 2EI1
MZdx + ~
Vydx + mgY + pgLY,
(10)
where the first two terms are due to the bending moment and shear force, respectively. The last two terms are due to the gravitational effects on the tip mass and the arm itself. For the system shown in Fig. 1, the non-conservative forces are the input torque T and the external forces [F~]. In joint space, the work done by these non-conservative forces is
We = T O - [Jc]r[Fe]O,
(11)
where [Je] r denotes the transpose of the Jacobian matrix. For the two-dimensional
Comprehensive dynamic modeling of flexible manipulators
135
Yo
y, {F.}
X
L
xcosO
Fig. 2. Deflection of a point on a flexible arm described in terms of inertia frame.
case, the Jacobian is defined as
-Lsin(@+ O) [Jc]
(12)
=
However, the angle of deflection o~(x, t) at location x along the arm and time t is equal to the linear combination of angle of rotation due to the bending moment ~(x, t) and angle of distortion due to the shear y(x, t). It is measured from the equilibrium state and is given by aD cv(x, t) - ax - ~p(x, t) + y(x, t).
(13)
The bending moment M and shear force V are related to the ~p(x,t) and y(x, t) [14] by the following relationships:
a~
M ( x , t) = E1 5x
and
V(x, t) = k ' G A y ,
where E G A k' I
= = = =
Young's modulus of the material of the arm, shear modulus of the material of the arm, cross-sectional area of the arm, numerical factor, and second area moment of inertia of the arm.
(14)
136
Y.-J. LIN and T.-S. LEE
From the Hamilton's principle, we obtain dynamic equations of the arm with four boundary conditions as follows: -o
+
ato ! +
'CA
-
= 0,
(15)
and
-J~ 32~p + E132--~ + k'GA 3t 2
-
~X 2
x=0,
D =0,
x=0,
~p=0, a~p
x = L,
E I ~ x =0,
=
-~
+m
~p = 0,
(16)
~2D]
Ld20
+ ~-t2 ] = 0,
(17)
as well as dynamic equation of motion of the joint L [x2d20 32D]dx ( d20 ~2 D ] fo - P ~ dt 2 + x 3 t 2 ] m L2 - dt - 2 + L ~t 2 x=L]
- ( m g L c o s O - rngsinOD[~_L._
d20
d20
- Ja L - -
Jh--
dt 2
dt 2
"}-
pgL2cos02
-- pgLsinODIx=L/2
+
[Jc]T[Fe])
+ T = 0.
(18)
Notice that the dynamic equations of the arm and the joint, Equations (15), (16) and (18), are coupled equations. However, Eqns (15) and (16) can be combined into one equation by using the calculus concept that a tangent of a curve is equal to the first spatial derivative of the function of that curve. Thus, one obtains the relationship of 3DB ~P- 3 x ' (19) where DB indicates the deflection due to the bending moment, Substituting Eqn (19) into Eqns (15) and (16), the deflection of the arm can be written in terms of DB as
--
[ D2D13 1
D = Dn + k'GA
Ja
~t 2
32DB] E1
(20)
~X 2 ]'
and ~4DB - -
E1 3x 4
S2DB +
p - - -
3t 2
34 DB J.
-
-
~X29t 2
pEl ~4 D13 k'GA ~x2~t 2
+
PJa ~4DB k'GA ~ t 4
px
d20 dt 2
(21) Equation (21) is the dynamic equation of motion of the arm subjected to the forcing function - p x d 20/dt 2.
Comprehensive dynamic modeling of flexible manipulators
137
Similarly, the boundary conditions of Eqn (21) will be transformed into the following forms 3DB x=0, =0 ax
x = O, x-- L, x=
L,
1 ja 32DB - at - -2 E , DB + k'G-----A
~2DB] a-Tr-j = o
~2Dl~
E1 - -0 ax 2 53D B Ja ~t2a x
m[ a4O.
33D B E1 ax----T- + -7 J" %xZSt 2
a4 .l E1
9x 4 j = 0.
(22)
Because the closed form solution of Eqn (21) is difficult to obtain, the approximate solution will be sought. Galerkin's method [15], which is a two-step numerical method, is used to find the approximate solution. At first, it assumes that the solution of the dynamic equations of motion of this continuous system can be approximated by the linear summation of the solution of each discretized system. Each of these discretized system solutions consists of the spatial-dependent comparison function ~i(x) multiplied by the corresponding time-dependent generalized coordinate qj(t). n
Du = ~
(&(x)qj(t).
j=l
(23)
Then, it assumes that the difference between the approximate and the exact solution is small such that the accumulation of these weighted errors along the length of the beam is considered to be zero.
L
fo ~(X, t)q)r(x)dx = O,
r = 1, 2, 3 . . . . .
n.
(24)
in which r is a dummy variable and (Pr(X) are called weighted functions. Through a laborious calculation, the comparison functions are obtained as
/3,x
. I~jx
_~_
QliSm--L
~bj(x) = cos T -
- Q2jcosh , _ .
+ Qljsinh
/3jx
L '
(25)
where /3j = jzr j =
Q lj 1 Q2j =
1,2,3 ..... =
Clj C2j
(26)
,
(27)
k ' G A + k'G---~ 2 k'GA
Clj = E I
n,
(28)
k'GA
sinh/3j - sinflj) +
+ Ja)~}(Q2jsinhflj+sinflj)+
,
0
mE-~I(~)3(Q2ciosh[3-jcos3ip)j.
1
(Q2jcoshfij + cosflj)
(29)
138
Y.-J. LIN and T.-S. LEE C2j = E I
(cosh/3j + cos/3j) + - -
( sinh/3i + sin/3j)
p
+ JaZ~(cosh/3/
-
c o s / ~ , ) + r -e E- l P
"--'[~ZI3( sinh/3j \c!
-
(30)
sin/3j).
However, following the standard procedures of Galerkin's method, the forward dynamics of the generalized coordinates and joint coordinates are given by [ql = ([M] + P [ S I [ U I T ) - I [ ( [ K ] + p[S][H]r)[q]] + ~
P ([M] + p [ S ] [ U ] r ) - I ( T
- R 2 c o s 0 - [Jc]r[Fe])[Sl,
(31) and 0 = ( T - e 2 c o s 0 - [Jc]T[Fe]) + ( [ H ] r - [ U ] r [ M ] - a [ K ] ) [ q ]
(1 +
RI(1 + p [ U ] T [ M ] - I [ S ] )
p[CqT[Ml-l[S])
'
(32)
where the parameters in these two equations are as follows: K,j = [K] = E1
k' G A
~1L
d4~] dx,
~ q)r ~
dx - p
S,.= [S] = I,Lxq),.dx, L 3
R1 = ~ -
(33)
¢p,. dx 4
+ m L 2 + Ja L + Jh,
pgL 2 R~ = m g L + - - 2
I q)'qJj d x ,
(34)
(35)
(36) (37)
Notice that in the derivation of forward dynamics, the fourth order time derivative term is ignored. This is due to the fact that the contribution of this term to the final solution is small. The inverse dynamics can also be obtained from Eqns (31) and (32), which will be used for designing the controller.
CONTROLLER
DESIGN
The main problem facing us here is how to design a controller that can apply a desired force to the environment while simultaneously attaining a desired position at the influence of the flexibility of the manipulators. Since the generalized coordinates used in Eqn (31) are unmeasurable in a physical system, the control law will be derived solely from the inverse dynamics equation deduced from Eqn (32). In order to ease the derivation of the control law, the inverse dynamics equation based on Eqn (32) will be expressed in the compact form.
Comprehensive dynamic modeling of flexible manipulators
139
Introducing the notations W, Cg and f into the inverse dynamics equation from Eqn (31) yields W O + Cg + f + [Jc]T[Fe] = T, (38) where W = RI(1 + p [ U ] r [ M ] - I [ S ] ) , (39) Cg = R2 cos 0 f = -RI([H]
T -
[u]T[M]-~[K])[q].
(40)
(41)
In view of Eqn (38), the notations W can be interpreted as the inertia term, Cg as the gravitational influence term, and f as the flexible effect term. The control law conceived here has two parts. The first part of this control law is called a model base. It is used to linearize the non-linearity of the system and to reduce the system to an unit mass system. The other part of the control law is called the servo part. This portion of the control law is error driven. Since the model base has reduced the system to an unit mass system, the servo portion is designed to regulate this unit mass system. In view of Eqn (38), the control law has the following form: T = r/Uf + /3,
(42)
where r / a n d / 3 are the model bases. In view of Eqns (38) and (42), it is clearly shown that in order to make this system a unit mass system, r/and/3 should be chosen as r/=W, /3 = Cg + f + [Jc]T[Fe].
(43)
Thus, the dynamic equation of the unit mass system is given by Uf = 0.
(44)
The servo part is chosen such that the position and force errors of this unit mass system can be regulated. If the servo part is chosen as a linear combination of PD-controller operating on the position and velocity errors and P-controller operating on the force error, then we obtain Uf = kp(O d - O) + k,,(O d - O) + kfp[Jc]T([Fe] d - [Fe] ) + 0 a,
(45)
where kp and k,, are the proportional and differential gain constants for the position error and velocity error, respectively, kfp is the proportional gain constant for the force error. Combining Eqns (43) and (45) with Eqn (42), the closed-loop control law is therefore given by T = W[O d + kv(O d - O) + kp(O d - O) + kfp[Jc]T([Fe] d - [Fe])] + Cg + f + [Jc]T[Fe].
(46)
Equation (46) shows that the controller indirectly controls the Cartesian space position and force by joint space position and force. It also shows that the dynamics of gravitational effect, flexible effect and the contact force are compensated for to reduce their influences on the accuracy of the controller.
140
Y.-J. LIN and T.-S. LEE
Depending on the physical condition of an environment in which a robot is working, the contact force can be modeled as a damper, spring, inertia, or any combination of these three models. It is assumed that the robots are doing assembly tasks. Thus, the hand of the robot will come in contact with the workpieces. It is also assumed that the workpieces are composed of hard material. Hence, the contact force can be modeled as a high stiffness linear spring model. In joint space this contact force is expressed as (
[Fo] -- [k0ql[Jc] A0 +
AD'x=L)
L
(47)
'
where A0 = ( 0 - 00) and ADI,=L/L = (Dk=L -- Do]x=L)/L. In the two-dimensional case, the equivalent stiffness matrix
[keq ] is given by
in which the subscripts x and y refer to the x and y directions in the end-effector frame. An element in the equivalent stiffness matrix is given by 1
keqi =
1 +__+__
,
(49)
kei ksi kbi where ke~ and ks~ are referred to as environment stiffness and sensor stiffness in the direction i, respectively. Similarly, kbi is referred to as the flexible arm stiffness in direction i. An element of the equivalent stiffness matrix [keq ] is equal to zero when in that particular direction the motion of the manipulator is unconstrained. Hence, the control law becomes 7 = w[O ~ +
k,,(0 ~ - 0) + k~(0 ~ - 0)]
+ W[kyp[Jc]r([Fe] d -
[Keq][Jc]( A O + ADL[X=L))]
+ Cg + f + [Jc]r[Keq][Jc] ( AO + ADIx=L) L "
(50)
The block diagram of this closed-loop system is given in Fig. 3. Generally, the force and position controller gains can be obtained from the characteristic equation of the closed-loop system, or they can be obtained by tuning them experimentally. The characteristic equation of our system is obtained by combining Eqn (50) with Eqn (38) first. It is then linearized at certain [Fe] a and 0 d and expressed in s-domain as follows:
s 2 + k,,s + (kp + kfp[Jc]r[keq][Jc]) = 0.
(51)
The controller gains appear in the above equation require to be tuned properly to result in suitable control actions. It is proposed here that the controller gains be tuned by using the root contour method. It can be shown that the rules to construct the
Comprehensive dynamic modeling of flexible manipulators
141
model based
servO J
Delm.t "--/-- + eO
(
~
Equations (3~) and (32)
I
IF I I Fig. 3. Closed-loop control for a single link flexible arm manipulator using extended partitional control law.
locus of a single variable gain system are also applicable to the construction of the locus of a multiple-gains system. The first step of root contour method involves setting two of the gains in Eqn (51) equal to zero. Let us set k,, and kfp be zero. Then, Eqn (51) becomes s 2 + k s = 0.
(52)
The root locus of this equation can be obtained by dividing both sides of the equation by s 2. Therefore, the equation becomes
ks
1 +-7=0.
(53)
The root locus of Eqn (53) can be constructed through the use of the open-loop transfer function, ks/s 2. If the k s is selected a s k s 1 , then the next step is to restore one of these two gains, kfs, or k,,. In the case of restoring the gain k];., Eqn (51) becomes S 2 "~ kpl + kfp[JclT[keq][Jc] = O.
(54)
Dividing both sides of Eqn (54) by s 2 + kp~ results in 1 +
kfp[ Jc] T [ keql[ Jc] = o. S 2 + kpl
(55)
Again the root locus can be drawn from this equation. Similarly, if be kfpl, then Eqn (51) can be written as
1 +
kvs s~ + (/~si + krpl[&]T[/%q][4])
= 0.
kys is chosen to (56)
The root locus can be plotted using this equation. In addition, k~ is selected such that the desired behavior of the manipulator is attained.
142
Y.-J. LIN and T.-S. LEE C O N T R O L L E D M O T I O N / F O R C E SIMULATION
To test the effectiveness of the control law derived in the previous section, a motion simulation program has been developed to simulate the dynamic behavior of a single link flexible manipulator. In the simulation we consider the manipulator to be a Timoshenko beam model which is made of aluminum. The cross-sectional shape of the manipulator arm is assumed to be thin hollow round. Thus, from Table 1, the correction factor is 2. The gravitational constant used in this simulation is 9.81 m s -2. In addition, the mass moment of inertia of the actuator and the tip mass are given as 0.1 kgm 2 and 0.03 kg, respectively. The physical parameters as well as mechanical properties of the link used are tabulated in Table 2. Now, the first two modes of the approximate solution are used in the simulation to model the dynamic behavior of the link. This approximation is justified by the fact that the natural frequencies of the arm beyond second mode are laid outside the bandwidth of the actuator. The sampling time of the simulation is taken to be 0.01 s. The tasks to be performed by the controllers are maneuvering the manipulator's arm to the position of 20 ° and at the same time applying a step change force with magnitude 10 N to the workpiece. The desired direction of the step change force is in the y-direction of the hand space. In other words, the x-direction of the hand space is unconstrained. Since the workpiece is assumed to be hard, the stiffness of the workpiece and that of the sensor should be much higher than that of the flexible arm. Thus, the stiffness of the workpiece and the sensor are assumed to be [kw] =
[0°
0 ]
1.5 x 104
N m -1,
(57)
1.1 x 103
N m -1,
(58)
and [ks] =
respectively. However, the stiffness of the flexible arm can be calculated from
[0°
[kb] =
3E!
N m -t.
L 3~
Table 2. The physical and mechanical properties of the link Parameter
Value
Link length Outer radius Inner radius Flexural stiffness (El) Mass moment of inertia Density Modulus of elasticity (E) Modulus of rigidity (G)
2m 5.6 mm 4 mm 40.5651 Nm 2 0.1002 kg m2 0.2712 × 104 kgm2 7.1 × 10mNm -2 2.62 x 10mNm <
(59)
Comprehensive dynamic modeling of flexible manipulators
143
In view of these three stiffness matrices, the equivalent stiffness matrix is given by [keq] ~--
~
0-] Nm-1.
(60)
15
In order to achieve the specified goals, the controller gains should be chosen carefully. The root contour method is used to design the controller gains and through the assistance of the CC software the controller gains will be tuned accordingly. With this approach the locus of the controller gains will be traced out by varying one gain at a time. Then the gains are selected through the analysis of the loci. Figure 4 shows the root locus of varying the position proportional gain kp. Figure 5 shows the locus of varying kyp with a fixed kp = 75. In addition, Figs 6 and 7 indicate the loci of varying ko for two different kp values. As shown in Fig. 4, the locus of varying kp is on the imaginary axis of the s-plane. This means that the larger kp is chosen, the more oscillatory the response of the system will become. Therefore, kp should be assumed to be a reasonably small value to avoid oscillations. In our simulation, kp is chosen to be 75. As illustrated earlier, the root of the characteristic equation of kp is the open-loop pole of the locus of kfp. Thus, as seen in Fig. 5, kfp has the effect of adding oscillation to the system. In other words, kfp should be chosen as small as possible. In our case, we select kfp to be 0.045. In order to suppress the oscillatory effect of kp and kfp as well as the flexible mode of the arm while simultaneously to achieve the specified objectives, k,, is chosen to be 10 which is corresponding to the underdamped case (see Fig. 6). Thus, it can be expected that there will be an overshoot in the system response and that the response is very fast. Indeed, the real time simulations reveal that the overshoots of the position and applied force are about 12.5% and 10% from their desired values, respectively (see Figs 8 and 9). However, the rise time of this underdamped system is very short (~ 0.35 s). On the other hand, when the value of k v = 18.11 is chosen, Figs 8 and 9 show that 13 1197 531-
i
kp = 7 5
-1
-9--
] "~'~"-- kp - 7,.q
-11 -13
I ",60
,
I ",2g
I ,20
I
I ,68
x 1 Real s
Fig. 4. Root locus for varying
kp.
144
Y.-J. LIN and T.-S. LEE 16 12-
t
k #p -
0.045
I
4-
I,.. ~" k/p
-12 -
-lt2
I -11
i
~
0.0,t~
I
m 11
O
I
21 Real s
Fig. 5. R o o t locus for varying
16I 12
kt), with
k s, = 75.
kv~ 10
,e,
,-.
04 1
~=
//\ -4
-0!
-
= 25
~ i i -'----
k,= 1o
-12 -10
\\ __
k . = i8.11 " ~ k .
i
I
-11
I
I
I
i
10
i
20 Real s
Fig. 6. R o o t locus for varying k,, with
kp
= 75 and
kip =
0.045.
the overshoots of the system are reduced to zero and the position and force responses of the system are slowing down. In fact these two phenomena correspond to the system that behaves critically. These phenomena are due to the fact that by increasing the value of kv the damping effect of the system will increase. Referring back to Fig. 6, the locations of the closed-loop poles that correspond to k,, = 40 have fallen into the overdamped region. The other two closed-loop poles are not shown in the figure because they are far away from the imaginary axis of the
Comprehensive dynamic modeling of flexible manipulators
145
20
16
m 12 ~.
k
=
10
I
8 4
)/ 8 -4
k, - 25
Ii ~
[-
k, = 20.69 ~
\\
-
=
-8 -12 k,=
-16
-2o.2o
i
l0
I -19
I
I 10
0
J 20 Real s
Fig. 7. Root locus for varying k,, with kp
=
]00 and k~ = 0.045.
40 35
--
k~= 18.11
k~- io
k,=40
30
20
..........
0
0.4
0.8
. .............
1.2
1.6
2
F
r
2.4
2.8
Time (s) Fig. 8. Closed-loop responses of joint angle with kp = 75, kyp = 0.045 and k,, = 10, 18.11, 40.
s - p l a n e . Since the s y s t e m t a k e s l o n g e r to d e c a y to its d e s i r e d v a l u e w h e n the p o l e s m o v e c l o s e r to the i m a g i n a r y axis, it can b e e x p e c t e d t h a t the r e s p o n s e s o f the s y s t e m will be q u i t e sluggish. This p r e d i c t i o n is v e r i f i e d by viewing t h e r e s p o n s e curves in Figs 8 a n d 9. M o r e o v e r , let kp = 100 while k e e p i n g k,, -- 10 a n d kfp = 0.045 as b e f o r e . It can b e e x p e c t e d t h a t k~ will n o t g e n e r a t e e n o u g h d a m p i n g to d a m p o u t the i n c r e a s e
146
Y.-J. LIN and T.-S. LEE
20]
,81
--k,-
10
k~- 18,11
- k~-40
16 14 '~12
~10
,/
~u 8 6
0
I
]
0.4
0.8
]
1.2
I
I
1.6
2
- - - T
2.4
]
2.8
T i m e (s)
Fig. 9. Closed-loop responses of applied force in the y-direction with kp = 75, k~p= 0.045 and k,, = 10, 18.11, 40.
oscillation generated by this kp value. Thus the percentage overshoot of the system has increased from 12.5% to about 25% for the position response while that of the force response has increased from 10% to 20% (see Figs 10 and 11). When the value of k,, is increased from 10, the two roots move towards the real axis along a circular arc as shown in Fig. 7. When k,, = 20.69, the roots are real and equal and the system is in the state of critically damped. In this case, the overshoots of the system should be reduced to zero. Figures 10 and 11 verify the claim. However, when k,, is increased beyond 20.69, the two roots become real and unequal, and the system is overdamped. Thus, the responses of the system are sluggish. This prediction is proved by inspecting the corresponding response curves in Figs 10 and 11.
CONCLUSION A theoretical study for the closed-loop position and force control of a flexible link robot with an end-point payload, maneuvered in an constraint environment, has been presented. The flexible arm's dynamic formulation is modeled as a Timoshenko beam. In addition, the gravitational effect is included in the modeling since the manipulator's workspace is three-dimensional. The comprehensive dynamic model along with the extended partitional control law architecture are then used to derive the control algorithm which deals with the manipulator's position and contact force simultaneously. The proposed controller gains tuning method is then applied to tune these controller's gains analytically. From the simulation results we observe that the performance of a flexible link robotic system can be predicted by investigating the roots of the system's characteristic equations. Hence, the total reliance on the
Comprehensive dynamic modeling of flexible manipulators
147
40 --k~'10
35
.....
k. = 20.69
...... k~-40
30 25 20 15 /
10 5 0
-5
I
I
I
I
I
I
I
0.4
0.8
1.2
1.6
2
2.4
2.8
Time
(s)
Fig. 10. Closed-loop responses of joint angle with kp = 100, kfp = 0.045 and k~ =
10, 20.69, 40.
20 18
--k~=lO
...... k o - 2 0 . 6 9
......
k~=40
16
"d
~
14
8 6 !
4
2 0 0
0.4
I
I
I
I
T
I
0.8
1,2
1.6
2
2.4
2.8
Time (s) Fig. 11. Closed-loop responses of applied force in the y-direction with kp = 100, krp = 0.045 and ko = 10, 20.69, 40. trial-and-error c o m p u t e r simulation to evaluate its p e r f o r m a n c e can be avoided or reduced. T h r o u g h the force and position controller design process it is seen that it is m u c h easier to design a multi-gains controller of a non-linear robotic system analytically using the p r o p o s e d root c o n t o u r m e t h o d than to design it by a trial-anderror a p p r o a c h for gains selection. In addition, since the position and force controller are affecting each o t h e r during m o t i o n , the root c o n t o u r m e t h o d gives better insight
148
Y.-J. LIN and T.-S. LEE
for controller gains selection. And from the motion control simulation results, the performance of the position and force controller designed analytically using the proposed method is quite satisfactory. This indicates that it is possible to extend the uses of the root contour method from designing controller gains of a linear system to tuning multi-controller-gains of a highly non-linear system such as a flexible link robot. Acknowledgements--The authors would like to acknowledge the partial support to this research from the Ohio Board of Regents through Research Challenge Grant # 5-34325 and the support obtained from the Office of Research Services and Sponsored Programs at The University of Akron.
REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.
Whitney D. E., J. Dynamic Systems, Measurement and Control 99 (1977). Hogan N., Proceedings of 1EEE Int. Conf. on Robotics and Automation 3 (1987). Hogan N., J. Dynamic Systems, Measurement and Control 107 (1985). Salisbury J. K., Proceedings of IEEE Conf. on Decision and Control (1980). Raibert M. and Craig J., J. Dynamic Systems, Measurement and Control 102 (1981). Shin K. G. and Lee C. P., Proceedings of 24th Conf. on Decision and Control (1985). Luh J. Y., Walker M. W. and Paul R. P., IEEE Trans on Automatic Control 25, 3 (1980). Khatib O. and Burdick J., IEEE Int. Conf. on Robotics and Automation 2 (1986). Khatib O., IEEE Trans on Robotics and Automation 3, 1 (1986). Lin Y. J. and Lee T. S., Proceedings of the 2nd National Conf. on Applied Mechanisms and Robotics, to appear (1992). Dorf R. C., Modern Control Systems, fifth Edition. Addison-Wesley, Reading, Massachusetts (1989). Craig J. J., Introduction to Robotics, second Edition. Addison-Wesley, Reading, Massachusetts (1989). Bayo E., lnt J. Robotics and Automation 4, 53-56 (1989). Timoshenko S. and Lessells J. M., Applied Elasticity, first Edition. Westinghouse Technical Night School Press, East Pittsburgh, Pensylvania (1925). Meirovitch L., Analytical Methods in Vibrations, first Edition. Macmillan, New York (1967).
0957-4158/92 $5.00+0.00 © 1992 Pergamon Press plc
COMPREHENSIVE DYNAMIC MODELING AND MOTION/FORCE CONTROL OF FLEXIBLE MANIPULATORS YUEH-JAw LIN and TIAN-SOoN LEE Department of Mechanical Enginccring, The University of Akron, Akron, OH 44325, U.S.A.
(Received 16 September 1991; accepted 30 October 1991) Abstract--This paper presents a comprehensive dynamic modeling of a flexible link manipulator by using Hamilton's principle. The model is formulated by considering the flexible arm as a Timoshenko beam model. Hence, the shear, bending as well as rotational inertia effect of the manipulator are all taken into account in the dynamic formulation. In order to obtain a comprehensive dynamic model, the gravitational effect, which is missing in most of the flexible manipulator dynamic models, is also included in the formulation. Then, an efficient motion and force controller design method utilizing root contour analysis is proposed. With the proposed method the multi controller gains of the position/force controller for the flexible manipulator can be chosen analytically, as opposed to the method of trial and error used conventionally for selecting manipulator controller gains. The tuning process of these controllers' gains proves the simplicity of gains selection. And the motion/force simulation results verify the effectiveness of the controllers. The simulation results also show that the derived control scheme, which is based on a linear control system analysis, can drive the highly non-linear flexible arm to achieve the desired position and force satisfactorily.
INTRODUCTION T h e main objective of this w o r k is to advance the state-of-the-art in the area of m o t i o n and force control of flexible a r m manipulators. T h e study is based on a flexible single-link m a n i p u l a t o r o p e r a t e d in a vertical plane as shown in Fig. 1. T h e driving t o r q u e is p r o v i d e d by an a c t u a t o r m o u n t e d on the m a n i p u l a t o r base joint. Position, velocity and force sensors are used to m e a s u r e joint angular position, angular velocity as well as wrist position, velocity and contact force. The study in this w o r k is aimed at developing an effective and systematic m e t h o d o l o g y to accurately control the m o t i o n o f flexible manipulators subjected to external contact forces. C u r r e n t research on flexible link manipulators has the following problems, namely: (1) most flexible m a n i p u l a t o r models u n d e r study are either considering the manipulator o p e r a t e d on a horizontal plane or the shear and gravitational effects are excluded in the d y n a m i c model, which m a k e s it difficult to control the m a n i p u l a t o r m o t i o n accurately; (2) the m o t i o n of the flexible manipulators being considered is usually a contact-free operation, i.e. without any contact forces, which m a k e s it impractical when p e r f o r m i n g industrial assembly tasks, and (3) the selection of controller gains is usually by empirical m e t h o d s or trial-and-error approaches, which m a k e s the controller design very time consuming. In view o f the a f o r e m e n t i o n e d p r o b l e m s the flexible link manipulator, u n d e r study in this paper, is constructed in such a way that it is o p e r a t e d vertically and in 129
130
Y.-J. LIN and T.-S. LEE (F.)
Y0
~+
T X0
Fig. 1. Configuration of a single link flexible arm manipulator.
addition, subjected to a contact force at the tip. The dynamic model is comprehensively formulated by considering shear, bending, rotational inertial and gravitational effects. The inclusion of gravitational effect is considered to be necessary since most industrial robots are operated in three-dimensional workspace. In the derivation of control algorithm of the flexible link manipulator, this paper also deals with external contact forces which cannot be avoided when performing assembly tasks. However, the research in this area is missing in the literature of flexible link manipulators as far as the authors are aware. Moreover, in this paper we propose a quick and systematic way of selecting good controller gains utilizing the root contour method. This method greatly reduces the laborious trial-and-error approach conventionally used in gains selection. Although the research in force control of flexible manipulators is rare, many research articles studying the force control problems of rigid link manipulators, which may be valuable to this study, can be found in the literature [1-9]. Whitney [1] in his early work proposed a feedback force control algorithm called accommodation control which modeled the end-point impedance as a dashpot. The controller is not considered accurate and is vulnerable to stability problem because it does not compensate for the manipulator dynamics. Hogan et al. [2, 3] developed a force control scheme which they called impedance control. This control method is evolved from observation of behavior of admittance and impedance subjected to an input. According to their definition, when a manipulator is interacting with an environment, the manipulator should assume the behavior of an impedance while the environment is an admittance. Another prominent force control method is the active stiffness control [4]. It started by defining a linear dynamic relationship between the contact
Comprehensive dynamic modeling of flexible manipulators
131
force and the deviation of the hand position from its desired position. This relationship in turn is transformed into joint space for implementation. In order to control position and force simultaneously, the hybrid position/force control technique [5] is developed. Since the controller of this method only consists of error driven servo, it is expected that the controller will not perform its tasks accurately. Recognizing the limitation of the hybrid position/force algorithm, Shin and Lee [6] have derived an extended resolved acceleration control scheme [7] which takes the manipulator dynamics into account. However, this method is devised for position control rather than for force control. Hence, to control position and force simultaneously, the extended resolved control method is proposed [8, 9]. The difference between the resolved acceleration and extended resolved acceleration control method is that in the latter, in addition to the position error, the force error term is included in the error driven servo. Since the controller implements its tasks in joint space, while the desired position and force are specified with respect to an external frame, considerable computational time will be consumed in calculating inverse Jacobian and homogeneous transformation matrices. As mentioned previously, the main objective of this paper is to develop an efficient controller which simultaneously deals with the flexible link manipulator's position and contact force and takes the manipulator dynamics into consideration to assure accuracy. A quick and systematic approach for tuning controller gains using the root contour method is proposed [10, 11]. Basically, the controller is developed based on the so-called partitional control law architecture [12]. However, it is extended in this paper to include both force and position errors driven servos. The artificial and natural constraint concept is used to specify and identify force and position control directions. This control scheme is different from others [1-9], which are used to control position and force of rigid link manipulators, in the sense that the Cartesian space position and force are indirectly controlled by joint space position and force. In addition, the dynamics of gravitational effect, flexible effect, and contact force are included in this controller design to reduce their influences on the accuracy of the motion of the flexible arm manipulator. In the following sections the detailed derivation of the control scheme based on the proposed methodology is described. Since the shear force depends on the shape of the cross-section of the arm, the numerical factors of three typical manipulator cross-sectional shapes, namely, hollow round, rectangular, and square, are derived in the next section. Following that, a complete set of dynamic equations of the flexible single-link manipulator with an endmass are formulated using the Hamilton's principle. The arm is considered to be a Timoshenko beam model when deriving the dynamic equations. The inclusion of shear force effect is required only in the cases of flexible robots with medium slenderness ratios subjected to very high speeds and with high precision requirements, as was addressed in [13]. For other cases, the shear force effect can usually be neglected for simplicity. The derived dynamic equations are solved by using the Galerkin's method. In Section 4, the control law that is used to control the motion of the flexible link manipulator subjected to external contact forces is developed. In addition, an efficient and systematic approach for tuning the controller gains is proposed. Finally, the motion simulation results of the flexible single link manipulator and the most
132
Y.-J. LIN and T.-S. LEE
important findings with the simulation are reported to verify the effectiveness of the proposed methodology.
N U M E R I C A L FACTORS As mentioned in the previous section, the upon the shape of the cross-section of the called numerical factors which represent the be needed in the dynamic formulation. Referring to [14], the numerical factor of a
shear force of the flexible arm depends arm. Therefore, the physical quantities characteristics of arm cross-sections, will flexible b e a m is defined as
Ib AQ
k' --
(1)
where I is defined as the second area m o m e n t of inertia of the cross-sectional shape of the b e a m computed with respect to its neutral axis. Q denotes the first m o m e n t about the neutral axis of the area contained between an edge of the cross-section of the b e a m parallel to the neutral axis and the surface at which the shear stress is to be computed. A and b represent the cross-sectional area and the width of the cross-sectional area at which the shear stress is required, respectively. Due to the fact that the numerical values of k' for most flexible beams are less than unity, there is a need to define the so-called correction factors k, which are inversely proportional to their corresponding numerical factors k' ]14]. Hence, k -
1
k'
-
AQ Ib
(2)
In this section we will derive the correction factors k for three typical robot arm cross-sectional shapes, namely, hollow round, hollow rectangular, and hollow square. Let us start with considering a b e a m with a hollow round cross-sectional shape. The inner and outer radius of this shape are given as r I and r 2, respectively. Thus, the approximate k for this shape is
4 (r~ + r2r 1 + r~) k
3
(r~ + r~)
(3)
for the thick type. However, for the thin type of this shape, the approximate k is obtained to be k = 2. Similarly, for a b e a m with a hollow rectangular cross-sectional shape, the height and width of the inner rectangle are given by h 1 and bl, respectively. The height of the outer rectangle is denoted by h2 and the width is represented by b2. The approximate ks for the thick and thin type of this shape are obtained as:
3 h l ( b 2 h ~ - b~h~) k - - ~ ( b 2 h ~ - blh 3) '
(4)
and k = 1, respectively. Notice that in obtaining these two ks the web area is used instead of the cross-sectional area. The hollow square cross-sectional b e a m means that the height and the width of the cross-section of the b e a m are equal. Thus, by setting
Comprehensive dynamic modeling of flexible manipulators
133
bl = h i and b 2 = h2, the correction factor k of the hollow square is therefore given by 3 hl(h~ + hlh2 + h~) k - - ~ (h2 + hl)(h~ + h12) "
(5)
If the web and the flange of the hollow square are thin, then it implies that ht ~- h 2 . Thus, k = 1.125. The values of the correction factors k of the hollow round, rectangular and square shapes are tabulated in Table 1.
COMPREHENSIVE
DYNAMIC
MODELING
Figure 1 shows a flexible arm with a concentrated mass at the tip of the arm. It also shows that the arm is subjected to an external force F e. On the one hand, the concentrated mass can be considered as the mass of a tool which is carried by the arm to p e r f o r m a constrained task. On the other hand, the external force can be considered as the contact force generated by interacting the arm with its environment. In the two-dimensional case, the contact force only consists of x and y components. The other end of the arm is clamped into the m o t o r shaft. The figure also shows that the length of the arm is L, the mass per unit length of the arm is p ( x ) , the mass m o m e n t of inertia of the actuator is Jh and the mass m o m e n t of inertia per unit length of the arm is Ja" The X o Y o Z o coordinates are chosen to denote the inertia frame of the system where the directions of X0-, Y0- and Z0-axis are the X0-axis points to the right, Y0-axis is normal to the X0-axis and Z0-axis is determined by completion of the right-hand rule. A rotational frame X ~ Y I Z 1 is attached and rotated simultaneously with the flexible arm. This rotational frame is used to describe the deflection of the arm. D ( x , t) denotes the transverse deflection of the arm from its equilibrium state at distance x along the Xl-axis and at time t. The assumptions of the dynamic derivation are as follows. (1) (2) (3) (4) (5)
Axial elongation, coriolis and torsional effects are assumed to be negligible. The frictional force and backlash of the m o t o r are ignored in modeling. The motion of the arm is in a vertical plane. The arm is considered to be a T i m o s h e n k o beam. The arm is considered to be uniform and prismatic. Table 1. Correction factor k for hollow round, rectangular and square Type
Thick
Thin
Hollow round
4 (r~ + r2rl + r~) 3 (r~ + r~)
2
Hollow rectangular
3 h~(b2h~- b,h~) 2 (b2h~- blh~)
1
Hollow square
3 h~(h~ + hlh2 + h~) (h2 + h,)(h~ + h~)
1.125
134
Y.-J. LIN and T.-S. LEE
Hamilton's method is used to derive dynamic equations of motion of the flexible manipulator. For a non-conservative system, Hamilton's principle is stated as
t2 f,
i
((~(T e - ge) -t-- a W e ) d t = 0,
(6)
where Te = kinetic energy of the system, V e = potential energy of the system, We = work done by the non-conservative forces to the system. Figure 2 illustrates the displacement of the arm from its equilibrium state at location x and the time t after the joint has rotated 0 degrees. In view of this figure, the kinematic relationships of any point on the flexible link described in terms of the inertial frame is given by, Y =xsin0+
Dcos0,
(7)
X = xcos0-
Dsin0.
(8)
The kinetic energy of the flexible arm system is the combination of the kinetic energy of the flexible arm, endmass and actuator. In view of these, the total kinetic energy of the manipulator T e expressed in terms of the inertial frame is therefore given by,
1 "2 + re= 1 f~ p(?2+22)dx+ m(X
?e) + 21 ft~
Ja~2dx + 1 fL
aO2dx +
] "2
(9) where the first term is the translational kinetic energy of the flexible arm. The second term is the translational kinetic energy of the endmass. The third term is the rotational kinetic energy due to the angular velocity caused by the bending moment. The fourth term is the rotational kinetic energy due to joint angular velocity. The last term is the kinetic energy of the joint itself. For simplicity, in Eqn (9) we neglect shear force effect and the Coriolis effect. The potential energy of the arm is caused by the internal bending moment and shear force as well as the gravity force. Therefore, the total potential energy of the system is
Ve - 2EI1
MZdx + ~
Vydx + mgY + pgLY,
(10)
where the first two terms are due to the bending moment and shear force, respectively. The last two terms are due to the gravitational effects on the tip mass and the arm itself. For the system shown in Fig. 1, the non-conservative forces are the input torque T and the external forces [F~]. In joint space, the work done by these non-conservative forces is
We = T O - [Jc]r[Fe]O,
(11)
where [Je] r denotes the transpose of the Jacobian matrix. For the two-dimensional
Comprehensive dynamic modeling of flexible manipulators
135
Yo
y, {F.}
X
L
xcosO
Fig. 2. Deflection of a point on a flexible arm described in terms of inertia frame.
case, the Jacobian is defined as
-Lsin(@+ O) [Jc]
(12)
=
However, the angle of deflection o~(x, t) at location x along the arm and time t is equal to the linear combination of angle of rotation due to the bending moment ~(x, t) and angle of distortion due to the shear y(x, t). It is measured from the equilibrium state and is given by aD cv(x, t) - ax - ~p(x, t) + y(x, t).
(13)
The bending moment M and shear force V are related to the ~p(x,t) and y(x, t) [14] by the following relationships:
a~
M ( x , t) = E1 5x
and
V(x, t) = k ' G A y ,
where E G A k' I
= = = =
Young's modulus of the material of the arm, shear modulus of the material of the arm, cross-sectional area of the arm, numerical factor, and second area moment of inertia of the arm.
(14)
136
Y.-J. LIN and T.-S. LEE
From the Hamilton's principle, we obtain dynamic equations of the arm with four boundary conditions as follows: -o
+
ato ! +
'CA
-
= 0,
(15)
and
-J~ 32~p + E132--~ + k'GA 3t 2
-
~X 2
x=0,
D =0,
x=0,
~p=0, a~p
x = L,
E I ~ x =0,
=
-~
+m
~p = 0,
(16)
~2D]
Ld20
+ ~-t2 ] = 0,
(17)
as well as dynamic equation of motion of the joint L [x2d20 32D]dx ( d20 ~2 D ] fo - P ~ dt 2 + x 3 t 2 ] m L2 - dt - 2 + L ~t 2 x=L]
- ( m g L c o s O - rngsinOD[~_L._
d20
d20
- Ja L - -
Jh--
dt 2
dt 2
"}-
pgL2cos02
-- pgLsinODIx=L/2
+
[Jc]T[Fe])
+ T = 0.
(18)
Notice that the dynamic equations of the arm and the joint, Equations (15), (16) and (18), are coupled equations. However, Eqns (15) and (16) can be combined into one equation by using the calculus concept that a tangent of a curve is equal to the first spatial derivative of the function of that curve. Thus, one obtains the relationship of 3DB ~P- 3 x ' (19) where DB indicates the deflection due to the bending moment, Substituting Eqn (19) into Eqns (15) and (16), the deflection of the arm can be written in terms of DB as
--
[ D2D13 1
D = Dn + k'GA
Ja
~t 2
32DB] E1
(20)
~X 2 ]'
and ~4DB - -
E1 3x 4
S2DB +
p - - -
3t 2
34 DB J.
-
-
~X29t 2
pEl ~4 D13 k'GA ~x2~t 2
+
PJa ~4DB k'GA ~ t 4
px
d20 dt 2
(21) Equation (21) is the dynamic equation of motion of the arm subjected to the forcing function - p x d 20/dt 2.
Comprehensive dynamic modeling of flexible manipulators
137
Similarly, the boundary conditions of Eqn (21) will be transformed into the following forms 3DB x=0, =0 ax
x = O, x-- L, x=
L,
1 ja 32DB - at - -2 E , DB + k'G-----A
~2DB] a-Tr-j = o
~2Dl~
E1 - -0 ax 2 53D B Ja ~t2a x
m[ a4O.
33D B E1 ax----T- + -7 J" %xZSt 2
a4 .l E1
9x 4 j = 0.
(22)
Because the closed form solution of Eqn (21) is difficult to obtain, the approximate solution will be sought. Galerkin's method [15], which is a two-step numerical method, is used to find the approximate solution. At first, it assumes that the solution of the dynamic equations of motion of this continuous system can be approximated by the linear summation of the solution of each discretized system. Each of these discretized system solutions consists of the spatial-dependent comparison function ~i(x) multiplied by the corresponding time-dependent generalized coordinate qj(t). n
Du = ~
(&(x)qj(t).
j=l
(23)
Then, it assumes that the difference between the approximate and the exact solution is small such that the accumulation of these weighted errors along the length of the beam is considered to be zero.
L
fo ~(X, t)q)r(x)dx = O,
r = 1, 2, 3 . . . . .
n.
(24)
in which r is a dummy variable and (Pr(X) are called weighted functions. Through a laborious calculation, the comparison functions are obtained as
/3,x
. I~jx
_~_
QliSm--L
~bj(x) = cos T -
- Q2jcosh , _ .
+ Qljsinh
/3jx
L '
(25)
where /3j = jzr j =
Q lj 1 Q2j =
1,2,3 ..... =
Clj C2j
(26)
,
(27)
k ' G A + k'G---~ 2 k'GA
Clj = E I
n,
(28)
k'GA
sinh/3j - sinflj) +
+ Ja)~}(Q2jsinhflj+sinflj)+
,
0
mE-~I(~)3(Q2ciosh[3-jcos3ip)j.
1
(Q2jcoshfij + cosflj)
(29)
138
Y.-J. LIN and T.-S. LEE C2j = E I
(cosh/3j + cos/3j) + - -
( sinh/3i + sin/3j)
p
+ JaZ~(cosh/3/
-
c o s / ~ , ) + r -e E- l P
"--'[~ZI3( sinh/3j \c!
-
(30)
sin/3j).
However, following the standard procedures of Galerkin's method, the forward dynamics of the generalized coordinates and joint coordinates are given by [ql = ([M] + P [ S I [ U I T ) - I [ ( [ K ] + p[S][H]r)[q]] + ~
P ([M] + p [ S ] [ U ] r ) - I ( T
- R 2 c o s 0 - [Jc]r[Fe])[Sl,
(31) and 0 = ( T - e 2 c o s 0 - [Jc]T[Fe]) + ( [ H ] r - [ U ] r [ M ] - a [ K ] ) [ q ]
(1 +
RI(1 + p [ U ] T [ M ] - I [ S ] )
p[CqT[Ml-l[S])
'
(32)
where the parameters in these two equations are as follows: K,j = [K] = E1
k' G A
~1L
d4~] dx,
~ q)r ~
dx - p
S,.= [S] = I,Lxq),.dx, L 3
R1 = ~ -
(33)
¢p,. dx 4
+ m L 2 + Ja L + Jh,
pgL 2 R~ = m g L + - - 2
I q)'qJj d x ,
(34)
(35)
(36) (37)
Notice that in the derivation of forward dynamics, the fourth order time derivative term is ignored. This is due to the fact that the contribution of this term to the final solution is small. The inverse dynamics can also be obtained from Eqns (31) and (32), which will be used for designing the controller.
CONTROLLER
DESIGN
The main problem facing us here is how to design a controller that can apply a desired force to the environment while simultaneously attaining a desired position at the influence of the flexibility of the manipulators. Since the generalized coordinates used in Eqn (31) are unmeasurable in a physical system, the control law will be derived solely from the inverse dynamics equation deduced from Eqn (32). In order to ease the derivation of the control law, the inverse dynamics equation based on Eqn (32) will be expressed in the compact form.
Comprehensive dynamic modeling of flexible manipulators
139
Introducing the notations W, Cg and f into the inverse dynamics equation from Eqn (31) yields W O + Cg + f + [Jc]T[Fe] = T, (38) where W = RI(1 + p [ U ] r [ M ] - I [ S ] ) , (39) Cg = R2 cos 0 f = -RI([H]
T -
[u]T[M]-~[K])[q].
(40)
(41)
In view of Eqn (38), the notations W can be interpreted as the inertia term, Cg as the gravitational influence term, and f as the flexible effect term. The control law conceived here has two parts. The first part of this control law is called a model base. It is used to linearize the non-linearity of the system and to reduce the system to an unit mass system. The other part of the control law is called the servo part. This portion of the control law is error driven. Since the model base has reduced the system to an unit mass system, the servo portion is designed to regulate this unit mass system. In view of Eqn (38), the control law has the following form: T = r/Uf + /3,
(42)
where r / a n d / 3 are the model bases. In view of Eqns (38) and (42), it is clearly shown that in order to make this system a unit mass system, r/and/3 should be chosen as r/=W, /3 = Cg + f + [Jc]T[Fe].
(43)
Thus, the dynamic equation of the unit mass system is given by Uf = 0.
(44)
The servo part is chosen such that the position and force errors of this unit mass system can be regulated. If the servo part is chosen as a linear combination of PD-controller operating on the position and velocity errors and P-controller operating on the force error, then we obtain Uf = kp(O d - O) + k,,(O d - O) + kfp[Jc]T([Fe] d - [Fe] ) + 0 a,
(45)
where kp and k,, are the proportional and differential gain constants for the position error and velocity error, respectively, kfp is the proportional gain constant for the force error. Combining Eqns (43) and (45) with Eqn (42), the closed-loop control law is therefore given by T = W[O d + kv(O d - O) + kp(O d - O) + kfp[Jc]T([Fe] d - [Fe])] + Cg + f + [Jc]T[Fe].
(46)
Equation (46) shows that the controller indirectly controls the Cartesian space position and force by joint space position and force. It also shows that the dynamics of gravitational effect, flexible effect and the contact force are compensated for to reduce their influences on the accuracy of the controller.
140
Y.-J. LIN and T.-S. LEE
Depending on the physical condition of an environment in which a robot is working, the contact force can be modeled as a damper, spring, inertia, or any combination of these three models. It is assumed that the robots are doing assembly tasks. Thus, the hand of the robot will come in contact with the workpieces. It is also assumed that the workpieces are composed of hard material. Hence, the contact force can be modeled as a high stiffness linear spring model. In joint space this contact force is expressed as (
[Fo] -- [k0ql[Jc] A0 +
AD'x=L)
L
(47)
'
where A0 = ( 0 - 00) and ADI,=L/L = (Dk=L -- Do]x=L)/L. In the two-dimensional case, the equivalent stiffness matrix
[keq ] is given by
in which the subscripts x and y refer to the x and y directions in the end-effector frame. An element in the equivalent stiffness matrix is given by 1
keqi =
1 +__+__
,
(49)
kei ksi kbi where ke~ and ks~ are referred to as environment stiffness and sensor stiffness in the direction i, respectively. Similarly, kbi is referred to as the flexible arm stiffness in direction i. An element of the equivalent stiffness matrix [keq ] is equal to zero when in that particular direction the motion of the manipulator is unconstrained. Hence, the control law becomes 7 = w[O ~ +
k,,(0 ~ - 0) + k~(0 ~ - 0)]
+ W[kyp[Jc]r([Fe] d -
[Keq][Jc]( A O + ADL[X=L))]
+ Cg + f + [Jc]r[Keq][Jc] ( AO + ADIx=L) L "
(50)
The block diagram of this closed-loop system is given in Fig. 3. Generally, the force and position controller gains can be obtained from the characteristic equation of the closed-loop system, or they can be obtained by tuning them experimentally. The characteristic equation of our system is obtained by combining Eqn (50) with Eqn (38) first. It is then linearized at certain [Fe] a and 0 d and expressed in s-domain as follows:
s 2 + k,,s + (kp + kfp[Jc]r[keq][Jc]) = 0.
(51)
The controller gains appear in the above equation require to be tuned properly to result in suitable control actions. It is proposed here that the controller gains be tuned by using the root contour method. It can be shown that the rules to construct the
Comprehensive dynamic modeling of flexible manipulators
141
model based
servO J
Delm.t "--/-- + eO
(
~
Equations (3~) and (32)
I
IF I I Fig. 3. Closed-loop control for a single link flexible arm manipulator using extended partitional control law.
locus of a single variable gain system are also applicable to the construction of the locus of a multiple-gains system. The first step of root contour method involves setting two of the gains in Eqn (51) equal to zero. Let us set k,, and kfp be zero. Then, Eqn (51) becomes s 2 + k s = 0.
(52)
The root locus of this equation can be obtained by dividing both sides of the equation by s 2. Therefore, the equation becomes
ks
1 +-7=0.
(53)
The root locus of Eqn (53) can be constructed through the use of the open-loop transfer function, ks/s 2. If the k s is selected a s k s 1 , then the next step is to restore one of these two gains, kfs, or k,,. In the case of restoring the gain k];., Eqn (51) becomes S 2 "~ kpl + kfp[JclT[keq][Jc] = O.
(54)
Dividing both sides of Eqn (54) by s 2 + kp~ results in 1 +
kfp[ Jc] T [ keql[ Jc] = o. S 2 + kpl
(55)
Again the root locus can be drawn from this equation. Similarly, if be kfpl, then Eqn (51) can be written as
1 +
kvs s~ + (/~si + krpl[&]T[/%q][4])
= 0.
kys is chosen to (56)
The root locus can be plotted using this equation. In addition, k~ is selected such that the desired behavior of the manipulator is attained.
142
Y.-J. LIN and T.-S. LEE C O N T R O L L E D M O T I O N / F O R C E SIMULATION
To test the effectiveness of the control law derived in the previous section, a motion simulation program has been developed to simulate the dynamic behavior of a single link flexible manipulator. In the simulation we consider the manipulator to be a Timoshenko beam model which is made of aluminum. The cross-sectional shape of the manipulator arm is assumed to be thin hollow round. Thus, from Table 1, the correction factor is 2. The gravitational constant used in this simulation is 9.81 m s -2. In addition, the mass moment of inertia of the actuator and the tip mass are given as 0.1 kgm 2 and 0.03 kg, respectively. The physical parameters as well as mechanical properties of the link used are tabulated in Table 2. Now, the first two modes of the approximate solution are used in the simulation to model the dynamic behavior of the link. This approximation is justified by the fact that the natural frequencies of the arm beyond second mode are laid outside the bandwidth of the actuator. The sampling time of the simulation is taken to be 0.01 s. The tasks to be performed by the controllers are maneuvering the manipulator's arm to the position of 20 ° and at the same time applying a step change force with magnitude 10 N to the workpiece. The desired direction of the step change force is in the y-direction of the hand space. In other words, the x-direction of the hand space is unconstrained. Since the workpiece is assumed to be hard, the stiffness of the workpiece and that of the sensor should be much higher than that of the flexible arm. Thus, the stiffness of the workpiece and the sensor are assumed to be [kw] =
[0°
0 ]
1.5 x 104
N m -1,
(57)
1.1 x 103
N m -1,
(58)
and [ks] =
respectively. However, the stiffness of the flexible arm can be calculated from
[0°
[kb] =
3E!
N m -t.
L 3~
Table 2. The physical and mechanical properties of the link Parameter
Value
Link length Outer radius Inner radius Flexural stiffness (El) Mass moment of inertia Density Modulus of elasticity (E) Modulus of rigidity (G)
2m 5.6 mm 4 mm 40.5651 Nm 2 0.1002 kg m2 0.2712 × 104 kgm2 7.1 × 10mNm -2 2.62 x 10mNm <
(59)
Comprehensive dynamic modeling of flexible manipulators
143
In view of these three stiffness matrices, the equivalent stiffness matrix is given by [keq] ~--
~
0-] Nm-1.
(60)
15
In order to achieve the specified goals, the controller gains should be chosen carefully. The root contour method is used to design the controller gains and through the assistance of the CC software the controller gains will be tuned accordingly. With this approach the locus of the controller gains will be traced out by varying one gain at a time. Then the gains are selected through the analysis of the loci. Figure 4 shows the root locus of varying the position proportional gain kp. Figure 5 shows the locus of varying kyp with a fixed kp = 75. In addition, Figs 6 and 7 indicate the loci of varying ko for two different kp values. As shown in Fig. 4, the locus of varying kp is on the imaginary axis of the s-plane. This means that the larger kp is chosen, the more oscillatory the response of the system will become. Therefore, kp should be assumed to be a reasonably small value to avoid oscillations. In our simulation, kp is chosen to be 75. As illustrated earlier, the root of the characteristic equation of kp is the open-loop pole of the locus of kfp. Thus, as seen in Fig. 5, kfp has the effect of adding oscillation to the system. In other words, kfp should be chosen as small as possible. In our case, we select kfp to be 0.045. In order to suppress the oscillatory effect of kp and kfp as well as the flexible mode of the arm while simultaneously to achieve the specified objectives, k,, is chosen to be 10 which is corresponding to the underdamped case (see Fig. 6). Thus, it can be expected that there will be an overshoot in the system response and that the response is very fast. Indeed, the real time simulations reveal that the overshoots of the position and applied force are about 12.5% and 10% from their desired values, respectively (see Figs 8 and 9). However, the rise time of this underdamped system is very short (~ 0.35 s). On the other hand, when the value of k v = 18.11 is chosen, Figs 8 and 9 show that 13 1197 531-
i
kp = 7 5
-1
-9--
] "~'~"-- kp - 7,.q
-11 -13
I ",60
,
I ",2g
I ,20
I
I ,68
x 1 Real s
Fig. 4. Root locus for varying
kp.
144
Y.-J. LIN and T.-S. LEE 16 12-
t
k #p -
0.045
I
4-
I,.. ~" k/p
-12 -
-lt2
I -11
i
~
0.0,t~
I
m 11
O
I
21 Real s
Fig. 5. R o o t locus for varying
16I 12
kt), with
k s, = 75.
kv~ 10
,e,
,-.
04 1
~=
//\ -4
-0!
-
= 25
~ i i -'----
k,= 1o
-12 -10
\\ __
k . = i8.11 " ~ k .
i
I
-11
I
I
I
i
10
i
20 Real s
Fig. 6. R o o t locus for varying k,, with
kp
= 75 and
kip =
0.045.
the overshoots of the system are reduced to zero and the position and force responses of the system are slowing down. In fact these two phenomena correspond to the system that behaves critically. These phenomena are due to the fact that by increasing the value of kv the damping effect of the system will increase. Referring back to Fig. 6, the locations of the closed-loop poles that correspond to k,, = 40 have fallen into the overdamped region. The other two closed-loop poles are not shown in the figure because they are far away from the imaginary axis of the
Comprehensive dynamic modeling of flexible manipulators
145
20
16
m 12 ~.
k
=
10
I
8 4
)/ 8 -4
k, - 25
Ii ~
[-
k, = 20.69 ~
\\
-
=
-8 -12 k,=
-16
-2o.2o
i
l0
I -19
I
I 10
0
J 20 Real s
Fig. 7. Root locus for varying k,, with kp
=
]00 and k~ = 0.045.
40 35
--
k~= 18.11
k~- io
k,=40
30
20
..........
0
0.4
0.8
. .............
1.2
1.6
2
F
r
2.4
2.8
Time (s) Fig. 8. Closed-loop responses of joint angle with kp = 75, kyp = 0.045 and k,, = 10, 18.11, 40.
s - p l a n e . Since the s y s t e m t a k e s l o n g e r to d e c a y to its d e s i r e d v a l u e w h e n the p o l e s m o v e c l o s e r to the i m a g i n a r y axis, it can b e e x p e c t e d t h a t the r e s p o n s e s o f the s y s t e m will be q u i t e sluggish. This p r e d i c t i o n is v e r i f i e d by viewing t h e r e s p o n s e curves in Figs 8 a n d 9. M o r e o v e r , let kp = 100 while k e e p i n g k,, -- 10 a n d kfp = 0.045 as b e f o r e . It can b e e x p e c t e d t h a t k~ will n o t g e n e r a t e e n o u g h d a m p i n g to d a m p o u t the i n c r e a s e
146
Y.-J. LIN and T.-S. LEE
20]
,81
--k,-
10
k~- 18,11
- k~-40
16 14 '~12
~10
,/
~u 8 6
0
I
]
0.4
0.8
]
1.2
I
I
1.6
2
- - - T
2.4
]
2.8
T i m e (s)
Fig. 9. Closed-loop responses of applied force in the y-direction with kp = 75, k~p= 0.045 and k,, = 10, 18.11, 40.
oscillation generated by this kp value. Thus the percentage overshoot of the system has increased from 12.5% to about 25% for the position response while that of the force response has increased from 10% to 20% (see Figs 10 and 11). When the value of k,, is increased from 10, the two roots move towards the real axis along a circular arc as shown in Fig. 7. When k,, = 20.69, the roots are real and equal and the system is in the state of critically damped. In this case, the overshoots of the system should be reduced to zero. Figures 10 and 11 verify the claim. However, when k,, is increased beyond 20.69, the two roots become real and unequal, and the system is overdamped. Thus, the responses of the system are sluggish. This prediction is proved by inspecting the corresponding response curves in Figs 10 and 11.
CONCLUSION A theoretical study for the closed-loop position and force control of a flexible link robot with an end-point payload, maneuvered in an constraint environment, has been presented. The flexible arm's dynamic formulation is modeled as a Timoshenko beam. In addition, the gravitational effect is included in the modeling since the manipulator's workspace is three-dimensional. The comprehensive dynamic model along with the extended partitional control law architecture are then used to derive the control algorithm which deals with the manipulator's position and contact force simultaneously. The proposed controller gains tuning method is then applied to tune these controller's gains analytically. From the simulation results we observe that the performance of a flexible link robotic system can be predicted by investigating the roots of the system's characteristic equations. Hence, the total reliance on the
Comprehensive dynamic modeling of flexible manipulators
147
40 --k~'10
35
.....
k. = 20.69
...... k~-40
30 25 20 15 /
10 5 0
-5
I
I
I
I
I
I
I
0.4
0.8
1.2
1.6
2
2.4
2.8
Time
(s)
Fig. 10. Closed-loop responses of joint angle with kp = 100, kfp = 0.045 and k~ =
10, 20.69, 40.
20 18
--k~=lO
...... k o - 2 0 . 6 9
......
k~=40
16
"d
~
14
8 6 !
4
2 0 0
0.4
I
I
I
I
T
I
0.8
1,2
1.6
2
2.4
2.8
Time (s) Fig. 11. Closed-loop responses of applied force in the y-direction with kp = 100, krp = 0.045 and ko = 10, 20.69, 40. trial-and-error c o m p u t e r simulation to evaluate its p e r f o r m a n c e can be avoided or reduced. T h r o u g h the force and position controller design process it is seen that it is m u c h easier to design a multi-gains controller of a non-linear robotic system analytically using the p r o p o s e d root c o n t o u r m e t h o d than to design it by a trial-anderror a p p r o a c h for gains selection. In addition, since the position and force controller are affecting each o t h e r during m o t i o n , the root c o n t o u r m e t h o d gives better insight
148
Y.-J. LIN and T.-S. LEE
for controller gains selection. And from the motion control simulation results, the performance of the position and force controller designed analytically using the proposed method is quite satisfactory. This indicates that it is possible to extend the uses of the root contour method from designing controller gains of a linear system to tuning multi-controller-gains of a highly non-linear system such as a flexible link robot. Acknowledgements--The authors would like to acknowledge the partial support to this research from the Ohio Board of Regents through Research Challenge Grant # 5-34325 and the support obtained from the Office of Research Services and Sponsored Programs at The University of Akron.
REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.
Whitney D. E., J. Dynamic Systems, Measurement and Control 99 (1977). Hogan N., Proceedings of 1EEE Int. Conf. on Robotics and Automation 3 (1987). Hogan N., J. Dynamic Systems, Measurement and Control 107 (1985). Salisbury J. K., Proceedings of IEEE Conf. on Decision and Control (1980). Raibert M. and Craig J., J. Dynamic Systems, Measurement and Control 102 (1981). Shin K. G. and Lee C. P., Proceedings of 24th Conf. on Decision and Control (1985). Luh J. Y., Walker M. W. and Paul R. P., IEEE Trans on Automatic Control 25, 3 (1980). Khatib O. and Burdick J., IEEE Int. Conf. on Robotics and Automation 2 (1986). Khatib O., IEEE Trans on Robotics and Automation 3, 1 (1986). Lin Y. J. and Lee T. S., Proceedings of the 2nd National Conf. on Applied Mechanisms and Robotics, to appear (1992). Dorf R. C., Modern Control Systems, fifth Edition. Addison-Wesley, Reading, Massachusetts (1989). Craig J. J., Introduction to Robotics, second Edition. Addison-Wesley, Reading, Massachusetts (1989). Bayo E., lnt J. Robotics and Automation 4, 53-56 (1989). Timoshenko S. and Lessells J. M., Applied Elasticity, first Edition. Westinghouse Technical Night School Press, East Pittsburgh, Pensylvania (1925). Meirovitch L., Analytical Methods in Vibrations, first Edition. Macmillan, New York (1967).