Estimation of leg stiffness using an approximation to the planar spring–mass system in high-speed running

Abstract

Leg stiffness plays a critical role in legged robots’ speed regulation. However, the analytic solutions to the differential equations of the stance phase do not exist, of course not for the exact analytical solution of stiffness. In view of the challenge in dealing with every circumstance by numerical methods, which have been adopted to tabulate approximate answers, the “harmonic motion model” was used as approximation of the stance phase. However, the wide range leg sweep angles and small fluctuations of the “center of mass” in fast movement were overlooked. In this article, we raise a “triangle motion model” with uniform forward speed, symmetric movement, and straight-line center of mass trajectory. The characters are then shifted to a quadratic equation by Taylor expansion and obtain an approximate analytical solution. Both the numerical simulation and ADAMS-Matlab co-simulation of the control system show the accuracy of the triangle motion model method in predicting leg stiffness even in the ultra-high-speed case, and it is also adaptable to low-speed cases. The study illuminates the relationship between leg stiffness and speed, and the approximation model of the planar spring–mass system may serve as an analytical tool for leg stiffness estimation in high-speed locomotion.

Introduction

The research on robot’s high-speed movement and operation ^{1 –4} is a significant topic, especially for legged robots, which could be used as an all-terrain platform in various missions. ⁵ In 1899, Muybridge used stop-motion photography to document animals’ locomotion, such as cats, dogs, camels, and horses. ⁶ Based on a large number of experimental biological data, researchers have summarized a spring-loaded inverted pendulum (SLIP) ⁷ model to describe the nature of animals’ dynamic motion, among which, spring represents the mechanical properties of animal tissues, such as bones, muscles, tendons, and ligaments. The SLIP model has become the most widely used model of legged robots in dynamic running. In this model, the centralized mass of the body and the stiffness of the legs have a decisive influence on energy storage and efficiency. ⁸ Based on the spring–mass model, Raibert and Tello had successfully developed a control scheme for legged robots from monopod to quadruped. ⁹ The iSprawl used tuned axial stiffness in conjunction with feed-forward position profiles of the prismatic leg. ¹⁰ The MIT Cheetah demonstrated a stable trot via the stiffness control of the stance legs. ¹¹ These studies confirm the importance of “leg stiffness” on locomotion behaviors. Legged robots’ high-speed locomotion is one of the challenging issues in robotic research. ^{11 –14} However, few studies covered the qualitative change of leg stiffness with the running speed increasing.

The use of mechanical impedance to control legged robots has been motivated by animals of different sizes and leg layouts using spring-like interactions with the ground. ^{15 –19} Many previous evidence have indicated that leg stiffness is independent of the forward speed, the stiffness of the leg (K _leg) changes slightly with speed in all observed animals (dogs, goats, horses, and kangaroos with a speed under 5 m/s), ^15,20,21 while other researchers reported a positive correlation between running speed and leg stiffness. ^{22 –26} The deviation comes from the different applied velocity, different leg sweeping angle, and the measurement accuracy. Animals are generally inclined to develop a lower leg stiffness if a larger sweep angle is available to support the velocity increment. ²⁴ Without significant changes in leg stiffness, a relatively stable unit-transportation-energy-consumption (COT) is achieved when speed increased gradually. ^21,27 However, the greatest touchdown angle by the leg spring during the ground contact phase was almost the same (around 28°) in animals of various body sizes, ¹⁵ indicating that all locomotion is restricted by the friction cone from increasing the sweep angle infinitely.

In 1990, a unique concept of “vertical leg stiffness” (K _vert) was raised by McMahon and Cheng to embrace the multivariables in running and it was quoted in many research works. ^15,24,25,28 The coupling mechanics show that the leg stiffness increases linearly and slowly with the horizontal velocity, while K _vert presents a significantly quadratic exponential growth when the other parameters remain unchanged. K _vert describes the vertical motion of the center of mass (COM) during the ground contact time and is essential in determining how long the spring–mass system remains on the ground. ²⁵ The vertical stiffness could be calculated by leg compression length and ground reaction force (GRF) (Figure 1(a)), where the maximum GRF can be measured by the force sensors installed on the feet/floor or by methods from TH Witte. ²⁹ K _vert is valuable in estimating the critical parameter, Stance Time. ³⁰ Since seeking the exact analytical solution of the motion during stance remains unsettled due to the nonintegral terms contained in the nonlinear dynamic equations, a common approach to tackle it is to simply the equations, ^{15,25,31 –34} resulting in an approximate solution of the state parameters. A mass–spring oscillator in the vertical direction was the most commonly used model to calculate the stance time as the gravitational force is ignored compared with the initial reaction force. In this article, we call this approximate model “harmonic motion model” (HMM) ¹⁵ . Under the approximation, the system will remain in contact with the ground for half of a resonant period of vertical vibration which is positively correlated with the reciprocal of K _vert, and K _vert is equal to K _leg when hopping vertically. ¹⁵

Figure 1.

Concept of the ‘Vertical stiffness’ and the ‘Triangle motion model’ in SLIP: (a) The ‘Vertical stiffness’ relies on the compression under the line between touchdown and liftoff and (b) the ‘triangle motion model’ assumes that the trajectory of the center of mass is a straight line. SLIP: spring-loaded inverted pendulum.

Nonetheless, the current definition of “vertical stiffness” still retains some ambiguities since it is difficult to distinguish the vertical displacement change which is caused by the tilting of the inverted pendulum or by the spring compression (view Figure 1(a)). In extreme cases, when an uncompressible inverted pendulum touches down with angle θ and is pushed up to the vertical position by inertia force, it will produce a negative △y, and △l equals to zero. Another case is that when the leg spring is hard enough, the compression of the leg in the vertical direction may be close to zero, which makes it difficult to calculate the actual value of △y. Therefore, it is hard to use the theory to estimate the stance time. In some instances of small leg sweep angle, researchers would use leg stiffness to replace vertical stiffness, and use the stance time of hopping in place to replace that with initial horizontal velocity and modify it step by step in the feedback process. ⁹ Since the Stance Time plays a key role in bridging stiffness selection and touchdown angle preparation on the premise of symmetric movement and neutral point theory, a more precise analytical equation of the stance time was adopted by taking gravity into account in McMahon and Chen’s paper. ^25,31,32

For the high-forward-velocity robots, more accurate analytical models were introduced. Geyer et al. derived analytical solutions to the stance map of the energy conservative SLIP model on symmetric stride. It was assumed that the sweep angle Δθ and the spring compression Δl during stance are sufficiently small (approach to zero). ³³ These remarkable assumptions guarantee the success of the approximating process. Similar assumptions were also used in Haitao Yu’s regular perturbation method, ³⁴ and both of them derived straightforward expansion solutions with high predicted precision. However, the approximations lose accuracy rapidly with an angular velocity at touchdown and radial spring compression rate increased, which are the classic characters of high-speed running.

With both the conciseness and accuracy in mind, and after extensive numerical simulations, a numerical fitting approach was taken to describe the relation between stance time and velocity. Then a modified stiffness with speed-related quadratic growth was applied to calculate the stance time. ³² Other empirical fitting formulas of stiffness from experiments were also accessible. ²⁵ However, the diversity of the actual robot system makes it laborious to fit the stiffness curve without pre-verified data under specific boundary conditions.

To sum up, accurate knowledge of leg stiffness is important as it is the basis for stabilizing control. It may lead to rapid failure in ultra-high-speed dynamic running if sustained errors are accumulated. ³¹ Although these studies have demonstrated the effectiveness of stiffness increment in high-speed running, yet almost none of them reported a rapid-response analytical approach for stiffness to handle the ultra-short contact phase. ³⁵ A better estimation of stiffness is required to avoid the wide range errors of the touchdown point as well as inadequate utilization of leg sweep angle. This work aims to pursue an accurate approximation of stance map validity in high-speed locomotion. The study illuminates the relationship between leg stiffness and speed, and the approximation model of the planar spring–mass system may serve as an analytical tool for estimating leg stiffness in high-speed locomotion, which is beneficial for making full use of the leg sweep angle and maintaining stability.

The remainder of this article is structured as follows: In the second section, the concise modeling of the SLIP while in high-speed locomotion is formulated in the “Problem formulation,” meanwhile a kind of “triangle motion model” (TMM) containing the simplest high-speed features is proposed. In the third section, the simplified solution of the stance phase is obtained by TMM and the stiffness strategy is introduced. In the fourth section, the method accuracy is analyzed by comparing the results between different analytical approaches and the results of numerical simulations. Multi-body dynamic simulations in ADAMS platform are outlined to consider the disturbance in the actual physical system and verify the theory. In the fifth section, limited conditions associated with incremental stiffness and the suitable range of the model are discussed. Finally, the sixth section presents the conclusion.

Problem formulation

As depicted in Figure 1(b), the conservative system consists of a point mass m and a massless prismatic leg connected to the hip with Hooke’s law. The spring–mass system lands on the ground with an initial touchdown angle θ ₀, forward velocity u, and downward velocity v_{y 0}. The velocity u takes the forward direction of bounce as positive and v_{y 0} is positive if it takes a vertical upward direction. The alternative touchdown and liftoff events contribute to the switching between the stance and flight phase. It is a passive model with no torque applied at the hip or active thrust in the elastic leg during the contact. The foothold keeps stationary without slipping. The system reaches the bottom as the vertical velocity decreases to zero. If the foothold falls on the “neutral point,” the trajectory of the COM will be symmetrical along the vertical line with a uniform touchdown and liftoff angle, and the forward velocity between the moment of touchdown and liftoff will be equal, while the symbol of the vertical velocity is reversed.

The incident point is taken as the origin of the horizontal motion. On the left side of the vertical line, the angle between the spring leg and the vertical line is negative, and on the right side of the vertical line, the angle is positive. As in high-speed locomotion, the touchdown angle is always swept to its maximum, which is represented as θ ₀. If the friction coefficient between the toe and the ground is represented by f, then the angle θ ₀ will be defined as a scalar that ${\theta }_0\in \left[0°,\text{ arctan }\left(f\right)\right)$ .

As mentioned above, leg stiffness acts as a critical bridge in speed regulation which determines the stance time and influences the leg sweep angle. The following equation describes the motion of the SLIP model in the vertical direction during the stance phase. It can be organized into a differential form

$m{a}_y=K_{\text{leg}}\text{Δ}l\text{cos}\hspace{0.17em}\theta -mg$

1

$\ddot{y}+{\omega }^{2}y={\omega }^2{l}_0\text{cos}\hspace{0.17em}\theta -g$

2

In this equation, ω denotes the undamped natural frequency. ${\omega }^2=K_{\text{leg}}/m$ . △l is the compression of the leg. Due to the nonlinear term ${\omega }^2{l}_0\text{cos}\hspace{0.17em}\theta$ in the right side of equation (2), the analytic solutions to the differential equations of the stance phase do not exist, of course not for the exact analytical solution of stiffness. Additionally, the wide range leg sweep angle which concerns cos θ is a major characteristic of high-speed locomotion to maximize stride length. It is important to note that, most research ignored this vital factor when simplifying the differential equations. ^{25,31 –34}

For the fast response control in ultra-high-speed locomotion, we have to obtain an approximated analytical solution of equation (2) by eliminating the nonlinear term while introducing the leg sweep angle θ and forward velocity u.

To simplify the equation, we reuse three similar assumptions as Marc Raibert used in the “neutral point” theory ⁹ which is accompanied by the feedback loop for fast error correction:

The forward velocity is nearly uniform in the contact phase, and the horizontal velocity u when the toe touches the ground is taken as the average velocity in the stance process.

The COM trajectory is symmetrical to the normal direction in steady running. Since it is unlikely to change stiffness frequently within one step, it is a good choice to maintain symmetry. The stance time error caused by the asymmetric motion is regarded as a disturbing factor which is corrected from step to step in the feedback process.

The COM trajectory is approximately a straight line connecting the incident point and the departure point.

Different from hopping in place, we call this “Triangle Motion Model” since the area in which the leg sweeps over forms a triangle. The kinematics and inertial navigation module can help obtain the horizontal velocity and apex height required in the application of these assumptions.

The TMM is a precise description of the SLIP model in high-speed motion. With these assumptions, we can revisit the dynamic equation of the contact phase to obtain a more accurate stiffness.

Method

Approximating the stance phase

As in the dynamics equation of the stance phase, we can replace cos (θ) in the nonhomogeneous terms as a square term Taylor expansion to realize linearization. The truncate error is no more than 0.0036 with a maximum touchdown angle θ ₀ = 30°, which covers almost the whole conventional friction cone, ¹⁵ thus equation (2) becomes

$\ddot{y}+{\omega }^{2}y={\omega }^2{l}_0\left(1-{\theta }^2/2\right)-g$

3

Combining the approximations of (a) and (b), we can obtain the stance time (T_{s 0}) required for the intended landing when the leg stretches to its maximum touchdown angle in high-speed running

$T_{s\text{0}}\approx 2l_0\text{sin}\hspace{0.17em}{\theta }_0/u(u>0)$

4

The constraint equation of real-time sweep angle is determined by the assumptions of (a) and (c) when timing begins at the moment of touching down

$\theta =\text{arctan}\left(\frac{ut-l_0\text{sin}\hspace{0.17em}{\theta }_0}{l_0\text{cos}\hspace{0.17em}{\theta }_0}\right)(0\le t\le T_{s0})$

5

We can expand arctan (x) to obtain a linearized expression of leg sweep angle as x because the maximum truncate error between the tangent angle and the real-time sweep angle is not more than 0.0538 within the contact phase

$\theta \approx \frac{ut}{l_0\text{cos}{\theta }_0}-\text{tan}{\theta }_0(\theta \in [-{\theta }_0,{\theta }_0])$

6

Thus, a time-varying nonhomogeneous term involved with velocity and touchdown angle elements is

$\left\{\begin{array}{c}\ddot{y}+{\omega }^{2}y=-{\omega }^2{l}_0\left(0.5A^2{t}^2+ABt+0.5B^2-1\right)-g\\ A=u/\left(l_0\text{cos}{\theta }_0\right) \\ B=-\text{tan}{\theta }_0 \end{array}\right.$

7

The general solution of equation (7) consists of a linear combination of the solution from the homogeneous differential equation and a special solution from the inhomogeneous equation generated by the undetermined coefficient method as equation (8)

$y(t)=C_1\text{cos}(\omega t)+C_2\text{sin}(\omega t)\text{+}a{t}^2+bt\text{+c}$

8

Also, the first-order derivative of equation (8) is

$\dot{y}\left(t\right)=-C_1\omega \text{sin}\left(\omega t\right)+C_2\omega \text{cos}\left(\omega t\right)+2at+b$

9

By combining the three equations above and substituting the initial boundary conditions below (equation (10)), we obtain the coefficients of the general solution in equation (11)

$\left\{\begin{array}{c}y\left(t=0\right)=l_{0}cos{\theta }_0\\ \dot{y}\left(t=0\right)=v_{y0}\end{array}\right.$

10

$\left\{\begin{array}{l}{C}_1=(0.5B^2+\text{cos}{\theta }_0-1)l_0+(g-A^2{l}_0)/{\omega }^2\\ C_2=(v_{y0}\text{+}AB{l}_0)/\omega \\ a=-0.5A^2{l}_0\\ b=-AB{l}_0\\ \text{c}=l_0-0.5B^2{l}_0+(A^2{l}_0-g)/{\omega }^2\end{array}\right.$

11

The symmetry assumption (b) helps to reduce the dimension of ω that the spring is deflected to its maximum and the vertical velocity of the COM is zero at mid-stance

$\dot{y}\left(t=T_{s0}/2\right)=0$

12

Substitute equation (12) into equation (9), we get the symmetric function

$C_1\omega \text{sin}(\omega T_{\text{s0}}/2)-C_2\omega \text{cos}(\omega T_{\text{s0}}/2)=a{T}_{\text{s0}}+b$

13

To solve ω, we introduced a new angular variable ϕ and scalar r. They satisfy the following definitions

$\left\{\begin{array}{l}r=\omega \sqrt{C_1^2+C_2^2}\\ \text{sin}\phi =C_2\omega /r\\ \text{cos}\phi =C_1\omega /r\end{array}\right.$

14

equation (13) can be rewritten as followed

$r[\text{sin}(\omega T_{s0}/2)\text{cos}\phi -\text{cos}(\omega T_{s0}/2)\text{sin}\phi ]=a{T}_{s0}+b$

15

thus

$\text{sin}(\omega T_{s0}/2-\phi )=(a{T}_s+b)/r$

16

With the equation (3), we were surprised to find that

$a{T}_{s0}+b\text{=}\frac{u(2l_0\text{sin}{\theta }_0-u{T}_s)}{2l_0{\text{cos}}^2{\theta }_0}\text{=0}$

17

which means that

$\left\{\begin{array}{l}\phi \hspace{0.17em}\text{=}\hspace{0.17em}\omega T_{s0}/2\\ \text{tan}\phi =C_2/C_1\end{array}\right.$

18

Unfortunately, C ₁, C ₂ are polynomial coefficients containing ω. We need to use Taylor expansion to solve the function since equation (18) is an implicit form. As we all know from traditional stance time estimation that T _s0 ≈ T/2, thus ωT _s0/2 ≈ π/2, which means that a duplication formula has to be used to achieve acceptable accuracy of Taylor expansion. As a trade-off, we used the triple-angle-formula

$\text{tan}\left(\frac{\omega T_{s0}}{2}\right)=\frac{3\text{tan}\left(\omega T_{s0}/\text{6}\right)-{\text{tan}}^3\left(\omega T_{s0}/\text{6}\right)}{1-3{\text{tan}}^2\left(\omega T_{s0}/\text{6}\right)}\approx \frac{108T_{s0}\omega -{\omega }^3{T}_{s0}^3}{216-18T_{s0}^2{\omega }^2}$

19

For the convenience of reading, we stipulate that

$\left\{\begin{array}{l}P=(0.5B^2+\text{cos}{\theta }_0-1)l_0\\ Q=g-A^2{l}_0\\ R=v_{y0}+AB{l}_0\end{array}\right.$

20

thus

$\frac{108T_{s0}\omega -{\omega }^3{T}_{s0}^3}{216-18T_{s0}{\hspace{0.17em}}^2{\omega }^2}=\frac{C_2}{C_1}\frac{R\omega }{P{\omega }^2+Q}$

21

Since ω ² = K _leg/m, the function can be simplified as

$T_{s0}^{3}P{({\omega }^2)}^2+(Q{T}_{s0}^3-18R{T}_{s0}^2-108P{T}_{s0})({\omega }^2)+(216R-108T_{s0}Q)=0({\theta }_0>0)$

22

With the computerized assistance we obtained the symbols of the coefficients in this quadratic equation and the solution

$\left\{\begin{array}{l}\Delta ={(Q{T}_{s0}^3-18R{T}_{s0}^3-108P{T}_{s0})}^2-\text{4}{T}_{s0}^{3}P(216R-108T_{s0}Q)\\ K_{\text{leg}}\text{=}m{\omega }^2=[\sqrt{\text{Δ}}-(Q{T}_{s0}^3-18R{T}_{s0}^2-108P{T}_{s0})]/(2T_{s0}^{3}P)\end{array}\right.$

23

So far, equation (23) constitutes the approximate analytical description of the leg stiffness of a high-speed system in passive motion. An accurate and stable stiffness could be obtained by one-step calculation compared with the multistep integrations and iterations in the numerical method.

It is important to note that, with the symmetric character, although the formula is derived on the scenario that the leg swings to its maximum touch angle θ ₀ in high-speed running, the approximation can also be applied to low-speed running while θ _td < θ ₀ based on the following definition

$\left\{\begin{array}{l}{T}_{s0}=T_{\text{s\_specify}}\\ {\theta }_0=\text{arc sin}(u{T}_{\text{s\_specif}y}/(2l_0))\end{array}\right.$

24

where T _{s_specify} is the stance time specified by the programmer.

Stiffness strategy

High-speed running animals have their own stiffness strategy. Step-length adjustment is the first choice of speed regulation in the low-speed movement. However, due to the limitation of structure and friction cone, further speed improvement can only be achieved by increasing leg stiffness and reducing the stance time. ³⁰ Based on the biological characteristics, a pre-calculated touchdown angle (θ _{0_pre}) deduced by the desired speed and the target stance time will be treated as the basis for strategy switching. As depicted in Figure 2.

Figure 2.

Process of stiffness selection in analytical solution.

1. When θ _{0_pre} < θ ₀, the touchdown angle increases with speed, while the stance time remains constant. θ _td = θ _{0_pre}.

2. When θ _{0_pre} ≥ θ ₀, the stance time decreases with speed, while the touchdown angle retains as θ ₀·θ _td = θ ₀.

Numerical simulation

We apply a numerical method in Matlab to evaluate the performance of the analytical solutions because the numerical method can obtain a sufficiently accurate solution satisfying locomotion details. The iterative process is shown in Figure 3.

Figure 3.

The monotonicity character and Newton method: (a) The monotonicity relationship between leg stiffness and stance time on the premise of a single variable (data from the simulation results in the next section) and (b) convergence function constructed by the Newton iteration method.

With the monotonicity ^{24 –26} of the function (see Figure 3(a)), we obtain a stance time with reasonable initial stiffness and state-values by Runge–Kutta method in the inner layer of the program and then construct an iteration formula in the outer layer to search the appropriate leg stiffness by the Newton-Iteration-method. The simulations are performed using the Runge–Kutta integrator, ode45 with event detection of leg touchdown and liftoff.

The analytical solutions of local linearization are used as the initial value of the search to accelerate the convergence speed since the Newton-iteration-method is sensitive to the initial values. The initial parameters of the SLIP model are set as in Table 1, which is close to the parameters of an adult cheetah. ³⁵

Table 1.

Model parameters for iteration.

Nomenclature	Symbol	Value
Total mass	m	30 kg
Initial height of the COM	h 0	0.8 m
Length of the leg	l 0	0.65 m
Original stance time	T s0	0.2 s
Maximum touchdown angle	θ 0	25°

COM: center of mass.

In regard to the parameters, we systematically change the input parameters relevant to the locomotion. It is important to note that, similar to animals, the maximum touchdown angle at low speed was restricted to 15° and increased gradually to 30° below 3 m/s and keeps this final maximum touchdown angle even in higher speed. Thus, preventing the calculated leg from being too soft.

Co-simulation platform

The Adams-Matlab co-simulation platform was also used to simulate the real physical environment of single leg hopping as well as to test the control performances. Comparisons are taken between the effects of TMM and HMM method on the accuracy of stiffness estimation and stable performance.

As seen in Figure 4, the robot consists of three parts: the body, the thigh, and the shank. A rotational joint connects the thigh to the body and a translational joint connects the shank to the thigh. It is constrained to move along a straight planar path. A body mass of 30 kg has been chosen. This is approximately half of the mass of an adult cheetah, allowing us to compare to this cursorial runner. The mass ratio between the body and the leg is set to 20:1 to avoid a significant influence on body posture when leg swings and the mass of the leg is evenly distributed. The static and dynamic friction coefficients were set to 0.8/0.75 because the maximum touchdown angles of most cursorial animals are around 28. ¹⁵

Figure 4.

The screenshot of co-simulation platform: Capture of the posture changes in one hopping cycle.

We adopted the well-regarded control strategy ‘Three-part Controller’ to manipulate rebound height, running speed and body attitude independently. ³⁶ The forward speed was controlled by the ‘neutral point’ theory with speed sampling at the beginning of the flight phase and the actual stance time of the final hopping cycle was selected during the calculation of the touchdown angle.

Within the boundaries of the stiffness estimation model, the desired touchdown angle and the desired stance time also include the PID value of velocity deviation during the calculation process because feedback was applied in the actual control framework and the estimation model should use data that is consistent with reality.

Since actual locomotion systems must have energy dissipation and injection, ³⁷ a damping ratio of 10 N·m/s was added to the translational joint to prevent leg tremor, also a rotational damping ratio of 10⁻⁶N·m·s/deg was added to the hip. A thrust was applied in the radius direction when the leg rebounds from the ground and the pitch angle was adjusted during the stance phase with a small hip torque.

The total hopping parameters used in the simulation are presented in Table 2.

Table 2.

Parameters in high-speed hopping.

Parameters	Value
Body mass	30 kg
Mass ratio between body and leg	20:1
Static friction coefficients	0.8
Dynamic friction coefficients	0.75
Maximum touchdown angle	25°
Translational damping ratio	10 N·m/s
Rotational damping ratio	5 × 10−7 N·m·s/deg
Acceleration	0.2 m/s2 /0.1 m/s2
Initial stance time	0.2 s
Minimum stance time	0.07 s

To avoid the influence of acceleration on observation, we set up the acceleration process with a multistep approach with a time interval for each speed. The predetermined initial stance time is 0.2 s, and the minimum stance time in super high-speed is set to 70 ms because the leg inertia is insurmountable.

Result

Numerical simulation result

Comparisons between the analytical and numerical methods for the required leg stiffness are illustrated in Figure 5. As in Figure 5(a), the TMM method (red solid line) fits well with the numerical result (blue solid line) because the touchdown angle and velocity have been taken into account within the model. The maximum stiffness error is less than 12% even in ultra-high-speed, while the method of HMM draws a significant gap as the speed increased to 10 m/s, almost twice as much as other methods. The vertical stiffness ‘K _vert’ (Pink dotted line) which is derived from the leg stiffness of numerical results is highly consistent with Thomas’s findings ¹⁵ when contrasting with K _leg.

Figure 5.

Stiffness increase with vertical speed: (a) Stiffness estimation results influenced by forwarding speed; (b) local enlargement of Figure 1(a); (c) stiffness estimation results influenced by the touchdown angle; and (d) stiffness estimation results influenced by the vertical velocity of the COM. COM: center of mass.

As seen in Figure 5(b), all these methods obtain a very similar stiffness prediction when forward velocity is less than 3 m/s. The HMM method is smaller than the numerical result at first and then grows significantly at the high-speed component, while the TMM method remains a little larger than the numerical result. Due to the limitation of Taylor’s expansion accuracy (Higher accuracy requires solving the complex quaternary equation), the TMM method could not fit with the numerical result completely. On the other hand, the symmetry of the numerical simulation is influenced by the assumption of constant forward velocity, which may cause the deviation of the touchdown angle from the theoretical value. It can also be seen in Figure 5(b) that leg stiffness increases dramatically with maximum touchdown angle θ ₀ becomes steeper since a smaller step length requires a much shorter stance time.

In Figure 5(c), K _leg is plotted vs the vertical velocity v_{y 0}, and it is apparent that K _leg increases approximately linearly with v_{y 0} as the forward velocity is set up as a constant value of 10 m/s.

By using different model results as the initial values of the numerical searching, the iteration time is presented in Figure 6(a). The TMM result is sufficiently close to the real value that most iterations can be completed within three steps, while HMM fits best at the intersection point but performs poorly in other points. As depicted in Figure 6(b), the stance time errors by substituting K _leg from TMM and HMM results also reveal the accuracy of the estimation, and the relative error trend in TMM is much smaller than HMM.

Figure 6.

The accuracy of the new analytical solution: (a) Iteration times by different initial values; the TMM method is sufficiently close to the real value that most iterations can be completed in three steps; (b) Error of stance time by different calculation method, the relative error in TMM is smaller than HMM. TMM: triangle motion model; HMM: harmonic motion model.

Adams-Matlab co-simulation

Despite the simplicity of the control law mentioned above, the simulation results show the stability of the forward speed, maximum height, and posture in both applications. However, detailed graphs in Figure 7 show a better velocity tracking with less fluctuation and higher stable speed by using the TMM method. Since the fluctuation in the stance time has a strong correlation with planned touchdown angle, the feedback system has to work consistently to compensate for the velocity changes during stance phase if the HMM method is adopted. The fluctuation is small in relation to velocity itself in high-speed locomotion, which is consistent with the assumptions of uniform horizontal velocity at stance phase.

Figure 7.

Velocity tracking: (a) velocity tracking by TMM method. The TMM method achieves a higher and more stable speed without oscillation; (b) velocity tracking by HMM method. TMM: triangle motion model; HMM: harmonic motion model.

Figure 8 (a) and (b) shows the leg sweep angle and its symmetricity during stance by TMM method. The magnitude of asymmetry during stance is relatively small when compared with the touchdown angle, which clarifies the validity of assumption (b). Due to the feedback within the physical system, the desired stance time in Figure 8(c) and (d) varied slightly with different methods. But the stance time in the TMM method is fundamentally consistent with the desired stance time, particularly during the high-speed interval (3–8 m/s). Yet the HMM method has considerable deviations at low-speed and high-speed sections and the minimum deviation occurs around 3 m/s, which is the intersection point of the two stiffness curves in Figure 6(b). What can be found particularly is the highly consistency of the stance time errors between Figure 6(b) and Figure 8. These phenomena were caused by the fact that HMM adopts lower stiffness at low speed, while much higher stiffness at high speed. Compared with the ultra-short stance time in ultra-high-speed locomotion, any small error will have an impressive effect on stability.

Figure 8.

Hip angle performance and stance time comparison: (a) Leg sweep angle in stance time by TMM method; (b) the representative $P_{\text{a}sym}={\theta }_{td}+{\theta }_{lo}$ to show the magnitude of symmetricity; (c) stance time obtained by TMM and HMM method; and (d) stance time obtained by the HMM method. TMM: triangle motion model; HMM: harmonic motion model.

The stiffness applied in the simulation and the corresponding effects on leg swing angles during the stance phase were presented in Figure 9. Due to the feedback loop, the control system can adjust the touchdown angle in real time to ensure the symmetry as long as possible. Although some disturbance exists in the actual physical system of TMM method such as the ground impact and motion delay caused by the inertia (Figure 9 (b)), the sweep angle in the TMM method is fully utilized and has a more stable performance than HMM. The HMM method generated a much narrow and unstable leg sweep angle in high speed owing to the great increment in leg stiffness.

Figure 9.

The stiffness adopted and the corresponding effect: (a) The HMM method adopts a much higher stiffness in high-speed section and (b) the sweep angles during stance by TMM and HMM show that the sweep angle in the TMM method is fully utilized. TMM: triangle motion model; HMM: harmonic motion model.

Discussion

A vital question regarding the stiffness selection in high-speed locomotion is how to limit maximum stiffness? The solution of leg stiffness in the TMM method mimics a similar result as speed increased to more than 10 m/s, even though it is not shown here because ultra-high-speed locomotion with a single leg is still quite difficult for a real robot limited by leg inertia, structural strength, and actuators.

Based on the consistent ground impulse received in a steady-state, the average GRF acting on the toe will increase significantly with the shortening of stance time. Although some methods can be taken to alleviate the ground impact, such as leg retraction before landing, the increase of impact force is inevitable and excessive impact may cause the fracture of mechanical structures. High stiffness also requires the actuators to generate a high force at a high frequency in active stiffness since the damper diminish significantly more mechanical energy as the impact force is increased. A mass–spring–damper (MSD) model developed by Alexander et al. ³⁸ also helps to highlight the relationship between leg stiffness and paw-pad stiffness in that it should not exceed five times of paw-pad stiffness in order to prevent chatting on the toe. These limitations make it impossible for us to increase leg stiffness infinitely.

Alternatively, if the compliance is too soft, it cannot provide a sufficient buffering force within the allowable compression range to balance the weight and resist the landing shock by squeezing energy into the spring. A softer leg may also degrade the positioning accuracy of the foot, and bigger amplitude of torso vibration caused by the compliance would add more noise and decrease sensing accuracy of the sensors on the back.

The TMM method for calculating leg stiffness is helpful to get real-time stiffness quickly through these complex boundaries in high-speed running. The accurate solution also helps to reduce the risk of using a much stiffer leg or much softer leg to harm the structure.

The Multi-leg cooperation with different gaits or spine movement could play a much more significant role in achieving a higher speed mission after the Limited Stance Time (LST) appeared.

Conclusions

Leg stiffness plays an essential role in touchdown angle and speed regulation. In relation to the HMM model, we propose a TMM model to solve the approximate analytical solution of stiffness in high-speed running. The TMM method accurately grasps the characteristics of high-speed running in that large range of leg sweep angle and small fluctuation of the COM occur. Both the numerical simulation and the multi-body dynamics simulation in ADAMS show the accuracy of the TMM method in the accurate prediction of leg stiffness even in ultra-high-speed. It is also adaptable to low speed and small sweep angle cases. The accurate solution of this model also helps the structure of the robots to avoid damage caused by a much stiffer leg or a much softer leg. Its accurate solution has a great advantage on velocity tracking and leg sweep angle utilizing as well as the system stability.

Nomenclatures

F	Friction	N	S	Step length	M
F	Ground reaction force	N	Ω	Natural frequency	rad/s
M	Body mass	kg	A	Acceleration	m·s− 2
Θ	Touchdown angle	rad	G	Gravitational acceleration	m·s− 2
L	Leg length	m	Subscript
U	Horizontal speed	m·s−1	max	The maximum value
vy	Vertical speed	m·s− 1	td	Leg touchdown
K	Stiffness	N·m− 1	leg	Leg
T	Time	s	S	Stance phase
H	Center of mass height	m	vert	The vertical direction
Δl	Compression length of the leg	m	0	Initial or desire situation
Δy	Distance from center of mass in stance	m	y	Vertical direction