3LP: a linear 3Dwalking model including torso and swing dynamics
Abstract
In this paper, we present a new model of biped locomotion which is composed of three linear pendulums (one per leg and one for the whole upper body) to describe stance, swing and torso dynamics. In addition to double support, this model has different actuation possibilities in the swing hip and stance ankle which could be widely used to produce different walking gaits. Without the need for numerical timeintegration, closedform solutions help finding periodic gaits which could be simply scaled in certain dimensions to modulate the motion online. Thanks to linearity properties, the proposed model can provide a computationally fast platform for model predictive controllers to predict the future and consider meaningful inequality constraints to ensure feasibility of the motion. Such property is coming from describing dynamics with joint torques directly and therefore, reflecting hardware limitations more precisely, even in the very abstract high level template space. The proposed model produces humanlike torque and ground reaction force profiles and thus, compared to pointmass models, it is more promising for precise control of humanoid robots. Despite being linear and lacking many other features of human walking like CoM excursion, knee flexion and ground clearance, we show that the proposed model can predict one of the main optimality trends in human walking, i.e. nonlinear speedfrequency relationship. In this paper, we mainly focus on describing the model and its capabilities, comparing it with human data and calculating optimal human gait variables. Setting up control problems and advanced biomechanical analysis still remain for future works.
I Introduction
Humanoid robots are challenging to control mainly due to their complex structure and floating base. During locomotion tasks, these systems introduce another complexity compared to wheeled or flying robots, which is the hybrid nature of stepping where the continuous model changes in each phase. Since the beginning, it has always been challenging to balance these robots considering the fact that they can only establish unilateral support with the environment. In addition, creating a sequence of motion, associated timing and the required control architecture in each phase of motion are other important topics in controlling humanoid robots. The main objectives are therefore being humanlike, energy efficient, versatile and of course agile like humans. In this paper, we are proposing a new template model that describes main aspects of walking while being computationally very efficient. Such model can be very useful in modern control architectures from the computational perspective. It can also go beyond conventional template models such as Linear Inverted Pendulum (LIP) by producing more natural motions in faster walking speeds, resembling human locomotion.
The design of controllers should address many concerns like fast implementation, stability, robustness to unknown model parameters and the degree of dependency on sensory data. It has always been appreciated if transition to different speeds and large disturbance rejection are handled in the same control framework. Therefore, candidate methods are normally coming with proper identification of the basin of attraction regarding system states, actuator limitations and violation of model assumptions. This identification is not always straight forward due to its nonlinear nature, though it has been postulated that two steps are enough to stabilize in almost all conditions [1]. In this regard, Model Predictive Control (MPC) is a powerful framework as it can find optimal policies constrained to certain actuation and state limitations. It can also predict if there is no feasible solution, in order to let the algorithm take a different decision.
Ia Hierarchical controllers
Recently, hierarchical control approaches are becoming popular, where a simple template model determines the overall dynamics in an abstract way and then, a detailed fullbody inverse dynamics controller converts this behavior to individual actuator inputs [2, 3, 4]. In dynamical systems, prediction of future evolution is mainly sensitive to the model and sensory data precision [5]. In hierarchical approaches similarly, dynamical matching between the template and the full model is crucial to ensure precise execution of the abstract plan.
[6]  [7]  [8]  [9]  [8]  [10]  [11]  
name  SLIP  IP          LIP 
knee          x     
steering              x 
3D            x  x 
smooth  x            x 
torso               
swing      x  x  x  x   
linear              x 
[12]  [13]  [14]  [15]  [16]  [17]  proposed  
name  BSLIP  FMCH          3LP 
knee          x  x   
steering            x  (x) 
3D            x  x 
smooth  x  x    x      x 
torso  rotating  rotating  rotating  rotating  rotating  rotating  fixed 
swing      x  x  x  x  x 
linear              x 
IB Inverted Pendulums
One of the earliest template models that roughly described human in single support was Inverted Pendulum (IP) [18] with fixed leg length. In this model, a single mass rolls over a contact point, established by a massless leg. IP is widely used to analyze passive walkers [18] and human motion [19]. Inspired by this idea, there are various simple robots built to walk naturally with minimal energy, pumped either in pushoff or swing hip or both [20]. Later this model was simplified to Linear Inverted Pendulum (LIP) [21], mainly favoring availability of analytical solutions instead of numerical integration. With proper modulation of Zero Moment Point (ZMP) [22], many position controlled robots like ASIMO [23] can perform walking through inverse kinematics methods. These algorithms are normally able to produce slow to moderate walking speeds. However complex robots using the LIP method usually walk with crouched knees to keep the Center of Mass (CoM) at constant height. In addition to increasing energy consumption, it is harmful for the robot in long term and less human like, though providing full controllability.
IC Multilink pendulums
Apart from these two models, there are other nonlinear extensions solved numerically. In [9, 8] the IP model is extended to have two separate masses for each leg as well as a single mass at hip. Using similar actuation schemes, this model produces compass gaits on 2Dconstrained simple robots. In [8], same model is modified to have another Degree of Freedom (DoF) in the knee for the swing leg to avoid foot scuffing. The stance leg however always stays straight. In [24], this advanced model is augmented with a torso and later, it is used also by [16] to perform natural walking on uneven terrain using a library of motion primitives. Another model with four masses in legs, hip and torso is proposed in [25] without any DoF in the knees, trying to slightly turn in 3D.
Targeting impactless walking, a simpler model with two passive springs in the hips is proposed by Gomes [15]. Addition of these springs is mainly motivated by elastic properties of human muscles. By exploiting torso motions, Gomes can find zero energy gaits that means impactlosses are less important in energetics of walking. Another very interesting extension of IPbased models is proposed in [17]. In this new model, the pelvis has a width in 3D with a mass in the center. The swing leg is also having another DoF in the knee. This new 3D model is allowed to take advantage of a limited transversal wrench in the contact point for better turning. However finding a periodic gait for such complicated model is difficult and computationally expensive.
ID SpringLoaded Pendulums
It is always questionable which template model produces more realistic motion. This can be inspected from the viewpoint of geometry, torques or energy. The aforementioned categories mainly address energetic and geometric similarities. However ground reaction forces and elastic behaviors are other favorite aspects addressed by another category based on Spring Loaded Inverted Pendulum (SLIP). This model is composed of two massless springs (legs) connected to a single mass. One can expect better description of energy exchange in this model over faster walking speeds and running, mimicking compliant properties of human tendons. Based on this model, Iida [6] built a hipactuated robot walking in 2D with various springs, similar to human muscles. Properties of the passive version without hip actuation were later widely explored in [26]. Using the concept of Virtual Pivot Point (VPP) to stabilize the torso, the model was also extended to have an upper body [13]. In Figure.1, we have briefly shown key models proposed for walking in the literature. Note that in some of them, passive springs are added to the hip actuators for energy storage, similar to humans. In this paper however, we do not investigate elastic behaviors and energy storage.
IE Control difficulty
Except LIP, all other models presented before require numerical integration to obtain time trajectories. Therefore, the Jacobian around a nominal solution linearizes the model and provides the framework for Floquet analysis or discrete controller design [26]. This approach can be used to create an optimal library of primitives [27, 16, 17]. However online reaction to disturbances as well as inclusion of other inequality constraints that are often ignored in calculating a stable basin of attraction, limit the generality of this framework. MPC on the other hand is powerful in this regard, however it requires simple and possibly linear models to provide online performance, depending on the hardware platform. LIP model therefore fits best in the MPC framework [2, 28]. MPC, its simpler version, LQR and sometimes Discretized LQR (DLQR) [29] are used a lot on nonlinear models to stabilize walking gaits and recover pushes [27, 9]. Either a library of optimal policies is generated offline or a discrete transition model is considered at specific events like CoM apex or heelstrike.
IF Why 3LP?
In this paper, we propose a more general version of LIP [21] with three linear pendulums (called 3LP) that captures torso and swing dynamics in 3D. This model allows prediction of future at any time in closed form which is favorable by limited computational resources and MPC. Compared to LIP, the planned motion with 3LP has instantaneous accelerations that are easier to track by inverse dynamics block in the hierarchical control. In other words, the motion of CoM is more natural for the humanoid robot as swing and torso dynamics are taken into account. Additionally, the swing trajectory is more natural compared to many other template models that track an imposed angle of attack with a stiff controller [27, 9, 20]. It should be noted however that the CoM height in our model will still be constant similar to LIP.
The 3LP model supports many different inputs, i.e. hip and ankle actuation authorities. These input dimensions let us find various types of gaits with simple algebraic routines, not optimizations. 3LP as a template model is in fact very useful for motion planning. Besides describing falling dynamics like IP [19] and LIP [30], hip torques required to keep the torso upright are also part of the model like [12]. The most outstanding feature of 3LP is considering swing dynamics that allow us to calculate natural cycles for the model. Considering Figure.1 again, there are few models in the literature that consider this integral part of walking. In these models due to nonlinearity, numerical integration is always needed to search for periodic gaits. In 3LP however we do not need to impose the timing, nor optimizing hip torque trajectories.
This paper merely focuses on introducing 3LP and its capabilities, while setting up control problems remain for future work. The main motivation behind 3LP is for control however, as it can provide more natural CoM and swing motions. It suits MPC control better than LIP in our previous work [2], because we do not need to define a desired footstep plan. The plan comes out of 3LP just by modulating hip torques. Closed form formulations also let the MPC controller change phase timings in case of extreme push recovery for example. Although the MPC problem becomes nonlinear in this case, one can define meaningful hip torque magnitude and rate limitations instead of putting vague timing boundaries on the step time which do not precisely reflect physical facts about the real hardware.
In the next section, we will explain the model details and assumptions behind as well as compare it to the existing models. Next, a method to find different periodic gaits is explained based on symmetry ideas. In the end we show that for a humanlike gait, actuation profiles of 3LP are similar to those of human. Additionally, we also show quantitatively that 3LP can explain main optimality criteria of human walking, i.e. speedfrequency relation, despite being linear.
Ii 3LP dynamics
Motivated by the fact that CoM motion is influenced by swing and torso dynamics, we have added two other pendulums to the normal linear inverted pendulum model, connected together with a pelvis of certain width. In this model, as shown in Figure.2, there are 2DoF actuators in each hip and ankle. We demonstrate feet of limited size in Figure.2 merely to mention availability of ankle actuation. The upper body is represented by a single mass referred to as torso and each leg is represented by a single inertialess mass. By construction (assuming ideal controllers), masses stay in horizontal planes of constant height and the torso is always upright without transversal rotation. These assumptions are used in [17] as well to decouple sagittal and lateral dynamics. Since the torso is connected to an accelerated frame (i.e. pelvis), the torque is always found to keep the torso upright. Such hip torques in literature are alternatively calculated by virtual pendulum concept in the models which let the torso pitch or roll freely during the motion [13]. Inspired by the fact that rolling contact constraint produces more human like walking [31], we also allow for transversal wrenches at stance contact to keep the pelvis orientation fixed. In 3LP model, we do not consider turning properties as they make the model nonlinear. It is however practically easy to turn the robot using inverse dynamics layer, as we showed in [2] using a simple LIP model. There is no need to add heelstrike and pushoff impulses since our double support phase smoothly takes care of transition. In SLIP based models [13], the compliant springs automatically produce a double support phase and perform weight transition without impulses. In 3LP however, switching to double support is triggered when both horizontal components of the swing foot velocity are zero. This assumption is typically used in models with swing dynamics and double support phase like [15], though unlike 3LP, Gomes removes double support for more simplicity.
In Figure.2, all external/internal forces and torques as well as positions are shown for each interesting point of the model, i.e. contacts, hips and pelvis center. These vector variables in our model are expressed in Cartesian frame. One can easily write geometric relations as:
(1) 
Where stands for pelvis width and , depending on left or right support. We can write total force equations for each mass i=1,2,3:
(2) 
and total moment as:
(3) 
Finally, we can write these equations for our massless pelvis around the center point:
(4) 
which link all other variables together. In these equations, we consider as independent variables and solve for others as dependent variables. Actuation possibilities for 3LP are selected as stance foot and swing hip torques in sagittal and lateral directions. Note that the variables stand for disturbances added to our model to simulate external forces. Now, the differential equations governing the system behavior could be obtained for different phases. Features of a full stride, consisting of a double support followed by a single support are defined in Table.I. In double support, the weight is transfered from to and in single support, the leg with variables of subscripts 2 will swing forward. We are going to find transfer matrices that relate the states of the system together at any instance of time over a stride phase.
Full stride =  double support +  single support 

duration  
timing order  1  2 
stance leg  subscripts 3  subscripts 3 
swing leg  X  subscripts 2 
control input  ankle  hip/ankle 
controllability  redundant  full 
Iia Single support
In this phase, the swing foot does not have any external forces, meaning:
(5) 
Also stance foot is fixed on the ground, meaning . The position variable , disturbance vector and contact position are defined as:
(6) 
And the inputs to the system are defined as:
(7) 
Note that we consider two modes of input torques, constant and linearly increasing with time (ramp profile). Therefore we have eight inputs to the system. Now for single support phase, we write linear differential equations and using Maple software [32], we solve them analytically to find the state evolution after arbitrary time with respect to the initial state.
(8) 
where is a constant matrix and is a linear function of time. Here the vector contains both inputs and states together:
(9) 
Note that inputs and are constant during a single support phase. As expected, could be written as:
(10) 
where and variables are complicated functions of system model parameters, though very sparse. Once these individual matrices are calculated offline, can be easily calculated online by few arithmetic operations. It is also worth mentioning that lateral and sagittal motions are decoupled, meaning that in , only elements with both odd or both even column and row index are nonzero. We also like to emphasize that the variables are assumed to be constant during the whole single support phase and therefore, the corresponding rows on only have on the diagonal. We keep everything packed for simpler manipulation of matrices in future. Note that in , initial foot velocities are zero, starting to swing from rest conditions. It should also be noted that the switch between left and right support only changes the variable and contact positions , not the system transition matrix .
IiB Double support
In this phase, the two feet are fixed and contact forces are being transfered from to . Note that although there are two constraints added regarding the new fixed foot, 6 other forces and torques are nonezero and should be found. Therefore we need 8 additional equations to replace (5) and (7). In our double support, we decide to linearly transfer the weight from one leg to another. This means the vertical component of the Ground Reaction Force (GRF) in previous stance leg (which is constant during single support) will go linearly to zero during double support time. Such policy will ensure smoother torque profiles which are easier to track on the real robot. The linear policy also makes simpler equations in analytic form (compared to quadratic or other forms).
Note that all contact forces are calculable at any time including right at the end of single support (as a function of state variables). If we want to preserve continuity of other horizontal components of GRF in phase transitions as well, in the right hand side of the differential equations (8), terms like will appear. Although these forms are still linear, they make finding analytical solutions difficult. Instead, we keep the Center of Pressure (CoP) for each foot constant during double support. Note that in single support, ground reaction torques in the stance foot are determined by (7) and the vertical GRF is constant. The CoP position can be then simply preserved in double support by linearly decreasing contact reaction moment in stance leg (and increasing accordingly). Overall, equations being considered additionally for the CoP constraints are:
(11)  
(12) 
Linear transfer of weight and transversal contact torques as well as continuity of hip torque profiles require 4 more equations. We solve all equations except 4, and replace variables in these 4 equations. This results in a set of 4 equations with 8 unknown variables, 6 hip torques and 2 vertical components of GRFs. Inspired by linear weight transfer idea, we obtain the following other 4 equations in (13), which transfer torques and forces uniformly during the double support:
(13) 
Here we calculate the Jacobian of these equations with respect to and , multiply them by these variables again and divide each side by proper time to induce linear force transition. Combining the last equation set in (13) with (to have 8 equations for 8 variables) and general equations of the robot mentioned before, similar to single support case, Maple finds analytical solutions of the form:
(14) 
Note that variables and (foot positions) are constant during double support phase and therefore, they are correctly linked to the righthand side vector in this equation. Again as before:
(15) 
Which ensures fast implementation of this matrix. At this stage, we would also like to define row selection matrices that if multiplied from left, they output corresponding rows. If their transpose is multiplied from right, they output corresponding columns alternatively, if dimensions match. These selection matrices are listed in Table.II.
Matrix  selected variables 

all states and except  
sagittal component of  
all inputs  
constant hip torques  
constant contact torques  
timeincreasing contact torques  
IiC Full stride
Once we have matrices for both phases, we can link them together to find the transfer matrix for the full stride phase:
Note that although we have used parameters and in calculating transfer matrices, they are defined by other methods explained later. The variable is crucial for double support calculations, though only determines the rate of timeincreasing input components in single support. The duration of a stride phase is therefore defined as . Apart from the forward transfer matrix , we can also define another back transfer matrix in a similar way:
(16) 
One can easily show that the matrices and could be written as:
(17) 
which favor fast implementation. Now having and matrices, we can define a back transfer matrix as:
Which always satisfies:
(18) 
This equation indicates that any intermediate state could be evolved to the end of the phase simply by multiplying the matrix . For the purpose of visualization or time integration, one can also use individual matrices and with sufficiently small . Figure.3 gives a better understanding of all transition matrices introduced so far as well as the composition of individual single and double support phases to make the full stride phase. Note that for an arbitrary state vector, it is straightforward to calculate all internal forces and torques in closed form as linear and quadratic functions of the state vector respectively.
The matrix is used to find evolution of the system state over time, specially at the end of the full stride where . Although the fact that the two feet are fixed in double support is already encoded in the matrix, one should consider another constraint for the system which forces the foot velocity to zero at the end of the stride. This means two input dimensions should be dedicated all the time and actively being regulated to ensure satisfaction of foot velocity constraint. For this purpose, we dedicate constant hip torque components for example:
(19) 
This equation in fact selects foot velocity rows, calculates the hip torques in terms of other variables and then replaces hip torques in other equations. Now using , foot velocity is always zero at the end of the stride, though two input dimensions are lost as expected. One can do the same to find which is used to find the evolution of current state (measured at time ) until the end of the phase. In next sections we will use these new augmented matrices for simplicity.
So far in this section, we have found transition matrices for the system and defined a full phase, consisting of a double support followed by a single support. We have also formulated matrices such that there are very fast for online calculation. It should be noted that for the purpose of MPC, only a single offline computation of the matrix is enough. In next section, we are going to find different openloop periodic gaits based on the type of actuation desired.
Iii Finding periodic gaits
All we need in this section is the matrix as the system is linear and fully expressed with this matrix. For the purpose of this paper unlike [26], we only focus on symmetric gaits observed in normal human walking. We find different classes of vectors that produce symmetric gaits. The concept of symmetry could be encoded in a single matrix along with the constraint of zero foot velocity at the end of the stride. Consider foot position , pelvis position , pelvis velocity and stance foot position packed in a vector:
(20) 
This vector can be of course extracted from the full vector (9) with the selection matrix . After a stride, we define relative vectors between base, swing foot and stance foot which are linked to those in the beginning in a certain way. These vectors could be defined in the following matrix , if multiplied by (20):
Comparing these quantities before and after a symmetric stride, in sagittal plane, components are equal and in lateral plane, they are linked with a negative sign. This fact could be explained in a matrix :
(21) 
Now we should also consider the transition matrix that exchanges the swing and support contact points after a stride as:
With these matrices, we can define the matrix which enfolds valid periodic gaits in its nullspace:
(22) 
This matrix in fact compares aforementioned vector quantities before and after a stride. Note that apart from state vectors, the hip and contact torques and , the disturbance vector and the variable could also be considered in calculating the nullspace. However, it is meaningless in practice to consider a periodic gait with constant external disturbance. It should also be mentioned that for a valid solution vector, initial swing velocities are zero (impactless model assumption) and contact positions are also set to zero to avoid redundant nullspace dimensions.
Remember that in the previous section, there were two variables to decide: and . In this section, we find various types of gaits by selecting different combinations of actuation and timing. We do so by considering the matrix which is in fact a function of timing variables. Any periodic solution (a vector containing initial states and actuation inputs) should lie in the nullspace of matrix. Note that only columns attributed to nonzero values in the solution vector are selected. We basically exclude columns related to initial foot velocity, contact position and disturbances to obtain the reduced matrix . We then find solution manifolds by combining different actuation dimensions. Indeed, possible actuations are swing hip and stance contact torques (in sagittal and lateral directions) and also their timeincreasing modes. In the rest of this paper, we use humanlike body parameters for numerical simulations, where mass distributions and geometries are taken from [33]. TableIII lists these parameters for two adultsize and kidsize models, used further in this paper.
Model  adultsize  kidsize  unit 

Total mass  70  30  kg 
Body length  1.7  1.0  m 
z1  0.89  0.52  m 
z2  0.32  0.19  m 
z3  0.36  0.22  m 
m1  45.7  19.6  kg 
m2=m3  12.15  5.2  kg 
0.1  0.06  m 
Iiia Pseudopassive gaits manifold
The first choice is to see whether the system has any pseudopassive walking pattern or not. By pseudopassive we mean a gait in which swing hip and stance contact torques are zero. The term pseudo is indeed indicating that in stance hip, the actuators are producing or dissipating power. It also refers to the fact that in our linear model, the legs are stretched or shortened by prismatic actuators, as part of model construction. To find pseudopassive gaits, we look at the reduced matrix (not necessarily square) extracted from by excluding also the 8 columns attributed to hip/ankle torques. Since any valid solution should lie in the nullspace of this matrix, we are interested in inspecting singular values. We normally calculate singular values of , since any vector in the nullspace of this new matrix lies in the nullspace of as well. The resulting singular values are shown in Figure.4 over time.
One can clearly see that there is a time where the system shows two zero singular values. can be found by a simple rootfinding algorithm. It is obvious that one of these singular values refers to the sagittal and the other to lateral dynamics. We can simply calculate corresponding eigenvectors of and find the nullspace manifold, composed of and . Note that there is only one lateral solution as the value of should only be . However the solution in the sagittal plane can be scaled by any arbitrary positive or negative value to obtain different modulated speeds. Therefore the manifold of pseudopassive compass gaits in this case is only 1dimensional. From such analysis, we can also conclude that if any other stride time is chosen, the robot cannot demonstrate pseudopassive forward progression and only steps in place. A demonstration of normal pseudopassive compass gaits can be found in Figure.6.
IiiB Actuated gaits manifold
In this part we are going to find manifolds of motion which can benefit from swing hip actuation and CoP modulation as well. These inputs are of course containing a constant and a time varying components for both sagittal and lateral dynamics, as discussed in the previous section. With this ability, we can pump energy to the swing leg and brake at the end to produce faster swing motions. We can also apply contact torques which modulate the CoP and resemble the fact the CoP in human goes forward from the heel to the toes over a swing phase. We simply perform all calculations on itself. The resulting singular values using the same method described earlier are shown in Figure.5 over time.
It is surprising that the system does not have a distinct zero singular value at like before. However, it has 7 singular values equal to zero that produce a nullspace at any given stride time. Each of the corresponding singular vectors where have similar dimensions with in (9), though with , and to be zero. These are nominal initial states with contact point at origin, resting swing foot and no disturbance of course. This nullspace is not 7dimensional however. The variable which should always be reduces the total dimensions to 6 like before. Now one can simply decide the timing and active actuators in order to reduce this high dimensionality and find a unique solution.
Now for any desired speed, we have the possibility to choose , and a linear combination of resulting 7 nullspace vectors. As demonstrated in pseudopassive gaits, at , a certain combination of these vectors can results in zerotorque gaits. Now what if we calculate the 7dimensional nullspace using ? Could we still find a pseudopassive linear combination? In fact it is possible, although no distinct zero singular value is observed in Figure.5. The reason is that the rank of actuation space at is equal to 5 in the 7dimensional nullspace manifold. This means if we constrain all of them to zero for pseudopassive walking, we only loose 5 ranks. The other 2 ranks are therefore dedicated to the variable and the desired speed, like before. So the nullspace manifold of actuated gaits already encompasses the one for pseudopassive gaits and we do not need to calculate them separately. In the next subsection, we are going to show a few examples of walking gaits using the nullspace calculated.
IiiC Numerical examples
Apart from pseudopassive walking which has a certain timing , we are going to show same speed walking solutions with different timing, once using only hip and once using hip and ankle torques. From singular value analysis, we have 7 singular vectors that can produce motion. We pack them together columnwise in a matrix . We also select a more humanlike choice of and calculated from pseudopassive walking. The choice of speed will be .
We setup an optimization problem to find linear coefficients which produce walking gaits of desired speed with minimal input torques. As mentioned before, this is just a demonstration and does not have any precise meaning in terms of energy. The purpose of this paper is not to find a cost function for energy that produce humanlike motions. It is part of our future work to compare this model with human gaits and their associated timing. In [8, 14] however, limit cycles for the real robot are found through dynamic energy minimization which is interesting and inspiring for future works. By considering a vector of , our optimization is formulated as:
(23) 
Now consider the following scenarios:

Pseudopassive walk: which is being calculated as mentioned before. Note that it can be shown if is chosen, the result of our optimization is the same pseudopassive gait, as obtained before.

Long double support: in this case we enforce ankle torques to zero by adding another constraint to the optimization:
(24) Keeping the same stride time, we double and decrease accordingly. Note that now the walking cannot be pseudopassive anymore, so the optimization will find minimal hip torques to produce the same speed and stride length.

Stage walk: here we constrain ankle torques to zero like before. Now instead of optimizing torques, we optimize the lateral velocity in the cost function. The optimization will give hip torques that produce a motion with minimal lateral bounce. In this case, the biped walks on a straight line without lateral bounce.

CoP modulated: given the length of human foot, the total weight and the timing of single support, we can calculate a constantly increasing ankle torque acting in sagittal plane to move the CoP to the toes gradually during single support. In this scenario, we force other components of ankle torques to zero and the timeincreasing component to by adding the following constraint to the optimization:
(25) The result is a gait with timeincreasing ankle torque profile and proper hip actuation, walking at the same frequency and speed like before.

LIP like: in this case, keeping the original timing, we change the model of the robot. We move most of the weight of each leg to the torso, and also move the 3 masses closer to the pelvis by decreasing and . The idea is to see a behavior similar to LIP. Again we disable all ankle torques as well.
It should be noted that these optimizations are very fast, in the order of microseconds. It is also always possible to find closedform solutions as the optimizations are equalityconstrained quadratic optimizations. The accompanied Multimedia Extension demonstrates different features of 3LP while Multimedia Extension shows movies of five previously mentioned scenarios. The 3D geometry of resulting gaits are shown in Figure.7 while a detailed diagram of each stride is shown in Figure.6. Our flexible model can produce periodic gaits with different actuation schemes. We consider piecewise linear profiles for each actuator to produce more humanlike torque and ground reaction profiles, investigated in the final section.
It can be concluded from Figure.7 that changing different parameters does not have major effect on the overall geometry of walking. However, since our target is to develop a template model that is easy for inverse dynamics to track, we are interested to investigate dynamic properties of these walking scenarios. For this purpose, we have shown CoM velocities for all scenarios in Figure.8.
Although CoM trajectories look similar in Figure.7, they have very different characteristics in terms of velocity variations. The LIP like model shows a large variation in sagittal velocity. It is not so obvious how swing and torso dynamics affect this motion at first glance. Remember that by torso dynamics, we mean the torque required by the stance hip to keep the torso always upright. This torque is not necessarily zero, since the pelvis is an accelerated frame. Therefore, hip torque can affect CoM motion considerably, specially since the torso is relatively heavy. Moreover, although the swing foot has smaller weight compared to other parts of the body, it can be seen from Figure.7 that swing leg has a faster motion during single support. Such fast motion quadratically increases kinetic energy and therefore results in a considerable work flow. In our model, we have described these effects in a simplified and linear fashion, but capturing important couplings between the 3 pendulums.
Taking a closer look at Figure.8 reveals that even maximal CoP modulation still does not change velocity profile considerably. This basically means the difference between pseudopassive walking and LIP is way larger than that between pseudopassive and CoP modulated walking. In other words, CoP authority can at most convert the pseudopassive gait to the CoP modulated gait. The available CoP authority is hardly enough to convert the pseudopassive gait to LIP gait and this gap increases mainly in faster walking speeds. Although here the speed is moderate, we can easily infer that LIP as a template model can only operate in a very limited range of walking speeds. Remember that in fact an inverse dynamics or kinematics approach eventually synthesizes the template motion with the full model by exploiting all control authorities of the robot (including CoP). This motivates therefore not to modulate the CoP in template level and leave the control authority free for fullbody controllers to mimic the template motion as precisely as possible.
In this section we discussed an easy method to find manifolds of periodic motions without any numerical forward simulation of the system. Once these manifolds were found, we also showed how to find individual solutions, based on the type of actuation and timing desired. We only considered gaits with minimal hip torques here. However to go further, we would like to investigate the effect of timing and walking speed as well. Such investigation reveals interesting energetic properties of 3LP, discussed in the next section.
Iv Comparison with human data
Compared to LIP, the 3LP model is much more similar to human locomotion, because it describes falling dynamics, swing motion, torso balance, lateral stepping and double support features all together. In addition to geometric similarities, we plot ground reaction forces and joint torques to compare the underlying dynamics that result in such geometric similarity. For this purpose, regarding the available data from human subjects [34], we select similar model parameters and timing, calculate periodic manifolds and find solutions with the same speed and CoP modulation pattern. Such comparison is demonstrated in Figure.9. In the following, we discuss various similarities observed in this figure.
Iva Sagittal dynamics
From the last column of Figure.9, one can observe a good match of hip extensor and ankle plantarflexor torques as well as AnteriorPosterior ground reaction forces. This is despite the fact that walking speed is relatively fast and the step length is about of the leg length. Note also that the constant and timeincreasing components for hip/ankle torques are roughly enough to describe major trends in human curves. Nonlinear profiles of 3LP are however related to the stance leg and generally those degrees of freedom which are not directly controlled by desirable input torques. The LIP model does not have hip torques and produces larger A/P GRF, because of different CoM trajectories shown in Figure.8.
IvB Vertical GRF
By model design, the CoM height is constant and we do not expect two peaks in the model profiles, similar to [26, 13]. But the general trapezoidal shape is preserved, thanks to our double support phase. The main consequence of such profile is walking with crouched knees which looks less humanlike compared to other template models. Note that LIP model produces the same profile as expected.
IvC Transversal rotation torques
The transversal rotation torques are preserving the general trend of the human data, but not matching very well, specially in the ankle. One major reason is that arm motions and pelvic rotations are not considered in the model. The other important reason is that here we just minimized all hip/ankle torques to find a unique solution out of the manifold of all available solutions. This minimization is not necessarily realistic and humanlike, as it causes wider lateral stepping, larger lateral motion and therefore more ground moment. A better cost function on energy might produce more humanlike gaits, although it is doubted in [35] that optimal gaits are defined merely by energy terms. They might be highly influenced by posture balance, at least in lower speeds. In faster walking speeds like the human data demonstrated here [34], humans tend to take laterally closer steps, compared to our model. Note that the transversal moment is required to keep the torso upright and straight ahead during the swing phase, compensating the moment produced by the swing leg. In LIP however, since there is no pelvis and inertia around the yaw axis, we do not expect transversal torques.
IvD Lateral dynamics
Again due to previously mentioned problem in optimizing lateral motion, we see that our model keeps general trends, like double peaks in the MedioLateral GRF, but cannot precisely describe other torque profiles. Humans normally tend to step as close as possible to minimize lateral motions and energy [10] (remember stage walking in Figure.7). However humans swing their foot over an arc shape to avoid scuffing as well. Such fine motion requires better objective functions to find proper solutions from the available manifolds. Note that the sagittal swing motion can influence lateral dynamics as well [10, 20], possibly through transversal moments. This might be another reason for the discrepancy observed between different lateral curves. Our model however completely decouples lateral and sagittal dynamics. The only linking parameter is stride time which relates to the zerovelocity assumption for the swing foot. In the LIP model, there is no pelvis modeled. However we consider a gait with the same lateral footstepping for better comparison. We can observe that although LIP does not explain hip torques, M/D GRF forces are yet similar to 3LP.
IvE Energetics
Apart from dynamics profiles, it is always interesting to investigate energy flow in the model and compare it with the real human. Remember that one of the main motivations behind developing 3LP is to match humanoid dynamics as precise as possible in the template space. Such matching could be viewed form the powerflow perspective, inspired by the fact that humans can walk very efficiently. Compared to LIP, 3LP can additionally describe swing and torso dynamics which are important aspects of locomotion, specially in faster speeds. Since the pelvis has accelerations in the sagittal and lateral directions, the whole upperbody requires a hip torque to remain upright during locomotion. Our model considers stance hip torques to fulfill this requirement and of course describes the influence of these torques on horizontal motions. We do not model torso movements in 3LP and assume they are negligible, but they might induce additional hip torques too. Although the weight of swing leg is relatively small, its peak velocity is about two times larger than the torso which has a heavier mass. Therefore, the peak kinetic energy of the swing leg is quite comparable with the torso. Such energy comes from both accelerated pelvis where the swing leg is attached to and the swing hip torques that can change swing dynamics. In trade off with the work required for accelerations and decelerations of the torso, this phenomenon explains optimal speedfrequency relations observed in metabolic cost of human walking [36].
Although 3LP does not describe many important walking features like heeltoe motions, knee flexions and CoM excursions, it is still very surprising that it can capture the main optimality trend in human walking. To demonstrate this capability, we consider two main parameters of locomotion, stepping frequency and forward speed. Inspired by [36], we calculate the economy of walking for different speedfrequency combinations. Such Economy is in fact the inverse of cost of transport, which is the total energy consumed per unit mass, per unit distance traveled. There are many ways to calculate such mechanical power in our model. In fact, 3LP does not simulate muscles and exact geometry of a human and is thus, unable to precisely reconstruct the human economy surface based on the metabolic cost. However by integrating the mechanical power on CoM, we are able to approximate a portion of this energy which still plays an important role in the overall energy, according to [37]. We are not able to model other costs like muscle activation, maintenance and shortening heat rates though, since we do not have muscles in 3LP.
By doing a systematic search over a grid of different speedfrequency combinations, we calculate the net positive work performed on the CoM per unit distance, divided by the total mass. The calculation of such energy is rather easy in our model, as it corresponds to the difference between minimum and maximum kinetic energy. Since the velocity of the CoM is a linear function of the state variable, it can be simply related to the known initial state through matrix which is a combination of exponential functions of time. (10,15). The problem reduces to solving two simple maximization and minimization problems which are fast, thanks to the simplicity of closed form solutions. Figure.10.A demonstrates economy contours for the choice of . It is surprising that the model is showing a peak line within the desired range of frequencyspeeds. Such peak in fact demonstrates the tradeoff between hip torques and contact switching acceleration/deceleration costs. The former increases in higher frequencies as energy needs to be pumped into the swing leg and taken out by braking at the end of swing [38]. The latter however increases by steplength which results more variations in the CoM velocity.
We have also repeated the same process for two other double support time choices of and . The resulting peak lines are demonstrated on Figure.10.A as well as the optimal trend of human data, taken from [36] (the case of constrained speed). Here we use the same average weight () and height () of human subjects in [36] as well as average weight distribution reported in [33]. Although 3LP is successful in identifying the tradeoff, using the choice of constant double support time ratio, it fails to match the human data.
It can be postulated that in slow speeds, humans have larger double support time ratios compared to higher speeds (refer to Figure.10.A). Experiments on real humans actually verify this postulation, suggesting a linear relation between double support time and the walking speed [39]. Here we implement similar simple law, making a linear function of walking velocity :
(26) 
This relation gives a ratio of at and at , close to our three initial conjectures. Repeating the same search process with this particular choice of double support time, we obtain Figure.10.C with the actual economy surface shown in Figure.10.D. It is surprising that 3LP is able to predict almost precisely the energyoptimal relation between frequency and speed in human walking. Although 3LP itself cannot probably explain the optimal double support time, it already shows that minimal knowledge of human data is enough. The dependence of (26) on body weight and geometry is yet to be explored in future works. If not dependent, the equation (26) and 3LP seem to be enough to predict optimal frequencyspeed relation, given specific body parameters.
Figure.10 is very promising and important, despite the fact that:

No impact or pushoff is considered.

CoM height is constant.

No foot clearance is modeled.

Heeltoe motions and knee flexion are not included.

Upper body is assumed to be fixed.

The torso has no rotation in any direction.

Weight loading cost on stance leg is not considered.

Walkingindependent metabolic cost is not modeled.
Therefore, although the overall frequencyspeed relation is correctly predicted, the 3D surface of economy is not matching the human data shown in Figure.10.B, taken from [36]. The correct prediction in highfrequency and highspeed walking conditions is surprising. Linear models are conventionally expected to be a linearization of the actual nonlinear system, thus performing well in nearstance (small stridelength) conditions where the geometry of 3LP is more close to the actual system. However here, although 3LP does not model healtoeknee motions and intrinsically simplifies them with an extensible prismatic actuator, it appears to demonstrate the same mechanical energy flow, even with large stride length. In future works, it is still interesting to explore other subcomponents of walking energy, keeping 3LP as a core model. Inclusion of other costs might reconstruct the actual human economy surface more precisely.
Overall, compared to LIP, while being still linear and slightly more complicated, the proposed model is able to describe much more features of human walking. It is very important in hierarchical control architectures that template models produce as accurate motions as possible to make tracking more precise. Despite being linear and of course operating in a more limited region of feasible states, this model is better than LIP to describe human dynamics and consequently, more precise for controlling humanoids. The energy flow of 3LP is also more similar to human than LIP, making the planned motion more natural for the humanoid robot which has similar body features.
V Conclusion
Compared to most of other template models listed in Figure.1, our proposed model considers swing and torso dynamics in a linear formulation. On the other hand, it is computationally similar to LIP which is vastly used in the literature to control real robots over relatively slow walking speeds [23]. Other nonlinear models are also used in simpler robots [20], but over a limited range of speeds. Template models try to describe major dynamics of the robot in an abstract way. This can be used either for analysis of human motion or control of a complex robot, probably in a hierarchy with more complex controllers. In such control paradigm, it is important to keep computational costs as minimal as possible, favoring future prediction.
In 3LP model, the pelvis width links linearly with the lateral motion. This parameter can be used to find lateral ankle/hip torques required for balancing and to determine natural lateral foot placement. This is compared to most of other methods like [2] where the two feet are enforced to be apart to avoid scuffing. In literature, the timing and footstep locations are mainly enforced by other desired trajectories. In our model however since swing dynamics is included, we introduce zero foot velocity assumption and let natural periodic gaits come out of equations. This assumption which relieves the need to calculate impact forces actually determines the timing and periodicity conditions as well.
The proposed model can predict human walking profiles quite well, even for relatively fast speeds where the linearity assumption might be quite limiting for operation on the real robot. Although IPbased models can demonstrate CoM excursions quite well, their nonlinear nature makes them less suitable for highly complex robots that require online planning. There are more advanced versions of IPbased models, including torso and swing dynamics. However again, nonlinear equations cannot be used for per timestep future prediction over a wide range of speeds, although they might better describe energetics of human walking. The proposed model keeps a tradeoff between geometric and dynamic similarities, favoring fast computation properties.
We also showed that 3LP can approximately describe the exchange of energy despite keeping the CoM height constant which is a less realistic assumption in fast speeds. We would like to mention that the vertical excursion of CoM in human locomotion, even in very fast walking at is still about [40] which is quite negligible compared to the leg length of about (pelvis excursion is about however). CoM excursion is mainly determined by the stride length which is not much larger in faster speeds. Humans in fact increase stepping frequency together with step length to walk faster, which highly affects swing dynamics and requires more energy pumped by hip torques. Therefore in faster speeds, apart from the doublepeaked vertical GRF, one should consider swing dynamics as well due to considerable exchange of energy. LIP of course fails to describe swing dynamics, but 3LP successfully captures them.
We have not used it for extensive biomechanical analysis in this paper. Rather, we focus on the fact that such similarity can be inspiring for generating more precise abstract plans, used to control humanoid robots. In brief, advantages of the 3LP can be listed as following:
 +

Swing dynamics.
 +

Torso balancing.
 +

Double support.
 +

Optimal frequencyspeed relation similar to human.
 +

Hip/ankle actuation possibilities.
 +

Computationally fast.
 +

Possibility to consider hip torque limits.
 +

Natural lateral motion.
 +

Natural periodic gaits.
 +

Pseudopassive compass gait.
While disadvantages are:
 

Flat vertical GRF profiles.
 

Stretching legs.
 

No steering capability yet.
 

No arm motion.
 

No torso pitch/roll DoF.
It should be noted that without pelvis width like LIP, steering is possible as demonstrated in our previous work [2]. Steering makes 3LP nonlinear, but one can compromise the lateral motion and let inverse dynamics find proper actuation patterns to turn around the yaw axis. In future works, we are going to replace the LIP with the proposed model and design better controllers to improve the performance on a full humanoid robot. We would also like to setup MPC with all aforementioned policies and inspirations to produce feasible plans for faster walking. 3LP can be extended to have two more degrees of freedom for the torso to describe asymmetries observed in human over fast walking speeds. It can also include actuation inputs proportional to the square of time to produce more precise torque profiles. Among many different advantages of 3LP, we favor its capability to produce more natural motions. In this way, we can expect our inverse dynamics layer to track the template model more precisely and therefore, being able to produce more humanlike motions. This paper is accompanied with a multimedia extension, demonstrating general features of 3LP and the different gaits it can produce. All the codes used in this article as well as the multimedia extension are available online at http://biorob.epfl.ch/page99800en.html.
Vi Acknowledgments
This work was funded by the WALKMAN project (European Community’s 7th Framework Programme: FP7ICT 611832).
References
 [1] P. Zaytsev, S. Hasaneini, and A. Ruina, “Two steps is enough: No need to plan far ahead for walking balance,” in Robotics and Automation (ICRA), 2015 IEEE International Conference on, May 2015, pp. 6295–6300.
 [2] S. Faraji, S. Pouya, and A. Ijspeert, “Robust and agile 3d biped walking with steering capability using a footstep predictive approach,” in Robotics Science and Systems (RSS), 2014.
 [3] S. Feng, X. Xinjilefu, W. Huang, and C. G. Atkeson, “3d walking based on online optimization,” in Humanoid Robots (Humanoids), 2013 13th IEEERAS International Conference on. IEEE, 2013, pp. 21–27.
 [4] S. Kuindersma, F. Permenter, and R. Tedrake, “An efficiently solvable quadratic program for stabilizing dynamic locomotion,” in Robotics and Automation (ICRA), 2014 IEEE International Conference on. IEEE, 2014, pp. 2589–2594.
 [5] P. A. Bhounsule, A. Ruina, and G. Stiesberg, “Discretedecision continuousactuation control: balance of an inverted pendulum and pumping a pendulum swing,” Journal of Dynamic Systems, Measurement, and Control, vol. 137, no. 5, p. 051012, 2015.
 [6] F. Iida, Y. Minekawa, J. Rummel, and A. Seyfarth, “Toward a humanlike biped robot with compliant legs,” Robotics and Autonomous Systems, vol. 57, no. 2, pp. 139–144, 2009.
 [7] H. Hemami and C. Golliday, “The inverted pendulum and biped stability,” Mathematical Biosciences, vol. 34, no. 1, pp. 95–110, 1977.
 [8] F. Asano, M. Yamakita, N. Kamamichi, and Z.W. Luo, “A novel gait generation for biped walking robots based on mechanical energy constraint,” Robotics and Automation, IEEE Transactions on, vol. 20, no. 3, pp. 565–573, 2004.
 [9] K. Byl and R. Tedrake, “Approximate optimal control of the compass gait on rough terrain,” in Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on. IEEE, 2008, pp. 1258–1263.
 [10] A. D. Kuo, “Stabilization of lateral motion in passive dynamic walking,” The International journal of robotics research, vol. 18, no. 9, pp. 917–930, 1999.
 [11] S. Kajita and K. Tan, “Study of dynamic biped locomotion on rugged terrainderivation and application of the linear inverted pendulum mode,” in Robotics and Automation, 1991. Proceedings., 1991 IEEE International Conference on. IEEE, 1991, pp. 1405–1411.
 [12] C. Maufroy, H. M. Maus, and A. Seyfarth, “Simplified control of upright walking by exploring asymmetric gaits induced by leg damping,” in Robotics and Biomimetics (ROBIO), 2011 IEEE International Conference on. IEEE, 2011, pp. 491–496.
 [13] M. Sharbafi and A. Seyfarth, “Fmch: a new model for humanlike postural control in walking,” in Intelligent Robots and Systems, 2015.(IROS 2015). Proceedings. 2015 IEEE/RSJ International Conference on, vol. 1. IEEE, 2015, pp. 5742–5747.
 [14] S. J. Hasaneini, C. Macnab, J. E. Bertram, and H. Leung, “The dynamic optimization approach to locomotion dynamics: humanlike gaits from a minimallyconstrained biped model,” Advanced Robotics, vol. 27, no. 11, pp. 845–859, 2013.
 [15] M. Gomes and A. Ruina, “Walking model with no energy cost,” Physical Review E, vol. 83, no. 3, p. 032901, 2011.
 [16] I. R. Manchester and J. Umenberger, “Realtime planning with primitives for dynamic walking over uneven terrain,” in Robotics and Automation (ICRA), 2014 IEEE International Conference on. IEEE, 2014, pp. 4639–4646.
 [17] R. D. Gregg, A. K. Tilton, S. Candido, T. Bretl, and M. W. Spong, “Control and planning of 3d dynamic walking with asymptotically stable gait primitives,” Robotics, IEEE Transactions on, vol. 28, no. 6, pp. 1415–1423, 2012.
 [18] T. McGeer, “Passive dynamic walking,” the international journal of robotics research, vol. 9, no. 2, pp. 62–82, 1990.
 [19] A. D. Kuo, J. M. Donelan, and A. Ruina, “Energetic consequences of walking like an inverted pendulum: steptostep transitions,” Exercise and sport sciences reviews, vol. 33, no. 2, pp. 88–97, 2005.
 [20] S. Collins, A. Ruina, R. Tedrake, and M. Wisse, “Efficient bipedal robots based on passivedynamic walkers,” Science, vol. 307, no. 5712, pp. 1082–1085, 2005.
 [21] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, and H. Hirukawa, “Biped walking pattern generation by using preview control of zeromoment point,” in Robotics and Automation, 2003. Proceedings. ICRA’03. IEEE International Conference on, vol. 2. IEEE, 2003, pp. 1620–1626.
 [22] P. Sardain and G. Bessonnet, “Forces acting on a biped robot. center of pressurezero moment point,” Systems, Man and Cybernetics, Part A: Systems and Humans, IEEE Transactions on, vol. 34, no. 5, pp. 630–637, 2004.
 [23] Y. Sakagami, R. Watanabe, C. Aoyama, S. Matsunaga, N. Higaki, and K. Fujimura, “The intelligent asimo: System overview and integration,” in Intelligent Robots and Systems, 2002. IEEE/RSJ International Conference on, vol. 3. IEEE, 2002, pp. 2478–2483.
 [24] E. R. Westervelt, J. W. Grizzle, C. Chevallereau, J. H. Choi, and B. Morris, Feedback control of dynamic bipedal robot locomotion. CRC press, 2007, vol. 28.
 [25] R. D. Gregg and M. W. Spong, “Bringing the compassgait bipedal walker to three dimensions,” in Intelligent Robots and Systems, 2009. IROS 2009. IEEE/RSJ International Conference on. IEEE, 2009, pp. 4469–4474.
 [26] J. Rummel, Y. Blum, H. M. Maus, C. Rode, and A. Seyfarth, “Stable and robust walking with compliant legs,” in Robotics and Automation (ICRA), 2010 IEEE International Conference on. IEEE, 2010, pp. 5250–5255.
 [27] M. Kelly and A. Ruina, “Nonlinear robust control for invertedpendulum 2d walking,” in Robotics and Automation (ICRA), 2015 IEEE International Conference on. IEEE, 2015, pp. 4353–4358.
 [28] A. Herdt, N. Perrin, and P.B. Wieber, “Walking without thinking about it,” in Intelligent Robots and Systems (IROS), 2010 IEEE/RSJ International Conference on. IEEE, 2010, pp. 190–195.
 [29] K. Ogata, Discretetime control systems. Prentice Hall Englewood Cliffs, NJ, 1995, vol. 2.
 [30] S. Kajita, F. Kanehiro, K. Kaneko, K. Yokoi, and H. Hirukawa, “The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation,” in Intelligent Robots and Systems, 2001. Proceedings. 2001 IEEE/RSJ International Conference on, vol. 1. IEEE, 2001, pp. 239–246.
 [31] S. R. Hamner, A. Seth, K. M. Steele, and S. L. Delp, “A rolling constraint reproduces ground reaction forces and moments in dynamic simulations of walking, running, and crouch gait,” Journal of biomechanics, vol. 46, no. 10, pp. 1772–1776, 2013.
 [32] M. B. Monagan, K. O. Geddes, K. M. Heal, G. Labahn, S. M. Vorkoetter, J. McCarron, and P. DeMarco, Maple 10 Programming Guide. Waterloo ON, Canada: Maplesoft, 2005.
 [33] P. De Leva, “Adjustments to zatsiorskyseluyanov’s segment inertia parameters,” Journal of biomechanics, vol. 29, no. 9, pp. 1223–1230, 1996.
 [34] J. J. Eng and D. A. Winter, “Kinetic analysis of the lower limbs during walking: what information can be gained from a threedimensional model?” Journal of biomechanics, vol. 28, no. 6, pp. 753–758, 1995.
 [35] J. M. Workman and B. W. Armstrong, “Metabolic cost of walking: equation and model,” Journal of Applied Physiology, vol. 61, no. 4, pp. 1369–1374, 1986.
 [36] J. E. Bertram, “Constrained optimization in human walking: cost minimization and gait plasticity,” Journal of experimental biology, vol. 208, no. 6, pp. 979–991, 2005.
 [37] F. C. Anderson and M. G. Pandy, “Dynamic optimization of human walking,” Journal of biomechanical engineering, vol. 123, no. 5, pp. 381–390, 2001.
 [38] J. Doke, J. M. Donelan, and A. D. Kuo, “Mechanics and energetics of swinging the human leg,” The Journal of Experimental Biology, vol. 208, no. 3, pp. 439–445, 2005.
 [39] G. Cappellini, Y. P. Ivanenko, R. E. Poppele, and F. Lacquaniti, “Motor patterns in human walking and running,” Journal of neurophysiology, vol. 95, no. 6, pp. 3426–3437, 2006.
 [40] S. A. Gard, S. C. Miff, and A. D. Kuo, “Comparison of kinematic and kinetic methods for computing the vertical motion of the body center of mass during walking,” Human movement science, vol. 22, no. 6, pp. 597–610, 2004.