Applied Bionics and Biomechanics
Vol. 6, No. 1, March 2009, 57–72

On control of reaching movements for musculo-skeletal redundant arm model
Kenji Taharaa,c,∗ , Suguru Arimotob,c , Masahiro Sekimotob and Zhi-Wei Luoda Organization for the Promotion of Advanced Research, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan;b Research Organization of Science and Engineering, Ritsumeikan University, 1-1-1Nojihigashi, Kusatsu 525-8577, Japan; c RIKEN-TRI

Collaboration Centre for Human-Interactive Robot Research, Riken, 2271-130 Anagahora, Shimoshidami, Moriyama-ku, Nagoya
463-0003, Japan; d Department of Computer Science and Systems Engineering, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe
525-8577, Japan

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

(Received 24 June 2008; final version received 1 February 2009)This paper focuses on a dynamic sensory-motor control mechanism of reaching movements for a musculo-skeletal redundant arm model. The formulation of a musculo-skeletal redundant arm system, which takes into account non-linear muscle
properties obtained by some physiological understandings, is introduced and numerical simulations are perfomed. The nonlinear properties of muscle dynamics make it possible to modulate the viscosity of the joints, and the end point of the arm
converges to the desired point with a simple task-space feedback when adequate internal forces are chosen, regardless of
the redundancy of the joint. Numerical simulations were performed and the effectiveness of our control scheme is discussed
through these results. The results suggest that the reaching movements can be achieved using only a simple task-space feedback scheme together with the internal force effect that comes from non-linear properties of skeletal muscles without any
complex mathematical computation such as an inverse dynamics or optimal trajectory derivation. In addition, the dynamic
damping ellipsoid for evaluating how the internal forces can be determined is introduced. The task-space feedback is extended to the ‘virtual spring-damper hypothesis’ based on the research by Arimoto et al. (2006) to reduce the muscle output
forces and heterogeneity of convergence depending on the initial state and desired position. The research suggests a new
direction for studies of brain-motor control mechanism of human movements.
Keywords: redundancy; musculo-skeletal model; sensory-motor control mechanism

1. Introduction
The natural movements of a human seem to be smooth,dexterous and sophisticated compared with today’s robots.
Until now, these excellent movements have attracted many
physiologists, kinesiologists and robotics researchers who
hope to understand and extract them. In paticular, obtaining and mimicking these natural movements in robotics has
become a hot topic, especially for robots that will be usedin close proximity to humans. Obtaining the fundamental
elements of how the central nervous system (CNS) generates control strategies to perform these kinds of movements is one of the ultimate objectives and dreams for understanding the nature of human movements. In robotics,
many human movements such as object manipulations,
walking, etc., have been mimicked. However, achieving
highly natural movements, like human behaviour, is still
difficult. There are many reasons for these difficulties. One
is the existence of ill-posed problems induced by the joint
and muscle redundancies that a human body intrinsicallypossesses (Bernstein 1935). Until now, many solutions forovercoming these ill-posed problems by introducing various optimisation criteria have been proposed not only inphysiology but also in robotics (Hogan 1984a; Hollerbach∗

Corresponding author. Email: tahara@ieee.org

ISSN: 1176-2322 print / 1754-2103 online
C 2009 Taylor & Francis
Copyright 
DOI: 10.1080/11762320902789848
http://www.informaworld.com

and Suh 1987; Katayama and Kawato 1993; Kawato et al.
1987, 1990; Nakamura 1991). In general, these optimisations require large calculation costs to solve complicatedoptimal functions as well as an inverse dynamics, which
would be very difficult for a human to perform in real time.
Indeed, do humans perform actually these complex calculations every time they move?
In this paper, we suggest that the sensory-motor control mechanism, which enables many kinds of dexterous
human movements by the CNS, may not be so complicatedand may even use simple task-space feedback by utilising
intrinsic non-linear mechanical characteristics of skeletal
muscles. Morasso (1981) observed human reaching
movements and suggested that the central commands inhuman movements may be coordinated not in joint space
but in task space. We focus on a simple reaching task
by using a three-link planar arm model driven by nine
redundant muscles so as to mimic the construction ofa human arm in a simple way. We introduce non-linear
muscle characteristics that correspond to the findings ofphysiological studies (Hill 1938; Mashima et al. 1972)
and propose a simple task-space feedback control scheme
for the arm model taking into account internal forcesDownloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

58

K. Tahara et al.

induced by redundant muscles. These internal forces canbe basically generated by the co-contraction of agonist andantagonist muscles. Arimoto et al. (2005) suggested thatthe performence of the convergence to the desired position
for human-like reaching movements by using a joint
redundant robotic arm depends upon the damping shapingof each joint. Unfortunately, humans cannot directly
modulate the rotational joint viscosity because each joint
is driven by linear movement of antagonistic muscles. How
do humans modulate the rotational joint viscosity? Hogan(1984b) remarked and Gribble et al. (2003) observedthat the co-contraction of agonist and antagonist muscles
contributes to the modulation of the joint impedance. In
our control scheme presented here, the co-contraction ofantagonistic muscles to invoke the internal forces comingfrom the redundant muscles can modulate the damping
shaping of each joint if we introduce non-linear muscle
dynamics based on Mashima et al.’s model (1972).
Basically, a redundant actuation system such as a
musculo-skeletal arm can generate internal forces. These
internal forces belong to the null-space of the muscle spacewith respect to the joint space and thereby they are independent of a principal part of the joint torques in relation
to generation of movements of the arm end point towards
the target. However, it causes an ill-posed problem of how
to determine the internal forces uniquely. Until now, many
optimisational methods to determine the internal forces
uniquely have been proposed by introducing some sort ofartificial indices (Delp and Loan 2000; Rasmussen et al.
2003). Unlike these methods, our proposed method does
not need any optimisation criteria, because the internal
forces can affect only part of the joint torques irrelevant
to motion of the arm end point towards the target in thetask space. Nevertheless, human-like slightly curved trajectories of the end point can be realised by adequately
adjusting the internal forces through the null space withrespect to the joint space. This point quite differs from thetraditional optimisation methods.
In our previous work (Tahara et al. 2005), we treated
the two-joint six-muscle model which has only muscle redundancy, that is, joint redundancy was not taken into consideration. In this paper, we take into account not only
the muscle redundancy but also the joint redundancy byintroducing a three-joint nine-muscle planar arm model.
The crucial difference between these two cases is the convergence analysis that uses Lyapunov’s direct method toprove that the closed-loop dynamics of the joint redundant
model cannot be applied. To prove the stability and convergence of the closed-loop dynamics under the existences ofboth muscle and joint redundancies, we introduce the novel
concept of ‘stability on a manifold’ originated by Arimoto(2005).
In what follows, we formulate the kinematics and dynamics of the three-joint planar arm model driven by six
mono-articular muscles and three bi-articular muscles governed by non-linear behaviour. The stability and conver-gence of the overall system is given by introducing a stability theory on a manifold. We report some numerical simulations to demonstrate that the simple sensory-motor control scheme makes it possible to make reaching movementseven though both joint and muscle redundancies are taken
into consideration. To discuss and evaluate how to choose
the internal forces, we introduce a dynamic damping ellipsoid. In addition, the task-space feedback manner is extended to the ‘virtual spring-damper hypothesis’ proposed
by Arimoto et al. (2005) to reduce the muscle output forces
and the heterogeneity of convergence, which depends on
the initial state and desired position of the arm. Finally,
these results suggest a direction for studies of the brainmotor control mechanism of human movements.
2.

Musculo-skeletal redundant system

We introduce a musculo-skeletal redundant arm model thatconsists of three serial links, six mono-articular muscles,
and three bi-articular muscles to imitate the configurationof a human arm. This model is shown in Figure 1. The threebi-articular muscles enhance the end-point stiffness.

2.1. Kinematics of muscle-joint space
Assume that the overall movements of the system are limited to a horizontal plane so that gravity does not affectFigure 1. Three-link musculo-skeletal arm model with six
mono-articular muscles and three bi-articular muscles.

Applied Bionics and Biomechanics

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

movements at all. Any friction such as Coulomb’s friction
and static friction is ignored in our model. All skeletal muscles included in the arm can only be linearly contracted
and are not curved on the way. Also the mass transfer ofthe muscles during contraction is omitted because its effect on arm movements is not so crucial. The physical parameters and variables are shown in Figure 1. The lengthsof the muscles l i (θ) ∈ R9 are formulated as the distances
between the insertion points of the muscles. They can be
expressed as follows:l(θ) = (l1 , l2 , l3 , l4 , l5 , l6 , l7 , l8 , l9 )T⎞⎛1(r1 2 + s1 2 + 2r1 s1 cos θ1 ) 2
⎟
⎜
1
⎟
⎜
(r2 2 + s2 2 − 2r2 s2 cos θ1 ) 2
⎟
⎜
1
⎟
⎜
(r3 2 + s3 2 + 2r3 s3 cos θ2 ) 2
⎟
⎜
⎟
⎜
1
⎟
⎜
2
2
(r4 + s4 − 2r4 s4 cos θ2 ) 2
⎟
⎜
⎟
⎜
⎟
⎜ (w51 2 + w52 2 + L1 2 + 2w51 L1 cos θ1
⎟
⎜
1 ⎟
⎜
⎜ +2w52 L1 cos θ2 + 2w51 w52 cos (θ1 + θ2 )) 2 ⎟⎟ ,(1)=⎜⎟⎜ (w61 2 + w62 2 + L1 2 − 2w61 L1 cos θ1
⎟
⎜
⎜
1 ⎟
⎜ −2w62 L1 cos θ2 + 2w61 w62 cos (θ1 + θ2 )) 2 ⎟
⎟
⎜
1
⎟
⎜
⎟
⎜
(r7 2 + s7 2 + 2r7 s7 cos θ3 ) 2
⎟
⎜
1
⎟
⎜
2
2
(r8 + s8 − 2r8 s8 cos θ3 ) 2
⎟
⎜
⎟
⎜
⎟
⎜ (w91 2 + w92 2 + L2 2 + 2w91 L2 cos θ2⎠⎝
1+2w92 L2 cos θ3 + 2w91 w92 cos (θ2 + θ3 )) 2
where li (i = 1−9) are the lengths of the muscles, andθi (i = 1−3) are the joint angles. Also, ri (i = 1−4, 7−8), and si (i = 1−4, 7−8) are the distances between
the centre of each joint and the insertion point of eachmono-articular muscle, and wi (i = 51, 52, 61, 62, 91, 92)
are those of bi-articular muscles. Taking the time derivative of Equation (1) yieldsl̇ = W T θ̇,(2)where W T ∈ R9×3 stands for the Jacobian matrix fromjoint space to muscle space. It is known that the relation between contractile forces of the muscles Fm ∈ R9 and joint
torques τ ∈ R3 can be expressed from the principle of virtual work as
τ = W Fm .

(3)

In this study, we assume that the Jacobian matrix W is ofrow full-rank during movement. Therefore, we can take theinverse of Equation (3) as
Fm = W + τ + (I 9 − W + W )ke ,

(4)

T −1
9×3
where W + = W T (W W
the pseudo ) ∈ +R  signifies
inverse matrix of W , I 9 − W W ke ∈ R9 stands for the

59

null-space of W , ke ∈ R9 is an arbitrary vector, and I 9 ∈R9×9 denotes an identity matrix. The physical meaning of
the second term of the right-hand side is the internal force
generated by the redundant muscles. This relation does notmean an optimisation, and it expresses only decompositionof the muscle force space with respect to the joint space
into two spaces, one is the image space and the other is thekernel space.

2.2. Kinematics of joint-task space
The end-point position of the arm model given by x =(x, y)T in Cartesian coordinates is expressed as
L1 cos θ1 +L2 cos(θ1 + θ2 ) + L3 cos(θ1 + θ2 + θ3 )x=
. (5)L1 sin θ1 +L2 sin(θ1 + θ2 ) + L3 sin(θ1 + θ2 + θ3 )Taking the time derivative of Equation (5) yieldsẋ = J θ̇ ,(6)

where J ∈ R2×3 stands for the Jacobian matrix from joint
space to task space. The relation between joint torquesτ ∈ R3 and output forces of the end point in Cartesian coordinates Fx ∈ R2 can be represented asτ = J T Fx .

(7)

Substituting Equation (7) into Equation (4) yields the relation between contractile forces of the muscles Fm andoutput forces of the end-point Fx . It is given by
Fm = W + J T Fx + I 9 − W + W ke .

(8)

2.3. Non-linear muscle model
It is known that the skeletal muscles are governed bystrong non-linearity. Here, we introduce a tangible nonlinear muscle model that is based on several physiological results. In the modeling of the muscles, we assume thatthe masses of all muscles are included in the masses ofeach link. The model of muscle is basically according to
the two-element model which is composed of contractile
element (CE) and series elastic element (SE). In this paper, we want to focus on the dumping effect on each joint
generated by CE, and thereby we introduce quite strong assumption that the SE part can be approximated as a rigid
element which indicate that its elastic coefficient becomes
infinity. Therefore, hereafter, we only consider the CE partin the muscle model. In physiology, the force–velocity relation of the skeletal muscle proposed by Hill (1938), whichis expressed as a simple hyperbolic equation, is known as

60

K. Tahara et al.

one of the most fundamental muscle features. This equation is(9)

where f is the output tensile force of the muscle, l˙ isthe contractile velocity of the muscle, f0 is the maximum
isometric tensile force of the muscle which depends on
its length, a is the heat constant and b is the rate constant of energy liberation. The shortening directions off and l˙ are defined as positive. In addition, it is knownin physiology that a skeletal muscle exerts more tensile
force in the lengthening phase than in the shortening phase.
Mashima et al. (1972) determined the above feature specifically through some experiments. By considering this feature, we can modify Equation (8) as follows:
bf0 − a l˙
α
˙ +b
|l|
˙ =
f (α, l)⎪bf0 − (2f0 + a)l˙
⎪
⎪
α
⎩
˙ +b
|l|

if l˙ ≥ 0(10)
if l˙ < 0,

⎧

0.9l0
0.25 ˙
⎪
⎪
ᾱ −
ᾱ l
if l˙ ≥ 0⎪⎪˙
⎨ 0.9l0 + |l|
0.9l0
⎪
⎪
⎪
⎪
⎩

(11)

2.25 ˙
0.9l0
ᾱ −
ᾱ l
if l˙ < 0˙
0.9l0
0.9l0 + |l|

Since a skeletal muscle intrinsically owns some damping factor, which comes from its component and material,
we introduce a muscle intrinsic viscosity c0 > 0 independently from ᾱ. Hence, the muscle dynamic model used inthis study can be expressed as˙ = p ᾱ − (ᾱc + c0 ) l˙ ,f (ᾱ, l)

100
80
60
40
20
0
-0.1

-0.05
0
0.05
Muscle contractile velocity (m/s)(12)

as the contractile velocity of the muscle increases and viceversa. Now, it should be remarked that in Equation (12) the
parameter p depending on l˙ satisfies the inequality Equation 0 < p ≤ 1 as long as l˙ is upper bounded. Eventually,
the overall muscle dynamics in this study can be expressed
as follows:
Fm = P ᾱ − P { A(ᾱ)C + C0 } l̇⎡
Fm = (f1 , f2 , . . . , f9 ) ∈ R9
⎢ P = diag(p , p , . . . , p ) ∈ R9×9
1
2
9
⎢
⎢
⎢ ᾱ = (ᾱ1 , ᾱ2 , . . . , ᾱ9 )T ∈ R9⎢
⎢ A(ᾱ) = diag(ᾱ , ᾱ , . . . , ᾱ ) ∈ R9×9
1
2
9
⎢
⎢
⎣ C = diag(c1 , c2 , . . . , c9 ) ∈ R9×9

p=

0.9l0
,
˙
0.9l0 + |l|c=

if l˙ ≥ 0

⎪
⎪
2.25
⎪
⎩
> 0,
0.9l0

if l˙ < 0.

The force–velocity curve of the muscle property expressedby Equation (12) is shown in Figure 2. As the figure shows,the contractile force of the muscle model reduces gradually(13)

C 0 = diag(c01 , c02 , . . . , c09 ) ∈ R9×9
It is remarkable that in Equation (13), the diagonal matrixA(ᾱ) is constructed by the elements of the control input
vector ᾱ (Tahara et al. 2005, 2006). This comes from thenon-linearity of the muscle model. Also, the skeletal muscle generates only a contractile force, so we define the control signal ᾱ, which is a saturated function so as to keep asemi-positive value, as shown in Figure 3. Therefore, A(ᾱ)
is semi-positive definite. Both damping matrices C and C 0

where
⎧
0.25
⎪
⎪
> 0,
⎪
⎨ 0.9l0

0.1

Figure 2. Force-velocity curve of the modified Hill’s muscle
property model, which is expressed by Equation (12).

where 0 ≤ α ≤ 1 is the muscle activation level, the coefficients a and b are determined experimentally in Mashima
et al. (1972) as a = |0.25 × f0 |, b = |0.9 × l0 | and l0 is theintrinsic rest length of the muscles. Now, we define the control input to the muscles as ᾱ = f0 × α, in order to set thedimension of the input as that of the muscle force. Therefore, Equation (10) can be represented as follows:˙ =f (ᾱ, l)120

Muscle input (N)Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

⎧
⎪
⎪
⎪
⎨

α = 20.0 (N)α = 60.0 (N)α = 100.0 (N)140
Muscle output force (N)(f + a)(l˙ + b) = b(f0 + a)

160

+

α

0
Time (s)–Figure 3. Saturated function of the input ᾱ.

Applied Bionics and Biomechanics

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

are also positive definite. In this paper, the positive coefficient c0 , which is the intrinsic viscosity of muscle, is defined as c0 = 0.2 by considering physiological aspects.

2.4. Control input
In physiological studies of the brain-motor control mechanism in human movements, several control strategies forskeletal muscles have been proposed since Fel’dman first
proposed a spring-like model called the λ-model (Fel’dman
1966, 1986) which exerts muscle force by modulating
the threshold of tonic stretch reflex about muscle length.
Mclntyre and Bizzi (1993) also proposed the α-model
which modulates spring-like characteristics of the muscles
by α neuron. These models are based on the equilibrium
point (EP) hypothesis proposed by Fel’dman and observed
by Bizzi et al. Many studies concerned with multijoint
reaching movement, which were based on the EP hypothesis, have been done (Bizzi et al. 1976, 1992; Flash 1987;
Flash and Hogan 1985). Recently, Arimoto et al. (2005)
have proposed a virtual spring hypothesis that a given
potential energy in task space can generate a spring-likeforce on the end point and distribute it to each joint torquethrough the transpose of the Jacobian matrix J T . Thistask-space control scheme originated from Takegaki and
Arimoto (1981). In developmental psychology, it is said
that a potential function that emanates neuro-motor signals
emerges from intentional incentive of a voluntary movement (Killeen 1992). The virtual spring hypothesis has emphasised that such a potential-generating neuro-motor signals should be coordinated not in the muscle or joint spaces
but in the task space in terms of the distance between apresent position of the end point and a desired one such as
1
x T K x. This potential is eventually distributed to each2
joint through the transposed Jacobian matrix J T which implies the joint configuration of a human body. In the present
study, we introduce a simple task-space feedback control
scheme with internal forces, which basically follows thevirtual spring hypothesis.
Now, let us consider the diagonal matrix P expressedby Equation (13). Since the boundaries of its components
pi are given by 0 < pi ≤ 1, P is positive definite. Hence,+
we can obtain a pseudo-inverse matrix W̄ = (W P)+ as
long as the muscle Jacobian matrix W does not degenerate
during movement. The control input to the muscles ᾱ isexpressed by++ᾱ = − W̄ J K x + (I 9 − W̄ W̄ )ke ,T(14)

where W̄ = W P, and x = x − xd is the end point position error between the present position x ∈ R2 and thedesired position xd ∈ R2 in the task space. The positive
definite diagonal matrix K = k I 2 ∈ R2×2 denotes the position feedback gain, and k signifies the positive stiffness61

coefficient of the virtual spring. It should be noted thatin Equation (14), the first term of the right-hand side denotes an artificial potential energy in task-space and produces spring-like forces to each muscle through two Ja+
cobian matrices J T and W̄ and the second term stands
for the internal forces generated by the redundant muscles,
where ke ∈ R9 is an arbitrary vector. This control signals
expressed by Equation (14) can be interpreted more specifically as that K x comes from intentional incentive of amovement which is generated in the motor cortex, and J T
implies the joint configuration of the arm relative to the arm++
end-point, W̄ and (I 9 − W̄ W̄ )ke come from the configuration of the related muscles relative to the joints.

2.5.

Dynamics of three-link planar arm model

The dynamics of a three-link planar arm model can be described by Lagrange’s equation of motion as
1H(θ )θ̈ +Ḣ(θ) + S(θ̇ , θ ) θ̇ = τ ,2(15)

where H(θ) ∈ R3×3 is the inertia matrix, which is symmetric and positive definite, and θ̈ ∈ R3 , θ̇ ∈ R3 , θ ∈ R3
are angular accelerations, angular velocities and angles ofjoints, respectively. S(θ̇, θ ) ∈ R3×3 is a skew-symmetricmatrix coming from Coriolis and centrifugal forces, andτ ∈ R3 are the input torques to the joints. Substituting
Equations (2), (3), (13) and (14) into Equation (15) yields
H(θ )θ̈ +
1
Ḣ(θ ) + S(θ, θ̇ ) θ̇2

+ J T K x + W P( AC + C 0 )W T θ̇ = 0. (16)Equation (16) shows the overall closed-loop dynamics,
which is expressed in the joint coordination system.

3. Stability of the closed-loop dynamics
This section illustrates the convergence of the overall system by introducing the concepts of ‘stability on a manifold’
and ‘transferability to a submanifold’. Now, it is convenient
to rewrite Equation (16) as
H(θ )θ̈ +
1
Ḣ(θ) + S(θ, θ̇ ) θ̇2

+ J T K x + W P ( AC + C 0 )W T θ̇ = u,(17)where u = 0 are the input to the closed-loop dynamics ofEquation (16). Taking the inner product between the input
as Equation (17) and the output as the joint angular velocity

62

K. Tahara et al.

θ̇ yields

Definition 2 (Transferability to a submanifold)Tθ̇ u =dT
E(t) + θ̇ W P ( AC + C 0 )W T θ̇,
dt

(18)

1 T
k
θ̇ H θ̇ + x2 ≥ 0.
2
2

(19)

where

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

E(t) =

Here, we should consider the damping matrix W P( AC +C0 )W T ∈ R3×3 such that the control inputs ᾱi , which arecomponents of the diagonal matrix A, are defined as a saturated function (see Figure 3) to be semi-positive. Therefore, the matrix ( AC + C0 ) is positive definite because
A ≥ 0, C > 0, C0 > 0 and P > 0. Hence, the damping
matrix W P ( AC + C0 )W T is positive definite as long as
the Jacobian matrix W is of row full-rank.
Now, integrating Equation (18) over time t ∈ [0, T )yields t

If for the reference state (θ 0 , 0) ∈ M1 , there exist a constant ε1 > 0 and another constant r1 > 0 (r1 < r0 ) suchthat any solution trajectory starting from an arbitrary initialstate contained in N 6 (δ(ε), r1 ) stays on N 6 (ε1 , r0 ) and converges asymptotically as t → ∞ to a point on a submanifold M1 ∩ N 6 (ε1 , r0 ), the reference state (θ 0 , 0) is said tobe transferable to a submanifold of M1 .
The physical meaning of Definition 1 is that the overall
energy of the system does not increase during movement.
In addition, the meaning of Definition 2 is that the desired
state (x = 0, θ̇ = 0) is guaranteed while the system isstable.
Taking the inner product between Equation (16) andθ̇ + βk J T x yields

1 T
k
T
2
θ̇ H θ̇ + x + θ̇ W P( AC + C 0 )W T θ̇
2
2


1
T
T
2
T
2
Ḣ + S θ̇
+βk  J x +βkx J H θ̈ + βkx J
2

d
dt



+βkx T J W P( AC + C 0 )W T θ̇ = 0,θ̇ udσ = E(t) − E(0)

0

+

(21)

T

 t

where β is a positive number. We introduce the following
equation.

Tθ̇ W P ( AC + C 0 )W T θ̇ dσ ≥ −E(0). (20)0This inequality means that the closed-loop dynamics ofEquation (16) satisfies the passivity condition. However,it should be noted that this system is a joint redundant system because the control input of task-space feedback requires only an end-point position (x, y) although the arm
consists of three joints. Then, the scalar function E(t) is
non-negative definite, but not positive definite for all statevariables of (θ, θ̇ ). Therefore, to discuss the convergenceof this redundant system, we introduce the novel concept of‘stability theory on a manifold’ proposed by Arimoto et al.
(2005). In this proof, we consider that the one-dimensional
submanifold M1 on the spectrum S 6 ∈ R6 such that thetransposed Jacobian matrix J T ∈ R3×2 is of column fullrank.
Let us introduce the following definitions to prove thestability and convergence of the closed-loop dynamics.

x T J H θ̈ =dT
x T J H θ̇ − θ̇ J T J H θ̇dt
− x T J̇ H θ̇ − x T J Ḣ θ̇ .

Substituting Equation (22) into Equation (21) yields
d
dt



1 T
k
θ̇ H θ̇ + ||x||2 + βkx T J H θ̇
2
2



T

+ θ̇ W P (AC + C0 )W T θ̇
+ βk 2 || J T x||2 + βkx T J W P( AC + C0 )W T θ̇
+ βkh(θ̇ , x) = 0,

(23)

where


1
h(θ̇, x) = x J − Ḣ + S θ̇2
TT− θ̇ J T J H θ̇ − x T J̇ H θ̇ .
Definition 1 (Stability on a manifold)
For the reference state (θ = θ 0 , θ̇ = 0) ∈ M1 , if for an arbitrarily given ε > 0, there exist a constant δ(ε) > 0 andanother constant r1 > 0 (r1 < r0 ) independent of ε suchthat a solution trajectory starting from any initial state(θ(0), θ̇ (0)) contained in N 6 (δ(ε), r1 ) in the neighbourhood of δ stays on N 6 (ε, r0 ) in the neighbourhood of ε, thereference state (θ 0 , 0) is said to be stable on a manifold.

(22)

(24)

The second, third and sixth terms of the left-hand side ofEquation (21) can be rewritten as follows:T−θ̇ W P( AC + C0 )W T θ̇ − βk 2 || J T x||2−βkx T J W P( AC + C0 )W T θ̇


1 T
k
θ̇ H θ̇ + ||x||2 + βkx T J H θ̇
= −γ
2
2

Applied Bionics and Biomechanicsγ−θ̇ W ( AC + C0 )W − H θ̇2γ− kx T βk J T J − I 2 x
2
T

where

T

where γ is a positive number. Also, the following inequality can be obtained by considering the positive definitenessof the damping matrix W P( AC + C0 )W T .
−βkx T J W P( AC + C0 )W T θ̇
≤ βkx T J W P( AC + C0 )W T J T x
βk Tθ̇ W P( AC + C0 )W T θ̇ .
4

(26)

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

Similarly, we can obtain the following inequality:
γβk T
γ αkx J H θ̇ ≤ γβkx J H J x +θ̇ H θ̇ . (27)4TT−θ̇ W P( AC + C0 )W θ̇ − βk || J x||T

In Equation (25), since Ḣ and S are in linear homogeneous with respect to the angular velocities θ̇ , and the orderof their coefficients is less than or equal to the maximum
eigenvalue of H. Also, J consists of a trigonometric function and physical parameters, and they are thereby totallybounded. Thus, as long as x(t) < x(0), there exists
a positive number ζ0 for h to satisfy the following inequality (Tahara et al. 2005, 2006):|h(θ̇, x)| ≤ ζ θ̇ 2 .

(32)

T

By substituting Equations (26) and (27) into Equation (25),
we obtain the following inequality:Tβk
W P( AC + C 0 )W T
1−f (θ̇, x) = θ̇
4
 

βk
γ
γ
1+
H θ̇ + kx T βk J TJ − I 2
−
2
2
2



−β J γ H + W P( AC + C 0 )W T J T x. (31)T

+γβkx T J H θ̇ − βkx T J W P( AC + C0 )W T θ̇ , (25)

+

63

2

T

2

−βkx T J W P( AC + C0 )W T θ̇


1 T
k
T
T
2
θ̇ H θ̇ + ||x|| + βkx J H θ̇ − θ̇
≤ −γ
2
2

 

βk
βk
γ
T
W P( AC + C0 )W −
1+
H θ̇
× 1−
4
2
2

γ
−kx T βk J T J − I 2
2
 
(28)
− β J γ H + W P (AC + C0 )W T J T x.

Now, assume that we can choose the numbers k, β, γ andthe internal force vector ke to satisfy the following inequalities.
I 2 − β 2k J T H J > 0

βk
γ
1−
W P( AC+C 0 )W T −
4
2

1+

βk
2

(33)



H > βkζ(34)γβk J T J− I 2 −β J γ H+W P( AC+C 0 )W2

T

J T >0.
(35)

By considering these inequalities, we see that the scalarfunction V (t) expressed by Equation (29) is positive, andthe inequality−βkh(θ̇ , x) − f (θ̇ , x) ≤ 0Here, we introduce the scalar function V (t) such that
1 Tkθ̇ H θ̇ + ||x||2 + βkx T J H θ̇
2
2
T 

1
=
θ̇ + βk J T x H θ̇ + βk J T x
2


k
(29)
+ x T I 2 − β 2 k J T H J x.
2

V (t) =

Substituting Equations (28) and (29) into Equation (23)yieldsd
V (t) ≤ −γ V (t) − βkh(θ̇, x) − f (θ̇ , x), (30)
dt

(36)

is satisfied. Hence, Equation (23) can be reduced todV (t) ≤ −γ V (t).
dt

(37)

Therefore, from Equation (37) we obtain
V (t) ≤ e−γ t V (0)(38)

and conclude that the scalar function V (t) converges tozero exponentially. The scalar function V (t) plays the roleof a modified energy function of E(t) in Equation (19).
Therefore, θ̇ → 0 and x → 0 when t → ∞.
Q.E.D.

64

K. Tahara et al.

Table 1. Physical parameters of the arm.

Upper arm
Forearm
Wrist

Length
(m)

Mass
(kg)

Inertia(kg·m2 )Mass centre
(m)

0.31
0.27
0.15

1.93
1.32
0.35

0.01410.01200.0010

0.165
0.135
0.075

Numerical simulation

Table 1 shows the physical parameters of the arm model,
Table 2 shows the insertions of the muscles and Table 3
shows the initial and desired positions of the end point.
In this paper, we treat two types of reaching movements:
short-range and middle-range reaching. We also determine
the internal force vector as ke = ke e in order to see theinfluence of ke on the dynamics, where ke is a positive
number and e = (1, 1, . . . , 1)T ∈ R9 . Simulation results
for the short-range reaching movement (|x| ≤ 0.3[m])with virtual spring coefficient k = 10.0 and internal force
coefficient ke = 5.0 are shown in Figure 4. The end-point
trajectory moves along a very winding path even though it
converges to the desired point. Simulation results for the
short-range reaching movement with k = 10.0 and ke =25.0 are shown in Figure 5. Here, in contrast to Figure 4,
the end-point converges to the desired point smoothly andits trajectory is close to a straight line. Figure 6 shows
the trajectories of the end-point position (x, y) of a shortrange reaching movement with the virtual spring coefficient k = 10.0 when ke is shifted from 5.0 to 30.0 by
5.0 units. These trajectories are made close to a straightline by increasing the internal force coefficient ke . Figure 7 shows the transient responses of the end-point velocity |v| = ẋ 2 + ẏ 2 of short-range reaching movementwith k = 10.0 when ke is shifted from 5.0 to 30.0 by 5.0
units. The velocity profiles are made smooth by increasing
the internal force coefficient ke . Figures 8 and 9 show thetransient responses of the end-point x- and y-position withk = 10.0 when ke is shifted from 5.0 to 30.0 by 5.0 units.
We see from these figures that the shortest settling timeof the end-point x-position is about 1.5 s when ke = 30.0,
Table 2. Insertion of each muscle.
Muscle
l1
l2
l3
l4
l5
l6
l7
l8
l9

Value (m)r1 = 0.055
r2 = 0.055
r3 = 0.030
r4 = 0.030
w51 = 0.040
w61 = 0.040
r7 = 0.035
r8 = 0.050
w91 = 0.040

Initial point (m)(x, y) = (0.0, 0.60)

Desired point (m)Short range
Middle range(xd , yd ) = (−0.20, 0.40)(xd , yd ) = (−0.40, 0.20)

and that of the y-position is about 2.0 s when ke = 20.0.
On the other hand, Figure 10 shows the simulation resultof the middle-range reaching movement with k = 10.0 andke = 5.0. The trajectory of the end point is not smooth anddescribe a more winding route than that of the short-range
reaching, and the final posture of the arm seems to be impossible for a real human. In contrast, Figure 11 shows thesimulation result of the middle-range reaching movementwith k = 10.0 and ke = 25.0. The end-point converges tothe desired point smoothly, and its trajectory becomes aquasi-straight line even though x ≤ 0.6 m holds forin middle-range reaching. Figure 12 shows trajectories ofthe end-point position (x, y) of the middle-range reaching
movement with k = 10.0 when ke is shifted from ke = 5.0 to
ke = 30.0 by 5.0 units. Like the short-range reaching movements shown in Figure 6, the end-point trajectories becomesmooth and close to straight lines as ke increases. Figure 13
shows transient responses of the end-point velocity |v| of
the middle-range reaching movement with k = 10.0 when
ke is shifted from ke = 5.0 to ke = 30.0 by 5.0 units. The
velocity profiles when ke = 5.0 and ke = 10.0 are oscillatory, and each settling time is more than 3.0 s. However,those in the case of ke = 25.0 and ke = 30.0 are smooth because of the higher ke , and also each settling time is about
1.5 s. Figures 14 and 15 show the transient responses ofthe end-point x- and y-position with k = 10.0 when ke isshifted from 5.0 to 30.0 by 5.0 units. We see from thesefigures that the fastest convergence time of the end-point

0.8
0.6
y-position (m)Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

4.

Table 3. Initial and desired points.

s1 = 0.080
s2 = 0.080
s3 = 0.120
s4 = 0.120
w52 = 0.045
w62 = 0.045
s7 = 0.220
s8 = 0.250
w92 = 0.030

0.4
0.2
0
–0.2

–0.4

–0.2
0
x-position (m)

0.2

0.4

Figure 4. Short-range reaching movement with k = 10.0 andke = 5.0.

Applied Bionics and Biomechanics

65

End-point velocity |v| (m/s)

0.8

y-position (m)

0.6
0.4
0.2
0
–0.2

–0.4

–0.2
0
x-position (m)0.2ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

1

0.5

0

0.4

0

1

2

3

4

5

Time (s)0.5ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.00ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

-0.1

-0.2

0

1

2

3

4

5

Time (s)Figure 8. Transient responses of the end-point x-position ofshort-range reaching movements with k = 10.0.

ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

0.6

y-position (m)0.6

Figure 7. Transient responses of the end-point velocity of shortrange reaching movements with k = 10.0.

x-position (m)x-position is about 2.0 s when ke = 15.0, and that of the
y-position is about 2.0 s when ke = 25.0. These simulation results suggest that the internal force coefficieint ke
are chosen according to the virtual spring coefficient k, theend point of the musculo-skeletal redundant arm converges
to the desired point smoothly without any complex mathematical optimization. The benefit of our proposed modelis that co-contraction of antagonistic muscles can modulate the damping effect in the joint space. Until now, many
control schemes for muscle redundant arm model are usedby some optimisation method to determine muscle output
force uniquely because the internal forces cannot affect tothe joint torques. Unlike these methods, in our model, internal forces can affect the joint torques, and thereby the
trajectory of the end point is changed by varying the internal force vector ke .

y-position (m)Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

Figure 5. Short-range reaching movement with k = 10.0 andke = 25.0.

0.5

0.4

0.4

0.3

-0.2

-0.1

0

x-position (m)Figure 6. Trajectories of the end-point position of short-range
reaching movements with k = 10.0.

0

1

2

3

4

5

Time (s)Figure 9. Transient responses of the end-point y-position ofshort-range reaching movements with k = 10.0.

66

K. Tahara et al.
0.8
End-point velocity |v| (m/s)

2y-position (m)

0.6
0.4
0.2
0

ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

1.5

1

0.5

0

–0.2

–0.4

–0.2
0
x-position (m)

0.2

0.4

0

2

3

5

Figure 13. Transient responses of the end-point velocity of
middle-range reaching movements with k = 10.0.

0.8

ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.00x-position (m)0.6
y-position (m)4Time (s)Figure 10. Middle-range reaching movement with k = 10.0 andke = 5.0.

0.4
0.2

–0.2

–0.4

0
–0.2

–0.6

–0.4

–0.2
0
x-position (m)

0.2

0.4

ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

0.4

1

2

3

4

5

Figure 14. Transient response of the end-point x-position of
middle-range reaching movements with k = 10.0.

ke = 5.0
ke = 10.0
ke = 15.0
ke = 20.0
ke = 25.0
ke = 30.0

0.6

y-position (m)0.60Time (s)Figure 11. Middle-range reaching movement with k = 10.0 andke = 25.0.

y-position (m)Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

1

0.4

0.2

0.2

–0.4

–0.2

0

x-position (m)Figure 12. Trajectories of the end-point position of middlerange reaching movements with k = 10.0.

0

1

2

3

4

5

Time (s)Figure 15. Transient response of the end-point y-position of
middle-range reaching movements with k = 10.0.

Applied Bionics and Biomechanics
5.

Dynamic damping ellipsoid

We introduce a dynamic damping ellipsoid, in order to discuss how the internal force vector ke is determined, and toevaluate it. The damping ellipsoid at the end point of thehuman arm has been proposed by Tsuji et al. (1994). The
ellipsoid can express the damping effect at the end pointof a human arm when the arm’s dynamics can be approximated as being almost governed by only a term related to
the velocity. It was derived from the result of measurement
and the evaluation of a real human hand’s damping effectin a quasi-static situation. Unlike the above ellipsoid, ourdynamic damping ellipsoid expresses the damping effectthat can be given to the end point of the arm by modulating
the internal force vector ke included in the control input tothe muscles. Figure 16 illustrates the dynamic damping ellipsoid at the end point of a human arm. We consider partof the control input in Equation (14), which can generate
the joint torque for only the damping effect as
Downloaded By: [Tahara, Kenji] At: 05:27 7 May 200967

τ D = −W P( AC + C 0 )W T θ̇ ,(39)

where τ D is the vector of the input torques for generating
the damping effect. Since the damping force FD at the end
point in task space can be expressed as
F D = − D ẋ,

(40)

where
D = −( J T )+ W P( AC + C 0 )W T J + ,Figure 16. Dynamic damping ellipsoid at the end point.

(41)

Figure 17. Dynamic damping ellipsoid at the end point whenk = 10.0 and ke = 25.0.

( J T )+ is the pseudo-inverse transposed Jacobian matrix of
J T , and J + means the pseudo-inverse Jacobian matrix of
J. Therefore, the exertable end-point damping force whenthe end-point velocity ẋ satisfies the inequality of  ẋ ≤ 1
can be expressed as
FD T ( D−1 )T D−1 FD ≤ 1.

(42)

This ellipsoid shows the end-point damping forces generated by the control input during movement. By considering
this ellipsoid, we can evaluate the internal forces required
to modulate the damping effect. Figure 17 shows the damping ellipsoid at the end point of the middle-range reaching movement with ke = 25.0, and Figure 18 shows thatwith ke = 5.0. We see from these figures that in the caseof a suitable internal force (ke = 25.0), the long axis of

Figure 18. Dynamic damping ellipsoid at the end point whenk = 10.0 and ke = 10.0.

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

68

K. Tahara et al.

the damping ellipsoid turns on the shoulder during movements. Therefore, the short axis of the damping ellipsoid
turns on the desired point, and the damping effect at theend-point is relatively small in the direction of the desired point during movement. Also at the beginning of
the movement, the ellipsoid is relatively large. It becomes
gradually small as the end-point velocity increases, and finally the shape of the ellipsoid becomes large again likeat the beginning shape. However, in the case of an unsuitable internal force (ke = 5.0), the damping ellipsoid isrelatively small: the long and short axes of the ellipsoid
are nearly equal. By comparing these results we can conclude that the end point can converge to the desired pointsmoothly by considering the size and direction of the dynamic damping ellipsoid during movement. The size anddirection of dynamic damping ellipsoid are governed by
the internal force vector ke and muscle contractile velocity. However in this paper, we only pointed out a relation
between the internal forces and the ellipsoid from the viewpoint of muscle and joint configurations as a posteriori reasoning, and could not disclose any dynamic relation between the internal force vector and the ellipsoid explicitlyfrom the dynamics viewpoint, because an explicit dynamics of co-activation model is not presented here. Therefore,

Table 4. Parameters for a virtual spring-damper.
k = 10.0

Virtual spring
Virtual damperζ = 1.0
ζ = 0.0(ke = 5.0)
(ke = 25.0)

a more systematic paradigm for determining ke according to the dynamic damping ellipsoid or a totally different modeling of muscle dynamics that may naturally prevent involvement of internal forces is needed in our nextwork.
6. Extension to ‘virtual spring-damper hypothesis’We extend our control stragegy to the ‘virtual springdamper hypothesis’ proposed by Arimoto et al. (2006).
This hypothesis made the point that in reaching movements, if we assume the existence of a virtual spring and avirtual damper between the end point and the desired point,
the end point of a human arm can converge to the desired
point smoothly according to this virtual spring-damper (seeFigure 19). Also in physiology, it is known that the cocontraction of agonist and antagonist muscles is gradually

Figure 19. Virtual spring-damper between the end point and the desired position.

Applied Bionics and Biomechanics

69

Output force of muscle f (N)30y-position (m)

0.6

0.4

0.2

f
1

f

25

2

f
3

20

f
4

f
5

15

f
6

f

10

7

f
8

f

5

9

0
0

0
–0.4

–0.2
0
x-position (m)12Time (s)34

0.2
Figure 22. Transient responses of the muscle’s output forces of
middle-range reaching with a virtual spring-damper.

√+ᾱ = − W̄ J {kx − ζ k ẋ} + (I 9 − W̄ W̄ )ke . (43)+60
Output force of muscle f (N)reduced as the process of learning a movement progresses(Osu et al. 1997). Now, we hypothesise that the internal
forces generated by the co-contraction of agonist and antagonist muscles are reduced, which enables an adequate
damping shape in the task space to be gradually obtained
as the learning process for reaching movements progresses.
We modify the control input to the virtual spring-damper
model. It is given by

f
1

f

50

2

f
3

f

40

4

f
5

30

f
6

f

20

7

f
8

f

10

9

T

0
0

1

2
Time (s)3

4Figure 23. Transient responses of the muscle’s output forces of
middle-range reaching with only a virtual spring.

0.6

E

0.6

F

D

0.4

y-position (m)y-position (m)Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

Figure 20. Middle-range reaching movement with a virtual
spring-damper when k = 10.0, ke = 5.0 and ζ = 1.0.

0.2

0

0.5
0.4

–0.2
0
x-position (m)0.2

Figure 21. Middle-range reaching movement with only a virtual
spring when k = 10.0, ke = 25.0 and ζ = 0.0.

G

0.3
B
0.2

–0.4

C

H
A

–0.3 –0.2 –0.1

0

0.1

0.2

x-position (m)Figure 24. End-point trajectories of all range-reaching movements with a virtual spring-damper.

70

K. Tahara et al.

y-position (m)End-point velocity (m/s)

E

0.6

F

D
0.5

G

0.4 C
0.3
B

H

0.2

A

–0.3 –0.2 –0.1

0

0.1

A
B
C
D
E
F
G
H

0.4
0.3
0.2
0.1
0

0.2

0

x-position (m)

1

2

3

4

5

Time (s)Figure 25. End-point trajectories of all range-reaching movements with only a virtual spring.

Figure 27. End-point velocities of all range-reaching movements with only a virtual spring.

Also, the closed-loop dynamics can be expressed as
1Ḣ(θ ) + S(θ, θ̇ ) θ̇ + J T kx
H(θ)θ̈ +2√+ { J T Jζ k + W P( AC +C 0 )W T }θ̇ = u = 0, (44)where the first term of the right-hand side of Equation (43)represents the virtual spring-damper and the second term
represents the internal forces generated by redundant muscles.
Table 4 shows each parameter for a virtual springdamper k and ζ used in the numerical simulations, andthese results are shown below. Figure 20 shows the simulation results with a virtual spring-damper, and Figure 21
shows those results without spring-damper. We see fromthese figures that the trajectories of the end point in thesetwo cases converge to the desired point smoothly. The final orientation of the wrist joint in the case of using a virtual spring-damper is more natural than that without one.
In contrast, Figures 22 and 23 show the transient responses

End-point velocity (m/s)

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009



A
B
C
D
E
F
G
H

0.4
0.3
0.2

of the muscle’s output forces. We see from these figuresthat the output forces in the case of using a virtual springdamper are obviously reduced compared with those without a virtual spring-damper, even though both end points
converge to the desired point smoothly. Figures 24 and 25
show the simulation results for reaching movements of all
ranges. We see from these figures that when the desired position is in the extension direction from the initial position(especially F, G, H in Figure 25), there is excessive overshoot, which depends on the configuration of the muscles.
However, this excessive overshoot can be reduced by introducing a virtual damper (see Figure 24). Figures 26 and 27
show the end-point velocities of all range-reaching movements with and without a virtual spring-damper. We see
from these figures that the velocity profiles become moresmooth when a virtual spring-damper is used. In paticular,
F, G, H in Figure 26, the oscillation is decreased. Thesesimulation results led us to suppose that reaching movements made before the movements were learned need acertain level of internal force to generate the joint damping
effect because there is no adequate virtual damper at theend point. After learning, the internal forces can be reducedby obtaining pro-prioceptive information in the CNS, so an
adequate virtual damper can be generated at the end point.
However, in order to discuss this more, we need to compare
a human’s movements and our scheme.

7.

0.1
0
0

1

2

3

4

5

Time (s)Figure 26. End-point velocities of all range-reaching movements with a virtual spring-damper.

Conclusion

This paper treated reaching movements by using themusculo-skeletal redundant arm model composed of three
joints and nine muscles with physiological non-linear muscle properties to mimic a human arm. Asymptotic convergence of the closed-loop dynamic system was obtained
based on the concept of stability theory on a manifold. We
finally confirm through numerical simulation results thatthe end point of the musculo-skeletal redundant arm couldDownloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

Applied Bionics and Biomechanics
converge to the desired point using only task-space feedback control with internal force adjustment under the existence of both joint and muscle redundancies, without any
complex mathematical computation. In addition, we introduced the dynamic damping elipse to determine the internal force vector and evaluated it. We extended our control
scheme to the virtual spring-damper hypothesis proposed
by Arimoto et al. (2006). By introducing a virtual damper
in addition to the virtual spring at the end point, we canreduce the internal forces. Moreover, the heterogeneity of
convergence depending on the initial state and desired position is also reduced. These results lead us to the understanding that learning of human movements might be interpreted as acquisition of adequate coordination of a virtual spring and damper in the task coordinates. However,
we have not treated the problem of learning of movementsin this paper yet. Therefore, in order to claim that our control scheme is reasonable, we should gain a more physical insight into determination of the internal force vector ke . In addition, we have only modeled co-contractionof agonist and antagonist muscles from the viewpoint of
their configurations and contractile velocities, but did notpresent dynamics of muscle co-activation which might begiven as a dynamic system of more complex differential
equations. By introducing such a dynamic model of coactivation, the structure of the internal force space described by (I 9 − W + W )ke and the mechanisum of cocontraction would be disclosed in a more explicit form.
Fortunately, since modeling of skeletal muscles is a hot
topic in physiology (Houdijk et al. 2006; Perreault et al.
2003) and many new models will be proposed hereafter,
we would like to introduce such skeletal muscle models
into our control strategy for redundant joint and musculoskeletal systems. On the other hand, one of the advantagesof our proposed method is that any planning of end-point
trajectories is unnecessary in advance. Similarly, Todorov
and Jordan (2002) have also proposed an optimal feedback
control theory that does not need end-point trajectories.
Our control method presented here may be interpreted as
one of such optimal feedback control schemes by regarding an artificial potential energy produced by the virtual
spring as one of the possible cost functions. The equivalence of both methods will be treated in our future works.
We would also like to introduce a learning strategy for theinternal force to investigate human-like multi-body control
mechanisms based on our control strategy.

Acknowledgements
This work was partially supported by the Kyushu University Research Superstar Program (SSP), based on the budget of Kyushu
University allocated under the President’s initiative, and by the
Ministry of Education, Science, Sports and Culture througha Grant-in-Aid for Young Scientists (B), 18760205, 2006–
2007.

71

References
Arimoto S, Sekimoto M, Hashiguchi H, Ozawa R. 2005. Natural
resolution of ill-posedness of inverse kinematics for redundant robots: a challenge to Bernstein’s degrees-of-freedom
problem. Adv Robot. 19(4):401–434.
Arimoto S, Sekimoto M. 2006. Human-like movements of roboticarms with redundant DOFs: virtual spring-damper hypothesis to tackle the Bernstein problem. In: Proceedings of IEEE
International Conference Robotics and Automation 1860–
1866.
Bernstein N. 1935. The problem of the interrelation of coordination and localization. Arch Biol Sci. 38:15–59 (reprinted
as Bernstein NA. 1967. The coordination and regulation of
movements. London: Pergamon.
Bizzi E, Hogan N, Mussa-Ivaldi F, Giszter S. 1992. Does the nervous system use equilibrium-point control to guide single andmultiple joint movements? Behav Brain Sci. 15:603–613.
Bizzi E, Polit A, Morasso P. 1976. Mechanism underlying achievement of final head position. J Neurophysiol.
3(9):435–444.
Delp S, Loan J. 2000. A computational framework for simulating
and analyzing human and animal movement. IEEE Comput
Sci Eng. 2:46–55.
Fel’dman A. 1966. Functional tuning of the nervous system with
control of movement or maintenance of steady posture. iii.
mechanographic analysis of the execution by man of the simplest robot tasks. Biophysics. 11:766–775.
Fel’dman A. 1986. Once more on the equilibrium-point hypothesis (λ model) for motor control. J Mot Behav. 18:17–54.
Flash T. 1987. The control of hand equilibrium trajectories in
multi-joint arm movements. Biol Cybern. 57(4–5):257–274.
Flash T, Hogan N. 1985. The coordination of arm movements: anexperimentally confirmed mathematical model. J Neurosci.
5(7):1688–1703.
Gribble P, Mullin L, Cothros N, Mattar A. 2003. Role of cocontraction in arm movement accuracy. J Neurophysiol.
89:2396–2405.
Hill A. 1938. The heat of shortening and the dynamic constantsof muscles. Proc R Soc B. (975):136–195.
Hogan N. 1984a. An organizing principle for a class of voluntary
movements. J Neurosci. 4(11):2745–2754.
Hogan N. 1984b. Adaptive control of mechanical impedanceby coactivation of antagonist muscles. IEEE Trans Automat
Contr AC-29(8):681–690.
Hollerbach J, Suh K. 1987. Redundancy resolution of manipulators through torque optimisation. IEEE Trans Robot Autom.
RA-3(4):308–316.
Houdijk H, Bobbert M, De Haan A. 2006. Evaluation of a Hill
based muscle model for the energy cost and efficiency ofmuscular contraction. J Biomech. 39(3):536–543.
Katayama M, Kawato M. 1993. Virtual trajectory and stiffness
ellipse during multijoint arm movement predicted by neural
inverse model. Biol Cybern. 69(5–6):353–362.
Kawato M, Furukawa K, Suzuki R. 1987. A hierarchical neuralnetwork model for control and learning of voluntary movement. Biol Cybern. 57(3):169–185.
Kawato M, Maeda Y, Uno Y, Suzuki R. 1990. Trajectory formation of arm movements by cascade neural network model
based on minimum torque-change criterion. Biol Cybern.
62(4):275–288.
Killeen P. 1992. Mechanics of the animate. J Exp Anal Behav.
57(3):429–463.
Mashima H, Akazawa K, Kushima H, Fujii K. 1972. The forceload-velocity relation and the viscous-like force in the frog
skeletal muscle. Jap J Physiol. 22:103–120.

72

K. Tahara et al.

Downloaded By: [Tahara, Kenji] At: 05:27 7 May 2009

Mclntyre J, Bizzi E. 1993. Servo hypotheses for the biological control of movement. J Mot Behav. 25(3):193–
202.
Morraso P. (1981) Spatial control of arm movements. Exp Brain
Res. 42:223–227.
Nakamura Y. 1991. Advanced robotics: redundancy and optimisation. Reading (MA): Addison-Wesley.
Osu R, Gomi H, Domen K, Yoshioka T, Kawato M. 1997. Decrease in joint stiffness during motor learning revealed frommeasured EMG signals. Paper presented at: The Society for
Neuroscience 27th Annual Meeting. 23:202.
Perreault E, Heckman C, Sandercock T. 2003. Hill muscle model
errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J Biomech.
36(2):211–218.
Rasmussen J, Damsgaard M, Surma E, Christensen S, deZee M, Vondrak V. 2003. AnyBody—a software system for ergonomic optimisation. Paper presented at: 5thWorld Congress on Structural and Multidisciplinary Op-

timisation. May 19–23, 2003; Lido di Jesolo, Venice,
Italy.
Tahara K, Luo ZW, Arimoto S, Kino H. 2005. Sensory-motor
control mechanism for reaching movements of a redundantmusculo-skeletal arm. J Robot Sys. 22(11):639–651.
Tahara K, Luo ZW, Arimotos. 2006. On control mechanism of
human-like reaching movements with musculo-skeletal redundancy. In: Proceedings of IEEE/RSJ International Conference Intellerence Robots and Systems 1402–1409; Edmonton, Canada.
Takegaki M, Arimoto S. 1981. A new feedback method for dynamic control of manipulators. ASME J Dynam Syst Meas
Control. 102:119–125.
Tsuji T, Morasso P, Goto K, Ito K. 1994. Human hand
impedance characteristics during maintained posture. Biol
Cybern. 72(6):475–485.
Todorov E, Jordan M. 2002. Optimal feedback control as a theoryof motor coordination. Nature Neuroscience. 5(11):1226–
1235.