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

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. 5 In 1899, Muybridge used stopmotion photography to document animals' locomotion, such as cats, dogs, camels, and horses. 6 Based on a large number of experimental biological data, researchers have summarized a spring-loaded inverted pendulum (SLIP) 7 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. 8 Based on the spring-mass model, Raibert and Tello had successfully developed a control scheme for legged robots from monopod to quadruped. 9 The iSprawl used tuned axial stiffness in conjunction with feed-forward position profiles of the prismatic leg. 10 The MIT Cheetah demonstrated a stable trot via the stiffness control of the stance legs. 11 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][12][13][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][16][17][18][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][23][24][25][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. 24 Without significant changes in leg stiffness, a relatively stable unittransportation-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, 15 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. 25 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. 29 K vert is valuable in estimating the critical parameter, Stance Time. 30 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][32][33][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) 15 . 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. 15 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 q and is pushed up to the vertical position by inertia force, it will produce a negative 4y, and 4l 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 4y. 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. 9 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 Dq and the spring compression Dl during stance are sufficiently small (approach to zero). 33 These remarkable assumptions guarantee the success of the approximating process. Similar assumptions were also used in Haitao Yu's regular perturbation method, 34 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 speedrelated quadratic growth was applied to calculate the stance time. 32 Other empirical fitting formulas of stiffness from experiments were also accessible. 25 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. 31 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. 35 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 springmass 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 q 0 , forward velocity u, and downward velocity v y0 . The velocity u takes the forward direction of bounce as positive and v y0 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 q 0 . If the friction coefficient between the toe and the ground is represented by f, then the angle q 0 will be defined as a scalar that q 0 2 0 ; arctan f ð Þ ½ Þ . 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 In this equation, o denotes the undamped natural frequency. ! 2 ¼ K leg =m. 4l is the compression of the leg. Due to the nonlinear term ! 2 l 0 cos q 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 q 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][32][33][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 q and forward velocity u.
To simplify the equation, we reuse three similar assumptions as Marc Raibert used in the "neutral point" theory 9 which is accompanied by the feedback loop for fast error correction: (a) 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. (b) 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. (c) 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 (q) 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 q 0 ¼ 30 , which covers almost the whole conventional friction cone, 15 thus equation (2) becomes Combining the approximations of (a) and (b), we can obtain the stance time (T s0 ) required for the intended landing when the leg stretches to its maximum touchdown angle in high-speed running 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 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 Thus, a time-varying nonhomogeneous term involved with velocity and touchdown angle elements is 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 cosð!tÞ þ C 2 sinð!tÞþat 2 þ btþc ð8Þ Also, the first-order derivative of equation (8) is 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) The symmetry assumption (b) helps to reduce the dimension of o that the spring is deflected to its maximum and the vertical velocity of the COM is zero at mid-stance Substitute equation (12) into equation (9), we get the symmetric function To solve o, we introduced a new angular variable and scalar r. They satisfy the following definitions Equation (13) can be rewritten as followed With the equation (3), we were surprised to find that which means that Unfortunately, C 1 , C 2 are polynomial coefficients containing o. 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 oT s0 /2 % p/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-angleformula For the convenience of reading, we stipulate that thus Since o 2 ¼ K leg /m, the function can be simplified as With the computerized assistance we obtained the symbols of the coefficients in this quadratic equation and the solution 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 q 0 in high-speed running, the approximation can also be applied to low-speed running while q td < q 0 based on the following definition 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. 30 Based on the biological characteristics, a pre-calculated touchdown angle (q 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. 1. When q 0_pre < q 0 , the touchdown angle increases with speed, while the stance time remains constant. q td ¼ q 0_pre . 2. When q 0_pre ! q 0 , the stance time decreases with speed, while the touchdown angle retains as q 0 Áq td ¼ q 0 .

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.
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. 35 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. 15  We adopted the well-regarded control strategy 'Three-part Controller' to manipulate rebound height, running speed and body attitude independently. 36 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, 37 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 À6 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.
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.

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 15 when contrasting with K leg .
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 q 0 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 y0 , and it is apparent that K leg increases approximately linearly with v y0 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.

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

Discussion
A vital question regarding the stiffness selection in highspeed 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 massspring-damper (MSD) model developed by Alexander et al. 38 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 highspeed 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-highspeed. 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 Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by Natural Science Foundation of China