Research on the biomimetic quadruped jumping robot based on an efficient energy storage structure

Abstract

The ability of quadruped robots to overcome obstacles is a critical factor that limits their practical application. Here, a design concept and a control algorithm are presented that aim at enhancing the explosive force of quadruped robots during jumping by utilizing elastic energy storage components. The hind legs of the quadruped robot are designed as energy storage units. Tension springs are utilized as components for storing energy and are installed in a parallel structure on the hind leg. Energy is stored during the compression process of the robot’s torso and released during the jumping phase. The optimal foot force is calculated using a single rigid body model. The mapping relationship between the force applied to the foot and the resulting joint torque is established by developing a dynamic model of the hind legs. Simulation experiments were conducted using the Webots physics engine to compare the impact of varying spring stiffness on joint torque during the jumping process. This study determined the optimal spring stiffness under specific conditions. The hind legs’ torque saving ratio reaches 19%, and the energy-saving ratio reaches 13%, which validates the effectiveness and feasibility of integrating elastic energy storage components.

1. Introduction

The forms and functions of robots have become diverse due to the continuous iteration of their structures and control algorithms [1–4]. The robots have become lighter, stronger, and more intelligent, which has improved their agility. In addition, the robots’ ability to adapt to the terrain was further enhanced through the use of deep reinforcement learning [5–7]. However, when a robot is confronted with a large and discontinuous terrain, it must resort to jumping to achieve spatial movement [8–14]. Imitating the behavior of animals during the jumping process can also yield certain effects [15–18]. During the jumping process, the robot’s body experiences high energy demands in a short time. The primary factors that hinder quadruped robots from successfully navigating vertical obstacles in natural environments are the need to increase joint burst torque during the jumping process and to reduce the torque output pressure.

Optimizing the trajectory to achieve the optimal joint torque for the jumping process is an attractive option. Park et al. [19] propose an optimization-based jumping method for a quadruped robot with elasticity. They plan the optimal launching area and the positions of the center of mass (CoM) and feet using Bézier curves. Zhang et al. [20] optimize the efficiency of the jump by minimizing the ground reaction forces (GRF) in the area of the CoM space and the joint torques in the joint space.

Other approaches, such as reducing the weight and improving the efficiency of the joint motor, can increase the efficiency of the jump. A low-cost, lightweight quadruped robot capable of a 41 cm vertical jump has been proposed [21]. The structure has been optimized, and materials have been changed to reduce the pressure exerted on the joint during jumping.

Additional energy storage components were incorporated to enhance joint torque during jumping. Jumping behavior is ubiquitous in the natural world, and grasshoppers are particularly notable for their ability to jump, relying on their unique energy storage mechanism to increase their power output during the jumping phase [22,23]. Based on this, promising jumping performance has been demonstrated by small jumping robots with energy storage components [24–27]. Haldane and Yim et al. [28,29] demonstrated the jumping capability of the Salto robot, which weighed 100 g and achieved a maximum jump height of approximately 1.007 m. They subsequently proposed the Salto-1P, which increased the jumping capability by approximately 25% with almost no change in weight. The jumping performance of the Salto-1P robot was further enhanced [30]. The application of elastic energy storage components in small robots is already well established. In the case of legged robots, the assistance of elastic energy storage components in jumping can effectively increase the explosive power of the legged robot launch, resulting in better jumping performance.

Although the jumping algorithm has improved the efficiency of jumping through trajectory planning, weight reduction, and increased hardware performance, there has been little optimization in terms of hardware structure. The algorithm cannot provide additional power beyond the motors for the robot’s jumping. Here, we developed a parallel-legged robot with elastic energy storage components and established the robot’s control system. To achieve higher jumping efficiency and energy utilization, we optimized the stiffness coefficient of the energy storage components, taking into account stability, and determined the optimal stiffness coefficient for the energy storage spring. Furthermore, we analyzed the energy changes of the quadruped robot during the jumping process, validating, and quantifying the energy storage effect of the elastic energy storage components.

The remainder of this article is organized as follows. Sect. 2 introduces the overall energy storage structure design and control system construction of the quadruped robot. Sect. 3 presents the dynamic model of the quadruped robot, including the modeling and parameter optimization of the elastic energy storage components. Sect. 4 describes the simulation experimental results and comparative analysis. Sect. 5 concludes and suggests potential research trajectories.

2. Structure design

2.1. Energy storage structure

The quadruped robot utilizes an energy storage structure to enhance its jumping ability. To achieve this, an elastic energy storage structure, as depicted in Fig. 1, is designed based on leg bending and energy storage, which is driven by dual motors.

Figure 1. Elastic energy storage structure. (a) Initial status. (b) Energy storage status.

The installation position is shown in Fig. 1(a). The hind legs of the robot adopt a diamond-shaped parallel structure. The length of the spring in the compressed leg changes from its original length when the robot transitions from free-standing to the energy storage state, as shown in Fig. 1(b). The initial spring length is maintained when the robot’s legs are in the initial state. Replace the original leg structure of the standard quadruped robot with the structure in Fig. 1.

Figure 2. Simplified model of the hind legs. (a) The left hind leg consists of four rods, with lengths a, b, c, and d, respectively (a = b = c = d^′). The joint between rod 1 and the torso is called joint 1, and the joint between rod 2 and the torso is called joint 2. Similarly, the joints for the right hind leg are joint 3 and joint 4. (b) The left front leg consists of two rods with lengths e, f. The joint between rod 5 and the torso is called joint 5, and the joint between rod 5 and rod 6 is called joint 6. Similarly, the joints for the right hind leg are joint 7 and joint 8.

2.2. Dynamics

The dynamic equations are established using the general form of the Lagrange dynamic equations, as the rear leg has a parallel leg structure.

(1) \begin{equation} \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_{j}}\right)-\frac{\partial L}{\partial q_{j}}=Q_{j} \end{equation}

As shown in Fig. 2(a), a coordinate system is established using the left leg as an example, with the hip joint as the origin. The center of mass of each rod should be considered at the center of its form. The velocities of the center of mass of the four rods are:

(2) \begin{equation} \mathrm{rod}\ 1: \left[\begin{array}{c} \dot{x}_{O1}\\[4pt] \dot{z}_{O1} \end{array}\right]=\left[\begin{array}{c} \dfrac{1}{2}a\cos \!\left(\theta _{1}\right)\\[9pt] \dfrac{1}{2}a\sin \!\left(\theta _{1}\right) \end{array}\right]\dot{\theta }_{1} \end{equation}

(3) \begin{equation} \mathrm{rod}\ 2:\left[\begin{array}{c} \dot{x}_{O2}\\[4pt] \dot{z}_{O2} \end{array}\right]=\left[\begin{array}{c} \dfrac{1}{2}b\cos \!\left(\theta _{2}\right)\\[9pt] \dfrac{1}{2}b\sin \!\left(\theta _{2}\right) \end{array}\right]\dot{\theta }_{2} \end{equation}

(4) \begin{equation} \mathrm{rod}\ 3:\left[\begin{array}{c} \dot{x}_{O3}\\[4pt] \dot{z}_{O3} \end{array}\right]=\left[\begin{array}{c@{\quad}c} b\cos \!\left(\theta _{2}\right) & \dfrac{1}{2}c\cos \!\left(\theta _{1}\right)\\[9pt] b\sin \!\left(\theta _{2}\right) & \dfrac{1}{2}c\sin \!\left(\theta _{1}\right) \end{array}\right]\left[\begin{array}{c} \dot{\theta }_{2}\\[4pt] \dot{\theta }_{1} \end{array}\right] \end{equation}

(5) \begin{equation} \mathrm{rod}\ 4:\left[\begin{array}{c} \dot{x}_{O4}\\[4pt] \dot{z}_{O4} \end{array}\right]=\left[\begin{array}{c@{\quad}c} a\cos \!\left(\theta _{1}\right) & \dfrac{1}{2}d\cos \!\left(\theta _{2}\right)\\[9pt] a\sin \!\left(\theta _{1}\right) & \dfrac{1}{2}d\sin \!\left(\theta _{2}\right) \end{array}\right]\left[\begin{array}{c} \dot{\theta }_{1}\\[4pt] \dot{\theta }_{2} \end{array}\right] \end{equation}

where $\dot{x}_{Oi}$ (i = 1 to 4) and $\dot{z}_{Oi}$ are the velocities along the x-axis and z-axis. The four rear leg rods have lengths a, b, c, and d, respectively. $\theta _{1}$ represents the angle between the rod 1 and the vertical direction, typically with a positive value. Similarly, $\theta _{2}$ represents the angle between rod 2 and the vertical direction, typically with a negative value. The system’s total translational kinetic energy is

(6) \begin{equation} K_{1}=\sum _{i=1}^{4}\frac{1}{2}m_{i}\!\left(\dot{x}_{oi}^{2}+\dot{z}_{oi}^{2}\right) \end{equation}

The angular velocity of the four rods $\omega _{Oi}$ can be expressed as

(7) \begin{equation} \mathrm{rod}\ 1:\ \omega _{O1}=\dot{\theta }_{1} \end{equation}

(8) \begin{equation} \mathrm{rod}\ 2:\ \omega _{O2}=\dot{\theta }_{2} \end{equation}

(9) \begin{equation} \mathrm{rod}\ 3:\ \omega _{O3}=\dot{\theta }_{1} \end{equation}

(10) \begin{equation} \mathrm{rod}\ 4:\ \omega _{O4}=\dot{\theta }_{2} \end{equation}

The system’s total rotational kinetic energy is

(11) \begin{equation} K_{2}=\sum _{i=1}^{4}\frac{1}{2}I_{i}\omega _{Oi}^{2} \end{equation}

where $I_{i}$ is the moment of inertia of the $i$ -th rod about the y-axis with respect to the Centroid frame, the Centroid frame represents a coordinate system with the robot’s center of mass as the origin, where the y-axis aligns with the {H} coordinate system, and the x-axis follows a right-handed orientation along the rod. $m_{i}$ is the mass of the $i$ -th rod, and $\omega _{Oi}$ is the angular velocity of the $i$ -th rod with respect to the frame{H}.

The gravitational potential energy of the system is calculated with respect to the z = 0 plane of the frame{H}, which is taken as the point of zero potential energy.

(12) \begin{equation} P_{1}=-\frac{1}{2}g\!\left(m_{1}a\cos \!\left(\theta _{1}\right)+m_{2}b\cos \!\left(\theta _{2}\right)+m_{3}\!\left(c\cos \!\left(\theta _{1}\right)+2b\cos \!\left(\theta _{2}\right)\right)+m_{4}\!\left(d\cos \!\left(\theta _{2}\right)+2a\cos \!\left(\theta _{1}\right)\right)\right) \end{equation}

Set the coefficient of elasticity of the spring to be $k$ , and establish the relationship between the spring length and the angle of the joint.

(13) \begin{equation} l=\left(a^{2}+b^{2}-2ab\cos \!\left(\theta _{1}-\theta _{2}\right)\right)^{\frac{1}{2}} \end{equation}

The elastic potential energy of a spring in a system can be calculated using the following formula.

(14) \begin{equation} P_{2}=\frac{1}{2}k\!\left(l-l_{0}\right)^{2} \end{equation}

where l is the length of the spring, and l ₀ is the original length of the spring. The Lagrangian equations with a nonpotential force are expressed as

(15) \begin{equation} Q_{j}=\sum _{i=1}^{n}\boldsymbol{F}_{i}\cdot \frac{\partial \boldsymbol{r}_{i}}{\partial q_{j}} \end{equation}

The GRF is transformed into a nonpotential force.

(16) \begin{equation} \boldsymbol{Q}=\left[\begin{array}{c} a\cos \!\left(\theta _{1}\right)F_{x} + a\sin \!\left(\theta _{1}\right)F_{z}\\[4pt] d\cos \!\left(\theta _{2}\right)F_{x} + d\cos \!\left(\theta _{2}\right)F_{z} \end{array}\right] \end{equation}

where F _x and F _y are the horizontal and vertical components of the foot force calculated from the single rigid body model, respectively.

By substituting Equations (6), (11)–(15) into the general form of the Lagrange equation, we can obtain

(17) \begin{equation} \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_{1}}\right);-\frac{\partial L}{\partial q_{1}}=a\cos \!\left(\theta _{1}\right)F_{x}+a\sin \!\left(\theta _{1}\right)F_{z}+\tau _{1} \end{equation}

(18) \begin{equation} \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_{2}}\right);-\frac{\partial L}{\partial q_{2}}=d\cos \!\left(\theta _{2}\right)F_{x}+d\cos \!\left(\theta _{2}\right)F_{z}+\tau _{2} \end{equation}

The solution provides the joint torques

(19) \begin{equation} \boldsymbol{\tau }=\boldsymbol{M}\!\left(\theta \right)\ddot{\theta }+\boldsymbol{c}\!\left(\theta,\dot{\theta }\right)+\boldsymbol{g}\!\left(\theta \right)+\boldsymbol{k}\!\left(\theta \right)+\boldsymbol{f}_{ext} \end{equation}

where

(20) \begin{equation} \boldsymbol{\tau }\boldsymbol{=}\left[\begin{array}{c} \tau _{1}\\[4pt] \tau _{2} \end{array}\right] \end{equation}

(21) \begin{equation} \boldsymbol{M}\!\left(\unicode{x03B8} \right)=\left[\begin{array}{c@{\quad}c} \left(m_{\mathrm{4}}+\dfrac{m_{{1}}}{4}\right)a^{\mathrm{2}}+\dfrac{m_{\mathrm{3}}c^{\mathrm{2}}}{4}+\mathrm{I}_{{1}}+\mathrm{I}_{\mathrm{3}} & {0}.\mathrm{5}\!\left(adm_{\mathrm{4}}+bcm_{\mathrm{3}}\right)\cos \!\left(\theta _{{1}}-\theta _{\mathrm{2}}\right)\\[8pt] {0}.\mathrm{5}\!\left(adm_{\mathrm{4}}+bcm_{\mathrm{3}}\right)\cos \!\left(\theta _{{1}}-\theta _{\mathrm{2}}\right) & \left(m_{\mathrm{3}}+\dfrac{m_{\mathrm{2}}}{4}\right)b^{\mathrm{2}}+\dfrac{m_{\mathrm{3}}c^{\mathrm{2}}}{4}+\mathrm{I}_{\mathrm{2}}+\mathrm{I}_{\mathrm{4}} \end{array}\right] \end{equation}

(22) \begin{equation} \boldsymbol{c}\!\left(\theta,\dot{\theta }\right)=\left[\begin{array}{c} \dfrac{1}{2}\left(adm_{4}+bcm_{3}\right)\dot{\theta }_{2}^{2}\sin \!\left(\theta _{1}-\theta _{2}\right)\\[9pt] -\dfrac{1}{2}\left(adm_{4}+bcm_{3}\right)\dot{\theta }_{1}^{2}\sin \!\left(\theta _{1}-\theta _{2}\right) \end{array}\right] \end{equation}

(23) \begin{equation} \boldsymbol{g}\!\left(\theta \right)=\left[\begin{array}{c} \dfrac{1}{2}\sin \!\left(\theta _{1}\right)\left(am_{1}+2am_{4}+cm_{3}\right)g\\[9pt] \dfrac{1}{2}\sin \!\left(\theta _{2}\right)\left(bm_{2}+2bm_{3}+dm_{4}\right)g \end{array}\right] \end{equation}

(24) \begin{equation} \boldsymbol{k}\!\left(\theta \right)=\left[\begin{array}{c} kab\sin \!\left(\theta _{1}-\theta _{2}\right)-kl_{0}ab\sin \!\left(\theta _{1}-\theta _{2}\right)\left(a^{2}+b^{2}-2ab\cos \!\left(\theta _{1}-\theta _{2}\right)\right)^{-\frac{1}{2}}\\[4pt] -kab\sin \!\left(\theta _{1}-\theta _{2}\right)+kl_{0}ab\sin \!\left(\theta _{1}-\theta _{2}\right)\left(a^{2}+b^{2}-2ab\cos \!\left(\theta _{1}-\theta _{2}\right)\right)^{-\frac{1}{2}} \end{array}\right] \end{equation}

(25) \begin{equation} \boldsymbol{f}_{ext}=\left[\begin{array}{c} -a\cos \!\left(\theta _{1}\right)F_{x}-a\sin \!\left(\theta _{1}\right)F_{z}\\[4pt] -d\cos \!\left(\theta _{2}\right)F_{x}-d\cos \!\left(\theta _{2}\right)F_{z} \end{array}\right] \end{equation}

2.3. Overall structure

The quadruped robot has a total of eight degrees of freedom. Each leg has two degrees of freedom. The front legs have a series of segments with a concave elbow design, while the hind legs have a parallel leg structure. The hind legs have a larger number of links, which shifts the overall center of mass rearward. As a result, the front legs require relatively lower torque. The quadruped robot’s front legs use a synchronous belt drive, while the hind legs utilize energy storage springs and a pulley drive with a reduction ratio to increase their maximum torque output. Figure 3 illustrates the overall structure of the robot.

Figure 3. The overall structure of a quadruped robot.

2.4. Jump controller

This paper focuses on the role of the energy storage element in the jumping process and briefly describes the control framework depicted in Fig. 4. The motion scheduler receives input velocity instructions, schedules joint positions and velocities, and determines the contact time of each foot. The model predictive control utilizes a three-dimensional single rigid body model to calculate ground reaction forces, and the simplified model enhances computational efficiency. The leg dynamics model calculates feedforward torques for the joints, which are then combined with PD (Proportional-Differential) control to achieve motor control for the joints. The orientation of the body is estimated using readings from the IMU (Inertial Measurement Unit) gyroscope and accelerometer. Additionally, the velocity and position of the body are estimated using a linear Kalman filter.

Figure 4. The control framework used in this work.

Figure 5. Schematic of the quadruped robot.

3. Dynamics model

After determining the structure of the quadruped robot, the single rigid body dynamics model of the quadruped robot was established using model predictive control to generate the corresponding GRF. The mapping torque of each joint was calculated from the kinetic models of the front and back legs.

3.1. Model predictive control

Figure 5 shows the global coordinate {W} and body coordinate {B}. The translational and rotational acceleration of the rigid body are obtained using Newton’s Euler equations [31]

(26) \begin{equation} {}^{W}{\ddot{\boldsymbol{p}}}{_{com}^{}}=\frac{\sum _{n=0}^{3}{}^{W}{\boldsymbol{f}}{_{n}^{}}}{m_{\textit{total}}}+{}^{W}{\boldsymbol{g}}{} \end{equation}

(27) \begin{equation} \frac{d}{dt}\left(\boldsymbol{I}\cdot {}^{W}{\boldsymbol{\omega }}{_{WB}^{}}\right)=\sum _{i=0}^{3}{}^{W}{\boldsymbol{r}}{_{i}^{}}\times {}^{W}{\boldsymbol{f}}{_{i}^{}}+\mathbf{0}_{3\times 1}\times \left(m_{\textit{total}}\cdot {}^{W}{\boldsymbol{g}}{}\right) \end{equation}

where ${}^{W}{\ddot{\boldsymbol{p}}}{_{com}^{}}$ and ${}^{W}{\boldsymbol{f}}{_{n}^{}}$ (n = 0 to 3) are three-dimensional vectors representing the robot’s acceleration and GRF with respect to the global frame. ^W g is the gravitational acceleration. $\boldsymbol{I}\in \boldsymbol{R}^{3\times 3}$ is the rotational inertia tensor and ${}^{W}{\boldsymbol{\omega }}{_{WB}^{}}$ is the angular velocity of the body with respect to the global frame. ${}^{W}{\boldsymbol{r}}{_{i}^{}}$ is the position of $i$ -th contact point with respect to the CoM of the robot. $m_{\textit{total}}$ is the robot’s body mass.

The conversion of the moment of inertia of rotation from the body frame to the oriented body frame {P} is as follows.

(28) \begin{equation} {}^{P}{\boldsymbol{I}}{}={}^{W}{\boldsymbol{R}}{_{B}^{}} \bullet {}^{B}{\boldsymbol{I}}{}\bullet {}^{W}{\boldsymbol{R}}{_{B}^{T}} \end{equation}

The oriented body frame aligns with the world coordinate system in direction, with the origin coinciding with the body coordinate system. $^{W}\boldsymbol{R}_{B}$ is a rotation matrix. Set the system state variables $x$ :

(29) \begin{equation} \boldsymbol{x}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \boldsymbol{\varTheta }^{T} & {}^{W}{{\boldsymbol{p}}}{_{com}^{T}} & {}^{W}{\boldsymbol{\omega }}{_{WB}^{T}} & {}^{W}{\dot{\boldsymbol{p}}}{_{com}^{T}} & {}^{W}{\boldsymbol{g}}{} \end{array}\right]^{T} \end{equation}

where $\boldsymbol{\varTheta }$ is the Euler angle to represent body orientation. ${}^{W}{\boldsymbol{p}}{_{com}^{}}$ and ${}^{W}{\dot{\boldsymbol{p}}}{_{com}^{}}$ are the robot’s position and robot’s velocity with respect to the global frame. We have made some assumptions because the jump process is a two-dimensional plane. Since the yaw and roll angles vary very little, we can assume that the differential of Euler angles maps to angular velocity as follows:

(30) \begin{equation} \left[\begin{array}{c} \begin{array}{c} \dot{\phi }\\[4pt] \dot{\Phi } \end{array}\\[4pt] \dot{\psi } \end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & \tan \boldsymbol{\varPhi} \\[4pt] 0 & 1 & 0\\[4pt] 1 & 0 & \cos \boldsymbol{\varPhi} ^{-1} \end{array}\right]{}^{O}{\omega }{_{OB}^{}}=\boldsymbol{R}_{y}^{T}\cdot {}^{O}{\boldsymbol{\omega }}{_{OB}^{}} \end{equation}

where $\dot{\phi }, \dot{\Phi }$ and $\dot{\psi }$ represent the yaw, pitch, and roll of the body orientation. The dynamics of the system can be expressed as

(31) \begin{equation} \dot{\boldsymbol{x}}=\boldsymbol{Ax}+\boldsymbol{Bu} \end{equation}

where

(32) \begin{equation} \boldsymbol{A}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \boldsymbol{R}_{y}^{\boldsymbol{T}}\left(\boldsymbol{\varPhi} \right) & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}\\[4pt] \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathrm{II}_{\mathrm{3}} & \mathbf{0}\\[4pt] \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}\\[4pt] \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & {1}\\[4pt] \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right] \end{equation}

(33) \begin{equation} \boldsymbol{B}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}}\\[4pt] \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}}\\[4pt] {}^{P}{\mathrm{I}}{_{}^{-{1}}}\boldsymbol{r}_{0\times } & {}^{P}{\mathrm{I}}{_{}^{-{1}}}\boldsymbol{r}_{1\times } & {}^{P}{\mathrm{I}}{_{}^{-{1}}}\boldsymbol{r}_{2\times } & {}^{P}{\mathrm{I}}{_{}^{-{1}}}\boldsymbol{r}_{3\times }\\[4pt] \mathrm{II}_{3}/\mathrm{m}_{\mathrm{total}} & \mathrm{II}_{3}/\mathrm{m}_{\mathrm{total}} & \mathrm{II}_{3}/\mathrm{m}_{\mathrm{total}} & \mathrm{II}_{3}/\mathrm{m}_{\mathrm{total}}\\[4pt] \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} & \mathbf{0}_{\mathrm{3}\times \mathrm{3}} \end{array}\right] \end{equation}

(34) \begin{equation} \boldsymbol{u}=\left[\begin{array}{c} {}^{W}{\boldsymbol{f}}{_{0}^{}}\\[4pt] {}^{W}{\boldsymbol{f}}{_{1}^{}}\\[4pt] {}^{W}{\boldsymbol{f}}{_{2}^{}}\\[4pt] {}^{W}{\boldsymbol{f}}{_{3}^{}} \end{array}\right] \end{equation}

where $\Pi$ ₃∈R³ is the identity matrix.

Discretize the process:

(35) \begin{equation} \boldsymbol{x}\!\left(k+{1}\right)=\boldsymbol{\Pi}+\Delta T\times \boldsymbol{A}\!\left(\boldsymbol\varPhi \right))\boldsymbol{x}(k)+\Delta T\times \boldsymbol{Bu}(k)=\boldsymbol{A}_{k}\boldsymbol{x}(k)+\boldsymbol{B}_{k}\boldsymbol{u}(k) \end{equation}

Using a discretized equation of state, the GRF is predicted for the subsequent n cycles of the jumping phase:

(36) \begin{align} &\boldsymbol{x}(1) =\boldsymbol{A}_{0}\boldsymbol{x}(0)+\boldsymbol{B}_{0}\boldsymbol{u}(0)\nonumber\\[4pt]& \boldsymbol{x}(2) =\boldsymbol{A}_{1}\boldsymbol{A}_{0}\boldsymbol{x}(0)+\boldsymbol{A}_{1}\boldsymbol{B}_{0}\boldsymbol{u}(0)+\boldsymbol{B}_{1}\boldsymbol{u}(1)\nonumber\\[4pt]& \boldsymbol{x}(3) =\boldsymbol{A}_{2}\boldsymbol{A}_{1}\boldsymbol{A}_{0}\boldsymbol{x}(0)+\boldsymbol{A}_{2}\boldsymbol{A}_{1}\boldsymbol{B}_{0}\boldsymbol{u}(0)+\boldsymbol{A}_{2}\boldsymbol{B}_{1}\boldsymbol{u}(1)+\boldsymbol{B}_{2}\boldsymbol{u}(2)\nonumber\\[4pt]& \vdots \nonumber\\[4pt]& \boldsymbol{x}(h) =\left(\prod\nolimits _{i=n-1}^{0}\boldsymbol{A}_{i}\right)\boldsymbol{x}(0)+\sum _{i=0}^{n-2}\left[\left(\prod\nolimits _{j=n-1}^{i+1}A_{i}\right)\boldsymbol{B}_{i}\boldsymbol{u}(i)\right]+\boldsymbol{B}_{n-1}\boldsymbol{u}\left(n-1\right) \end{align}

We construct a QP that minimizes

(37) \begin{equation} \min J\!\left(\boldsymbol{U}\right)=\left(\boldsymbol{X}-\boldsymbol{D}\right)^{T}\boldsymbol{Q}\!\left(\boldsymbol{X}-\boldsymbol{D}\right)+\boldsymbol{U}^{T}\boldsymbol{RU} \end{equation}

where

(38) \begin{equation} \boldsymbol{X}=\boldsymbol{A}_{qp}\boldsymbol{x}(0)+\boldsymbol{B}_{qp}\boldsymbol{U} \end{equation}

(39) \begin{equation} \boldsymbol{U}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \boldsymbol{u}(0) & \boldsymbol{u}(1) & \cdots & \boldsymbol{u}\!\left(n-1\right) \end{array}\right] \end{equation}

$\boldsymbol{D}$ is the desired state of the system. The following constraints are presented below. The limits of the friction cone are as follows:

(40) \begin{equation} \underline{\boldsymbol{c}}_{i}\leq \boldsymbol{C}_{i}\cdot {}^{W}{\boldsymbol{f}}{_{c}^{i}}\leq \overline{\boldsymbol{c}}_{i} \end{equation}

S _i denotes the contact of leg i. The contacted constraints are expressed as

(41) \begin{equation} {}^{W}{\boldsymbol{f}}{_{i}^{}}=0_{3\times 1,}\forall S_{n}=0 \end{equation}

Quadratic programming can be utilized to calculate the optimal sequence of foot end forces. The first set of foot end forces that are optimized to produce this sequence should be used as the desired foot end forces.

(42) \begin{equation} {}^{W}{\boldsymbol{f}}{_{}^{MPC}}=\boldsymbol{U}(0) \end{equation}

The method for joint control is as follows:

(43) \begin{equation} \tau =\tau _{ff}+k_{p}\left(q_{ref}-q\right)+k_{d}\left(\dot{q}_{ref}-\dot{q}\right) \end{equation}

where $\tau _{ff}$ is the feedforward torque and $k_{p}$ and $k_{d}$ represent the position stiffness and the velocity stiffness. The simulation controls the motor using five parameters: $\tau _{ff}, k_{p}, k_{d}, q_{ref}, \dot{q}_{ref}$ .

The command for the acceleration of the center of mass of a quadruped robot involves the stiffness coefficient and damping coefficient. This command is given by the following equations:

(44) \begin{equation} {}^{W}{\ddot{\boldsymbol{p}}}{_{COM}^{}}=\boldsymbol{K}_{p}\cdot {}^{W}{\boldsymbol{p}}{_{COM}^{}}+\boldsymbol{D}_{p}\cdot {}^{W}{\dot{\boldsymbol{p}}}{_{COM}^{}} \end{equation}

(45) \begin{equation} {}^{W}{\dot{\boldsymbol{\omega }}}{_{WB}^{}}=\boldsymbol{K}_{\omega }\hat{\boldsymbol{\omega }}\theta +\boldsymbol{D}_{\omega }{}^{W}{\boldsymbol{\omega }}{_{WB}^{}} \end{equation}

where $\hat{\boldsymbol{\omega }}\theta$ is the exponential coordinate of the rotation matrix.

3.2. Optimization of the spring coefficient

The spring coefficient k of the rear leg joint has a significant impact on the joint torque. A larger spring constant results in greater energy being stored in the system, but it also demands more torque from the motor for energy storage. During the energy storage phase, the reverse torque required by the motor increases with the spring constant k. It is important to ensure that the reverse torque of a single motor does not exceed the rated torque. When the energy storage state of the system reaches a predetermined maximum state, the maximum torque of the motor should not exceed the maximum torque of the motor without springs.

In the energy storage state, the elongation of the spring is determined without considering the velocities and accelerations of the system’s various components.

(46) \begin{equation} \Delta l=\left(a^{2}+b^{2}-2ab\cos \!\left(\theta _{1}-\theta _{2}\right)\right)^{\frac{1}{2}}-l_{0} \end{equation}

The dynamic model of the hind leg is used to determine the maximum torque. It represents the state of maximum energy storage during the energy storage phase.

(47) \begin{equation} \boldsymbol{\tau }=\boldsymbol{g}\!\left(\theta \right)+\boldsymbol{k}\!\left(\theta \right)+\boldsymbol{f}_{ext} \end{equation}

The variables Fx and Fy represent the horizontal and vertical components of the force on the foot end of the back leg when the robot is at rest. Additionally, $\theta$ represents the joint angles after the energy storage in the energy storage phase. All moments must satisfy the following constraints.

(48) \begin{equation} \boldsymbol{\tau }_{1}\prec \tau _{1-std\max } \end{equation}

(49) \begin{equation} \boldsymbol{\tau }_{2}\prec \tau _{2-std\max } \end{equation}

(50) \begin{equation} \boldsymbol{\tau }_{1}\prec T_{\max } \end{equation}

(51) \begin{equation} \boldsymbol{\tau }_{2}\prec T_{\max } \end{equation}

where $\tau _{1-std\max }$ and $\tau _{2-std\max }$ represent the maximum torque of joint 1 and joint 2 during the storage energy and jumping process without the elastic energy storage component (k = 0). $T_{\max }$ stands for the peak torque of the joint motor.

The Equations (48–51) show that $\tau _{1}$ and $\tau _{2}$ should be less than the maximum moment of the quadruped robot during the process of completing the predetermined jumping target velocity, including the energy storage phase and the jumping phase.

To determine the optimal spring coefficient, in the case where the quadrupedal robot is in the state of maximum stored energy and the constraints of Equations (48) and (49) are satisfied, and the relevant parameters are brought into the Equaions (47). We then solved for the optimal spring coefficient.

3.3. Analysis of system energy efficiency

After determining the optimal spring coefficient for the motor, the hardware parameters for the entire system can be established. Next, an energy analysis of the entire quadruped robot system is conducted.

During the energy storage phase, the quadruped robot transitions from a freely standing state to the posture where energy is stored. The velocity of the CoM of the robot is very slow during this process, and the total energy of the system does not account for the robot’s kinetic energy. The ground where the robot is located is taken as the zero potential energy point. The system’s total energy is given as

(52) \begin{align} E_{1\_ k} & =m_{\text{total}}gh_{1}+2\times \frac{1}{2}k_{op}\!\left(l-l_{0}\right)^{2}\nonumber\\[4pt] & =m_{\textit{total}}gh_{1}+k_{op}\!\left(l-l_{0}\right)^{2} \end{align}

where h ₁ is the height of the center of mass during the energy storage phase. The energy output of the motor is :

(53) \begin{equation} W_{{1}\_ \mathrm{k}}=\sum _{j={1}}^{\mathrm{8}}\left(\int \left| \tau _{j}\right| d\theta \right) \end{equation}

where τ _j denotes the moment of the j-th joint. During the jumping phase, the quadruped robot initiates a jump from a posture where energy is stored, resulting in both feet leaving the ground. After this phase, the robot’s total energy can be calculated.

(54) \begin{equation} E_{2\_ k}=m_{\textit{total}}gh_{2}+\frac{1}{2}m_{\textit{total}}\nu ^{2}+k\!\left(l-l_{0}\right)^{2} \end{equation}

where h ₂ is the height of the center of mass during the jumping phase. The energy output of the motor is as follows:

(55) \begin{equation} W_{2\_ k}=\sum _{j=1}^{8}\left(\int \!\left| \tau _{j}\right| d\theta \right) \end{equation}

Similarly, the energy for each phase of the system can be calculated without the presence of a stored energy spring ( $k=0$ ). Details are given below.

(56) \begin{equation} E_{1\_ std}=m_{\textit{total}}gh_{1} \end{equation}

(57) \begin{equation} W_{1\_ std}=\sum _{j=1}^{8}\left(\int \left| \tau _{j}\right| d\theta \right) \end{equation}

(58) \begin{equation} E_{2\_ std}=m_{\text{total}}gh_{2}+\frac{1}{2}m_{\textit{total}}\nu ^{2} \end{equation}

(59) \begin{equation} W_{2\_ std}=\sum _{j=1}^{8}\left(\int \left| \tau _{j}\right| d\theta \right) \end{equation}

To describe torque optimization more objectively, we use the Equation (60) to define the torque saving factor.

(60) \begin{equation} \eta =0.5\times \frac{|\boldsymbol{\tau }_{1-\textit{stdmax}}-\boldsymbol{\tau }_{1-\textit{sprmax}}|}{|\boldsymbol{\tau }_{1-\textit{stdmax}}|}+0.5\times \frac{|\boldsymbol{\tau }_{2-\textit{stdmax}}-\boldsymbol{\tau }_{2-\textit{sprmax}}|}{|\boldsymbol{\tau }_{2-\textit{stdmax}}|} \end{equation}

The maximum torque of the two motors in the back leg (k = 600) is $\tau _{1-\textit{sprmax}}$ and $\tau _{2-\textit{sprmax}}$ . Additionally, the energy-saving factor can be determined as

(61) \begin{equation} \xi =\frac{W_{1\_ std}+W_{2\_ std}-W_{1\_ spr}-W_{2\_ spr}}{W_{1\_ std}+W_{2\_ std}} \end{equation}

The motor’s energy outputs at standard state ( $k=0$ ) are represented by $W_{1\_ std}$ and $W_{2\_ std}$ , while they are represented by $W_{1\_ spr}$ and $W_{2\_ spr}$ at the added spring state ( $k=k_{op}$ ). This allows the calculation of the robot’s energy change over time.

4. Experiment and analysis

4.1. Parameter determination

The simulation experiment for this robot was conducted on the Webots platform using a self-designed quadruped robot model. The simulation parameters are detailed in Table I. The stiffness coefficient and damping coefficient of centroid position, angle, and eight motors are detailed in Table II.

Table I. Parameters of the robot.

Table II. Parameters of K and D.

4.2. Jumping process

The quadruped robot stores energy through downward pressure, as shown in Fig. 6. Once the center of mass reaches the predetermined position, it initiates the jumping phase. The quadruped robot undergoes a power accumulation phase (0–300 ms), followed by a takeoff phase (300 ms–400 ms), and finally an air phase (400 ms–720 ms).

Figure 6. Quadrupedal robot jumping process.

Figure 7. Joint torque without energy storage components.

4.3. Results

To compare the subsequent addition of an elastic energy storage component, we first tested the quadruped robot without the energy storage component (i.e. k = 0) during both the downward pressure and jumping phases. To reduce the impact of jumping velocity on the elastic energy stored during jumping, we used a target velocity of 2.5 m/s vertically in this study. Figure 7 illustrates the variation in joint torque during the takeoff phase for the hind limb joints. The feedforward torques, τ_1,ff and τ_2,ff, were calculated based on the dynamic model of the energy storage structure mentioned earlier for joints 1 to 2, while τ₁ and τ₂ represent the actual torque changes for the corresponding joints. The same applies to joints 3 and 4. The graph indicates that each joint exhibits a peak in output torque at approximately 310 ms. Figure 8 displays the angular variation of the hind limb-driven joints during the takeoff phase. The desired joint angles for joints 1 to 4 are represented by θ _1,des to θ _4,des, while the actual joint angles are represented by θ ₁ to θ ₄. Both graphs demonstrate that the error between the actual joint torques and angles is less than 10%, thus validating the correctness of the hind limb dynamics.

Figure 8. Joint angle.

Figure 9 shows the variation of joint torque over time under different elastic coefficients at the same target velocity. It has been observed that within a certain range, increasing the elastic coefficient k results in a relatively smaller overall torque sequence for the hind leg motors, without considering the motor torque limit. However, a large elastic coefficient leads to significant torque fluctuations, indicating an unstable system state.

Figure 9. Joint torques with different elastic coefficients. (l₀= 206 mm).

Figure 10. Torque under maximum energy storage.

The relationship between joint motor torque and various elastic coefficients is derived by transforming Equation (47), as shown in Fig. 10. τ _1,ff and τ _2,ff represent the relationship between the actual feedforward torques of joint 1 and joint 2 under maximum energy storage conditions and the stiffness coefficient of the spring. As the stiffness coefficient increases, the required torque to drive the joint decreases, as shown in Fig 9. To ensure that the maximum torque of the motor is smaller than the maximum torque without the energy storage structure installed, the optimal elasticity coefficient is calculated as

(62) \begin{equation} \left[\begin{array}{c} k_{op}\\[4pt] k_{op} \end{array}\right]=\left(\left[\begin{array}{c} \tau _{1-std\max }\\[4pt] \tau _{2-std\max } \end{array}\right]-\boldsymbol{f}_{ext}-\boldsymbol{g}\!\left(\theta \right)\right)\left(\boldsymbol{k}\!\left(\theta \right)/k\right)^{-1} \end{equation}

According to Fig. 7, the maximum torques of the motor are 15 Nm and –17 Nm, respectively. To ensure system stability, minimize spring torque fluctuation, and maintain a sufficient safety margin, the optimal elastic coefficient is set at 600 while meeting this condition. Figure 11 shows the variation in motor joint torque from the energy storage phase to the jumping phase. During the transition from the energy storage phase to the takeoff phase, there is a sudden surge in motor torque. Figure 12 depicts the variation in the robot’s CoM height and velocity. The energy storage phase is represented by the period from 0 to 300 ms, during which the robot’s center of mass height decreases from 0 ms and reaches the predetermined height by 300 ms, marking the end of the energy storage phase. At 301 ms, the robot’s joint torque sharply increases, initiating the jumping phase. The jumping phase concludes at 400 ms when the robot reaches the predetermined center of mass velocity, signifying the end of the jumping phase. Figure 13 illustrates the variation in spring length. The spring length gradually increases during the energy storage phase and reaches its maximum length upon completion of this phase. During the jumping phase, the length of the spring decreases rapidly, providing an additional explosive force for the jumping process of the robot.

Figure 11. Torque with optimal elastic coefficients.

Figure 12. State of center of mass of the z-axis ( $k=k_{op}$ ).

Figure 13. Spring length ( $k=k_{op}$ ).

Equations (48)–(55) can be used to calculate the energy variation of the robot at different times. Figure 14 illustrates the total energy variation of the quadruped robot without elastic energy storage. The change in gravitation potential energy while crouching is stored by the springs as elastic potential energy, and this energy in turn is used in the jumping phase. This increases additional output power.

Figure 14. The total energy of the robot.

Equation (60) shows that at the target velocity, the torque reduction factors for the left and right hind legs of the robot are 0.1821 and 0.1978, respectively. The addition of elastic energy storage components enables the robot to significantly reduce motor torque output pressure during energy storage and jumping phases. Calculations based on Equation (61) demonstrate that the energy storage components enhance the energy utilization of the robot by collecting the reduced gravitational potential energy during the energy storage phase. Also, this validates the feasibility and efficiency of installing energy storage components by reducing the robot’s energy output by approximately 13%.

5. Conclusion

Here we present a storage structure and jumping algorithm for quadruped robots based on elastic energy storage components. The storage structure adopts a parallel mechanism with added energy storage elements, completing the construction of a quadruped robot with a storage structure by replacing the original hind leg structure. The jumping algorithm is based on model predictive control of a single rigid body, to optimize the elastic coefficient. This reduces the torque required by the hind leg joints by 19% while meeting the same jumping requirements, resulting in a saving of approximately 13% of energy output. The study confirms that the robot system can collect and release gravitational potential energy by adding elastic energy storage components before the jumping phase. This leads to an additional burst of torque during the jumping phase through energy release. The feasibility and effectiveness of the storage structure under the balanced jumping algorithm framework were verified through calculation and simulation.

This study investigates the use of elastic energy storage devices in high-performance jumping of quadruped robots. The research demonstrates that this strategy can offer supplementary energy in extreme jumping scenarios, enhance energy utilization during the jumping process, and increase the maximum jumping height under the same torque conditions. Based on prior research, this study presents a practical approach for future scenarios involving modular elastic energy storage components.