APCOM’07 in conjunction with EPMESC XI, December 3-6, 2007, Kyoto, JAPANDevelopment of a Unified Feed-Forward Control System for
Robotic Mechanisms using Finite Element Approach
Daigoro Isobe1, Kouji Yamanaka2*
1

Department of Engineering Mechanics and Energy, University of Tsukuba, 1-1-1 Tennodai Tsukuba-shi, Ibaraki
305-8573, Japan
2
Graduate School, University of Tsukuba, 1-1-1 Tennodai Tsukuba-shi, Ibaraki 305-8573, Japane-mail: isobe@kz.tsukuba.ac.jp, e0201364@edu.esys.tsukuba.ac.jp
Abstract
Dynamic equations used for feed-forward control of robotic mechanisms include
interdependent variables between the constituting links, since they are normally evaluated in relative polar
coordinates and in the dimension of torque. Accordingly, it will become highly complicated to derive
inverse dynamics of closed-loop link systems, continuously transforming systems, or of flexible linksystems. Consideration of dynamics is required to realize stable control of robotic systems, and many
researchers have tried to deal with the dynamics by improving theories and methods against each system.
Isobe, on the other hand, developed a completely new solution scheme for inverse dynamics called theparallel solution scheme, which can be commonly applied in different types of link systems such as openor closed-loop mechanisms, or ones constituted with rigid or flexible link members. The scheme isdeveloped using a finite element approach, handling the entire system as a continuum. By taking
advantage of natural characteristics of the finite element method (FEM), i.e., the capability of expressing
the behavior of each discrete element as well as that of the entire continuous system, local informationsuch as nodal forces and displacements can be calculated in parallel. It evaluates the analyzed model inabsolute Cartesian coordinates with the equation of motion expressed in dimension of force. The inverse
dynamics is calculated using a matrix form relation to the nodal forces obtained by the finite element
calculation. The matrix-form equations are divided individually into terms of force, transformation
between coordinates, and length, which makes the scheme potentially higher in applicability andexpansibility.
The scheme can not only deal with open- and closed-loop link systems independently, but it can also deal
seamlessly with those that gradually change their forms and dynamics. There is also no need to revise thebasic numerical algorithm of the scheme, regardless of the stiffness of the constituting link member.
Particularly, it is considered to be valid for link systems with elastic members, since the calculation
process of the scheme is based upon the finite element approach.
The main objective of this study is to verify the extensive ability of the scheme as a unified scheme, by
carrying out inverse dynamics calculations on several types of rigid and flexible manipulators, along withapplications to feed-forward control of various types of link systems and robotic mechanisms.
Key words: feed-forward control, link systems, robotic mechanisms, parallel solution scheme, finite
element method
INTRODUCTION
Dynamic equations derived by general schemes, such as the Newton-Euler method or the Lagrangian
method, include interdependent variables between the constituent links, since they are evaluated in relative
polar coordinates and in the dimension of torque. Accordingly, it is highly complicated to derive theinverse dynamics of closed-loop link systems, continuously transforming systems, and flexible linksystems. Nakamura and Ghodoussi proposed a systematic computational scheme of the inverse dynamics
specified for closed-loop link systems, which was derived using d'Alembert's principle [1]. Sugimoto
derived an equation of motion for closed-loop link systems by the motor algebra [2]. Nakamura andYamane developed a computational algorithm for the inverse and forward dynamics of open and closed
kinematic chains, which can be applied seamlessly to the motion of any rigid link system withoutswitching algorithms [3]. However, these methods are specified for rigid-body link systems, and users
occasionally need to revise the methods or dynamic equations when the form or elasticity of the hardware
system is changed.
An approach with a completely different viewpoint is described in this paper, which enables not only the
inverse dynamics computation of any rigid link system, but also that of flexible link systems by the sameprocedure. Isobe and Nakagawa first attempted an application of the finite element method (FEM), a
widely used computational tool for analyzing structures and fluids, to a control system of connectedpiezoelectric actuators, and achieved good control not only of the actuator itself but also of the entire
system [4]. By taking advantage of the natural characteristics of the FEM, i.e., the capability of expressing
the behaviour of each discrete element as well as that of the entire continuous system, local informationsuch as nodal forces and displacements can be calculated in parallel. After finding that the FEM can beused as a control scheme of a continuum, Isobe implemented the FEM in a solution scheme of inverse
dynamics for various types of link system [5-7]. The scheme is called the parallel solution scheme, since
the nodal forces in discrete models are calculated in parallel from a mechanical point of view. In thescheme, the analyzed model is evaluated in absolute Cartesian coordinates with the equation of motion
expressed in the dimension of force. The nodal forces are calculated incrementally in matrix form, which
does not require any revision of the overall frame, and the variables forming the frame can be revised bysimply changing the input data in the case of a physical change in the hardware system. The calculated
nodal forces are then converted into joint torques using a matrix-form equation divided into terms of force,the transformation between the coordinates, and length. The structure of the algorithm makes it seamlessin its application to different types of link system under various boundary conditions such as open- or
closed-loop link systems. The scheme can also be applied to link systems with different values of flexural
stiffness. For many years, considerable effort has been made to model [8-11] or even calculate the inverse
dynamics [12] of robotic arms with elastic members. In contrast, the parallel solution scheme only requires
a precise target trajectory that takes into account stiffness and damping.
There are two different versions of the parallel solution scheme for calculating the inverse dynamics oflink systems. One version uses linear Timoshenko beam elements with a single integration point to
explicitly express the center of gravity of links, and the other uses Bernoulli-Euler beam elements with two
integration points to express the consistent mass distribution of the links. The former version is easier tocompare with conventional dynamic equations since it explicitly expresses the center of gravity. The
latter, on the other hand, has merits in accuracy, calculation time and stability. The latter version is
described in this paper, and some numerical examples are shown for several types of link system. The
scheme is also implemented in a control system to evaluate the performance under actual control withdynamics compensation, and some control applications are shown for several types of link system and arobotic arm.
PARALLEL SOLUTION SCHEME
In the parallel solution scheme, a link system consisting of motor joints and links, as shown in Fig. 1(a), ismodeled and subdivided using finite elements. One method of modeling the system is to subdivide a link
into two linear Timoshenko beam elements, as shown in Fig. 1(b), with an intermediate node at the centerof gravity [6]. Another method is to substitute a link with a Bernoulli-Euler beam element, as shown in Fig.
1(c), considering the consistent mass distribution along the link. In this section, a type of modeling shownin Fig. 1(c) is described.
1. Formulation using Bernoulli-Euler beam elements The order of the displacement function used for
the Bernoulli-Euler beam element is higher than that of the linear Timoshenko beam element, and thus, it
can express deformation more accurately with fewer elements. In the model, the link mass is distributed
consistently, in the same manner as the displacement function, and does not require an expression for thecenter of gravity. Note that, since the deformation of the element is defined using a higher-order function,
it requires only one-element subdivision per member for infinitesimal deformation cases.

Joint (Servo motor)

Center ofgravity

Nodes

Nodes

+

+

(a) Link system

Bernoulli-Euler
beam elements

Linear Timoshenko
beam elements

Link: Motor mass

: Motor mass

: Lumped mass

: Consistent mass(b) Linear Timoshenko beam elements(c) Bernoulli-Euler beam elements

Fig. 1 Modeling of link systems using finite elements
Pn

Pn

j=i+1 FjZFiûZ
Link ij=i+1 FjY

FiûY

FiûX

Z

Pn

Link i + 1

yj=i+1 FjXz
Yx

Elemental coordinate

X：NodesGlobal coordinate

Fig. 2 Nodal forces acting in a link system with consistent mass distribution
Figure 2 shows the nodal forces (based on global coordinates) acting on the ith link in a three-dimensional
open-loop n-link system with a consistent mass distribution. The joint torque τix required around thex-elemental axis on the ith link is determined by adding the i+1th joint torque τ(i+1)x to the sum of the
moments of inertia acting on this link, and is expressed by the nodal forces Fiy and Fiφx based on elemental(or link) coordinates asúix = li (nXFj )y + Fiûx + ú(i+1)x .

(1)

j=i+1

By considering the other components around the y- and z-axes and arranging them into global coordinates(X, Y, Z) in matrix form, the joint torque vector is expressed as(2)

fún g = [Ln ][T n ]fP n g ,where {Pn} is a 6n×1 vector related to nodal force, defined as
9
8 P
n
9
>
8
>
>
j=i+1 FjX >
>
>
Pn
>
>
>
P1 >
>
>
>
>
>
>
>
>
F
jY
>
>
>
>
j=i+1
>
>
>
>
>
>
>
=
< P2 >=< Pnnj=i+1 FjZ
fP g =; (i = 1; ÅÅÅ; n)Å ; where fPi g =
>
>
>
>
FiûX
>
>
>
>
>
>
>
>
>
>
>
>
> Å >
>
>
>
>
>
>
FiûY
;
:
>
>
>
>
Pn
>
>
;
:
FiûZ

(3)
.

[Tn] is a 6n×6n transformation matrix, which is defined as(4)n[T n ] = [hn ][TGE].

where [hn] is a correction matrix between the x-y and z-x coordinate systems, which simply inverts thesigns of the components in the y-axis direction. [TnGE] is a transformation matrix between the global andelemental coordinates, which is expressed as
2

6
6
6
6
h
i 6
6
n
=6
TGE
6
6
6
6
4

where
h

Ti

i

and
h

Ai

i

=

3

T1
T2

0
T3
Å

0

"

Ai

0

0

Ai

2

cosûiXx#6= 4 cosûiXy
cosûiXz

Å

Å

Tn

7
7
7
7
7
7
7
7
7
7
7
5

(5)

,

(6a)
,
3
cosûiZx
7cosûiZy 5cosûiZzcosûiY xcosûiY ycosûiY z(6b),where φiXx, for example, represents the rotational angle between the X (global) and x (elemental)
coordinates.
[Ln] is a 3n×6n matrix related to link length, which is expressed as[Ln ] = [TÉn ][Én ] ,(7)where [TnΛ] is a transformation matrix between each elemental coordinate, and [Λn] is a matrix expressed
as
2

6
6
6
6
h
i 6
6
n
=6
É
6
6
6
6
4

where
h

Éi

i

2

3

É1
É2

0
É3

Å
0

6
= 4 li
0

0

li

0 1

0

0 0

0

0 0

Å

Å

3

7
1 0 5
0 1
0 0

Én

7
7
7
7
7
7
7
7
7
7
7
5

(8)

,

(9)
.

2. Calculation of kinematics for flexible models To calculate the inverse dynamics for flexible models,
information on target trajectories that compensate for the inertial forces acting at the links and thestiffnesses of the links is required. Therefore, a solution scheme of kinematics is developed on the basis of
the FEM. Moreover, it is combined with the parallel solution scheme to handle analyzed models
comprehensively in a single calculation process.
When inertia caused by the overall system motion is taken into account, the equation of motion at timet+∆t is derived from the principle of virtual work as follows:[M]fu°
m gt+Å t +[M]fu°
d gt+Å t +[C]fu_d gt+Å t +[K]fÅud g = fF gt+Å t ÄfRgt,(10)

where [M] is the total mass matrix, [C] is the total damping matrix, [K] is the total stiffness matrix, {F} isthe external force vector, and {R} is the internal force vector. The vectors fu°
m g and fu°
d g denote theacceleration vectors for the overall system motion components and material deformation components,
respectively, fu_d g denotes the velocity vector for material deformation components, and {∆ud} is theincremental material deformation vector at time t+∆t. By applying Newmark’s β method (integration
parameter: δ=1/2) as a time integration scheme, the acceleration and velocity vectors are calculated as1
1
1
fu_mgt Ä( Ä1)fu°fÅumgÄ
m gt
2
åÅt
åÅt
2å
,
1
1
1
fu_d gt Ä( Ä1)fu°
fu°fÅud gÄ
d gt+Å t =
d gt
åÅt2
åÅt
2å
,
1
1
1 Ä 4å
fu_mgt+Å t =fÅumgÄ( Ä1)fu_mgt Ä()fu°m gt Åt
2åÅt
2å
4å
,
1
1
1 Ä 4å
fu_d gt+Å t =fÅud gÄ( Ä1)fu_d gt Ä()fu°
d gt Åt
2åÅt
2å
4å
,
fu°
m gt+Å t =

(11a)
(11b)
(11c)
(11d)

where fu_m g is the velocity vector for the overall system motion components and {∆um} is the incremental
motion vector at time t+∆t. β is another integration parameter and a value of 1/4 is used in the examples.
Substituting Eq. (11) into Eq. (10) yields
([K]+

1
1
1
1
[M]+
[C])fÅud g = fF gt+Å t ÄfRgt +[M](fu_d gt +( Ä1)fu°d gt
åÅt2
2åÅt
åÅt
2å

1
1
1
1
1 Ä 4å
fu_mgt Ä( Ä1)fu°
Ä1)fu_d gt +(
)fu°
Ä[M](
fÅumgÄ
m gt )+[C]((d gt Åt)
åÅt2
åÅt
2å
2å
4å
.

(12)

By applying {∆um} as the input in a time integration loop of Eq. (12) and using the vectors of Eq. (11) attime t, {∆ud} at each time step can be successively obtained. The displacement vector {um} for the overallsystem motion components, and the displacement vector {ud} for the material deformation components arecalculated incrementally as
fumgt+Å t = fumgt + fÅumg ,fudgt+Å t = fudgt + fÅudg .

(13a)
(13b)

The total displacement {u} is obtained by summing both the overall system motion and material
deformation components as
fugt+Å t = fumgt+Å t + fud gt+Å t.

(14)

Final target trajectories considering the effects of stiffness and damping are obtained using Eq. (14). The
resultant forces acting on the elements can also be calculated using the obtained displacements.

NUMERICAL EXAMPLES OF INVERSE DYNAMICS CALCULATIONSIn this section, some numerical examples are shown on various types of link system to confirm theapplicability and expansibility of the parallel solution scheme.
1. Closed-loop, structure-varying, and multibranch link systems Although special attention must be
paid to the number of incremental steps, one can obtain torque curves for a closed-loop link system (seeFig. 3(a)), a structure-varying type (see Fig. 3(b)), or even for a multibranch link system (see Fig. 3(c))
without revising any part of the numerical algorithm in the parallel solution scheme [6]. This is one of themost significant advantages of using the scheme.
2. Flexible link systems A numerical test is carried out on a 5-joint link system that consists of fivedifferent stiffness links, as shown in Table 1. Damping is considered in the calculations. In Table 1, α andβ are Rayleigh’s damping coefficients. An initial trajectory for a rigid-body model is given in Fig. 4(a).

4m= 1.0 kg

Torque [N m]

2

3

0

6

m
5

2
1
4
1
Mass ○：
Joint
●：

-2
-4

Joint 1
Joint 2
Joint 3
Joint 4
Joint 5
Joint 6

-6
-8
-10

0

0.2

0.4

0.6

0.8

1

Time [s](a) Closed loop(b) Structure-varying type(c) multibranch

Fig. 3 Torque curves for various types of link system
Table 1. Link parameters of a five-joint flexible link system

Link 1
Link 2
Link 3
Link 4
Link 5

1

0.200
0.200
0.200
0.200
0.200

Y [m]

-0.5

5

Torque [N*m]

0.5

J1 2

Extra
mass [kg]

0.021
0.021
0.021
0.021
0.021

0.500
0.200
0.100
0.050
0.010

1
0.4

0

Link mass[kg]

Joint 1
Joint 2
Joint 3
Joint 4
Joint 5

0.2

α

β

0.100
0.100
0.100
0.100
0.100

2.0×10-4
2.0×10-4
2.0×10-4
2.0×10-4
2.0×10-4

Y [m]
0.4

0.5

0

0

Damping coefficients

Torque [N*m]Flexural
stiffness[Nm2]
0.53840.4376
0.29170.14590.02917

Link
length [m]Joint 1
Joint 2
Joint 3
Joint 4
Joint 5

0.2

0

3
4
-0.2

X [m]
0 0.2 0.4 0.6 0.8

-0.5
0

1

2

3

Time [s]

4

5

6

-0.2

X [m]
0 0.2 0.4 0.6 0.8 1

0

1

2

3

Time [s]456(a) Initial trajectory (b) Torque curves for rigid-body model (c) Trajectory for flexible model obtained by kinematics
computation (d) Torque curves for flexible modelFig. 4 Inverse dynamics calculation of a five-joint flexible link system with different stiffness values

The same trajectory is used and is substituted directly into the parallel solution scheme to obtain the joint
torque values shown in Fig. 4(b). A smooth joint torque curve is obtained for the rigid-body model. If
stiffness must be considered, the initial trajectory is substituted into the kinematics computation, and a new
trajectory for the flexible model is obtained, as shown in Fig. 4(c). Some bending deformation can beobserved in the figure. The new trajectory is then substituted into the parallel solution scheme in the samemanner as that for the rigid-body model, and torque curves for the flexible model are obtained, as shown in
Fig. 4(d). In contrast to the rigid-body model, vibration can be observed in the torque curve of the flexible
model owing to the flexural vibration that occurs in the model. It is confirmed that the torque curves
vibrate along with those of the rigid-body model at the natural frequency of the whole system. The results
lead us to conclude that the solution scheme of kinematics developed by the FEM can handle systems withlinks of various stiffness values.
APPLICATIONS TO FEEDFORWARD CONTROLThe parallel solution scheme is implemented in the control system shown in Fig. 5, and it is applied tofeedforward control experiments on a flexible link system and a robotic arm.
1. Flexible link system Some feedforward control experiments are carried out on a two-joint flexible linksystem, as shown in Fig. 6, which contains two links made of different materials: Link 1: stainless steel,
Link 2: acryl (see Table 2 for details). Direct-drive DC motors (Maxon F2260 for Joint 1, RE40 for Joint 2)
are used at link joints to maximize the effect of the dynamics. Other details of the experimental setup are as
follows: servo amplifiers: Maxon 4-DC-Q, encoders: Maxon HEDL554, D/A board: Interface PCI-3340,
encoder board: Interface PCI-6201E, control PC: Pentium III 600 MHz processor, 512 MB memory, OS:
Windows 2000, and sampling rate: 10 ms. The control torque for feedforward and feedback compensation(τFF+FB) is calculated using(15)úF F +F B = úmotor + úI:D: + úfeedback ,

where
(16)

úmotor = J q°d + Dq_d + fc .

Here, J is the moment of inertia of the motors, D is the viscosity coefficient of the motors, fc is the dynamic
friction force, τI.D. is the computed inverse dynamics obtained by the parallel solution scheme, and q_d andq°
d are the target angular velocity and angular acceleration, respectively. D and fc are identified beforehandby simple experiments. τfeedback is the proportional-integral (PI) feedback torque, in this case, which is
obtained using(17)úfeedback = KP (qd Ä q) + KI Ü(qd Ä q) ,Joint 1

Joint 2
Link 2
Link 1

Fig. 5 Outline of the control system

Fig. 6 Two-joint flexible link system

Table 2. Link parameters of a two-joint flexible link system
Parameter
Link length [m]Young’s modulus [GPa]Moment of inertia [m4]Link mass [kg]Extra mass [kg]Link 1
0.25
200
2.92×10-12
6.83×10-2
5.81×10-1

Link 2
0.25
4.38
7.72×10-12
6.20×10-3
4.20×10-2

Fig. 7 Control angle for vibration control

where q and qd are the actual angle acquired from the attached encoders and the target angle, respectively.
KP and KI are the feedback gain for the angle and the integrated value, respectively. Values of 3.5 and0.017 are used in the feedforward and feedback control, respectively.
Sensorless, model-based vibration control is adopted to reduce the vibration of the link head and to verify
the accuracy of the calculated trajectories and torque values. The general concept of the vibration control
is to rotate the base joint by a control angle θ, as shown in Fig. 7, at the instant when the link is maximallydeformed after stopping. The timing of the rotation is determined by predicting the time when the speed ofthe link head becomes zero from the calculated trajectory. The base joint is rotated until the amplitude issufficiently reduced. In this case, the vibration control is adopted only for Joint 2 for simplicity. Joint 1 iscontrolled by the trajectory without considering any oscillating motion.
The given task is a simple rotational motion, as shown in Fig. 8(a), which stops after 2 s. A trajectory forthe flexible model with vibration control is obtained as shown in Fig. 8(b), from the kinematics
computation. The vibration that occurs after the stop is markedly reduced. Figures 9(a) and 9(b) show acomparison of the computed and actual torques. The computed torques (τI.D.) are in good agreement withthe actual control torques using feedforward and feedback torques (τFF+FB) compared with the control

0.4

0.4

0.2

0.2

0.05

J2

0 J1

0
-0.2

-0.2

-0.4

-0.4

-0.6

Y [m]

0.2

0

-0.5

0

0.4

0.2

0

1

0.4

3

4

Time [s](a) Joint 1

Fig. 8 Initial and obtained trajectories-0.05

Computed torque (τI.D. )Control torque (τFB )Control torque (τFF+FB )
0

5

1

0.4

X [m]

2

0.2

1

0

Joint 20
1

23Time [s]4(a) Rotational angle

5

0

1

2

3

Time [s]4(b) X-coordinate of link headFig. 10 Control results

3

Time [s]Computation
Experiment (FF+FB)
Joint 102(b) Joint 2

Fig. 9 Torque curves

TargetFB
FF+FB

3

Rotational angle [rad]

2

0

-0.1

Computed torque (τI.D. )Control torque (τFB )Control torque (τFF+FB )

X [m]

X [m]
0

Torque [N*m]

0.6

Y [m]

Torque [N*m]

0.6

5

4

5

torques using only feedback torques (τFB). The complicated motion required to restrain vibration isrealized in the rotational angle (Fig. 10(a)), and the computed location of the link head also agrees wellwith the experimental result (Fig. 10(b)). These results indicate the high accuracy of the trajectories andtorque values obtained using the parallel solution scheme.
2. Seven-DOF robotic arm system Finally, some feedforward control experiments using the parallel
solution scheme are carried out on a robotic arm, as shown in Fig. 11, to verify the scheme’s practicability
for actual robotic systems. The robotic arm [13] consists of 7-DOF joints actuated using an air cylinder,
and is very light compared with an arm using 7 DC motors but has sufficient stiffness to act as anendoskeleton. The linear motion of the air cylinder is transformed to rotational motion using a screw
mechanism attached to the inside of the upper arm and forearm. Input torque values are converted to airpressure values by a linear relationship. The assumed mass distribution of the mechanism is shown in Fig.
12, and it is modeled and computed explicitly using the parallel solution scheme considering the center ofgravity [6]. PID feedback control is adopted, and each feedback gain is fixed to specific values throughout
the experiments. One of the advantages of adopting the feedforward control is that one can restrain the
feedback gains to small values. A simple 2-DOF motion is given to Joints 2 and 3, and the control results
obtained are shown in Fig. 13. The computed torque curves are shown in Fig. 14. Slight delays comparedwith the target motion can be observed in the results, which originate from the mechanical delay caused bycompressing the air cylinder. However, the positive effects of the feedforward control are apparent,
particularly at the initial stage and in the overall stability. Note that although the motions are very simple,
the dynamics of the system is regarded as that of a 7-DOF system, which may be very complicated tohandle using conventional dynamic equations.

Fig. 11 Robotic arm with air cylinder usedas endoskeleton [13]5050target
FB
FF+FB

20

J2

10
0

30

J3
20
10

0.5

1
Time[s](a) Joint 2

1.5

2

1

0.5

target
FB
FF+FB

0
0

Torque[N*m]30

Joint2
Joint3

1.5

40
Rotational angle[degree]

40
Rotational angle[degree]

Fig. 12 Modeling of the robotic arm

0

0.5

1
Time[s](b) Joint 3

Fig. 13 Control results for 2-DOF motion

1.5

2

0

0

0.2

0.4
0.6
Time[s]0.8

Fig. 14 Torque curves

1CONCLUSION
A completely new solution scheme for calculating the inverse dynamics of robotic systems is described inthis paper with various numerical examples and experimental results. The scheme derives nodal forces inparallel and converts them to a joint torque, and can easily be applied to many types of link system undervarious boundary conditions or with various stiffness values. The basic numerical algorithm of the scheme
requires no revision for different cases. The simple algorithm enables the scheme to be more easily
implemented in a program or for its functionality to be enhanced, even for nonexpert users. Theseadvantages cannot be realized by conventional schemes based upon generally used dynamic equations.
The parallel solution scheme may also be used to achieve stability and smoothness in robotic motion
control. The numerical examples and control results reveal the possibility of using the scheme as a unified
scheme, independent of the system configuration and physical parameters of link systems.
Acknowledgements This work was supported by Grant-in-Aid for Scientific Research (C) MEXT and theElectro-Mechanic Technology Advancing Foundation.
REFERENCES[1] Y. Nakamura, M. Ghodoussi, Dynamics Computation of Closed-Link Robot Mechanisms withNonredundant and Redundant Actuators, IEEE Trans. Robotics and Automation, 5(3), (1989),
294-302.
[2] K. Sugimoto, Dynamic Analysis of Closed Loop Mechanisms on the Basis Vectors of Passive Joint
Axes, J. Robotic Systems, 20(8), (2003), 501-508.
[3] Y. Nakamura, K. Yamane, Dynamics Computation of Structure-Varying Kinematic Chains and Its
Application to Human Figures, IEEE Trans. Robotics and Automation, 16(2), (2000), 124-134.
[4] D. Isobe, H. Nakagawa,. A Parallel Control System for Continuous Architecture Using Finite Element
Method, J. Intel. Mat. Syst. Structures, 9(12), (1998), 1038-1045.
[5] D. Isobe, D. Imaizumi, Y. Chikugo, S. Sato, A Parallel Solution Scheme for Inverse Dynamics and Its
Application in Feed-forward Control of Link Mechanisms, J. Robotics and Mechatronics, 15(1),(2003), 1-7.
[6] D. Isobe, A Unified Solution Scheme for Inverse Dynamics, Advanced Robotics, 18(9), (2004),
859-880.
[7] D. Isobe, A. Yagi, S. Sato, General-Purpose Expression of Structural Connectivity in the Parallel
Solution Scheme and Its Application, JSME Int. J. Series C, 49(3), (2006), 789-798.
[8] W. H. Sunada, S. Dubowsky, On the Dynamic Analysis and Behavior of Industrial Robotic
Manipulators with Elastic Members, ASME J. Mechanisms, Transmissions, and Automation in
Design, 105, (1983), 42-51.
[9] E. Bayo, A Finite-Element Approach to Control the End-Point Motion of a Single-Link Flexible Robot,
J. Robotic Systems, 4(1), (1987), 63-75.
[10] A. A. Shabana, Dynamics of Flexible Bodies Using Generalized Newton-Euler Equations, ASME J.
Dyn. Sys. Measurement Control, 112, (1990), 496-503.
[11] R. J. Theodore, A. Ghosal, Comparison of the Assumed Modes and Finite Element Models for
Flexible Multilink Manipulators, Int. J. Robotics Research, 14(2), (1995), 91-111.
[12] H. Asada, Z. D. Ma, H. Tokumaru, Inverse Dynamics of Flexible Robot Arms: Modeling andComputation for Trajectory Control, ASME J. Dyn. Sys. Measurement Control, 112, (1990),
177-185.
[13] Y. Koyabu, I. Kawabuchi, K. Hoshino, Control of Humanoid Robot Arm with Air Cylinders as
Endoskeletons, Trans. Inst. Electronics, Information and Communication Engineers A, J88-A(11),(2005), 1318-1325, in Japanese.