Kinematics of Spherical Robots Rolling over 3D Terrains

Although the kinematics and dynamics of spherical robots (SRs) on at horizontal and inclined 2D surfaces are thoroughly investigated, their rolling behavior on generic 3D terrains has remained unexplored. is paper derives the kinematics equations of the most common SR congurations rolling over 3D surfaces. First, the kinematics equations for a geometrical sphere rolling over a 3D surface are derived along with the characterization of the modeling method. Next, a brief review of current mechanical congurations of SRs is presented as well as a novel classication for SRs based on their kinematics. en, considering the mechanical constraints of each category, the kinematics equations for each group of SRs are derived. Afterward, a path-tracking method is utilized for a desired 3D trajectory. Finally, simulations are carried out to validate the developed models and the eectiveness of the proposed control scheme.


INTRODUCTION
Spherical robots are a class of mobile robots that are generally recognized by their ball-shaped shell and internal driving components that provide torques required for their rolling motion. From the ball-shaped exterior, SRs inherit multiple advantages over other types of mobile robots, in which skidding, tipping over, falling, or friction with the surface makes them vulnerable or inefficient [1,2]. Despite all the unique features, the complicated nonlinear behavior of SRs has remained as a hurdle to fully comprehend their dynamics, motion kinematics, and unveil their maneuverability capabilities.
Although the earliest efforts to analytically capture the kinematics and dynamics of rolling geometrical balls on mathematical surfaces back to more than a century ago e.g. in E. Routh and S.A. Chaplygin's works, this topic is still an open discussion that is being investigated in more recent papers such as [3][4][5]. Numerous researches have been conducted on the mathematical modeling of kinematics and dynamics of SRs with a variety of mechanical configurations [6,7]. Among them, a widespread assumption is that SRs roll over an ideally flat horizontal plane. In a number of available works, where SRs are studied on non-horizontal surfaces, simplifying assumptions are made. For instance, in [2,8], two SR designs are investigated that the kinematics model's accuracy and its controller's tracking efficiency in different models along with comparison between their behavior.
The remainder of this paper is organized as follows: the next section is devoted to the model description.
The third section elaborates the problem of kinematics modeling of a sphere rolling over a 3D surface. In section four, the objectives of the paper are addressed, including presenting a new classification for SRs by their kinematics properties and then deriving kinematics equations of SRs based on specifications of each selected SR category. In section five, a 3D path tracking kinematic control scheme is presented. The simulation results in MATLAB are represented in section six, followed by the conclusion and future work.

MODEL DESCRIPTION
Consider a geometrical sphere rolling over a 3D surface  as shown in Fig. 1. To derive the kinematics equations of the rolling motion, the first step is to define the required reference frames, as  , defined as world, translated, surface-side tangent, robot-side tangent, and local reference frames respectively. In each reference frame set, the first element shows its origin and the rest are the reference frame's axes. Accordingly, as it is illustrated in Fig. 1 (1) Through utilizing transformation matrices, a fourth element of 0 is added to velocity vectors as translation does not apply to them. A fourth element of 1 in position vectors guarantees both rotation and translation.
In Fig. 1, the position of the contact point of the sphere with the 3D surface  is represented as 0 P :  In   L  , L x  and L y  are also tangent to T  . L x  is along the sphere's longitudinal direction. L z  is normal to T  coinciding with T z  and normal axis of the sphere, and L y  is mutually perpendicular to L x  and L z  following the right-hand rule. One can note that L y  can be considered as the projection of the sphere's lateral axis on The next step is to define rotation angles of the sphere. According to Fig. 1 implying that for any x and y , the elevation of the surface  along the vertical axis of W Z  is Assumptions: The mathematical model is derived based on the following assumptions: 1) The sphere rolls over a 3D surface without slipping, skidding or falling, i.e. it does not lose traction with the surface.
( , ) z f x y  is a continuous function of class 1 C or higher [26].
3) The sphere is rigid, and the radius of curvature of the surface is never smaller than the radius of the sphere. Therefore, the sphere always remains in a single point contact with the surface.

KINEMATICS MODELING OF SPHERES ON 3D SURFACES
The next phase to derive the kinematics equations of the sphere rolling over  is to calculate the utilized reference frames' unit vectors analytically along with deriving transformation matrices between reference frames. In this process, first, we will determine the vector quantities in   can be written as follows: where, 3 3  І is the identity matrix. From (3), the normal vector of surface  at position 0 P is given by: where, Defining n s as the following: we can write: Seemingly, the direction of n  is perpendicular to the surface and always downward, i.e., ˆ1 n Z    . As robotsite normal, n , is in the opposite direction to the surface-side normal and based on assumption (1), there is no relative velocity between   0 T  and   T  , therefore, we can write: In fact, n and ˆT k can be calculated by rotating n   by 180 degrees about 0 T i , therefore: (10) Next, we calculate the unit vectors of   T  axes. Let us consider a sphere that rolls over  in a way that 0   , that means only  and  contribute in the rolling action. Due to this rolling, the direction of sphere's longitudinal and lateral axes changes accordingly to match the geometry of  . Based on the Euler's rotation theorem, any Cartesian coordinate system in 3D space with a common origin are related by a rotation about a unique axis called Euler axis denoted as ê . Intuitively, this rotation maps axes of   T  to   Tr  . We know that the rotation maps n to ˆT r k and it is assumed that 0   , therefore, as it is depicted in Fig. 2, ê should lie in the tangent plane and be perpendicular to both n and ˆT r k , leading to: Considering ˆ[ 0,0,1] T Tr k  andn from (10), e  is calculated as: Consequently: Moreover, the angle of rotation can be calculated as: where atan2 is four-quadrant inverse tangent, where atan2 is four-quadrant inverse tangent, Therefore, to avoid singularity, from (14) and (15),  is written in the following conditional statement form: Equation (17) can be written in matrix form to calculate TTr  using Rodrigues' rotation formula. Let E  denote cross-product matrix for ê , i.e., a matrix when multiplied from the left with a vector, gives the cross product. It can be shown that: Now, it can be shown that rotation matrix TTr  can be written in the following form: that, using (18), can be written as the following: where, for the sake of brevity, sin( )  and os( ) c  are abbreviated to S  and C  respectively. Considering Also, form (13) and (15) it can be concluded that: As (3 1) 0 TTr    , and by substituting (21) and (22) in (20), we can write: Alternatively, to calculate TTr  , rotation quaternions can be utilized which is explained in details in appendix A. Depending on the application TTr  can be calculated from (19), (23), or through quaternions. To calculate ˆT j , we can write: Next, as   L  is resulted from the rotation of   T  through  about T z  , then, and LT  can be written as follows: Finally, from (4) , (23), and (25) it can be concluded that: 11 where, Consequently, the position of the sphere can be calculated by integrating (28) Obviously, 1 g and 2 g functions vary depending on the configuration of the SR. Transforming L V   to   W  by plugging (29) into (28) and using terms presented in (27) we have: Resulting in: 11 1 12 2 x a g a g    The next section, first, presents a classification based on kinematics specifications of SRs, then, provides calculation steps of L V   for each class of SR.

KINEMATICS MODELING OF SPHERICAL ROBOTS ROLLING OVER 3D SURFACES
This section is dedicated to deriving the kinematics of SRs rolling over 3D surfaces. First, we present a brief classification in which the most popular configurations of SRs are divided into two main categories and a number of subcategories. Secondly, by applying constraints of each category (or type), the kinematics of SRs on 3D surfaces is derived.
A. Continuous Rolling Spherical Robots (CR-SR): CR-SRs are capable of performing rolling maneuvers continuously, i.e., when the robot rolls in any direction there is no limitation for its rotation angle. However, based on their mechanical configurations, these robots can roll in different directions based on which this category can be divided into three subcategories explained in the following:

1) Triple-Axis Rolling Spherical Robot (3R-SR).
This subcategory is the most general type of CR-SR. As illustrated in Fig. 3, in 3R-SRs the internal components are mechanically capable of providing torques required for the robot to rotate about all of its three main axes, namely angular velocities of   ,   , and  . Therefore, 3R-SRs are omni-directional. Several mechanical configurations are available in the literature that rely on this actuation method [28][29][30].  Fig. 4 illustrates differential rolling lengths for a 3R-SR in an infinitesimally small time dt along ĥ and l directions. Based on assumption (1) and using equal-arc-length rule for rolling without slipping [31], we have: where R denotes the sphere's radius. Integrating (32) over time, L V   can be calculated as the following: Now, from (28) , (31) and (33)   2) Dual Axis Rolling Spherical Robot (2R-SR).
In this subcategory of SRs, the robot provides torques and consequently the ability to roll about its two main axes which are parallel to T  . There are a number of SR designs that relate to this type [23,[32][33][34][35][36][37] such as hamster wheel SRs. Comparing to 3R-SRs, similar steps should be taken in order to derive kinematics equation of 2R-SRs. The only distinction is that, in 2R-SRs, the robot is not capable of rotating about its vertical axis, i.e., it is assumed that 0    , therefore we can write: 11 12 x a R a R        , 21 22 31 32

3) Rolling and Turning Spherical Robot (RT-SR).
In addition to rolling action, RT-SRs are able to turn about their vertical axis in order to change their moving direction. In other words, internal parts can provide angular velocities of   and  , while it is assumed that 0    . RT-SRs are relatively less common but still some designs can be found in the literature [38][39][40]

B. Rolling and Steering Spherical Robots (RS-SR):
As it is shown in Fig. 5, consider as angular velocity of a RS-SR around its transverse axis that provides forward and backward rolling motion of the robot. In addition, the robot can tilt about its longitudinal axis as its steering mechanism providing turning ability for the robot through an arc shape rolling motion. Unlike CR-SRs, in RS-SRs tilting angle  is limited. Another difference is that in RS-SRs the transverse axis is defined in a way that tilts with the robot, so it is not always parallel to the tangent plane. This type of SRs uses an underactuated mechanical propulsion mechanism that allows the robot to roll in a nonholonomic manner [41][42][43][44]. Fig. 5 illustrates a RS-SR with tilting angle of  with respect to the tangent plane. The turning action of the sphere can be modeled as a cone with apex angle of 2 , purely rolling over the surface  , for which the instantaneous axis of rotation is the imaginary line passing through cone's vertex V , and 0 P that lies in the tangent plane T  . Considering the instantaneous differential rolling arc of the cone's base circle, we have: where, c r is the cone's base circle radius, d and d are differential rotation of the sphere during differential of time, dt , about the transverse and longitudinal axis respectively.  is the instantaneous radius of curvature of the turning motion that is the distance between the vertex of the cone and the contact point. These two quantities can be calculated as: and, Substituting (41) and (42) into (40), results in: Dividing both sides of (43) by dt , we have: Moving on, the position vector of the center point of the sphere in   L  in accordance with 0 P can be written as: and angular velocity of the sphere in   L  can be written as: Then, the linear velocity of the sphere's center point can be calculated as follows: Forthwith, we are able to form L V   as the following: Eventually, from (31) and (48) we can calculate components of W V   as follows: 11 12 cos( ) 21 22 cos( ) Equations in (49), using terms presented in (27), can be utilized to calculate the position of the SR, 0 P , based on input parameters of robot,  and  as the following: In all types of SR, to calculate the position of the center point of the sphere denoted as [ , , ] T Additionally, defining deviation angle  as the angle between the heading direction, ĥ , and e , we have: The objective of this algorithm is to converge the tracking error e along with  to zero for all four proposed categories of SRs. To that end, in each category the path planning algorithm computes values for actuated kinematics states, namely   ,   ,  , and  .

A. 3R-SR Kinematics Control
As 3R-SRs are capable of providing   ,   , and  , the controller is designed to provide required values accordingly as follows: is the angular velocity that is required for the sphere to pursue the target based on its speed. Therefore, this term is used as the bias value for   . For   , if 0   , i.e., the robot is not heading towards the target point, the robot tries to compensate its lateral distance to the target by rolling laterally using   in parallel to other actuators. If the robot is heading towards the target, 0   , the robot relies mostly on   to approach the target.
Finally,  is designed to be directly proportional to  as a measure of required turning to compensate the deviation angle.

B. 2R-SR Kinematics Control
Considering the fact that 2R-SRs are not able to turn, i.e. 0    , it can be implied that a 2R-SR is not capable of following the target by changing its heading direction and consequently the robot can only use   and   to move towards the target. Therefore, the controller is designed such that:   (56)

C. RT-SR Kinematics Control
RT-SRs can move towards the target point by compensating their deviation error,  through turning action, and then using   to approach and follow the target. This method is used to design the controller to provide required angular velocities as the following: (57)

D. RS-SR Kinematics Control
Calculated values for  and   are as the following for RS-SRs: It should be noted that  is driven and cannot be used as actuated value in RS-SRs.

SIMULATION RESULT
This section provides the simulation results through MATLAB/Simulink to verify the presented method and the controllers that are proposed in section 5. The simulation is carried out, where for the terrain surface,  is defined as the following:

CONCLUSION
This paper addressed the problem of modeling the kinematics and path tracking control of popular SR configurations rolling over 3D terrains. It is noteworthy that the aim of this paper is not to compare the performance of different types of spherical robots. In fact, it provides a classification of SRs based on their kinematics behavior to derive kinematics equations on 3D surfaces. Utilizing the proposed kinematics, one can use the Euler-Lagrange method to derive the dynamics of robots on uneven terrains by keeping in mind that the gravity direction is not always perpendicular to the tangent plane.

APPENDIX A
In this appendix, calculation steps of utilizing rotation quaternions instead of Rodriguez rotation method to calculate TTr  are presented. Generally, if we write unit rotation axis vector ê and ˆT TTr i j k r i k j k i r i k j r j k i r i j q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q , it can be shown that, plugging these terms in (A6) results in (20).