The generalized spring-loaded inverted pendulum model for analysis of various planar reduced-order models and for optimal robot leg design

This paper proposes a generalized spring-loaded inverted pendulum (G-SLIP) model to explore various popular reduced-order dynamic models’ characteristics and suggest a better robot leg design under specified performance indices. The G-SLIP model’s composition can be varied by changing the model’s parameters, such as ground contacting type and spring property. It can be transformed into four widely used models: the spring-loaded inverted pendulum (SLIP) model, the two-segment leg model, the SLIP with rolling foot model, and the rolling SLIP model. The effects of rolling contact and spring configuration on the dynamic behavior and fixed-point distribution of the G-SLIP model were analyzed, and the basins of attraction of the four described models were studied. By varying the parameters of the G-SLIP model, the dynamic behavior of the model can be optimized. Optimized for general locomotion running at various speeds, the model provided leg design guidelines. The leg was empirically fabricated and installed on the hexapod for experimental evaluation. The results indicated that the robot with a designed leg runs faster and is more power-efficient.


Introduction
In nature, most creatures can move agilely in their habitat after evolving for millions of years, even when their surroundings are uneven, complex, or threatening.Researchers are inspired by the fantastic mobility and adaptability of legged animals to improve the motion of robots.For running gait, scientists found that a creature's motion can be widely described by the spring-loaded inverted pendulum (SLIP) model [1][2][3][4].Full and Koditschek then proposed the 'template and anchor' concept, which stated that the complicated motion of an animal (i.e. the anchor) can be extracted by a reduced-order model (i.e. the template), and the model can suggest the control strategy and excite the locomotion of the anchor system [5].Since then, researchers have developed various dynamic models for legged locomotion [6,7].
The proposed reduced-order models have distinct configurations and dynamic properties.The SLIP model [1,2], composed of a point mass and a linear spring, is the most widely used.Another well-known one is a two-segment leg (TSL) model consisting of two straight bars connected by a torsion spring, which yields nonlinear force and displacement relations between the mass point and the ground [8].Based on the SLIP model, Seipel and Holmes considered clock-torqued control on the leg and developed a clock-torqued SLIP model [9].In the meantime, researchers tried to add torque and damping terms to the SLIP model (such as the TD-SLIP model [10,11] and the Hip-SLIP model [12]), and the model was utilized on compliant terrain as well [13].Jun et al also reported a torque-driven and damped TSL (TD-TSL) model [14].In addition, researchers have tried to address the effects of rolling contact.Ankarali et al proposed a C-Pod model [15].Jun and Clark studied the rolling effect by adding a circular foot to the SLIP (SLIP-RF) model as well as torque-driven and damped components (i.e.TD-SLIP-RF model) [16].The two-segment form was also adopted to simulate the rolling effect (i.e.TD-HCL model) [17].With this series of reduced-order models, Jun and Clark further studied the dynamic properties of the models [18].Rao et al compared the actuated two-segment model with a knee spring and the actuated SLIP model and found that they had similar movements in which the two-segment model had slightly higher stability and the SLIP model had a slightly lower energy cost.However, the force acting on the mass of the two-segment model is more complicated-it differs from that of the SLIP model because the hip torque affects the force along the line from the hip to the ground contact point and the force perpendicular to it [19].Descriptions of different extended SLIP models can be found in [20].The use of models to control legged robots is also reported [21][22][23][24][25][26].In addition to single-leg models, bipedal models and their controls have been widely studied [6,7,27].
Previously, we studied reduced-order models and their relation to the dynamic motions of empirical robots on several occasions.A rolling SLIP (R-SLIP) model [28] was developed to serve as the motion template of an RHex-style [29] robot so that the robot runs stably according to the passive dynamics of the R-SLIP model.Later, we proposed a clocktorqued R-SLIP model [30], which added an active leg control to match the empirical control scheme of the robot.The model successfully served as a template for transient, steady-state, and variable-speed running by a robot.A torque-actuated, dissipative R-SLIP model was developed to analyze a robot's energy flow and transformation [31].A physics-data hybrid motion template was designed to bridge the unmodeled discrepancy between the R-SLIP model and the robot [32].In addition, an eccentric SLIP model was developed to model the differentiated leg morphology and function of the fore and hind legs [33], and a three-dimensional rolling template was designed to model the rolling dynamics of a hexapod robot [34].
While various dynamic models have been proposed and analyzed, the models' behaviors and characteristics have not been studied and compared at the quantitative level, especially their fundamental differences in morphology and composition, such as ground contact style (rolling vs. point contact) and leg morphology (single vs. multiple segments).This work proposes a new generalized spring-loaded inverted pendulum (G-SLIP) model parametrized by seven parameters.By adjusting the parameter values, the morphology of the G-SLIP model can be transformed or approximated into four distinct and widely discussed models, as shown in figure 1, including the traditional SLIP model, the SLIP-RF model, the R-SLIP model, and the TSL model.Thus, the characteristics of these models can be analyzed and compared using the same methods and approaches, yielding unified validation results.To the best of our knowledge, this is the first work to utilize a single model to simulate the behavior of multiple-legged models and compare their performance on the quantitative level.As the first exploration of this approach, this work focuses on how intrinsic morphology differences affect the model's dynamics.By varying the parameters of the G-SLIP model, we sought to find the optimal model morphology for a motion template and suggest a suitable leg morphology for the robot.Furthermore, a novel robot leg based on the optimal model was fabricated and tested to validate the model's performance in the real world.
The remainder of the paper is organized as follows: section 2 reports on the construction of the G-SLIP model method; section 3 presents the dynamic characteristics of the G-SLIP model, including those of the four widely-used cases-the SLIP model, the TSL model, the SLIP-RF model, and the R-SLIP model; section 4 describes the optimization of the G-SLIP model; and section 5 discusses the experimental validation of the optimal design.Section 6 concludes the work.

The G-SLIP model and the analysis methods
The analysis of the model in this work involves two steps.In the first step, the conservative G-SLIP model is utilized to find the fixed points, because we would like to use the motion of the conservative model as a template-a reference to initiate the passivedynamic motion of the empirical robot as the anchor.However, because empirically the robot leg has inertia and may be subjected to energy loss, the robot leg is usually controlled.In the second step, to take empirical conditions into consideration, the G-SLIP model is revised to add clocked torque.The basins of attraction (BOA) is analyzed using the clocked-torque G-SLIP model because its motion would be closer to the empirical robot's behavior.

The G-SLIP model
The G-SLIP model was developed to explore the general dynamics of the legged model, such as the effect of rolling contact and spring configuration on dynamic motion.The model has seven intrinsic parameters, including mass (m), upper bar length (l 1 ), lower bar length (l 2 ), radius of the circular rim (r), torsion spring stiffness (k t ), natural angle of the torsion spring (ϕ 0 ), and the angle between the lower bar and the rim (ψ ), as shown in figure 2(a).Details to be noted are: (i) the angle ψ is fixed during motion since the lower bar and rim form a rigid body.Thus, when the circular rim rolls, the lower bar rolls at the same angular speed.(ii) The point mass can be out or in the abstracted circle formed by the rim.(iii) The line connecting the other ends of the upper and lower bar (i.e.l 4 ) does not necessarily pass the center of the circular rim.(iv) The upper bar, the lower bar, and the circular rim are massless.In other words, the model only contains a point mass.As shown in figure 2(a), the model has three touchdown states (i.e. the initial conditions), including touchdown speed (v), touchdown angle (α), and landing angle (β).The touchdown angle (α) is defined as the angle included by the horizontal line and the direction of the touchdown velocity vector.The landing angle (β) is the angle included by the horizontal line and the line segment from the point mass to the center of the circular rim.Thus, the β defines the initial and relative configuration between the point mass and the center of the circular rim.
The motion of the G-SLIP model in stance phase can be described by two generalized coordinates, q = θ ϕ T , as shown in figure 2(a), representing the configurations of the leg and of the torsion spring.
Their initial values (θ 0 ,ϕ 0 ) are the angles when the model touches the ground.The ϕ 0 is the natural angle of the torsion spring as described before.The θ 0 is determined by the model parameters and the landing angle (β), where various assisting parameters are utilized to ease the derivation of θ 0 as shown in figure 2(b).Based on the trigonometric relations, the parameters l 3 and η can be derived as Since the intrinsic parameters (l 2 , r, ψ ) are constant, l 3 and η are constants as well.The θ 0 can then be represented as: where l 0 can be computed as: The l 0 represents the distance between the point mass and the center of the circular rim when the torsion spring is at its natural configuration (i.e.ϕ 0 ), so the subscript 0 is utilized to emphasize this fact.The l 0 is a constant as well.
The equations of motion (EOMs) of the G-SLIP model was derived using Lagrangian method.In stance phase, the fore/aft and vertical displacements of the point mass (x, z) can be expressed as functions of the generalized coordinates (θ, ϕ ): The Lagrangian of the G-SLIP model, L, can be expressed as a function of the generalized coordinates: Thus, the EOMs of the G-SLIP model in stance phase can be derived by the Lagrangian method: where Π and ∆ represent the energy supplied into the system and the Rayleigh dissipation function, respectively.Note that for the G-SLIP model, Π and ∆ are both zero due to the conservative condition.Given touchdown states, (v, α, β) as shown in figure 2(a), the motion of the G-SLIP model in stance phase can be numerically computed.After the model takes off, it enters its flight phase.The EOM of the G-SLIP model in flight phase is a simple ballistic motion equation.The switching events of the stance and flight phases of the G-SLIP model are the touchdown and liftoff events, as shown in figure 1(b).The touchdown happens when the model flies downward till the height of the point mass equals the vertical length of the leg: After the touchdown, the model enters the stance phase.The point mass is going forward and downward, and the torsion spring starts to be compressed.After the maximum compression, the spring gradually returns to its natural configuration (i.e.ϕ 0 ), and the model lifts off: Note that the spring would not exceed the natural configuration (i.e.extension) because the lower bar and the rim are massless.In addition, the model usually has upward velocity when lifting off.In some extreme cases, the model may move forward and downward when the leg rotates and lifts off, so the upward velocity (i.e.ż > 0) is not the necessary condition for the liftoff event.The model in these extreme cases is most likely to fail in the next touchdown, and the failure check of the model can screen these unsuitable gaits.
Finally, with the EOMs of the G-SLIP model in both stance and flight phases, as well as the defined switching events, the dynamic motion of the model can be continuously simulated.Because the stance dynamics is nonintegrable, the analytical solution of ( 7) is unavailable.Thus, a numerical solution is derived.

Fixed points of the G-SLIP model
Similar to the analysis reported in our previous work [28], this work also utilizes a Poincaré map to analyze the G-LSIP model's stability.The touchdown event is selected as the Poincaré section.The touchdown conditions include three states (v, α, β), so in general the map is three-dimensional.Here, the landing angle (β) is set the same for all touchdowns because only the gait with period one is considered in the work.This setting is also commonly utilized in the analysis of SLIP-like models [3,9,12].In this case, the touchdown speed (v) also remains the same because the G-SLIP model is conservative.Therefore, only the one-dimensional Poincaré map of α is required to be analyzed.
The G-SLIP model is considered to operate at a fixed point α * if In this one-dimension case, the stability of the fixed point can be evaluated by the slope.If the absolute value of the slope is smaller than 1, the fixed point is stable.In contrast, it is unstable if the value is larger than 1.Thus, by this method, the distribution of the stable/unstable fixed points of the G-SLIP model with respect to the touchdown conditions (v, α, β) can be obtained.
When running at a fixed point described by a specific set of touchdown states (v, α, β) fp , the model exhibits periodic and alternating stance motion and flight motion.This dynamic motion can be numerically found using the stance EOMs and flight EOMs.More specifically, the generalized coordinates of the model running at the fixed point versus time, θ fp (t) , ϕ fp (t) stance , can be computed.Note that θ (t) represents the trajectory of the model leg, and in empirical implementation, this trajectory is controllable.For easy implementation effort, a fifth-order polynomial is utilized to parametrize the fixed-point stance trajectory of the model leg, and a trapezoid trajectory is implemented while the model leg is in the aerial phase.The combined leg trajectory in one period (i.e. one stance and one aerial) is hereafter referred to as θ fp (t).The parametrized fixed-point trajectory of the conservative model, as the template, can be served as the motion reference of the empirical robot, as the anchor, to initiate its model-like dynamic motion.The methodology can be found in our previous work [28].

The clocked-torque G-SLIP model
Because this work also aims to implement the optimized G-SLIP model on the robot, empirical conditions must be considered.While the massless model leg can move passively in stance phase and can be instantly reposed in aerial phase without effort, the robot leg needs to be controlled to compensate for the leg inertia effect and the energy loss during the motion.Thus, in empirical implementation a clocked-torque controller was added at the hip of the model to control the leg motion, the same strategy as reported in [30].In this case, the EOMs of the clocked-torque G-SLIP model in stance phase can be rewritten as where k P and k D are the proportional and derivative gains.Because of the added clocked torque, the clocked-torque G-SLIP model is not conservative.It does not contain damping terms, as reported in [9,10,14,[16][17][18]35].The energy flowing in and out of the model is utilized to accelerate and decelerate leg motion.
The clocked-torque G-SLIP model regulates leg motion in both the stance and aerial phases.The model leg continuously moves, using θ fp (t) as a reference.If the model leg moves perfectly according to the fixed-point motion, the leg will touch down with the designated landing angle β.The controller has no effect if the model leg moves perfectly according to the fixed-point motion.In contrast, if the leg does not touch down at the desired timing or is subjected to disturbance, the landing angle differs from β.To accurately parametrize the touchdown states with this possible variation of landing angle, the relative phase (ρ) is defined and utilized: where ∆t TD and T represent the time offset of the touchdown timing and the period of the fixedpoint trajectory, respectively.In such situations, the clocked-torque G-SLIP model moves differently from the fixed-point motion of the G-SLIP model.In these cases, additional torques are applied, and the model generally converges its motion to the fixed-point motion.
In summary, the clocked-torque G-SLIP model can be regarded as a modified version of the conservative G-SLIP model.The clocked-torque model utilizes the realistic leg trajectory as a reference (i.e.fixed-point trajectory in stance and time-required aerial trajectory), as opposed to the traditional setting where the leg in the aerial phase can be reposed instantly.In addition, the added clocked torque resembles the setting in the empirical robot's controlled leg.As for the fixed-point trajectory, using the one from the conservative G-SLIP model is desired because the idea of a template and anchor is to make the robot move like a passive model.

The basin of attraction of the clocked-torque G-SLIP model
This work also investigated the BOA to further study the performance of the model.The BOA of the model at a specific operation point is the collection of all the touchdown conditions of the model, whose motion converges to that fixed point.If the model touches down exactly according to the conditions of the fixed point, then the model will follow the motion of the fixed point periodically and precisely, without any deviation.In contrast, if the model touches down according to conditions that differ from those of the fixed point but within the BOA of that fixed point, the model's motion will gradually converge to the motion of the fixed point.Thus, the plots of BOAs shown in this paper illustrate the collection of touchdown conditions that are different from that of the specific fixed point but in which the model's motion can converge to the motion of the fixed point.Both the G-SLIP model and the clocked-torque G-SLIP model have BOAs.Because the added clocked-torque controller allows for more realistic mapping between the model and the empirical robot, the clocked-torque G-SLIP model was utilized for BOA analysis.Our previous results reported in [30] show that the added clocked-torque controller enlarges the BOA in comparison with the original conservative model, which is beneficial to the model's stability.
The construction of the model's BOA at a specific fixed point requires investigating all neighborhood touchdown conditions.The touchdown of the model is parametrized by three states (v, α, β).As mentioned in section 2.3, the leg of the clocked-torque G-SLIP model is controlled to follow the fixed-point trajectory θ fp (t).When the leg does not touch down perfectly, all three states likely differ from those of the fixed point.Thus, the BOA should be parametrized by three variables, and it can be represented using a three-dimensional plot.Because the relative phase (ρ) can present deviation better than the landing angle (β), the touchdown states (v, α, ρ) are utilized.
The BOA of the clocked-torque G-SLIP model at a specific fixed point is computed using the step-to-fall method.If a specific set of touchdown states (v, α, ρ) converges to the designated specific fixed point within specified steps and with an error (ε) less than 1%, the touchdown state is regarded as in the BOA:

Four special cases of the G-SLIP model
As shown in figure 1, the four models-SLIP, SLIP-RF, TSL, and R-SLIP-can be regarded as special cases of the G-SLIP model.Two issues should be addressed for this shape adaptation: first, transforming rolling contact to point contact, and second, simulating linear spring using the original torsion spring.
The transformation from rolling contact to point contact is straightforward.Because the rolling effect is parametrized by the parameter r, when r = 0, the rolling contact becomes point contact.The maximum r = 0.075 m is set according to the dimension of the robot leg.Two more values are evaluated in the following analysis, r = 0.025, 0.05 m.Note that pure rolling is assumed, so the contact point does not slip.
Simulating a linear spring effect using the original torsion spring of the G-SLIP model requires some effort.As shown in figure 3, when the lengths of the bars connecting the torsion spring (l 1 , l 2 ) are great in comparison to L, the spring effect between the other ends of the bars approximates to a linear spring.The relation can be presented as: where the subscript 0 indicates the initial/natural condition.The equation (14) indicates that the linear spring is a function of (l 1 , l 2 , k t , ϕ 0 ).Empirically, the maximum deformation of the torsion spring in the operating range of the running hexapod robot is about 20%-30% of its leg length.In addition, other literature also reveals that the maximum deformation of animals' legs in running and hopping is about  20%-30% [4,36].As a result, the equivalent linear stiffness was approximated using regression, where 0%-20% deformation and corresponding force data were taken.Thus, in addition to the original torsion spring behavior, this setting allows the G-SLIP model to simulate the linear spring behaviors.
To parametrize the variations between the linear spring and the torsion spring, a parameter σ is defined, as shown in figure 4. The quantitative definition of sigma is the relative position of the torsion spring on the red line vector between the point on the rim (i.e.ϕ 0 = 66 • ) and the point formed by two line segments l 1 = l 2 = 1 m, as shown in figure 4. When σ = 0, the torsion spring is located on the rim of the leg, and this is the configuration for the R-SLIP and TSL models.The elastic behavior of the torsion spring included by two short line segments yields nonlinear translational elastic behavior between the mass and the ground, as shown in figure 3.In contrast, when σ = 1, the torsion spring is included by two long bars with lengths l 1 = l 2 = 1 m, which yields approximately linear translational elastic behavior.Therefore, this setup is suitable for simulating the linear spring effect of the SLIP and SLIP-RF models.In addition to σ = 0 and σ = 1, σ = 0.33, 0.66 are also chosen for analysis.The parameter σ can also be regarded as the linearity of the elasticity between the mass point and the end point of the lower bar.The spring is more linear as σ moves closer to 1, and it is more nonlinear as σ moves closer to 0. Note that for each spring position, the stiffness k t should be recomputed so that the equivalent linear spring (or the empirical leg compliance) is set in a similar range for comparison.
With the described equivalent methods, the four models can be successfully transformed from the G-SLIP model.That is, the four models can be regarded as the special cases of the G-SLIP model with specific parameter settings.Table 1 lists the quantitative parameter values of these four models using the parameters of the G-SLIP model.The derivations of the EOMS of these models are included in the appendix.

Dimensionless settings of the G-SLIP model
To discuss the model's characteristics more generally, a dimensionless G-SLIP model was utilized in the analysis.The dimensionless parameters and touchdown states are listed below.
where mass (m), the maximum leg length from mass point to lower bar edge (l 0 + r), and gravity acceleration (g) are used as the standard.Table 2 lists the dimensionless parameters of four special cases of the G-SLIP model: SLIP, SLIP-RF, TSL, and R-SLIP.

Dynamic characteristics of the models
In this paper, we explore and compare the dynamic characteristics of models in different morphologies.
The G-SLIP model offers a standard model structure as the basis for transformation into four commonly used models (i.e.SLIP, SLIP-RF, TSL, and R-SLIP models, as shown in figure 1) using parameter adjustment, as described in section 2. In addition, by interpolating the parameter space, the change in dynamic characteristics owing to the parameter change can be clearly observed, which greatly helps to determine the effect of each parameter.As described in section 2, the fixed point and the BOA are utilized to evaluate the model's performance.

The fixed-point distribution of the model
The fixed points of the G-SLIP models with r = 0, 0.025, 0.05, 0.075 m and σ = 0, 0.33, 0.67, 1 were evaluated, and figure 5 shows the fixed-point distribution of the 16 models evaluated.Based on the parameter setting, the models at the four corners of the figure represent SLIP, SLIP-RF, TSL, and R-SLIP, respectively.Note that the fixed-point distribution of the original SLIP model and that of the G-SLIP model still have some discrepancy caused by spring equivalence, and this discrepancy is further amplified in the SLIP-RF model owing to rolling motion.The discrepancy occurs in large α and v regions, where the model rolls for longer.Because the discrepancy between the original SLIP/SLIP-RF model and the G-SLIP model in the commonly used region is tolerable, the G-SLIP model in the following contents is utilized for formulation consistency.
The models analyzed here are dimensionless, so the results are informative for models with the same dimensionless parameters.Several trends can be observed: (i) As the parameter r increases, the fixed points gradually extend to a low landing angle area and cross over to the right side of the α = β line.This is because rolling backward slightly can

The BOA of the clock-torqued models
To further explore the dynamic characteristics of the G-SLIP model, the BOAs of the model parametrized in four commonly used models shown in figure 1    cooperate closely, so the BOAs are smaller.In contrast, the BOAs of the other three models running at 25 fixed points are more evenly distributed.
Note that in this analysis, the G-SLIP model was utilized, rather than the four original models.Thus, the linear spring was approximated by the torsion spring with large arms.In this case, the two DOFs are not completely independent.This helps to improve the performance of the model.Even with clocktorqued hip actuation, the SLIP model can barely survive the disturbance, and a damper parallel to the spring is required to stabilize the nonpassive motion, as reported in [9].In addition, an actual robot leg of the G-SLIP type can be reproduced by implementing a torsion spring between two leg segments, which would be similar to [31].

Model optimization and results
In addition to comparing the dynamic characteristics of various popular reduced-order models on a consistent basis, the other purpose of constructing the G-SLIP model was to evaluate the optimized form of the model that yielded the desired dynamic response.Owing to the parametrized and dimensionless model structure, optimization can be achieved by searching for the best parameter set of the G-SLIP model.To fairly justify the parameter space (i.e.dimensionless setting), the models in optimization have the same mass m and maximum leg length l 0 + r.Thus, the five DOFs of the dimensionless G-SLIP model can be varied.Here, let κ be the angle between upper bar l 1 and l 0 , and ξ be the angle between l 3 and r.Five parameters (k t , l 1 , κ, ξ , r) were selected as the variables of optimization, as shown in figure 11.The searching space was formed by these five parameters, and their ranges were bounded based on the reasonable physical and geometric ranges of the model and the robot.Table 3 lists the upper and lower limits of these five variables.Note that these parameters are transformed into dimensionless ones in the optimizing process.Moreover, the cost function can be defined based on the different demands described in the following subsections.The optimization in this paper was executed using the fmincon function in MATLAB®, where an interior-point algorithm was utilized.To yield a global-like optimal solution, 3000 initial points were randomly and evenly generated in the searching space, and 100 of them with the lowest costs were selected as the first-step screening.Next, these 100 points were utilized as the starting points for the optimization process, and the one with the lowest cost was regarded as the global optimal solution.

Optimization of the model running at a single fixed point
In general, the model can run at many different fixed points.As a first-step optimization, running at a single fixed point is utilized, and the fixed point (v, α) = 2 m s −1 , 10 • is selected because it was widely utilized in our past robot experiment.The optimization is set toward two goals: faster forward speed and larger vertical displacement, and the cost functions are: cost maximum forward speed = −ṽ x avg (17) cost maximum vertical displacement = −d H (18) where ṽx avg is the mean of dimensionless forward speeds in one step and d H is the difference between maximum and minimum vertical displacement in one step.In addition, the stance phase and gait period of the model are set in the range of 0.1-0.5 s and 0.3-0.7 s, respectively, matching the time spans of running animals with similar weights [37].
Figures 12(a) and (b) plot the running motions and configurations of the top five models using speed cost and displacement cost, respectively.Except for the first result in 12(a) and the fourth one in 12(b), the other eight models have similar configurations and dynamic motions.For instance, the models tend to have a similar angle of the upper bar at touchdown.
Since in each cost setting the costs of the top five optimal models are close to each other, here the second model with the fastest forward speed cost and the first model with the largest vertical displacement cost are taken for further discussion.Table 4 lists the optimized parameters of these two models, and figure 13 plots the running motions and configurations of these two models.Several observations can be derived: (i) Note that by changing the touchdown conditions, the model can exhibit different flight behaviors for a long stride length or large ground clearance.Here, the common touchdown condition for running (v, α) = 2 m s −1 , 10 • as the fixed point (i.e. the constraint) was utilized intentionally for optimization.In this case, how high the model can fly or how fast it can run were determined by merely changing its geometrical and elastic parameters.Different parameter settings yield different liftoff conditions, so the resulting flight behavior differs.(ii) As shown in both figures 12 and 13, regardless of whether the model is optimized for the fastest speed or the highest displacement, the models tend to have long arms for torsion spring.This indicates that a linear spring is preferable.(iii) As shown in both figures 12 and 13, regardless of whether the model is optimized for the fastest speed or the highest displacement, the models tend to have large circular rims.This indicates that rolling contact is preferable.It is speculated that rolling can not only assist in increasing forward displacement in stance phase but also mitigate the vertical impulse.(iv) As shown in both figures 12 and 13, regardless of whether the model is optimized for the fastest speed or the highest displacement, the trajectories of the mass seem straight compared with those of the standard SLIP or R-SLIP models.We believe the straight-line-like stance trajectory of the optimized model results from the combination of geometrical configuration, rolling effect, and close-to-linear spring effect.
(v) In the case of fastest forward speed, the trajectory of the mass point has relatively smaller vertical oscillation and energy conversion is less.This is because most of the energy is utilized to run forward, and less energy is stored in the torsion spring.To achieve these dynamics, the model has larger spring stiffness and slightly shorter bar length than the model optimized for the maximum vertical displacement.(vi) In contrast, in the case of largest vertical displacement, the model has larger energy transfer among kinetic energy, gravitational energy, and spring potential energy.Smaller spring stiffness is required to provide larger displacement and energy storage.(vii) The period of the model optimized for the fastest forward speed is relatively short, while that for largest vertical displacement is long.Furthermore, the ratio of the flight phase to period of the former model is about 0.52, larger than that of the latter case, 0.47.The results verify that most progress occurs in flight phase.

Optimization with combined indicator
For empirical implementation, the robot is expected to run at various speeds (i.e.different fixed points), so the optimization of the model should consider using multiple fixed-point motions to yield the optimal model parameters.Thus, a model optimized for a combination of three indicators was developed.The first indicator is gait variability, which utilizes the number of fixed points as the indicator.The 25 commonly used touchdown states v = 1.25, 1.5, 1.75, 2, 2.25 m s −1 and α = 5, 10, 15, 20, 25 • are imported into the model to check whether the model has fixed points at these touchdown states.The second indicator is the mean of the forward speeds of the model running at the fixed points.Forward speed is an important indicator of the running model.In addition, to avoid large oscillations, the vertical displacement of the mass point should be taken into account.As a result, the cost function of the combined indicator can be formulated as follows: where N fixed point represents the number of existing fixed points out of 25 fixed points for a certain model.The subscripts mean and std present the mean value and standard deviation, respectively.The symbols (w 1 , w 2 , w 3 ) are weightings.The three indicators  shown in (19) are standardized and weighted by the weightings, and table 5 summarizes the associated values.
Figure 12(c) plots the running motions and configurations of the top five models using the combined indicator.Table 4 and figure 14 show the associated results of the optimized model.Several trends can be observed in the results: (i) As shown in figure 12(c), the optimal models optimized for the combined indicator have a similar configuration to those optimized for the fastest forward speed and largest vertical displacement, as shown in figures 12(a) and (b): The linear spring and rolling contact are preferable.(ii) Table 4 reveals that the optimal model performs well in forward velocity, has no exaggerating vertical displacement, and has all fixed points assigned, showing well-performed locomotion and gait variability.(iii) Owing to a small vertical displacement setting in the cost function, the spring stiffness is larger, and less compression occurs.Moreover, to optimize the motion at all 25 fixed points, the bar length of the optimal model is much shorter than those in other cases because the shorter bar results in a larger tolerance for different directions of touchdown impact.This is because radial and angular DOF are more dependent when the springs are relatively nonlinear.(iv) The optimal model has fixed points on the given 25 touchdown states.The landing angles of these fixed points are located in the range of about 18 • .This phenomenon is reasonable because all 25 fixed-point motions are optimized for similar criteria (i.e.fast forward speed and less vertical displacement), resulting in similar landing angles.
Note that the original G-SLIP model utilizes the infinite lengths of l 1 and l 2 to approximate the linear spring behavior, so using only one model can simulate behaviors of both linear and torsion springs for analysis.In an empirical implementation, if the optimized model has long l 1 and l 2 , a linear spring can be utilized directly.
The elastic characteristics of the optimized model resemble those of biological creatures and other dynamic robots, where relative leg stiffness is utilized as the index.According to [4], relative leg stiffness is defined as the normalized ratio of peak groundreaction force to peak leg compression.Studies [38][39][40] found that the relative leg stiffness of animals and  robots is within the range of 7-27.Our optimal leg design using the combined indicator has a relative leg stiffness of 18, which is within the observed range.
In this paper, single-objective optimization is utilized for either the fastest forward speed, the largest vertical displacement, or the combined indicator.The Pareto optimal set is a set of optimal solutions for multi-objective optimization; thus, we do not have a separately as objectives, a multi-objective optimization can be conducted, and then the Pareto optimal set and Pareto front can be identified.This is another approach to optimization.

Experimental implementation and results
The analysis of the G-SLIP model not only provides a comparison basis for existing models as described in sections 2-4, but also provides a suggestion for the leg design of the robot.Following this, the G-SLIP leg with optimal configuration was fabricated and implemented on TWIX, the RHex-style hexapod robot [29] in the lab, for experimental evaluation.Note that the G-SLIP model only has a point mass and one elastic leg, but the robot has a rigid body and six legs.Thus, the single-legged model and the hexapod robot need certain mapping, so the template-andanchor strategy can be realized.The model's mass is identical to the robot's mass, but the model's leg stiffness is triple that of the robot's leg stiffness because the robot runs in an alternating tripod gait.Other mapping details can be found in our previous work [28].The TWIX moves its legs just like the RHex by monotonically rotating the legs in the same direction.
To make the leg realization feasible, some modifications were required.Because the relatively long, straight bars of the optimal G-SLIP model were originally designed to approximate the linear spring and the torsion spring compression angles are small (about 3 • ), the linear spring was utilized.The dynamic behavior of the original optimal G-SLIP model shown in figure 15(a) and its approximated version shown in figure 15(b) were compared.Their responses were close to identical, as shown in figure 15(c), and this result confirms that the discrepancy from the approximation was negligible.Note that the configuration of the approximate model is similar to that of the SLIP-RF model but, unlike the latter, the linear spring is not perpendicular to the rim.The G-SLIP leg can also be fabricated with other configurations using the same approach reported in this work.
A set of legs was fabricated based on the leg morphology of the approximated G-SLIP model, as shown in figure 16(a), and installed on the robot TWIX for experiments, as shown in figure 16(b).The leg was made up of a linear spring, a linear slider leg structure made of Polylactic acid, and tire tread, which provides the leg traction on the ground.The robot runs with an alternating tripod gait, the same as in other previous works [28,[30][31][32].The motion trajectory of the robot leg in stance phase follows that of the optimal G-SLIP model.To be more specific, the optimal G-SLIP model serves as the motion template of the robot as the anchor [5].
In the experiments, the robot ran straightly through the seven-meter runway, and the motion capture system, which consists of 10 high-speed  Figure 17 shows the experimental results of the robot running at the selected fixed points.The results are summarized in a statistical manner.The robot's motion fits the trend of the dynamics of the optimal G-SLIP model, indicating that the model (i.e. the template) effectively excites the stable running gait of the robot (i.e. the anchor) using legs designed by the approximated G-SLIP model.Moreover, the goal of the model optimization using the combined indicator is to find the model suitable for faster forward speed and avoid large motion oscillations at 25 commonly used fixed points.To verify the performance of the optimization, the experimental results were compared to those reported in our previous work [30], where the robot ran with a C-shaped leg made with ultra-high-molecular-weight polyethylene (UHMWPE).While the magnitudes of vertical body oscillation are nearly the same, the horizontal speeds of the robot with a G-SLIP leg are higher than those with UHMWPE legs.That is, the optimal G-SLIP leg indeed improves the performance of the running gait as the goal is set.
Other performance metrics were also evaluated.Table 6 lists the root mean squared errors of the robot's pitch and roll, which on some levels maps to the robot's motion stability.The table shows that the values are less than 5 • , enabling the robot to perform stable running.In addition, specific resistance (SR), which represents the averaged energy cost of the system per unit weight moving at a unit distance, is also analyzed: where P is the power and ∆x is the horizontal displacement.Here, the SR is computed based on the   one-step motion of the robot, so ∆x and ∆t represent the displacement and period of one-step motion.Table 6 lists the SR values of the experiments.The SR of the robot with a UHMWPE leg is about 0.6 and 0.4 for low-speed fixed points v = 1.5 m s −1 and high-speed fixed points v = 2.0 m s −1 , respectively.Thus, it is obvious that the SR of a robot with a G-SLIP leg has better motion efficiency because of the cost of fast forward speed in optimization.In summary, the generalized G-SLIP model formulation provides an opportunity to compare the performance of various popular reduced-order models.This formulation also allows the optimization of the model morphology to improve the model's dynamic performance.The optimal model can then serve as the design reference for the robot to yield better performance.Finally, the experimental validation provided in this section confirms the advantages of this design approach.
It is possible to optimize the model for minimum SR using the proposed methodology, but the model itself would require significant rework.The G-SLIP model is conservative, so there is no energy movement in or out of the model.In this work, optimization equates to finding the model's optimal geometrical configuration and elastic property so that it can perform certain explicit behaviors, such as fast running or a continuous forward high jump.If optimizing the SR of the model is desired, the model should contain energy-in-and-out terms.

Conclusion
This paper reports on a generalized model, the G-SLIP model, which was utilized to explore the characteristics of various dynamic models and to suggest a robot leg design.The G-SLIP model can achieve different morphologies by changing its model parameters.The popular SLIP model, the TSL model, the SLIP-RF model, and the R-SLIP model are four special cases in the discussion that focus on the effect of rolling contact and spring configuration on dynamic motion.The analysis reveals that rolling contact can help keep locomotion from failing, and the model with medium or large σ parameter has more stable fixed points.Moreover, the BOAs of the SLIP model, the TSL model, the SLIP-RF model, and the R-SLIP model were analyzed.The model with rolling contact has larger BOAs among fixed points, and the model with two independent degree of freedom has significantly larger BOAs at extreme fixed points than at nonextreme ones.
The parameter optimization of the G-SLIP model yields better locomotion performance given particular goals.The fast forward speed, or large vertical displacement, is the cost for optimization of the model running at a specific fixed point, and the combined indicator is the cost for optimization of the model running at general multiple speeds.The process initially utilizes random deployment and then optimization.The results show that for each cost function, the top optimal models have similar configurations.As optimized for the fastest forward speed, the model's energy conversion to the spring is less than that for the largest vertical displacement.As optimized for the combined indicator, where the model is suitable for running under various conditions, the optimal result provides fast running gaits at commonly used fixed points and enough gait variability.Generally, a model with a medium or long bar length and rolling contact is preferable, and this informative result can serve as the leg design guideline.
This work was also dedicated to manufacturing the robot's leg based on the optimal model using the combined indicator and then performing experiments on an empirical robot with the designed legs.The torsion spring with two long bars of the optimal G-SLIP model is implemented by a linear spring, and both models performed similarly.The experimental results indicate that the forward speed of the robot with the G-SLIP leg is higher than that of the UHMWPE legs.In addition, the former has better energy efficiency (i.e.low SR) than the latter.
In the next step, the stability and BOA of the model can be considered in the optimization.In this case, the computation cost would dramatically increase, and a more efficient optimal method must be developed.In addition, the current G-SLIP model has a parameter space large enough to compose various popular models, but it is not sufficiently generalized to cover the complete morphology space of the planar reduced-order dynamic models.This would require the model to be composed in a general manner.

C. The SLIP-RF model
The four intrinsic parameters of the SLIP-RF model are mass (m), spring length (l 0 ), spring stiffness (k), and radius of the circular rim (r).The relationship between generalized (l, ϕ ) and Cartesian coordinates (x, z) can be expressed as: Thus, the Lagrangian L of the SLIP-RF model would be: Finally, the switching events of phases are listed as follows: where the landing angle β = π /2 + ϕ 0 .touchdown event: z = r + (l 0 − r) sin β (33) liftoff event: l = l 0 (34)

Figure 1 .
Figure 1.The G-SLIP model: (a) its four transformations using model parameter adjustment: SLIP, TSL, SLIP-RF, and R-SLIP model.(b) Running of the G-SLIP model, which is composed of alternating stance phase and flight phase in each stride.

Figure 3 .
Figure 3.The illustrative drawing that shows the approximation of the torsion spring to the linear spring.

Figure 4 .
Figure 4.The setup to parametrize the positions of the torsion spring in the analysis using σ.When σ changes from 0 to 1, the spring configuration and stiffness change from the R-SLIP or TSL model to the SLIP or SLIP-RF model.The red dot represents the position of the torsion spring.

Figure 5 .
Figure 5.The fixed-point distribution of 16 models with various parameters r and σ.In each subfigure, a different color represents a different touchdown speed (v).The circular and cross markers indicate stable and unstable fixed points, respectively.
were investigated.Here, a clock-torqued control was added to the hip of the model.The BOAs of each commonly used model were investigated in 25 variations parameterized by two parameters: v = 1.25, 1.5, 1.75, 2, 2.25 m s −1 and α = 5, 10, 15, 20, 25 • , which are common fixed points for empirical implementations on the robot.The results

Figure 6 .
Figure 6.The 25 BOAs of the TSL model with commonly used fixed points.The x, y, and z axes of each subfigure represent touchdown speed (v), touchdown angle (α), and relative phase (ρ), respectively.The colors represent the number of steps for the model to converge.

Figure 7 .
Figure 7.The 25 BOAs of the SLIP model with commonly used fixed points.The x, y, and z axes of each subfigure represent touchdown speed (v), touchdown angle (α), and relative phase (ρ), respectively.The colors represent the number of steps for the model to converge.

Figure 8 .
Figure 8.The 25 BOAs of the R-SLIP model with commonly used fixed points.The x, y, and z axes of each subfigure represent touchdown speed (v), touchdown angle (α), and relative phase (ρ), respectively.The colors represent the number of steps for the model to converge.

Figure 9 .
Figure 9.The 25 BOAs of the SLIP-RF model with commonly used fixed points.The x, y, and z axes of each subfigure represent touchdown speed (v), touchdown angle (α), and relative phase (ρ), respectively.The colors represent the number of steps for the model to converge.

Figure 10 .
Figure 10.The volumes of the BOAs of the models running at 25 different fixed points, where the models include (a) the TSL model, (b) the SLIP model, (c) the R-SLIP model, and (d) the SLIP-RF model.Each bar represents the volume of the BOA of the model running at a specific fixed point.The bars in different colors represent different touchdown speeds (v) of the model.

Figure 11 .
Figure 11.The five parameters (red) of the G-SLIP model utilized in optimization.The two parameters (black), m and l0 + r, are the same for all models and serve as a dimensional standard in optimization (m = 7.78 kg and l0 + r = 0.15 m).

Figure 12 .
Figure 12.The top five models, which are optimized using (a) the fastest forward speed, (b) the largest vertical displacement, and (c) the combined indicator.Each subfigure presents the optimized configuration of the model at the touchdown moment.The solid red curve and blue dash curves represent the stance and flight phases of the model, respectively.

Figure 13 .
Figure 13.The running motions and configurations of the model optimized for (a) the fastest forward speed and (b) the largest vertical displacement.The corresponding dimensionless energy in one period are shown in (c) and (d), respectively.The solid and dotted curves represent the model in stance phase and flight phase, respectively.

Figure 14 .
Figure 14.The model was optimized using the combined indicator: (a) The configuration of the model, (b) the dimensionless energy of the model in one period, where the model runs at the fixed point (v, α) = ( 2 m s −1 , 10 • ) , (c) the distribution of 25 commonly used fixed points of this model, and (d) the trajectories of the mass point of the model running at these 25 fixed points.

Figure 15 .
Figure 15.The approximation of the optimal G-SLIP model for empirical implementation: (a) The optimal G-SLIP model with a torsion spring and two long bars; (b) the approximated model with a linear spring; and (c) the motion profile of the optimal G-SLIP model (blue) and the approximated model (red).

Figure 16 .
Figure 16.The robot leg was designed based on an approximation of the optimal G-SLIP model: (a) a photo of the leg, and (b) the hexapod robot, TWIX, with the designed legs.

Figure 17 .
Figure 17.The body states of the COM of the robot running at fixed points: (a) (v, α) = ( 1.5 m s −1 , 10 • ) and (b)(v, α) = ( 2 m s −1 , 10 • ).The blue curves and the light blue regions represent the mean and standard deviation of the experimental data.The red curves represent the motion of the optimal G-SLIP model.

Table 1 .
Parameters of four special cases of the G-SLIP Model.

Table 2 .
Dimensionless parameters of four special cases.

Table 3 .
Upper and lower boundaries of parameters of the model in optimization.

Table 4 .
The parameters of the models optimized using different costs.

Table 5 .
The parameters utilized in the optimization process of the model using the combined indicator.
Pareto front in this case.If −

Table 6 .
Averaged pitch, roll, and SR of the robot in experiments.