~
Pergamon
00945765(94)000816
Acta Astronautica Vol. 32, No. 9, pp, 561576, 1994 Copyright ,~ 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 00945765/94 $7.00 + 0.00
CLOSEDLOOP VIBRATION CONTROL OF FLEXIBLE SPACE SHUTTLE MANIPULATOR COLIN L. KIRK College of Aeronautics, Cranfield University, Cranfield, Bedford MK43 0AL, U.K. and FRANK L. DOENGI Institute of Statics and Dynamics (ISD), University of Stuttgart, Pfaffenwaldring 27, W7000 Stuttgart 80, Germany (Received 10 February 1993: revised 18 March 1993; received for publication 26 May 1994)
AbstractThipaper s develops a new approach to dynamic modelling and closedloop vibration damping
control of a Hermes type multilink flexible manipulator for twodimensional vibratory motion of the vehicle/manipulator system in the pitch plane. The flexible links are modelled by replacing them by rigid members and fictitious rotational springs at the servo motor joints to give a dynamically equivalent system. Linear equations of motion are derived by Lagrangian formulation for small motions relative to an arbitrary baseline configuration and eigenvalue analysis shows that mode shapes and frequencies are relatively insensitive to the angle between the two main manipulator arms. Closedloop modal control is then performed by velocity feedback from the active servo joints to the shoulder joint and a simple, stable damping control law is derived for the noncollocated sensoractuator system. The emphasis is on a simplified approach which facilitates physical interpretation of vehicle/manipulator dynamics and control, and which can be used to verify more comprehensive multibody computer codes.
i. INTRODUCTION In the design and operation of large space manipulators arm flexibility has a significant effect on system performance in terms of elastic vibrations induced by payload manoeuvres. In [1], Scott et al. indicate that longperiod oscillations persist in the U.S. Space Shuttle remote manipulator system (RMS) after routine operational manoeuvres. In the construction of the Space Station Freedom it has also been estimated that as much as 10h could be spent waiting for oscillations to decay to a 2.5 cm amplitude but this time would be reduced to 4 h by doubling system damping. It is therefore clear that closedloop damping and stiffness control of a manipulator system is essential to enhance its performance. A considerable volume of literature has been devoted to the study of stiff robotic arms which is primarily related to small and robustly constructed systems. A useful selection of material covering the fundamentals and applications in this field are given by Refs [2] [6]. It is difficult to specify the size of a manipulator at which the effect of flexibility becomes important but it is known to be significant for the manipulator systems of the U.S. Space Shuttle, the Space Station Freedom[l] and for similar systems such as the Hermes robotic arm [7].
The dynamics of a chain of flexible bodies was studied by Hughes in 1979 [8] with an application to payload manoeuvring in the Shuttle. The dynamics of a cantilever attached to a moving base and its relationship to certain multibody computer programs was considered by Kane et al. in 1987 [9]. Nonlinear equations of motion formed by recursive Lagrangian methods for flexible manipulator arms with rotary joints were formulated by Book in 1984110] for comparison of the computational time for rigid and flexible arms. Schmitz and Ramey [11] formulated the nonlinear equations of motion of a large flexible space manipulator using Kane's dynamical equations. A 2D experimental model was tested and compared with theoretical results for closedloop control, and implementation of various digital compensators for active control of a single elastic beam were considered. A study of the effect of joint flexibility and damping due to position and rate feedback for the Shuttle RMS was given by Barrows in 1984112]. Emphasis was on a large payload, linearization of nonlinear gearbox stiffness with brakes on and system response to RCS jet firing. A more comprehensive report on RMS dynamics is given in [I] where the nonlinear Draper RMS simulator code is used to predict arm dynamic response, including joint servo motors, gearboxes, arm flexibility and active augmentation. 561
562
(.'OLIN L. KIRK a n d FRANK L. DOENGI
The effect of flexibility on multibody dynamics has been considered by Vukasovic el a/. [13] using a fully Cartesian system of coordinates with an example of the Hermes spaceplane/manipulator/satellite payload system for large angle rotations, but closedloop control is not considered. The subject of elastic manipulator dynamics and control is closely related to the more general topic of the dynamics and control of large space structures (LSS) and in ,~articular to the Control of Flexible Structures (COFS) programme. The reader is referred to the extensive survey on structural dynamics and control research by Hallauer [14]. An example of the COFS programme is COFSII in which the dynamics and control of a large flexible antenna astromast system on the Space Shuttle is studied by Fay et al. [15]. Kirk and Moch [16] have also analysed this system for openloop pitch slewing and closedloop pole placement control using fullstate feedback of antenna and mast response into the antenna servo motor. The core of the present problem is that of an elastic fourarm manipulator driven by joint electric servo motors in openloop and closedloop operation. For large slew angles the system inertia matrix is timedependent and nonlinear, hence linear control theory is inapplicable. The aim of this paper is to study variation of system eigenvalues and mode shapes for small perturbations about an arbitrary fixed manipulator configuration. By this means the inertia matrix is linearized enabling modal analysis to be applied to transient and closedloop damped response. The simplification in the maniplalor arm elastic modelling assumes as in [7] that a single flexural shape function or degree of freedom for each arm adequately represents the dynamics, thereby eliminating the possibility of higher mode response involving two or more degrees of freedom per arm. The danger in this assumption is the possibility of spillover problems as pointed out by Hung and Weng[17]. Instead of following the usual approach of first producing a rational simplified structural model and then coupling it to the control system, they model the structure and controllers simultaneously using the Finite Element Transfer Function Matrix Method. One justification for neglecting higher order arm modes is that typical actuator torque pulse durations cxcecd 5s and have a Fourier translbrm given by r(m)\(I cos~,)v),,~,). Since the highest mode modelled has a natural frequency ~4 ~ 80 rad/s it is considered that higher mode response will be negligible unless much shorter pulses are used. 2. CLOSEI)I.OOP C O N T R O L OF FI.EXIBLE MANIPULATOR
A study of the literature reveals that a considerable amount of progress has been achieved in developing multipurpose computer codes for flexible multibody
systems which have been applied to the Shuttle RMS[I]. The present paper provides a nc~ and simplified approach to RMS closedloop dynamics for small amplitude payload oscillations in the shuttle pitch plane only, as I'ollo~s: (1) Arm flexibility is taken into account by replacing the elastic system by rigid arms and fictitious rotational springs at the joints which alh)w for servo stiffness and arm cantilever flexure. (2) The system inertia matrix is linearized b) evaluating it for small angular displacements relative to an arbitrary fixed arm configuration with closedloop stiffness in the joints. (3) The inertia matrix includes the payload, the shuttle and the manipulator moments of inertia. The distributed mass of the arms in the assumed rigid configuration is corrected to account for the actual elastic flexure. (4) The stiffness matrix yielded by the equivalent torsion springs at the joints takes into account the cantilever arm flexure or slope at the point of connection with the next arm. (5) An eigenvalue analysis is performed for the system for a range of included angles between the two long manipulator arms using joint stiffness corresponding to closedloop displacement feedback. This stiffness is also assumed to be equal to the case for brakes on.
(6) Modal analysis is applied to closedloop control in which ,joint relative angular velocities are fed back to the shoulder joint providing a generalized modal damping force which effectively couples the modes. (7) For quite large modal damping it is shown that modal responses are approximately uncoupled so that the feedback velocity gain constant can be directly selected to give specified modal damping. (8) The advantage of the modal approach is that resort to the much more complex pole placement method[16], applied to the original equations of motion, is avoided. A somewhat inconclusivc attempt at representing arm flexibility by means of fictitious elastic springs at the joints of the Hermes robotic ann {HERA) was presented by Woerkom [7] who refers to tile significant complex variations in system inertia matrix ik~r large angle motions and the development of an accurate or "truth" model of the system. The complexity of the model arises primaril,,, due to the 3D nature of the vehicle/arm/payhmd system, the elastic deformation and the various rigid bod.,, modes. 1! is hoped thai the new approach m the present paper using, albeit, a simpler 2D pitch planc analysis will help remove some o1 the conceptual difficulties in [7].
563
Flexible manipulator
........
Fig. 1. Schematic of Hermes Space Shuttle and robotic manipulator. 3. S I M P L I F I E D E I G E N V A L U E A N A L Y S I S
The aim of this section is to present a simple dynamic formulation of the Hermes robotic manipulator or similar system for operations in the pitch plane, without recourse to complex rigid or flexible multibody analysis. Emphasis is placed on designing a closedloop system in which rate sensors (tachometers) in the joints provide velocity feedback signals which are fed into the shoulder joint in order to damp out transient oscillations arising from manoeuvres or slews. The equations of motion are formulated for small perturbations relative to an arbitrary arm configuration that is assumed fixed. This approach yields a linear eigenvalue problem which permits the natural frequencies and mode shapes to be determined for a range of included angles between the main arms of the manipulator for various masses and moments of inertia. By this means a physical understanding of manipulator dynamics is obtained enabling a clear and simple closedloop system to be designed with specified modal damping. The decay time following transient disturbances can then be easily determined.
3. I. Formulation r?["the equations of motion A schematic of the Hermes robotic arm (HERA) is shown in Fig. 1 consisting of a 1.54m deployment arm, two 4.3 m main arms and the 1.3 m end effector. Each arm has a specified flexibility denoted by the El, where E is Young's modulus and I is the second moment of area of tubular arm in bending. Brushless d.c. serve motors exist at the joints and use reduction gear boxes to reduce the motor speed to an appropriate slow rotational rate, typically 1.66 rpm maximum. Although the gear teeth are known to have nonlinear stiffness (zero torsional stiffness for very small angular displacements and stiffening characteristics for increasing rotations) it will be assumed that the gear box stiffness for the brakeson situation is linear. In the case of brakes off the only joint stiffness present in closedloop operation is due to joint rotational displacement feedback into the serve motors. It can be seen that in order to have a positive definite system
stiffness matrix, joint stiffness must exist and is obtained by selecting feedback gains which give the same stiffness as the brakeson case. The effect of arm and joint flexibility will be modelled by assuming a rigid arm and an equivalent joint stiffness which gives the same tip deflection as the flexible arm. Thus, referring to Fig. 2 the deflection of the arm consists of a rigid arm rotation ,9 against the torsion spring of stiffness 2 and a cantilever beam flexure with tip displacement q. Under torque M the deflection q is given by q = ML2/2EI and the rotation by ,9 = M/2 where L is the link length. Introducing an equivalent rotation qJ = q/L an equivalent torsional spring stiffness 2~ is defined by 2c(0 + ~) = 2c~ = M which yields
M 2ZEI )~ (M +~IML)  2EI + L2"
(I)
Lagrange's equations are used to formulate the equations of motion of the coupled robotic arm/payload/vehicle system. We thus calculate the system kinetic and strain energy as follows. Inertia matrix. To derive the system inertia matrix using Lagrange's equations we refer to Fig. 3 showing the robotic arm in an arbitrary position defined by fixed angles 0o, 0~, 02 and 0~. The elastic degrees of freedom are defined by ~0, ~ , ~2, ~3 and the rigid shuttle degree of freedom by ~. The Hermes robotic arm (HERA) is stowed during launch at the back of the resources module (HRM) shown in Fig. 1, after which it can be relocated to various baseplates. Here, for simplicity, we assume M
I Fig. 2. Modelling of elastic manipulator arm by rigid link and fictitious torsion spring at serve motor joint.
564
( o l IN [., KIRK a n d FRANK L. DOEN(;I
Payload
," ~
~"
%\
• LI
I ~ O( 0 Jr" ~
~
~C
\ \o
3\'\ _ /
%,
I"v~ "'~
b"
'
. . . . . C
t
.
S
Shuttle Fig. 3. Coordinate system for dynamic analysis of shuttle/manipulator/payload system in the pitch plane. that the deployment arm is located above the Hermes mass centre. Vehicle rotation will not therefore give a velocity component of the deployment link in the Y direction. Also due to the large inertia of the vehicle and the small dynamic torque at the base rotary joint, point 0, the X component of velocity is ignored. This leads to considerable simplification in evaluating the kinetic energy and hence mass matrix of H E R A and the payload. Payload. We require the absolute translational and rotational velocity of the payload of mass Mp and moment of inertia Ip. The translational velocity of point P is found by summing the x and y coordinate velocity components of the end points of the previous links. Thus from Fig. 4 for the ith link, assuming
where for i = 3 L; = L* = L 3 + e, e being the eccentricity of the payload mass centre. Noting that for 0o = g/2, sin 00 = 1 and cos0o = 0 we have Tp = ½lp0i~ + ~
[(0~0L 0 + zi, L, sin 0, + ~2L~_sin 02
+ ~3L~ sin 03)' + (0~I L I cos 0~ + ~2L~ cos 0: + 0~3L* cos 03):]
(4)
The system matrix is then obtained from d/dt(tIT/&i;) which for the payload terms gives %:
MpLo[Lo~i o + L~ ~iI sin 01 + L2~i2 sin 02 + L*&"3 sin O~]
~:
MpLl[Lo6gosin 01 + Li6il + L2~': cos(01  02) + L~ ~cos(01  03)]
t,~.; =  L,~; sin 0, r~.i = L;~; cos 0,
(2)
~, :
M p L 2[Lo6io sin 0~ + Ll ~1 cos(01  02 ) + L26:2 + L* ~3 cos(0:  03 )]
Thus the payload translational velocity components %:
are
MvL't[Lo6iosinO~ + Li6:lcos(O~03) + L2 ~2 cos(0~  03 ) + L ~ ~3] + lp &'x
V,.p =  ~ L,~;sin 0, i 0
t'~.p = f
L~, cos 0,
(3)
I=(]
V. , Y
Vi
(5)
Link mass terms. Although the manipulator arms are light compared to the payloads normally considered, in the case of small payloads or no payload, arm inertia terms are required. For the deployment link L 0 the inertia term is simply Io~0 where Io = moL~/3. To illustrate the general procedure we consider links Lt and L2 as shown in Fig. 5(a) and (b). The mass per unit length of the link, assumed to be rigid, will be modified as shown subsequently, to allow for the actual cantilever flexure. For L~, see Fig. 5(a), the x and y velocity components of a point on the arm distance x from the lefthand joint are t', =  L0~o  .v:~lsin 01 r, = x ~ l cos 01
0i X Fig. 4. Velocity components for extremity of ith manipulator arm.
(6)
thus
T, = !m,
f
Li
W:, + ;,~)d.v )
(7)
Flexible manipulator
Inert& matr&. Finally assembling all the inertia terms the symmetric inertia matrix M = [M,;] = [M;;] is obtained with the following elements:
" LI &l
•
~o~
(a,
~
~
o ~t
(xv:,~
MII = I~ Mi2 = Mla = Mj4 = MI5 = 0
aolLo Ll&lcos Ol
565
M22 = I0 + Lg(ml Lt + m2L2 + m3L3 + M o) M23 = LoLI sin 01(½m I Lj + m2L 2 + m3L 3 + Mp)
~ L~2
(b)
M:4 = LoL 2 sin 02(½m2L2 + m3L 3 + Mp) M2s=LoL3sinO 3 m 3 L 3 + M p l + M33 = I 1 + L~(m2L2 + m3L 3 + M n)
Lo~ o + Li sin 0i fil Fig. 5. Velocity components for determination of kinetic energy of manipulator links L~ and L2.
M34 = L 1L2 cos(01  02)1½m2 L2 + m3 t3 +Mp
1+
which after inserting (6) and integrating yields
T, = ½m, (C~oL,no + Lot~ sin 0, a0~, + ~Lfaf)
(8)
where m~ is the uncorrected mass per unit length of the link. The terms for the inertia matrix are obtained by differentiation as %:
ml(L2Ll6~o + ½LoL~ sin 01~1)
~,:
m,(½LoL~sinO,~o+½L~6i~)
~ = ~ m2 I0 L2 [(aoLo+~L~sinOt+~2xsin02) 2 (10)
To avoid the algebra of squaring and integrating for Lagrange's equation we differentiate under the integral sign with respect to cq, i = 0, 1,2 and t and then integrate, obtaining
~o
m2[LoL,.6~o+ LoLtL2sinO, 6i, + ½LoL~sin02~:]
~i
m2[LoLiLzsinOl£o + L~L26il
+~t~2]
T = 5m I
f0 fo[
(x0 + 3))2 dx
x~+qtz)
But q = (L22/XEI)O =pO hence
m2[½LoL~ sin 02~0 + ½Lt L~ cos(0,  02)~'~
(11)
The terms corresponding to link L 3 are derived analogously. Shuttle inertia torque. The shuttle inertia torque l~a~ is added, where I~ is the moment of inertia about the vehicle mass centre G, Fig. 1, giving a system with the 5d.o.f. ct~,ct0, ~q, ~2, ~3.
(12)
Effective mass distribution offlexible arm. To allow for the effect of the actual cantilever flexure of the arms on the distributed mass we consider Fig. 2 and assume an elastic deflection shape function of the form y = q(x/L) 3 which allows both shear force and bending moment at the righthand end. Curve (1) denotes the actual elastic deflection and line (2) denotes the position of the assumed rigid arm. By comparing the kinetic energies of the two configurations we obtain an effective mass m e as follows. As discussed previously we note that q/9 = L22/2EI thus for configuration (1) the kinetic energy is
=½m
+ ½L, L~ cos(01  02)&'2]
~3
M,5 = 6 + (L~ + e):M~ + Ip I; = m;L3/3, i = 1,2, 3
(9)
For link L 2, Fig. 5(b) shows that the velocity of a point at distance x from joint 2 is found by considering the x and y velocity components of joint 2 due to link L 1 and the velocity due to rotation of L 2. Following (7) this gives
+ (~l LI cos 01 + ~2x cos 02)2] dx
M,~ = 6 + L~(m3L3 + Mp)
J dx
where
(13)
p = (L2212EI),
+'if) J dx =
i .q2rirS+~pL2+~p2L] ~m~ t~
(14)
For configuration (2) x2(,9+ ~)2 dx
T = 7m~ do
05)
566
('OHN L. KIRK and FRANK L. D()r~:sca
but O = #L). ~2El which on substitution in ( 15), after integration gives
T:,m~,,4"L'(I+' "~

•
3\
M xI
Ol = _ ,
hence
01 O1

V
/~'12 =  ")~e0ro
K22 = 2~or~ + ,:.~l(I + RII)2
_2EI~,~
L121 /
K23 = ,~lrl(1 + R0)
buttp,=q,/L,
K33 = )o~tr~ + 2~2(1 + Rt) 2
q~ 
2El
K34 = 2~2r2(1 + RI)
(19)
K44 = ).c2r~ + 2¢~(1 + R2) 2
Lt;.i
2EIqJ 1 2Elql LI2~  L~2 I
2EI L~2~+2EiOq
K45 =  t~.~r3 ( 1 + R2 ) (20)
or in general for the ith arm rotation against the rotation spring ,9 i =
rig * where
r; =
2El Li2 i + 2El
d I (x),]
, = , . ~  L,
L,).,+2E/al
(22)
Tangent at j o i n t 2
l
],
__
KI3 = KL4 = K]5 = K24 = K25 = K3~ = 0
(26)
3.2. Eigenanalysis of the equivalent rigid system (21)
Again from Fig. 6 the slope of L~ at joint 2 due to flexure only is q S , = d x q' L
(25)
K~] = ).~o
LI).I 2El
l+
,9~
r 3 9 3 ] 2]
and using (?V/g'.g, the terms of the symmetric stiffness matrix K = [K,I = [K,I are obtained as
L~ :q
Finally
(24)
5lz,0(zqrog0)2+2~l[%(l+Ro)rlgl] z
  Q g2] 2 + }.e319~2(1 I R 2 ) 
From Fig. 2. using (18).
hence
and i + I is
+ ,~.~2[gi(l + RI)
(18)
0q
(23)
The strain energy is then
Hence the distributed mass terms in (12) are modified using (17). St(if}tess matrix. To evaluate the strain energy in the joint torsion springs it is necessary to consider relative rotation between the adjacent arms. We therefore refer to Fig. 6 and note the following relationships:
(
3L,2
L, 2. + 2El
(17)
I+
qt=Ll(g]at)=Li
R.
7, = qS, + [~1, 0; ~ l = (1 + R i ) g ,  r, ig,~j
l + g Z +7\L/
and
where
The angle of rotation between arms therefore
Ol c = D~
ql MLI L i 2EI
0 , = R, gi
(16)
Equating (14) and (15) gives
01
In general
The equations of motion for the system are written in matrix form as M ~ + Kx = bz~
(27)
where x = [cq, %, gE, g2, g3] r, r~ is the control torque at joint 1, and b = [0,  1, 1,0, O]r. The eigenanalysis
to L I
~ / "
~
.""
L',o,
I~01
L2
I
.T
Fig. 6. Link geometry for determination of strain energy in fictitious torsion springs at manipulator drive joints.
567
Flexible manipulator for a fixed configuration of the system yields the modal matrix X of eigenvectors normalized such that XrMXii + XTKXq = Xrbt
mine ct and c2 but for three modes, ct, c2 and c3 will be required. The equations for the first two modes are then written #1 + og~ql = Rl[(Cl l/ill "~ C2 I/J21)ql
(28)
+ (g ~J2 + c2~22)q2]
gives the decoupled modal matrix equation + f 2 q = Rt.
q2 + o9~q2 =  R2[(cl ~H + c2 ~P21)ql
(29)
Here, q = [q~, q2, q3, q4, qrb]t is the vector of modal coordinates, f~ the diagonal matrix of squared eigenvalues and R = Xtb is the modal input vector, q~b in vector q is the uncontrolled zero frequency rigidbody mode of the coupled vehicle/robotic arm/payload system. Numerical results for the first four eigenvalues are given in Figs 9 and 10. 4. POLE PLACEMENT MODAL CONTROL
The determination of the eigenvalues and eigenvectors of the coupled Hermes/robotic arm system in the pitch plane for small amplitude oscillations relative to any fixed arm configuration provides useful data for the design and analysis of both open and closedloop control of flexible arms. It is also essential to minimise orbit operational time by effective damping of the payload oscillations caused by manoeuvres (see Scott et HI. [1]). This paper assesses the problem of closedloop damping of the HERA system by rotational velocity feedback from joints l, 2 and 3, assuming that the deployment joint 0 is inactive. The feedback rotational rates are obtained from sensors mounted in the geared servo motors and are fed back to the shoulder joint 1 to control vibration of the system following a disturbance. To illustrate the procedure we consider the summation of the modal rates from the three joints in the first two modes giving a feedback signal  G = c, [(Xll  X01)q I Jr (X12  X02)q2]
(32)
+ (C I I//12 "F C2 IlO22)q2]
(33)
where Ri are the corresponding elements of the modal input vector R. Equations (32) and (33) are now coupled by rate feedback terms and combine to form a fourthorder characteristic equation containing the unknown gain constants c~ and c2. It is the aim of the analysis to select these gain constants to give specified closedloop poles in the complex splane. With the abbreviations gt = c~ ~P. + c2 ~v2~ and g2 = c~ ~J12 F C2I//22 we obtain ijj + o9~ql + Rlgl ql "F Rig2t]2 = 0
(34a)
i12 + o9~q2 + R2gl (1! + R2g2q2 = 0
(34b)
from which we recover ct and c 2 as g2 tP2t  gl ~P22 cl = ~'12~2j  ~ . ~P22' gl ~12  g2 ~PII c2 = ~12 ~21  ~P. ~22"
(35)
Taking Laplace transforms of (34) and eliminating ql, q2 yields the characteristic equation s 4 + s3(Rlgl + REg2) + s2(o9~ + o9~) + s(R~g~ o92 + R2g2og~) + o92o9~ = 0.
(36)
The characteristic equation is now factorized into (s + a~)(s + a2)(s + a3)(s + a4) = O, a, = fl, ogl + iogl X / 1  fl~,
a2 = fl,ogl  i o g , x / 1  fl~,
a3 = fl2o92 + io92x/1  fl~,
a4 = fl2o92  io92\/I  f12 (37)
[ C2 [(X21  Xll )ql t ( X 2 2  ~12)q2] I c 3 [(.](31  X21 )ql F (X32  ¥22)q2]
(30)
The sensors will yield a signal equal to superposition of the modal velocities which can be decomposed into the respective qi by D / A conversion and use of the fast Fourier transform algorithm. Here, the X u are the elements of matrix X, and the c~ are gain constants. We introduce ~li] ~ X ( j   X i
i./
(31)
as the rotation of link i minus rotation of link i  1 in mode j for joint i. Note that we use the modal rotations of the equivalent rigid system which take into account the arm end slope. In (30) it is seen that we have deliberately included feedback from all three joints with unknown gain constants c~, c2 and c 3. However, as only two modes have been considered it is clear that only feedback from joints 1 and 2 is permissible in order to deterAA ~2 9B
where the 1%are the modal damping ratios. The above equation can be written as s4 + s3(al + a2 + a3 + a4)
+ s2[al a 2 + (Hi + a2)(a3 + a4) + a3a4] + s[a3a4(al + a2) + al a2(a3 + a4)] + aj a2aaa 4 = 0
or S 4 "Jr 2S3(fllCOl "F fl2o92)
+ s2(o92 + 4fl, f12o9, o92 + o9~) +2S((D~O92f12+O92Oglfll) + O)10322 2=0
(38)
Equating the coefficients in eqns (36) and (38) gives the system of equations RIgl + R2g2 = 2(fllogj + fl2o92)
(39a)
o9~ + o92 = o9~ + o9~ + 4fl, fl2o9,o92
(39b)
R, og~gl + R2og~g 2 = 2(o9~o92fl2+ og~og~fl~)
(39c)
568
COLIN L. KIRKand FRANKL. DOENG1
It is obvious that there can be no exact solution to eqn (39b) except fl~ or f12 equal to zero. This is equivalent to the situation of uncoupled modes. However, for small fl, we can neglect the nonlinear term in (39b) and obtain the approximate solution for the unknown gains gj and g2 from eqns (39a) and (39c), thus:
2fll~l gl = T i
2fl2co2 ,
g2 
R2
(40)
The error involved in neglecting the term 4flj fl2co1092 in eqn (39b) is small even for say fl~ = f12 = 0.5. For example with col = 0.6 rad/s and co2 = 8 rad/s we have 4ill f12col co2/(co~ + (o~) = 0.075. Hence eqn (40) can be used for larger values of ft. It is necessary, however, to ensure that the gains g] and g2 do not lead to excessive control torque demand on the servo motor. Substitution of eqn (40) in (35) yields the gain constants c~ and c2. It is seen from eqn (40) that for given eigenfrequencies cot, co2 it is only necessary to specify the appropriate damping ratios fl~, f12 for the closedloop system. Thus far we have a simple result for determination of the velocity feedback gains to give pole placements for the first two modes for small modal damping ratios. For design purposes this may be adequate if only the first two modes are assumed to give the most significant response during manoeuvring the robotic arm. Furthermore, it can be shown for the case with three modes that the simple result [eqn (40)] can be generalized to
g,
2fl, co,
i = 1.2, 3.
(41)
Ri
Since our approach is equal to neglecting the coupling via damping terms in eqn (34) for small damping ratios fi. eqn (41) can be assumed valid for an arbitrary number of modes. 5. DYNAMIC RESPONSE
For a fixed arm configuration with closedloop stiffness in the joints the eigenvalues are constant and closedloop feedback damping can be used in a linear control system. Similarly for a slewing system the eigenvalues will change slowly with time for small accelerations and decelerations and the modal method can also be used but with the joints having zero stiffness. Variations of natural frequencies with the elbow angle 0~02 are shown in Figs 9 and 10 for joints having closed loop stiffness ). =4.5 × 105 Nm/rad. It is thus possible to apply the modal method to determine the openloop response during a slew and then use the closedloop damping system to eliminate the residual vibrations at the end of the slew. To simplify the procedure it is assumed that only the first two modes are excited by the application of a torque applied by the shoulder servo motor, thus z~0= 450 Nm as shown in Fig. 8.
5.1. Transient response The modal equations for transient response caused by a disturbing torque r~(t) at joint 1, neglecting joint friction and structural damping are ~/,+co2q~=R~rc(t )
i=1,2
(42)
where L(t) = _+r~oand R i is the appropriate element of the modal input vector R. The displacement time history of the payload can then be obtained by modal superposition. This procedure yields the modal initial conditions q~(t), (h(t) following the disturbance, to be used subsequently for the study of closed loop damping. Control criteria. In controlling the payload dynamic response the following factors are important: (I) The translational and rotational displacement of the payload relative to the shuttle. (2) The absolute translational and rotational motion of the payload in relation for example to an adjacent satellite or the Space Station. (3) Effective damping of the payload.
Static relative payload displacement. Referring to Fig. 7 the static displacement coordinates of payload P relative to point 0 are xs = L~ cos 0is + L2 cos 02~+ L* cos 03~ Ys = L0 + LI sin 01s + L2 sin 02s + L* sin 03s (43)
Dynamic payload displacement. The modal rotation vectors of the robot arms from the eigenvalue analysis are defined as
xi = [x0. x~,, x2~, x3,] j = l. 2
(44)
where X0 is the rotation of link i in modej. Assuming that the elastic modal response rotations of the links are small, the dynamic displacement components of P are
x ( t ) =  {L0lX01ql (t) + Xo2q2(t)] + L 1[Xilq] (t) + Xi2q2 (t)]sin 01 + L2[X21 ql (t) + X22q2(t)]sin 02
+ L3 [X31q~(t) + Xs2q2(t)]sin 03 }
y ( t ) = L 1[Xuq I (t) + X12qz(t)]cos Oi + L2[X21 ql (t) + X22q2(t)]cos 02 + L3 [X31ql (t) + X32q2(t)]cos 03
(45)
or
x ( t ) =  Fl qj (t)  F2q2(t) y(t ) = Gi ql (t ) + Gzq2(t )
(46)
Flexible manipulator
569
L~ L2 ~
I
i LI
T= ///

i
Xs
~[
//
Fig. 7. Geometry for static displacement components of payload relative to shuttle.
R i"cco qi(t) = ~ [2 cos ogi(t  To)
where
F l = LoXol + LiXll sin 01 + LzX21 sin 02
 cos wit  1]
+ L3X31 sin 0 3
(49)
Rj%0(1 _ 2 e  r ~ + e 2Td,) phase 3:
qi(s) =
s(s2 +~o2 )
F: = LoX02 F Li Xl2 sin 01 + L:X:2 sin 02
qi(t) = Riz~,° [2 cos ogi(t  Td)
} L 3X32 sin 03
 cos ~oi(t  2T0 )  cos o~t]
GI = LiXll COS 01 + L2X2j cos 02 + L3X31 cos 03 G2 = LI XI2 cos 0j + LI X2., cos 02 +
(50) Ri r~o
L3X32COS 03
qi (2 Td ) = ~
(2 COS ~i Td  COS 2oJi Td  1)
(47) Ri Tco q/(2Td) =    (2 sin o9i To  sin 2o~, T~)
(Oi
Response functions. F r o m Fig. 8 we obtain the Laplace transforms and the respective inverses of eqn (42) for the acceleration and deceleration phases I and 2 under step torques %0 and the residual free vibration phase 3 as follows: Ri zco
phase 1:
qi(S)S(S2q_O)2) Ri zeo qi(t) = o9~'(1  cos og,Q R~z¢0(l  2e r~.~.)
phase 2:
qi(s)
s(s2 + ~o2g)
l
Phase1
"[co ~
Td
(48)
(50 The previous equation gives the initial conditions for the closedloop control analysis at the end of the disturbance torques Zoo, i.e. for t = 2T0: Maximum response. The maximum response may occur during the acceleration or the deceleration pulse phases. For example the response excited during acceleration may be enhanced by the addition of the response due to the deceleration phase. On the other hand the responses from the two phases may partially cancel each other depending on the time at which the m o t o r torque is reversed. Acceleration phase. Substituting eqn (48) we find from eqn (46) that
Phase 3
Td 
'[co
Phase 2
Fig. 8. Step function disturbance torques r¢o applied to shoulder joint for T0s. Phase l: accelerating torque; Phase 2: decelerating torque; Phase 3:%0 = 0 with closedloop feedback velocity control of residual vibration.
570 x(t) "Qtl 1'(t) Zc0
COLIN U KIRK and FRANK L. DOENGI x(t)
F IR I /'; R~ ~ (1cos091t)z~(Icos092t) 09 ~ (O

& R~ v [2 cos .,, (t  T d )  cos '.')1 t .... I i

7"c0
U) 1
G,R~
GIRl
~ (Icos091t)+ 09 ~
F:v212cos092( . T,I) . . c. o s 0. 9 , t
" "(1cos092t) ~O"~ (52)
y(t)
_
Gt R
_1 [2 cos 091 (t  To)  cos 091t  I] 09T
Zoo T h e m a x i m u m values o f x ( t ) a n d y(t) are obtained f r o m dx/dt = 0 a n d dy/dt = 0 . Since the m o t i o n is u n d a m p e d a n d periodic the m a x i m a will repeat periodically d u r i n g the acceleration slew. D i f f e r e n t i a t i n g the a b o v e e q u a t i o n s gives the times at w h i c h Xmax a n d Ym,x o c c u r w h i c h are t h e n b a c k s u b stituted into e q n (52). T h e c o n d i t i o n s for m a x i m a read
1 [DI
Fi RI sin 09, t +
1 

F 2 R 2 sin
+
Differentiating the above equations with respect to t gives the following c o n d i t i o n s for the times at which the m a x i m a occur: 1 
091
O,I2
1
I
091
092
"
G2 R: 092 [ 2 c ° s 0 9 2 ( t  T o )  c ° s 0 9 "  1]
(54)
092 t = 0
GjR~sin09~t +   G . R . sin092t=O
FI Rl [2 sin 091 (t  To )  sin 091 t]
1
+   F 2 R 2 1 2 s i n 0 9 2 ( t  To)sin092t]=O
(53)
092
~
1 
T h e s e e q u a t i o n s are solved n u m e r i c a l l y for tmax. It is also n o t e d t h a t the m a x i m u m m a y o c c u r at the e n d o f the first phase, i.e. at t = T d w h i c h can be inserted into e q n (52) as a check• Deceleration phase. T h e p a y l o a d r e s p o n s e d u r i n g d e c e l e r a t i o n is f o u n d by s u b s t i t u t i n g e q n (49) into e q n (46) to give after t r a n s f o r m a t i o n
091
GIRi [2 sin 09~ (t  Td)  sin col t]
1
+   G2 R212 sin 092(t  Td)  sin 09:t] = 0 092
..)
~4
;
12.5
1.00 10.0 0.75 v
7.5
~.~
7.
0.50 "''
f2
''"
5.O 0.25 ~ \
/ /, N
fj
I I.5
1.0
(55)
A s b e f o r e the values /max are inserted in e q n (54)•
Mp = 3500 kg
1.25
0
I1
09
/
I
I
0.5
0
I 0.5
1.0
f3
2.5 1.5
Angle between main arms, 0 (rad) Fig. 9. Variation of the first four natural frequencies of the manipulator with the included angle 0 between the two main arms, for M p = 3500kg. Joint stiffness ). =4.5 × 105 Nm/rad.
571
Flexible manipulator
The constants A~, A 2, B~ and B 2 depend on the initial conditions and are given by
5.2. Closedloop damping The damped free vibrations of the two coupled modes are obtained from eqn (34) subject to the initial conditions of eqn (51). Taking Laplace transforms of eqn (34) gives
Ai = ql(0),
BI = ql(O) + Rlg2q2(O) + Rlglql(O),
[s2q, (s)  sq, (0)  q, (0)] + og~q I (s) + R~g, [sq~(s)
B2 = q2(0) + R2glql(O) + R2g2q2(O).
 ql (0)] + R~g2[sq2(s)  q2(0)] = 0
[s2q2(s)  sq2 (0)

#2(0)]
"l
(60)
Inversion. The Heaviside inversion formula gives the modal time histories as
o)2q2(s) d R2g2[sq2(s)
 q2(0)] + R2gl[sql(s)  ql(0)] = 0
A2 = q2(0),
(56)
qi(t) = ~ N,(ar)
r=l
Solving for ql(s) and q2(s) yields
e.,,
i=
q, (s) 
(Ajs + BI)(s2 + co~ + R 2 g 2 s )  R I g z A 2 s 2  RIg2B2 s Nl(s)  O(s) O(s)
q2(s) 
(A2s + B2)(s2 + o92 + R I g l s )  R 2 g l A 2 s 2  R2glBIs N2(s) _ __ D(s) D(s)
1,2
(61)
(57)
where
D(s) = s 4 + s3(Rlgl + R2g2) + s2(o9~ + o~)
where the N~ are obtained from eqn (57) and the poles a r from eqn (59). The derivative of D(s) is
+ s(Rlgloo~ + R2g2og~ ) + tolo2 2 2 (58)
D'(s) = 4s 3 + 3s2(Rlgl + R2g2) + 2s(to~ + ~ )
which has been previously factorised approximately, based on the assumption of small fl~ and f12, as
d(g,glo922 q R2g2o~)
D(s) = (s + a I )(s + a2)(s + a~)(s + a4)
which is evaluated for s = at. Evaluation of eqn (61) yields the damped closedloop responses ql (t), q2 (t). Alternatively, the q~(t) can be obtained from numerical integration of the equations of motion• The damped payload response is then obtained by insertion of the modal time histories into eqn (46).
where al. 2 = fll(Di
___
io9,~/1
 flf
and a3., = fl2co2 _ i~o2x / l  fl]
(59)
Mp
0.225
=
20,000 kg
m "''...°
\
m
&
. . . .
0.220
12.5
. j .°•,
10.0
"•''"••'""'''. • .. f2. .
0.175
.....
°.
..° ''' 7.5
v
0.150
5.0 i
0.125

I
2.5 f3 //
O. I00
\
I I .5
f
F .... I .0
t 0.5
(62)
"
q 0
I ..... 0.5
I 1.0
I 1.5
Angle between main arms, 0 (rad) Fig. I0. Variation of the first four natural frequencies of the manipulator with the included angle 0 between the two main arms, for Mp = 20,000 kg.
572
C o l IN L. KIRK a n d FRANK L. DOEN(;t Table I. Payload displacement maxima during transient rcspt)nsc Payload
0~
O.
0~
Small Sinai] Small Small Small Small Large Large Large Large Large Large
9(1 90 90 90 9(I 90 90 90 90 90 90 90
90 90 45 45 (I 0 90 90 45 45 0 0
90 90 45 45 0 0 90 90 45 45 0 0
Coordinate
15.01 15.01 19.51 15.113 14.95 14.95 14.95 17.58 14.97 12.57 15.41 18.45
6. R E S U L T S
Extensive numerical investigations have been performed using the analytical results from the previous sections for three configurations with two different payloads. In this section some o f the numerical results are presented and discussed. Referring to Fig. 3, the following data have been used for the geometrical and physical configuration of the system: Link lengths: Link masses: F o r all joints a n d links: Spacecraft inertia: (a) Small payload: (b) Large payload:
/,~,,, [s]
Maxim um (m) 0 lob5 • 08485 ~ 0 1588 (11436. 0 182X • 01042. 01640, 0.o891 • tl.2873. n 3362 ,
I0 I0 I11 I0 IO 10 l0 l(I I0 lu
!L4292
I0
.
11.4805 ,, 10
simultaneously varied from 90 to + 9 0 . At a n u m b e r of intermediate positions 0 = 0~  0: eigenanalyses have been made resuhing in Figs 9 and 10 for the small a n d large payloads, respectively. It becomes a p p a r e n t that lk~r large angular ranges, especially a b o u t 0 t  (L = t) . the eigenfrequencies of the system remain fairly constant, a fact that was used above for simplification of thc dynamic response studies.
L0 = 1.54 m, L~ = 4.3 m, L, = 4.3 m , L 4 = 1.3 m m0 = ml = m~ = m3 = 30 kg/m 2 = 4.5 x 105 N m / r a d , E1 = 3.0 x 105 Nm" /, = 300,000 kg m' M v = 3500 kg lp = 3430 kg m" e = 1.5 m M 0=20000kg Ip=600000kgm 2 e =3.4m
F o r b o t h cases (a) a n d (b) with the small a n d large payloads eigenanalyses of the linearized systems have been performed. The angle 00 was m a i n t a i n e d at 9 0 , whereas b o t h angles 0~ and 02 have been
Transient response calculations are based on application of positive and negative torque pulses L,, (Fig. 8) of m a g n i t u d e 450 N m and d u r a t i o n T d = 10 s. The m a x i m a of the payload displacement co
x 10 .2
Modal coordinate ql
5.0
M p = 3500 kg
2.5
=2 0 ~g
2.5
5.0
I 0
25
5O
I 75
I I00
I 125
k 150
I 175
i 200
t (s) Fig. I1. Closedloop response for first mode generalized coordinate q~ for fl~ = 0.03, [/_,= 0.01 and Mv = 3500 kg.
573
Flexible manipulator Modal coordinate q2
0.02
Mp = 3500 kg
0.01
0.01
0.02 0
50
I
[
I
I
I
100
150
200
250
300
t (s) Fig. 12. Closedloop response for second mode generalized coordinate ql for fll = 0.03, f12= 0.01 and Mp=3500 kg. ordinates have been determined using eqns (52)(55). It was found that all maxima occur during the deceleration phase as in Table 1. Using the initial values given by eqn (51) numerical integration of the closedloop time histories has been performed. For the two cases (small/large payload) with a straight manipulator (0~ = 90 °, 02 = 03 = 0) the results are given in Figs 1 !16. For the selected cases the payload displacement y coordinate is virtually zero. We therefore confine the discussion to the modal coordinates q~, q2 and the payload displace
ment x coordinate. Since the damping increases with the eigenfrequency of the mode, the modal damping ratio has been chosen highest for the first mode. For the computed cases we have set the modal damping ratios to fl~ = 0.03 and f12 = 0.01. The time histories of the modal coordinates q~ and q2 for the case with small payload are shown in Figs l l and 12, respectively. The appropriate eigenfrequencies are to~ = 0.6285 rad/s and ~o2 = 8.053 rad/s. In Fig. 11, the coupling of the modes via the damping terms in eqn (34) becomes
x 10 4 3
Payload
displacement:
x
2
1
~
0
I
2
3 II 0
I
I
I
50
I O0
150
L 200
I
I
250
300
t (s) Fig. 13. Closedloop payload response in the x direction for Mp = 3500 kg.
574
(~tJ:,; L, KIRK and FRANK L. DOENGI x It)
Modal coordinate ql Mp = 20,000 kg
5.0
25
2
0
2.5
5.0 I 50
0
L Ioo
t 150
L 200
I 250
I 300
I 350
I 40O
t ts) Fig. 14. Closedloop response for first mode generalized coordinate q~ for fl~ =0.03, //2 = 0.01 and Mp = 20,000 kg.
a p p a r e n t in the hightYequency ripples in the time history q~ (t) for small t. Since ~o2 is m u c h higher t h a n ogj the second m o d e is d a m p e d out m u c h sooner t h a n the first. A n interesting coincidence occurs when c o m p u t i n g the payload displacement in this special case: the oscillation period of the first mode is virtually equal to the acceleration a n d deceleration pulse lengths r j = 10 s. Therefore, the c o n t r i b u t i o n of this mode to the subsequent closedloop m o t i o n vanishes, i.e. the a p p r o p r i a t e initial conditions are
equal to zero. F o r this reason, the x coordinate time history s h o w n in Fig. 13 is completely determined by the second mode. F o r the case with large payload, presented in Figs 1416, the situation is different. Both modes, with ~oj = 0.6176 rad/s and to 2 = 1.246 rad/s, are d a m p e d out in similar times. This is due to co2 being m u c h lower t h a n in the smallpayload case. The coupling between the modes via the d a m p i n g terms is best seen in the time history of qz(t) in Fig. 15 when looking
x 10
2
Modal coordinate q2 Mp = 20,000 kg
i
1
l
0
50
i 00
150
200
250
300
350
400
t (s) Fig. 15. Closedloop response for second mode generalized coordinate qj for fit = 0.03, fl~ = 0.01 and Mp = 20,000 kg.
Flexible manipulator xl0 5 2
575
Payload displacement: x Mp = 20,000 kg
1
~
0
I 2
0
50
I 100
I 150
I 200
I 250
I 300
L 350
I 400
t (s) Fig. 16. Closedloop payload response in the x direction for Me = 20,000 kg. at the negative peaks for small t. The period of the first mode being sufficiently different from 10 s, the payload displacement x coordinate shown in Fig. 16 clearly exhibits the contributions of both modes. The maximum feedback damping control torque has been checked by substituting eqn (40) into eqn (34a) giving re(t) = 2 I ~ t
1q. + ~ 0 2 ] R ,
(63)
The feedback control torque required to provide stiffness at the joints is not considered specifically in this paper. Substitution of typical peak values of q~ and q2 shows that re(t) is small compared to 450 Nm. The gain constants have been obtained from (35) for the small payload case as c~= 1.14x 105(l/s), c2 =  1.17 x 104 ( 1 / s ) , and for the large payload case as c~ = 1.20 x 10S(l/s), c2 =  6 . 4 2 x 104(l/s). 7. CONCLUSIONS A new method of formulating linear equations of motion for a flexible multilink space manipulator/ shuttle system undergoing pitch plain disturbances has been presented, in which the elastic links are replaced by rigid members and equivalent fictitious torsion springs at the joints. For an arbitrary fixed baseline configuration of the links an eigenvalue analysis is performed enabling openloop equations of motion to be obtained in uncoupled modal form. Linearization of the equations of motion is justified on the basis that the natural frequencies and mode shapes remain sensibly constant for relatively large included angles between the long arms of the manipulator. AA 32 9 (
Closedloop modal damping control of the system is obtained by velocity feedback from the various joints into the shoulder joint and it is shown that velocity coupling effects between various modes are negligible for even high damping ratios. The analysis presented gives a clear physical insight into twodimensional manipulator motion and can be used to verify the results obtained from other methods (e.g. D I S C O S [18]). It is suggested that the present method could be developed for more compled threedimensional analysis of shuttle/manipulator systems.
REFERENCES
1. M. A. Scott, M. G. Gilbert and M. E. Demeo, Active vibration damping of the space shuttle remote manipulator system• AIAA912621CP, Proceedings of the AIAA Guidance, Navigation and Control Conference, New Orleans, 12 14 Aug., Vol. 1, pp. 194204 (1991). 2. J. Denavit and R. Hartenburg, Kinematic notation for lowerpair mechanisms based on matrices• A S M E J. Appl. Mech. 77, 215221 (1985). 3. D. Whitney, Resolved motion rate control of manipulators. IEEE Trans. Man Mach. Syst. 10, 47 53 (1969). 4. J. Luh, M. Walker and R. Paul, Online computational scheme for mechanical manipulators. IEEE Trans. Syst. Man Cybernet. 14, 444~450 (1980). 5. J. Hollerbach, Recursive Lagrangian formulation of manipulator dynamics and a comparative study of dynamic formulation complexity. IEEE Trans. Syst. Man Cybernet. 10, 730 736 (1980). 6. K. S. Fu, R. C. Gonzales and C. S. G. Lee, Robotics Control, Sensing, Vision and Intelligence. McGraw Hill, New York (1987). 7. P. Th. L. M. Woerkom, On fictitious joints modelling of manipulator link flexibility for the HERA simulation facility plot. National Aerospace Laboratory NLR, TR 88086 U, The Netherlands (1988).
576
('~t I': L. KIRK and FRANK L. DOENt;I
:~. P. C. Hughes, D~namics of a chain of flexible bodies. ,l, .4~'lrommt..S'ci. 27, 359 380 (1979), 9. T. R. Kane, R. R, Ryan and A. K. Banerjee, Dynamics of a cantilever beam attached to a moving base. AIAA J. Guidance Control Dvnam. 10, 139 287 (1987). 10. W. J. Book, Recursive Lagrangian dynamics of flexible manipulator arms. hll. J. Robot. Res. 3, 87 101 (1984). 11. E. Schmitz and M. Ramey, Large space manipulators stud~. Martin Marietta Corporation MCR90513: AFOSRTR901031, June (19901, 12. T+ Barrows, A simplified simulation for the analysis of RMS motion with heavy payloads. Charles Stark Draper Laboratory Inc., Intralab Memorandum No. TMB842, July [1984). 13. N. Vukasovic, J. T. Celigueta, J. Garcia de Jalon and E. Bayo, Flexible multibody dynamics based on a fully cartesian system of support coordinates. Proc. 1st Int. Con/Z Dym#nics
O/ Fh'xihh,
Structures
in
Space,
Cranfield Institute of Technology, U K . , 15 [8 May 1990 (Edited by C. L. Kirk and J. L. Junkins), pp. 353 365. Computational Mechanics Publications, U.K. Springer, Berlin (1990).
14. W. L. Hallauer, Recent literature on cxpernncntai structural dynamics and control research. In ~lechanic.~ and Control +d' Lar~,,e Fh,~ihle Structures, Pro.~res,~ m Astronautics and Aeronautic.v (Edited b; J. L. JunkinsL
Vol. 129, pp. 465 489 (1990). 15. S. Fay et al., Control of flexible structures II I ( O F S ll) flight control, structure and gimbal system interaction study. Charles Stark Draper LaboratoD Inc., NASA CR172095, September (1988J. 16. C. L. Kirk and M. Moch, Dynamics and control of shuttle mast/hoopcolumn antenna system COFSII in Dynamics and Control o/'Structures #~ Space 11. Proc. of the Second Int. Conf, on Dynamics and Control of Structures in Space, Cranfield Inst. of Technol.. ILK. (6 10 September 1993). 17. S. C, C. Ht.ng and C I. Weng, Modal analysis of controlled multilink systems with flexible links and joints..41AA J. Guidance Control Dv,,:,am. 15, 634 641 (1992). 18. C. S. Bodley, A. D. Dever, A. C. Park and H. P. Frisch, A digital computer program for the dynamic interaction simulation of controls and structures (DISCOS), Vols I 4. NASA Technical Paper 1219, May (1978).