Design and verification of a parallel elastic robotic leg

This paper presents the design and experimental verification of a parallel elastic robotic leg mechanism that aims to capture the dynamics of the linear mass-spring-damper model. The mechanism utilizes a wrapping cam mechanism to linearize the non-linear force resulting from the elongation of the parallel elastic element. Firstly, we explain the desired dynamics of the mass-spring-damper model, including the impact transitions, and the design of the wrapping cam mechanism. We then introduce a system identification procedure to estimate the parameters of the leg mechanism corresponding to the dynamic model. The estimated parameters are tested with a cross-validation approach to evaluate the mechanism’s performance in tracking the desired model. The experimental results show that the passive dynamics of the mechanism resemble the linear model as intended. Thus, the robot provides a basis for using parallel elastic actuation while using model-based controllers that benefit the analytic solutions of the linear model.


Introduction
As evident from nature, compliant actuation and control principles form a basis for highly dynamic and agile locomotion.Recent legged platforms have utilized either high-torque driving units that provide lossless transmission enabling the implementation of virtual impedance control [1][2][3], or passive elastic elements that provide greater energy efficiency at the cost of introduced control and engineering difficulties [4,5].The spring-loaded inverted pendulum model (SLIP) and its extensions have long been accepted as a basis model explaining animals' efficient and stable running dynamics [6][7][8].The model ′ s simplicity, along with its energetic and stability characteristics, make it a control objective beyond an analysis tool for legged locomotion [9][10][11][12].
Controlling a legged robot to adopt the dynamics of a simpler model like SLIP is desirable; however, it is not straightforward, as physical realizations of robots have complex mechanical structures.The robot's mechanics should be designed to be suitable for the target model to some extent, and further advanced control actions may be required.In this context, our study aims to bridge the gap between theoretical models and practical implementations, particularly for SLIP-based robotic systems.By focusing on the design of a parallel elastically actuated robotic leg, we address the challenge of translating the simplistic yet effective dynamics of the SLIP model into a physically realizable and efficient robotic structure.This endeavor is not only about achieving a technical feat, but also about paving the way for more energy-efficient, agile, and dynamically stable robots, drawing inspiration from the natural efficiency observed in animal locomotion.
Raibert [13] adopted a simplified approach where the non-linear leg dynamics are neglected, and the desired spring forces are exerted virtually on the body of a quadruped.Later, De and Koditschek [14] adopted a similar approach on the Minitaur robot platform to realize decoupled springy hoppers as a basis for the running control of quadrupeds.Pratt et al [15] followed a similar approach on the 'spring turkey' biped to achieve virtual compliance on the robot.Hutter et al [16] considered full-body dynamics of an articulated leg and used operational space control methods to compensate for the non-linear effects.Direct implementation of passive elastic elements to obtain desired virtual compliance requires careful design of the geometric placement and dimensioning of the elastic components.Early examples [13,17,18] mostly use prismatic actuation accompanied by loaded springs or have compliance mostly in limited degrees of freedom restraining the versatility of motion.Abate et al [19] provided the compliance and inertia analysis of a segmented leg as a design framework for deciding the locations of the elastic elements and the mass distribution of the leg to capture the mass-spring dynamics.The ATRIAS robot [4] is a notable example employing a lightweight parallelogram leg to constrain the ground reaction forces to be on the radial axis as in the SLIP model.
Drawing inspiration from nature, where we observe parallel and series elasticity in the musculoskeletal system of animals [20], an increasing number of robots have been developed with parallel elastic actuation (PEA) in their designs.In principle, PEA is favorable over series elastic actuation (SEA) since the parallel elastic element provides additive forces to the actuator, which enables the usage of relatively small actuation units by reducing the required peak torques of the actuator and energy consumption.Yesilevskiy et al [21] compared the efficiency of monopedal hopping with series and parallel elastically actuated configurations, and reported that PEA performs better when the actuation unit has a small transmission ratio.Another advantage is that the PEA provides a higher control bandwidth compared to SEA, since the direct regulation of the forces at the end-effector is possible [22][23][24].The Spacebok robot [25] uses parallel springs in its parallel leg mechanism.It benefits from the high output forces from the parallel spring to reach high jumping heights, which makes it favorable in unstructured environments.Liu et al [26] and Nan et al [27] employed PEA in segmented robotic leg mechanisms with a focus on the compliance adjustment mechanisms to obtain efficient motion while eliminating the joint dexterity problem.
The majority of work on passive elastically actuated robots evaluates the impact mitigation, the energetic efficiency, and the torque capacity of the robots.Having integrated compliance in the joint level, due to the non-linear mapping between joint and task space velocities, most of the robots obtain a nonlinear compliant behavior in the task space [4,26,27].The purpose of this study is to exploit the efficiency of PEA while preserving SLIP dynamics to utilize simple model-based control.In our previous work [28], we introduced the synthesis of a wrapping CAM mechanism for the linearization of the virtual spring forces and validated the concept design with static force measurements.In this study, we extend these static measurements to dynamic hopping experiments on a monopedal robot with an added impact transition map to provide comprehensive proof of the design concept.
The following section details the target springmass model dynamics and the leg mechanism's design principles.Section 3 firstly introduces the experimental setup and the parameter estimation process of the mechanism.Then, the experimental results evaluating the passive dynamic structure of the mechanism are presented.Section 4 presents our conclusions to the experimental results and a discussion on possible future directions.

Description of the mechanism
The design goal of the leg mechanism is to achieve the linear stance dynamics of a vertical hopper depicted in figure 1(a).The equation of motion of the system in the stance phase is described by m is the point mass lumped at the body, k is the equivalent spring constant, and d s is the viscous damping constant.F in is the force created at the toe by the actuator and z 0 is the position where the spring has the rest length.
The considered leg mechanism is shown in figure 2. It is a four-link mechanism with passive compliance.The primary function of the mechanism is to realize a hopping motion along the vertical axis with the aid of the actuation input and the compliant element.The mechanism is constrained to move along the vertical axis with prismatic joints at the body and the toe.All links are interconnected with revolute joints.Only one actuator is placed on the body to apply torque input between the two upper links.The compliant element is a tension spring that stores energy while the leg shortens in the radial direction.As the upper links rotate relatively, the spring is wrapped on the surface of the cam.The system is parallel elastically actuated, considering the generated force at the end effector (toe).Because, the forces exerted on the toe by both the actuator torque from the electrical motor and the additional forces applied by the tension spring directly combine.As the central concern of the mechanism is to capture the linear dynamics in (1), the primary assumption for the mechanism is that the link masses are negligible, so we can disregard the inertial effects.With the neglected inertial effects, one can relate the actuation torque, τ ,  to the input force in (1) simply by using the virtual work equivalence.The relation between the vertical position and the generalized coordinate of the mechanism θ is obtained by where z c is the vertical position of the contact point and L is the constant link length.Then, using the relation between the infinitesimal change in position variables, the input force on the vertical axis becomes

Wrapping cam design
The key aspect we present to mimic SLIP compliance in the linkage-based leg is the design of the wrapping cam component.The cam is integrated into the system such that as the mechanism contracts the ground, the spring elongates and a string attached to the spring is wrapped around the cam to control the spring extension.The cam profile should be designed such that the spring extension results in an equivalent radial elastic force, which is a linear function of the radial displacement of the center of mass (COM).
The design methodology and the working mechanism are adopted from our previous work [28].
The cam is rigidly attached to upper link 2, as indicated in figure 2. To simplify the analysis, the principle of kinematic inversion is used, where the cam frame is fixed and all the motion is given to upper link 1.This approach, as well as the kinematic diagram, are summarized in figure 3. The synthesis is based on three equations.The first is a vector loop formed by the vector ⃗ L, which stands for the upper link 1 length and direction, vector ⃗ R, which is the tangent point on the cam where the string leaves, and vector ⃗ T, which goes along the cam spring and starts from the tangent point.The second is the relation between the required spring elongation and the amount of wrapping, and the third is the tangency relation, which connects the infinitesimal amount of wrapping and the tangent vector on the cam curve.
The amount of required wrapping can be derived from the equivalence we enforced between the SLIP and physical system dynamics.The input forces are equated through torque control, as discussed in the previous section.Kinetic energies and damping losses can be equated simply by correct parameter identification.The cam's aim is to equate elastic potential energies: where k c represents the physical cam spring, k represents the virtual SLIP spring, and w represents the total elastic element elongation.To solve for w(θ), a solution with a minus sign is preferred to achieve an increasing function because, as θ increases, the spring should elongate and the potential energy should increase as in the SLIP: Here, L 0 can be defined in terms of an initial angle θ 0 , such that L 0 = 2L cos(θ 0 ).Also, the ratio k/k c can be parameterized with a stiffness ratio K.This is an important design parameter as it determines which spring stiffness and cam size should be selected.The vector loop, the wrap and elongation relation, and the differential tangency condition can be written with complex notation as where the following vector definitions hold: in which 2θ represents the rotation of the link with respect to the fixed cam frame.To utilize the tangent condition, equation ( 6) is differentiated with respect to θ and vector definitions are substituted, while we know that L is constant as it is the link length but β changes as a function of θ: noting that e (−i π/2) = −i.Substituting equation (8) gives us and this can be simplified by differentiating and substituting equation (7) as dw = ds + dT, and dividing both sides by e iβ : equating the real and imaginary parts gives us The intermediate variable β can be solved by using equation ( 14): and dβ/dθ can be solved by differentiating (14) and knowing β: then the intermediate variable T can be solved by using equation ( 15), and utilizing (16), and ( 18): then, the cam profile can be solved as Note that in CAD design, ⃗ R should be rotated by −(π/2 − θ 0 ) or the link should be tilted with θ 0 in the drawings.
CAM surface synthesis has some limitations.These either arise mathematically during derivation steps or physically exist but are mathematically expressable.For example, arccos is only defined between −1 and 1, which puts a limit on the demanded wrapping function.These constraints, existence, and uniqueness issues are discussed in the literature.Numerical optimization approaches satisfying these constraints are widely used to relax the analytical limitations for cam surface generation [29,30].
One should observe the inverse relation between the cam size and the stiffness ratio K.The physical spring force decreases with the decreasing physical spring constant; however, for a given equivalent stiffness, the ratio K increases and the amount of wrapping required increases.Therefore, if one wants to use a small cam to decrease leg inertia, one should use a stiffer spring, which will increase the loads in the system components.This is discussed in detail in our previous work [28], with differences in the coordinate system definitions.

State transitions and impact
Although we have described the equation of motion for the stance phase, the complete hopping motion consists of cyclic changes between the stance and flight phases.Therefore, one should also describe the flight dynamics and the transition between the two states for the complete analysis.In the flight phase, the mechanism simply undergoes a free-fall motion.However, due to the prismatic joints in the body and toe, a more realistic model also includes a lumped viscous damping loss in the flight phase.The descriptive model for the flight is illustrated in figure 1(b).The corresponding equation of motion is We assume ideal plastic collision between the toe and the ground at the touchdown and lift-off instants.While an obvious collision between the toe and the ground occurs at touchdown, the other collision occurs at lift-off due to the mechanical limit of the mechanism.While plastic collisions cause an instant stop for the toe, the rest of the mechanism undergoes a sudden change in velocity.The impact is simply modeled as an abrupt velocity change of the COM of the system: The term β is a parametric expression of the abrupt velocity change at the phase transitions.The subscripts td and lo refer to the touchdown and liftoff instants, respectively.We assume that the flightto-stance transition occurs when the COM position becomes equal to the free length of the leg, z = 2Lcos(20 • ).A stance-to-flight transition occurs when the ground reaction force becomes zero, F GRF = 0, or the COM position becomes equal to the free length, as in the flight-to-stance transition.The ground reaction force is estimated as follows:   and the height encoder.Communication is done through the EtherCat protocol with a 1 kHz cycle frequency.The Simulink target PC is the main controller and acts as an EtherCat master.The data acquisition device and the motor driver act as EtherCat slaves to provide sensory feedback to the controller and to receive torque/position commands for the motor.

Hopping control
In our experiments, we update the control commands based on an event-triggered predefined cyclic openloop torque input scheme.As the parallel spring creates a dominant force along the motion trajectory, the input torque should be in sync with the linear velocity of the system to obtain a smooth hopping motion.This can also be interpreted as exciting a system with its natural frequency.Since the input torque can change the system ′ s momentum only in the stance phase, the input torque takes the form of where T is the amplitude of the torque in N m −1 , P is the actuation period in seconds, t is the time, and every cycle starts with the touchdown instant, t d .Figure 6 shows one example of system state evolution and applied input torque under this control policy.

Parameter estimation
Our system identification method in this study relies on fitting the physical system parameters to the adopted dynamical model for the systematically obtained experimental data.In this regard, the objective is to find the lumped mass of the system m, the linear spring coefficient k, the damping coefficients in the stance and flight phases d s and d f , and the impact parameters β td and β lo .The main idea is to find a set of these physical parameters that minimizes the difference between the simulated and measured output states.With the discretization of the dynamic formulation given in ( 1), (23), and ( 24), the state evaluation can be formulated as a function of the parameters as where α is a vector of parameters, and the state vector is defined as x = [z, ż] T .Note that the direct measurements of both states are available.Only apex-to-apex single-stride hopping trajectories are considered.The initial conditions of the simulations are selected directly from the experimental trajectories, and the other states are obtained with (27).We evaluate the accuracy of the simulated hopping trajectory with a cost function that considers the normalized mean square error (NMSE) in the COM position and velocity as follows: Here, z and ż are the vectors of the sampled COM positions and velocities, respectively.n is the index of the sample points and N is the total number of samples of the vector.The subscripts 'exp' and 'sim' refer to the sampled experimental and the simulated trajectories, respectively.The mean values of the vectors are defined as The cost function aims to minimize the tracking error both in the position and the velocity.As we use the squared error, the cost function is of the quadratic form, and the weights for the position and velocity should affect only the convergence rate.Therefore, we simply use the unit weight for each component.Parameters related to the different phases of locomotion are estimated separately using distinct initial conditions of the states from the given experimental trajectory.For each phase, multiple singlestride hopping experiments are conducted, and the overall objective function is set to find the parameters that result in the best fit for all of the experimental trajectories.Hence, the overall objective is to find a set of parameters that minimizes the sum of the costs of all trajectories: where the subscript i is the index of the experimental trajectory.We use the fmincon function of MATLAB [31] to solve (29) with physically consistent parameter constraints.The initial point of the optimization is selected randomly.Figure 7 shows an example single-stride experiment.The initial position of the simulation is selected near the transition point of the related phase.The preliminary results indicate that even if a sensible set of parameters is estimated, when we select a noisy initial condition, the simulated trajectories deflect significantly from the experiments.
To eliminate this, the noise in the velocity profile is filtered out with a 15 Hz cutoff frequency with a third-order zero-phase-lag Butterworth filter (using the filtfilt method of MATLAB).This method eliminates the phase shift resulting from the filtering by processing the data in both the forward and backward directions.The cutoff frequency of the filter is chosen by empirical iterations with a focus on not distorting the main character of the trajectory.As the nominal stride frequency of the leg is between 1.5 and 2.5 Hz, the cutoff frequency is 7-10 times the operational range.Therefore, we simply filter out the noise of the encoder while preserving the natural hopping dynamics of the mechanism.The noisy transition in velocity just after the impact instants is encircled and the filtered velocity profile is presented in figure 7(b).Although it is not feasible to directly obtain all physical parameters from the measurements, using the measured parameters wherever possible is beneficial to decrease the number of variables to be estimated.The total mass of the system and the spring constant are the two lumped quantities that can be measured with relatively high accuracy.As the natural frequency of the system is a distinctive dynamic property, we accept the measured quantity of either the lumped mass or the spring coefficient and estimate the other from the experiment 7 .The lumped mass of the system is simply measured to be 1.62 kg.The equivalent spring constant is measured using a force plate following the procedure in [28].The ground reaction force is measured at discrete deflection points with a AMTI he6x6 force plate, and a linear fit to the measurements is obtained as shown in figure 8.The equation for the best fit line is 605 N m −1 +2.7 N. The offset force is associated with the pretension force of the tension spring [28].

Experimental results
During the estimation process explained in section 3, different excitation trajectories are obtained by changing the amplitude and period of actuation in (26).Firstly, the stance and flight phase parameters are estimated.Then, fixing the mean values for the estimated parameters, the impact parameters are estimated to obtain a complete hopping trajectory.
For all the physical parameters, the estimation and validation sets are formed using the k-fold crossvalidation scheme.For the flight phase, 20 experiments are conducted, and fivefold cross-validation is used for validation.For the stance phase and the impact parameters 60 and 52, valid experiments are used with tenfold validation.The distributions of the NMSE, mean absolute error, and maximum absolute error for each validation run of the stance and flight phases are presented in table 1.The estimated parameters result in a NMSE of less than 2% in both the position and velocity of the stance phase.While the percentage error is also small in the flight velocity, one can notice that it is relatively large in the flight position.This is due to the fact that the flight phase consists of a smaller motion range, and the normalization we used in (28) focuses on the dispersion of the trajectory elements around their mean value, rather than explicitly considering the extent of the motion range.One can see that even the normalized error is relatively high, and the mean absolute error is smaller than 2 mm.Confirming the performance of the isolated stance and flight phases, we fix the mean values to estimate the impact parameters.Apex-to-apex single-stride trajectories are used to find the touchdown and lift-off coefficients.The mean values of the estimated parameters are presented in table 2. As seen in the table, the standard deviation of the damping coefficient in the flight phase is notably larger than that of the other estimated parameters.Our observation is that the linear guide used in the setup provides a frictionless motion, so that the inertial forces become dominant against the damping losses in the flight dynamics.Therefore, even with the relatively large variation in the estimated damping coefficient, the validation results do not differ significantly.Single-stride experiments that are used to estimate the impact coefficients are used to validate the tracking performance of the overall parameter set. Figure 9 shows the distribution of the NMSE for the complete trajectory.One can see that with the estimated set of parameters, we can predict the overall trajectory of both the position and velocity of the robot with an NMSE of less than 5%.Two example simulations showing the tracking performance of the overall estimated parameters are presented in figure 10.One can see the effect of the impact as an instantaneous velocity change at the touchdown and lift-off instances.Note that the velocity change in the experiments at the impact seems to be smooth as the trajectory is filtered.Figure 10(b) is elaborately chosen as one of the worst tracking cases in all 52 experiments.The resulting lift-off parameter larger than 1 seems to be contradictory, as the impact is commonly associated with the energy loss.This result is related to the change of the velocity profile between the flight and stance phases.In the stance phase, the kinetic energy of the leg has both rotational and translational velocity components.When the system jumps to the flight phase, the mechanism has pure translational velocity along the vertical axis; thus, we observe an increase in the velocity regardless of the net energy loss of the system.
After obtaining the physical parameters defining the stance and flight dynamics, the main contribution of the impact model is to carry the system to the correct initial conditions at the phase transition instants.To evaluate the effect of the impact, we compare the apex-to-apex prediction performance  of our model to the performance of a model that neglects the impact.In the impact-free model, we simply take the impact coefficients in (24) as unity, while we take the mean values of the estimates given in table 2 for the model with impact.Figures 11(a

Conclusion
This paper presents the detailed design and dynamic verification of a parallel elastically actuated robotic leg mechanism.The system's passive (non-actuated) dynamics are tailored to realize the mass-springdamper dynamics by utilizing a novel wrapping cam mechanism.In the paper, we first introduced the dynamical model of the system together with the hybrid transition dynamics.Then, we validated the approach using a prototype hopper system through systematic experimental analysis.In our empirical analysis, we mainly focused on estimating the physical parameters of the prototype system and the associated dynamic model, and the cross-validation of these parameters and the model itself by utilizing categorically different experiments, which effectively elicit the generalization ability and potential applicability on the various platforms.
More specifically, the validation results show that the piece-wise linear (though still hybrid since transitions depend on state variables) model closely tracks the experimental single-stride jumping trajectories.The results also indicate the importance of incorporating touchdown and lift-off dynamics since even simple (linear) transition models significantly increase the apex-to-apex prediction performance, which is critical for adopting return-mapbased control policies.Accuracy in predicting the robot's apex position and timing leads to the usage of simple model-based controllers.This study forms the basis for the utilization of highly efficient PEA in robotic systems and/or model-based locomotion control policies.
Looking ahead, our future work will build upon the foundations laid in this study.We aim to extend this monopedal test platform to a floating-base mobile robot, employ contact estimation algorithms to eliminate the need for external contact sensing measurements, and assess the performance of the SLIP model-based controller with a robot that has optimal passive dynamics specifically targeting the SLIP model.
These future endeavors align with our vision of advancing the field of robotics by exploring innovative mechanisms, control strategies, and applicability in various domains.We believe that the integration of elastic actuation and optimization of passive dynamics will lead to more efficient and versatile robotic systems with potential applications in areas such as search and rescue, exploration, and assistive technologies.

Figure 2 .
Figure 2. Description of the mechanism.

Figure 3 .
Figure 3. Wrapping cam in action, and the definition of variables used.

3. 1 .
Hardware descriptionThe leg mechanism based on a four-bar linkage contains an actuation unit of a GhostRobotics Minotaur robot on the body.It is a T-Motor U10 80 kV 1200 W brushless DC motor with an embedded motor driver unit.The motor driver has an EtherCat communication interface, which can be controlled in torque and position control modes.The leg links are 3D printed with ABS plastic material with an FDM printer.The links have an eye-fork type end structure, interconnected with bearings from two sides to construct revolute joints with low friction and minimized moment loads.The mechanism contacts the ground via radial ball bearings fitted on the leg to have a point contact convenient to the considered dynamic models.The body of the mechanism is constrained in the vertical direction by a linear bearing that connects the body to a guide rail, and the toe is constrained with a 3D-printed guide located on the ground.The mechanism utilizes a wrapping cam designed with the previously explained method to have a stiffness ratio of K = 0.1.The linear position and velocity are obtained by a 3D-printed rack and pinion mechanism equipped with a 2000 CPR encoder that provides 0.09 mm resolution.A contact switch is placed at the bottom to identify the stance-flight transitions.Images of the leg and the experimental platform are shown in figure4.The overall communication structure is depicted in figure5.The system is controlled with Simulink Real-time [31] running on a x86 target computer.A Copley-Accelnet servo motor driver module is used as the data acquisition device by collecting auxiliary sensory input from the contact switch

Figure 4 .
Figure 4.The top image shows the experimental platform, while the bottom one shows the leg mechanism.The red arrows indicate the important parts.

Figure 5 .
Figure 5. Schematic of overall communication structure.

Figure 6 .
Figure 6.Example of continuous hopping trajectory obtained with harmonic forcing with T = 2 Nm and P = 0.25 s.

Figure 7 .
Figure 7. Single-stride data for parameter estimation.Dashed lines separate different phases.(a) COM position trajectory.(b) COM velocity trajectory with a focus on impact transitions.Raw and filtered velocity profiles are shown in blue and red, respectively.

Figure 8 .
Figure 8.The static force measurements of equivalent vertical spring along with the best fit line.The blue circles show measured force, and the red line shows the best fit line.R 2 = 0.9957.

Figure 9 .
Figure 9. Validation of overall trajectory through NMSE of the states.Bottom and top edges of the box indicate the 25th and 75th percentiles, respectively.The whiskers extend to the most extreme data points.The red '+'s show the outliers.Horizontal red line indicates the median of the distribution.

Figure 10 .
Figure 10.Sample trajectories with different validation error for the overall motion trajectory of the linear model.Blue: experiment, red: simulation.(a) Mean Absolute Error (MAE) in COM position: 2.7 mm; MAE in COM velocity: 49.6 mm s −1 .(b) MAE in COM position: 4.7 mm; MAE in COM velocity: 89.1 mm s −1 .

Figure 11 .
Figure 11.Comparison of apex-to-apex performance of the models with and without the impact dynamics.MWI: Model with impact, MW/OI: Model without impact.Bottom and top edges of the box indicate the 25th and 75th percentiles, respectively.The whiskers extend to the most extreme data points.The red '+'s show the outliers.Horizontal red line indicates the median of the distribution.
) and (b) show the distribution of the mean absolute errors of these two models predicting the apex height and apex time, respectively, in 60 singlestride experiments.The results show that accurately estimated impact parameters significantly improve the predictions of both the apex position and apex timing.

Table 1 .
Validation of results for each phase using different error metrics.σ is the standard deviation.