Research on Dynamic Modeling Method and Flying Gait Characteristics of Quadruped Robots with Flexible Spines

1. Introduction

The high adaptability and agile movement characteristics of quadruped robots in unstructured environments make them significant in specific applications and research areas, such as earthquake or fire rescue operations, wilderness military support, indoor household services, and educational entertainment. Consequently, quadruped robots have become a focal point in the field of robotics. Currently, research on quadruped robots mostly revolves around rigid-bodied quadruped robots [1-5], with a primary focus on bio-inspired leg structures and the study of gait generation and coordinated control among the four limbs. In practical quadruped mammals, however, the flexible spine joints play a crucial role in achieving high-speed, agile, and stable motion.

In recent years, there have been two main types of spine structures adopted for quadruped robots, namely articulated spines [6-8] and integral compliant spines[9-14]. Although the mechanisms by which a flexible spine influences dynamic performance remain elusive, domestic and international research suggests that a flexible spine can enhance a robot's movement speed, increase stride length, enable elastic energy storage, and reduce energy consumption [6-8, 10-14]. Stoch2, a quadruped robot designed by Shounak Bhattacharya and colleagues at the Indian Institute of Science, features a dual-joint spine. Drawing inspiration from the running motion of cheetahs, this robot utilizes deep reinforcement learning to generate joint control signals. The results indicate a significant improvement in the robot's running speed, and stride length, and a reduction in locomotion costs facilitated by the flexible spine [6]. Sugoi-Neco, a passive single-joint spine quadruped robot designed by Ryosuke Kawasaki and his team at the University of Electro-Communications in Japan, employs a spine joint with torsional springs for storing elastic potential energy. The findings suggest that the energy stored in the flexible spine allows the robot to generate larger horizontal and vertical impulses during the support phase, enabling it to achieve greater leaps [8]. The MIT Cheetah I robot, designed by Sangbae Kim's team at the Massachusetts Institute of Technology, features a spine composed of rigid vertebrae interspersed with polyurethane rubber rings. The spine is coupled to the hind legs by connecting ropes, utilizing a differential mechanism linked to the hip joint motors. This coupling achieves synchronized motion between the spine and hind legs. The robot ultimately achieved a running speed of 22 km/h in a trot gait on a treadmill, with the flexible spine contributing to reducing the robot's locomotion costs [9, 10]. The Kitty robot, designed by Zhao Qian and colleagues at the University of Zurich, has a spine consisting of ABS-made cross-shaped vertebrae alternately connected with silicone blocks. Four servos drive ropes in four directions to bend the spine. It was found that the robot exhibits the fastest forward speed when the virtual spine joint is located on the posterior side of the spine [11]. Canid and Inu are two integral compliant spine quadruped robots developed with the involvement of Daniel E. Koditschek's team at the University of Pennsylvania. The spines of these robots directly use carbon fiber plate springs and glass fiber leaf spring plates, respectively, as the main body of the spine. Motors rotate to actively control the spine by pulling ropes or belts passing through it. Experimental results indicate that a flexible spine can increase the average speed of quadruped robots by 0.4 m/s, decrease the body pitch angle by 17%, and reduce the average vertical height of the center of mass by approximately 3 cm [12, 13]. The passive integral compliant spine quadruped robot, Fanari, designed by Hasan H. Kani and colleagues at the University of Tehran in Iran, features a spine composed of multiple organic glass segments parallelly connected with linear springs. In experiments where the robot sprints along a slope without any external energy, it was observed that the spine, capable of bending in both upward and downward directions, enables the robot to achieve faster and more stable movement compared to spines that can only bend upward or rigid spines [14].

Due to the additional degrees of freedom introduced by a flexible spine in quadruped robots, and the necessity for the motion of the spine to be coordinated with the movement of the remaining degrees of freedom, without mutual inhibition, numerous researchers worldwide are currently exploring the mechanisms underlying the impact of a flexible spine on various performance metrics of quadruped robots. They are establishing corresponding dynamic models and conducting motion characteristic analyses [15-22]. However, there is limited research investigating the influence of a flexible spine on the flight phase motion of quadruped robots. In 2015, Soroush Maleki and colleagues from the University of Tehran addressed quadruped robots with flexible spines. Using the Lagrange method, they formulated motion equations and, based on the model, employed a genetic algorithm for gait design optimization, ultimately reducing the system's energy consumption [16]. In 2018, Yesilevskiy proposed a quadruped robot model with distributed mass. By comparing the dynamic characteristics of rigid-spine and flexible-spine quadruped robots, it was found that, for asymmetric gaits, introducing a flexible spine can increase the maximum running speed, extend the stride, and improve energy utilization efficiency through the energy storage and conversion functions of the spine [21]. In 2020, Zeng Shun and colleagues from Wuhan University of Technology investigated the influence of the center of mass position of spine segments on the dynamic performance of the bound gait in quadruped robots. They established a passive simplified model in the sagittal plane, and treated the center of mass position of the spine segment as an independent variable. The results indicated that when the center of mass of the spine segment is closer to the hip joint, the bound gait cycle and stride will be longer. Conversely, when the center of mass is closer to the midpoint of the spine segment, the fluctuation of the robot's center of gravity and the ground reaction force will be smaller [22].

In this paper, we concentrate on quadruped robots with a flexible spine and formulate dynamic models with varying degrees of simplification. Drawing inspiration from the observed flight phase motion in cheetahs during high-speed galloping, we introduce parameterized descriptions of the flight phase process for both rigid-bodied and flexible-spine quadruped robots. To assess the impact of a flexible spine on the flight phase, we assign different initial values to the spine joint and compute the flight phase outcomes for both types of quadruped robots at various initial speeds. Through comparative experiments, we affirm the influence of a flexible spine during the flight phase. The primary contributions of this paper include the establishment of dynamic models of different degrees and the parameterized description of the flight phase for quadruped robots. The study validates the influence of a flexible spine during the flight phase, inspired by observations of cheetahs in high-speed galloping.

The remainder of the paper is organized as follows: In Section 2, the dynamic models of quadruped robots with different degrees of simplification are introduced. Section 3 provides a parameterized description of the flight phase process for both rigid-bodied and flexible-spine quadruped robots. Section 4 compares the impact of rigid-bodied and flexible-spine configurations on the flight phase of quadruped robots and analyzes the results of this comparison. Finally, Section 5 concludes the paper and discusses future work.

2. Quadruped Robot Model

2.1. Dynamic Models of Varying Degrees of Simplification

We have established a quadruped robot model with a spine, as illustrated in Figure 1. This model consists of a body with a spine and four legs, each leg comprising a thigh and a shank. The model has 9 degrees of freedom, with each leg having two degrees of freedom: knee joint pitch and hip joint pitch. There is no hip roll degree of freedom, and there is one spine joint degree of freedom. The parameters of this model are presented in Table 1.

2.1.1. Mass-Point Model

Due to the prevalent design of low inertia in the leg structure and the centralized arrangement of drive motors at the hip joints in contemporary quadruped robots [2], this simplified model considers only the combined center of mass of the two trunk segments. The model depicted in Figure 1 is thus simplified to a point mass model, where the mass possesses no volume or shape considerations, as shown in Figure 2(a). The external forces acting on it include constraint forces from the ground and the gravitational forces of the two trunk segments. The forces can be entirely translated to a single point, forming a concurrent force system. The position vector of the modeled point mass relative to the world coordinate system can be expressed as:

$${\mathit{r}}_{\mathrm{c}}=\frac{m_{\mathrm{bB}}{\mathit{r}}_{\mathrm{bB}}+m_{\mathrm{bF}}{\mathit{r}}_{\mathrm{bF}}}{m_{\mathrm{bB}}+m_{\mathrm{bF}}}=\frac{m_{\mathrm{bB}}{\mathit{r}}_{\mathrm{bB}}+m_{\mathrm{bF}}{\mathit{r}}_{\mathrm{bF}}}{m_0}$$

(1)

Where, $m_0$ represents the mass of the robot's trunk, and ${\mathit{r}}_c$ is the position vector of the point mass in the world coordinate system. Since forces only influence changes in the position of the point mass, based on Newton's Second Law, the dynamic equation for this point mass can be derived as:

$$m_0\mathit{g}+{\displaystyle \sum _{i=1}^4{\mathsf{\lambda }}_i}=m_0 · {\ddot{\mathit{r}}}_{\mathrm{c}}$$

(2)

Where, ${\lambda }_i$ represents the constraint force applied by the ground at the end of each leg of the robot.

2.1.2. Single Rigid Body Model

The simplification in the single rigid body model involves neglecting the dynamic characteristics of the robot's legs, considering the robot's torso as a single rigid body, and assuming the robot has only a torso and foot ends. The model is illustrated in Figure 2(b), where the coordinate system $\left\{B\right\}$ represents the body-attached coordinate system.

The external forces acting on the four-legged robot with a single rigid body consist of the body gravity and the ground reaction forces at the foot ends. For the translational motion of the rigid body, the translational motion equation of this model can be derived from Euler's first law of motion as follows:

$$m_0\mathit{g}+{\displaystyle \sum _{i=1}^4{}_W\mathsf{\lambda }{ }_i}=m_0·{}_W\ddot{\mathit{r}}{ }_{\mathrm{c}}$$

(3)

For the rotational motion of the rigid body, choosing the body's center of mass as the reference point, the external torque on the rigid body can be obtained from Euler's second law as:

$${}_W\mathit{M}={\displaystyle \sum _{i=1}^4{}_W\mathit{r}{ }_{cti}\times {}_W\mathit{\lambda }{ }_i}+0_{3\times 1}\times \left(m_0\cdot {}_W\mathit{g}\right)$$

(4)

Where, $M$ represents the external torque acting on the rigid body.

The angular momentum $L$ of the rigid body, rotating about an axis passing through its center of mass and perpendicular to the plane of the paper, is equal to the product of the rotational inertia $J$ and the angular velocity ${\omega }_{WB}$ of the rigid body relative to the world coordinate system. This relationship is expressed in Equation (5). Here, the rotational inertia $J$ about the axis passing through the center of mass and perpendicular to the paper can be calculated as $J=\frac{1}{12}\cdot m_0\cdot {\left(l_{bB}+l_{bF}\right)}^2$ ，where $l_{bB}$ is the distance from the center of mass to the back foot, and $l_{bF}$ is the distance from the center of mass to the front foot.

$${}_W\mathit{L}=J·{}_W\mathsf{\omega }{ }_{WB}$$

(5)

Thus, Equations (3) and (6) together form the dynamic equations of the single rigid body model.

According to Euler's second law, combining Equations (4) and (5), we can express the equations of motion for the rigid body undergoing rotational motion about a fixed axis in the plane as:

$${\displaystyle \sum _{i=1}^4{}_W\mathit{r}{ }_{cti}\times {}_W\mathsf{\lambda }{ }_i}+0_{3\times 1}\times \left(m_0\cdot {}_W\mathit{g}\right)=\frac{\mathrm{d}{}_W\mathit{L}}{\mathrm{d}t}=J·\frac{\mathrm{d}{}_W\mathsf{\omega }{ }_{WB}}{\mathrm{d}t}$$

(6)

2.1.3. Planar Multi-Rigid Body Model

Unlike the point mass model and the single rigid body model, the planar multi-rigid body model considers the legs explicitly and precisely accounts for the motion of each body in the robot. As shown in Figure 3, a planar multi-rigid body model is established, treating the quadruped robot as a floating base object. The motion of the floating base with respect to the world coordinate system is achieved through interactions between the foot ends and the ground, without direct actuation. This model is formulated using the first kind of Lagrange equations, and according to the first kind of Lagrange formalism, dynamic equations are derived for the abstracted 12-DOF model depicted in Figure 3.

$$\frac{d}{dt}(\frac{\partial L}{\partial {\dot{q}}_k})-\frac{\partial L}{\partial q_k}=Q_k+{\Gamma }_k\begin{array}{cc} & \end{array} \left(k=1,2,\cdots ,12\right)$$

(7)

To express each term of Equation (7), write the equation for each generalized coordinate's motion. Then, organize these 12 equations into matrix form, which can be found in the appendix A. Finally, provide a concise representation of the matrix form.

$$\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})+\mathit{G}(\mathit{q})=\mathit{B}\cdot \mathit{u}+{\mathit{J}}^{\mathrm{T}}\mathsf{\lambda }$$

(8)

Where, the coefficient matrix $M\left(q\right)$ is referred to as the system's inertia matrix, $C\left(q,\dot{q}\right)$ is the Coriolis and centrifugal vector matrix of the robot, $G\left(q\right)$ is the gravity vector matrix of the robot, and $B$ is the selection matrix. The 8-dimensional input torque vector $u$ ,representing the torques applied at each joint, is transformed into a 12-dimensional column vector. The first four elements of the vector are zeros, and the last eight elements are equal to the input torque vector $u$ . $J$ is the stacked Jacobian matrix of the ground-contacted foot ends.

3. Quadruped robot flight phase motion parameterization

3.1. Flight phase motion simplification

To simplify the analysis, this paper chooses the process starting from the cheetah's hind legs touching the ground, through the extension of the trunk, to the realization of the large flight phase, and finally the front legs touching the ground as a reference [23]. According to the force and motion characteristics of the cheetah, the whole process is divided into take-off motion and flight phase motion, which is subjected to gravity and ground reaction force during the take-off motion, and only gravity during the flight phase motion, as shown in Figure 4.

Referring to the motion characteristics of the cheetah's large flight phase process in Figure 4, the flight phase process of the quadruped robot with rigid torso and flexible spine is simplified as Figure 5and Figure 6. In the simplified process, the left and right pairs of legs of the robot are set to have the same motion phase, and in the take-off motion, only the robot's back leg exerts force, and the angle of the joints of the front leg does not change, and the angle of each joint of the front and back legs is set to remain constant until the end of the flight phase motion. The angles of the front and back legs of the robot remain unchanged until the end of the flight phase motion, and the angles of the spine joints are set to change only during the takeoff motion, while the angles remain fixed during the flight phase motion.

3.2. Rigid torso quadruped robot flight phase motion parameterization

The flight phase motion of a rigid torso quadruped robot can be decomposed into translational motion of the center of mass and fixed-axis rotation over the center of mass and perpendicular to the plane, as shown in Figure 7.

In the take-off motion process, only by the vertical direction of gravity, the whole process is equivalent to the oblique throw movement, to take off the point as a benchmark, combined with the expected to achieve the flight phase trajectory height of $h$ and the landing point (LA) coordinates for $\left(x_{LA}, y_{LA}\right)$ can be known as the shape of the parabola, and can be obtained center of mass in the horizontal and vertical direction of the speed and position change trajectory, the time of the fight phase, the angle of the landing point:

$$\left\{\begin{array}{l}{v}_{cx\_FL}\left(t\right)=v_{LOx},v_{cy\_FL}\left(t\right)=v_{LOy}-\mathit{g}t\\ x_{c\_FL}\left(t\right)=v_{LOx}·t ,y_{c\_FL}\left(t\right)=v_{LOy}t-\frac{1}{2}\mathit{g}{t}^2\\ \mathsf{\Delta }{t}_{\mathrm{FL}}=\frac{\sqrt{2\mathit{g}h}+\sqrt{2\mathit{g}h-2\mathit{g}{y}_{\mathrm{LA}}}}{\mathit{g}}\\ {\theta }_{\mathrm{LA}}={\theta }_{\mathrm{LO}}+{\omega }_{\mathrm{LO}}·\mathsf{\Delta }{t}_{\mathrm{FL}}\end{array}\right.$$

(9)

In the flight phase motion process, according to Ludovic D. Maes et al. on the biology of dogs, it was shown that there is a certain pattern between the gait frequency $f$ as well as the loading factor $\beta$ and the horizontal velocity v_x, which corresponds to the mathematical relationship [24]:

$$\left\{\begin{array}{l}f=1.314+0.762\ln v_{\mathrm{x}}\\ \beta =0.652-0.21\ln v_{\mathrm{x}}\end{array}\right.$$

(10)

The horizontal instantaneous velocity at the onset point is approximately equal to its average velocity during a gait cycle[23]. The duration $△t_{\mathrm{TO}}$ of the takeoff motion can be obtained as:

$$△t_{\mathrm{TO}}=\frac{\beta }{f}=\frac{0.652-0.21\ln v_{\mathrm{STx}}}{1.314+0.762\ln v_{\mathrm{STx}}}$$

(11)

According to the momentum theorem, the horizontal and vertical ground reaction forces ${\lambda }_x\left(t\right)$ and ${\lambda }_y\left(t\right)$ applied to the robot during the take-off motion should satisfy the following conditions:

$${\displaystyle \underset{0}{\overset{\frac{\beta }{f}}{\int }}{\lambda }_{\mathrm{x}}(t)} dt=m(v_{\mathrm{LOx}}-v_{\mathrm{STx}}) \Rightarrow {\displaystyle \underset{0}{\overset{\frac{\beta }{f}}{\int }}{\lambda }_{\mathrm{x}}(t)} dt=\frac{m\mathit{g}{x}_{\mathrm{LA}}}{\sqrt{2\mathit{g}h}+\sqrt{2\mathit{g}h-2\mathit{g}{y}_{\mathrm{LA}}}}-m{v}_{\mathrm{STx}}$$

(12)

$${\displaystyle \underset{0}{\overset{\frac{\beta }{f}}{\int }}{\lambda }_{\mathrm{y}}(t)}-m_0\mathit{g} dt=m(v_{\mathrm{LOy}}-v_{\mathrm{STy}}) \Rightarrow {\displaystyle \underset{0}{\overset{\frac{\beta }{f}}{\int }}{\lambda }_{\mathrm{y}}(t)} dt=m\mathit{g}\frac{\beta }{f}+m(\sqrt{2\mathit{g}h}-v_{\mathrm{STy}})$$

(13)

Based on the above conditions, it is still not possible to obtain specific trajectories of ground reaction forces. Referring to the trajectories of ground reaction force in gallop gait measured by Penny E. Hudson et al. for cheetahs and greyhounds running at high speed[25] and Rebecca M. Walter and David R. Carrier for adult dogs running at high speed[26], Fourier series were used to fit the trajectories of ground reaction forces over time.

For the vertical ground reaction force ${\lambda }_y\left(t\right)$ , $\left(0,0\right), \left(0.4\beta /f, {\lambda }_y^{*}\right), \left(\beta /f, 0\right)$ three points are chosen as fitting points, and at the same time, the impulse constraints of equation (13) should be satisfied, so the vertical ground reaction force trajectory ${\lambda }_y\left(t\right)$ meets the bionic characteristics, and also meets the conditions of impulse in the take-off motion. As shown in equation (14).

$${\lambda }_{\mathrm{y}}(t)=a_1\sin (2\pi ft)+a_2\sin 2(2\pi ft)$$

(14)

Where, $a_1,a_2$ are both Fourier series coefficients.

Similarly, on the horizontal direction ground reaction force ${\lambda }_x\left(t\right)$ , five points $\left(0,0\right)$ , $\left(0.14\beta /f,-{\lambda }_x^{*}\right), \left(0.28\beta /f, 0\right), \left(0.8\beta /f, 0.3{\lambda }_x^{*}\right), \left(\beta /f, 0\right)$ are selected as the fitting points, and at the same time satisfy Equation (12), then the trajectory of the ground reaction force ${\lambda }_x\left(t\right)$ in the horizontal direction can be expressed as:

$$\begin{array}{c}{\lambda }_{\mathrm{x}}(t)=b_1\sin (2\pi ft)+b_2\sin 2(2\pi ft)+b_3\sin 3(2\pi ft)\\ +b_4\sin 4(2\pi ft)+b_5\cos 5(2\pi ft)\end{array}$$

(15)

Where, $b_1~b_5$ are all Fourier series coefficients.

After determining the trajectory of the ground reaction force in the takeoff motion, the trajectory of the center of mass in the process in terms of velocity and position over time can be determined, as well as the angular trajectory of the torso, as shown in Equation (16).

$$\left\{\begin{array}{l}{v}_{\mathrm{cy}\text{\_}\mathrm{TO}}(t)=v_{\mathrm{STy}}+{\displaystyle \underset{0}{\overset{t}{\int }}{a}_{\mathrm{y}\text{\_}\mathrm{TO}}(t) dt}=v_{\mathrm{STy}}+{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{y}}(t)-m_0\mathit{g}}{m_0} dt}\\ y{}_{\mathrm{c}\text{\_}\mathrm{TO}}(_t^{)}=y_{\mathrm{ST}}+v_{\mathrm{STy}}·t+{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{y}}(t)-m_0\mathit{g}}{m_0} dtdt}}\\ v_{\mathrm{cx}\text{\_}\mathrm{TO}}(t)=v_{\mathrm{STx}}+{\displaystyle \underset{0}{\overset{t}{\int }}{a}_{\mathrm{x}\text{\_}\mathrm{TO}}(t) dt}=v_{\mathrm{STx}}+{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{x}}(t)}{m_0} dt}\\ x{}_{\mathrm{c}\text{\_}\mathrm{TO}}(_t^{)}=x_{\mathrm{ST}}+v_{\mathrm{ST}\mathrm{x}}·t+{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{x}}(t)}{m_0} dtdt}}\\ {\theta }_{\mathrm{TO}}\left(t\right)={\theta }_{\mathrm{ST}}+{\omega }_{\mathrm{ST}}·t+{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{M\left(t\right)}{J}} dtdt}\end{array}\right.$$

(16)

In summary, the parameterization of the flying phase process of the rigid torso quadruped robot is completed. As shown in Table 2.

3.3. Flexible spine quadruped robot flight phase motion parameterization

The overall center of mass of the torso composed of the centers of mass of the anterior and posterior torsos is selected as the object of the analyzed mass point, and the motion process of this mass point in the two-dimensional plane under a specific force is analyzed and described parametrically to represent the flight phase motion process of the flexible spine quadruped robot in a similar manner, which is shown in the schematic diagram in Figure 8(a). In order to describe the parameterized flight phase motion process of the flexible spine quadruped robot, the trajectories of the spine joints in the whole motion process and the ground reaction force in the take-off motion process need to be clarified first.

3.3.1. Spinal joint trajectory planning

Referring to the changing rule of flexible spine in cheetah running set the spine joints only in the take-off motion occurs obvious stretching action, while in the flight phase segment spine joint angle remains unchanged. The flexible spine joints are set to have the locked/unlocked function, and this locked/unlocked function of the flexible spine has been applied to the physical prototype in the research of other domestic counterparts, which is feasible [19]. The motion of the spinal joints is shown schematically in Figure 8(b).

Let the angle of the spinal joints at the start point be ${\phi }_{s0}$ and the stiffness of the joint torsion springs be. the anterior trunk will be subjected to a joint torque $M_s\left(t\right)$ due to the bending of the torsion springs, and $M_s\left(t\right)$ can be expressed as:

$$M_{\mathrm{s}}\left(t\right)=-k_{\mathrm{s}}·{\phi }_s\left(t\right)=J_{\mathrm{bF}}{\ddot{\phi }}_s\left(t\right)=\frac{1}{3}{m}_{\mathrm{bF}}{l}^2{ }_{\mathrm{bF}}·{\ddot{\phi }}_s\left(t\right)$$

(17)

The solution corresponding to the spinal joint angle ${\phi }_s\left(t\right)$ can be found as:

$${\phi }_s\left(t\right)=A\sin \left({\omega }_{\mathrm{n}}t+\phi \right)$$

(18)

where $A$ is the amplitude of the vibration, ${\omega }_n$ is the circular frequency, and it is easy to know that $\phi$ is the initial phase of the vibration. According to the initial condition, ${\phi }_s\left(0\right)={\phi }_{s0}$ ,it can be solved that the amplitude of the vibration $A={\phi }_{s0}$ and the initial phase $\phi =\pi /2$ .

From the circular frequency ${\omega }_n$ of the system, the period $T$ of the vibration can be written:

$$T=\frac{2\pi }{{\omega }_{\mathrm{n}}}$$

(19)

The time of the angular motion of the joint is half of its vibration period and at the same time is equal to the time of the take-off motion $\mathsf{\Delta }{t}_{TO\_flex}$ . Let the initial speeds of the center of mass at its starting point in the horizontal and vertical directions be $v_{STx\_flex}$ , $v_{STy\_flex}$ . Based on Equation (10), the load factor ${\beta }_{flex}$ as well as the gait frequency $f_{flex}$ can be solved for the flexible spine quadruped robot under this motion. Then：

$$\frac{T}{2}=\sqrt{\frac{{\pi }^2{m}_{\mathrm{bF}}{l}_{\mathrm{bF}}{ }^2}{3k_{\mathrm{s}}}}=△t_{\mathrm{TO}\text{\_}\mathrm{flex}}=\frac{{\beta }_{\mathrm{flex}}}{f_{\mathrm{flex}}}$$

(20)

The stiffness $k_s$ of the torsion spring of the spinal joint is:

$$k_{\mathrm{s}}=\frac{{\pi }^2{m}_{\mathrm{bF}}{l}_{\mathrm{bF}}{ }^2{f}_{\mathrm{flex}}{ }^2}{3{\beta }_{\mathrm{flex}}{ }^2}$$

(21)

The obtained stiffness $k_s$ , amplitude $A$ , and initial phase $\phi$ of the torsion spring of the spine joint are brought into Equation (18), and the derivation of the time is made, which means that the angular position of the spine joint and the trajectory of the angular velocity with time can be obtained:

$$\left\{\begin{array}{l}{\phi }_s\left(t\right)=\left\{\begin{array}{c}{\phi }_s{ }_0\sin \left(\frac{f_{\mathrm{flex}}}{{\beta }_{\mathrm{flex}}}\pi t+\frac{\pi }{2}\right)t\in \left[0,{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}}\right] \\ -{\phi }_s{ }_{0}t\in \left[{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}},{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}}+△t_{\mathrm{FL}\text{\_}\mathrm{flex}}\right]\end{array}\right.\\ {\dot{\phi }}_s\left(t\right)=\left\{\begin{array}{c}\frac{{\phi }_s{ }_0\pi f_{\mathrm{flex}}}{{\beta }_{\mathrm{flex}}}\cos \left(\frac{f_{\mathrm{flex}}}{{\beta }_{\mathrm{flex}}}\pi t+\frac{\pi }{2}\right)t\in \left[0,{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}}\right]\\ 0t\in \left[{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}},{\beta }_{\mathrm{flex}}/f_{\mathrm{flex}}+△t_{\mathrm{FL}\text{\_}\mathrm{flex}}\right]\end{array}\right.\end{array}\right.$$

(22)

3.3.2. Ground Reaction Force

In order to obtain the ground reaction force which is only affected by the flexible spine compared to the rigid torso quadruped robot, the solution process is based on the ground reaction force of the rigid torso quadruped robot. Firstly, according to the described flying phase process of the rigid torso quadruped robot, the final obtained trajectories of the center of mass position, velocity, torso attitude angle, and angular velocity of the rigid torso over time are converted into the trajectories of the generalized position and generalized velocity variables over time in the complete dynamics model based on the geometric relationships. The geometric relationship between the above several physical quantities is schematically shown in Figure 9(a) set the hip and knee angles of the front leg to maintain fixed values throughout the process. The specific values of both are set here as shown in equation (23). After determining the trajectories of the variables of the generalized position, the generalized velocity is then derived for the obtained generalized position with respect to time.

$$\left\{\begin{array}{l}{\alpha }_{\mathrm{hFL}\text{\_}\mathrm{rigid}\left(\mathrm{hFR}\text{\_}\mathrm{rigid}\right)}=\frac{10}{9}\pi \\ {\alpha }_{\mathrm{kFL}\text{\_}\mathrm{rigid}\left(\mathrm{kFR}\text{\_}\mathrm{rigid}\right)}=\frac{7}{9}\pi \end{array}\right.$$

(23)

The input torque $u_{\mathit{r}i\mathit{g}id}\left(t\right)$ for each joint of the rear leg of the transition segment can be found by bringing the ground reaction force of the rigid torso quadruped robot based on the Fourier series fit with the derived generalized position and velocity trajectories into the planar multi-rigid body model (Equation (18)).

3.3.3. Flight Phase Motion Parameterization

Since the position of the total center of mass is related to the position of the centers of mass of the front and rear trunks, the motions of the respective centers of mass of the front and rear trunks under the corresponding forces are first analyzed in the take-off motion. Subsequently, the motions of both are synthesized by mass distribution to obtain the motion of the total center of mass in the takeoff motion. A schematic diagram of the force and the corresponding velocity analysis for the front and back trunks is shown in Figure 9(b). The velocity $v_{bB}$ of the center of mass of the posterior torso is only affected by its own gravity $G_{bB}$ and the ground reaction force ${\lambda }_{flex}\left(t\right)$ , while the velocity of the center of mass of the anterior torso is not only affected by its own gravity $G_{bF}$ and the ground reaction force ${\lambda }_{flex}\left(t\right)$ , but also by the elastic force $F$ corresponding to the release of flexible spine, and the two kinds of velocities are expressed as and , respectively. The final combined velocity of and is the velocity of the anterior trunk center of mass in takeoff motion.

Let the horizontal and vertical initial velocities of the robot's total center of mass as well as the front and rear torsos at the starting point be $v_{STx\_flex}$ , $v_{STy\_flex}$ . then, according to the ground reaction force ${\lambda }_{flex\_x}\left(t\right)$ , ${\lambda }_{flex\_y}\left(t\right)$ of the flexible spine object obtained in the previous section, we can obtain the trajectory of the center of mass velocity of the front and rear torsos' centers of mass in the presence of the ground reaction force with time:

$$\left\{\begin{array}{l}{v}_{\mathrm{bBx}}\left(t\right)=v_{\mathrm{STx}\text{\_}\mathrm{flex}}+{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{flex}}{ }_{\_\mathrm{x}}(t)}{m_0} }dt=v{`}_{\mathrm{bFx}}\left(t\right)\\ v_{\mathrm{bBy}}\left(t\right)=v_{\mathrm{STy}\text{\_}\mathrm{flex}}+{\displaystyle \underset{0}{\overset{t}{\int }}\frac{{\lambda }_{\mathrm{flex}\text{\_}\mathrm{y}}(t)-m_0\mathit{g}}{m_0} }dt=v{`}_{\mathrm{bFy}}\left(t\right)\end{array}\right.$$

(24)

Since the quadrupedal robot with a flexible torso is reduced to a prime model, the angle between its rear side torso and the horizontal direction cannot be obtained, so this quantity is set equal to the torso attitude angle of the quadrupedal robot with a rigid torso. It can be obtained：

$$\left\{\begin{array}{l}v{*}_{\mathrm{bFx}}\left(t\right)=-\frac{1}{2}{\dot{\phi }}_{\mathrm{s}}(t)l_{\mathrm{bF}}·\sin \left[{\phi }_{\mathrm{s}}\left(t\right)+{\phi }_{\mathrm{b}\text{\_}\mathrm{rigid}}\left(t\right)\right]\\ v{*}_{\mathrm{bFy}}\left(t\right)=\frac{1}{2}{\dot{\phi }}_{\mathrm{s}}(t)l_{\mathrm{bF}}·\cos \left[{\phi }_{\mathrm{s}}\left(t\right)+{\phi }_{\mathrm{b}\text{\_}\mathrm{rigid}}\left(t\right)\right]\end{array}\right.$$

(25)

Velocity and positional trajectories of the total center of mass of the torso in the horizontal and vertical directions can be obtained：

$$\left\{\begin{array}{l}{v}_{\mathrm{cx}\text{\_}\mathrm{flex}}\left(t\right)=\frac{m_{\mathrm{bB}}·v_{\mathrm{bBx}}\left(t\right)+m_{\mathrm{bF}}·v_{\mathrm{bFx}}\left(t\right)}{m_0}{v}_{\mathrm{cy}\text{\_}\mathrm{flex}}\left(t\right)=\frac{m_{\mathrm{bB}}·v_{\mathrm{bBy}}\left(t\right)+m_{\mathrm{bF}}·v_{\mathrm{bFy}}\left(t\right)}{m_0}\\ x_{\mathrm{c}\text{\_}\mathrm{flex}\text{\_}\mathrm{tran}}\left(t\right)=x_{\mathrm{ST}\text{\_}\mathrm{flex}}+{\displaystyle \underset{0}{\overset{t}{\int }}{v}_{\mathrm{cx}\text{\_}\mathrm{flex}}\left(t\right) }dt{y}_{\mathrm{c}\text{\_}\mathrm{flex}\text{\_}\mathrm{tran}}\left(t\right)=y_{\mathrm{ST}\text{\_}\mathrm{flex}}+{\displaystyle \underset{0}{\overset{t}{\int }}{v}_{\mathrm{cy}\text{\_}\mathrm{flex}}\left(t\right) }dt\end{array}\right.$$

(26)

Since the flexible spinal joints are in the "locked" state in the flight phase, the motion of the flexible torso center of mass in this process is the same as that of the rigid torso center of mass in the flight phase.

In summary, the parameterization of the flying phase process of the rigid torso quadruped robot is completed. As shown in Table 3.

4. Flight Phase Motion Experiment

We compare the torso center of mass trajectory, torso center of mass velocity, ground reaction force in the transition section, and motion energy during the flying phase motion of the rigid torso quadruped robot and the flexible spine torso quadruped robot, and then we analyze the effect of the flexible spine on the flying phase motion. The BQR3 quadruped robot is taken as the research object, as shown in Figure 10. The length and mass of the torso of the BQR3 robot are corresponded to the simplified model according to the treatment of front: back=3:2. The parameters of each size and mass of the BQR3 robot itself are shown in Table 4. and to determine the specific initial value of the flight phase motion process of the individual physical quantities and their initial values are shown in Table 5.

4.1. Trunk center of mass trajectory

A comparison of the trajectories of the torso center of mass motion trajectories over time for a rigid torso quadruped robot and a flexible spine quadruped robot with different spine initial joint angles at different initial velocities is shown in Figure 11. The dashed line in the figure indicates the take-off motion and the solid line indicates the flying phase motion.

From a single graph, the initial joint angle of the flexible spine affects the position of the total center of mass of the torso; the larger the absolute value of the initial angle, the lower the position of the center of mass at the take-off point. On the other hand, the larger the absolute value of the initial angle, the faster the rise of the center of mass position in the take-off motion, the higher the position of the takeoff point, and the higher the highest position of the center of mass during the whole motion.

Combining the multiple graphs, the horizontal distance moved by the center of mass in the takeoff motion does not change significantly with increasing speed for the same spine angle, but the horizontal distance moved by the center of mass in the flying phase motion increases significantly. The flexible spine can not only enhance the height of the center of mass movement, but also enhance the horizontal distance moved by the center of mass when the velocity is larger, which can help the robot to achieve a higher and farther flying phase trajectory.

4.2. Trunk center of mass Velocity

The trajectories of the horizontal and vertical velocities of the torso center of mass of the rigid torso quadruped robot compared with those of the flexible spine quadruped robot with different initial joint angles of the spine at different initial velocities over time are shown in Figure 12. From the single plot, the appropriate initial spine joint angle can improve the horizontal velocity fluctuation due to touchdown in the transition section. From the results of multiple plots, compared with the rigid torso, in the case of low speed and low stiffness of the spine joints, the appropriate initial joint angle of the spine can improve the horizontal velocity of the center of mass at the moment of takeoff. And in the case of high speed and high stiffness, the larger the initial joint angle of the spine, the greater the enhancement effect on the horizontal and vertical velocities at the moment of take-off.

4.3. Ground Reaction Force

A comparison of the ground reaction force on the rigid torso quadruped robot and the flexible spine quadruped robot with different initial joint angles of the spine at different initial velocities in the transition segment is shown in Figure 13. In order to quantify the magnitude of the ground reaction force more intuitively, the normalized ground reaction force is obtained by dividing the obtained ground reaction force with the gravity of the robot torso, and the graph is plotted with this quantity as the vertical coordinate. From a single graph, for the ground reaction force in the horizontal direction, the larger the initial joint angle of the flexible spine, the larger the peak force of the ground reaction force in both positive and negative directions. From the results of multiple plots, for the ground reaction force in the vertical direction, the ground reaction force in the vertical direction fluctuates based on the ground reaction force curve corresponding to the rigid torso as the spine joint angle increases at low speed. In the high-speed case, the impulse of the vertical ground reaction force increases with the increase of the initial spine joint angle. The same is true for the ground reaction force in the horizontal direction. Therefore, the spinal joints will have an effect on the ground reaction force. At low speeds and low stiffness, the effect on the horizontal ground reaction force is more pronounced; at high speeds and high stiffness, it has an effect on the ground reaction force in both directions, and a flexible spine will result in a greater value of the foot-ground contact force relative to a rigid spine, and the greater the angle of the spinal joint, the greater the effect.

5. Discussion and Conclusions

In this paper, we focused on the investigation of the impact of a flexible spine on the aerial phase motion of quadruped robots. We conducted dynamic modeling and motion characteristic analysis on the quadruped robot. Inspired by the aerial phase motion in biological systems, we used Fourier series to fit the ground reaction forces acting on the robot and described the parameterized aerial phase motion for both rigid and flexible spine quadruped robots. Employing a comparative analysis approach, we set different initial values for the parameterized aerial phase motion of the robot and compared the motion data between rigid and flexible spine quadruped robots. We observed differences in the robot's performance in terms of the trajectory of the trunk's center of mass, the velocity trajectory of the trunk's center of mass, and the ground reaction forces during the transition phase. We analyzed the impact of a flexible spine on the aerial phase motion of quadruped robots. The results showed that a flexible spine could increase the upward speed during the transition phase, achieve a higher vertical position at the takeoff point, and result in a higher maximum height during the aerial phase. In high-speed scenarios with high spinal stiffness, the flexible spine enabled the robot to cover a longer horizontal distance during the aerial phase. For initial joint angles of moderate size, the flexible spine improved the negative work done on horizontal velocity caused by ground reaction forces at the beginning of the touchdown, reducing the decrease in horizontal velocity and minimizing its fluctuation. The flexible spine had an impact on ground reaction forces, generally leading to larger values in the corresponding positive or negative direction. In the future, we plan to use optimization methods to consider the position and attitude trajectories of the center of mass and plan the ground reaction forces during the transition phase comprehensively. Additionally, we intend to analyze the phase from the front leg touchdown, short aerial phase, to the rear leg touchdown in a similar manner, and integrate the results of both analyses to make the analysis of the motion process more complete and the results more convincing. In the next steps of our work, we aim to validate the results using simulations or physical experiments.