A simulation framework for simultaneous design and control of passivity based walkers

In this paper, we propose a simulation framework which simultaneously computes both the design and the control of bipedal walkers. The problem of computing a design and a control is formulated as a single large-scale parametric optimal control problem on hybrid dynamics with path constraints (e.g. non sliding and non slipping contact constraints). Our framework relies on state-of-the-art numerical optimal control techniques and efficient computation of the multi-body rigid dynamics. It allows to compute both the parametrized model and the control of passive walkers on different scenarios, in only few seconds on a standard computer. The framework is illustrated by several examples which highlight the interest of the approach.


I. INTRODUCTION
Passive walkers are bipedal robots essentially powered by gravity.They exploit their natural dynamics to move forward, but in the meanwhile they are unable to exhibit quasi-static behaviors.Such mechanical systems show an excellent cost of transport (CoT).They are effective platforms, both for biomechanicians to better understand the essence of bipedal walk, and for robot designers to build efficient humanoid robots.
The purpose of this paper is to propose a generic dynamic simulation framework for optimizing both the design and the controller of such robots with respect to a given cost function.An overview of the framework is given in Fig. 1. [1], many passive walkers have been designed and crafted.A great introduction to this field is given in Collins et al. [2].A complete methodology to build incrementally complex passive walkers is described by Wisse et al. [3].It starts from a simple compass-like model and goes to Denise, a 3D dynamic walking robot with a bisecting mechanism for the torso, two arms rigidly coupled to the hip angle, knees that are unlocked during the swing phase, and ankles engineered to steer in the direction that prevents falling.Such incremental approach reveals the critical issues in the design of walker with an anthropomorphic shape.The first issue is the transition from the simplest compass model evolving in the sagittal plane to 3D system [4], [5], [6].A second issue is the extension of the compass model to articulated legs [1], [4], [6], [7], [8].The addition of ankles is addressed in [6], [9], arms in [4], [6] and even neck and head in [10], [11].The mass distribution {saurel,carpentier,mansard,jpl}@laas.fr for a given kinematic structure is also an important question which is addressed in [12], [13].
By opposition to non-passive preview control approaches (e.g.Kajita [18], capture-point [19]) that require footstep planning, the control of passive walkers aims at maintaining at the best the natural steady gait that originates from the mechanical design.Passive walkers do not compute their footsteps in advance.Mechanical design and control are deeply connected and the challenge is to consider both simultaneously.The natural dynamics of passive walkers is periodic.The mechanical design tends to create limit cycles at the origin of the steady gait.The role of the controller is to maintain the system around this limit cycle in spite of environment perturbation.
In model-based optimization, the problem of finding a good controller is expressed as a constrained numerical optimization problem that benefits from solid theoretical analysis of stability [20].This general numerical framework has been applied to passive bipedal walkers for the optimization of mass distribution [12].It has been also instantiated to the study of passive quadrupeds [21] and then to passive walkers for the creation and analysis of efficient gaits [13].In this last paper, the authors propose a dynamic simulation package whose scope is illustrated by various examples including passive walkers and running robots.Our paper falls in the same framework.

B. Problem statement
In this paper, we want to simultaneously solve the two following objectives: 1) for a given kinematic chain architecture, we want to optimize the parameters of a robot (e.g.mass distribution, body segments lengths, slope of the ground, mean forward velocity, etc.), 2) for a given robot, we want to find the best controller with respect to a cost function (e.g. the cost of transport, the minimal time, etc.).In short, for a given kinematic structure, we propose a generic approach to design all robot parameters together with the controller that optimizes a given cost function.

C. Contribution
The first originality of this approach is its ability to deal with complex architectures like human-like walkers, both in 2D and 3D.We do not restrict the motion to only lie in the sagittal plane.Moreover and contrary to similar works on this topic, our framework automatically computes the full dynamics.Therefore, there is no need of writing down the complex dynamic equations of polyarticulated systems.Hence, many passive walkers can be efficiently designed, optimized and compared.
Secondly, control can be chosen to be active or passive.We also handle periodic as well as non-periodic gaits.
Finally, we can optimize various parameters of a given walker (slope, lengths, masses, speeds, etc.) with respect to a given cost function.Fig. 1 illustrates the global architecture of our framework.

D. Outline of the paper
In Sec.II, we first introduce the dynamic contact simulator at the root of the framework.We highlight how the various problem parameters are distributed either in the mechanical model or in the controller.In Sec.III, we set up the generic optimal control formulation which allows the computation of both design and control parameters.Finally, the setup is introduced in Sec.IV and experimental results are presented in Sec.V.

II. DYNAMIC CONTACT SIMULATION
Passive walkers are intrinsically hybrid systems.They are submitted to a continuous dynamics when the stance leg is in contact with the ground, and they are also subject to impacts when the swing foot hits the ground.In addition to this hybrid dynamics, some contact constraints must be satisfied to ensure the feasibility of the entire motion.
In the following of this section, we expose the notations, the parametric model of the walkers and the contact formulation used inside the framework.

A. Notations
In this paper, we assimilate a passive walker to a free-floating base system.We denote by q ∈ SE(3) × R n its configuration vector, with SE(3) the special Euclidian group of dimension 3 encoding the placement of the robot base and n the number of degrees of freedom (DoF).The tangent velocity and acceleration of the configuration vector are denoted by q and q respectively and live in R 6+n .Finally, the torque applied at each joint is denoted by τ ∈ R n .

B. The parametric model
A passive walker is primarily a kinematic tree, i.e. a tree of joints where each joint has its own topology (e.g.revolute, free-flyer, spherical, etc.) and a particular placement regarding to its parent joint.The joints are the nodes of the tree.In addition, each joint supports a body, which is defined by its mass, its center of mass (CoM) and its inertia matrix.All the bodies together define the mass distribution of the model.The tree structure with the mass distribution correspond to the structural parameters of the system.The model of the passive walker is then parametrized by those two sets of parameters: model(tree, mass_distribution)

C. The parametric controller
A gait is characterized by its controller which is represented by a set of real parameters.For instance, a controller may be a set of splines which encode the torque trajectories or just the PID gains (stiffness and damping of the spring attached to the joint) in the case of pure passive controller.

D. Continuous contact dynamics
For the continuous dynamics, we make the hypothesis of punctual rigid contact with Coulomb friction cone.The dynamic equation of the polyarticulated system with constraints can be stated as: J c (q)q + Jc (q, q) q = 0 where M (q) is the joint space inertia matrix, b(q, q) corresponds to the Coriolis, centrifugal and gravitational effects, S is a selection matrix encoding the underactutation, J c (q) is the contact Jacobian with f c the contact forces and .denotes the transpose operator.
A necessary and sufficient condition for non-sliding and non-slipping of the contact is described by the constraint that f c must remain inside the Coulomb friction cone K c .This cone reflects the fact that the normal component of the contact force is positive (the ground cannot pull), the norm of its tangential components and the normal torque are limited by the normal component.
We make the choice to solve (1)-( 2) together, and to add the conic constraint inside the optimal control problem to enforce the Coulomb contact model [22].This implies that we have to impose contact phases where (1)-( 2) are enforced.Other approaches have tried to get rid of this extra hypothesis [23], [22] leading to difficult and still incompletely solved problems for trajectory optimisation [24], [25].In the context of this paper, selecting in advance the contact phases is not a limitation.
The joint acceleration and contact force vectors are then given by: with Λ c def = J c M −1 J c the so-called operational-space inertia matrix and I n the identity matrix of dimension n.The dependences on q and q have been omitted for simplicity of notation.

E. Impact dynamics
Passive walkers are also subject to impacts.Here, we make the assumption of instantaneous inelastic impacts with a restitution coefficient sets to zero, i.e. the post-impact velocity of the contact point is null.The impact dynamics then leads to a discontinuity in the joint velocity space, which is described by the two following equations: with q− , q+ the pre-impact and post-impact generalized velocities and λ c is the impulse resulting from the impact [22].Other impact model (e.g.elastic) could also be introduced without loss of generality.This impact model is frequently used in the litterature [26], even though it is contested for its physical consistency [27].

F. Efficient rigid-body dynamic computation
The optimal control solver must evaluate thousands of times the multi-body dynamics either inside the numerical integration procedure or for the evaluation of the dynamic sensitivities regarding to the model and control parameters.For that purpose, we used Pinocchio [28], a whole new, open-source and efficient C++ library to model and compute the forward and the inverse dynamics of polyarticulated systems in contact.Pinocchio is written on top of the efficient Eigen C++ library [29] for Linear Algebra.It is based on Featherstone's algorithms [30] but they have been implemented in a way to take profit of branch prediction and cache mechanisms of modern processors [31].

CONTROL
Optimal control is a powerful and generic mathematical tool which allows the exhibition of a particular solution among an infinite number of candidates.While geometric results only exist for a very limited class of systems, the past few years have seen the emergence of efficient and reliable numerical optimal control frameworks working on high dimensional and complex systems [32], [33].
In our framework, the simultaneous search of model and control parameters is set up as an optimal control problem with a prescribed cost function.This cost function reflects the objective of the gait and can be any real value function.Some examples of cost functions in the context of passive walkers are the cost of transport or even the minimal time.In addition, we add the possibility to set the duration of the motion as a free parameter of the problem.Other free variables, like the gait stride length for example or the slope of the ground, can be stacked to the list of parameters.

A. Notations
We denote by x def = (q, q) the state of the systems, u is the control vector and p is the vector of parameters, composed of both the model parameters and the aforementioned free variables.State and control trajectories are denoted by x and u respectively.The cost function and dynamics of the system are then written as (t, x, u, p) and dx dt = f (t, x, u, p) respectively.Here, we use a slight abuse of notations to denote with the same notation both the continuous and the impact dynamics.Finally, g(t, x, u, p) corresponds to the inequality constraints that must be satisfied along the path (equality constraint can be also considered without loss of generality).As mentioned in Sec.II-D, the function g is first and foremost composed of the Coulomb conic constraints.

B. The Optimal Control Problem formulation
The hybrid dynamics of passive walkers can be seen as a multi-phase system, each stage corresponding either to the single, double support or impact phases.Thereafter, the integer s refers to the index of the s th stage.
The generic optimal control problem for simultaneously computing model and control parameters with multi-phase dynamics can be written as: where π in (7d) is a function which acts both on the state and control trajectories to enforce periodicity constraints.∆t s is the duration of the phase s, and T = ∆t s is the total duration of the motion.In case of impacts, the phase duration is reduced to 0.

C. Solving the Optimal Control Problem
Two major directions exist to solve the infinite dimensional problem (7).The first direction belongs to the so-called indirect methods.It consists in exploiting the necessary conditions for optimality, namely the Pontryagin's maximum principle [34], which transforms the problem (7) into a boundary value problem working on ordinary differential equations.However, such methods are currently unable to track path constraints.
The second direction corresponds to direct methods.Direct methods first discretize the original problem into a finite dimensional nonlinear programming problem (NLP), which is then solved with standard NLP strategies.Among NLP strategies, three of them are now popular: (i) single shooting, (ii) collocation and (iii) multiple shooting.In what follows, we briefly survey these three methods.For further details, we refer the reader to the general overview on numerical optimal control methods written by Diehl et al. [35].
1) Single shooting: discretizes the control and constraints according to a temporal grid.The state trajectory is recovered by integration of the discrete control trajectory along this grid.As single shooting method reduces the NLP to the search of a control trajectory, the optimization problem is of low dimension.However, the solver is hard to initialize if only an initial guess on the state trajectory is available or it may not converge at all in the context of unstable systems.
2) Collocation: disctrizes both the control and the state trajectories according to a temporal grid.In addition to the classic discritized constraints, the state trajectory is enforce to match the dynamics equation (7b) at each grid node.The problem can then be easily initialized from a given state trajectory and collocation handles well unstable dynamics.However, a very fine grid is required to make the state trajectory closer to the true dynamics of the system.
3) Multiple shooting: takes profit of both previous methods.It works on a coarser time grid, which are called multiple shooting intervals.On each interval, the control is discretized as well as the initial state value.The final state value on the interval is then obtained by integration of the system dynamics (7b).In this way, each interval is set independent from its neighbours.And the dependencies between successive intervals is shifted as equality constraints of the NLP.The NLP remains a low dimensional problem and can be easily warm-started with an initial guess on the state trajectory.Furthermore, multiple shooting is really suited for multi-phase dynamics, as each phase is set independent to the others.
Multiple shooting has been successfully applied for the modelling of human running [26].Following several advantages listed in Sec.III-C3, we chose this strategy to solve the problem (7).

D. Efficient optimal control solver
Our framework relies on the 20 years old optimization package MUSCOD-II [32], a multiple shooting solver for highly nonlinear systems submitted to path equality and inequality constraints, developed inside the Optimization and Simulation group at the University of Heidelberg.MUSCOD-II handles multi-phase systems with discontinue dynamics and periodic constraints with efficiency.While MUSCOD-II is a closed-source framework, one can depend on ACADO [33] which implements similar features.

IV. EXPERIMENTAL SETUP
The generality of our simulation framework is illustrated by various examples depicted in Tab.I.Those examples include different kinematic structures and different control schemes.Doing so, we may compare different polyarticulated topologies and, for a same topology, different control schemes, including either active or passive actuators.
Once a topology is chosen, we show that it is possible to replace the active actuation by a passive spring damper system.The motions of the first five robots are constrained to lie in the sagittal plane.Such restriction is removed for the robot with arms which is in 3D.

A. Inputs and Outputs
Fig. 1 shows the general inputs and outputs of our framework.In the experiments of Tab.I, the inputs of the system are the kinematic structure of the robot with the anthropometric parameters of the body segments (see Sec. IV-E), and the cost function, constraints and actuation type (see resp.Sections IV-B to IV-D).
In those experiments, we choose to set at a fixed value both the step duration (0.8 s) and the slope (0.05 rad) for each scenario, and let the solver optimize the step length.We limit the step length to the bounds [0.4; 1.] m.
The outputs of the system are the optimal cost of transport, together with the associated step length of the gait and the state and control trajectories of the joints during one step.
We also collect the number of iterations needed by the solver to converge and the total computation time for each experiment.

B. Cost function
We use the same cost function in each scenario.This cost function corresponds to the classic cost of transport (CoT).The CoT is a non-dimensional quantity which reflects the energy efficiency of a locomotion pattern.By definition, the CoT is the ratio between the energy consumed by the system (E) and its weight (mg) multiplied by the traveled distance (d): with g the gravity field value and m the mass of the system.The energy (E) consumed by the system is the sum of the potential energy mgh (with h the variation in altitude of the center of mass) with the integral of the power input In this last formula, δ is the Dirac impulsion corresponding to the impact instant, • is the dot product operator and x M def = √ x M x.Finally, on a slope with angle α, the CoT is given by: This cost function is only used as an example for the different case studies in Sec.V; in real life examples, the choice of a better suited cost function is still an open problem.Moreover, the CoT we give are not meant to be compared with examples outside of this paper.

C. Constraints
To reduce the dimension of the NLP, we compute a solution only for a half-step due to the periodicity and the symmetry between right and left segments.We then constrain the swing leg to be in contact with the ground only at the beginning and the end of the simulation.Then, the cyclicality  Cases studies are based on six walkers models to measure the impact of the knees, torso, neck, arms and actuation type.For each example, the considered output are the cost of transport and the step length.The two last lines give the performance of the algorithm.In all examples the duration of a step is fixed (0.8s) as well as the slope of the ground (0.05rad).
and the symmetry of the motion is enforced by periodic constraints on both the initial and final configuration, velocity and torque of each joint trajectory.

D. Actuation
In a first stage, we propose an active actuation.The torque trajectory of each joint is modelized by piecewise cubic polynomials.An hardware implementation would then need an external source of power for the system, as a battery or an air canister, and an intelligent controller to deliver the desired torque.
In a second stage, we compare this active actuation with a passive one, which is simply a proportional derivative (PD) controller.In this case, we apply the torque τ j to the joint j: τ j = −K pj (q j − q 0j ) − K dj qj (10) where q j is the position of the joint and qj its velocity.K pj is the spring constant, q 0j the free length of the spring, and K dj the damping coefficient.Those three parameters are optimized by the solver.From those spring-damper parameters, one may craft a purely passive walker.In this case, the single source of power is the gravity field.Of course, another constraint for the numerical solver is to use the same spring and damper for each symmetrical joint.

E. Body parameters
All the robots have an anthropometric mass distribution, i.e. segment length, mass, center of mass position and inertia tensor are scaled according to a reference height and a reference mass.Those inertial parameters follow the anthropometric table proposed in [36].
Our minimal model is composed of a pelvis, two thighs and two legs.Ankles (and knees on the model M B ) are actuated.On top of that, in the models M C , M D and M E , we add a torso and a head.This head is only actuated in the model M D .At last, we add two arms and forearms in the model M E with actuated shoulders.All actuated joints are revolute along parallel axes.

V. EXPERIMENTAL RESULTS
In the following, all the experiments are composed of 9 nodes for the single support phase and 1 single node for the impact phase.

A. Influence of the knees
The influence of the knees is highlighted by comparing models M A and M B in Tab.I.
The experimental results show that adding knees to a compass walker roughly divides by two the CoT while increases the optimal step length.The computational cost of adding knees is negligible.

B. Influence of the torso
Here, we study the impact of adding a fixed torso-head segment attached to the pelvis.This situation corresponds to the models M A and M C of Tab.I respectively, with a non-passive actuation.
Then, adding a fixed torso on top of the pelvis leads to a 40 percent decrease of the CoT and also reduces the optimal step length down to 40 cm, which is the lower bound we allowed.In the meanwhile, the computation times remain the same.

C. Influence of the neck
We consider the influence of actuating the neck by comparing the models M C and M D while keeping the actuation type.
The experiments on those models show that the increase of the CoT is lower than one percent even if we add an actuator.However, the solver needs more time to converge.

D. Comparison of actuation type
This experiment considers the model M C only and study the differences between active and passive controllers.
If we confront the results of the third and fourth columns of Tab.I, it appears that the cost is higher for the passive walker with respect to the choosen cost function, while the computational time is lower in the passive case.This can be explained by the dimensionality of the problems: in the first case, the dimensionality of the joint actuation is defined by four polynomial parameters times the number of shooting nodes while it is defined by three scalar parameters corresponding to the spring-damper model in the other case.

E. Beyond the 2D sagittal plane
Our simulation framework is general and it allows to address 3D model.In order to let the robot keep its balance, we add two actuated arms on the model M E .
In comparison to the 2D model without arms M C , the CoT is a bit higher, but remains lower than the first compass walker M A .We also notice that both the computation time and the number of iterations are higher, but the solver still converges in a similar time.

VI. DISCUSSIONS AND FUTURE WORKS
While our framework shares some common features with the Matlab package developed by Remy et al. [13], its differs on various aspects.The first difference lies at the optimal control method level.In [13], the authors use either collocation or shooting methods.Still, as mentioned in [35] and recalls in Sec.III-C, the multiple shooting approach is the most suited in the context of multi-phase dynamical systems submitted to impacts and constraints.In addition, the choice of the C++ programming language to write both our dynamic simulation and optimal control problem is fundamental for efficient and fast computation.
Currently, we do not check the stability of the resulting motions, but only its balance at a dynamic level.As future extension, we plan to implement an approach similar to [37] to converge on optimal and stable open-loop motions.
An other extension of this framework can be achieved at the level of contact modelling.Currently, punctual contacts are the standard for passive walkers.They are easy to simulate and to mechanically integrate as feet for passive walkers.However, rolling contacts may lead to more efficient gaits, as suggested by Kuo et al. [38].We recently proposed in [39] an analytical formulation of the rolling contact which is compatible with our dynamic formulation presented in Sec.II-D.The addition of this contact model inside our framework can largely improve the quality of the generated motions and enable us to deal with more complex scenarios.
Our framework is not restricted to passive walkers.On a wider scale, it can evaluate some paradigms in the design and the control of new humanoid robots [40], [41] and exoskeletons [42].In the near future, we will exploit this framework to evaluate and design a bipedal robot inspired from passive walkers.

Fig. 1 :
Fig.1: Overview of the simulation framework.The simulator is described in Sec.II, and the solver is detailed in Sec.III.