A Novel Model for 3D Motion Planning for a Generalized Dubins Vehicle with Pitch and Yaw Rate Constraints
^††thanks: ¹ Deepak Prakash Kumar is a Postdoctoral Scholar in the Center of Resilient Autonomous Systems at University of California, Irvine, CA 92697, USA (e-mail: deepakprakash1997@gmail.com).
² Swaroop Darbha is with the Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843, USA (e-mail: dswaroop@tamu.edu).
³ Satyanarayana Gupta Manyam is with DCS Corporation, Dayton, OH 45431 USA (e-mail: msngupta@gmail.com)
⁴ David Casbeer is with the Control Science Center, Air Force Research Laboratory, Wright-Patterson Air Force Base, OH 45433 USA (e-mail: david.casbeer@us.af.mil).
DISTRIBUTION STATEMENT A. Approved for public release. Distribution is unlimited. AFRL-2025-3035; Cleared 06/17/2025.

The use of Unmanned Aerial Vehicles (UAVs) is rapidly growing in civilian and military applications, including search and rescue and surveillance. Fixed-wing UAVs are of particular interest due to longer flight times, larger payload capacity, and the ability to fly at higher altitudes[1, 2]. However, they are persistently in motion, i.e., cannot stop or hover mid-air, and cannot change their heading angle instantaneously. Hence, they have a bound on the rate of change of their heading/orientation, which manifests itself as curvature constraints on the path. Motion planning is important for these vehicles, in which the goal is to plan the optimal path to travel from one configuration (i.e., position and orientation together) to another. The objective of interest in this paper is to obtain the minimum-time (or distance) path(s). We seek a finite set of candidate paths that includes the optimal path for any boundary condition. These candidate paths are suitable for constructing paths for fixed-wing aircraft or yaw rate-constrained vehicles.

Motion planning for yaw rate-constrained vehicles is typically addressed by considering a simplified kinematic model, called the Dubins model. This models a vehicle traveling at a constant speed and has a minimum turning radius constraint, which is suitable for UAVs traveling at constant altitude (in 2D). Dubins [3] solved the problem of the shortest path between a pair of configurations on a plane. It was shown that the optimal path is of type $CSC,CCC,$ or a degenerate path of the same, where $C=L,R$ denotes a left or a right turn of minimum turning radius, and $S$ denotes a straight line segment. Although Dubins showed this result using geometric techniques, the same result was later derived using Pontryagin’s Minimum Principle [4] in [5] using simpler proofs. Various variants of the planar path planning problem have been explored with variations in the model and/or the objective, such as in [6], where different left and right turning radius was considered.

Motion planning for such vehicles in 3D has also been an area of interest, where the shortest path to travel from one configuration to another, considering their motion constraints, is sought. The 3D problem applies not only to fixed-wing UAVs but also to underwater gliders and robots [7, 8]. Although specifying the heading angle alone uniquely defines the orientation of the vehicle in the 2D problem, the 3D problem requires both the heading angle and the plane in which the UAV lies. The UAV’s plane can be uniquely described using two additional angles (pitch and roll) or by a vector along its lateral or normal direction.

Simple kinematic models have also been used to tackle the 3D problem, where, similar to the 2D problem, the generated path can be tracked using a lower-level controller [9]. Because these paths can be computed very quickly, they can be combined with algorithms such as Rapidly Exploring Random Tree (RRT) to find feasible, obstacle-avoiding routes [10]. This method has also been experimentally demonstrated using an ARDrone in [11].

To our knowledge, the first exploration of the 3D problem was by Sussman [12]. The author showed that the optimal path is of type $CSC,$ $CCC,$ or a degenerate¹¹1Degenerate paths of $CSC$ and $CCC$ paths are $CS,$ $SC,$ $C,$ and $S$. version of these paths, or a helicoidal arc to connect a given location and heading direction²²2We note that in 3D, specifying only the location and heading direction does not fully determine the vehicle’s configuration, since the orientation of the plane containing the vehicle is not uniquely defined. In contrast, for 2D motion planning, location and heading are sufficient to uniquely specify the configuration.. Unlike the 2D problem, there may exist infinitely many $CSC$ paths between the initial and final configurations since the plane containing the $C$ segments can be freely chosen. Hence, efficient construction of $CSC$ paths has been explored in the literature. In [13], a $CSC$ path was constructed for instances where the initial and final locations are spaced sufficiently far apart using geometric and numerical approaches. In a later work [14], the authors adopted this path as an initial guess for a nonlinear optimization problem that was solved using a multiple-shooting method to improve the solution. The $CSC$ path construction was addressed in [15] as an inverse kinematics problem for a five degrees-of-freedom robotic manipulator. Analytical solutions were obtained for the path parameters to improve computational efficiency.

In [16] the 3D problem was addressed with a model that has two controllable inputs - one for yaw rate, and another for the rate of change of the altitude. In this work, a “Dubins airplane model” was proposed, and it was shown that the optimal trajectory comprises segments with arcs of minimum turning radius, straight lines, or a Dubins path of a certain length. Depending on the altitude difference between the initial and final locations, feasible solutions were generated by introducing additional segments to attain the desired altitude. This model was modified in [17], wherein a pitch angle constraint was introduced. Contrary to [16], this paper generates trajectories by incorporating helicoidal segments to attain the desired altitude when necessary. Additionally, the constructed paths were validated by simulations using a six-degree-of-freedom model (given in [18]) and a vector-field-based guidance law, inspired by [19], for tracking.

The 3D motion planning problem has also been addressed using the 2D Dubins result in [20]. The authors consider a total curvature constraint and a bounded pitch angle for the vehicle. They decouple the $CSC$ path to connect the given initial and final locations and heading vector into horizontal and vertical components. In this regard, they use the horizontal projection of the configurations to connect the locations and heading angles. The obtained path length is utilized as the $x$ coordinate in the vertical plane, with the desired altitude difference to be attained serving as the $z$ coordinate. Using iterative optimization of the horizontal and vertical turning radii, a feasible solution is constructed such that the total curvature bound and pitch angle bounds are satisfied. The authors used this solution to provide an initial guess to a non-linear optimization problem to further refine the path in [21].

In the existing literature, we observe that the complete orientation of the vehicle has not been considered for 3D motion planning. This is because the description of the pitch and heading angles alone does not uniquely describe its orientation. For example, Fig. 1 shows two orientations for the same vehicle moving along a straight line trajectory with a pitch angle of $0^{\circ}$ and a heading angle of $45^{\circ}:$ one where the roll angle is $0^{\circ}$ and another where the roll angle is $45^{\circ}.$ Infinitely many orientations exist for the same trajectory. However, from these figures, we can observe that prescribing the longitudinal direction and the lateral or normal direction of the vehicle would uniquely describe the orientation.³³3We remark here that in Fig. 1, unit vectors along the longitudinal, lateral, and normal directions are referred to as tangent, tangent normal, and surface normal vectors, respectively. The latter notation would be used later in the paper to describe the rotation minimizing frame model.  While the study in [22] models the complete configuration of a general robot, it is unclear how their model relates to the kinematics of a fixed-wing UAV.

The literature on path planning for aerial robots has primarily focused on models with a single control input, such as yaw rate, or a single constraint on the path’s curvature. In the 3D generalization of path planning, a single control input alone cannot capture the range of motions. This is because there are two elementary motions of interest for the motion planning of aerial vehicles: pitch and yaw. For fixed-wing UAVs, pitch motion is achieved using elevators, and yaw motion is achieved with rudder and/or aileron⁴⁴4Ailerons control the roll of the UAV, and this roll motion leads to a corresponding yaw motion for the vehicle. [18]. Since yaw and pitch motion are controlled by separate actuators, considering a single control is not sufficient. Hence, it is crucial to consider two control inputs: a bounded pitch rate and a bounded yaw rate. These constraints correspond to the minimum turning radii, $R_{pitch}$ and $R_{yaw}$, of the path’s curvature. The bounds on the pitch and yaw rates manifest as locally inaccessible spherical regions of radii $R_{pitch}$ and $R_{yaw}$, respectively. ⁵⁵5Imagine a vehicle moving in 3D space with a maximum allowable yaw rate (i.e., how quickly it can turn left or right). Due to this constraint, the vehicle cannot immediately change direction; it needs a certain minimum turning radius to execute a yaw maneuver. A yaw motion sphere is a spherical region around the vehicle that it cannot enter directly unless it has traveled a sufficient distance to allow for the required turning maneuver. This is analogous to the turning circles on the left and right in the 2D Dubins path problem. Similar to the yaw motion spheres that define inaccessible regions in the horizontal plane, pitch constraints define vertical maneuverability limits. Though two control inputs are considered in [16], the second control input is the rate of change of the altitude of the vehicle; furthermore, the pitch angle is not considered in this model, which make it tmore appropriate for a quadcopter.

A demonstration of the limitations of state-of-the-art approaches that rely on a single control input is presented in Fig. 2 for two instances. In both of these cases, paths were generated using the method from [20] - for which the code is publicly available. The minimum turning radius for this model was set to 40 meters to enforce the curvature constraint. In the following, we analyze these examples to explain why this path planning method is either infeasible or inefficient for the proposed (3D) model.

1. 1.

   The proposed model incorporates two independent control inputs for the yaw and pitch rates. We set the minimum pitch turning radius, $R_{pitch}$, to 40 meters⁶⁶6In the paper, we will interchangeably use “m” to denote “meters”. for the first example. This choice is consistent with the parameter used in [20]. Since the yaw turning radius, $R_{yaw}$, can be chosen independently, we set it to 50 meters. The solution generated by [20], shown in Fig. 2a, violates the yaw rate constraint by entering the yaw motion sphere. Such a maneuver is infeasible as the vehicle must travel a sufficient distance outside the sphere before entering it. This behavior is analogous to the infeasibility of entering the left or right turning circles in the 2D Dubins problem [3].

2. 2.

   In the second example, we alternately pick $R_{yaw}$ in our model to be the same as the parameter in [20], which is $40$ meters. Since $R_{pitch}$ is free to choose, we set $R_{pitch}$ to be equal to $60$ meters. The path obtained using the algorithm from [20] is shown in Fig. 2b. We observe that the path enters one of the pitch motion spheres (which lies at the top of the vehicle), and hence violates the pitch rate constraint.

Alternatively, one could argue that the maximum of $R_{yaw}$ and $R_{pitch}$ can be chosen as the minimum turning radius in [20]. However, this would lead to inefficient motion planning, since the vehicle would take larger-than-necessary turns in some instances.

The presented issues were addressed by the authors in [23], where a special case of motion planning on the surface of a sphere was studied. This paper provided insights for the 3D motion planning by considering a vehicle model with bounded yaw rate and pitch rate. Furthermore, the spherical motion planning problem was shown to arise as an intermediary problem to be solved for the 3D problem.

Having identified two major issues in the literature, the contributions of this paper are as follows:

1. 1.

   We present a novel model using a rotation minimizing frame, also called the Bishop frame, to obtain the shortest path for a vehicle subject to pitch rate and yaw rate constraints. We build on the insights provided in [23] for the 3D problem.

2. 2.

   We prove that the pitch rate and yaw rate constraints manifest as four spheres around the vehicle that represent temporarily inaccessible regions, thereby appropriately generalizing the 2D Dubins model to 3D.

3. 3.

   We propose a path construction algorithm that consists of three classes of paths. The main idea is to build path segments on spherical surfaces that are tangent to the initial and final configurations, and these segments are connected by an intermediary surface. This surface could be a cylindrical envelope, a cross-tangent plane, or another spherical surface.

4. 4.

   We pose and solve a Dubins-type path planning problem subject to curvature constraints on a cylindrical surface. To the best of our knowledge, the cylindrical motion planning problem has not been addressed in the literature. The proposed solution method involves unwrapping the surface to a two-dimensional plane, computing the optimal Dubins path, and then mapping back onto the cylindrical surface.

5. 5.

   We present extensive numerical results on several instances. We show the effect of $(i)$ model that defines the complete configuration of the vehicle, (i.e., heading and lateral orientation and $(ii)$ the impact of minimum turning radii on the best feasible path. We also observe that our algorithm can produce a high-quality feasible solution within 10 seconds. Additionally, we provide the code in a publicly available repository.⁷⁷7The code is available at https://github.com/DeepakPrakashKumar/3D-Motion-Planning-for-Generalized-Dubins-with-Pitch-Yaw-constraints.

Let $t$ and $s$ denote the time and arc length, respectively, and $\mathbf{X}(s)$ denote the instantaneous location of the vehicle. We consider a Rotation-Minimizing frame, also called a Bishop frame [24], attached to the center of mass of the UAV. Let $\mathbf{T},\mathbf{Y},\mathbf{U}$ denote the unit vectors of the Bishop frame with $\mathbf{T},\mathbf{Y}$ directed along the longitudinal and lateral directions of the vehicle, respectively. The vector $\mathbf{U}:=\mathbf{T}\times\mathbf{Y}$ is along the normal direction of the vehicle. Fig. 3 shows the vehicle configuration with the vectors $\mathbf{T},\mathbf{Y}$ and $\mathbf{U}$.

The instantaneous angular velocity of the frame can be written as

$\displaystyle\omega(t)=\omega_{x}(t){\bf T}(t)+\omega_{y}(t){\bf Y}(t)+\omega_{z}(t){\bf U}(t),$

with $\omega_{x},\omega_{y},\omega_{z}$ denoting the components in the ${\bf T},{\bf Y}$, and ${\bf U}$ directions, respectively. One may think of them as the roll, pitch, and yaw rates of the body, respectively. Notably, the Rotation Minimizing Frame (RMF) is constructed so that $\omega_{x}(t)=0$: there is no rotation about the tangent, minimizing frame twisting (essentially, the roll rate is set to zero). The kinematics of the frame then satisfy:

	$\displaystyle\frac{d{\bf T}}{dt}$	$\displaystyle=\omega(t)\times{\bf T}(t)=\omega_{z}(t){\bf Y}(t)-\omega_{y}(t){\bf U}(t),$
	$\displaystyle\frac{d{\bf Y}}{dt}$	$\displaystyle=\omega(t)\times{\bf Y}(t)=-\omega_{z}(t){\bf T}(t),$
	$\displaystyle\frac{d{\bf U}}{dt}$	$\displaystyle=\omega(t)\times{\bf U}(t)=\omega_{y}(t){\bf T}(t).$

A key property of an RMF is that ${\bf Y}(t)$ and ${\bf U}(t)$ change only in the direction of ${\bf T}(t)$.

Assume that the vehicle moves at a constant, nonzero longitudinal speed $V_{0}$. Defining

$\displaystyle\kappa_{g}(s):=\frac{\omega_{z}(s)}{V_{0}},\qquad\kappa_{n}(s):=\frac{\omega_{y}(s)}{V_{0}},$

(1)

the kinematic equations parameterized in terms of $s$ become

	$\displaystyle\frac{d\mathbf{X}(s)}{ds}$	$\displaystyle=\frac{1}{V_{0}}\frac{d{\bf X}}{dt}=\mathbf{T}(s),\quad$		(2)
	$\displaystyle\frac{d\mathbf{T}(s)}{ds}$	$\displaystyle=\kappa_{g}(s)\mathbf{Y}(s)+\kappa_{n}(s)\mathbf{U}(s),$		(3)
	$\displaystyle\frac{d\mathbf{Y}(s)}{ds}$	$\displaystyle=-\kappa_{g}(s)\mathbf{T}(s),\quad$		(4)
	$\displaystyle\frac{d\mathbf{U}(s)}{ds}$	$\displaystyle=-\kappa_{n}(s)\mathbf{T}(s).$		(5)

The bounds for control inputs $\kappa_{g}$ and $\kappa_{n}$, which we refer to as geodesic curvature and normal curvature, respectively, are stated as shown below,⁸⁸8We use the notation $\kappa_{g}$ and $\kappa_{n}$, since they have the same form as that of the geodesic and normal curvatures in the Darboux frame model, a differential geometric model.

$\displaystyle|\kappa_{n}|\leq\frac{1}{R_{pitch}},\quad|\kappa_{g}|\leq\frac{1}{R_{yaw}}.$

(6)

We shall later show that $R_{pitch}$ is the minimum turning radius corresponding to pitch motion and $R_{yaw}$ is the minimum turning radius corresponding to the yaw motion. The objective is to compute the minimum distance trajectory from the initial to final configuration, defined by $\mathbf{X},$ $\mathbf{T},$ $\mathbf{Y},$ and $\mathbf{U}$. Hence, the cost to minimize is $J=\int 1ds.$ Geometrically, the path planning problem can be depicted as shown in Fig. 4. A detailed description of this figure follows.

In Fig. 4, it can be observed that four spheres are constructed, which are tangential to the vehicle configuration. To understand how they appear, we need to understand the geometric impact of the curvatures, $\kappa_{n}$ and $\kappa_{g}.$ To this end, we derive the closed-form expression for $\mathbf{X},$ $\mathbf{T},$ $\mathbf{Y},$ and $\mathbf{U}$ over an interval wherein $\kappa_{g}$ and $\kappa_{n}$ are constants. The obtained expressions are shown in Appendix -A.

We remark here that since the control inputs appear linearly in the differential equations (3)-(5), and a minimum time problem is considered, the optimal control actions are expected to be bang-bang from Pontryagin’s minimum principle [4]. Therefore, it is sufficient to consider intervals in which $\kappa_{g}$ and $\kappa_{n}$ are constant. Furthermore, it suffices to consider $\kappa_{g}\in\{-\frac{1}{R_{yaw}},0,\frac{1}{R_{yaw}}\}$ and $\kappa_{n}\in\{-\frac{1}{R_{pitch}},0,\frac{1}{R_{pitch}}\}$.

Using the closed-form expressions derived in Appendix -A, we state and prove the following two Lemmas.

The following lemma states a similar result in a different axis, and the proof follows the same reasoning.

From these two lemmas, we see that the normal curvature $\kappa_{n}$ governs the pitch motion, while the geodesic curvature $\kappa_{g}$ governs the yaw motion. In fact, these curvatures directly correspond to the vehicle’s pitch rate and yaw rate, respectively (which is expected based on the Bishop frame setup and the definitions in (1)). When the vehicle moves with its maximum pitch rate and zero yaw rate, it follows a great circle of radius $R_{pitch}$ on the orange or purple sphere shown in Fig. 4. This result comes from Lemma 1. Since the vehicle travels at unit speed, the time to complete the circle is $t_{pitch}=2\pi R_{pitch}.$ Over this time, the pitch angle changes by $2\pi$, so the pitch rate is $\frac{2\pi}{t_{pitch}}=\frac{1}{R_{pitch}}.$ A similar argument holds for the yaw rate, giving a maximum value of $\frac{1}{R_{yaw}}$. Therefore, $\kappa_{n}$ and $\kappa_{g}$ represent the vehicle’s pitch and yaw rates, respectively.

By varying $\kappa_{n}\in\{-\frac{1}{R_{pitch}},0,\frac{1}{R_{pitch}}\}$ and $\kappa_{g}\in\{-\frac{1}{R_{yaw}},0,\frac{1}{R_{yaw}}\}$, we obtain nine distinct motion primitives, shown in Fig. 5. These were generated using the closed-form expressions, presented in Appendix -A. Using Lemma 1, we find that the segments $L_{si},$ $R_{si},$ $L_{so},$ and $R_{so}$ have radius $\frac{1}{\sqrt{\frac{1}{R_{yaw}^{2}}+\frac{1}{R_{pitch}^{2}}}}$, corresponding to motion with maximum absolute pitch and yaw rates. Here, $L$ and $R$ denote a left turn and right turn, respectively, which correspond to $\kappa_{g}=\frac{1}{R_{yaw}}$ and $\kappa_{g}=-\frac{1}{R_{yaw}},$ respectively. Additionally, subscripts “$si$” and “$so$” are used to refer to the segments that lie on the “inner” sphere and “outer” sphere, respectively; the inner sphere corresponds to $\kappa_{n}=\frac{1}{R_{pitch}},$ and the outer sphere corresponds to $\kappa_{n}=-\frac{1}{R_{pitch}}.$ The segments $G_{si}$ and $G_{so}$ result from pure pitch motion ($\kappa_{g}=0$), while $L_{p}$ and $R_{p}$ result from pure yaw motion. When both curvatures are zero ($\kappa_{n}=\kappa_{g}=0$), the vehicle moves in a straight line segment $S$.

Using the obtained motion primitives and the observation that $\kappa_{n}$ and $\kappa_{g}$ attaining values of $\pm\frac{1}{R_{pitch}}$ and $\pm\frac{1}{R_{yaw}}$ yields two spheres each (a pair along $\mathbf{U}$ and a pair along $\mathbf{Y}$), we can observe that at both the initial and final configuration, four spheres exist around the vehicle. Additionally, portions of the optimal path will lie on one of the four spheres at the initial configuration and one of the four spheres at the final configuration⁹⁹9The only motion primitive that does not lie on a sphere is a straight line segment $S.$ However, even in this case, a portion of the optimal path can be modeled to lie on one of the four spheres; however, the path will be trivial, i.e., of zero length.. Hence, we propose three classes of paths to construct a feasible path connecting one of the initial spheres with one of the final spheres. We construct the path using three types of intermediary surfaces (or classes): a cylindrical envelope, a planar surface, or a spherical surface. These three classes of paths are a generalization of the classical $CSC$ and $CCC$ paths for the planar Dubins problems. In our algorithm, we consider a sphere at the initial or final configuration to serve as a generalization of the turn segment ($C$) in a plane; furthermore, we consider the cylindrical envelope and planar surface to generalize the $S$ segment. In the following section, we describe the three classes of paths in more detail.

Remark: In general, the yaw and pitch rates for different UAVs may vary and can be coupled. The problem posed here is still of significant interest in obtaining lower and upper bounds for the shortest path length. One such case is illustrated by the region within the boundaries, shown in brown, in Fig. 6. However, by replacing the boundary with a rectangular region inscribed within this area, one can derive an upper bound that is a feasible solution. Similarly, by outer approximating the allowable region with a larger rectangular region, shown in green in Fig. 6, a lower bound for the optimal path length can be obtained.

The constructed paths start and end on a spherical surface at the initial and final configurations. Three distinct classes of paths are presented, each using a different intermediary surface for the sub-path between the spheres, which can be cylindrical, planar, or spherical. We refer to the spheres centered along the $\mathbf{U}$-axis as the inner (orange, along $\mathbf{U}$) and outer (purple, along $-\mathbf{U}$) spheres, as shown in Fig. 4. Similarly, the spheres centered along the $\mathbf{Y}$-axis are called the left (green, along $\mathbf{Y}$) and right (blue, along $-\mathbf{Y}$) spheres.

The intermediary sub-paths considered are as follows:

1. 1.

   In the first class, the sub-path is constructed using a cylindrical envelope. In this case, we connect the pair of spheres of the same type at the initial and final configurations. There are four such pairs: inner-to-inner, outer-to-outer, left-to-left, and right-to-right. Fig. 7a illustrates the cylindrical envelope between inner and outer spheres. We will later show that these paths satisfy the curvature constraints (6). The full construction is detailed in Section IV.

2. 2.

   In the second class of paths, a sub-path between a pair of spheres is constructed on a cross-tangent plane. For spheres of opposite type, such as inner-to-outer, outer-to-inner, left-to-right, or right-to-left, a path through a cylindrical envelope is not feasible. This is because the normal vector of the cylindrical surface, $\mathbf{U}$ for pitch spheres and $\mathbf{Y}$ for yaw spheres, remains constant along the envelope and does not support a continuous feasible orientation between opposing directions. An example for the inner-to-outer case is shown in Fig. 7b. Details of this construction are provided in Section V.

3. 3.

   In the third class of paths, we construct sub-paths between pairs of spheres of the same type using an intermediary sphere. There are four such configurations: inner–outer–inner, outer–inner–outer, left–right–left, and right–left–right. These paths are designed for initial and final locations that are close to each other. An example of this type is illustrated in Fig. 7c, and the full construction is described in Section VI.

To generate a feasible path connecting two spheres of the same radius and type via a cylindrical envelope (as shown in Fig. 7a), the vehicle follows this sequence:

• •

  Step 1: The path starts on the initial sphere and transitions to the cylindrical surface. The transition point, $\mathbf{X}_{ic}$, lies on the boundary circle formed by the intersection of the sphere and the cylindrical envelope. At this point, its longitudinal direction aligns with $\mathbf{T}_{ic},$ as illustrated in Fig. 8.

• •

  Step 2: The path exits the cylindrical envelope at $\mathbf{X}_{oc}$, with longitudinal direction $\mathbf{T}_{oc}$.

• •

  Step 3: Finally, the path continues on the final sphere to reach the desired final configuration.

Note that the entry and exit points on the cylinder, along with their corresponding tangent directions, are not fixed and can be freely chosen. We parameterize these directions using four angles: $\theta_{ic},$ $\phi_{ic},$ $\theta_{oc},$ and $\phi_{oc}.$ The following section describes how the path on the spheres and the cylinder is constructed based on these four parameters.

In this section, we describe the second class of paths where the sub-path between the initial and final spheres is constructed on a cross-tangent plane. There exist infinitely many cross-tangent planes between these two spheres; the locus of the point of intersection of these cross-tangent planes with the initial/final sphere will be a circle, as shown in Fig. 7b. To uniquely define a cross-tangent plane, we use angle $\theta$ as a parameter (see Fig. 13).

Remark: From Fig. 13, we can observe that the considered cross-tangent plane exists when $\|\mathbf{r}_{f}-\mathbf{r}_{i}\|_{2}\geq 2\overline{R},$ i.e., when the spheres at the initial and final configurations do not intersect.

We denote the center of the locus at the initial and final spheres by $\mathbf{A}$ and $\mathbf{B}$, respectively, as shown in Fig. 13. The location of $\mathbf{A}$ and $\mathbf{B}$ can be obtained as

	$\displaystyle\mathbf{A}$	$\displaystyle=\mathbf{r}_{i}+\overline{R}\cos{\alpha}\left(\frac{\mathbf{r}_{f}-\mathbf{r}_{i}}{\|\mathbf{r}_{f}-\mathbf{r}_{i}\|_{2}}\right),$		(21)
	$\displaystyle\mathbf{B}$	$\displaystyle=\mathbf{r}_{f}-\overline{R}\cos{\alpha}\left(\frac{\mathbf{r}_{f}-\mathbf{r}_{i}}{\|\mathbf{r}_{f}-\mathbf{r}_{i}\|_{2}}\right),$		(22)

where $\alpha:=\cos^{-1}\left(\frac{2\overline{R}}{\|\mathbf{r}_{f}-\mathbf{r}_{i}\|_{2}}\right)$. Here, the expressions for $\mathbf{r}_{i}$ and $\mathbf{r}_{f}$ are given in (7) and (8), respectively, and $\overline{R}$ is given in (11). The cross-tangent plane between the initial and final spheres is needed for inner-outer, outer-inner, left-right, and right-left pairings. It follows that $\delta^{initial}_{i,o}=-\delta^{final}_{i,o}$ and $\delta^{initial}_{l,r}=-\delta^{final}_{l,r}.$

To define the parameter $\theta$, we first designate a unit vector $\mathbf{x}$ perpendicular to $\mathbf{r}_{f}-\mathbf{r}_{i}$, as shown in Fig. 13.¹⁵¹⁵15The generation of $\mathbf{x}$ is similar to the procedure described for the cylindrical envelope. The angle $\theta$ specifies the point of tangency on the initial sphere, with respect to $\mathbf{x}$. We then define another unit vector $\mathbf{y}:=\left(\frac{\mathbf{r}_{f}-\mathbf{r}_{i}}{\|\mathbf{r}_{f}-\mathbf{r}_{i}\|_{2}}\right)\times\mathbf{x}$, which is orthogonal to both $\mathbf{x}$ and the axis ($\mathbf{k}$, defined in (9)) between the spheres. Using $\mathbf{x}$ and $\mathbf{y}$, we now express the entry and exit points $\mathbf{X}_{ic}$ and $\mathbf{X}_{oc}$, which are the points of tangency (see Fig. 13).

	$\displaystyle\mathbf{X}_{ic}(\theta)$	$\displaystyle=\mathbf{A}+\overline{R}\sin{\alpha}\cos{\theta}\mathbf{x}+\overline{R}\sin{\alpha}\sin{\theta}\mathbf{y},$
	$\displaystyle\mathbf{X}_{oc}(\theta)$	$\displaystyle=\mathbf{B}+\overline{R}\sin{\alpha}\cos{(\theta+\pi)}\mathbf{x}+\overline{R}\sin{\alpha}\sin{(\theta+\pi)}\mathbf{y}.$

Next, we parameterize the tangent vectors at the entry and exit points using angles $\phi_{ic}$ and $\phi_{oc}$, defined relative to the axis $\mathbf{t}$, where $\mathbf{t}$ is a unit vector pointing from $\mathbf{X}_{ic}(\theta)$ to $\mathbf{X}_{oc}(\theta)$ (see Fig. 13). The tangent vectors $\mathbf{T}_{ic}$ and $\mathbf{T}_{oc}$ at the entry and exit points are derived as below:

	$\displaystyle\mathbf{T}_{ic}(\phi_{ic})$	$\displaystyle=\cos{\phi_{ic}}\mathbf{t}(\theta)+\sin{\phi_{ic}}\left(\frac{\mathbf{X}_{ic}(\theta)-\mathbf{r}_{i}}{\overline{R}}\times\mathbf{t}(\theta)\right),$
	$\displaystyle\mathbf{T}_{oc}(\phi_{oc})$	$\displaystyle=\cos{\phi_{oc}}\mathbf{t}(\theta)+\sin{\phi_{oc}}\left(\frac{\mathbf{X}_{ic}(\theta)-\mathbf{r}_{i}}{\overline{R}}\times\mathbf{t}(\theta)\right).$

Given the location and tangent vectors at the exit point from the initial sphere and the entry point of the final sphere, we can construct the optimal path over each sphere using the approach in Section IV-C. We can also construct the path on the cross-tangent plane using the 2D Dubins result [3, 28], illustrated in Fig. 14. In this figure, the minimum turning radius $R_{plane}$ is dictated by the type of spheres considered at the initial and final configurations. If we are considering inner-outer or outer-inner connections, $R_{plane}=R_{yaw}$ since the vehicle moves in the $\mathbf{T}-\mathbf{Y}$ plane (as observed from Fig. 13). Alternatively, $R_{plane}=R_{pitch}$ if left-right or right-left sphere connections are considered.