# Motion Planning for a Spin-Rolling Sphere on a Plane

Seyed Amir Tafrishi<EMAIL_ADDRESS>Mikhail Svinin  Motoji
Yamamoto  Yasuhisa Hirata  Department of Robotics, Tohoku University, Sendai,
Japan  Department of Information Science and Engineering, Ritsumeikan
University, Kyoto, Japan  Department of Mechanical Engineering, Kyushu
University, Kyushu, Japan

###### Abstract

The paper deals with motion planning for a spin-rolling sphere when the sphere
follows a straight path on a plane. Since the motion of the sphere is
constrained by the straight line, the control of the sphere’s spin motion is
essential to converge to a desired configuration of the sphere. In this paper,
we show a new geometric-based planning approach that is based on a full-state
description of this nonlinear system. First, the problem statement of the
motion planning is posed. Next, we develop a geometric controller implemented
as a virtual surface by using the Darboux frame kinematics. This virtual
surface generates arc-length-based inputs for controlling the trajectories of
the sphere. Then, an iterative algorithm is designed to tune these inputs for
the desired configurations. The feasibility of the proposed approach is
verified by simulations.

###### keywords:

Path planning , spin-rolling , Darboux frame , shortest path , smooth
trajectory , geometric control

††journal: Latex

##  1 Introduction

Manipulation of a spin-rolling sphere to a desired configuration is an
important research problem that can find many applications. Different
mechanisms can realize a spin-rolling motion of the sphere on a plane (see
Fig. 1). In the field of dexterous manipulation, rolling fingertips with
convex surfaces can be utilized to grasp and manipulate different objects [1,
2, 3]. For example, Yuan et al. developed rolling spherical fingertips for
grasping mechanisms [3]. There are propulsion systems that can realize spin-
rolling motions with multiple-degrees-of-freedom (DoFs) actuators. For
instance, spherical mobile robots with mass-imbalance propulsion [4, 5], with
cart-based actuators inside their spherical shell [6, 7] and internal motor
drivers [8, 9, 10] can have a similar feature with simple increase in number
of actuators. Additionally, the ballbot robots with drivers on top of the
rolling ball can be included in the kinematic model of spin-rolling motion
[11, 12]. Furthermore, nano/micro spherical particle manipulation on plane
surfaces is finding important applications [13, 14, 15]. Thus, the problem of
path planning for a spin-rolling sphere is applicable to different mechanisms.

Figure 1: Realization of the spin-rolling motion of the sphere, e.g., the
micro/nano particle manipulation [13], the spherical rolling robots [6],
ballbots [12] and the Dexterous manipulation [3].

The analysis of the spin-rolling motion can be done by using the conventional
ball-plate system kinematics [16, 17] but we cannot spin the ball directly by
the rotation of the upper sandwiched plane [18]. Different kinematic
parametrization of rolling contact have been considered in [18, 19, 20, 21,
22, 23]. To conduct motion planning, Kiss et al. [1] used three independent
planes to manipulate the sphere at the kinematic level and controlled the
relative angles of the sphere without considering desired plane configuration.
Date et al. used the advantage of spinning indirectly to control the ball-
plate system [24]. Their motion planning algorithm was based on shifting the
coordinate of the actuating plane in iterations with respect to different
reference frames. This change can be looked like a rotational virtual center.
However, the time scaling with the included coordinate transformation of the
kinematic model can bring the system to an uncontrollable state [25].

In the literature on the pure rolling, different feed-forward [16, 26, 27, 10]
and feedback control [24, 25, 28, 29] methods were developed. Since the change
of the spin angle affects the orientation of the sphere, the conventional
planning approaches, such as the ones based on the geometric phase shifting
[30, 31, 32, 33, 34], cannot be adjusted naturally for the spin-rolling motion
planning. Arthurs et al. [35] and Jurdjevic [16] proposed planning approaches
based on the optimal control theory. This research direction was developed
later on by Sachkov [36, 37]. As an alternative to feedforward planning, to
stabilize the non-differentially flat ball-plate system, an iterative steering
was formulated by Oriolo and Vendittelli[25]. However, the solutions had
serious fluctuations as trajectories converge to the desired states. Thus,
their created trajectories were not easy to be realized by dynamical systems
[38] and there were singularities in different regions of spherical manifold
due to the locality of the solutions. Beschastnyi [39] studied the optimality
problem of the spin-rolling sphere. Extremal trajectories were parameterized,
and their cut times were estimated for optimality. It was proposed that
Maxwell time can be determined while the sphere follows the straight path for
the optimal solution. However, this research did not cover the control problem
for arbitrary desired states.

In this paper, we approach the motion planning problem with two main
motivations: considering a spin-rolling sphere on the plane and developing
optimal smooth trajectories for the arbitrary desired states. Our approach is
based on parameterizing the spin-rolling sphere motion by using the Darboux-
frame-based kinematics of the ball-plate system originally proposed by Cui and
Dai [21, 2, 40]. In our previous work[18], we transformed the underactuated
ball-plate system to a fully-actuated one with the arc-length-based control
inputs. In this paper, after explaining the motion planning problem, we
introduce a virtual surface to manipulate the control inputs of the Darboux-
frame-based kinematic model to a desired configuration. To the best of our
knowledge, this is the first proposed geometric arc-length-based control
strategy for the motion planning of a spin-rolling sphere. This control
strategy separates the time scale from the kinematic equations, which allows
constructing motion with different convergence rates in a given time. Finally,
we propose a tuning algorithm to iteratively solve the differential kinematic
model until achieving a successful final full configuration. It is important
to note that the designed Darboux-frame parametrization increases the number
of planning parameters that simplify the developed algorithm.

This paper is organized as follows. In Section 2, the motion planning problem
is stated and a distance constraint due to the no-sliding constraint is
explained. Next, Section 3 includes the design of the geometric controller for
a Darboux frame kinematics. In Section 4, the iterative algorithm with
included tuning variables is described. Simulations of this new approach are
demonstrated and discussed in Section 5. Finally, we conclude our findings in
Section 6.

##  2 Problem Statement

It is assumed that the sphere that rolls and spins along a straight path on
the plane, and the the path on the sphere, connecting the initial and final
configurations, must be established. Because of the no-sliding constraint the
length of the curve on the sphere ${\bm{\mathbf{L}}}_{o}$ is equal to that on
the the plane ${\bm{\mathbf{L}}}_{s}$, the curve may not be long enough to
reach certain configurations on the sphere. This limitation is mainly due to
our considered approach with a restricted path on the plane. To deal with this
issue, we find a distance constraint that determines the minimum length of the
curve for reaching the sphere to its desired states. Note that this limitation
only appears when the length between the initial and final states of the plane
is less than the circumference of the sphere. We calculate this distance
constraint for the plane by applying the Gauss-Bonnet theorem [41] to the
created circular path on the desired angular states of the sphere.

Figure 2: Statement of planning problem while the sphere with contact path of
${\bm{\mathbf{L}}}_{o}$ follows a straight path ${\bm{\mathbf{L}}}_{s}$ on
plane.

Configuration of a sphere on a plane (see Fig. 2) is described with the sphere
position on plane ${\bm{\mathbf{P}}}$=$(u_{s},v_{s})$ and its local
orientation $\bm{\Psi}=(u_{o},v_{o})$ which $\psi$ is the spin angle between
the sphere and plane. We consider our initial and final states with following
notations $\\{{\bm{\mathbf{P}}}_{0}$, $\bm{\Psi}_{0}$,
$\psi_{0}\\}$$=\\{u_{s,0}$, $v_{s,0}$, $u_{o,0}$, $v_{o,0}$, $\psi_{0}\\}$ and
$\\{{\bm{\mathbf{P}}}_{f}$, $\bm{\Psi}_{f}$, $\psi_{f}\\}$$=\\{u_{s,f}$,
$v_{s,f}$, $u_{o,f}$, $v_{o,f}$, $\psi_{f}\\}$, respectively. The coordinate
system and preliminary information of this ball-plate system is described in
A. In this planning, the state equation is solved in $k$ iterations till
finding admissible paths, where the traveled paths on the sphere and plane are
noted as ${\bm{\mathbf{L}}}_{o}$ and ${\bm{\mathbf{L}}}_{s}$. We plan the
sphere to reach its final configuration $C_{0}\rightarrow C_{1}$ with a
straight trajectory on ${\bm{\mathbf{L}}}_{s}$ in the given $t_{f}$ time.

It is assumed that the initial states
$\\{{\bm{\mathbf{P}}}_{0},\bm{\Psi}_{0},\psi_{0}\\}$ and the final desired
states of the sphere $\\{\bm{\Psi}_{f},\psi_{f}\\}$ are given. However, the
final desired position on the plane ${\bm{\mathbf{P}}}_{f}$ has to be chosen
under the consideration that the length between initial and final position
$||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}=[(u_{s,f}-u_{s,0})^{2}+(v_{s,f}-v_{s,0})^{2}]^{\frac{1}{2}}$
is greater than a minimum distance variable $d$ where this constraint is
calculated by using the rest of desired and initial states. This length
limitation appears because the total curve length of ${\bm{\mathbf{L}}}_{o}$
on the sphere is always the same (no-sliding constraint) as
${\bm{\mathbf{L}}}_{s}$ that goes in a straight trajectory. Also, the sphere
has to arrive at desired local coordinate $\bm{\Psi}_{f}$ with different
approaching angles, as the desired spin $\psi_{f}$, while it follows a
straight optimal line rather than maneuvering freely through the plane
$U_{S}$. Hence, by knowing the sphere desired states
$\\{\bm{\Psi}_{f},\psi_{f}\\}$, plane’s desired states ${\bm{\mathbf{P}}}_{f}$
are chosen in any desired direction $\\{u_{s,f},v_{s,f}\\}$ with the length
$||{\bm{\mathbf{L}}}_{s}||=||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$
larger than the minimum distance variable $d$ as

$d<||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2},$ (1)

Note that this constraint is not important and the sphere can reach all
possible configurations, if the length of the desired position
$||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$ is set larger than
sphere circumference $2\pi R_{o}$.

In order to find the minimum distance variable $d$, a path is constructed as a
circular segment on the sphere that passes $\bm{\Psi}_{f}$, and then the cap
area is changed toward true desired $\psi_{f}$ by Gauss-Bonnet theorem [41].
Then, by the trigonometric relations in the new under-cap area, we find the
length of ${\bm{\mathbf{L}}}_{o}$ as the minimum distance, where
$||{\bm{\mathbf{L}}}_{s}||=||{\bm{\mathbf{L}}}_{o}||$.

Figure 3: Limit of minimum distance $d$ that is required to reach desired
$\psi_{f}$. Note that $G_{f}$ shows the desired goal angle on plane.

Let the minimum distance $d=2\pi a_{c}\alpha$ be the circumference of the
orange cap-area that sphere travels, red line ${\bm{\mathbf{L}}}_{o}$ on the
right sphere at Fig. 3, where $a_{c}$ and $\alpha$ are the base diameter of
this circumference and the angle of sector from $\bm{\Psi}_{0}$ to
$\bm{\Psi}_{f}$ on the sphere $U_{C}$, respectively. First, we find the base
diameter $a_{c}$ of the cap with the following formulation

$\displaystyle a_{c}=2\sqrt{(S_{t}/\pi)-h^{2}_{c}},\;\;\;h_{c}=S_{t}/(2\pi
R_{o}),$ (2)

where $S_{t}$ and $h_{c}$ are the total area of the cap-shaped region [orange
and blue regions in Fig. 3] and height of the cap area. Notice that the
rotation of the sphere along the red path ${\bm{\mathbf{L}}}_{o}$ doesn’t
cover the whole cap-shaped region. Thus, we find the total area as
$S_{t}=2S^{\prime}_{c}+\Delta S$ where $S^{\prime}_{c}$ is the constructed cap
area of the circular path by $\bm{\Psi}_{f}$ [see the blue cap-shaped region
in Fig. 3] and $\Delta S$ is the area change [see the orange region in Fig. 3]
for reaching from the spin angle $\psi^{\prime}$ with cap area of
$S^{\prime}_{c}$ to desired spin angle $\psi_{f}$. Here, $S_{t}$ has two parts
which combining them gives the full segment of the circular cap (dashed and
solid red line ${\bm{\mathbf{L}}}_{o}$). First, the area of the cap-shaped
region of $S^{\prime}_{c}$ is determined from the closed simple circle that
passes $\bm{\Psi}_{f}$

$\displaystyle
S^{\prime}_{c}=\left(\alpha^{\prime}/2\pi\right)\left[({a}^{\prime}_{c}/2)^{2}+{h^{\prime}}_{c}^{2}\right],$
(3)

where ${a^{\prime}}_{c}$, ${h^{\prime}}_{c}$ and $\alpha^{\prime}$ are the
diameter of cap’s base, height of the cap and the sector angle from sphere
initial to final configuration $\bm{\Psi}_{f}$ on the cap’s base,
$\alpha^{\prime}$ equals to $\pi$. Following parameters in (3) are calculated
with the help of Eq. (29) and trigonometric relations, shown in Fig. 3,

$\begin{split}&h=R_{o}\left[1-\cos u_{o,f}\cos
v_{o,f}\right],\;{a^{\prime}}_{c}=\left[h^{2}+R_{o}^{2}\left(\sin^{2}v_{o,f}+\sin^{2}u_{o,f}\cos^{2}v_{o,f}\right)\right]^{1/2},\;\\\
&{h^{\prime}}_{c}=\begin{cases}&R_{o}\left[1-\cos\left(Q_{o,f}/2\right)\right],\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{a}^{\prime}_{c}\leq
2R_{o}\\\
&R_{o}\left[1-{a}^{\prime}_{c}\cos\left(Q_{o,f}/2\right)/2R_{o}\right],\;\;\;\;\;\;{a}^{\prime}_{c}>2R_{o}\end{cases}\end{split}$

where $Q_{o,f}$ is the angle from initial contact point
$\bm{\Psi}_{0}=\\{0,0\\}$ to desired local coordinates $\bm{\Psi}_{f}$ as

$Q_{o,f}=\begin{cases}&\pi-2\cos^{-1}\left(h/a^{\prime}_{c}\right),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;h\leq
R_{o}\\\
&\frac{\pi}{2}-\sin^{-1}\left((h-R_{o})/R_{o}\right),\;\;\;\;\;\;\;\;\hfill
h>R_{o}\end{cases}$

Also, by knowing the $a_{c}$ and $a^{\prime}_{c}$, the new sector angle from
$\bm{\Psi}_{0}$ to $\\{\bm{\Psi}_{f},\psi_{f}\\}$ ( ${\bm{\mathbf{L}}}_{o}$ as
the red line in Fig. 3) becomes

$\alpha=\frac{1}{2\pi}\cdot\begin{cases}&1-2\sin^{-1}(a^{\prime}_{c}/a_{c}),\;\;\;\;\;\;\;\;\;\;\;\Delta
S\geq 0\\\
&2\sin^{-1}(a^{\prime}_{c}/a_{c}),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Delta
S<0\end{cases}$ (4)

To find the achieved spin angle $\psi^{\prime}$ by the contacted
$S^{\prime}_{c}$ cap-area, we use the Gauss-Bonnet theorem [41] while the
sphere touches the circular closed blue path

$\Delta\psi=\psi^{\prime}-\psi_{0}=\iint_{{S}^{\prime}_{c}}\kappa_{o}\;dS=S^{\prime}_{c}/R_{o}^{2},$
(5)

Figure 4: Normalized minimum distance $d/R_{o}$ for different desired angles
where $u_{o,f}\in[0,\pi]$ and $v_{o,f}=0.01$.

where $\kappa_{o}=1/R^{2}_{o}$ is the Gaussian curvature and the initial spin
angle is assumed $\psi_{0}=0$. Finding the spin angle $\psi^{\prime}$ also let
us calculate $\Delta S$ in $S_{t}$. The required area-change $\Delta S$ for
achieving $\psi_{f}$ is calculated by the same Gauss-Bonnet theorem as follows

$\Delta
S=(\psi_{f}-{\psi^{\prime}})/\kappa_{o}=R^{2}_{o}(\psi_{f}-{\psi^{\prime}}).$
(6)

After obtaining the total area $S_{t}$ for reaching $\psi_{f}$, $a_{c}$ in (2)
and $\alpha$ in (4) give us the value of minimum distance $d$. Fig. 4 depicts
the normalized minimum distance for different desired angles which shows how
the small distance $||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$ can
be for the example final configurations. It is clear that as the desired local
coordinate $\bm{\Psi}_{f}$ moves to upper-hemisphere, ${\bm{\mathbf{P}}}_{f}$
has to set for the longer length from ${\bm{\mathbf{P}}}_{0}$. Also. the same
property is true for rising the spin angle $\psi_{f}$. Moreover, the largest
distance requirement happens at lower points of $\bm{\Psi}_{f}$
($Q_{o,f}\leq\pi/2$) with larger desired spin angles $\psi_{f}$. Note that
this constraint (1), with the assumption of a simple circular curve, is looked
as the shortest length under the isoperimetric inequality, $4\pi S_{t}\leq
d^{2}$.

##  3 Geometric Controller for Darboux-Frame-Based Kinematics

We utilize the model of the Darboux frame at the contact point of the spin-
rolling sphere and plane presented in [18]. The kinematics readily suits both
contact trajectories and arbitrary parameters of the surfaces [2, 21, 42].
Also, this transformation provides two significant benefits through our
planning approach. First, the spin-rolling angular rotations explicitly appear
on the relative curvature and torsion [21] that makes it easier for
manipulation. Second, the Darboux frame separates the time variable from the
planning due to its time- and coordinate-invariance.

As show in [18], the transformed fully-actuated Darboux-frame-based kinematics
model of a spin-rolling sphere on a plane in the time domain is presented as

$\displaystyle\begin{split}&\left[\begin{array}[]{c}\dot{u}_{s}(t)\\\
\dot{v}_{s}(t)\\\ \dot{u}_{o}(t)\\\ \dot{v}_{o}(t)\\\
\dot{\psi}(t)\end{array}\right]=\delta\Bigg{(}\left[\begin{array}[]{c}\sin(\theta+\varphi)\\\
\sin(\theta+\varphi)\\\
\frac{\sin(\theta+\varphi)[\sin{\psi}-\cos{\psi}]}{R_{o}\cos{v_{o}}}\\\
\frac{\sin(\theta+\varphi)[\cos{\psi}+\sin{\psi}]}{R_{o}}\\\
\frac{\tan{{v}_{o}}[\sin(\theta+\varphi)(\sin\psi-\cos\psi)+\cos\varphi]}{R_{o}}\end{array}\right]+\left[\begin{array}[]{c}-R_{o}\sin(\theta+\varphi)\\\
-R_{o}\sin(\theta+\varphi)\\\
\frac{\sin(\theta+\varphi)[\cos{\psi}-\sin{\psi}]}{\cos{v_{o}}}\\\
-\sin(\theta+\varphi)[\sin{\psi}+\cos{\psi}]\\\
\tan{{v}_{o}}[\sin(\theta+\varphi)(\cos\psi-\sin\psi)]\end{array}\right]\;\gamma_{s}\\\
&+\left[\begin{array}[]{c}R_{o}\sin(\theta+\varphi)\\\
-R_{o}\cos(\theta+\varphi)\\\
\frac{-\sin{(\psi+\theta+\varphi)}}{\cos{v_{o}}}\\\
-\cos{(\psi+\theta+\varphi)}\\\
-\tan{{v}_{o}}\sin{(\psi+\theta+\varphi)}\end{array}\right]\;\beta_{s}+\left[\begin{array}[]{c}0\\\
0\\\ 0\\\ 0\\\ -1\end{array}\right]\;\alpha_{s}\Bigg{)}.\end{split}$ (7)

where five state system $\\{{\bm{\mathbf{P}}}(t),\bm{\Psi}(t),\psi(t)\\}$ has
$\\{\alpha_{s},\beta_{s},\gamma_{s},G_{f},\delta\\}$ inputs which

$\displaystyle\theta(\beta_{s},\gamma_{s},G_{f})=\cot^{-1}\Big{[}\frac{1}{\beta_{s}}\Big{(}\frac{1}{R_{o}}(1-\tan
G_{f})+\gamma_{s}(-1+\tan G_{f})-\beta_{s}\tan{G_{f}}\Big{)}\Big{]}-\psi_{q},$
$\displaystyle\varphi(G_{f})=\psi_{q}+\begin{cases}&\pi,\;\;\;\;\;\;\;\;\;\;\;-\frac{3\pi}{4}<G_{f}<0\;\;\;\;\&\;\;\;\;\;\;0\leq
G_{f}<\frac{\pi}{4}\\\
&0,\;\;\;\;\;\;\;\;\;\;\;\;-\pi<G_{f}<-\frac{3\pi}{4}\;\&\;-\pi\leq
G_{f}<\frac{\pi}{4}\end{cases}.$ (8)

where $\psi_{q}$ is the spin angle deviation. The three arc-length-based
control inputs $\\{\alpha_{s},\beta_{s},\gamma_{s}\\}$ correspond to the
angular velocities of the sphere in time domain. The remaining angular $G_{f}$
and rolling rate $\delta$ inputs are for directing sphere on fixed surface
(plane) while they are constraining virtual surface inputs
$\\{\alpha_{s},\beta_{s},\gamma_{s}\\}$. The controllability of this geometric
model (7) is in [18]. This geometric-based controller converges the rotating
object (sphere) trajectory to the desired angular states of the sphere while
$G_{f}$ keeps the sphere on suitable direction with rolling rate $\delta$. In
here, we define the goal angle by
$G_{f}$=$\tan^{-1}\big{[}(v_{s,f}-v_{s,0})$$/(u_{s,f}-u_{s,0})\big{]}$; hence,
$\theta$ and $\varphi$ will change relative to any values of arc-length
inputs, for keeping the sphere always along $G_{f}$ angle to reach
${\bm{\mathbf{P}}}_{f}$. Therefore, the sphere stays in the prescribed
straight direction by substituting constant value to $G_{f}$ from the desired
position on the plane ${\bm{\mathbf{P}}}_{f}$.

The Darboux-frame-based kinematics (7)-(8) with the inputs in the arc-length
domain requires a geometric control for converging spin-rolling sphere to
desired angular configuration $\bm{\Psi}_{f}$. We introduce a virtual surface
$U_{V}$ to define these arc-length-based inputs
$\\{\alpha_{s},\beta_{s},\gamma_{s}\\}$. The virtual surface is a surface
sandwiched between sphere and plane at the contact frame. From a physical
point of view, deformation (curvature changes) of this virtual surface is
projected onto both sphere and plane trajectories [18]. These changes in the
curvature manipulates the curve on the sphere ${\bm{\mathbf{L}}}_{o}$ like a
flexible rope. To bend this rope-like curve toward the desired full-
configuration, we propose a tuning algorithm to update the variables of this
virtual surface in the incoming section.

a)  b)

Figure 5: a) Geometric shape of the used virtual Surface $U_{V}$. Note: Green
and blue surfaces are for $R_{t}=0$ as the sphere and $R_{t}=R_{v}$ as the
helicoid, b) Design of the arc-length-based inputs by using desired virtual
surface. Note that schematic is drawn for $u_{o}=0$.

Here, we design the virtual surface for controlling the curve
${\bm{\mathbf{L}}}_{o}$ on the sphere $U_{C}$ when ${\bm{\mathbf{L}}}_{s}$ is
specified as a straight path. In order to develop a geometric controller based
on the virtual surface, we use the induced curvature (differential
characteristics), with the geodesic curvature $k_{g}$, the normal curvature
$k_{n}$ and the geodesic torsion $\tau_{g}$, that consists of the rolling
convex object (sphere) on a plane with sandwiched virtual surface as follows
[18]

$\displaystyle
k_{g}=k^{o}_{g}-k^{s}_{g}-\alpha_{s},\;\tau_{g}=\tau_{g}^{o}-\tau^{s}_{g}-\beta_{s},\;k_{n}=k^{o}_{n}-k^{s}_{n}-\gamma_{s},$
(9)

where $\\{k^{o}_{g},\tau^{o}_{g},k^{o}_{n}\\}$,
$\\{k^{s}_{g},\tau^{s}_{g},k^{s}_{n}\\}$ and
$\\{\alpha_{s},\beta_{s},\gamma_{s}\\}$ are the geodesic curvature, geodesic
torsion and normal curvature for, respectively, the rolling object (sphere),
plane and virtual surface. In the geometric controller, the curvature
properties are designed with similarity to rotating object (sphere) manifold
$U_{C}$ (30) for planning on the virtual surface. However, the moving object
is a sphere with $\tau^{o}_{g}=0$. This cause the model (7) to be unstable
[18] because $\beta_{s}=0$. As a solution, geodesic torsion of a helicoid
shape is applied where the geometric shape is shown in Fig. 5-a (see B for
details of derivation)

$\tau^{v}_{g}=\frac{1}{R_{v}^{2}}|R_{v}^{2}\cos^{2}{v_{v}(t)}-R^{2}_{t}|^{\frac{1}{2}},$
(10)

where $R_{v}$, $R_{t}$ and $v_{v}(t)$ are defined by main spherical radius,
sum of spherical and torsion radii and $v$-curve angle of helicoid surface,
respectively. As a keynote, geometric surface shown in Fig. 5-a is able to
transform from the sphere to the helicoid surface by changing $R_{t}$.

Next, we utilize (9) to construct the desired virtual surface (see Fig. 5-b)
with the normal and geodesic curvature of the spherical surface
$\\{k_{n}^{d}(s),k_{g}^{d}(s)\\}$ in (30)-(31) and the geodesic torsion of the
helicoid surface $\tau^{d}_{g}(s)$ (10) as follows

$\displaystyle
k_{g}^{d}(s)=\tan{\zeta}/R_{t}=\tan{v_{o,f}}/R_{o}-\alpha_{s},\;\tau^{d}_{g}(s)=\frac{1}{R^{2}_{o}}\left|R^{2}_{o}\cos^{2}{v^{\prime}_{o}}-R^{2}_{t}\right|^{\frac{1}{2}}=\beta_{s},\;k_{n}^{d}(s)=1/R_{n}=1/R_{o}-\gamma_{s}.$
(11)

where $R_{t}=R_{n}+R_{g}$ and $v^{\prime}_{o}=v_{o,f}-v_{o}(t)$ are the total
radius and the angle feed of $v$-curve in which $R_{n}$ and $R_{g}$ are the
desired normal curvature and geodesic torsion radii. Also, we assume $\zeta$
is the desired stereographic projection angle for the spherical object along
$k^{d}_{g}\bm{e}_{3}$ vector. By considering that the diameter of projected
curve is same (red and blue curves in Fig. 5-b ) and $R_{t}\tan\zeta=R_{o}\tan
v_{o}$, $\zeta$ is determined

$\zeta=\tan^{-1}\left[R_{o}\tan\left(v_{o,f}+\zeta^{\prime}\right)/R_{t}\right],$
(12)

where $\zeta^{\prime}$ is the constant angle shift that will be used during
planning updates. Also, $R_{n}$ and $R_{g}$ are defined as

$\displaystyle R_{n}(t)=R_{g}(t)=\left[R_{i}(t)+R_{a}\right]/2,$ (13)

Figure 6: $R_{i}$ determination by incircle of $u^{\prime}_{o}$.

where $R_{i}$ and $R_{a}$ are the built-in circle radius of $u$-curve and
tuning constants that happens by algorithm, respectively. By using schematic
of Fig. 6, built-in radius $R_{i}$ is computed with the change of the designed
incircle radius as follows

$R_{i}(t)=\begin{cases}\begin{split}&\left[\frac{(S_{i}-r_{i})^{2}(S_{i}-l_{i})}{S_{i}}\right]^{\frac{1}{2}},\;\;\;\;\;\;\;\;\;\;\;0\leq
u^{\prime}_{o}<\frac{\pi}{2}\\\
&\frac{R_{o}}{\mu_{r}}+\left[\frac{(S_{i}-r_{i})^{2}(S_{i}-l_{i})}{S_{i}}\right]^{\frac{1}{2}},\;\;\frac{\pi}{2}\leq
u^{\prime}_{o}\leq\pi\end{split}\end{cases}$ (14)

where $S_{i}=(2r_{i}+l_{i})/2$ is the area of encompassed triangle of the
incricle, $u^{\prime}_{o}=u_{o,f}-u_{o}(t)$ is the convergence of $u$-curve
angle, $\mu_{r}$ is the scaler to limit the maximum built-in radius $R_{i}$,
and also $r_{i}=R_{o}/\cos u^{\prime}_{o}$ and
$l_{i}=2R_{o}\tan{u^{\prime}_{o}}$ are adjacent and hypotenuse sides of
isosceles triangle (see Fig. 6).

######  Remark 1.

The sphere curvature in (11) is included to cancel out the existing sphere
$U_{C}$ properties that drift terms is presenting in the kinematic model (7).
Thus, this design gives a direct manipulation of virtual surface curvatures
(s-domain) on ${\bm{\mathbf{L}}}_{o}$ trajectory for $\bm{\Psi}_{f}$.

Now, we can find the geometric control inputs from (11) as

$\displaystyle\alpha_{s}=\tan
v_{o,f}/R_{o}-\tan\zeta/R_{t},\;\beta_{s}=\frac{1}{R^{2}_{o}}|R^{2}_{o}\cos^{2}{v^{\prime}_{o}}-R^{2}_{t}|^{\frac{1}{2}},\;\gamma_{s}=\frac{R_{n}-R_{o}}{R_{n}R_{o}}.$
(15)

By substituting the desired values $(u_{o,f},v_{o,f})$, the curves of
${\bm{\mathbf{L}}}_{o}$ are constrained on the virtual surface with the
projection on the rotating object, as the example of Fig. 5. This constrained
curve is created from a conservative vector field in (15) where independent
inputs act as force fields
${\bm{\mathbf{F}}}=\nabla{\bm{\mathbf{E}}}_{in}(u_{o}(t),v_{o}(t),u_{o,f},v_{o,f})=\\{\alpha_{s},\beta_{s},\gamma_{s}\\}$
to manipulate at contact frame on $U_{C}$ manifold [43]. By using kinematics
(7)-(8), it is interpreted that the desired virtual surface bends this vector
field towards the desired local coordinate $\bm{\Psi}_{f}$. Note that our
curve manipulation on the sphere manifold $U_{C}$ can be imagined as a
flexible rope with length of
$||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$ that $\\{R_{n},R_{g}\\}$
and $\zeta$ terms change the curvature radius and angle of $v$-curve for
${\bm{\mathbf{L}}}_{o}$ trajectory on the spherical surface $U_{C}$.

It should also be noted that the proposed virtual surface is designed in a way
that eases our planning problem with spherical curvature properties
$k^{v}_{n}$ and $k^{v}_{g}$. However, more complicated virtual surfaces with
existing geodesic torsion $\tau^{v}_{g}$ can be applied to this planning
approach for creating different contact paths but it can complicate planning
problem, e.g., algorithm tuning, which requires a separate study.

##  4 Motion Planning Algorithm

###  4.1 Iterative Tuning Algorithm

a)  b)

Figure 7: a) The flowchart of planning algorithm, b) Manipulation of the rope-
like ${\bm{\mathbf{L}}}_{o}$ curve with the designed virtual surface.Note that
following curves are obtained by using formulas in Phase I for $\psi_{q}$ and
$\delta$ in Eq. (16) and (17).

We plan the motion of the sphere spinning on the plane along the straight line
${\bm{\mathbf{L}}}_{s}$ with reaching a desired orientation. The flowchart in
Fig. 7-a shows our algorithm that uses the control inputs (15). Geometrically,
it can be considered as manipulating curvature radii $\\{R_{n},R_{g}\\}$ and
angular location $\zeta^{\prime}$ to bend the created rope-like curve
${\bm{\mathbf{L}}}_{o}$ from (15) on the sphere. For example, constant raise
in $\\{R_{n},R_{g}\\}$ expand the curve ${\bm{\mathbf{L}}}_{o}$ like enlarging
the loops of the rope while its length is same [see Fig. 7-b]. Also, constant
changes on $\zeta^{\prime}$ moves the curve ${\bm{\mathbf{L}}}_{o}$ on the
cutting plane of $u_{o,f}$. To do the convergence of this rope-like model, the
proposed algorithm will solve (7)-(8) and (15) in iterations with the re-tuned
constants in $\\{R_{n},R_{g}\\}$ and $\zeta^{\prime}$.

In general, to shift the curve of the sphere ${\bm{\mathbf{L}}}_{o}$ by
$\zeta^{\prime}(k)=\zeta_{q}(k)+\zeta_{u}(k)$ on the cutting plane $u_{o,f}$,
updates happen in Phase I and II by $\zeta_{q}$ and $\zeta_{u}$, respectively.
Also, constant change in the desired radius $R_{a}$ is defined by the
inclusion of all phases as $R_{a}(k)=R_{q}(k)+R_{u}(k)$, where $R_{q}$ and
$R_{u}$ are tuning constants of radii in the steps of the Phase II. Our tuning
algorithm consists of three primary steps:

The first step, Phase I, is the main part of the algorithm where the kinematic
model (7)-(8) with designed arc-length-based inputs (15) is numerically solved
in the time domain. After obtaining ${\bm{\mathbf{L}}}_{o}$ and
${\bm{\mathbf{L}}}_{s}$ curves, Phase I checks whether there is a point on
${\bm{\mathbf{L}}}_{o}$ that passes the desired local coordinates on the
sphere $\bm{\Psi}_{f}$. If ${\bm{\mathbf{L}}}_{o}$ fails to reach
$\bm{\Psi}_{f}$, $\zeta_{q}$ is re-tuned and kinematic model is resolved
numerically. After succeeding Phase I, the curve is checked whether the final
point of local coordinate $\bm{\Psi}_{l}$ is at desired values of
$\bm{\Psi}_{f}$ as Fig. 8. If the curve ${\bm{\mathbf{L}}}_{o}$ fails from
$\bm{\Psi}_{l}=\bm{\Psi}_{f}$, the radius of desired virtual surface $R_{q}$
is increased by obtained error. Also, Phase II has a second extra step to
reach the sphere final distance ${\bm{\mathbf{P}}}_{l}$ to the exact desired
plane configuration ${\bm{\mathbf{P}}}_{l}={\bm{\mathbf{P}}}_{f}$ by tuning
$R_{u}$ and $\zeta_{u}$ parameters. In the final step, the final state as the
desired spin angle $\psi_{f}$ is achieved. Note that this algorithm requires
to numerically solve differential equations in $k$ iterations until full
convergence of the chosen configuration.

####  4.1.1 Phase I

This phase solves the kinematic model (7)-(8) first. Then, it finds the curve
${\bm{\mathbf{L}}}_{o}$ on the rotating sphere that passes $\bm{\Psi}_{f}$
while the sphere moves toward desired position ${\bm{\mathbf{P}}}_{f}$.

Figure 8: Parametrization of the curve ${\bm{\mathbf{L}}}_{o}$ by cutting
planes of $u_{o,f}$ and $v_{o,f}$.

Before defining the tuning functions, we have to develop a converging state
for the spin angle. The spin motion can make tracking unstable when the sphere
passes the same desired spin angle many times during its rotation. Thus, we
develop an indirect convergence by using spinning angle deviation $\psi_{q}$
to converge the desired spin angle $\psi_{f}$. During calculations of minimum
distance $d$ at (1), we assumed that under-cap area $S_{t}$ of traversed
distance ${\bm{\mathbf{L}}}_{o}$ is constant while $C_{0}\rightarrow C_{1}$.
Thus, the desired under-cap area is defined by
$S_{t}=S^{\prime}_{c}+R^{2}_{o}(\psi_{f}+\psi_{u}(k)-{\psi^{\prime}})$ (only
orange and blue cap areas excluding dashed-line part as Fig. 3), where
$\psi_{u}(k)$ is the constant re-tuning spin angle for Phase III. Next, we
extend the Gauss-Bonnet theorem between two points of
$\\{u_{o}(t),v_{o}(t),\psi(t)\\}$ and $\\{u_{o,f},v_{o,f},\psi_{f}\\}$ in time
$t$, as

$\psi_{q}(t)=\frac{1}{R^{2}_{o}}\left[S_{t}-S_{i}(t)\right]=\frac{1}{R^{2}_{o}}\left[S_{t}-R^{2}_{o}\left(\psi(t)-\psi_{0}\right)\right],$
(16)

where $S_{i}(t)$ and $\psi(t)$ are the changing cap area from the initial
state till current time $t$ and the current spin angle.

As the last parameter in the kinematic model, the rolling rate (arc-length
derivation relative to the time $\delta(t)$) is the parameter that is
multiplied to all three inputs and drift term. This variable is related to the
time ($ds/dt$) and it varies the arc-length step of curve
${\bm{\mathbf{L}}}_{o}$ and ${\bm{\mathbf{L}}}_{s}$ in the given time $t$.
Therefore, $\delta$ is defined to have a rest-to-rest motion by

$\delta(t,u_{o},u_{s},v_{s})=||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}(t)||_{2}\cdot\Big{|}\frac{v_{o,f}\cdot
u^{\prime}_{o}}{T}\Big{|}=\left[(u_{s,f}-u_{s}(t))^{2}+(v_{s,f}-v_{s}(t))^{2}\right]^{\frac{1}{2}}\cdot\Big{|}\frac{v_{o,f}\cdot
u^{\prime}_{o}}{T}\Big{|},$ (17)

where $T$ is the time scaling.

Algorithm 1 Phase I

1:procedure
ConfigSolve($\bm{\Psi}_{f},{\bm{\mathbf{P}}}_{f},\psi_{f},\bm{\Psi}_{0},{\bm{\mathbf{P}}}_{0},\psi_{0}$)

2: while $e_{n}\geq\epsilon_{n}$ do

3: Solve the kinematic model (7) numerically

4: Collect $\bm{\Psi}_{n}$, $\bm{\Psi}_{l}$, $\bm{\Psi}_{v}$ and
$\bm{\Psi}_{u}$

5: Compute $e_{n}$ and $e^{\prime}_{n}$ by (18) $\triangleright$ Error of
nearest point to $\bm{\Psi}_{f}$

6: if $e_{n}(k)>\epsilon_{n}$ then

7: Apply the directional comparison by $\bm{\Psi}_{n}$, $\bm{\Psi}_{v}$ and
$\bm{\Psi}_{u}$ $\triangleright$ Details in C

8: Calculate $\zeta_{q}(k)$ by (19)

9: end if

10: end while

11:end procedure

Algorithm 1 shows the computations of Phase I. By solving the Darboux-based
kinematic model, the geometric parameters ([see Fig. 8) are determined from
obtained ${\bm{\mathbf{L}}}_{o}$ and ${\bm{\mathbf{L}}}_{s}$ trajectories.
$\bm{\Psi}_{l}$ and ${\bm{\mathbf{P}}}_{l}$ are the final arrived
configuration by the numerical solution of the found trajectories. Also,
$\bm{\Psi}_{n}$, $\bm{\Psi}_{u}$ and $\bm{\Psi}_{v}$ are the nearest point on
${\bm{\mathbf{L}}}_{o}$ to $\bm{\Psi}_{f}$ and existing nearest points on
${\bm{\mathbf{L}}}_{o}$ that are obtained by cutting planes of $u_{o,f}$ and
$v_{o,f}$ on the sphere. Note that all points are collected as the final
existing values since the curve can passes the nearest point or cutting planes
more than once.

In Algorithm 1, $\epsilon_{n}$ is a small value for breaking while loop, when
required accuracy is achieved for Phase I. Also, the error $e_{n}$ is the
Euclidean distance of nearest $\bm{\Psi}_{n}$ and desired $\bm{\Psi}_{f}$
points in $\mathbb{R}^{3}$

$\displaystyle\begin{split}&e_{n}=||\bm{\Psi}_{n}-\bm{\Psi}_{f}||_{3}=R_{o}\big{[}(\sin
v_{o,f}-\sin v_{o,n})^{2}+(\cos u_{o,f}\cos v_{o,f}-\cos u_{o,n}\cos
v_{o,n})^{2}\\\ &+(\sin u_{o,f}\cos v_{o,f}-\sin u_{o,n}\cos
v_{o,n})^{2}\big{]}^{\frac{1}{2}}.\end{split}$ (18)

Next, the absolute angle difference between $\bm{\Psi}_{f}$ and
$\bm{\Psi}_{n}$ are utilized to update the $\zeta_{q}(k)$ angle in k-th
iteration of Phase I as follows

$\begin{split}\zeta_{q}(k)=\zeta_{q}(k-1)\pm\begin{cases}&e^{\prime}_{n}\cdot|Q^{zx}_{f}-Q^{zx}_{n}|,\;\;\;\;\;Q^{zx}_{f}=Q^{zx}_{n}\\\
&e^{\prime}_{n}\cdot|Q^{zy}_{f}-Q^{zy}_{n}|,\;\;\;\;\;Q^{zx}_{f}\neq
Q^{zx}_{n}\end{cases}\end{split}$ (19)

where $e^{\prime}_{n}=\min\\{e_{n}\\}$ is smallest error of $e_{n}$ till
iteration $k$ and the angle differences with respect to $Z-X$ plane $Q^{zx}$
and $Z-Y$ plane $Q^{zy}$ for $\bm{\Psi}_{f}$/$\bm{\Psi}_{n}$ are obtained by
Eq. (29) as

$\displaystyle Q^{zx}=\left|\tan^{-1}\left(\frac{\sin u_{o}\cos v_{o}}{\cos
u_{o}\cos v_{o}}\right)\right|=|u_{o}|,Q^{zy}=\left|\tan^{-1}\left(\frac{\sin
u_{o}\cos v_{o}}{\sin v_{o}}\right)\right|.$ (20)

The sign of each update in (19) is chosen from the directional updates in C.
The goal of the directional update is to always move the curve
${\bm{\mathbf{L}}}_{o}$ towards $\bm{\Psi}_{f}$ on the spherical surface
$U_{C}$.

####  4.1.2 Phase II and III

Algorithm 2 Complete Computation

1:procedure
ConfigCom($\bm{\Psi}_{f},{\bm{\mathbf{P}}}_{f},\psi_{f},\bm{\Psi}_{0},{\bm{\mathbf{P}}}_{0},\psi_{0}$)

2: while $e_{s}\geq\epsilon_{s}$ do $\triangleright$ Phase III

3: while
$\left(e_{n}\geq\epsilon_{n}\right)\;\&\;\left(e_{r}\geq\epsilon_{r}\right)\;\&\;{\bm{\mathbf{P}}}_{l}\neq{\bm{\mathbf{P}}}_{f}\;$
do $\triangleright$ Phase II

4: while $e_{r}\geq\epsilon_{r}$ do $\triangleright$ Phase II

5: CONFIGSOLVE(.) $\triangleright$ Phase I

6: Compute $e_{r}$, $e^{\prime}_{r}$, $n_{s}$ and $d_{s}$

7: if $\bm{\Psi}_{l}\neq\bm{\Psi}_{n}$ then $\triangleright$ Updates $R_{q}$

8: $R_{q}(k)\leftarrow+R_{q}(k-1)$ according to (22)

9: else

10: $R_{q}(k)\leftarrow-R_{q}(k-1)$ by $n_{s}\leq 1$ case according to (22)

11: end if

12: end while

13: if $||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{l}||_{2}>\epsilon_{p}$ then
$\triangleright$ Phase II

14: Update $R_{u}$ and $\zeta_{u}$

15: end if

16: end while

17: if $\psi_{l}\neq\psi_{f}$ then $\triangleright$ Phase III

18: Update $\psi_{u}$

19: end if

20: Reset $e_{n}$, $e^{\prime}_{n}$, $e_{r}$ and $e^{\prime}_{r}$

21: end while

22:return ${\bm{\mathbf{L}}}_{o}$ and ${\bm{\mathbf{L}}}_{s}$

23:end procedure

After finding the suitable curve ${\bm{\mathbf{L}}}_{o}$ from Phase I that
passes $\bm{\Psi}_{f}$, Phase II tunes variable for converging
${\bm{\mathbf{L}}}_{o}$ and ${\bm{\mathbf{L}}}_{s}$ final arrived points
($\bm{\Psi}_{l}$ and ${\bm{\mathbf{P}}}_{l}$) to the desired configurations on
$\bm{\Psi}_{f}$ and ${\bm{\mathbf{P}}}_{f}$. Algorithm 2 depicts the overall
computations, including all the phases using the flowchart of Fig. 7-a.

In Phase II, the error $e_{r}$ of final arrival point $\bm{\Psi}_{l}$ with
respect to desire values of $\bm{\Psi}_{f}$ is first determined by

$\displaystyle\begin{split}&e_{r}=||\bm{\Psi}_{l}-\bm{\Psi}_{f}||_{3}=R_{o}\big{[}(\sin
v_{o,f}-\sin v_{o,l})^{2}+(\cos u_{o,f}\cos v_{o,f}-\cos u_{o,l}\cos
v_{o,l})^{2}\\\ &+(\sin u_{o,f}\cos v_{o,f}-\sin u_{o,l}\cos
v_{o,l})^{2}\big{]}^{\frac{1}{2}}.\end{split}$ (21)

Next, the update of $R_{q}$ is to enlarge the loops of ${\bm{\mathbf{L}}}_{o}$
for shorting the distance between nearest $\bm{\Psi}_{n}$ and final arrived
$\bm{\Psi}_{l}$ points; hence, the algorithm computes the $R_{q}$ update for
$k$-th iteration by using the error $e_{r}$ as

$\begin{split}R_{q}(k)=R_{q}(k-1)+R_{o}\cdot\begin{cases}&d_{s}\cdot
e^{\prime}_{r},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;n_{s}\leq 1\\\
&n_{s}\cdot d_{s}\cdot e^{\prime}_{r},\;\;\;\;\;\;\;1<n_{s}\leq 2\\\
&\left(d_{s}\cdot
e^{\prime}_{r}\right)/n_{s},\;\;\;\;\;\;\;\;\;\;\;n_{s}>2\end{cases}\end{split}$
(22)

where $e^{\prime}_{r}=\min\\{e_{r}\\}$, $d_{s}$ and $n_{s}$ are the smallest
error of $e_{r}$ till iteration $k$, distance ratio of ${\bm{\mathbf{L}}}_{s}$
for the loop on the curve ${\bm{\mathbf{L}}}_{o}$ and the numbers of created
full loops by spinning $\psi$, respectively. The distance ratio $d_{s}$ and
the number of created loops $n_{s}\in\mathbb{R}$ are calculated as follows

$d_{s}=\frac{||{\bm{\mathbf{P}}}_{n}-{\bm{\mathbf{P}}}_{f}||_{2}}{||{\bm{\mathbf{P}}}_{0}-{\bm{\mathbf{P}}}_{f}||_{2}}=\left[\frac{(u_{s,n}-u_{s,f})^{2}+(v_{s,n}-v_{s,f})^{2}}{(u_{s,0}-u_{s,f})^{2}+(v_{s,0}-v_{s,f})^{2}}\right]^{\frac{1}{2}},n_{s}=|\psi_{l}/\psi_{n}|.$
(23)

Note that $n_{s}$ is computed with the assumption that $\psi$ value can be
greater than $2\pi$ which every $2\pi$ orientation represents a loop on spin
angle [43].

As the final step of Phase II, because our arc-length-based inputs create
curve ${\bm{\mathbf{L}}}_{o}$ with multiple loops (specially for
$d\ll||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$), we use $R_{u}$ and
$\zeta_{u}$ variable for shifting ${\bm{\mathbf{L}}}_{o}$ away from
$\bm{\Psi}_{0}$ on the surface of the sphere $U_{C}$ (can be looked as
decreasing the size of the enlarged $R_{q}$) with same expressed flexible rope
model. This computation step with updating $\zeta_{u}$ and $R_{u}$ tries to
reach $\bm{\Psi}_{l}\rightarrow{\bm{\mathbf{P}}}_{f}$ as previous loops try to
keep $\bm{\Psi}_{l}=\bm{\Psi}_{f}$ condition true,

$\displaystyle\begin{split}&R_{u}(k)=R_{u}(k-1)-R_{q}(k)\cdot\frac{||{\bm{\mathbf{P}}}_{l}-{\bm{\mathbf{P}}}_{f}||_{2}}{||{\bm{\mathbf{P}}}_{l}-{\bm{\mathbf{P}}}_{0}||_{2}},\\\
&\zeta_{u}(k)=\zeta_{u}(k-1)+\begin{cases}&-\tan^{-1}\left(\frac{R_{u}(k)\tan\zeta_{q}(k-1)}{R_{o}}\right),\;\;\;\;\;\textnormal{
for }v_{o,f}\geq 0\\\
&\tan^{-1}\left(\frac{R_{u}(k)\tan\zeta_{q}(k-1)}{R_{o}}\right),\;\;\;\;\;\;\;\;\textnormal{
for }v_{o,f}<0\end{cases}\end{split}$ (24)

Final phase is to update the tuning variable of the spinning angle $\psi_{u}$
in Eq. (16) with

$\psi_{u}(k)=\psi_{u}(k-1)-e^{\prime}_{s}\operatorname{sgn}{\left[\psi_{f}-\psi_{l}\right]},$
(25)

where $e^{\prime}_{s}=\min\\{\psi_{f}-\psi_{l}\\}$ is the smallest error till
iteration $k$. Also, $\epsilon_{s}$, $\epsilon_{r}$ and $\epsilon_{p}$ are the
accuracies of each designated step in Algorithm 2.

##  5 Global Convergence of The Planning Algorithm

To prove the convergence of proposed iterative Algorithm 2, one can utilize
Zangwill’s convergence theorem [44].

######  Theorem 1.

Let the iterative algorithm
${\bm{\mathbf{A}}}:\;\mathbb{R}^{5}\rightarrow\mathbb{R}^{3}$, showed in the
flowchart at Fig. 7, be on
${\bm{\mathbf{X}}}={\bm{\mathbf{X}}}_{I}\circ{\bm{\mathbf{X}}}_{II}\circ{\bm{\mathbf{X}}}_{III}$
compact set where ${\bm{\mathbf{X}}}_{I}$, ${\bm{\mathbf{X}}}_{II}$ and
${\bm{\mathbf{X}}}_{III}$ are sets of phase I, II and phase III computation
steps, respectively. Given ${\bm{\mathbf{x}}}_{0}\in{\bm{\mathbf{X}}}$, the
created sequence $\\{{\bm{\mathbf{x}}}_{k}\\}^{\infty}_{k=1}$ satisfies
${\bm{\mathbf{x}}}_{k+1}\in{\bm{\mathbf{A}}}({\bm{\mathbf{x}}}_{k})$. Then,
${\bm{\mathbf{A}}}$ is globally convergent to a solution
$\bm{\Gamma}\subset{\bm{\mathbf{X}}}$ with following conditions

$\displaystyle\begin{split}&\bullet\exists\textit{ a descent function
${\bm{\mathbf{z}}}$ for }\bm{\Gamma}\textit{ and }{\bm{\mathbf{A}}},\\\
&\bullet\textit{The
squence}\;\\{{\bm{\mathbf{x}}}_{k}\\}^{\infty}_{k=0}\subset
S\;\textnormal{for}\;S\subset{\bm{\mathbf{X}}}\textit{ is a compact set,}\\\
&\bullet\textit{The mapping }{\bm{\mathbf{A}}}\textit{ is closed at all points
of }{\bm{\mathbf{X}}}/\bm{\Gamma}.\end{split}$

Regarding the first condition in Theorem 1, the designed line search algorithm
${\bm{\mathbf{A}}}$ have descent function
${\bm{\mathbf{z}}}=\\{R_{q},\zeta_{q},R_{u},\zeta_{u},\psi_{u}\\}$ which
roughly can be presented as

${\bm{\mathbf{z}}}(x)={\bm{\mathbf{x}}}_{k}+{\bm{\mathbf{h}}}_{k}{\bm{\mathbf{d}}}_{k}$
(26)

where ${\bm{\mathbf{d}}}_{k}$ and
${\bm{\mathbf{h}}}_{k}=\textnormal{diag}\\{h(1),h(2),h(3),h(4),h(5)\\}$ are
the direction of iteration and step size of iterations. Note that $\zeta_{u}$
is updated in the same way with others, only because $R_{u}$ changes the curve
${\bm{\mathbf{L}}}_{o}$, the new angular location shift happens by
$\zeta_{u}$. Next, the given descent function ${\bm{\mathbf{z}}}(x)$ has
following properties:

$\displaystyle\begin{cases}&\textnormal{If
}{\bm{\mathbf{x}}}\notin\bm{\Gamma}\textnormal{ and
}{\bm{\mathbf{y}}}\in{\bm{\mathbf{X}}},\;{\bm{\mathbf{z}}}(y)<{\bm{\mathbf{z}}}(x),\\\
&\textnormal{If }{\bm{\mathbf{x}}}\in\bm{\Gamma}\textnormal{ and
}{\bm{\mathbf{y}}}\in{\bm{\mathbf{X}}},\;{\bm{\mathbf{z}}}(y)\leq{\bm{\mathbf{z}}}(x),\end{cases}$

hence, the first condition of Theorem 1 is satisfied. The algorithm
${\bm{\mathbf{A}}}$ is on the compact manifold of the sphere $U_{C}$ and the
plane $U_{S}$ so ${\bm{\mathbf{h}}}_{k}$ norms as the error e.g., $e_{r}$,
$\mathbb{R}^{n}\rightarrow\mathbb{R}$ become a compact set as well. This can
be clearly extended for the second condition about all produced sequences
$\\{{\bm{\mathbf{x}}}_{k}\\}^{\infty}_{k=0}$ on ${\bm{\mathbf{X}}}$. Now, we
provide a proposition to satisfy final condition as

######  Proposition 1.

Let $f$ be a real continues function on ${\bm{\mathbf{X}}}$. Then, Algorithm
${\bm{\mathbf{A}}}$ with solution set $\bm{\Gamma}$ of

$\bm{\Gamma}({\bm{\mathbf{x}}},{\bm{\mathbf{d}}})=\\{{\bm{\mathbf{y}}}\in\mathbb{R}^{n}\;|\;{\bm{\mathbf{y}}}={\bm{\mathbf{x}}}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}},\;{\bm{\mathbf{h}}}>0,\;\textnormal{and
}f({\bm{\mathbf{y}}})=\min{\;f({\bm{\mathbf{x}}}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}})}\\}$

is closed at any point $({\bm{\mathbf{x}}},{\bm{\mathbf{d}}})$ at which always
${\bm{\mathbf{d}}}\neq 0$.

D presents the proof of the given proposition. By satisfying all the three
conditions in Theorem 1, the proposed line search algorithm is always
convergent to a solution for the desired configuration.

Figure 9: Simulation results for a case study with the final configuration
$\\{u_{s,f},v_{s,f},u_{o,f},v_{o,f},\psi_{f}\\}=\\{3,3.2,-\frac{\pi}{2}-0.8,0.8,0.8\\}$.

##  6 Results and Discussion

The operation and performance of the proposed planning approach are tested
under simulations. First, we check how phases of the algorithm work to achieve
a successful convergence. Next, we analyze the simulation results for
different desired spin angles while the rest of the configuration is the same.

1\.
2\.
3\.
4\.

Figure 10: Collected iterations of the given final configuration: 1. The first
solution of system ($k=0$), 2. The first successful results of Phase I
($k=11$), 3. Converging Phase II ($k=14$), 4. The first results before
applying Phase III ($k=41$).

In the simulations, the kinematic model (7)-(8) is solved together with (15)
and (17) by using Matlab’s ODE45 routine. We set the initial and final
configurations as
$\\{u_{s,0}$,$v_{s,0}$,$u_{o,0}$,$v_{o,0}$,$\psi_{0}\\}=\\{0$,$0$,$0$,$0$,$0\\}$
and $\\{u_{s,f}$, $v_{s,f}$, $u_{o,f}$, $v_{o,f}$, $\psi_{f}$$\\}=\\{3$,
$3.2$, $-\frac{\pi}{2}-0.8$, $0.8$ ,$0.8\\}$. Note that our plane’s final
configuration satisfies distance constraint (1) where $d$ for our case is
$2.15\leq||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}$, and to see
curve with multiple loops, we have
$||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}_{0}||_{2}=4.38$ m. For this case,
the simulation time is set to $t_{f}=15$ s and time constant $T$ in (17) is 1.
The accuracy of phases in Algorithm 2 are $\epsilon_{n}=\epsilon_{r}=0.07$,
$\epsilon_{p}=0.12$ and $\epsilon_{s}=0.05$. As expressed in the
controllability analysis in Ref. [18], there is a uncontrollable point at
$\pi/4$ on the plane. For this case, we chose it away from point with
$G_{f}=0.754$ rad to see our approach abilities. The initial value of $R_{q}$
is 0.005 to prevent any computation singularities. Also, we set the sphere
radius and maximum divider as $R_{o}=0.5$ m and $\mu_{r}=4$.

Figure 11:  Angular velocities of the sphere with the designed arc-length-
based based controller.

After running the simulation, the final results are shown as Fig. 9. The video
of the simulation results for this section is available on Youtube [45]. The
calculation steps $k$ for this case is achieved with 48 iterations. The
computation process takes about 7.2 s. By the achieved successful final
iteration, we can interpret that while the sphere is rolling along the given
desired final position ${\bm{\mathbf{P}}}_{f}$, it spins with smooth
trajectories toward its final configuration for
$\\{\bm{\Psi}_{f},\psi_{f}\\}$.

a)  b)

Figure 12:  a) The designed function $\omega_{s}$ in time, b) Smoothed angular
velocities of the sphere in two different cases of time scaling $T$.

To understand better how the algorithm works by utilizing the derived virtual
surface, we plot four collected iteration steps [see Fig. 10]. From the
initial iteration ($k=0$), it is clear that ${\bm{\mathbf{L}}}_{o}$ is created
by the designed virtual surface toward the desired states $\bm{\Psi}_{f}$ but
because ${\bm{\mathbf{P}}}_{f}$ is far, ${\bm{\mathbf{L}}}_{o}$ curve moves
along the $u_{o,f}$ cutting plane. A rope on a sphere explains the way that
the achieved curve ${\bm{\mathbf{L}}}_{o}$ on $U_{C}$ can be imagined. So, the
algorithm tunes the virtual surface and this bends the obtained curve
${\bm{\mathbf{L}}}_{o}$ in loops, like a rope, toward the desired
$\bm{\Psi}_{f}$. Here, Phase I rises the ${\bm{\mathbf{L}}}_{o}$ to move
$\bm{\Psi}_{n}$ near to $\bm{\Psi}_{f}$ as the result of succeeded iteration
$k=11$. To move $\bm{\Psi}_{l}$ on $\bm{\Psi}_{f}$, Phase II enlarges the
radius of the rope-like curve ${\bm{\mathbf{L}}}_{o}$, as the example
iteration $k=14$ shows. The iteration $k=41$ in Fig. 10 satisfies
$\bm{\Psi}_{l}=\bm{\Psi}_{f}$ condition while the next step is tuning the
plain configuration to achieve ${\bm{\mathbf{P}}}_{l}={\bm{\mathbf{P}}}_{f}$
using Eq. (24). Finally, the Phase III condition is satisfied for the desired
$\psi_{f}$ value shown in Fig. 9.

Figure 13:  Comparing angular convergence of $\bm{\Psi}_{f}$ at two cases.

Figure 14: Simulation results for different final spin angles
$\psi_{f}=\\{-1.7,1.3,2.3\\}$.

The angular velocities of the sphere collected by search algorithm are
presented in Fig. 11 but there are certain fast accelerations in the angular
velocities. Note that the problem is considered in kinematics level without
any slipping; Hence, arriving longer distances in a short amount of time
results in faster velocities. However, we can use the aforementioned advantage
in the derived kinematic model (7) where the time domain is separated from the
kinematics by the rolling rate $\delta$. Thus, We alter the sphere angular
velocities in two cases:

In the initial case, to decrease the overall velocity, the time span can be
enlarged without affecting the achieved paths ${\bm{\mathbf{L}}}_{o}$ and
${\bm{\mathbf{L}}}_{s}$. For the considered example, we increase the time
scaling constant $T$ (it is at $\delta$ in Eq. (17) ) and simulation time
$t_{f}$ to 160 s while the kinematic model is solved with the succeeded
parameters in $\zeta^{\prime}$ and $R_{q}+R_{u}$ from the algorithm. Clearly,
the overall velocities decreased by expanding them in time as Fig. 12-b
(dashed line velocities). However, there are fast accelerations (around 3 s
and 21 s ) that makes it harder to be applied practically. To solve this
issue, we keep the simulation time $t_{f}$ same as 160 s but a smooth function
is chosen for $T$ as the second case study. With the known relation of
$\omega^{o}_{z}=\delta k^{*}_{g}\overset{\Delta}{=}\delta\alpha_{s}$ [18] and
Eq. (17), we define $T$ as

$\displaystyle T\overset{\Delta}{=}(c\cdot\omega_{s})/\alpha_{s},$ (27)

where $c=\left|v_{o,f}\cdot
u^{\prime}_{o}\right|\cdot||{\bm{\mathbf{P}}}_{f}-{\bm{\mathbf{P}}}(t)||_{2}$
and $\omega_{s}$ is the desired smooth angular velocity. The desired smooth
velocity $\omega_{s}$ is created with a symmetric time function as

$\displaystyle\omega_{s}=a\left(-\frac{140}{T_{s}^{7}}t^{7}+\frac{420}{T_{s}^{6}}t^{6}-\frac{420}{T_{s}^{5}}t^{5}+\frac{140}{T_{s}^{4}}t^{4}\right),$
(28)

where $a=12.91$ and $T_{s}=160$ are the time constant and the amplitude. Fig.
12-a shows the rest-to-rest function behavior of Eq. (28). After running the
same simulation with our defined $T$, we get the results in Fig. 12-b (solid
lines velocities). It is clear that the second case has smooth angular
velocities which are practical without large accelerations at the same time
span of $t_{f}$. The angular displacement in Fig. 13 also demonstrate how the
desired angles of the sphere get expanded with the same convergence trend
(${\bm{\mathbf{L}}}_{o}$ trajectory is the same as Fig. 9 ). Note that more
complex functions can be defined for $T$ to converge with different angular
velocities without affecting the trajectories of ${\bm{\mathbf{L}}}_{o}$ and
${\bm{\mathbf{L}}}_{s}$.

Next, we analyze the patterns of the created curves ${\bm{\mathbf{L}}}_{o}$
for different desired spin angles with the same final configuration. In here,
the final configuration is considered $\\{u_{s,f}$, $v_{s,f}$, $u_{o,f}$,
$v_{o,f}\\}=\\{3,3.2,0.6,0.7\\}$ besides different desired spin angles
$\psi_{f}=\\{-1.7,1.3,2.3\\}$. To have these different curves approximately on
same configuration, the $\epsilon_{n}$ and $\epsilon_{r}$ accuracies are
chosen $\\{0.045,0.05\\}$, $\\{0.048,0.034\\}$ and $\\{0.048,0.05\\}$ for
$\psi_{f}=\\{-1.7,1.3,2.3\\}$ cases, in the given order. Also, simulation time
is set for $t_{f}=20$ s with $T=1$. Note that the rest of the initial
conditions are the same as the previous case study.

Simulation results for different final spin angles are presented in Fig. 14.
As can be seen, the algorithm converges successfully to desired states with
different trajectories. The average process time for these simulations was
about 40 s. Note that the rise in the average process time is due to
exceptional computation for $\psi_{f}=-1.7$ rad. This occurs depending on the
accuracy and how near the desired $\bm{\Psi}_{f}$ angles are to
$\bm{\Psi}_{0}$ which causes larger iterations (110 s process time).

Also, for the case of $\psi_{f}=-1.7$ rad, the curve creates multiple loops to
reach the desired goal. As mentioned before, the reason is that some spin
angles exhibit slower convergence by the designed approach. To achieve a
faster convergence, we have to change the region of attraction for the desired
configuration. This can be done by adding $+\frac{\pi}{2}$ to $v$-curve for
the $\alpha_{s}(v_{o,f}+\frac{\pi}{2})$ arc-length-based input for this
configuration. Since we have already guarantied the convergence of the desired
goal by the Theorem 1, we can have regional shifts in the $v_{o,f}$ and
$u_{o,f}$. Note that this is an interesting property of proposed algorithm
since multiple solutions would be achievable for the same configuration with
changing the time constant $T$, resulting to different regions of attraction
on $u$\- and $v$-curves.

##  7 Conclusion

In this paper, we presented a new approach for path planning of the spin-
rolling sphere on a straight plane trajectory. The planning was mainly
constructed by the introduced geometric control. In the beginning, the path
planning problem is explained. Next, a virtual surface was designed to produce
arc-length-based control inputs for the fully-actuated transformed model in
[18]. The problem is imagined as manipulating a flexible rope-like curve on
the sphere with a constant length. Thus, we proposed an iterative algorithm to
tune the traversed curve toward desired configurations. The achieved
simulations clarified the performance of the proposed planning approach.

The proposed algorithm produces smooth trajectories of the spin-rolling
sphere. In this connection, it should be noted that the feedback-based
planning algorithms [24, 25, 26, 28] results in piece-wise smooth
trajectories. Also, our approach does not require decomposition of the
planning strategy into motion steps as featured in geometric phase shift
methods [30, 28, 24, 38]. Another advantage is that the resulting angular
velocities of the spin-rolling sphere can shaped to different desired
convergence rates.

## References

  * [1] B. Kiss, J. Lévine, B. Lantos, On motion planning for robotic manipulation with permanent rolling contacts, Int. J. Robot. Res. 21 (5-6) (2002) 443–461.
  * [2] L. Cui, J. Sun, J. Dai, In-hand forward and inverse kinematics with rolling contact, Robotica 35 (12) (2017) 2381–2399.
  * [3] S. Yuan, L. Shao, C. L. Yako, A. Gruebele, J. K. Salisbury, Design and control of roller grasper v2 for in-hand manipulation, in: IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2020, pp. 9151–9158.
  * [4] S. A. Tafrishi, ”RollRoller” novel spherical mobile robot basic dynamical analysis and motion simulations, Master’s thesis, University of Sheffield, Sheffield, UK (2014).
  * [5] S. A. Tafrishi, M. Svinin, E. Esmaeilzadeh, M. Yamamoto, Design, modeling, and motion analysis of a novel fluid actuated spherical rolling robot, ASME J. Mech. Robot. 11 (4) (2019) 041010.
  * [6] M. Ishikawa, R. Kitayoshi, T. Sugie, Volvot: A spherical mobile robot with eccentric twin rotors, in: Proc. IEEE Int. Conf. Robot. Biomimetics, 2011, pp. 1462–1467.
  * [7] Y. Karavaev, I. Mamaev, A. Kilin, E. Pivovarova, Spherical rolling robots: Different designs and control algorithms, in: Robots in Human Life: Proc. of the 23rd Int. Conf. on Climbing and Walking Robots and the Support Technologies for Mobile Machines (CLAWAR 2020, Moscow, Aug 24–26 2020), 2020, pp. 195–202.
  * [8] A. Borisov, A. Kilin, I. Mamaev, How to control chaplygin sphere using rotors, Regular and Chaotic Dynamics 17 (3) (2012) 258–272.
  * [9] A. Borisov, A. Kilin, I. Mamaev, How to control the chaplygin ball using rotors. ii, Regular and Chaotic Dynamics 18 (1) (2013) 144–158.
  * [10] M. Svinin, A. Morinaga, M. Yamamoto, On the dynamic model and motion planning for a spherical rolling robot actuated by orthogonal internal rotors, Regular and Chaotic Dynamics 18 (1) (2013) 126–143.
  * [11] P. Fankhauser, C. Gwerder, Modeling and control of a ballbot, B.S. thesis, Eidgenössische Technische Hochschule Zürich (2010).
  * [12] J. Johnson, R. Senthilnathan, M. Negi, R. G. Patel, A. Bhattacherjee, A fuzzy logic-in-loop control for a novel reduced height ballbot prototype, Procedia Comput. Sci. 133 (2018) 960–967.
  * [13] B. Sümer, M. Sitti, Rolling and spinning friction characterization of fine particles using lateral force microscopy based contact pushing, J. Adhes Sci. Technol. 22 (5-6) (2008) 481–506.
  * [14] E. Diller, M. Sitti, et al., Micro-scale mobile robotics, Found. Trends Robot., 2 (3) (2013) 143–259.
  * [15] A. Fernández-Pacheco, R. Streubel, O. Fruchart, R. Hertel, P. Fischer, R. P. Cowburn, Three-dimensional nanomagnetism, Nat. Commun. 8 (2017) 15756.
  * [16] V. Jurdjevic, The geometry of the plate-ball problem, Arch. Ratio. Mech. Anal. 124 (4) (1993) 305–328.
  * [17] A. Marigo, A. Bicchi, Rolling bodies with regular surface: Controllability theory and applications, IEEE Trans. Autom. Control 45 (9) (2000) 1586–1599.
  * [18] S. A. Tafrishi, M. Svinin, M. Yamamoto, Darboux-frame-based parametrization for a spin-rolling sphere on a plane: A nonlinear transformation of underactuated system to fully-actuated model, Mechanism and Machine Theory 164 (2021) 104415\.
  * [19] N. Sankar, V. Kumar, X. Yun, Velocity and acceleration analysis of contact between three-dimensional rigid bodies.
  * [20] D. J. Montana, The kinematics of contact and grasp, Int. J. Robot. Res. 7 (3) (1988) 17–32.
  * [21] L. Cui, J. S. Dai, A darboux-frame-based formulation of spin-rolling motion of rigid objects with point contact, IEEE Trans. Robot. 26 (2) (2010) 383–388.
  * [22] I. Bizyaev, A. Borisov, I. Mamaev, Different models of rolling for a robot ball on a plane as a generalization of the chaplygin ball problem, Regular and Chaotic Dynamics 24 (5) (2019) 560–582.
  * [23] Z. Woodruff, K. Lynch, Second-order contact kinematics between three-dimensional rigid bodies, Journal of Applied Mechanics 86 (8).
  * [24] H. Date, M. Sampei, M. Ishikawa, M. Koga, Simultaneous control of position and orientation for ball-plate manipulation problem based on time-state control form, IEEE Trans. Robot. 20 (3) (2004) 465–480.
  * [25] G. Oriolo, M. Vendittelli, A framework for the stabilization of general nonholonomic systems with an application to the plate-ball mechanism, IEEE Trans. Robot. 21 (2) (2005) 162–175.
  * [26] T. Das, R. Mukherjee, Exponential stabilization of the rolling sphere, Automatica 40 (11) (2004) 1877–1889.
  * [27] F. Alouges, Y. Chitour, R. Long, A motion-planning algorithm for the rolling-body problem, IEEE Trans. Robot. 5 (26) (2010) 827–836.
  * [28] R. Mukherjee, T. Das, Feedback stabilization of a spherical mobile robot, in: Proc. IEEE/RSJ Int. Conf. on Intell. Robots Sys., Vol. 3, 2002, pp. 2154–2162.
  * [29] Z. Woodruff, S. Ren, K. Lynch, Motion planning and feedback control of rolling bodies, IEEE Access 8 (2020) 31780–31791.
  * [30] Z. Li, J. Canny, Motion of two rigid bodies with rolling constraint, IEEE Trans. Robot. Autom. 6 (1) (1990) 62–72.
  * [31] R. Mukherjee, M. A. Minor, J. T. Pukrushpan, Motion planning for a spherical mobile robot: Revisiting the classical ball-plate problem, ASME Journal of Dynamic Systems, Measurement, and Control 124 (4) (2002) 502–511.
  * [32] M. Svinin, S. Hosoe, Motion planning algorithms for a rolling sphere with limited contact area, IEEE Trans. Robot. 24 (3) (2008) 612–625.
  * [33] A. Kilin, E. Pivovarova, T. Ivanova, Spherical robot of combined type: Dynamics and control, Regular and chaotic dynamics 20 (6) (2015) 716–728.
  * [34] Y. Bai, M. Svinin, M. Yamamoto, Dynamics-based motion planning for a pendulum-actuated spherical rolling robot, Regular and Chaotic Dynamics 23 (4) (2018) 243–259.
  * [35] A. Arthurs, G. Walsh, On hammersley’s minimum problem for a rolling sphere, in: Math. Proc. Cambridge Philos. Soc., Vol. 99, Cambridge University Press, 1986, pp. 529–534.
  * [36] Y. Sachkov, Maxwell strata and symmetries in the problem of optimal rolling of a sphere over a plane, Sb. Math 201 (7) (2010) 1029–1051.
  * [37] A. Mashtakov, Y. Sachkov, Extremal trajectories and the asymptotics of the maxwell time in the problem of the optimal rolling of a sphere on a plane, Sb. Math 202 (9) (2011) 1347–1371.
  * [38] A. Morinaga, M. Svinin, M. Yamamoto, A motion planning strategy for a spherical rolling robot driven by two internal rotors, IEEE Trans. Robot. 30 (4) (2014) 993–1002.
  * [39] I. Beschatnyi, The optimal rolling of a sphere, with twisting but without slipping, Mat. Sb. 205 (2) (2014) 157.
  * [40] L. Cui, J. Dai, Sliding-rolling Contact and In-hand Manipulation, World Scientific, 2020.
  * [41] M. P. do Carmo, Differential Geometry of Curves and Surfaces, 2nd Edition, Prentice-Hall, 1976.
  * [42] E. Cartan, Riemannian Geometry in an Orthogonal Frame, 1st Edition, World Scientific Pub Co Inc, 2002.
  * [43] K. Tapp, Differential geometry of curves and surfaces, Springer, New York, 2016\.
  * [44] W. Zangwill, Nonlinear programming: a unified approach, Vol. 196, Prentice-Hall Englewood Cliffs, NJ, 1969.
  * [45] A. Tafrishi, M. Svinin, M. Yamamoto, Y. Hirata, Path planning of spin-rolling sphere on a plane (November 2021).
URL {https://www.youtube.com/watch?v=7JK3VtT0aqs}

##  Appendix A Ball-Plate System Preliminaries

Figure 15: Kinematic model of rotating sphere. Note: $\psi$ is the spin angle
between sphere and plane surfaces.

Fig. 15 illustrates the rotating object and surface coordinates in the ball-
plate system. Here, we have $\Sigma_{o}$ and $\Sigma_{s}$ as the fixed frames
on the rolling object (red sphere) and the plane. The frame $\Sigma_{s}$ is
fixed relative to other coordinate frames. It is assumed that the sphere with
radius $R_{o}$ is rotating with no sliding constraint. The local coordinate
systems for the sphere and plane are considered

$\begin{split}&f_{o}:U_{C}\rightarrow\mathbb{R}^{3}:c(u_{o},v_{o})\mapsto(-R_{o}\sin{u_{o}}\cos{v_{o}},R_{o}\sin{v_{o}},-R_{o}\cos{u_{o}}\cos{v_{o}}),\\\
&f_{s}:U_{S}\rightarrow\mathbb{R}^{3}:c(u_{s},v_{s})\mapsto(u_{s},v_{s},0),\end{split}$
(29)

where $c(u_{o},v_{o})\in[-\pi,\pi]$ and $c(u_{s},v_{s})$ are contact
parameters of the sphere and plane. For the sphere we have [30]:

$k^{o}_{nu}=k^{o}_{nv}=1/R_{o},\;\tau^{o}_{gu}=\tau^{o}_{gv}=0,\;k^{o}_{gu}=\tan(v_{o})/R_{o},\;k^{o}_{gv}=0,$
(30)

where $k^{o}_{nu}$, $k^{o}_{nv}$, $\tau^{o}_{gu}$, $\tau^{o}_{gv}$,
$k^{o}_{gu}$, $k^{o}_{gv}$ are the normal curvature, geodesic torsion and
geodesic curvature of the rotating body respect to $u_{o}$ and $v_{o}$
principle angle. The curvature terms of the plane surface $U_{S}$ are

$\displaystyle\begin{split}&k^{s}_{nu}=k^{s}_{nv}=\tau^{s}_{gu}=\tau^{s}_{gv}=k^{s}_{gu}=k^{s}_{gv}=0,\end{split}$
(31)

where $k^{s}_{nu}$, $k^{s}_{nv}$, $\tau^{s}_{gu}$, $\tau^{s}_{gv}$,
$k^{s}_{gu}$, $k^{s}_{gv}$ are the normal curvature, geodesic torsion and the
geodesic curvature of the fixed surface (plane) respect to $u_{s}$ and $v_{s}$
principle angles.

##  Appendix B Geodesic Torsion Design of the Virtual Surface

The rotating object and plane do not have the geodesic torsion,
$\tau^{*}_{g}=-\beta_{s}$. It makes the kinematic model (7) uncontrollable
[18]. Thus, a helicoid virtual surface with similar curvature properties with
rotating object (see Fig. 5) is proposed as

$\displaystyle\begin{split}&f_{v}:U_{V}\rightarrow
R^{3}:c(u_{v},v_{v})\mapsto(-R_{v}\sin{u_{v}}\cos{v_{v}},R_{v}\sin{v_{v}}+R_{t}u_{v},-R_{v}\cos{u_{v}}\cos{v_{v}}),\end{split}$
(32)

where $R_{v}$ and $R_{t}$ are defined by main spherical and sum of spherical
and torsion radii, respectively. The curvature properties can be obtained as

$\displaystyle
k^{v}_{n}=\frac{1}{R_{v}},\;k^{v}_{g}=\frac{R_{v}\cos{v_{v}}\sin{v_{v}}}{R^{2}_{v}\cos^{2}{v_{v}}+R^{2}_{t}},\tau^{v}_{g}=\frac{1}{R_{v}^{2}}(R_{v}^{2}\cos^{2}{v_{v}}+R^{2}_{t})^{\frac{1}{2}}$
(33)

where $k^{v}_{n},\;k^{v}_{g},\;\tau^{v}_{g}$ are normal curvature, geodesic
curvature and geodesic torsion. We here mainly care about
$\tau^{v}_{g}(R_{v},R_{t},v_{v})$ and $k^{v}_{n}(R_{v})$ to understand
relation of geodesic torsion design. Thats why we related them with separate
corresponding radii. To prove that this surface let us manipulate $\tau_{g}$
by using $R_{t}$, we consider $R_{t}=nR_{v}$ where $n>1$, which results

$\displaystyle\tau^{v}_{g}=\frac{1}{R_{v}}(\cos^{2}{v_{v}}+n^{2})^{\frac{1}{2}}.$
(34)

######  Theorem 2.

At an arbitrary point $P$ on a surface, geodesic torsion relation with normal
curvature is defined as [42]

$\displaystyle\tau_{g}=\frac{1}{2}(k_{nu}-k_{nv})\sin 2\Omega$ (35)

where $k_{nu}$ and $k_{nv}$ are the principle curvatures. Also, $\Omega$ is
the counterclockwise angle from the direction of minimum curvature $k_{nu}$ on
the tangent plane.

By relying on the Theorem 2, it is clear that sphere tangent plane angle
$\Omega$ is $\pi/4$. Now, we can do some algebraic operations on (35) with
considering $k_{nu}=1/R_{v}$ and Eq. (34),

$k_{nv}=1/R_{v}\left[1+(\cos^{2}{v_{v}}+n^{2})^{\frac{1}{2}}\right]$ (36)

where $k_{nv}$ is the second principle curvature of the virtual surface. This
supports our assumption that designed geodesic torsion is inverse of sphere
radius with similar unity of normal curvature $R_{t}=nR_{v}$. We do this for
sake of separating the radius of geodesic torsion and normal curvature in Eq.
(11) as we use two separate surfaces (sphere and helicoid) in arc-length-based
inputs.

##  Appendix C Phase I Directional Update

We can separate the sphere regions to four by $c(u_{o,f},v_{o,f})$ cutting
planes (see example in Fig. 16) for the desired configuration where
$\zeta_{q}$ is calculated for its directional updates.

Figure 16: The sectioned sphere based on the $\bm{\Psi}_{v}$ and
$\bm{\Psi}_{u}$ cutting planes for directional updates of $\zeta_{q}$. Note
that $\bm{\Psi}_{v}$ and $\bm{\Psi}_{u}$ are on $\varrho_{2}$ and
$\varrho_{3}$ as an example.

Because we need to find directional updates in different $G_{f}$ angles on
plane, we apply a spherical rotational transformation to ease our regional
computation. To find corresponding local points’ right location for the
directional update from plane angle $G_{f}$ respect to $v_{s}$, we use an
arbitrary local coordinate $c(u_{o},v_{o})$ rotation respect to $-\pi/4$ base
angle, $G^{\prime}_{f}=G_{f}-\pi/4$, where gives new coordinate
$c(u^{r}_{o},v^{r}_{o})$ as

$\begin{split}&v^{r}_{o}=\sin^{-1}\big{[}-\sin G^{\prime}_{f}\sin u_{o}\cos
v_{o}+\cos G^{\prime}_{f}\sin v_{o}\big{]},\;u^{r}_{o}=\sin^{-1}\big{[}(\cos
G^{\prime}_{f}\sin u_{o}\cos v_{o}+\sin G^{\prime}_{f}\sin v_{o})/\cos
v^{r}_{o}]\\\ \end{split}$ (37)

Next, to avoid the numerical solution of following Eqs. (37), we present it in
algebraic form depending on $c(u_{o},v_{o})$ location by

$\begin{cases}&\begin{cases}&u^{r}_{o}\leftarrow-u^{r}_{o},\\\
&v^{r}_{o}\leftarrow\pi-v^{r}_{o},\end{cases}\;\;\textnormal{for}\;\Big{[}\left(0\leq|v_{o}|\leq\frac{\pi}{2}\;\&\;\frac{\pi}{2}\leq|u_{o}|\leq\pi\right)||\left(\frac{\pi}{2}<|v_{o}|\leq\pi\;\&\;0\leq|u_{o}|<\frac{\pi}{2}\right)\Big{]}\\\
&\begin{cases}&u^{r}_{o}\leftarrow u^{r}_{o},\\\ &v^{r}_{o}\leftarrow
v^{r}_{o},\end{cases}\;\;\;\;\;\;\;\;\textnormal{for}\;\Big{[}\left(0<|v_{o}|<\frac{\pi}{2}\;\&\;0<|u_{o}|<\frac{\pi}{2}\right)||\left(\frac{\pi}{2}<|v_{o}|<\pi\;\&\;\frac{\pi}{2}<|u_{o}|<\pi\right)\Big{]}\\\
\end{cases}$ (38)

The operations (37)-(38) are applied to rotate $\bm{\Psi}_{f}$,
$\bm{\Psi}_{u}$, $\bm{\Psi}_{v}$ and $\bm{\Psi}_{n}$ relative to angle $G_{f}$
to find them always in one $G^{\prime}_{f}$ direction. Now, the sectioned
sphere as Fig. 16 is obtained by cutting planes $\bm{\Psi}_{v}$ and
$\bm{\Psi}_{u}$ where they create $\\{\varrho_{1},\varrho_{2}\\}$ and
$\\{\varrho_{3},\varrho_{4}\\}$ regions, respectively. Note that for nearest
point $\bm{\Psi}_{n}$ on ${\bm{\mathbf{L}}}_{o}$, the corresponding $\varrho$
are shown as $\varrho_{n1}$, $\varrho_{n2}$, $\varrho_{n3}$, $\varrho_{n4}$.
We design our directional updates by comparing $\bm{\Psi}_{n}$ with
$\bm{\Psi}_{u}$, $\bm{\Psi}_{v}$ as following computation

Calculate $G^{\prime}_{f}$ depending on the $\\{u_{s,f},v_{s,f}\\}$

Calculate rotated coordinates of $\bm{\Psi}_{f}$, $\bm{\Psi}_{u}$,
$\bm{\Psi}_{v}$ and $\bm{\Psi}_{n}$ by (38)

Calculate $Q^{zy}_{f}$, $Q^{zx}_{f}$, $Q^{zy}_{n}$ and $Q^{zx}_{n}$ according
to (20)

if $\varrho_{2}=1\;\&\;\varrho_{4}=1$ then $\triangleright$ Exceptional
regional updates

if only $\varrho_{n2}=1$ then

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

else if $\left(\varrho_{n2}=1\;\&\;\varrho_{n4}=1\right)$ $||$ $($only
$\varrho_{n4}=1)$ then

if $v^{r}_{o,f}\geq 0$ then

if $u^{r}_{o,f}\geq 0$ then

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

else

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

end if

else

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

end if

end if

else if
$\left(\varrho_{n2}=1\;\&\;\varrho_{n3}=1\right)||\left(\varrho_{3}=1\;\&\;\varrho_{n2}=1\;\&\;\varrho_{n4}=1\right)$
then

if $u^{r}_{o,f}\geq 0$ then

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

else

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

end if

else if $\varrho_{1}=\varrho_{2}=\varrho_{3}=\varrho_{4}=0$ then

if $v^{r}_{o,f}\geq 0$ then

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

else

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

end if

else$\triangleright$ Normal regional updates

if only $\left(\varrho_{n1}=1\right)$ $||$ $\left(\varrho_{n4}=1\right)$ then

if $v^{r}_{o,f}\geq 0$ then

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

else

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

end if

else if only $\left(\varrho_{n2}=1\right)$ $||$ $\left(\varrho_{n3}=1\right)$
then

if $v^{r}_{o,f}\geq 0$ then

$\zeta_{q}(k)\leftarrow-\zeta_{q}(k-1)$

else

$\zeta_{q}(k)\leftarrow+\zeta_{q}(k-1)$

end if

end if

end if

where $\varrho_{ni}$ presents the nearest point $\bm{\Psi}_{n}$ location on
$i$-th section of the spherical surface ($i\in[1,4]$) with respect to
$\bm{\Psi}_{f}$. Also, $\varrho_{ni}=1$ means the following argument is true.
Note that this complexity in sectioned spherical surface, as shown in Fig. 16,
is due to rolling surface that creates different signs of updates of
$\zeta_{q}(k)$ in (19) based on
$\\{\bm{\Psi}_{0},\bm{\Psi}_{f},\bm{\Psi}_{n}\\}$. In additions, the normal
updates in the algorithm find $\bm{\Psi}_{n}$ location on the spherical
surface and compare it to the location of desired goal $\bm{\Psi}_{f}$ to give
correct sign update for $\zeta_{q}(k)$. However, exceptional updates are
designed to exist from trajectory loops of ${\bm{\mathbf{L}}}_{o}$ for
$\bm{\Psi}_{n}$ and bring the ${\bm{\mathbf{L}}}_{o}$ curve toward the
location of $\bm{\Psi}_{0}$.

##  Appendix D Proof of Closed Set

Let ${\bm{\mathbf{x}}}_{k}\rightarrow{\bm{\mathbf{x}}}_{0}$ and
${\bm{\mathbf{d}}}_{k}\rightarrow{\bm{\mathbf{d}}}_{0}$ be the convergence to
desired values as $k\rightarrow\infty$ via sequences of
$\\{{\bm{\mathbf{x}}}_{k}\\}^{\infty}_{k=1}$ and
$\\{{\bm{\mathbf{d}}}_{k}\\}^{\infty}_{k=1}$. Also, suppose
$\\{{\bm{\mathbf{y}}}_{k}\\}^{\infty}_{k=1}$ is a sequence where
${\bm{\mathbf{y}}}\in\Gamma({\bm{\mathbf{x}}}_{k},{\bm{\mathbf{d}}}_{k})$ for
all $k$ while it has convergence of
${\bm{\mathbf{y}}}_{k}\rightarrow{\bm{\mathbf{y}}}_{0}$ as
$k\rightarrow\infty$. We want to illustrate that
${\bm{\mathbf{y}}}_{0}\in\Gamma({\bm{\mathbf{x}}}_{0},{\bm{\mathbf{d}}}_{0})$
for having a closed set.

In each iteration $k$,
${\bm{\mathbf{y}}}_{k}={\bm{\mathbf{x}}}_{k}+{\bm{\mathbf{h}}}_{k}{\bm{\mathbf{d}}}_{k}$
for ${\bm{\mathbf{h}}}_{k}>0$. Thus, by knowing fact that
$||{\bm{\mathbf{d}}}_{k}||=[1\;..\;1]^{T}$, there is

$\small{\begin{split}&{\bm{\mathbf{h}}}_{k}=||{\bm{\mathbf{y}}}_{k}-{\bm{\mathbf{x}}}_{k}||=\left[\begin{array}[]{cccc}||{\bm{\mathbf{y}}}_{k}(1)-{\bm{\mathbf{x}}}_{k}(1)||&...\\\
...&||{\bm{\mathbf{y}}}_{k}(n)-{\bm{\mathbf{x}}}_{k}(n)||\end{array}\right]\rightarrow{\bm{\mathbf{h}}}^{*}=||{\bm{\mathbf{y}}}_{0}-{\bm{\mathbf{x}}}_{0}||,\end{split}}$

hence it implies that
${\bm{\mathbf{y}}}_{0}={\bm{\mathbf{x}}}_{0}+{\bm{\mathbf{h}}}^{*}{\bm{\mathbf{d}}}_{0}$.

Now, to show ${\bm{\mathbf{y}}}_{0}$ minimizes the $f$ along the
${\bm{\mathbf{x}}}_{0}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}}_{0}$, we know for
each k and ${\bm{\mathbf{h}}}$, $0<{\bm{\mathbf{h}}}<\infty$, there is

$f({\bm{\mathbf{y}}}_{k})\leq
f({\bm{\mathbf{x}}}_{k}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}}_{k}).$

Therefore, by continuity of $f$, $k\rightarrow\infty$ leads to
$f({\bm{\mathbf{y}}}_{0})\leq
f({\bm{\mathbf{x}}}_{0}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}}_{0})$ for all
${\bm{\mathbf{h}}}$ in which shows

$f({\bm{\mathbf{y}}}_{k})\leq\min{\;f({\bm{\mathbf{x}}}_{k}+{\bm{\mathbf{h}}}{\bm{\mathbf{d}}}_{k})}.$

This proves the definition of
${\bm{\mathbf{y}}}_{0}\in\Gamma({\bm{\mathbf{x}}}_{0},{\bm{\mathbf{d}}}_{0})$
for having a closed set.