How a serpentine tail assists agile motions of kangaroo rats: a dynamics and control approach

Kangaroo rat is a good representative for general bipedalism with a serpentine tail. Modeling and analyzing the kangaroo rat motion helps to understand the serpentine tail functionalities in agile motions of bipedal mobile platforms, and this understanding is expected to lay the foundation for the future development of such robotic systems. This paper analyzes the kangaroo rat motions through dynamic modeling and control. The system dynamic model is established using the inertia matrix method, and two typical serpentine tail models are considered: a continuum tail model where the tail is modeled as several constant curvature arcs, and an articulated tail model where the tail is discretized into rigid links. Regularized contact model is used to compute the ground reaction force (GRF). To automatically plan the tail motion, numerical optimal control techniques (i.e., direct collocation method) are utilized. Partial feedback linearization is then used to track the designed tail trajectory. Based on the formulated dynamic model and motion controller, two representative tail functions (airborne righting and supporting) were simulated and analyzed. The results validated the proposed modeling and control framework and showed the nontrivial functionalities of the serpentine tail in helping the kangaroo rat to achieve agile motions. Moreover, comparative studies on the two tail models and the tail segmentations were performed to analyze the model differences. The results demonstrated that the articulated tail model is a good approximation of the continuum tail model, and more tail segments and links enhance the kangaroo rat’s ability to deliberately adjust its motion.

From a mechanics perspective, the tail functionalities could be summarized into two categories: static and dynamic functions. The static functions of the tail mean that the tail is only used to inject static forces into the system. The system motion is governed by statics, and only position information is used. Such functions include adjusting the system center of mass (COM), supporting the body, grasping, manipulation, etc. The dynamic functions of the tail mean that the tail is used to inject inertial forces into the system. The system motion is governed by the full-level dynamics, and kinematic information including acceleration is necessary. Such functions include airborne righting, maneuvering, stabilization, etc.
Kangaroo rat is an animal that relies heavily on its tail to assist its agile motions. For instance, Fig. 1 shows that a kangaroo rat flips over its body in the air using its tail after being attacked by a rattlesnake [29]. Further observations imply that kangaroo rats may use most of the functionalities mentioned above, that is, using the tail as an additional appendage (e.g., leg) to support its locomotion on the ground, using the tail in the air to change its torso orientation, and using the tail to adjust the ground reaction force (GRF) during ground contact events (indirectly affecting the locomotion). Therefore, the kangaroo rat is thought to be a great benchmark model to investigate the functions of a tail for bipedal locomotion. This idea recently started attracting attention and a series of research efforts were carried out, from both the biological side [30,31] and the robotic side [10,32,33]. However, the existing efforts on the biological side mainly focus on data interpretation using statistical analysis, and the efforts on the robotic side mainly focus on using the pendulum tail abstraction. Both approaches use a simplified dynamic model, which is either in 2D or using a single-segment singlelink tail model. This brings up an important missing part for this research-a high-fidelity 3D dynamic model that captures the serpentine nature of the tail. Therefore, in this paper, we aim to investigate how the kangaroo rat uses its serpentine tail to achieve agile motions from a dynamics and control perspective, that is, establishing a 3D dynamic model of the kangaroo rat with a serpentine tail and analyzing the tail's functionalities with helping the kangaroo rat to achieve agile motions. This investigation is expected to provide a solid understanding (mechanically and quantitatively) of the kangaroo rat motions and what role the serpentine tail plays in achieving these motions. These understandings could further benefit Fig. 1 A kangaroo rat flips itself in the air by swinging its serpentine tail. Snapshots are captured from the supplementary video 1 in [30], which is also available online at [29] the development and control of such robotic systems, for instance, developing novel tailed biped robots for executing highly agile motions (military technology) or developing tail-assisted fast maneuvering techniques for space robots (space technology). In control theory, solving the unique kangaroo rat tail-assisted locomotion problem could further enrich the existing whole-body control theory.
The major contributions of this work are summarized as follows. Firstly, a general 3D bipedal model with multi-contact feet and a general serpentine tail is established. Secondly, two representative and general serpentine tail structures, namely the articulated tail model and the continuum tail model, are formulated. Thirdly, numerical optimal control methods (direct collocation) are used to automatically plan the tail motion, which allows incorporating advanced input and state constraints. Fourthly, through the established dynamic model and controller, agile motions, tail functionalities, and tail model comparative studies of the kangaroo rat are analyzed. Each of these four points is proposed in the literature for the first time. The combination of these four points provides a highfidelity kangaroo rat model and explains the serpentine tail's functionalities with assisting the kangaroo rat to achieve agile motions from a mechanics perspective for the first time.
The overview of the modeling and control approaches in this work is shown in Fig. 2 where the entire process consists of three parts: dynamic modeling, motion planning, and simulation. The dynamic modeling part establishes the equation of motion (EOM) of a kangaroo rat model considering two types of serpentine tail models and the ground contact events. This model is used in both the motion planning part and the simulation part. The difference is that in motion planning, the dynamic model is used as path constraints for the tail trajectory optimization tasks while in simulation, the dynamic model is numerically solved to verify the motion planning effectiveness. To convert the designed trajectory to actual joint inputs, a feedback-based trajectory tracking controller was implemented in the simulation.
The following sections of this paper are organized as follows. Section 2 formulates the dynamic model of the kangaroo rat. Section 3 derives the motion planning algorithm and the trajectory tracking controller for the serpentine tail motion control problem. Section 4 implements the numerical experiments and explores the serpentine tail functionalities in kangaroo rat locomotion. Section 5 summarizes the main results of this paper.

Dynamic modelling
The dynamic model of the kangaroo rat is shown in Fig. 3 where the body is abstracted into an ellipsoid  Fig. 3 Kangaroo rat kinematic model. Note that the animal on the top is for morphological reference only. It is a jerboa-a different species with similar morphology and behaviors and the leg is modeled as a planar three-parallelogram mechanism. The foot is modeled with a flat sole such that it allows stable standing. The inertial frame RS :¼ ðS;x s ;y s ;z s Þ is attached on the ground with z s pointing vertically up and y s pointing forward. The body-fixed frame of the torso RB :¼ ðB;x b ;y b ;z b Þ is attached at the torso COM with its basis vector aligning with the principal axes of the torso moment of inertia (MOI). Note that the torso COM is not necessarily located at the ellipsoid center or its focus, but the MOI is always computed with respect to the ellipsoid center.
The EOM is established using the inertia matrix method in [34], as shown in Eq. (1). The leg inertial loading is neglected due to its small influence (lightweight legs) on the overall motion. However, its kinematics is still needed to generate the GRF.
where H is the system inertia matrix, C is the bias force containing all the gravitational, centrifugal, and Coriolis forces, s ta is the actuation input from the tail, f c is the GRF, J ta and J f are the corresponding Jacobians for the actuation and GRF, respectively, where p b is the position of point B, / b ¼ ½/ x / y / z T is the torso orientation, and q t is the generalized coordinates that describe the tail configuration. The serpentine tail consists of m serially connected segments. Each segment can bend spatially (i.e., bending to left and right or bending up and down), and is fully defined by these two mobilities. Therefore, the tail has 2m generalized coordinates and the system has d ¼ 6 þ 2m DOF. The components in Eq. (1) are obtained as follows: where m b , I b , H i are the torso mass, torso MOI (measured in the inertial frame), and the inertia matrix for the i-th segment of the tail, respectively. J b;x and J b;x are the trivial Jacobians for p b and / b , respectively. g ¼ ½0 0 g T is the negative gravitational acceleration.
is the angular velocity of the torso. The tilde above a vector indicates its skewsymmetric expansion. The ''ID(tail)'' function represents the inverse dynamics of the tail subsystem that contains all the non-actuation loading (inertial, gravitational, frictional, etc.): where s i is the inertial loadings for the i-th tail segment.
The critical kinematic and dynamic information required by Eqs. (1)-(4) is presented in the following subsections in detail. To compare different serpentine tail structures, a continuum tail model and a discrete tail model are formulated.

Leg kinematics
For simplicity, the leg mechanism is set to move only in the plane that is parallel to the y b Bz b plane, i.e., there is no abduction DOF in the leg. Due to the threeparallelogram mechanism, all the hip (H), knee (K), ankle (A), metatarsophalangeal (M), and toe (T) positions in the leg plane could be described by two hip angles c l and d l , where the subscript l ¼ f1; 2g indicates the right (1) or the left (2) leg, respectively. Therefore, the metatarsophalangeal and toe positions are calculated using Eqs. (5)- (9), where the left superscript is used to indicate the measuring frame other than the inertial frame. Other subscripts are labeled to indicate a point or a line segment. For instance, p m;l stands for the position of point M l , p h2m;l means the vector from H l to M l , and L h2d is the linkage length from H l to D l .
The velocities could be calculated by differentiating the above equations directly and the corresponding Jacobians are obtained by factoring out the generalized velocity _ q. Then, the generalized GRF is obtained as

Articulated tail model
The articulated tail assumes that each tail segment consists of n serially connected links and adjacent links are connected by universal joints. All joints have the same rotation in the same segment. To compute the kinematics of this tail, a similar approach as in [35] is used. As shown in Fig. 4, the body-fixed frames are defined as follows. Frame RJ j :¼ ðJ j ;x j ;y j ;z j Þ is attached to the j-th link and is located at the j-th joint center J j . y j points along with the link and to the tail root direction. z j aligns with the yaw axis of the j-th universal joint. Note that the first joint coincides with the tail mounting point T.
With the above definitions, the orientation of each link (R j ) could be obtained recursively, as shown in Eq. (11) and (12) where R x and R z are the principal rotation matrix functions with respect to the x-axis and z-axis, respectively. a i and b i are the pitch and yaw rotation angles for the universal joints in the i-th segment, respectively. Therefore, the tail generalized coordinates for the articulated tail case is The joint positions p j;jnt and the COM position p j;com of each link could then be calculated recursively too, as shown in Eqs. (13)- (15), where L j2c and L j2j are the joint-to-COM and joint-to-joint distance, respectively. p j;com ¼ p j;jnt þ p j;j2c ð13Þ The corresponding angular velocities, linear velocities, Jacobians, and accelerations are calculated by differentiating the position relationships directly, as outlined in Appendix A. Therefore, with the above kinematic information, the segment inertia and nonactuation loading could be computed as where m at and I j;at are the mass and the MOI (measured in the inertial frame) for the j-th link, respectively. Note that Eq. (17) only accounts for the inertial and gravitational loadings. The internal joint frictions are neglected to simplify the modeling. The reasoning for this is that previous modeling and simulations [35,38] showed that the internal joint frictions mainly affect the control effort (requiring larger motor torques) and do not change the overall motion significantly. Since the focus of this paper is on the overall motions instead of the joint torques, the joint frictions were neglected for simplicity.

Continuum tail model
The continuum tail model assumes that each tail segment is a non-extensible constant curvature rod and there is no twist motion inside the rod. Mass is uniformly distributed along the rod. Due to the results in [36] that the angular kinetic energy constitutes only 5% of the total kinetic energy, the angular inertial loading for an infinitesimal slice (referring to Fig. 5) is neglected, i.e., the mass is assumed to concentrate on the neutral line. Similar to the articulated tail model, all frictions and elasticities are also neglected. Only the inertial loading and the gravitational loading are computed. Based on these assumptions, the shape of Yaw βi  Figure 5 illustrates the kinematic configuration of the continuum tail model. Note that since the tail is not extensible, h i r i ¼ L ct where the L ct is the segment length.
There are several existing methods to formulate the continuum robot dynamics, such as the classic [37] or the new [36] modal approaches. The main differences between these methods are their kinematic representations, i.e., how to describe the position and orientation of the slice x. We apply the same approach as in [38] due to its advantage to represent the kinematic integral explicitly. However, since the continuum tail model is non-extensible and non-elastic, the formulation is slightly modified. That is, the inertial matrix and the non-actuation loading of the i-th segment are computed as integrals along the neutral line: where m ct is the segment mass, x is the acceleration of O i;x , and x is the local parametric variable. Note that Eq. (19) neglects the angular inertial loading due to its small contribution to the total inertial loading. This is equivalent to concentrating the mass on the neutral line.
With the constant curvature assumption, Eq. (18) and (19) turn out to be integrable, i.e., the integral variable x could be separated from the kinematic terms. Therefore, expressing the O i;x position as an explicit function of x yields The velocity and Jacobian of O i;x are then obtained by directly differentiating Eq. (20): where J i;q , J i;a , J i;b , and J i;h are the corresponding Jacobians for p i;q , a i , b i , and h i , respectively. In Eq. (25), the Jacobian is written as a block-wise matrix multiplication (denoted by '''') form such that the integral terms are separated from the non-integral terms. This form makes the expression concise and facilitates the following integrations. Therefore, the inertia matrix for the i-th segment could be evaluated as where Adding up the gravitational acceleration and using the block-wise matrix multiplication notation, the nonactuation loading is computed as The orientation and position of each segment frame P O i could be calculated recursively, as shown in Eq. (33) and (34). The higher-order information (velocities, Jacobians, and accelerations) could be calculated by directly differentiating their position relationships.
Moreover, since the above kinematic representation is essentially a variant of the Frenet-Serret frame [37], it is not able to handle the case when the curvature equals zero, i.e., the segment becomes a straight line. For such cases (e.g., when h i \10 À6 rad), the segment dynamics are treated separately as a straight rigid link.

Ground contact model
The same ground contact model as in [39] is used, which is a regularized compliant contact model [40] that formulates the normal force as a nonlinear springdamper system and the friction force as a linear springdamper system. The friction force model obeys the Coulomb friction law, which states that the quotient of the friction force and the normal force cannot be greater than the friction coefficient. Therefore, for each contact point, the GRF is obtained as where z, K n , and D n are the penetration depth, ground stiffness, and ground damping coefficient, respectively. K f , D f , x, and l are the friction spring stiffness, friction spring damping coefficient, offset distance in x s direction away from the first contact point, and friction coefficient, respectively. f y is computed the same as the f x k k except replacing x with y. The contact detection problem (computing the contact points) for the legs and the articulated tail are straightforward, which are just the foot points M l , T l , A l , and the tail joint points J j . The contact point calculation for the continuum tail, however, requires special treatments since it changes continuously as the tail configuration changes. An easy way to find out these points is to calculate the lowest point of each tail segment, denoted as C i;ct ¼ fO i;x : x 2 0; 1 ½ minimizes the z component of p i;o;x g. Note that x here has the same sense as that in Sect. 2.3, instead of the x in Eq. (39). The necessary condition for x to minimize the z component of p i;o;x could be obtained by letting the derivative of Eq. (20) with respect to x vanish, which yields: where x is the optimal value of x, a i;z and b i;z are the z component of a i and b i , respectively. Therefore, the contact points C i;ct could be determined as For the cases that x is undetermined (e.g., when the bending plane of segment i is horizontal), x is set to 0.5. For the cases that the tail segment becomes a line segment, the contact points are manually set as its two endpoints.

Tail motion control
Controlling the motion of a 3D biped with a multisegment serpentine tail is more challenging than a 2D biped [5,10,33] or a 3D single body with a singlesegment single-link tail [7-9, 13, 32, 33] since the traditional intuition-based tail motion planning and control methods [5-20, 32, 33] are not applicable anymore. To solve this problem, we propose to use the numerical optimal control method [41] to automatically synthesize the tail motion for the kangaroo rat agile motions. The existing numerical optimal control techniques could be roughly divided into two classes: the indirect methods such as differential dynamic programming (DDP) [42], and the direct methods such as direct collocation [43] or orthogonal collocation [44]. In this paper, we chose to use the direct collocation method [43] due to its robustness for high-dimensionality systems and its flexibility to include advanced constraints (e.g., box constraints for joint motion range and joint torque limit). The motion planning will only focus on the airborne phase since the tail motion in this phase is more challenging to be planned manually. More importantly, contacts are not involved in the airborne phase and thus the system dynamics are smooth, which makes the optimization process efficient and robust. After the tail motion trajectory was generated, a partial feedback linearization (PFL) based controller [45] was implemented to track the trajectory. This motion control framework is illustrated in Fig. 6 where the active tail motion control starts from a randomly picked moment when the kangaroo rat is in the air and ends until the moment before the kangaroo rat lands on the ground. All the motions in other moments use a passive tail controller (simple damping) for which the details will be presented in Sect. 3.3.

Trajectory optimization through direct collocation
The core idea of the direct collocation technique is to treat both the system states and control inputs as decision variables for the nonlinear programming (NLP) problem [46,47]. For the motion planning problem of the kangaroo rat in the airborne phase, the objective is to design appropriate tail motion profiles q t;r ðtÞ, _ q t;r ðtÞ and torque profile s ta;r ðtÞ such that the kangaroo rat body could reach desired orientation before landing. Using the direct collocation technique, this problem could be formulated as an NLP problem: f bc q 1 ; _ q 1 ; u 1 ; q N ; _ q N ; s ta;N À Á 0; g pc q k ; _ q k ; s ta;k À Á 0; bc; pc 2 N where the decision variables are p d ¼ ft N ; q 1 ; _ q 1 ; s ta; 1; . . .; q N ; _ q N ; s ta;N g. ' 0 and ' k are the boundary objective and the path objective functions, respectively. f bc and g pc are the boundary and path constraints, respectively, including all the box constraints (e.g., joint motion range and joint torque limit) and more complicated kinematic constraints (e.g., equality constraints for the final state). Note that the dynamic constraints in Eqs. (41) and (42) do not include the GRF term (referring to Eq. (1) to see the difference) due to the airborne assumption. Moreover, Eqs. (41) and (42) implement the trapezoidal quadrature rule to enhance numerical stability. Other quadrature rules such as the Hermite-Simpson rule could be also used to increase the numerical integration accuracy with a cost of longer computation time [47].

Nontrivial optimization implementation details
Since the system has d ¼ 6 þ 2m DOF and 2m inputs (the leg inputs c l and d l are ignored since the leg motions do not affect the airborne dynamics), the NLP problem complexity could easily go up to hundreds of the decision variables and involve highly nonlinear constraints, which makes the solving process timeconsuming and susceptible to an infeasible local minimum. To ease these issues, several special optimization techniques are used. The first technique is a warm-start technique, which utilizes dynamic simulation to generate the initial trajectory guess. This way, the initial guess always satisfies the dynamics constraints in Eqs. (41) and (42), which helps the NLP to converge faster. The dynamic simulation is carried out by solving Eq. (1) using the given boundary value q 1 ¼ qðt 1 Þ and time span ½t 1 ; t N . The tail joint torques are computed based on a purely damping controller (Eq. (47)), which will be detailed in the next section. The second technique is to design the objectives or constraints according to the simulation results. For instance, for the airborne righting tasks, the boundary constraints f bc in Eq. (43) are usually set to make the body vertically up at the end of the trajectory, i.e., / x ðt N Þ ¼ 0, / y ðt N Þ ¼ 0. However, the body pitch angle / x and roll angle / y are both modulo 2p numbers. Therefore, strictly driving the body orientation back to / x ¼ 0 and / y ¼ 0 is not necessary. Instead, depending on the simulated results, / x ðt N Þ and / y ðt N Þ could be set as the closest multiples of 2p to the end body orientation. For instance, if the simulation shows that the body pitch angle reaches 15 radians and the body roll angle reaches -5 radians at the end, / x ðt N Þ and / y ðt N Þ constraints could be set to equal to 4p and À2p, respectively. This technique takes advantage of the existing system momentum and tends to make the system perform minimal work to reach the desired orientation.

Trajectory interpolation and joint-level control
The obtained tail trajectory from the NLP is a table of isolated values on collocation points. To construct a continuous trajectory for joint-level control, piecewise spline interpolation is utilized using the MATLAB function spline. To track this tail trajectory in the simulation, a PFL controller that was previously proposed by the authors [45] is used. That is, the tail control effort s ta is calculated to just linearize part of the system dynamics, which is, in this case, the dynamics corresponding to the tail part of the generalized coordinates q t . Therefore, to track the tail trajectory q t;r ðtÞ, the system output is constructed as y ¼ q t À q t;r ðtÞ ð 44Þ To linearize the system to asymptotically stabilize this output, a spring-damper system that is known to be stable is constructed as where K d ¼ K d I 2mÂ2m and K p ¼ K p I 2mÂ2m with K d ; K p [ 0 and I 2mÂ2m 2 R 2mÂ2m being the identity matrix. Substituting Eq. (44) into Eq. (55), using a selection matrix S 2 R 2mÂd such that q t ¼ Sq, and using Eq. (1), Eq. (45) could be solved as in which X ¼ SH À1 J T ta and X þ is the Moore-Penrose inverse of X. Again, the GRF term does not appear due to the airborne assumption. It is worth noting that in practice, we did not use the € q t;r term and it turned out that the controller still worked well. It is also worth noting that due to the coupled and highly nonlinear dynamics between the tail and the body, a simple model-free proportional-derivative (PD) controller turned out not to work well in this case, which is the major reason that we chose to use PFL, which empirically shows better tracking performances.
For the phase when the tail is not tracking a prescribed trajectory, the tail joint controller is set to a passive damping controller, which simply resists the tail motion and purely consumes the system energy. The control effort is calculated as

Numerical experiments
All computations are carried out using MATLAB. The simulation is performed using a variable step ordinary differential equation (ODE) solver ode45 with an absolute tolerance of 1e-6 and a relative tolerance of 1e-8 to guarantee the simulation accuracy. The NLP problem for the tail trajectory optimization is solved using the function fmincon with ''interior-point'' method. The computing environment consists of a workstation with an AMD Ryzen 9 3900X 12 Core CPU (3.793 GHz) and a MATLAB R2021b distribution running on Microsoft Windows 10 Education operating system. All the model properties are summarized in Table 1 where the critical kangaroo rat parameters (body mass, tail mass, tail length) are taken from the measurements in [31]. The detailed tail parameters such as the segment length of the continuum tail or the link mass of the articulated tail, are calculated by evenly distributing the total tail properties onto each segment/link, depending on the actual segment number m and link number n. For instance, the link mass of the articulated tail is determined to be m at ¼ tail mass=mn. The MOI of the body and each link in the articulated tail model are estimated by assuming that each rigid body is an ellipsoid. The COM of each articulated tail link locates at the geometric center of the link, i.e., L j2j ¼ 2L j2c . It is worth noting that the friction coefficient l is set to 1 in the airborne righting motion tests in Sects. 4.2 and 4.3 to guarantee a no-slip condition between the foot and the ground during jumping. This value is later set to 0.5 in the tail-ground contact motion tests in Sect. 4.4 to allow slippage between the tail and the ground. This change is accommodated to facilitate the numerical computations without affecting the motion conclusions.
Since this paper focuses on the tail behaviors, the leg motion is planned manually, and only simple trajectories are used. For instance, for the airborne simulations, the leg motions are designed to be which keeps c l constant and rotates d l in a constant speed from À15 to À75 . t kick is the moment to start the leg motion and in this paper, all the simulations use the same t kick ¼ 0:4s. This leg kicking motion drives the kangaroo rat to jump to a height of around 0.6 m. It is worth noting that for this kicking motion, the kangaroo rat model uses its ankle (A) point and metatarsophalangeal (M) point to push the ground, instead of the toe (T) point. This is similar to the jumping motion exhibited by the real kangaroo rat. For better illustration, animations were created for the simulations in this section, and they can be found online at https://youtu.be/wy067QQ0Cvs. This section investigates the performance difference between the continuum tail model and the articulated tail model. Since the fundamental difference between these two models is that the articulated model discretizes the continuum tail, this comparative study aims to investigate the effects of the number of links on the tail performance. Therefore, the tail segment number is set to 1 (m ¼ 1), and standard airborne yaw rotation motions are simulated. All other model parameters are set to be the same. The tail trajectories are planned so that the articulated tail has the same COM angle as the continuum tail, as shown in Fig. 7. This criterion was previously proposed by the authors to compare the dynamic differences for different robotic tail structures [35]. The tail COM angle h com is defined as the rotation angle of the line segment connecting the tail COM and the tail mounting point T.
To compute the h com of the continuum tail (which is an arc), Eq. (86) in [35] is used, which gives The h com trajectory is planned to be a simple pointto-point motion from 0 radians to atanðp=2Þ radians using cubic polynomial, which corresponds to h 1 from 0 to 180 .
The comparison results are presented in Fig. 8 where the Fig. 8a shows the body yaw responses / z t ð Þ for the continuum tail model and the articulated tail model (four links are used, i.e., n ¼ 4), respectively. Figure 8b plots the trend of the body yaw angle differences as the link number increases. From  Fig. 8b, it can be found that the overall model difference is small (less than 0:5 ) and this difference approaches zero as more links are used for the articulated model. The largest difference happens when the link number is four. However, even for this worst case, the absolute model difference is still small, which can be also observed from Fig. 8a.

Airborne righting tests
From the results of Sect. 4.1, the dynamic difference between the articulated tail model and the continuum tail model is not significant, especially when the link number of the articulated tail becomes large. Therefore, for this section and the next section, the experiments will only focus on the articulated tail structure. This section aims to verify the developed dynamic model and the tail motion controller, where the tail trajectory is generated by the direct collocation method in Sect. 3. The objective function and constraints settings are summarized as follows: . . .; N À 1) to minimize the control effort over a given time span.
from a given initial state qðt 1 Þ and reach a desired body orientation / b;d;N with a zero body angular velocity. (3) g 1 ¼ q t;k À b t;lim e 2m 0, g 2 ¼ Àq t;k À b t;lim e 2m 0, g 3 ¼ s ta;k À b ta;lim e 2m 0, g 4 ¼ Às ta;k À b ta;lim e 2m 0, g 5 ¼ ½t 1 À t N t N À t 1 À t s T 0 to add joint motion limit (b t;lim ¼ 180 =mn), joint torque limit (b ta;lim ¼ 0:02Nm), and motion planning time span t s [ 0. e 2m is a column vector of ones that has 2m elements.
The initial state qðt 1 Þ for the trajectory optimization is selected at a random moment t 1 after the kangaroo rat leaves the ground, although this ''random'' selection has a prerequisite that it must guarantee that there exists at least one feasible solution to the NLP problem. For instance, if t 1 is chosen to be too close to t N , it is possible that there does not exist a feasible tail trajectory that can drive the system to the desired state, i.e., no matter how the tail rotates, the kangaroo rat is not able to reach its desired posture in time. This arbitrary selection of starting moment is meant to show the effectiveness of the motion planning algorithm on handling the airborne righting tasks for any given state. For the experiments in this section, we used a one-segment tail with three links per segment (m ¼ 1, n ¼ 3), selected t 1 ¼ 0:5s (0.1 s after the leg yb xb T θcom Fig. 7 Tail trajectory setting for the comparative study kicking the ground), t s ¼ 0:58s (the moment just before the landing), used a desired final posture of / x ¼ 0, / y ¼ 0, and set a collocation rate of 100 Hz (step size of 10 ms). The snapshots for the entire jumping motion with active tail control are shown in Fig. 9b and c where the dashed lines show the body COM (point B) trajectories. To compare the tail motion effects, a simulation for the same jumping motion but without active tail motion control was also conducted and included in Fig. 9a. The entire motion in Fig. 9a was simulated using the full dynamics of Eq. (1) and the tail controller of Eq. (47). In Fig. 9b, the green segment (green dashed line) of the motion was simulated using Eqs. (1) and (47), but the red segment (red dashed line) was a direct plot of the NLP results in Sect. 3.1, which used all the designed motion information including both the body motion and the tail motion. To verify that the designed motion was actually feasible and can meet the desired goal, Fig. 9c simulated the red part using Eq. (1) and the tail controller in Eq. (46). Note the slight differences of how Fig. 9b used the designed whole-body motion for plotting while Fig. 9c used only the tail trajectory. This is because in simulation, the body motion is the under-actuated DOF and only the tail trajectories could be properly controlled. It is also worth noting that the tail contact model was not included for these experiments and therefore the tail could penetrate the ground in Fig. 9. Comparing Fig. 9a and c validates the tail's functionalities in helping the airborne righting maneuvers of kangaroo rats. Comparing Fig. 9b and c verifies the effectiveness of the motion planning algorithms in Sect. 3.1 although only a simple trapezoidal quadrature rule was used. Figure 10 shows the time domain details for the motions in Fig. 9, where the Fig. 10a plots the body COM trajectory p b ðtÞ (p b;x , p b;y , and p b;z are the x, y, z components, respectively), Fig. 10b and c show the body orientation trajectories, Fig. 10d plots the tail joint trajectories, Fig. 10e plots the trajectory tracking Fig. 9 Snapshots of a a simulated jumping motion without active tail control, b a designed whole-body motion using the proposed motion planning algorithm, and c the simulated motion by tracking the designed tail motion According to Fig. 10a, the matching trend among the three trails verifies the effectiveness of the motion planning algorithm. It also implies that the tail motion has marginal effects on the body COM motion. However, Fig. 10b and c demonstrate the tail's usefulness in affecting the body's orientation, although the mismatch between the planned (blue dot-dashed line) and simulated (green solid line) body orientation trajectories reflect the accuracy limit of using the trapezoidal quadrature rule to implement the dynamics constraints. Figure 10d and e show the effectiveness of the PFL controller for the tail's jointlevel trajectory tracking control. Figure 10f first verifies how the PFL controller can generate appropriate control effort. It also shows how the open-loop motion planning algorithm can generate a similar control effort as the closed-loop feedback controller.
To further investigate the numerical aspects (e.g., accuracy and efficiency) of the motion planning algorithm, different collocation frequencies are used, and the results are collected in Fig. 11. The dynamics error is measured based on the distance between the planned motion x nlp and the simulated motion x track with the planned tail trajectory being tracked, which is a L 2 -norm: From the results, it can be found that increasing the collocation frequency helps reduce the dynamics error but increased computation time. However, when the dynamics error is already reduced to a low level, this trend does not continue ([ 100 Hz in Fig. 11). The iteration number and objective value are shown not to be significantly affected by the collocation frequency.

Comparative study on tail segmentation
After verifying the modeling and control frameworks using a one-segment tail, this section investigates the differences of the tail with multiple segments, that is, how the tail's number of segments affects the airborne righting performance as well as the numerical aspects of the trajectory optimization algorithm. For this purpose, the total number of links of the articulated tail is set to 6 so that it could be divided by the segment numbers 1, 2, and 3. For instance, when the segment number is set to m ¼ 2, the link number is set to n ¼ 3. The remaining settings are the same as in Sect. 4.2 except for changing b t;lim ¼ 360 =mn to allow a larger workspace of the tail to find feasible solutions for the NLP problem.  Table 2 where the tracking error is computed according to maxfabsfe 1;a ; e 1;b ; . . .; e 3;a ; e 3;b gg. Table 2 shows similar observations as in the figures whereas more segments in the tail tend to save energy and reduce peak power by finding better tail trajectories (with lower optimal value). However, this trend is not significant between the two-segment case and the three-segment case. It is worth noting that Table 2 shows a noticeable dynamics error (the larger the error, the larger the difference between the planned motion and the simulated motion). However, based on the conclusions from Sect. 4.2, this issue could be solved by increasing the collocation frequency.

Tail-ground contact tests
When the kangaroo rat is on the ground, the tail could be used as an additional appendage to assist the locomotion. This section investigates this static tail functionality by conducting case studies of a jumping motion and a self-recovery motion. Both the articulated tail model and the continuum tail model are included to evaluate the effectiveness of the contact model established in Sect. 2.4. Due to the large number of contact events, the non-stiff ODE solver ode45 is replaced by a stiff ODE solver ode15s for faster evaluation of the contact events. The friction coefficient l is also adjusted to 0.5 to allow sliding motion on the contact surface.

Case study on jumping motion
Two sets of experiments were designed to compare the difference between the two tail types: one jumping motion using the articulated tail and one jumping motion using the continuum tail. Both cases utilized the tail contact event to stabilize locomotion. For comparison purposes, one additional jumping motion without using the tail for locomotion was also performed. It is worth noting that other than changing the ODE solver and the friction coefficient, the leg trajectory component d l ðtÞ was also changed to rotating from À40 to À90 for kicking and then rotating back to À30 for landing. Both tails used three segments (m ¼ 3) and the articulated tail used four links per segment (n ¼ 4).
The simulation results are demonstrated in Fig. 14, where Fig. 14a-d show the jumping motion without using the tail support on the ground, Fig. 14e-h show the jumping motion using the articulated tail support on the ground, and Fig. 14i-l show the jumping motion using the continuum tail support on the ground. From the case study results, it can be found that the tailground contact model is able to generate realistic contact phenomena, which enables future investigations on contact-rich motion planning and control of the kangaroo rat model. In addition, the case study shows that through environmental contacts, the tail Fig. 13 Detailed plots of the planned motion for one-(subfigures (a), (d), (g)), two-(subfigures (b), (e), (h)), and three-segment (subfigures (c), (f), (i)) tail case. The line style code for the subfigures a-c is that the ''red solid'', ''green dashed'', and ''blue dash-dotted'' lines correspond to / x , / y , and / z , respectively. The line style code for the subfigures d-f is that the ''red solid'' line, ''green dashed'' line, ''blue dash-dotted'' line, ''black dotted'' line, ''black solid'' line with circle markers, and the ''black solid'' line with square markers correspond to a 1 , b 1 , a 2 , b 2 , a 3 , and b 3 , respectively. The line style code for the subfigures g-i is that the ''red solid'' line, ''green dashed'' line, ''blue dash-dotted'' line, ''black dotted'' line, ''black solid'' line with circle markers, and the ''black solid'' line with square markers correspond to s ta;a;1 , s ta;b;1 , s ta;a;2 , s ta;b;2 , s ta;a;3 , and s ta;b;3 , respectively

Case study on self-recovery motion
For the self-recovery test, the robot was initialized laying down on the ground. The tail motion was manually planned to push the body to stand up. The results are presented in Fig. 15 where Fig. 15a-d show the snapshots of the self-recovery motion of the kangaroo rat using the articulated tail and Fig. 15e-h show the snapshots of the same motion using the continuum tail. Note that due to the soft contact model, the tail may penetrate the ground for transient motions. From Fig. 15, the tail-ground contact models are verified further, especially for the varying contact point problem (Eqs. 38 and 39) of the continuum tail.
The importance of this case study is that it illustrates how the tail could be used to help the manipulation tasks of the kangaroo rat through tail-environment contacts. In these scenarios, the serpentine tail structure helps to enlarge the tail workspace and enhance the tail dexterity significantly, which is a fundamental benefit in comparison with single-link tails.

Discussion
The tail controller in this paper mostly focuses on the airborne righting motions, for which the dynamics are continuous. In future work, we wish to extend the work into the non-smooth domain, i.e., designing the tail controller together with the leg motion planning, which involves tail-ground contact and leg-ground contact. Due to the rich contact type and large contact number, the traditional method based on mode schedule [48] is less attractive. To solve this problem, through-contact trajectory optimization may be required, such as those based on direct methods [49,50] and those based on indirect methods [51,52]. This work was motivated by understanding the tail's functionalities on agile motions of the kangaroo rat from dynamics and control perspectives. This understanding, together with the modeling and control framework developed in this paper, aims to lay the foundation for further developments of agile biped robots with a serpentine robotic tail. Therefore, one important future work is to develop such a robotic system that achieves the same level of agility and dexterity as the kangaroo rat, which may have unique applications as a highly agile terrestrial moving platform. Although the proposed dynamic model computes fast (less than 0.1 ms), the motion planning algorithm takes a considerable amount of time (from several minutes to a couple of hours) to converge. To mitigate this drawback in the future, an analytical model with automatic differentiation [53] may be used to fasten the NLP solving process. Moreover, due to the nonlinear nature of the system dynamics, the motion planning algorithm is only able to find a local minimum. To find the global minimum or a unique solution for every motion planning, the model complexity may need to be reduced so that the NLP problem could be formulated as a convex programming problem [54]. In addition, the more generic techniques (e.g., genetic algorithm, Monte-Carlo method) may also be used to find a global optimum or a suboptimum at the cost of a longer computation time. For the specific case in this paper, one practical way is to generate more initial guesses spanning the entire feasible region, and then use a small disturbance to reinitiate the optimization process after finding a local minimum. This technique may help to find a better solution but global optimum is still not guaranteed.

Conclusion
This work investigated the agile motions and tail functionalities of the kangaroo rat from the perspective of dynamics and control. Two representative, general serpentine tail models were proposed: one articulated tail model and one continuum tail model. The contact events among the feet, tail, and the ground were modeled using a regularized compliant contact model. To automatically design the tail motion, a direct collocation-based numerical optimal control method was utilized. The planned tail trajectories were then tracked using a partial feedback linearization controller. Using the established dynamic model and control framework, various numerical experiments were performed. A comparative study on the two serpentine tail models was conducted first and the results validated the common impression that the articulated tail model approaches the continuum tail model as the number of links increases. Two sets of airborne righting tests were then performed, focusing on verifying the dynamics and control framework and investigating the differences between tail segmentations, respectively. To test the tail contact model and explore the tail's functionalities in supporting the body, case studies on jumping and self-recovery motions were conducted. The simulation results demonstrated the unique functionalities of the serpentine tail with assisting the kangaroo rat's locomotion.
Data availability The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Code availability Custom code.

Declarations
Conflicts of interest The authors declare that they have no conflict of interest.

Appendix A: Articulated tail kinematics
The velocities, Jacobians, acceleration, and MOI for each link of the articulated tail model are computed recursively using Eqs. (A1-A13). Eq. (A14) computes the torso MOI. u x;y is an x-dimension unit column vector with 1 on the y-th entry.
J j;com ¼ J j;jnt þ J j;j2c ðA6Þ _ v j;j2c ¼ e _ x j p j;j2c þx 2 j p j;j2c _ v j;j2j ¼ e _ x j p j;j2j þx 2 j p j;j2j ( ðA12Þ I j;at ¼ R j j I j;at R T j ðA13Þ Appendix B: Continuum tail kinematics The detailed expression of matrix E i;v is given as