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.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 [Reference Kim, Di Carlo, Katz, Bledt and Kim31]

(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.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%.