Igor Goncharenko1
3D Incorporated,
Kanagawa-ku, Sakaecho 1-1, Urban Square,
Yokohama 221-0052, Japane-mail: igor@ddd.co.jp

Mikhail Svinin
Mem. ASME
e-mail: svinin@bmc.riken.jp

Shigeyuki Hosoee-mail: hosoe@bmc.riken.jp
Bio-Mimetic Control Research Center,
RIKEN,
Anagahora, Shimoshidami, Moriyama-ku,
Nagoya 463-0003, Japan1

Dynamic Model, Haptic Solution,and Human-Inspired Motion
Planning for Rolling-Based
Manipulation
A virtual reality haptic system for capturing skillful human movements in control of ahemisphere rolling on a plane without slipping is presented in this paper. A dynamic
model of this nonholonomic rolling system with configuration-dependent inertia andgravity is derived, and a solver, required for the real-time haptic interaction, is implemented. The performance of the haptic system is verified under experiments with humansubjects. Experimental data recorded by the haptic system are analyzed and some common features of human movements in the precession phase of the manipulation of therolling system are observed. Finally, a simple actuation scheme, capturing these features,is proposed and verified under simulation. 关DOI: 10.1115/1.3074282兴

Introduction

In recent years, there appears an interest in robotic systems
where nonholonomic rolling constraints are used not only for manipulation but also for locomotion. In such systems self-propelled
movements are usually generated by creating imbalance andchanging the system inertia. The applications of rolling-based locomotion systems vary from the planetary surface exploration 关1兴
to the toy industry 关2兴. Different designs of spherical mobile robotic systems are reported in Refs. 关3–6兴. These vehicles can be
useful when the usage of traditional machines is limited or undesirable.
It should be noted that the analysis of motion of a sphere rollingon a plane without slipping is a classical problem in rational mechanics. The motion of a dynamically symmetric ball with a gyroscope inside 共and this was perhaps the first spherical robot!兲 wasstudied in Ref. 关7兴. The theoretical solution obtained there wasverified under experiment where the ball rolled on a plane powdered with likopodium. Note that from the viewpoint of analyticalmechanics the main problem is the integrability of the motion
equations in quadratures, and no notion of control is involved in
the analysis. Whether or not the wealth of results obtained inanalytical mechanics for rolling systems 关8–10兴 can be utilized bycontrol theorists and engineers yet remains to be seen.
One of the key problems in the control of nonholonomic systems with rolling constraints is the construction of motion planning algorithms. Typically, this problem is considered in the kinematic formulation 关11–15兴, and in practical implementations thecontrol system operates in the quasistatic mode 关16–18兴. The extension of the proposed motion planning algorithms to the dynamic domain is a very challenging, if not formidable, research
problem. The main difficulty here is that in addition to the nonholonomic constraints one also has the configuration-dependent
inertia and the gravity drift. This combination is rarely consideredin the control literature, and the motion planning frameworks forsuch systems are not yet developed.2
1

Corresponding author.
To the best of our knowledge, the only exception is the study 关20兴 where a
motion planner, employing state-of-the-art algorithms of computer science 共rapidly
exploring random trees兲, is proposed.
Contributed by the Engineering Simulation and Visualization Committee of
ASME for publication in the JOURNAL OF COMPUTING AND INFORMATION SCIENCE IN
ENGINEERING. Manuscript received September 30, 2007; final manuscript received
June 18, 2008; published online February 20, 2009. Guest Editors: J. Oliver, M.
Omalley, and K. Kesavadas.
2

It is also worth noting that the controllability analysis of nonholonomic systems with rolling constraints in dynamic formulation is highly nontrivial at best 关19兴 and, in general, quite cumbersome 关20兴. In this situation, when the controllability analysis is
not conclusive 共or simply not feasible兲 and the controllability
property needs to be hypothesized 关20兴, the insight gained in theanalysis of human movements in rolling-based manipulation andlocomotion may prove helpful for attacking the motion planningproblem.
It is interesting to note that trained humans and animals exhibit
highly skillful movements when performing rolling-based locomotion tasks 共see Fig. 1, reproducing the famous painting of Pablo
Picasso “The Girl on a Ball”兲. Watching them performing locomotion feats on a rolling sphere in circuses, one can wonder about
how these skills can be formalized in the form of a control strategy. To imitate the artistic movements, one could place a robotic
mechanism on top of a spherical object. This system, however,would be extremely complex for an initial study, and it would bereasonable to resort to simplifications. In a simplified formulation,
one could deal with a hemisphere and replace the robot with its
center of masses projected on the main hemisphere plane. The
driving principle for such a system would then be based on controlling the position of the center of mass and exploiting the nonholonomic rolling constraint to propel the hemisphere.
To extract the patterns of human movements in the control ofthis simplified system, one can resort to the virtual reality technology 关21兴. The model of a virtual hemisphere with the driving
mass point is suitable for the implementation in haptic systems,
because the mass point can be directly controlled by a humanusing a standard point-force haptic device. In a nutshell, this virtual reality haptic system would resemble the real manipulations
shown in the right part of Fig. 1. Changing the position of themass point by the haptic stylus would result into the mass redistribution and lead to the motion of the hemisphere, which can berealized in a haptic simulator. The development of such a simulator constitutes the main goal of this paper. Note that apart fromdiscovering nontrivial motion strategies in the control of nonholonomic systems, the simulator itself can be used also for studying
human motor behavior and for the development of games basedon the physical modeling.
Several issues need to be addressed in the implementation of
the haptic simulator. For the high fidelity force feedback, thesimulation of the virtual rolling system needs to be synchronizedwith the haptic hardware frequency. This brings the requirement

Journal of Computing and Information Science in Engineering
Copyright © 2009 by ASMEMARCH 2009, Vol. 9 / 011004-1

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

z[x o , yo , zo ]Σoy[x p , yp ]xc (u o , vo )
Σ co

Σb

Σ ca

c (u a , va )ψΣaFig. 1 Motivating example „left… and the prototype of the virtual reality haptic simulator „right…

Fig. 2 System formalizationof high computational efficiency. For this purpose, the state-space
dynamics of the rolling system needs to be established in theminimal form, as the often used representation employing
Lagrange multipliers 关6兴 can be more costly for real-time computations. Next, the force feedback generation procedures, suitable
for the control of this system, need to be assessed, and the feasibility of the haptic rendering needs to be tested under experimentswith human subjects.
The paper is organized as follows. First, in Sec. 2 we establish
the kinematic and dynamic models of the rolling system underconsideration. The architecture of a real-time haptic simulator isdescribed in Sec. 3. Here, we present a solver, integrating the state
dynamics, and describe force generation procedures for the haptic
interaction. Initial experiments with human subjects, using thehaptic simulator developed, are reported in Sec. 4. In Sec. 5, we
discuss how the experimental patterns of human movements canbe captured using simple mathematical models. Finally, conclusions are drawn in Sec. 6.

obtained by rotating the vector 关0 , 0 , −R兴T first around the x-axisof ⌺o through the angle vo and then around the y-axis 共of theframe associated with the first rotation兲 through the angle uo. In
this parametrization the lower hemisphere is selected by imposing
the following constraints:− ␲/2 ⬍ uo ⬍ ␲/2,

Mathematical Model

This section describes the kinematics and dynamics of a hemisphere with a mass point moving on the equatorial plane of thehemisphere. It is assumed that the hemisphere can roll on thecontact plane without slipping. The effects of rolling and spinning
friction are ignored.
2.1 Kinematic Model. To describe the system under consideration, we introduce the following coordinate frames 共see Fig. 2兲:⌺b is an inertial frame fixed at the base, ⌺o is a frame fixed at thecenter of the equatorial plane of the hemisphere, and ⌺a is a frame
fixed at the contact plane. In addition, at the contact point we
introduce the contact frame of the object ⌺co, and the contact
frame of the plane ⌺ca. In what follows, we will assume that theframe ⌺a coincides with ⌺b and the frame ⌺ca is parallel to ⌺a.
The position of a point on the sphere is parametrized as

冤 冥 冤
0

− sin uo cos vo

sin voc共uo, vo兲  Ry共uo兲Rx共vo兲 0 = R−R− cos uo cos vo

冥共1兲

where R is the radius of the sphere. The elementary rotation matrices around the axes x, y, and z are denoted by the corresponding
indices. In the parametrization 共1兲 the origin is placed at the southpole of the sphere. The geometrical meaning of the angular coordinates uo and vo is established as follows. The vector c共uo , vo兲 is011004-2 / Vol. 9, MARCH 2009

共2兲The parametrization 共共1兲兲 defines the Gaussian frame ⌺co whoseorientation relative to ⌺o iso

Rco共uo, vo兲 =冋

cu cv cu ⫻ cv
兩cu兩 兩cv兩 兩cu ⫻ cv兩册

共3兲

where cu  ⳵c共uo , vo兲 / ⳵uo and cv  ⳵c共uo , vo兲 / ⳵vo. Having defined
cu and cv, one can write this matrix in the explicit formo2

− ␲/2 ⬍ vo ⬍ ␲/2

Rco共uo, vo兲 =冤− cos uo sin uo sin vo − sin uo cos vo
cos vo0
sin uosin vocos uo sin vo − cos uo cos vo

冥共4兲Note that this matrix can be composed of the elementary rotations
as oRco共uo , vo兲 = Ry共uo兲Rx共vo兲Ry共␲兲.
Define the contact angle ␺ as the angle between the x-axes of⌺co and ⌺ca 共the z-axes of these frames are aligned as depicted inFig. 2兲. The orientation matrix of the sphere 共the orientation of ⌺orelative to the base ⌺b兲 can be defined as follows. First, by superposing the axes of ⌺co with those of ⌺ca, we rotate ⌺co first by
Ry共␲兲 and then by Rz共␺兲. Thus, bRcoRy共␲兲Rz共␺兲 = bRca and theorientation of ⌺ca relative to ⌺co are defined asco

Rca = Ry共␲兲Rz共␺兲 =冤− cos ␺ sin ␺

0

sin ␺

cos ␺

0

0

0

−1

冥

共5兲

Then, from closing the kinematic chains associated with the coordinate frames, one obtains bRo oRcoRy共␲兲Rz共␺兲 = bRca  I, where I
stands
for
the
unit
3⫻3
matrix.
Thereforeb
RoRy共uo兲Rx共vo兲Rz共␺兲 = I and the orientation of the sphere aredefined as
R  bRo = RzT共␺兲RTx 共vo兲RTy 共uo兲

共6兲˜ = 关␻
˜ x,␻
˜ y,␻
˜ z兴 be, respectively, theLet ṽ = 关ṽx , ṽy , ṽz兴 and ␻translational and rotational velocities of ⌺co relative to ⌺ca. Define
T

TTransactions of the ASME

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

roo

is denoted by the vector r p, and its absolute velocity by the vector
v p. The position of the mass point with respect to the origin of ⌺ois denoted by ␳ p.
It is assumed that the particle moves on the main hemisphereplane, and the vector of the driving force F acting on the point
mass lies in the main hemisphere plane. Therefore, the relative
position ␳ p and the driving force F can be parametrized as␳ p  iox p + joy p = RNp, F  io f x + jo f y = RNf, where p  兵x p , y p其T,f  兵f x , f y其T, R is the matrix of orientation of the frame ⌺o withrespect to ⌺b, and the projection matrix

p

s

b

rc

rp

c

冤 冥
1 0

N 0 1
0 0

Fig. 3 Definition of basic vectors

uco = 关uo , vo兴 and uca = 关ua , va兴 , where ua and va are the coordinates of the contact point in the contact plane. The evolution of
the contact coordinates is represented by the Montana equations
关22,23兴TT−1u̇co = M−1
o 共Ko + K̃a兲冉冋 册 冋 册冊
冉冋 册 冋 册冊
˜y
−␻
˜x
␻

−1
u̇ca = M−1
a R␺共Ko + K̃a兲

˜y
−␻
˜x
␻

− K̃a

ṽx

ṽy

+ Ko

ṽx

ṽy

共7兲
共8兲

˜ z + ToMou̇co + TaMau̇ca␺˙ = ␻

共9兲ṽz = 0

共10兲

where Ko, Mo, To and Ka, Ma, Ta stand for the curvature form,
connection form, and the metric tensor of the surfaces o and a,
respectively, at the contact point. R␺ is the upper 2 ⫻ 2 diagonal
block of the orientation matrix 共共5兲兲. K̃a = R␺KaR␺ is the curvature of a at the point of contact relative to the x- and y-axes of⌺co. In our specific case, Ma = I, Ka = O, Ta = 关0 , 0兴, and Mo= R diag兵cos vo , 1其, Ko = I / R, To = 关tan vo / R , 0兴. Here, I and O are,respectively, the unit and zero 2 ⫻ 2 matrices. Note that Eq. 共10兲
corresponds to the holonomic constraint maintaining the contact
between the bodies, while the rolling conditions are defined by thenonholonomic constraints ṽx = 0 and ṽy = 0.
Let ␻o = 兵␻x , ␻y , ␻z其T be the angular velocity of the frame ⌺o,defined in projections onto the axis of the base frame ⌺b. The
˜ are related as ␻˜ = coRca␻o. Expressing the convectors ␻ and ␻
tact kinematic equations 共7兲–共9兲 in terms of the components of
␻o, one obtainsu̇a = R␻y共11兲v̇a = − R␻x共12兲u̇o = − 共␻x sin ␺ + ␻y cos ␺兲/cos vo

共13兲v̇o = − ␻x cos ␺ + ␻y sin ␺

共14兲␺˙ = − ␻z − 共␻x sin ␺ + ␻y cos ␺兲tan vo

共15兲

Note that we do not impose the so-called no-spinning constraint
共␻z = 0兲 because it is difficult to justify its validity in our actuation
scheme.
2.2 Dynamic Model. Let ro be the vector of the origin of theframe ⌺o, vo its absolute velocity, and ␻o is the angular velocityof the frame ⌺o. The unit base vectors of the frame ⌺o are denotedby io , jo , ko, and those of the inertial frame ⌺b are denoted by ib,
jb, and kb. The system under consideration is propelled by moving
the mass point P and thus changing the system inertia and thegravity center 共see Fig. 3兲. The absolute position of the mass point

共16兲This matrix is necessary to keep the z-components of the projections of the vectors ␳ p and F onto the axes of the frame ⌺ozero.
The derivation of the dynamic model will be based on the useof D’Alembert’s principle 关24兴. Let ␦ro and ␦r p be the virtual
displacements of the origin of ⌺o and the driving point, respectively, and ␦␲o be the vector of infinitesimal rotation of ⌺o. For
the unconstrained system composed of the hemisphere and thepoint mass, D’Alembert’s principle can be written down as follows:˙ o ⫻ ms␳s + ␻o ⫻ 共␻o ⫻ ms␳s兲 − msg其␦ro · 兵msv̇o + ␻˙ o + ␻o ⫻ Jo␻o − ␳s ⫻ msg其+ ␦␲o · 兵ms␳s ⫻ v̇o + Jo␻+ ␦r p · 兵m pv̇ p − m pg − F其 = 0

共17兲

where ms is the mass of the hemisphere, m p is the mass of theparticle, ␳s is the vector of the center of mass of the hemisphere
defined in the frame ⌺o, Jo is the inertia tensor of the hemispherewith respect to the origin of ⌺o, and g is the gravity vector. Under3
the assumption that the hemisphere is homogeneous, ␳s = − 8 R ko,
2
2
Jo = 5 msR I, where I stands for the unit tensor.
First, we need to transform ␦r p in Eq. 共17兲 to ␦ro, ␦␲o, and␦⬘␳ p, where ␦⬘␳ p  io␦x p + jo␦y p = RN␦p is the virtual displacement of the driving point relative to ⌺o. Taking into account thatr p = ro + ␳ p, ␦r p = ␦ro + ␦␲o ⫻ ␳ p + ␦⬘␳ p, and ␻o ⫻ Jo␻o = 0, we
transform Eq. 共17兲 to the following form:˙ o ⫻ ms␳s + ␻o ⫻ 共␻o ⫻ ms␳s兲 + m pv̇ p
␦ro · 兵msv̇o + ␻˙o− 共ms + m p兲g − F其 + ␦␲o · 兵ms␳s ⫻ v̇o + Jo␻+ ␳ p ⫻ m pv̇ p − 共ms␳s + m p␳ p兲 ⫻ g − ␳ p ⫻ F其
+ ␦⬘␳ p · 兵m pv̇ p − m pg − F其 = 0

共18兲

Here, the terms involving v p and v̇ p can be expanded further by˙ o ⫻ ␳p
the kinematic formulas v p = vo + ␻o ⫻ ␳ p + ␳˚ p and v̇ p = v̇o + ␻ⴰⴰ+ ␻o ⫻ 共␻o ⫻ ␳ p兲 + 2␻o ⫻ ␳˚ p + ␳ p, where the local derivatives areⴰⴰdefined as ␳˚ p  ioẋ p + joẏ p = RNṗ and ␳ p  ioẍ p + joÿ p = RNp̈.
Next, we transform the virtual displacement ␦ro in Eq. 共18兲 tothe point of contact, defined by the vector rc 共see Fig. 3兲. The
reason is to introduce variations in quasicoordinates that would
vanish by virtue of the nonholonomic constraints 关25兴. Taking into
account that rc = ro + ␳c, ␦rc = ␦ro + ␦␲o ⫻ ␳c, and ␳c ⫻ g = 0, where␳c = −R kb, we transform Eq. 共18兲 to the following form:˙ o ⫻ ms␳s + ␻o ⫻ 共␻o ⫻ ms␳s兲 + m pv̇ p
␦rc · 兵msv̇o + ␻− 共ms + m p兲g − F其 + ␦␲o · 兵ms共␳s − ␳c兲 ⫻ v̇o
˙ o − ␳c ⫻ 关␻˙ o ⫻ ms␳s + ␻o ⫻ 共␻o ⫻ ms␳s兲兴
+ J o␻
+ 共␳ p − ␳c兲 ⫻ m pv̇ p − 共ms␳s + m p␳ p兲 ⫻ g− 共␳ p − ␳c兲 ⫻ F其 + ␦⬘␳ p · 兵m pv̇ p − m pg − F其 = 0.

共19兲

The rolling constraints, formulated in terms of the variations in

Journal of Computing and Information Science in Engineering

MARCH 2009, Vol. 9 / 011004-3

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

INTERCHANGABLE ITEMSCOMMON FRAMEWORK
Haptic Rendering

History Unit

Configuration
Constraint Library
VRML Repository

Physical Simulator

Graphics

Haptic Interface

Input

OutputPhantom Manipulator

3D Window

Fig. 4 System architecture

the quasicoordinates, can be written down as ␦rc = 0, and the corresponding terms in Eq. 共19兲 do not enter the dynamic equations
共the standard procedure based on the introduction of the Lagrangian multipliers moves the terms in the corresponding brackets tothe constraint reactions 关26兴兲.
In deriving the motion equations from Eq. 共19兲 we can eliminate the terms with v̇o using the fact that vo + ␻o ⫻ ␳c = 0, and,˙ o ⫻ ␳c. Now, setting ␦x p , ␦y p and ␦␲o as the intherefore, v̇o = −␻dependent variations, we obtain from Eq. 共19兲 the dynamic equations for the mass point
NTRT兵m pv̇ p − m pg其 = f

共20兲

and for the hemisphere˙ o − ␳c ⫻ 关␻o ⫻ 共␻o ⫻ ms␳s兲兴 + 共␳ p − ␳c兲 ⫻ m pv̇ p − 共ms␳s
J c␻
+ m p␳ p兲 ⫻ g − 共 ␳ p − ␳ c兲 ⫻ F = 0

共21兲

where Jc = Jo + ms共␳ˆ s␳ˆ c + ␳ˆ c␳ˆ s − ␳ˆ c␳ˆ c兲 is the inertia tensor of thehemisphere at the contact point, and hat is the skew-symmetric
operator such that a ⫻ a ⬅ âa.

3

Design and Implementation of the Haptic Simulator

To record and analyze human movements in the manipulationof the hemisphere rolling on the plane, we built a virtual reality
haptic system. The system architecture is schematically shown in
Fig. 4. The hardware part of the system is based on a point-force
haptic device 共PHANToM Premium 1.5 manipulator, 3DOF,
maximum exertable force 8.5 N兲 connected to a PC 共dual coreCPU, Intel Pentium 4, 3.0 GHz兲 through a PCI 共Peripheral Component Interconnect兲 interface.
The software runs in two threads responsible for independent
haptic and graphical renderings. In the haptic rendering loop, thePhysical Simulator module integrates the system of first-order ordinary differential equations 共ODEs兲, describing the evolution of
the rolling system under a specified force input in the time domain. The state dynamics are simulated in real time 共using fourthorder Runge–Kutta method 关27兴 with constant step h = 0.001 s兲
and the haptic feedback is supplied to the human hand through the
PHANToM stylus. Thus, the simulation is synchronized with theconstant haptic hardware frequency 共1000 Hz兲.
Similar to the architecture presented in Ref. 关28兴, the ODE
solver and the haptic feedback generation procedures are dynamically linked to the library of external constraints. When the application starts, the Physical Simulator module acquires in real time
the position and velocity of the PHANToM stylus, calculates thehaptic feedback force and the force applied to the virtual reality
object, and integrates the ODEs to obtain the state of the system atthe next moment of time. The geometric output of this module, theposition and orientation of the hemisphere with the mass point, isrendered in the 3D Window by the graphical thread.
The History Unit module records all the time histories of thestate variables as well as the haptic feedback forces. These data
are recorded at 100 Hz, which is enough to analyze human movements in our manipulation task. The physical parameters of the011004-4 / Vol. 9, MARCH 2009

Fig. 5Graphical user interface

rolling system as well as the initial and the target values of thestate variables are defined in the Configuration module. The
graphical representation of the rolling system 共referred to as a
Virtual Reality Modeling Language 共VRML兲 scene兲 is also defined in this module. Different segments 共files兲 of the VRML
scene are collected in the configuration repository. The VRML
files define the geometry of the hemisphere and the mass point
drawn as a small sphere on the main hemisphere plane.
The initial position of the rolling system is shown in Fig. 4. The
graphical user interface, shown in Fig. 5, corresponds to the stateof the rolling system after 30 s of haptic interaction. In this figure,the real-time scrolling window displays the evolution of thex-coordinate of the contact trajectory.
The rest of this section describes in detail the structure of the
ODE solver and the haptic feedback force models in our simulator.
3.1 Solver. The state of the system under consideration 共rolling hemisphere with moving mass point兲 can be described interms of 12 variables. The 12-dimensional state vector
x = 兵pT,ṗT, ␻To ,uTca,uTco, ␺其T

共22兲

or x = 兵x p , y p , ẋ p , ẏ p , ␻x , ␻y , ␻z , ua , va , uo , vo , ␺其 in the expanded
form. The solver integrates the system of differential equationsẋ = X共x , f兲, describing the evolution of the state vector x under the
input force f. To form the right-hand part X共x , f兲 of the state
˙ o as functions of x anddynamics, we first need to define p̈ and ␻
f. Upon expressing v̇ p in Eqs. 共20兲 and 共21兲 through the statevariables, one can represent the dynamic equations 共20兲 and 共21兲in the following matrix form:T冋

H pp H ps
Hsp Hss册冋 册 冋 册
p̈
hp
=
˙o
hs
␻

共23兲

where
H pp = m pI

共24兲
H ps = − m pNTRT共␳
p − ␳ c兲共25兲T

Hsp = m p共␳
p − ␳c兲RN = H ps

共26兲

2

Hss = Jc − m p共␳
p − ␳ c兲共27兲h p = f − m pNTRT关␻o ⫻ 共␻o ⫻ ␳ p兲 + 2␻o ⫻ ␳˚ p − g兴

共28兲Transactions of the ASME

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

hs = 共␳ p − ␳c兲 ⫻ RNf − m p共␳ p − ␳c兲 ⫻ 关␻o ⫻ 共␻o ⫻ ␳ p兲
+ 2␻o ⫻ ␳˚ p兴 + ␳c ⫻ 关␻o ⫻ 共␻o ⫻ ms␳s兲兴 + 共ms␳s + m p␳ p兲 ⫻ g
共29兲
Having formed the system matrices, one can define the accelerations

冋 册冋
p̈

=

˙o
␻

H pp H ps
Hsp Hss册冋 册 冋 册
−1

hp
hs



P共x,f兲

Q共x,f兲

共30兲

and then compile the kinematic and dynamic equations 共11兲–共15兲
and 共30兲 into the state-space form

冤 冥冤
ṗ

ẋ 

ṗ

p̈

P共x,f兲

˙o␻

Q共x,f兲

u̇a

v̇a

R␻ y

=

− R␻x

u̇o

− 共␻x sin ␺ + ␻y cos ␺兲/cos vov̇o− ␻x cos ␺ + ␻y sin ␺␺˙− ␻z − 共␻x sin ␺ + ␻y cos ␺兲tan vo

冥 X共x,f兲Fh = − kr共r p − rh兲 − kv共v p − vh兲

共31兲

To be able to simulate the state dynamics 共共31兲兲 in real time, thenumerical integration of the state equations must be synchronizedwith the constant haptic hardware frequency 共1000 Hz兲. This constitutes higher requirements to the efficiency of the computationof the right-hand part of the state equations 共共31兲兲. To increase thecomputational efficiency, one can reduce the inversion of the 5⫻ 5 system matrix in Eq. 共30兲 to the inversion of a 3 ⫻ 3 matrix.
Let us define冋

H pp H ps
Hsp Hss册 冋−1

H pp H ps
Hsp Hss册

共32兲

where, in the general case, the blocks if the inverse matrix aredefined as 关29兴−1
Hss = 共Hss − HspH−1
pp H ps兲

共33兲

H ps = − H−1
pp H psHss

共34兲

Hsp = − HssHspH−1
pp

共35兲

−1
−1
H pp = H−1
pp + H pp H psHssHspH pp

共36兲

The straight computation of the matrix blocks for the system under consideration leads to the following formulas:
T  T −1
Hss = 共Jc + m p共␳
p − ␳c兲koko 共␳ p − ␳c兲 兲

H pp =共37兲

T
H ps = − N RT共␳p − ␳c兲 Hss

共38兲

T

Hsp = − Hss共␳
p − ␳c兲RN = H ps

共39兲

1
T


I + N TR T共 ␳
p − ␳c兲 Hss共␳ p − ␳c兲RN
mp

共40兲

T

tual spring-damper models 关30兴. These models generate the haptic
force Fh based on the position and velocity differences between
the haptic proxy 共the end-point of the PHANToM manipulatortransformed to the virtual space兲 and a virtual object 共the driven
mass point, in our case兲. The haptic force is generated in real timeand is supplied to the hardware interface.
Note that in programming the PHANToM manipulator the haptic force Fh is specified in the inertial 共base兲 frame 关31兴. Let rh andvh be, respectively, the absolute position and velocity of the haptic
point. The absolute position and velocity of the particle are defined as r p = ro + ␳ p, and v p = ␻o ⫻ 共␳ p − ␳c兲 + ␳˚ p , where coordinates of the center of the sphere are ro = 兵ua , va , R其T, and ␳ p= RNp, ␳˚ p = RNṗ. Having defined these quantities, one can calculate the haptic force in the inertial frame. The following threesolutions are implemented in our haptic simulator.
First method (M1). The haptic interaction between the human
hand and the mass point is modeled by the elastic force with thespring and damping matrices proportional to the unit matrix

Thus, in the calculation of the matrix 共共32兲兲 we need to inverse
only the 3 ⫻ 3 matrix in Eq. 共37兲. Note that this inverse matrix isnever singular 共unless the mass of the sphere ms = 0兲. Also notethat in the software coding the last term 共共40兲兲 can be computed as−1
H pp = 共1 / m p兲I + H psHss
Hsp.
3.2 Force Input. The input force f, driving the point mass inthe dynamic model, is to be defined from the force Fh developed
at the human hand through the haptic interface. Typically, in the
impedance-type haptic devices the force Fh is defined using vir-

共41兲

where kr and kv are the scalar stiffness and damping coefficients.
To define the driving force f, one needs to transform Fh into themoving frame ⌺o, and then project it onto the main hemisphereplain. As a result, one obtains
f = N TR TF h

共42兲This method has one disadvantage. When the haptic force Fh is
orthogonal to the main hemisphere plane it does not affect the
motion equation. This may produce some unrealistic feeling for
the user.
Second method (M2). This method is a simple modification of
the first one. As before, we define the haptic force by Eq. 共41兲.
Then, we replace f in Eq. 共28兲 with NTRTFh, and in Eq. 共29兲 wereplace the term RNf with Fh. From the perceptional point ofview this model can be more realistic. In this model, when thehaptic force Fh is orthogonal to the main hemisphere plane it does
not affect the motion of the particle but it does affect the motionof the hemisphere.
Third method (M3). This method is similar to the second one
but the haptic force here is defined differently
Fh = − Kr共r p − rh兲 − Kv共v p − vh兲

共43兲

where Kr and Kv are the stiffness and damping matrices expressed in the inertial 共base兲 frame. How do we specify them?
Consider first the stiffness matrix. When there is no damping, Eq.
共b兲
共b兲
共43兲 can be written as Fh = −Kr ⌬r共b兲. Here, by the upper superscript we denote the frame 共the base ⌺b or the object ⌺o兲 in which
共b兲
共o兲
a vector or tensor quantity is expressed. Since Fh = RFh and
共o兲
共o兲
共b兲
共o兲
T 共b兲
共o兲
共o兲
⌬r = R⌬r , we get Fh = −R Kr R⌬r  −Kr ⌬r , andthereforeKr  Kr共b兲 = RKr共o兲RT

共44兲Here, the stiffness matrix in the object frame can be specified as
共o兲
Kr = diag兵krx , kry , krz其, where krx = kry and krz Ⰷ krx. If the coefficients are specified in this way, the user will feel the “plane constraint.” It will be harder to move the particle perpendicular to themain hemisphere plane compared with the “in-plane” movements.
Thus, this method of modeling the haptic interaction can be more
realistic.
The damping matrix in Eq. 共43兲 is defined in the similar way.
Specifically, we assign
Kv  Kv共b兲 = RKv共o兲RT

共45兲
共o兲

where the damping matrix in the object frame is Kv= diag兵kvx , kvy , kvz其, and kvx = kvy and kvz Ⰷ kvx.
It is instructive to clarify the differences among the three methods from the standpoint of perception. In the first method, when
the haptic force Fh does not lie in the main hemisphere plane, its

Journal of Computing and Information Science in Engineering

MARCH 2009, Vol. 9 / 011004-5

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

out-of-plane component does not enter the motion equations.
When it happens, the user may feel strange sensations, resulting
from mismatching of the visual and haptic information and leading to unrealistic perception of the system behavior. This situation
is excluded in the second and third methods, as we allow for theout-of-plane component of Fh to influence the motion of the rolling system. Thus, these methods are expected to produce more
realistic perception of the system behavior.
The difference between the second and third methods is only inthe specification of the stiffness and damping tensors. Specifically,in the second method these tensors are isotropic 共given by Eq.
共41兲兲 while in the third method they can be anisotropic 共defined by
Eq. 共43兲兲. Setting the coefficients of the anisotropic tensors in aright way can allow the user to feel the main hemisphere plane as
the out-of-plane movement of the haptic proxy would be felt withgreater resistance. The idea of the anisotropic tensors is basicallysimilar to that of the compliance control of robotic mechanisms inconstrained movements 关32兴. In principle, the third method can
produce more realistic perception. However, the design of thestable haptic interaction in this method becomes complex and requires special techniques 关33,34兴.

4

Experiments

The experimental system, described in Sec. 3, was used foranalysis of human movements in the manipulation of the virtual
rolling system. In the experiments, the parameters of the systemwere configured as follows. The gravity acceleration and themasses of the driving point and the hemisphere were set to9.81 m / s2, 1 kg, and 0.1 kg, respectively. The radius of the hemisphere was set as R = 4 cm.
Four naïve subjects 共males, aged between 25 and 30 years old,three right-handed and one left-handed兲 participated in the experiments. The subjects were instructed to transport the hemispherefrom the initial configuration, where all the state variables 共Eq.
共22兲兲 are zero, to a final state and stop it there. At the initialconfiguration the system was at rest, with the main symmetry axisof the hemisphere vertical and the point mass being on the top of
the hemisphere, as shown in the 3D Window block of Fig. 4.
Before starting the movement the subject positioned the haptic
proxy at the point mass and switched on the haptic interaction bypressing a button on the computer keyboard.
In the target configuration the vector uca was set as 兵−4 ,−4其T cm while the vectors p, ṗ, ␻, and uco were set zero. The
component ␺ of the state vector x at the target configuration was
not specified.3 Geometrically, this means that at the target configuration the system is at rest, the main symmetry axis of thehemisphere is vertical, the balancing point mass is on the top of
the hemisphere, but the orientation of the hemisphere, defined at
this state only by the angle ␺, can be arbitrary.
The subjects were requested to produce the reaching movements in a natural way, on their own pace, trading off the speed
and the comfortability. A movement trial was considered to be
successful if the movement task was completed within 4 min andwithin certain position and velocity tolerances. The subject wasgiven an audio feedback, generated by the computer, if a trial was
successful. In the course of preliminary evaluations, the following
tolerances were specified: ⌬␻ = ⫾ 0.01 rad/ s for the componentsof the angular velocity of the hemisphere, ⌬uca = ⫾ 0.4 cm for thecomponents of the contact point on the plane, ⌬uco= ⫾ ␲ / 18 rad for the components of the contact point on thehemisphere, ⌬p = ⫾ 0.4 cm and ⌬ṗ = ⫾ 0.4 cm/ s for the components of, respectively, the position and velocity of the driving
mass point.
3
In the preliminary evaluations, we found that specifying the complete orientation
at the target configuration made the reaching task much more complex. One canspeculate that the acquisition of the necessary motor skills for such a task requiresconsiderably more learning time.

011004-6 / Vol. 9, MARCH 2009

4.1 Assessment of the Force Feedback Methods. The experiments with the human subjects were conducted in two days.
On the first day we conducted a general evaluation of the subject
performance. The subjects familiarized with the experimental
setup, comprehended the reaching task, and performed movement
trials. During these preliminary experiments we also conducted
the assessment of the force feedback methods described in Sec.
3.2.
For the methods M1 and M2, the PHANToM stiffness anddamping coefficients were set as, respectively, kr = 250 N / m andkv = 10 N s / m. It was found that when using the method M1, all
the subjects could manipulate the rolling system and transport it to
the target configuration. However, they mentioned strange sensations and inconveniences in the manipulation of the hemispherewhen the hand movement was in the direction orthogonal to theequatorial plane. Compare with the method M1, the subjects
found the method M2 easier and more comfortable.
In assessing method M3, we set krx = kry = 250 N / m, kvx = kvy= 10 N s / m, and varied the ratio ␣ = krz / krx = kvz / kvx from 1 to 4.
When the ratio factor ␣ was closed to or exceeded 4, the subjects
evaluated the control as bumpy, uncomfortable, and unstable.
While the exact stability analysis is out of the scope of this paper,
the subjects’ perceptions can be partially explained by the apparatus limitations—the PHANToM manufacturer 关7兴 does not recommend setting up the stiffness coefficients more than 800 Nm.
In the course of the experiments, it was found that when theratio factor ␣ varied from 1.5 to 2.5, the control of the hemisphere
was stable, with smooth force sensations. As a result, the subjects
manipulated the rolling system faster; two subjects could complete the reaching task within 40 s time interval several times.
Overall, among the three methods for the generation of the feedback force, the subjects preferred method M3 共with ␣ = 2兲 as morenatural, comfortable, and realistic in perception. They sensed better the object reaction force and the inertial properties of thesystem.
Finally, it should be mentioned that in the assessment of theforce feedback methods we also experimented with a special casewhen the driving force f, sent to the physical simulator, was calculated correctly but the haptic force Fh was zeroed when applied
to the hand through the PHANToM stylus. This corresponds to thevisual-based control without the haptic perception. In this case all
the subjects failed to produce the requested reaching movements.
Even though the subjects could see that the point mass followed
the hand movement, they could not control the system as thehemisphere motion became quite chaotic. They reasoned that thesystem became uncontrollable without the haptic feedback. This
fact may evidence that the central nervous system cannot control
constrained movements based only on visual information.
4.2 Experimental Results. On the second day of the experiments, the subjects performed 15 movement trials, using method
M3 with ␣ = 2, and the data were collected for analysis. In all thetrials the movement task was completed successfully. However,
the movement time was varied from 46 s to 170 s.
It is interesting to note that the subjects showed roughly similar
motion patterns resulting possibly from a similar control strategy.
When asked, after the experiments, about their perceptions, thesubjects reported they felt loosing control of the system if theyvisually concentrated on the position of the point mass and/or thatof the haptic proxy. It appears that they paid visual attention
mostly to the current and target contact points on the contact plane
while controlling the motion of the hemisphere through the haptic
feedback.
When asked about their movements, the subjects could not
identify the specific trajectory of the hand. However, the general
movement strategy, reconstructed from their comments, can besummarized as follows. First, they tried to slightly unbalance thehemisphere and then gently increase its precession. Having developed a certain level of the hemisphere precession, they tried toTransactions of the ASME

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

Contact point on the plane

2

3

1

2

0

1

−1

yp (cm)va (cm)Driving point on the main hemisphere plane

−2

0

−1

−3
−2

−4
−5
−6

−5

−4

−3

−2

−1

u (cm)

0

−3
−3

1

−2

−1

1

2

3

p

a

Contact point on the plane

10x (cm)

2.5

Driving point on the main hemisphere plane

2

0

1.5

−1

1

−2

y (cm)0pav (cm)

0.5

−3
−4

−0.5

−5

−1

−6

−1.5

−7

−2

−8
−5

−4

−3

−2

−1

u (cm)

0

1

2

−2.5
−3

3

keep it constant until the rolling system reached a vicinity of thetarget point at the contact plane. At the final stage of the maneuver, the subjects decreased the precession and tried to stabilize theangular movements of the hemisphere. This was produced bysomewhat jerky, if not chaotic, movements of the hand. All thesubjects mentioned that the choice of the precession axis and theprecession speed is important for the successful steering of therolling system
There can be two reasons why all the subjects preferred thisstrategy. First, it appears to be conceptually simple and, therefore,
it could be more comfortable. The second reason is not physical or
physiological but rather social. The subjects could observe thetrials of the current performer and discuss them between eachother. So, once a successful motion strategy was picked by someone it could be reproduced by the others.
Two examples of the trajectory of the contact point in the contact plane are shown in Fig. 6. The trajectories are circularly
curved and appear to be helical. The corresponding trajectories ofthe driving point mass in the equatorial plane of the hemisphere
are shown in Fig. 7. To keep the graphs clean and readable, these
trajectories are taken for 20 s intervals of time in the middle partof the movement depicted in Fig. 6. The trace of the curves in Fig.
7 is almost closed, and the curves themselves appear to be roselike. One may conjecture that one leaf of the rose corresponds toone spiral loop in Fig. 6. This roselike pattern was visible in all−10

x (cm)

1

2

3

p

a

Fig. 6 Contact point trajectories for the right-handed „top… andleft-handed „bottom… subjects−2Fig. 7 Driving point trajectories for the right-handed „top… andleft-handed „bottom… subjects

the subjects’ data. The number of leaves varied from 8 or 15 共onaverage 10兲. The leaf radius changed from 0.4 cm to 1.0 cm 共onaverage 0.7 cm兲, while the outer radius of the roselike curve varied from 1.2 cm to 2.4 cm 共on average 1.8 cm兲. Also, the timespan featuring the roselike pattern varied from 18 s to 32 s 共on
average 23 s兲.

5Simple Actuation Scheme

Having observed the pattern of human movements, reported inSec. 4, we would like to explain its formation using, as much as
possible, simple mathematical tools. In general, the hand trajectory is composed of the acceleration, the steady state precession,and the deceleration/stabilization phases. In what follows, we will
be interested in capturing only the precession phase of this trajectory.
Assume for simplicity that the human can perfectly control andspecify the trajectory p共t兲 of the mass point on the main hemisphere plane. From first equation in Eq. 共23兲 we get
˙ o = f + ĥ p
H ppp̈ + H ps␻

共46兲

where ĥ p  共h p兲f=0 is defined by Eq. 共28兲 with f = 0. Let alsoĥs  共hs兲f=0. It is defined by Eq. 共29兲 with f = 0. Upon substituting
Eq. 共46兲 into Eq, 共23兲, after some algebra one obtains

Journal of Computing and Information Science in Engineering

MARCH 2009, Vol. 9 / 011004-7

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

Driving point on the main hemisphere plane

2.5

r2
r1

2

-k t

1.5

t

1

y (cm)

0.5

p

0
−0.5
−1

Fig. 8 Double pendulum as a propelling mechanism

−1.5
−2

共48兲

where Q̂ is defined by the right-hand side of Eq. 共47兲. By integrating the state dynamics 共共48兲兲, one first finds the system evolution x共t兲 and then estimates from Eq. 共46兲 the driving force f共t兲
corresponding to the prescribed motion of the particle p共t兲.
To qualitatively capture the pattern of human movements,
shown in Fig. 7, one may look for a family of roselike curves withan adjustable number of leaves. There are many such curves
known in the literature 关35兴. In our study, we model the trajectory
p共t兲 of the mass point as
x p共t兲 = r1 cos ␸共t兲 + r2 cos共k ␸共t兲兲

共49兲

y p共t兲 = r1 sin ␸共t兲 − r2 sin共k ␸共t兲兲

共50兲This curve describes the trajectory of the end-point of a virtual
double pendulum with the link lengths r1 and r2, rotating in theequatorial plane of the hemisphere, as shown in Fig. 8. The linksof the pendulum are rotated in the opposite directions, with thefirst link circulating k + 1 times slower than the second one. If theparameter k is integer the curve is closed and it has k + 1 leaves. To
capture the precession phase of the human trajectory, one specifies
the motion of the joint angle as ␸共t兲 = ⍀t, where ⍀ is a constantparameter.
Note that the curve 共Eqs. 共49兲 and 共50兲兲 can be viewed as apartial case of epitrochoids. However, the double pendulum interpretation is more attractive. One can think of the double pendulum
as a propeller, and it gives a mechanical analog for our actuation
scheme.
The idea of the double pendulum actuation scheme can be supported by a simulation example. In the example, in accordance

−2

−1

0

x (cm)

12

3pFig. 9 Roselike trajectory of the mass point on the equatorial
planewith the experimental data reported in Sec. 4, we set r1 = 1.5 cm
and r2 = 0.8 cm. To obtain one closed cycle with nine leaves we
set k = 8. The precession constant is set as ⍀ = −␲ / 10 rad/ s,
which corresponds to the clockwise rotation of the first link andcounterclockwise rotation of the second link. The position of thedriving mass point in the equatorial plane of the hemisphere, defined by the curve 共Eqs. 共49兲 and 共50兲兲, is shown in Fig. 9. As canbe seen, qualitatively, the roselike pattern of human movements is
captured.
Having specified p共t兲, we then integrated the state dynamic
equations 共共48兲兲. The motion is simulated on the time interval 关0,35兴 s. Note that the initial position of the mass point is not zero,and in the simulation we start to move the hemisphere from theequilibrium state corresponding to p共0兲. This can be considered asa modeling uncertainty. Figure 10 shows the trajectory of the contact point on the contact plane. The trajectory resembles thatshown in Fig. 6 共top兲. As can be seen, qualitatively, the basic
features of this trajectory 共the curvature of the trace and the spirallike loops along the trace兲 are similar to those observed in thehuman movements.
While the comparison conducted may look optimistic andspeaks in favor of the actuation scheme given by Eqs. 共49兲 and
共50兲, the real story of human manipulations is, most probably, still4

Contact point on the plane

3
2
1
0
a

Note that under our assumption 共the direct specification of p兲 theacceleration of the mass point p̈ does not influence the angularacceleration of the sphere as it does not enter the sphere dynamics
共共47兲兲.
The dynamics of the system under consideration can be described now in terms of eight state variables. The eightdimensional state vector x = 兵␻To , uTca , uTco , ␺其T, or, in the expanded
form x = 兵␻x , ␻y , ␻z , ua , va , uo , vo , ␺其T. Since the trajectory of themass point is specified explicitly as a function of time, the state
dynamics ẋ = X共x , p共t兲 , ṗ共t兲兲 are defined as

011004-8 / Vol. 9, MARCH 2009

−2.5
−3

共47兲

v (cm)−1˙ o = 共Hss − 共␳␻p − ␳c兲RNH ps兲 共ĥs − 共␳ p − ␳c兲RNĥ p兲ua(cm)

0

1

2

Fig. 10 Trajectory of the contact point in the contact plane

Transactions of the ASME

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

Contact angle
80
70
60

ψ(rad)

50
40
30
20
10
0

0

5

10

15

20

25

30

35

time (sec)Fig. 11 Precession angle ␺

far from being that simple. In this connection, it should be notedthat in the course of simulation we found the motion of the hemisphere to be very sensitive to the choice of k, ⍀, r1, and r2. We
also found that the curvature of the arc trace observed in theexperimental data is achieved only for a relatively narrow range ofthese parameters. In this range the precession speed becomes almost constant 共see Fig. 11兲 after recovering from initial imbalancecaused by the nonzero values of p共0兲. The exact tuning of theseparameters 共or their possible self-synchronization in the course of
the production of comfortable periodic movements兲 remains an
open issue.
Note that the proposed mechanism for capturing the pattern of
human movements is applicable only for the middle part of the
trajectory. It does not work for the beginning 共acceleration phase兲and end 共deceleration/stabilization兲 parts where the patterns of thehand movements are quite different. In principle, the actuation
scheme given by Eqs. 共49兲 and 共50兲 can form a part of a motion
planner for the generation of humanlike movements for rollingbased locomotion and manipulation systems. However, it captures
only the general tendency of the rhythmic movement in the precession phase. Different mechanisms for the generation of such
movements can be realized using the approaches developed inbiological cybernetics. In particular, it is believed that in the central nervous system the generation of complex periodic movements is organized in the form of a central pattern generator 关36兴
that, computationally, can be implemented by a network ofcoupled oscillators 关37兴. Testing this approach for the generationof humanlike movements in control of nonholonomic systems
with rolling constraints would be an interesting subject for future
research.

6

Conclusions

A virtual reality haptic system for capturing skillful human
movements in control of a hemisphere rolling on a plane withoutslipping has been presented in this paper. To steer the hemisphere,a mass point, moving in the main hemisphere plane and changing
the system inertia and the gravity center, was introduced. The
driving principle for the rolling system is based on controlling theposition of the mass point by the human hand via a haptic device.
A dynamic model of the nonholonomic rolling system withconfiguration-dependent inertia and gravity has been derived, anda solver, required for the real-time haptic interaction, has been
implemented.
The performance of the haptic system has been verified underexperiments with human subjects. It was confirmed that the control of the rolling system is very hard, if not impossible, without

the kinesthetic sensations provided by the feedback forces. Whenestimating the haptic feedback, the subjects preferred the one thatprovides the better feeling of the reaction from the main hemisphere plane. Experimental data recorded by the haptic systemwere analyzed and some common features of the human movement in the manipulation of the rolling system have been observed. In particular, it was found that the trajectory of the contact
point in the contact plane appears to be helical, with its path being
circularly curved. At the same time the trajectory of the driving
mass point looks like a periodic rose curve. A simple actuation
scheme, capturing this pattern of human movements, has been
proposed and verified under simulation.
As we see it now, the haptic simulator presented in this paper
can be used mainly for discovering nontrivial motion strategies in
the control of nonholonomic systems with rolling constraints, as
the analytical motion planning techniques for such systems arevery hard to construct. Apart from the main purpose, we believe
the simulator developed can be used also for studying humanmotor behavior and in gaming applications based on the physicalmodeling.
Several issues should be addressed in future research. First, it
would be reasonable to include in our model the effects of rollingand spinning frictions ignored in the current study. Next, a rigorous analysis of the stability of the haptic interaction would benecessary for the finer tuning of the feedback force models withanisotropic stiffness and damping. Finally, it would be interesting
to analyze the formation of the human trajectory not only in the
middle 共precession phase兲 part but also in the acceleration anddeceleration phases. The prediction of the whole trajectory is, perhaps, an ambitious goal, but is worth trying.

Acknowledgment
The authors would like thank the anonymous referees for helpful comments and valuable suggestions.

References
关1兴 Michaud, F., de Lafontaine, J., and Caron, S., 2001, “A Spherical Robot forPlanetary Surface Exploration,” Proceedings of the Sixth International Symposium on Artificial Intelligence, Robotics and Automation in Space, CanadianSpace Agency, St-Hubert Quebec, Canada.
关2兴 Michaud, F., and Caron, S., 2002, “Roball, the Rolling Robot,” Auton. Rob.,
12共2兲, pp. 211–222.
关3兴 Bicchi, A., Balluchi, A., Prattichizzo, D., and Gorelli, A., 1997, “Introducing
the ‘Sphericle’: An Experimental Testbed for Research and Teaching in Nonholonomy,” Proceedings of the IEEE International Conference on Robotics
and Automation, Vol. 3, pp. 2620–2625.
关4兴 Bhattacharya, S., and Agrawal, S., 2000, “Spherical Rolling Robot: A Designand Motion Planning Studies,” IEEE Trans. Rob. Autom., 16共6兲, pp. 835–839.
关5兴 Javadi A., A. H., and Mojabi, P., 2004, “Introducing Glory: A Novel Strategy
for an Omnidirectional Spherical Rolling Robot,” ASME J. Dyn. Syst., Meas.,
Control, 126共3兲, pp. 678–683.
关6兴 Balasubramanian, R., 2006, “Modeling and Control Techniques for a Class ofMobile Robot Error Recovery Problems,” Ph.D. thesis, Carnegie Mellon University, Pittsburgh, PA.
关7兴 Bobylev, D., 1892, “On a Ball With a Gyroscope Inside Rolling Without
Sliding on a Horizontal Plane,” Mat. Sb., 16共3兲, pp. 544–581.
关8兴 Chaplygin, S., 1949, “On a Ball Rolling on the Horizontal Plane,” Investigations on the Dynamics of Nonholonomic Systems 共Classics of Natural Science兲,N. Semenova, ed., GITTL, Moscow, pp. 72–99, in Russian.
关9兴 Kilin, A., 2001, “The Dynamics of Chaplygin Ball: The Qualitative and Computer Analysis,” Regular Chaotic Dyn., 6共3兲, pp. 291–306.
关10兴 Borisov, A., and Mamaev, I., 2002, “Rolling of a Rigid Body on Plane andSphere: Hierarchy of Dynamics,” Regular Chaotic Dyn., 7共2兲, pp. 177–200.
关11兴 Li, Z., and Canny, J., 1990, “Motion of Two Rigid Bodies With Rolling Constraint,” IEEE Trans. Rob. Autom., 6共1兲, pp. 62–72.
关12兴 Jurdjevic, V., 1993, “The Geometry of the Plate-Ball Problem,” Arch. Ration.
Mech. Anal., 124共4兲, pp. 305–328.
关13兴 Marigo, A., and Bicchi, A., 2000, “Rolling Bodies With Regular Surface:Controllability Theory and Applications,” IEEE Trans. Autom. Control, 45共9兲,
pp. 1586–1599.
关14兴 Mukherjee, R., Minor, M., and Pukrushpan, J., 2002, “Motion Planning for a
Spherical Mobile Robot: Revisiting the Classical Ball-Plate Problem,” ASME
J. Dyn. Syst., Meas., Control, 124共4兲, pp. 502–511.
关15兴 Svinin, M., and Hosoe, S., 2008, “Motion Planning Algorithms for a Rolling
Sphere With Limited Contact Area,” IEEE Transactions on Robotics, 24共3兲,pp. 612–625.
关16兴 Camicia, C., Conticelli, F., and Bicchi, A., 2000, “Nonholonomic Kinematics

Journal of Computing and Information Science in Engineering

MARCH 2009, Vol. 9 / 011004-9

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm

and Dynamics of the Spherical,” Proceedings of the IEEE/RSJ International
Conference on Intelligent Robots and Systems, IROS’2000, pp. 805–810.
关17兴 Harada, K., Kawashima, T., and Kaneko, M., 2002, “Rolling Based Manipulation Under Neighborhood Equilibrium,” Int. J. Robot. Res., 21共5–6兲, pp.
463–474.
关18兴 Suzuki, K., Svinin, M., and Hosoe, S., 2005, “Motion Planning for Rolling
Based Locomotion,” Journal of Robotics and Mechatronics, 17共5兲, pp. 537–545.
关19兴 Hristu-Varsakelis, D., 2001, “The Dynamics of a Forced Sphere-Plate Mechanical System,” IEEE Trans. Autom. Control, 46共5兲, pp. 678–686.
关20兴 Choudhury, P., and Lynch, K., 2002, “Rolling Manipulation With a Single
Control,” Int. J. Robot. Res., 21共5–6兲, pp. 475–487.
关21兴 Burdea, G., 2003, Virtual Reality Technology, 2nd ed., Wiley, New York.
关22兴 Montana, D., 1988, “The Kinematics of Contact and Grasp,” Int. J. Robot.
Res., 7共3兲, pp. 17–32.
关23兴 Murray, R., Li, Z., and Sastry, S., 1994, A Mathematical Introduction to Robotic Manipulation, CRC, Boca Raton, FL.
关24兴 Wittenburg, J., 1977, Dynamics of Systems of Rigid Bodies, Teubner, Stuttgart,
Germany.
关25兴 Neimark, J., and Fufaev, N., 1972, “Dynamics of Nonholonomic Systems,”
Translation of Mathematical Monographs, Vol. 33, American MathematicalSociety, Providence, RI.
关26兴 Lurie, A., 2002, “Analytical Mechanics,” Foundations of Engineering Mechanics. Springer-Verlag, Berlin.
关27兴 Press, W., Teukolsky, S., Vetterling, W., and Flannery, B., 1992, Numerical

011004-10 / Vol. 9, MARCH 2009

Recipes in C: The Art of Scientific Computing, Cambridge University Press,
New York.
关28兴 Goncharenko, I., Svinin, M., Kanou, Y., and Hosoe, S., 2006, “Predictability of
Rest-to-Rest Movements in Haptic Environments With 3D Constraints,” Journal of Robotics and Mechatronics, 18共4兲, pp. 458–466.
关29兴 Horn, R., and Johnson, C., 1986, Matrix Analysis, Cambridge University
Press, Cambridge.
关30兴 Burdea, G., 1996, Force and Touch Feedback for Virtual Reality, Wiley, NewYork.
关31兴 SensAble Devices Inc., PHANToM Master User’s Manual, http://
sensable.com.
关32兴 Craig, J., 1989, Introduction to Robotics: Mechanics and Control, 2nd ed.,Addison-Wesley, Reading, MA.
关33兴 Hogan, N., and Colgate, J., 1989, “Stability Problems in Contact Tasks,” Robotics Review, Vol. 1, J. Craig, O. Khatib, and T. Lozano-Perez, eds., MIT,
Cambridge, MA, pp. 339–348.
关34兴 Miller, B., Colgate, J., and Freeman, R., 2000, “Guaranteed Stability of Haptic
Systems With Nonlinear Virtual Environments,” IEEE Trans. Rob. Autom.,
16共6兲, pp. 712–719.
关35兴 Gray, A., 1998, Modern Differential Geometry of Curves and Surfaces With
Mathematica, 2nd ed., CRC, Boca Raton, FL.
关36兴 Marder, E., and Calabrese, R., 1996, “Principles of Rhythmic Motor Pattern
Generation,” Physiol. Rev., 76, pp. 688–717.
关37兴 Matsuoka, K., 1987, “Mechanisms of Frequency and Pattern Control in the
Neural Rhythm Generators,” Biol. Cybern., 56, pp. 345–353.

Transactions of the ASME

Downloaded 22 Feb 2009 to 133.6.1.49. Redistribution subject to ASME license or copyright; see http://www.asme.org/terms/Terms_Use.cfm