vector of patient's support leg joint torques in R 6 [N.m] τ sw P D vector of exoskeleton swing leg high-gain PD joint torques in R 6 [N.m] τ sw GC vector of exoskeleton swing leg gravity compensation torques mapping joint position in R 6 to joint torques in R 6 [N.m] 1.1. EXOSKELETONS FOR WALKING REHABILITATION 17 To this day, REX (Rex Bionics) and Atalante [Gurriet2018] are the only selfbalanced exoskeletons on the market, while the XoMotion (Human in motion Robotics) [Shahrokhshahi2022] is still under development

Abstract

Self-balanced exoskeletons offer individuals with disabilities recovered ambulation capabilities thanks to their embedded control algorithms which realize autonomous and stable walking. This thesis is concerned with enabling patients to use these systems for rehabilitation purposes according to the concept of active gait rehabilitation. In our interpretation, this concept entails giving back some freedom of motion to the patients. This objective is conflicting with ensuring safety, which implies taking control over the motion. The main topic addressed in the thesis is the challenge of reconciling these two objectives.

Our solution is twofold: empower the patient and monitor safety. Firstly, we introduce a split in the control architecture, allowing patients to actively participate in the swing leg motion. Concurrently, the control algorithm of the support leg regulates the unilateral contact forces. Secondly, we solve online optimal control problems in free final time on a reduced model to generate stabilizable trajectories. The duration of these trajectories serves as a determinant for when the control system should take back control over the execution of the walking pattern.

The manuscript presents the necessary control architecture updates along with the development of the safety management algorithm. Experimental results are provided to assess the effectiveness of our approach in both flat-foot and foot-rolling walking exercises.

Résumé

Les exosquelettes auto-équilibrés, grâce à leurs algorithmes de contrôle qui permettent une marche autonome et stable, offrent aux personnes handicapées des capacités de marche retrouvées. Cette thèse vise à permettre aux patients d'utiliser ces systèmes à des fins de rééducation selon le concept de rééducation active de la marche. Selon notre interprétation, ce concept implique de redonner une certaine liberté de mouvement aux patients. Cet objectif est en conflit avec la garantie de la sécurité, qui implique de prendre le contrôle du mouvement. Le principal sujet abordé dans cette thèse est le défi que représente la conciliation de ces deux objectifs.

Notre solution est double : donner de la liberté au patient et en garantissant l'équilibre. Tout d'abord, nous introduisons une division dans l'architecture de contrôle existante, permettant aux patients de participer activement au mouvement de la jambe de vol. Parallèlement, l'algorithme de contrôle de la jambe de support régule les forces de contact. Deuxièmement, nous résolvons de manière embarquée des problèmes de contrôle optimal en temps final libre sur un modèle réduit afin de générer des trajectoires stabilisables. La durée de ces trajectoires sert à déterminer le moment où le système doit reprendre le contrôle de l'exécution de la marche.

Le manuscrit présente les mises à jour nécessaires de l'architecture de contrôle ainsi que le développement de l'algorithme de gestion de la sécurité. Des résultats expérimentaux sont fournis pour évaluer l'efficacité de notre approche lors d'exercices de rééducation de la marche en utilisant une trajectoire sans et une trajectoire avec déroulé du pied.

Je tiens à remercier tous les membres du jury d'avoir pris le temps de relire ce manuscrit et d'assister à ma soutenance. En particulier, merci à Katja Mombaur, pour tous vos commentaires constructifs sur mon manuscrit, merci à Aaron Young, pour toutes vos questions et d'avoir fait le voyage spécialement pour ma soutenance, à Heike Vallery, pour nos échanges et votre aide à ICRA, ainsi que d'avoir présidée ce jury, et merci à Pierre-Brice Wieber pour la clarté de vos travaux qui m'ont grandement éclairci les idées, ainsi que vos questions très intéressantes lors de ma soutenance.

Merci également à toi Marine pour ton soutient constant, ta patience et l'accueil challengeant mais toujours positif que tu réserves à mes idées toujours surprenantes voir extravagantes depuis plus de quatre années. Merci Florent pour la patience que tu fais preuves depuis mes premiers TD d'automatique avec toi et mes mille questions évidentes posées sans vergogne, mais aussi pour tous nos échanges informels, intéressants et tout autant instructifs. Merci Nicolas pour la confiance que vous m'accordez depuis nos premiers échanges lorsque je cherchais un stage de césure, pour votre disponibilité et votre soutien, tant technique que moral, ainsi que tous nos café-débats sans filtre enrichissants.

Je tiens évidemment à remercier mes parents, pour tout ce qu'ils ont fait pour moi, mais en particulier pour m'avoir toujours permis de travailler dans de très bonnes conditions et de m'avoir toujours poussé à donner le meilleur de moi-même. Je remercie également Fanny, pour m'avoir soutenu depuis le début de cette thèse, autant dans les bons moments que dans ceux plus difficiles, pour toujours avoir été de mon côté quand je me plaignais, peu importe ce que tu en pensais vraiment, pour ta patience quand tu écoutais mes anecdotes et élucubrations pas si (peu) intéressantes, et enfin d'avoir supporté quand je te réveillais avec mes insomnies. Je remercie également tous mes amis d'avoir supporté mon indisponibilité chronique et d'avoir pourtant été là pour me changer les idées à chaque fois que j'en avais besoin.

Enfin, je tiens à remercier tous mes collègues, tant ceux du CAS que 6 ceux de Wandercraft, pour toutes nos discussions enrichissantes et votre aide, avec une mention spéciale pour mes collègues de Wandercraft qui supportent avec patience d'écouter mes idées farfelues sans s'en plaindre, même si je ne choisissais pas toujours le meilleur moment pour vous embêter. to the CoM state in R 4 [m, m.s -1 ] q vector of generalized positions in R 3 × R 15 [m x rad] q sw vector of swing leg joint positions in R 6 [rad] q sp vector of support leg joint positions in R 6 [rad] q un vector of unactuated DoF in R 3 × R 3 [m x rad] q t vector of target joint positions in R 12 [rad] q sw t vector of swing leg target joint positions in R 6 [rad] q sp t vector of support leg target joint positions in R 6 [rad] q T R time-rescaled joint trajectory, computed from the nominal joint trajectory T , mapping time in R + to joint space in R 12 [rad] q vector of generalized velocities in R 3 × R 15 [m x rad] -1 qsw vector of swing leg joint velocities in R 6 [rad.s -1 ] qsp vector of support leg joint velocities in R 6 [rad.s -1 ] qun vector of generalized velocities in R 3 × R 15 [m x rad] -1 qt vector of target joint velocities in R 12 [rad.s -1 ] qsw t vector of swing leg target joint velocities in R 6 [rad.s -1 ] Glossary 15 qsp t vector of support leg target joint velocities in R 6 [rad.s -1 ] q vector of generalized accelerations in R 3 × R 15 [m x rad] -2 qsw vector of swing leg joint accelerations in R 6 [rad.s -2 ] qsp vector of support leg joint accelerations in R 6 [rad.s -2 ] qun vector of generalized accelerations in R 3 × R 15 [m x rad] -2 qt vector of target joint accelerations in R 12 [rad.s -2 ] qsp t vector of support leg target joint accelerations in R 6 [rad.s -2 ] τ vector of exoskeleton joint torques in R 1 2 [N.m] τ sw vector of exoskeleton swing leg joint torques in R 6 [N.m] τ sp vector of exoskeleton support leg joint torques in R 6 [N.m] τ u vector of patient's joint torques in R 12 [N.m] τ sw u vector of patient's swing leg joint torques in R 6 [N.m] τ sp Chapter 1 Introduction Chapitre 1-Introduction: Ce chapitre introductif décrit les exercices de rééducation de la marche, et en particulier s'attarde sur les principales motivations des développements robotiques pour la rééducation. Les différentes familles d'exosquelettes, sous-catégorie des robots de rééducation, sont présentées, ainsi que les différents types d'exercices de rééducation qu'ils permettent. Une caractéristique essentielle de ces exercices est mise en avant : la possibilité pour le patient de modifier le mouvement à l'aide de sa force musculaire. Nous appelons exercices actifs les exercices permettant une telle modification du mouvement par le patient, et, par opposition, exercices passifs ceux ne permettant pas de modification du mouvement. L'exosquelette auto-équilibré Atalante, cible des travaux présentés dans ce manuscrit, et son mode nominal de fonctionnement, permettant la rééducation passive de la marche, sont décrits. Ensuite, deux catégories de changement permis par les modes de controôle actifs d'exosquelettes sont mis en avant: les changements de la forme des mouvements et les changements de la durée des mouvements. Finalement, la principale contribution de ce manuscrit, un contrôleur offrant une liberté temporelle au patient tout en monitorant l'équilibre, est présentée, ainsi que l'organisation du reste du manuscrit. 1.1 Exoskeletons for walking rehabilitation

Acronyms

A concise view of walking rehabilitation

Human walking can be seen as repetitively achieving a gait 1 . Each leg has a stance phase and a swing phase. These define, for the human, single and double support phases. They are illustrated on Figure 1.1. At some instant of the double support phase, the weight is transferred from one leg to the other and the Center of Pressure (CoP) 2 shifts from one foot to the other. At all instants of single support phases, the human body is propelled forward. Each single support phase includes at least one controlled falling phase [Perry2010],

i.e. a phase during which the Center of Mass (CoM) of the human body is not directly above the base of support, the support foot, and the landing of the other foot is required to prevent an actual fall 3 . Walking requires two primary abilities: transferring the weight during double support phases and moving both the CoM and the swing foot forward during single support phases.

1 A gait is a cyclic motion pattern that produces locomotion through a sequence of foot contacts with the ground [Haynes2006]. 2 The Center of Pressure is the average of all the contact points weighted by the vertical component of the force pushing on the ground.

3 This is a particularity of human walking. Impaired patients following gait training programs are trained to regain the ability to walk (see Figure 1.2 for examples). Because of their neurological disorders 4 , they have strong difficulties controlling the weight transfer or producing the efforts required to generate healthy single support motion by themselves. Walking rehabilitation is exhausting for both the physiotherapist and the patient because they need to maintain the balance of the patient during the exercise by producing important external forces, in a well-coordinated manner. These efforts often even lead to lower exercise intensity and shorter durations of training sessions than what would be ideal. This is an important matter as strong evidence indicates that higher intensity (in the duration of exercises, number of repetitions and level of effort) leads to faster recovery [Lohse2014;Lang2015;Hornby2019;Moore2020].

With their powerful actuators, robotic devices have the potential to handle safety [Gassert2018], helping both the physiotherapist in their work and the patient in their recovery process. Robotic rehabilitation is especially promising for gait training as the potential load that could be lifted from physiotherapists and patients is heavier. For this reason, exoskeletons are becoming more and more popular for gait training [Doppmann2015;Young2017;Eveld2021;Shahrokhshahi2022;Kang2023].

Using exoskeletons during training

Exoskeletons are a special class of robots [Gassert2018]. They are meant to be worn by humans s.t. they can exert efforts on several body parts simultaneously. They vastly differ from robotic manipulators which interact with humans only at their end-effector (see [Maciejasz2014] and references therein). Exoskeletons are designed to align some of their joints with the morphological joints of the human wearer, s.t. they share some Degrees of Freedom (DoF), as opposed to manipulators (which rarely align their joints with humans').

Lower-limb exoskeletons (see examples on Figure 1.3) are a subclass

EXOSKELETONS FOR WALKING REHABILITATION

21

specifically designed for walking rehabilitation or ambulation5 . They are the subject of this thesis. From now on, we only discuss lower-limb exoskeletons, and we simply call them "exoskeletons".

Passive vs active gait training

With such systems, two modes of operation must be distinguished. Passive gait training is reported to be less effective than active gait training, as it leads to the phenomena of "learned helplessness" [Skinner1979;Wool1980]. Intuitively, one can easily understand that rewarding the proper efforts and penalizing the inappropriate ones should encourage the patients toward producing the right efforts. If the patient has no freedom of motion and senses no difference in the robot's motion regardless of their efforts, they will tend to let the robot produce all the efforts. This is a negative and undesired outcome.

Three families of exoskeletons

During gait training the exoskeleton must guarantee the safety of the patient. Chronologically, three types of solutions have been developed for this purpose. They have driven the development of three different types of exoskeletons: platform-based exoskeletons (for instance Lokomat pictured in Figure 1.3a), which are virtually fixed to the ground (either by being truly fixed to the ground or by being sufficiently heavy s.t. their base remains still), crutched exoskeletons (for instance Ekso pictured in Figure 1.3b), designed to be used with crutches or external help for stabilization, and, more recently, self-balanced exoskeletons (for instance REX pictured in Figure 1 These various types of exoskeletons constitute various trade-offs between the realism6 , safety7 and intensity8 of the walking exercises. Platform-based exoskeletons are, by definition, the most stable as they ensure the patients will never fall. They also allow high-intensity walking exercises. However, they reduce the realism of the walking exercise. Hence, they favor safety and intensity over the realism of the walking exercise.

The realism of the walking exercises is higher with crutched exoskeletons, as the legs of the patients have to support their weight. However, they leave the patients and/or physiotherapists with the exhausting walking stabilization task, which limits the intensity of the rehabilitation exercises. The probability of falling with these exoskeletons is not small. Hence, they favor realism over the safety and intensity of the walking exercise.

Self-balanced exoskeletons, while not suffering from realism or intensity issues, face a major challenge: ensuring the balance of the walk. Designing a control algorithm to address this issue is not straightforward. In principle, self-balanced exoskeletons could re-employ the vast humanoid and bipedal robots' literature to realize stable walking. However, the presence of a patient in the system implies specific needs, unanswered by the literature9 . First, some freedom of motion must be given to the patients to boost their rehabilitation, while humanoid robot controllers favor accurate tracking of the planned motion. Second, by producing inappropriate efforts, and because the behavior of the patient is only approximately known and difficult to model, patients may disturb the balance of the system during walking exercises, which implies the design of stabilization algorithms especially robust to disturbances and model errors. We are interested in the design of gait training control methods for self-balanced exoskeletons.

The work presented in this manuscript targets one of these self-balanced exoskeletons, Atalante. We now describe its hardware and the control algorithms designed for passive walking.

Description of the Atalante platform

The exoskeleton Atalante is pictured on Figure 1.4. It has been certified as a medical device for passive use by paraplegic patients in the European market (in 2019) and in the US market (in 2023). Atalante carries a patient standing on its feet s.t. the robot supports its weight. In addition, the patient is strapped to the exoskeleton at his feet, shins, thighs and trunk (see Figure 1.3d). To ease its use in medical care centers it is possible to adjust the length of the shins and thighs to match the patient's body parts length and align their joints with the exoskeleton. In total, the exoskeleton has 13 segments, and 7 of them are attached to the body parts of the patient. The 13 segments constitute a kinematic chain having 12 actuated revolute joints (6 for each leg: two joints at the ankle, one joint at the knee, and three joints at the hip). These joints are mechanically optimized to be best aligned with the morphological joints of the patients.

PASSIVE WALKING WITH A SELF-BALANCED EXOSKELETON 25

This mechanical structure is put into motion by 12 independent electrical actuators, each composed of an electrical motor and a gearbox. The whole exoskeleton is powered by two batteries, located on its thighs, which offer an autonomy of two hours under intensive use. Overall, the exoskeleton weighs approximately 80 kg. Except for the ankles, all the gearboxes feature a high reduction factor (approximately 100). The position of each joint of the exoskeleton is measured by encoders located on the motor side. Atalante also possesses 8 single-axis force sensors (4 under each of its feet) used to sense the contact forces. In addition, 5 Inertial Measurement Unit (IMU) are rigidly mounted on its feet, shins and back of the exoskeleton. They serve two purposes. First, they are used to estimate the underactuated position and orientation of the back of the robot, through the sensing of the proper acceleration, which includes the gravity vector defining the local frame, the direct kinematic model of the exoskeleton and the joint angles. Then, they are used for the estimation and control of the mechanical structure's deformations. Indeed, the Atalante mechanical structure undergoes deformations, similar to those of rotational springs located in its hips and ankles [Vigne2021] 10 .

Atalante's control algorithms are running on the onboard dual-core Intel Core i5-4300U @ 1.90GHz computer. The robot trajectories, tailored to each patient's measurements (their weight, height and length of the different parts of their legs), are generated offline for each patient. Various motion patterns are available: walking in a straight or curved line, walking sideways and backward, turning around, standing up, and sitting down. The modes of control and motions of the exoskeleton are commanded by a two-step procedure: first, the patient uses a remote to select the desired control mode, then the motion is triggered by the user when a specific movement of the patient's trunk is detected, using an IMU fitted within the jacket worn by the patient (visible on Figure 1.3d).

Below, we describe the passive walking controller of Atalante and give a brief overview of the humanoid robots' literature which inspired its development.

Passive walking control architecture: plan then execute

The primary control paradigm of Atalante is to plan then execute. It takes the form of a passive walking controller, inspired by the humanoid robots' literature, and is composed of three elements (illustrated on Figure 1.5):

10 The estimation and control of the deformations are beyond the scope of the presented work, and we refer the interested reader to the work in [Benallegue2014;Benallegue2015;Mifsud2016] for humanoid robots, and in [Vigne2021] for a thorough study on Atalante. • a trajectory of the exoskeleton's joints and underactuated DoFs, computed offline;

• a stabilization algorithm, adapting this joint trajectory online to track the underactuated DoFs trajectory;

• a low-level controller, composed of joint-independent high-gain PD controllers, tracking the modified joint trajectory.

To follow the paradigm of humanoid robots, the patient is viewed as a set of rigid bodies, rigidly fixed to the robot, adding mass and inertia to the exoskeleton. The passive walking controller of Atalante is inspired by two broad families of control methodologies used to realize walking motions with humanoid robots11 . A first version has appeared in [Gurriet2018] 12 . Concerning trajectory design, the main sources of inspiration are the works presented below. The method of virtual constraints and Hybrid Zero Dynamics (HZD) is exposed in [Westervelt2007; Ames2014; Reher2016; Grizzle2017]. This control methodology is based on the offline design of optimal walking trajectories [Hereid2016], which possess some unstable dynamics that are stabilized online through the feedback of their timescale. The generalization of this methodology is still the subject of some recent work [Castillo2019;Castillo2020;Galliker2022]. A second family of controllers generates the trajectories of the robot online. Because the trajectory is computed online with these methodologies, it is better suited to react to perturbations in a Model Predictive Control (MPC) fashion, execute various motions, and modify the plan online. A line of work relies on some approximation of 1.2.

PASSIVE WALKING WITH A SELF-BALANCED EXOSKELETON 27

the dynamics13 of the robot to generate trajectories of the reduced state of the robot online, either as feedback controllers [Kajita2003] or as the solution to an optimal trajectory design (in the form of a QP) over some horizon [Wieber2006b;Herdt2010;Caron2016;Griffin2017], at a frequency typically lower than 100 Hz. The controller of Atalante uses the same offline design of optimal walking trajectories as [Hereid2016] but does not rely on the HZD stabilization Concerning stabilization, we employ the whole-body admittance algorithm [Caron2019], which enables position-controlled robots to regulate the unilateral contact forces and track the CoM trajectory. In this paradigm, the state trajectories are converted at high-frequency (typically 1 kHz) into an instantaneous acceleration in the joint-space through Inverse Kinematics (IK) 14 . The next step is to achieve these joint accelerations. In principle, this can be done by torque control or position control. However, mostly due to the patient, the torques are vastly uncertain. This discards the torque control strategies found on Atlas [Johnson2015], TORO [Englsberger2014] or Talos [Villa2022], and favor instead position control ones found on HRP-4 [Kaneko2011].

Finally, low-level controllers are implemented to achieve disturbance rejection and to track the previously defined positions. They take the form of independent high-gain PD controllers (running at 1 kHz). This decentralized approach shows satisfactory results due to the high reduction factor gearboxes (see [Finet2017] for more details).

In this passive walking controller, the patient efforts are not accounted for: they are only seen as unmodeled disturbances to be compensated for by the low-level controllers. In other words, all the patient deviations from the trajectory are hampered by the exoskeleton. Despite its merits for passive gait training, this control architecture is not suited for active gait training, which is the topic of this thesis. Reference trajectory tracking can be relaxed in various ways, leading to different types of freedom: either temporal freedom, spatial freedom, or both at the same time can be offered. We consider that a control law offers temporal freedom (respectively spatial freedom) if, through their efforts, a patient can substantially change the duration of the steps (respectively the shape of the steps' trajectory i.e. the geometric path followed during the steps). These two freedoms are illustrated on Figure 1.6. Most of the active gait training controllers found in the literature offer some form of both temporal and spatial freedom. For instance, in [Aoyagi2007; Banala2009; Vallery2009; Duschau-Wicke2010], a nearest-neighbor algorithm is used to offer temporal freedom by finding the closest point, within a reference trajectory, to the current position of the legs of the patient (in joint-space) 16 . This methodology is extended in [Martínez2018] to provide step-to-step variation of the step length, through online modification of the reference trajectory between each step. Another control law is proposed in [Martínez2019], based on a flow field rather than a potential field, in order to combine in a more native way the freedom provided by the controller from [Duschau-Wicke2010] and some form of assistance to the user.

Altering the trajectories of the system is troublesome as it may jeopardize the balance of the exercises. This problem is not critical in the works cited above because the exoskeletons under consideration are platform-based exoskeletons, where the safety of the exercise is never at risk, or crutched exoskeletons, which leave the physiotherapist or patient in charge of the balance of the exercise. With self-balanced exoskeletons, maintaining balance while allowing freedom is a seldom studied problem. This specific problem lies at the intersection of the literature on humanoid walking robots (see Section 1.2.2) and the one on robotic gait training (see above). While these two fields of research have received numerous contributions over the past 20 years, their intersection has only drawn little attention. This is certainly caused by the small number of available self-balanced exoskeletons 17 .

Exploiting spatial freedom with self-balanced exoskeletons

To our knowledge, only one active walking controller for self-balanced exoskeletons, demonstrated with Atalante, has been proposed in [Gurriet2020]. This controller, illustrated on Figure 1.6a, offers spatial freedom to the patient by relaxing the tracking of the trajectory on the swing hip and knee, and uses Control Barrier Functions (CBF) to limit the deviation of these joints from a nominal walk pattern. Specifically, safe regions are designed in joint-space s.t. the freedom of the patient is saturated at the beginning and the end of each step, and, therefore, the exoskeleton feet lift and land correctly.

This controller constitutes the first demonstration of an active walking controller for self-balanced exoskeletons. It offers spatial freedom around the offline trajectory, but no temporal freedom. In this thesis, we propose a rehabilitation controller for Atalante which offers temporal freedom around the offline trajectory.

Contribution of the thesis: exploit temporal freedom and maintain balance

Temporal freedom has the potential to make the patient participate more significantly. This point is illustrated on Figure 1.6 which reports the swing hip position as a function of time for various degrees of spatial (left) and temporal (right) freedom, respectively. The envelope of Figure 1.6b is larger than Figure 1.6a, hinting that temporal freedom would be easier to perceive for patients. Our objective in this thesis is to develop a new feature for Atalante: a rehabilitation controller giving substantial temporal freedom to the patient without requiring the physiotherapist's assistance.

We propose a control system update. This update has a dual purpose: it empowers the patient and monitors safety.

Overview of the solution The modifications are schematically depicted on Figure 1.7. We choose to let the patient physically contribute to the motion of the swing leg and, consequently, allow modulation of the velocity at which a predefined gait is traveled. Because this directly impacts safety, a new reference trajectory is computed for the low-level controllers of the whole system. This is done online by an optimization-based trajectory planning algorithm.

Detailed view

The two legs of Atalante are the subjects of very distinct changes. On the swing leg, we proceed as follows. The nominal gait for the two-legged system, generated using the same offline optimization methodology 1.3.

ACTIVE GAIT TRAINING WITH SELF-BALANCED EX-OSKELETONS 31

as for the passive walking, serves to define a geometric path for the swing leg. The exoskeleton efforts in the longitudinal direction of the path are nullified and left for the patient to produce, while the robot motion is strictly controlled in the hyperplane orthogonal to the path. To this end, we rely on the Virtual Guides (VG) methodology which maps the high-dimensional user efforts to a one-dimensional quantity: the velocity at which the swing leg's geometric path is followed. Implicitly, this defines a new schedule for the nominal trajectory, the patient schedule, which is followed as long as safety is not threatened.

On the support leg, the control structure is also updated. Because the swing leg DoFs have been updated (see above), only the support leg can be used for the admittance controller. This admittance controller successfully tracks reference trajectories that satisfy the unilateral contact constraints.

To define a reference trajectory synchronous with the patient schedule, a simple time rescaling of the nominal gait would be possible but will most likely violate the aforementioned constraints. Instead, we chose to adopt an online planning strategy.

The online planning strategy we consider solves an Optimal Control Problem (OCP) over an unspecified horizon for a Linear Inverted Pendulum (LIP) model, which represents the overall balance dynamics of the system (patient+exoskeleton). The patient schedule defined above is treated as a penalty on the final time. The unilateral contact constraints is a constraint on the CoP. To ensure the feasibility of the next step, it is sufficient to require that the trajectory endpoint satisfies some geometric constraints. These are included in the OCP. The simplicity of the LIP model allows us to use a fast resolution method for this nonlinear OCP, relying on the theoretical study about the nature of the set of feasible times for the LIP dynamics (i.e. the set of final times for which a solution to the OCP exists).

Note that, although it is known the Atalante platform undergoes deformations, we do not try to compensate them in the swing or support leg control algorithms, nor the replanning. This choice is made in order to alleviate some complexity to the already complex problem at hand. However, these deformations might affect the performance of our algorithms when testing them on the actual platform Atalante.

Thesis outline This manuscript presents the development of the solution described above and experimental results stressing its performance.

CHAPTER 1. INTRODUCTION

Chapter 2: Atalante: models and passive control laws. This chapter presents in detail the experimental platform Atalante and its control methodology during passive walking.

Chapter 3: Empower the patient for gait training. This chapter presents the controller enabling temporal gait training. In particular, the split of the exoskeleton's DoF in half is detailed. The poor performance of a simple time rescaling strategy is exposed.

Chapter 4: Monitor safety. This chapter presents the replanning methodology developed to ensure the balance of the system. The method for fast solving of the OCP, compatible with real-time requirements on the Atalante's embedded controller, is provided. A theoretical study of the nature of the solution space is provided in Appendix A. Extensive experimental results are presented.

Chapter 5: Conclusion and perspectives. This chapter stresses the benefits and limitations of the proposed active gait training controller. Possible extensions are presented.

Publications in this thesis

The work done in this thesis led to the following published results: In this chapter, we describe the hardware and the passive walking control law of the exoskeleton Atalante. Contrary to walking rehabilitation, during which the patient efforts on the exoskeleton are substantially impacting the motion of the robot, in passive walking the controllers completely assign the motion of the robot, no matter the level of patient efforts.

• Maxime Brunet,
We start in Section 2.1 by describing the mechanical structure, actuators and sensors of Atalante. Then, we derive the equation of motions under the assumption that the patient motion is completely assigned by the exoskeleton. The governing equations are similar to those of a humanoid robot. Next, in Section 2.2, we describe how trajectory generation for the patient-exoskeleton system can be cast into a non-linear Optimal Control Problem (OCP). A direct collocation transcription of the OCP is described. It is used offline. We detail two specific instances of OCPs:

• a flat-foot trajectory, for which the supporting leg switches instantaneously from one leg to the other. There is no double support phase in between two single support phases, and the support foot is always resting flat on the ground.

• a foot-rolling trajectory, for which there exist double support phases. During these, the supporting feet are not resting still on the ground but are rolling on their front (for the trailing foot) or back edges (for the leading foot). The support foot is resting flat on the ground during single support phases.

The flat-foot trajectory is simpler to study but has a lower average translational velocity and is less anthropomorphic than the foot-rolling trajectory. It will serve as our use-case trajectory throughout Chapter 3 and Chapter 4. The foot-rolling trajectory will be used in the experiments discussed in Chapter 5. Finally, in Section 2.3, we expose how these offline trajectories are tracked by a stabilization control law based on a linearized model of the system, relying on position-based impedance control [Hogan1984;Heinrichs1996], also called admittance control.

2.1 Atalante: description and dynamical models

Definitions of the frames of reference

We consider three axes defined w.r.t. the patient, as depicted on Figure 2.1. The X-axis points forward, the Y-axis to the left, and the Z-axis is aligned with the local gravity vector, pointing upward. The local frame of reference is chosen lying on the ground, under one of the exoskeleton feet1 . The sagittal, frontal and transverse planes are the three orthogonal planes intersecting at the patient's CoM as depicted on Figure 2.1.

Description of Atalante

The Atalante exoskeleton is depicted on Figure 2.2. It is composed of 13 rigid bodies linked together by twelve actuated revolute joints (in blue). Except for the two Henke ankle axes, all of them are orthogonal to one of the planes defined on Figure 2.1 when the exoskeleton is standing still. The axis of each knee joint is orthogonal to the sagittal plane. The three axes of the three hip actuators of each hip are orthogonal to one of the sagittal, frontal or transverse planes. Each ankle possesses two joints: one orthogonal to the sagittal plane, and one aligned with the morphological Henke axis2 . The rigid body connecting the two legs of the exoskeleton is called the back body. The upper-leg part (between the sagittal hip joint and the knee joint) is called the thigh body. The lower-leg part (between the knee joint and the sagittal ankle joint) is called the shin body. The metal plate in contact with the ground which supports the patient's foot is called the foot body.

Atalante is equipped with several types of sensors (see on Figure 2.2). In details, encoders are sensing the angular positions3 of each joint of Atalante. Single-axis force sensors are located under each foot of the robot (4 sensors under each foot) and sense the unilateral contact forces with the ground. Several Inertial Measurement Unit (IMU), rigidly fixed to the exoskeleton, are sensing the proper acceleration and angular velocity of several bodies of the exoskeleton.

The patient is fastened using a special vest and straps for the feet, shins and thighs, visible on Figure 2.2.

Kinematic model

The fastening system of Atalante completely assigns the positions of the lower limbs of the patient w.r.t. the robot. The patient's torso is also fastened to the exoskeleton, but some freedom in the upper-body motion is allowed (in particular in their head and arms motion). We neglect the impact of these disturbing motions on the dynamics of the patient-exoskeleton system and simply consider that the patient is an additional rigid body attached to the exoskeleton.

According to the description above, we model the patient-exoskeleton system as an articulated rigid-body system of mass m (the patient mass plus the exoskeleton's) with 12 actuated Degrees of Freedom (DoF), corresponding to the joints of the exoskeleton, aligned with the patient's, and 6 unactuated DoF, defining the position and orientation of the back of the exoskeleton, w.r.t. the inertial frame W.

In details, following [Featherstone2008], we denote q ∈ R 18 the vector of generalized positions, (q sw , q sp ) ∈ R 6 × R 6 the vector of actuated positions for the swing leg and the support leg, and q un ∈ R 6 the vector of unactuated positions, s.t.

q =    q sw q sp q un    (2.1)
Likewise, we denote q ∈ R 18 (resp. q ∈ R 18 ) the vector of generalized velocities (resp. acceleration), ( qsw , qsp ) ∈ R 6 × R 6 (resp. (q sw , qsp ) ∈ R 6 × R 6 ) the vector of actuated velocities (resp. acceleration) for the swing leg and the

CHAPTER 2. ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER

support leg, and qun ∈ R 6 (resp. qun ∈ R 6 ) the vector of unactuated velocities (resp. acceleration), s.t.

q =    qsw qsp qun    q =    qsw qsp qun    (2.2)
We employ the general model of the displacement of a rigid body B w.r.t. to a frame W , living in the Special Euclidean Group SE(3). It is defined by the position O B ∈ R 3 of the origin of B, and the rotation matrix

R W B ∈ SO(3). This displacement is represented by a homogenous transformation matrix, denoted X W B ∈ SE(3), which writes X W B = R W B O B 0 1×3 1 (2.3)
In the rest of the thesis, X W B is referred to as a transform. The Special Euclidean Group SE(3) is equipped with the composition operator × • s.t.

× • : X W L × X L B ∈ SE(3)×SE(3) → X W L × • X L B ≜ X W B ∈ SE(3) (2.4)
For simplicity of notations, we omit this composition operator for the rest of the thesis.

From the vector of generalized positions, and a kinematic model defining the kinematic chain of the robot, we define the Forward Kinematics (FK) function as the map

FK : q ∈ R 18 → X W B (q) ∈ SE(3) (2.5)
We denote τ ∈ R 12 the vector of exoskeleton joint torques to be chosen by the controllers, (τ sw , τ sp ) ∈ R 6 × R 6 the vectors of exoskeleton swing and support joint torques, and (τ sw u , τ sp u ) ∈ R 6 × R 6 the vectors of swing and support joint torques created by the patient, s.t.

τ = τ sw τ sp τ u = τ sw u τ sp u (2.6)
The Center of Mass (CoM) of the patient-exoskeleton system is denoted c ∈ R 3 , and the total angular momentum of the system w.r.t. the CoM, in the inertial frame W, is denoted L s.t.

c = k m k O k k m k L = k (O k -c) × m k Ȯk + I k ω k (2.7)

ATALANTE: DESCRIPTION AND DYNAMICAL MODELS

39

with O k ∈ R 3 (resp. Ȯk ∈ R 3 ) the position (resp. velocity) of the origin of the frame F k attached to the part k ∈ [0 . . . 13] of the robot, w.r.t. the inertial frame W, ω k ∈ R 3 the rotation velocity of the frame F k w.r.t. the inertial frame W, m k ∈ R + the mass of the part k (including the patient mass), and I k ∈ M (3,3) the inertia matrix of the part k of the system (including the patient inertia matrix), and × the cross product of R 3 . We denote . x,y the 2D vector representing the position of 3D points projected onto the horizontal plane of reference. Analogously, we denote . z the height of 3D points, or the vertical component of 3D vectors.

Finally, we denote ξ the Divergent Component of Motion (DCM), defined as follows

ξ = c x,y + ċx,y ω ∈ R 2 (2.8)
with ω ≜ g c z , and p the Center of Pressure (CoP), defined as

p ≜ p x,y i f z i f z i ∈ R 2 (2.9)
with f i ∈ R 3 the external force vectors acting on the patient-exoskeleton system, in the inertial frame W, and p i ∈ R 3 the location of application of the external forces, in the inertial frame W. Note that, by definition, the CoP belongs to the convex hull of all the contact points, also called the support polygon SP. These quantities are of particular importance for the design of the stabilization control law detailed in Section 2.3.

Lagrangian dynamics

With the above notations, the Lagrangian dynamics write, following [Feath-erstone2008] 4 ,

M (q)q + C(q, q) = H(τ + τ u ) + i J i (p i , q) T f i (2.10)
with M (q) the generalized inertia matrix of the system, C(q, q) the combined gravity and inertia effects vector, H ≜ I 12 O 6 ⊤ the matrix mapping the joint torques to the actuated generalized coordinates, J i (p i , q) the associated Jacobian matrices at each contact point p i .

The Lagrangian dynamics (2.10) are identical to the one of walking robots, except for the vectors of patient joint torques τ u . These dynamics CHAPTER 2.

ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER

are instrumental in the offline trajectory generation framework described in Section 2.2, where patient torques are set to zero to model a paraplegic patient5 .

Newton-Euler dynamics

The Newton and Euler equations [Wieber2006a], in the inertial frame W, are, in accordance with (2.10),

m(c + g) = i f i , L = i (p i -c) × f i (2.11)

Offline trajectory generation

We now describe the Optimal Control Problem (OCP) formulation describing the walking of Atalante, the transcription of this OCP into a discrete optimization problem, and a numerical method to solve it. Then, two illustrative OCP examples are discussed: a simple flat-foot OCP and a more advanced foot-rolling OCP.

Generic walking motion optimal control problem

More details on the content of this section can be found in [Hereid2016;Gurriet2018].

OCP formulation

The problem of multi-contact trajectory generation is an OCP for hybrid dynamic. Each contact phase, defined by a different set of contact bodies (i.e. one or two feet resting flat or on one edge on the ground), defines distinct dynamics. These different contact phases, also called domains, are separated using several time instants t i . The OCP allows to optimize the vector τ of joint torques, the timing of each step t i , the contact forces f i and application point p i to minimize a general cost function subjected to constraints of various nature. The general OCP defining walking motion of humanoid robots is of the following form [Goswami2019, Section VII], which is also considered for Atalante, Problem 2.1 (Generic multi-step walking motion OCP).

min τ,t j ,f i ,p i n j=1 t j t j-1 Φ j (q(t), q(t), τ (t))dt + n j k=0 Φ j,k (t j , q(t j,k ), q(t j,k ), τ (t j,k )) s.t. for j = 1 . . . n, ∀t ∈ [t j-1 , t j ] q(t) = f (q(t), τ (t), p i , f i (t)) (2.12a) g j (t, q(t), q(t), τ (t), p i , f i (t)) ≥ 0 (2.12b) for j = 1 . . . n, q(t - j ) = q(t + j ) + k q (t - j , q(t - j ), q(t - j )) (2.12c) q(t - j ) = q(t + j ) + k q(t - j , q(t - j ), q(t - j )) (2.12d)
r cyc j (q(t j-1 ), q(t j-1 ), q(t j ), q(t j )) = 0 (2.12e)

for j = 1 . . . n, for k = 1 . . . n j , r eq j,k (t j , q(t j,k ), q(t j,k ), τ (t j,k ), p i , f i (t j,k )) = 0 (2.12f) r inj,k (t j , q(t j,k ), q(t j,k ), τ (t j,k ), p i , f i (t j,k )) ≥ 0 (2.12g) with f (q, τ, p i , f i ) = M -1 (q) -C(q, q) + Hτ + i J i (p i , q) T f i (2.13)
where (2.13) is obtained from (2.10), expressed in the forward sense (i.e. with torques as inputs), where the patient torques τ u are identically null.

The cost functions Φ i are running costs and are used, for instance, to guide the OCP toward solutions with the best energy efficiency, desired walking velocity, and those with the highest margins (for instance in the sense of penalizing solutions requiring contact forces generating large moments along the vertical axis of the inertial frame W). The cost functions Φ i,j favor solutions passing through waypoints at specific times6 and are used to produce trajectories with their feet rolling on edges between taking off, rising their feet high enough, or avoid known obstacles like stairs.

The first type of constraints included in Problem 2.1 are called running constraints and must be satisfied at all times. The running constraints (2.12a,2.12b) are used to express the physical constraints acting on the system, like the dynamics of the system, the friction cone constraints (which include the unilateral contact constraints), the joint position limits, or CHAPTER 2.

ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER

the actuation torque limits. They can also be used to incorporate constraints like minimal foot clearance, bounded CoM height, or bounded torso orientation. The second type of constraints included in Problem 2.1 are waypoints constraints and hold only at specific times 6 . Among them, constraints (2.12c) and (2.12d) can be incorporated to reflect the impact dynamics, discrete by nature. However, no impact map (k q , k q) has been implemented to generate trajectories for Atalante so far, and, simply,

k q ≜ 0 k q ≜ 0 (2.14)
A special constraint can be included to enforce the periodicity through (2.12e).

It is used to generate trajectories that can be infinitely looped onto themselves by symmetrization. Finally, equality and inequality waypoint constraints (2.12f) and (2.12g) are used to shape trajectories 6 by placing feet or CoM at specific locations, constraining the total footprint of trajectories, or bounding the velocity of some bodies of the system.

OCP transcription

The OCP defined in continuous time by Problem 2.1 is transcribed into an optimization problem using a finite number of variables using the Direct Collocation method [Hargraves1987;vStryk1993]. Direct multiple shooting methods [Bock1984; Mombaur2001; Hereid2015] have first been used to solve non-linear OCPs of the type of Problem 2.1. However, they have proved hard to scale with the increasing number of variables (arising for robots with a high number of DoF and/or with numerous domains). Instead, Direct Collocation has been found to be robust and scalable. The results reported in [Huynh2021] constitute an example of successful optimization of trajectories for Atalante involving more than 10 independent domains. The general idea behind Direct Collocation is to approximate the trajectory at a finite number of points x i , called collocation points, replacing the explicit forward integration of the dynamical systems with a series of integration constraints. For the trajectory represented by these points x i to be an approximation of the continuous solution of OCPs, a piece-wise polynomial interpolation is used to estimate the continuous solution between collocation points. Direct Collocation is generally used with defect constraints on defect variables δ. Defect variables δ are additional optimization variables that relax equality constraints at collocation points. Defect constraints are bounds on the defect variables δ (typically a small value). Introducing defect variables and constraints, instead of closed-form constraints (i.e. directly expressed as a function of the trajectory's optimization variable), has several advantages, like being faster to compute, decoupling the optimization into many independent subproblems, and allowing the solver to violate the equality constraints during the first iterations of the optimization. We refer the interested reader to [Hargraves1987;vStryk1993;Duburcq2022a] for more details about Direct Collocation and its implementation.

OCP resolution

The resulting discrete optimization problem is a relatively large Non Linear Program (NLP). Typical setup includes up to 5000 variables in [Huynh2021]. However, it is very sparse [Hereid2016], which allows it to be solved efficiently by solvers such as IPOPT [Wächter2006].

The output of the trajectory generation framework is a trajectory

Q Q : t ∈ [0, T f ] → Q(t) ∈ R 18 (2.15)

Flat-foot trajectory optimal control problem

The general multi-step OCP defined by Problem 2.1 can be further specialized to generate trajectories of interest. A first example of such specialization gives the flat-foot trajectory. It is composed of 5 different contact phases:

• STARTD: a starting double-support weight transfer phase, during which both feet are constrained to rest flat and still on the ground, the weight of the system is gradually shifted toward the starting stance foot, and the initial joint position is constrained (corresponding to a standing posture).

• STARTS: a starting single-support phase, during which only the starting stance foot is constrained to rest flat and still on the ground, and the CoM and flying foot are both moved forward.

• CS: a cyclic single-support phase, during which only the left foot is constrained to rest flat and still on the ground, the CoM and right flying foot are both moved forward and the initial and final states are constrained to be symmetric. Gathering this phase and its symmetric constitute a cyclic trajectory, which can be repeated as much as requested at run-time. • STOPS: a stopping single-support phase, during which only the starting stance foot is constrained to rest flat and still on the ground, and the CoM and flying foot are both moved forward.

• STOPD: a stopping double-support weight transfer phase, during which both feet are constrained to rest flat and still on the ground, the weight of the system is gradually shifted from the stopping stance foot to the middle of the feet, and the final posture is constrained (corresponding to a standing posture).

In addition, continuity constraints are added so that the whole motion is continuous, as illustrated on Figure 2.3. The particularities of this OCP are that the CoP is always constrained to lie within a box smaller than the actual support polygon SP (to add safety margins), the step duration is arbitrarily fixed to 1.0 s for single support phases and 0.5 s for double support ones, the step length is fixed (approximately to 15 cm), and the feet are constrained to take off and land horizontally to the ground. Thus, the average translational velocity is approximately 15 cm/s.

The CoM, CoP and DCM of the trajectory solution of the flat-foot OCP are reported on Figure 2.4, and the swing foot clearance during the (CS) phase is given on Figure 2.6 (Left).

Foot-rolling trajectory optimal control problem

A second example of specialization of Problem 2.1 gives the foot-rolling trajectory generation OCP, generating a more complex but more anthropomorphic walking trajectory. In detail, most of the OCP is similar to the one described above for the flat-foot trajectory. The major difference is the addition of a cyclic double support domain (CD), described below (and depicted on Figure 2.5):

• CD: a cyclic double-support weight transfer phase, during which the heel of the leading foot and the toe of the trailing foot are constrained to rest still on the ground, but their pitch angle is not constrained, s.t. the feet of the robot are rolling on their edges while the weight of the system is gradually shifted from the trailing foot to the leading foot, and the initial and final states are constrained to be symmetric.

For all single support domains defined in Section 2.2.2, the swing foot is no longer constrained to land flat on the ground at the end of the domain, but rather on its heel. The step timing is arbitrarily fixed to 0.81 s, for single support phases, to 0.23 s, for the (CD) phase, and to 0.

CENTER OF MASS STABILIZATION AND ADMITTANCE CONTROLLER47

Center of Mass stabilization and admittance controller

The two trajectories presented in the previous section satisfy the Lagrangian dynamics (2.10) and all the constraints acting on the system. However, due to the unstable nature of the dynamics (see the eigenvalues of (2.17) for instance), they can not be directly executed on the robot as pure open-loop joint torque trajectories, nor can they be naively followed by independent joint controllers, without leading to the divergence of the generalized coordinates q and, eventually, falling. For this reason, a stabilization algorithm, adapted from the humanoid robots' literature, ensures the correct execution of the trajectories presented above on Atalante7 .

In details, the stabilization algorithm tracks the horizontal CoM state trajectory. This stabilization algorithm is a feedback law designed using a linearized model of the patient-exoskeleton system, which performs a particular feedback on the CoM through the contact forces. The contact forces are in turn regulated by the so-called admittance control method.

CENTER OF MASS STABILIZATION AND ADMITTANCE CONTROLLER49

DCM control using admittance

The linearized model of the system from which the feedback law stems is derived for the Newton-Euler dynamics (2.11). Following [Wieber2016], dynamics (2.11) can be simplified into the Linear Inverted Pendulum (LIP) model. First, notice that (2.11) yields

mc × (c + g) + L m(c z + g z ) = p i f z i f z i (2.16)
Along gait patterns of moderate velocity, the angular momentum variations are small and can be neglected, i.e. L = 0. Further, assuming that the robot walks on horizontal ground and that the CoM remains at a constant height c z , the X and Y dynamics of (2.16) simplify to the following uncoupled LIP dynamics cx,y = ω 2 (c x,y -p) (2.17)

where ω= g c z is the natural frequency, and p is the CoP defined by (2.9). The CoP acts as bounded virtual8 actuator in (2.17) (recall that, by definition, the CoP is bounded inside the support polygon SP).

The dynamics (2.17) reproduce the unstable nature of the system. For stabilization, a stabilizing state feedback law on the CoM, such as the one detailed in [Caron2019], is used to drive c * -c to 0. First, one defines

p d = p * -(1 + k p ω )(ξ * -ξ) - k i ω (ξ * -ξ) + k d ( ξ * -ξ) (2.18)
with ξ the DCM, readily computed from the horizontal LIP state (c

Hierarchical Inverse Kinematic problems: the Stack-of-Tasks framework

Inverse Kinematics (IK) problems are ubiquitous in the humanoid robots' literature. They naturally arise because objectives are best formulated in spaces (like the Cartesian inertial frame W) which are different from the actuators' space (the joint space). Hence, the FK map (2.5) needs to be inverted to generate trajectories to be followed in the joint space. Mathematically, if the FK map were invertible, the solution q * would write

q * = FK -1 (X d ) (2.20)
with X d the desired transform of some rigid body of the robot. However, the FK map is not injective due to kinematic redundancy. A classical workaround is to define objectives as task errors and to ensure they are asymptotically driven to zero. In details, following the task-function approach [Samson1991], a task function e i is any twice differentiable function that maps the robot configuration space to R n . Its image space is called the task space. The function e i may define various objectives. We detail here a few important task functions.

The task function e q is expressed directly from the robot's generalized coordinates e q : q ∈ R 18 → q d -q ∈ R 18 (2.21) with q d ∈ R 18 the desired configuration, which is time-varying. The task function e x ∈ R 3 in the literature is written w.r.t. the position x of a point,

CENTER OF MASS STABILIZATION AND ADMITTANCE CONTROLLER51

like the CoM, in Cartesian inertial frame W

e x : q ∈ R 18 → x d -x(q) ∈ R 3 (2.22)
with x d ∈ R 3 the desired location, which is time-varying. Finally, the last important task function e X is expressed from the transform of one of the robot's rigid bodies, in the sense of (2.3), and using the matrix log map, which maps the Lie group SE(3) to the Lie algebra se(3

) = R 6 , e X : q ∈ R 18 → log(X d FK(q)) ∈ R 6 (2.23)
with X d ∈ SE(3) the desired transform, which is time-varying. For every task functions e i , we define its Jacobian J i as follows ėi = J i (p i , q) q (2.24) which yields ëi = Ji (q) q + J i (p i , q)q (2.25)

For each task function e i , the Jacobian J i and its time derivative Ji , defining the objective w i , are computed using the efficient library Pinocchio [Carpen-tier2019]. From the task acceleration (2.25), we define the objective w i for all task functions e i as follows

w i (q) = ||ë i -ë * i || 2 = ||J i (p i , q)q -(ë * i -Ji (q) q)|| 2 (2.26)
with a target task acceleration ë * i , illustrated on a Cartesian task, driving the task error e i exponentially to zero, such as

ë * i = ẍd -k p i e i -k di ėi (2.27)
with k p , k d ≻ 0 and ẍd denoting the feedforward in the Cartesian space. The desired motion of the robot is usually defined by several objectives w i , s.t. the total objective w writes

w(q) = i λ i w i (q)
(2.28) with λ i ∈ R + weights expressing the relative importance of the objectives w i .

The IK problem is to find optimal joint accelerations q which minimize w.

In addition to these objectives, the physical limitations of the robot, such as the joint position, velocity and acceleration bounds, must be embedded into the IK problem. Finally, it writes q * = arg min s.t.

Aq ≤ B i λ i w i (q) (2.29) CHAPTER 2

. ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER

with A (resp. B) a matrix (resp. vector) to represent the aforementioned bounds. While the IK problem (2.29) has been successfully used to realize stable walking of humanoid robots [Caron2019], we favor a different approach. Indeed, the tuning of the weights λ i can be cumbersome, and some objectives often have a strict priority over others (for instance, the robot's feet should always be on the ground, while the posture tasks are used only for regularization). An alternate formulation has been proposed, the Stack-of-Tasks (SoT) framework [Mansard2009], enabling one to solve IK problems composed of objectives listed in strict hierarchical order. A Stack-of-Tasks (SoT) with n ≥ 0 levels is an ordered list of objectives w j i , grouped in several levels j, with w j their cost defined as follows

w j (q) = i λ j i w j i (q) (2.30)
The levels are sorted in strict decreasing priority order, s.t. the cost w j can only be decreased if the value w k , of all the levels k < j ≤ n, remains constant (i.e. equal to the optimal cost w k * ). Mathematically, a Hierarchical IK problem with constraints writes ∀j = 0 . . . n, q * =arg min i

w j i (q) s.t. Aq ≤ B ∀k < j, w k (q) = w k * ≜ min q i w k i (q) (2.31)
On Atalante, the SoTs are solved in real-time using state-of-the-art hierarchical optimization routines, such as the one proposed in [Escande2014]. The key idea is to solve the Hierarchical IK problems (2.31) iteratively, starting with the highest priority objectives, and then minimizing the objectives of lower priority in the level set of the higher priority objectives already minimized. We briefly illustrate below this iterative process.

For ease of notations, let us define a simplified Hierarchical IK problem, with only two levels (n = 1), one task at each level, and no constraint. For example, one has q * = arg min w 1 0 (q) s.t.

w 0 (q) = w 0 * (2.32)

CENTER OF MASS STABILIZATION AND ADMITTANCE CONTROLLER53

with w 0 (q) = w 0 0 (q) = ||J 0 0 (q)q -b 0 || 2 b 0 = ë0 * 0 -J0 0 (q) q (2.33) and w 1 (q) = w 1 0 (q) = ||J 1 0 (q)q -b 1 || 2 b 1 = ë0 1 * -J1 0 (q) q (2.34)
As there is no constraint, in this simple case, all the solutions q0 * minimizing w 0 (q) write q0 * = J 0 0 † (q)b 0 + P 0 q1 (2.35)

with . † the matrix pseudo-inverse, P 0 ≜ (I -J 0 0 † (q)J 0 0 (q)) a projector on the null-space ker(J 0 0 (q)) of J 0 0 (q), and q1 an arbitrary vector. If the kernel ker(J 0 0 (q)) is not null, then one can attempt to minimize the cost of the level 1 problem, of lower priority.

The level 1 problem (2.34) is solved in level 0's kernel by selecting q1 only. A direct resolution of this simple case yields

q1 * = J 1 0 (q)P 0 † b 1 -J 1 0 (q)J 0 0 † (q)b 0 (2.36)
and the solution of problem (2.34) write q * = J 0 0 † (q)b 0 + P 0 q1 * (2.37) Equation (2.37) and (2.36) show how one can iteratively solve a simple Hierarchical IK problem such as (2.32). For more complex problems, we refer the interested reader to [Escande2014] to get efficient solvers addressing them.

The nominal SoT 1 of Atalante

The goal of Atalante's nominal SoT, denoted SoT 1 , is to stabilize the system around the nominal trajectory, computed offline as described in Section 2.2. To this end, most objectives included in SoT 1 are the tracking of quantities defined by the nominal trajectory, except for the task at level 1, which is defined by the admittance horizontal acceleration (2.19).

Namely, the level 0 of SoT 1 , hence the highest priority, includes two tasks on the robot's feet,

w 0 i (q) = ||J 0 i (q)q -(ë * 0 i -J0 i (q) q)|| 2 (2.
e 0 i = log(X 0 nom i X 0 i (q)) (2.39)
where i = 0 (resp. i = 1) denotes the left foot (resp. the right foot), X 0 i is the left (resp. right) foot transform w.r.t. the inertial frame W, J 0 i (q) is the left (resp. right) foot transform Jacobian, and K 0 p , K 0 d ≻ 0 are two matrices 10 . With these notations, level 0 writes

w 0 (q) = w 0 0 (q) + w 0 1 (q) (2.40)
The level 1 of SoT 1 has only one task, defined by (2.19)

w 1 (q) = ||J 1 c x,y (q)q -(c d -J1 c x,y (q) q)|| 2 (2.41)
where J 1 c x,y (q) is the horizontal CoM transform Jacobian. Recall, that

cd = c * + A(p -p d ) (2.42)
with c * is horizontal CoM acceleration of the nominal trajectory.

The level 2 of SoT 1 has only one task, concerning the robot's back orientation,

w 2 (q) = ||J 2 back (q)q -(ë * back -J2 back (q) q)|| 2 (2.43)
defined by the rotation task function

e back = log 3 (R 2 nom back R 2 back (q)) (2.44)
where R 2 back is the orientation of the robot's back w.r.t. the inertial frame W, J 2 back (q) is the Jacobian of the orientation of the robot's back, and K 2 p , K 2 d ≻ 0 are two high gain matrices.

Finally, the level 3 of SoT 1 is composed of one task, expressed on the robot's generalized position q, defined as follows

w 3 (q) = ||q -(q nom + K 2 p (q nom -q) + K 2 d ( qnom -q))|| 2 (2.45)
where (q nom , qnom , qnom ) are the generalized position, velocity and acceleration of the nominal trajectory. From all the objectives of SoT 1 described above, only the task at level 1 (see (2.42)) contains a term not only expressed using the nominal trajectory. This suggests that the stabilization will be obtained by updating the horizontal position of the CoM, as can be visible in experiments on the hip lateral motion. SoT 1 can be summarized as follows, with levels listed in decreasing order of priority 10 Chosen as high gain. On Atalante, the Hierarchical IK problem SoT 1 is solved online by a solver inspired from [Escande2014]. Its output is the optimal acceleration qt . This signal is integrated twice to get joint position and velocity targets (q t , qt ). Finally, these joint targets are tracked using a high-gain PD controller

CENTER OF MASS STABILIZATION AND ADMITTANCE CONTROLLER55

Stack of Tasks

τ sw τ sp = K p (q t - q sw q sp ) + K d ( qt - qsw qsp ) (2.46)
From the definition of SoT 1 and (2.46), the nominal controller is now fully described. It is schematically depicted on Figure 2 In this chapter, we propose a new controller for Atalante which enables the user to participate in the motion. More precisely, we adapt the existing controllers to 1. let the patient modify the swing leg's schedule 2. stabilize the CoM dynamics around its reference trajectory using the support leg.

We modify the existing controllers only during single support phases, and keep them unchanged during the double support phases 1 . We propose to use the Virtual Guides (VG) framework [Joly1995; Sanchez Restrepo2018], originally developed for teleoperation and co-manipulation, on the swing leg Degrees of Freedom (DoF). This requires some adaptations and yields a controller akin to a path control feedback [Banala2007; Duschau-Wicke2010] that aims at stabilizing a parametric curve in the articular space. Because the velocity along the curve is left free to choose by this controller, the patient is free to express an intent via the displacement of the system along this path. This displacement defines the so-called patient 3.1.

SATISFYING THE PATIENT INTENT WITH THE SWING LEG 59

schedule. On the other hand, we propose to use an adaptation of whole-body admittance, originally developed for the stabilization of humanoid robots during walking, for the support leg only. The purpose of the admittance task is to stabilize the CoM dynamics around a reference trajectory. To make this reference consistent with the swing leg motion, a natural idea is to apply a Time Rescaling (TR) of the nominal CoM trajectory following the patient schedule. The stabilizing properties of admittance around the nominal trajectory are not diminished when using it on the support leg only, as shown in Section 3.2.4. More precisely, its stabilization properties are preserved within a neighborhood of the nominal trajectory, when rescaled by the patient schedule, as shown in Section 3.3. The chapter is organized as follows. First, we describe the VG control law of the swing leg. Second, we describe our adaptation of the whole-body admittance control law, for its use on the support leg only, along with a natural synchronization strategy with the swing leg. Then, we prove that, provided the support leg has sufficient control authority (a sufficient number of DoF), and the motion of the swing leg is known, the use of the sole support leg is sufficient for the admittance task to succeed. Finally, we study, in simulation and experimentally, the stabilization property offered by our controller when the patient schedule differs from the nominal schedule.

Satisfying the patient intent with the swing leg

The VG methodology [Joly1995] allows a parametric curve P to be followed at a velocity prescribed by efforts not produced by the robot. This control methodology is classically used during teleoperation and co-manipulation in order to constrain a slave robot onto a desired manifold (e.g. in Cartesian space) and/or, more generally, help an operator guide the robot during a specific task.

The VG framework has two ingredients 1. a control law providing a contraction property in the direction orthogonal to the path. Usually, this control law is implemented as torques at the joint level.

2. a control law accounting for the patient efforts. This controller is implemented as a real-time re-scheduling of the nominal trajectory.

In a way, the VG methodology is a technological solution to realize a path controller, which can be compared to several alternative formulations [Ba-nala2007; Aoyagi2007; Vallery2009; Duschau-Wicke2010; Martínez2018; 60CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING Martínez2019]. All these path controllers are stabilizing around the closest point in the trajectory from the current state of the robot. We choose the VG methodology as it possesses the advantage of explicitly minimizing the robot forces acting on the patient, as will be detailed next. Also, its computation load is light, which is handy for embedded implementation.

Virtual Guides (VG) allows the patient's efforts onto the exoskeleton to drive the velocity of his swing leg in real-time within a predefined path (computed from the nominal swing leg trajectory). This approach enables us to satisfy the user intent without actually measuring the patient's efforts, which is convenient since Atalante does not possess torque or force sensors at the interfaces with the patient nor at the joint level.

Choice of the Virtual Guides definition space

Virtual Guides (VG) are defined in a certain task-space 2 , in which both internal efforts (robot efforts) and external efforts (in our case, patient efforts), are defined. In many occurrences, e.g. manipulators interacting with their environment at their end-effector, VG are defined in Cartesian task-space. This is not the case here, because exoskeletons interact with their users at every link. Indeed, the patient does not exert efforts only at the end-effector (the swing foot). This appears in the Lagrangian dynamics (2.10) of the patient-exoskeleton system, where the patient efforts appear at every joint. Besides, defining the control law of the exoskeleton in Cartesian space would not enforce inter-joint coordination (but joint space does) due to the noninjectivity of the joint to Cartesian position mapping. For these two reasons, we consider the VG approach in joint-space, similar to most of the pathcontrol laws for exoskeletons [Aoyagi2007; Vallery2009; Duschau-Wicke2010; Martínez2018; Martínez2019].

Construction of a parametric curve of reference

As already mentioned, Virtual Guides (VG) aim at constraining the motion of the robot along a parametric curve or manifold. The parametric curve used by our controller is built from a nominal joint trajectory T

T : t ∈ [0, T f ] → T (t) ∈ R 12 (3.1)
2 Here, any normed space enabling one to define error metrics (e.g. the Cartesian space R 3 , the Special Orthogonal group of rotation matrices SO(3), or the Special Euclidean group SE(3)).

SATISFYING THE PATIENT INTENT WITH THE SWING LEG 61

readily obtained from a trajectory Q computed by the offline trajectory generation pipeline presented in Section 2.2. The swing leg trajectory T sw is readily extracted from the complete trajectory T

T sw : t ∈ [0, T f ] → M i T sw (t) ∈ R 6 (3.2)
with M i ∈ R 6×12 , i ∈ {l, r} a masking rectangle matrix depending on the stance side (left or right). To obtain a parametric curve from the swing leg trajectory T sw , we reparametrize it w.r.t. its curvilinear abscissa s

s : τ ∈ [0, T f ] → τ 0 || Ṫ sw (t)|| 2 dt ∈ [0, L max ] (3.3)
with L max the total length (in joint space) of the swing leg trajectory. By this formula, the curvilinear abscissa s is monotonous, therefore, assuming further that the Euclidean norm || Ṫ sw (t)|| 2 is non-zero for all t, s can be inverted. Then, we define the parametric curve as

P ≜ s ∈ [0, L max ] → T sw • s -1 (s) ∈ R 6 (3.4)
To process this equation, one can numerically integrate the differential equation

ṡ(t) = || Ṫ sw (t)|| 2 (3.5)
Then, we gather the two grid vectors t = { i N T f , ∀i ∈ [0, N ]} and s = {s(t i ), ∀i ∈ [0, N ]}, and we fit 3 rd order Bézier curves (Splines) allowing us to evaluate s(t) for arbitrary values of t. The same fitting method is used to evaluate s -1 (σ) for arbitrary values of σ. An example of the result representing the nominal trajectory is depicted on Figure 3.1, which stresses that the mapping t → s is monotonous, invertible but non-linear.

Virtual Guides controller

In the following

σ : t → σ(t) ∈ [0, L max ] (3.6)
is now used as a control variable. It defines the current set-point P (σ(t)) along the swing leg path. Schematically, the VG methodology constrains the position x of the robot onto the parametric path P , as pictured on Figure 3.2. At the current point P (σ(t)), the Frenet-Serret unit tangent vector to the curve P , pointing in the direction of motion is The exoskeleton joint torques of the swing leg τ sw are computed as the sum of two torque vectors. First, a (vector-valued) torque τ sw P D aiming at providing a contraction property in the direction orthogonal to T , computed as a high-gain PD controller

T (σ) ≜ dP dσ (σ) (
τ sw P D (σ, σ, q sw , qsw ) = K sw p (P (σ) -q sw ) + K sw d (T (σ) σ -qsw ) (3.8)
with K sw p , K sw d ≻ 0 constant diagonal gain matrices. Second, a gravity compensation torque τ sw GC , aiming at compensating the gravity effects on the robot, computed by inverting the dynamics (2.10) of the exoskeleton (using a holonomic fixed support-foot constraint [Featherstone2008]), s.t.

τ sw GC (q) = H sw H M [I -J M ] H ⊤ † H M [I -J M ] C R (q, 0) (3.9) J M = J ⊤ (q)(J(q)M -1 R (q)J ⊤ (q)) -1 J(q)M -1 R (q) (3.10) H M = HM -1 R (q) (3.11)
with M R the generalized inertia matrix of the robot only and C R (q, 0) the gravity effects vector of the robot only. Finally, the swing leg τ sw exoskeleton 3.1.

SATISFYING THE PATIENT INTENT WITH THE SWING LEG 63

Figure 3.2: Schematic description of the VG approach.

torques are computed as

τ sw (σ, σ) = τ sw P D (σ, σ, q sw , qsw ) + τ sw GC (q) (3.12)
This PD plus gravity compensation controller is very similar to the pathcontrol controllers in [Banala2007; Duschau-Wicke2010], the main difference with these lies in the choice of the control variable. Following our VG approach, we define the control variable σ such that the projection of the efforts τ sw P D along the path P is nullified, which reads T (σ) ⊤ τ sw P D (σ, σ, q sw , qsw ) = 0 (3.13)

Notice that we only nullify the efforts of the PD controller, as we do not aim at canceling the gravity compensation efforts in any direction. This yields, by injecting (3.8) and solving for σ

[Sanchez Restrepo2018], σ(σ, q sw , qsw ) ≜ T (σ) ⊤ K sw p (q sw -P (σ)) + K sw d qsw T (σ) T K sw d T (σ) (3.14)
In summary, our VG controller boils down to implementing the joint torques (3.12) and integrating (3.14) to get σ from σ.

Behavior of Virtual Guides

In this section, we illustrate the difference between the above-defined VG controller and the path-control nearest-neighbor approach. This approach, detailed in [Duschau-Wicke2010], consists in solving the following minimization

64CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING problem σ nn = arg min σ ||P (σ) -q sw m || 2 (3.15)
where q sw m is the vector of measured positions of the swing leg's joints. We call "nearest neighbor" the position q sw nn = P (σ nn ) corresponding to σ nn defined by (3.15).

Almost a closest point algorithm: toy examples

In this section, we illustrate the asymptotic behavior of the VG dynamics (3.14). To this end, we study the behavior of σ defined by dynamics (3.14) with arbitrarily fixed values of q sw and qsw . We denote these fictitious joint positions and velocities as follows

∀t , q sw (t) = q sw f i , qsw (t) = qsw f i (3.16)
In addition, we choose K sw p , K sw d ≻ 0 as constant gain matrices (with K d the critically damped matrix)

K p = αI 6 , K d = 2 K p (3.17)
for some α ∈ R + . We study the convergence of (3.14) in the case where the measured position lies perfectly on the trajectory, and the measured velocity is null, i.e.

q sw f i = T sw (t f i ) , qsw f i = 0 6 (3.18)
with t f i ∈ [0, T f ], the fictitious time. The solutions of (3.14), integrated over 200 ms from several initial conditions σ 0 (s.t. t 0 ≜ s -1 (σ 0 )), and for several t f i , are reported on Figure 3.3. The final value of the solution σ is denoted σ f (i.e. t f ≜ s -1 (σ f )). These solutions illustrate that the convergence of σ to σ f i ≜ s(t f i ) is not guaranteed, which implies the existence of several equilibrium points to (3.14). By definition, the equilibrium points of the dynamical system defined by (3.14) are s.t.

T (σ) ⊤ K sw p (q sw -P (σ)) + K sw d qsw = 0 (3.19)
Interestingly, by using the weighted norm For an arbitrary function h defined over [0, L max ], we denote h its normalized version,

f ≜ ||K sw p (q sw -P (σ)) + K sw d qsw || 2 N (3.20) with N = K
h ≜ h -min [0,Lmax] h max [0,Lmax] h -min [0,Lmax] h (3.24)
We denote g(σ) ≜ ||P (σ) -q f i || 2 2 the Euclidean distance from q f i . On Figure 3.4, which depicts f , we verify that the local minima of f are indeed the equilibrium points of (3.14). In this example, with qsw f i = 0 and K p a constant diagonal matrix, Hence, ḡ is identical to f . Therefore, Figure 3.4 also reveals that not all the equilibrium points of (3.14) minimize the Euclidean distance, as the path-control methodology does. The dynamics (3.14) converge toward a local minimum, closer to the initial state.

f (σ) = ||K sw p (q sw -P (σ))|| 2 N = α 2 g (3.
The VG methodology also differs from the path-control approach because it is impacted by the measured velocity. In particular, a non-zero measured velocity qsw f i can shift its equilibrium points. This point is illustrated on Figure 3.5, where the VG equilibrium is slightly different from the closest point when the measured velocity qsw f i is not null, for instance when the velocity is nominal

qsw f i = Ṫ sw (t f i ) (3.26)
Further, in practice, the measured position is never precisely on the trajectory (i.e. ∄ t s.t. q sw f i = T sw (t)). We study the impact of this discrepancy by adding a component to q sw f i belonging to the hyperplane orthogonal to T (σ f i ), in the direction of the nominal acceleration projected in this hyperplane, i.e.

q sw f i = T sw (t f i ) + a f i T sw (t f i ) -T (σ f i ) ⊤ T sw (t f i )T (σ f i ) || T sw (t f i ) -T (σ f i ) ⊤ T sw (t f i )T (σ f i )|| 2 (3.27)
for some a f i ∈ R. Choosing a f i ̸ = 0 modifies both the VG equilibrium point and the nearest-neighbor, as depicted on Figure 3.6. Note that the impact on the equilibrium point is sensitive to the values of t f i , t 0 and a f i . On Figure 3.6a, the equilibrium is greatly modified, differing from both the Euclidean minimum distance and σ f i . However, on Figure 3.6b the impact of this change seems negligible, and the VG equilibrium is still equal to the nearest-neighbor, which is also equal to σ f i . Finally, all the previous analyses were conducted with isotropic matrices K p and K d (proportional to identity), which is never the case on the robot. We study the impact of anisotropic K p and K d with the example matrices below

K ′ p = Diag( 1e4 1e4 1e4 2e4 2e4 2e4 ) , K ′ d = 2 K ′ p (3.28)
The results for this choice of K ′ p and K ′ d are pictured on Figure 3.7. Interestingly, this modification changes the profile of f , but the nearest-neighbor is still an equilibrium of the VG, and the other VG equilibrium can be changed, as revealed on Figure 3.7b. Nonetheless, with the gains tuned for the exoskeleton K exo p and K exo d , all the VG equilibria are unchanged w.r.t. the isotropic matrices case, as shown on Figures 3.7c and 3.7d. In summary, the VG controller locally behaves similarly to a path control algorithm, when initialized close enough to the nearest neighbor. Additionally, it cancels efforts tangential to the path. Hence, as expected, Figure 3.8 shows that the reference set-points (P (σ), T (σ) σ) of the PD controller Eq. (3.8) closely follow the fictitious joint positions and velocities, even in the case when these positions and velocities are time-varying.

SATISFYING THE PATIENT INTENT WITH THE SWING LEG 69 (a)

t f i = 0.5 s, t 0 = 0.1 s, K ′ p , K ′ d (b) t f i = 0.8 s, t 0 = 0.5 s, K ′ p , K ′ d (c) t f i = 0.5 s, t 0 = 0.1 s, K exo p , K exo d (d) t f i = 0.

Virtual Guides behavior: experiment with a user

The VG behavior during an experiment with an able-bodied user inside the exoskeleton is depicted on Figure 3.9 and Figure 3.10. This experiment has been conducted using the complete algorithmic solution proposed in this chapter. An experimented user of the exoskeleton provides external help, holding the handles of the exoskeleton on its sides, to ensure the balance during the whole experiment and focus on the VG behavior. The support leg torques τ sp used during this experiment are yet to be specified, but we focus here on the behavior of the swing leg VG controller, which is already fully described by (3.12) and (3.14).

On

Stabilizing the Center of Mass with the support leg

In this section, we detail our approach to stabilize the Center of Mass (CoM) dynamics around a reference CoM trajectory. We adapt the existing admittance controller. Whole-body admittance control is used to stabilize the CoM around the trajectory during passive walking, as detailed in Section 2.3. However, during single support phases, the swing leg torques are entirely defined by the VG controller. Therefore, we propose below (in Section 3.2.1) an adaptation of the passive walking admittance SoT from Section 2.3 to compute the support leg joint targets taking into account the motion of the swing leg. In Section 3.2.4, we expose the rationale of our admittance SoT

STABILIZING THE CENTER OF MASS WITH THE SUPPORT LEG

73 modification and provide a numerical illustration of its performance for the admittance task.

Modification of admittance's SoT

As discussed earlier, the VG methodology implies that the swing leg target position q sw t follows a fixed path P . To prescribe the motion of the swing leg, we design a single support Stack-of-Tasks (SoT), denoted SoT 2 , with a task on the swing leg joints having higher priority than the CoM acceleration tracking task, as follows The swing leg task is to follow the swing leg target joint position. It writes e 0 0 = P (σ) -q sw (3.29)

We design the reference acceleration of this task s.t. it yields fast convergence to the same targets (P (σ), T (σ) σ) as the VG controller (3.8)

qsw * = K sw p (P (σ) -q sw t ) + K sw d (T (σ) σ -qsw t ) (3.30)
with (K sw p , K sw d ) high gains and (q sw t , qsw t ) the swing leg joint target computed during the previous control loop. This task is the one with the highest priority, therefore its associated cost w 0 0 = ||q sw * -qsw || 2 2 is always identically zero, yielding qsw t = qsw * (3.31) Thus, fast convergence of (q sw t , qsw t ) to (P (σ), T (σ) σ) is ensured (by the choice of high gains matrices (K sw p , K sw d )), which enables the IK solver to compensate for the bias introduced by the swing leg motion in the CoM acceleration task. This property is further derived in Section 3.2.4.

As in Section 2.3, the optimal acceleration of the support leg qsp t , resulting from tasks of level 1 -3, is integrated twice to define the target support leg position and velocity (q sp t , qsp t ). Finally, the support leg joint torques τ sp are computed using a high-gain PD controller

τ sp = K sp p (q sp t -q sp ) + K sp d ( qsp t -qsp ) (3.

Synchronization of legs by Time Rescaling (TR)

The SoT 2 described in the previous section incorporates the same CoM acceleration task as SoT 1 at level 1. The desired acceleration for this task is computed by (2.19) which, in turn, aims at tracking the desired CoP computed from (2.18). These two desired quantities are computed to stabilize the CoM around a state trajectory x, with corresponding input trajectory p. We now define these reference trajectories synchronized with the patient schedule σ.

A natural way to define two synchronized trajectories x and p from the nominal trajectory Q is to first time-rescale the nominal trajectory Q with the patient schedule σ, defining the time-rescaled trajectory

q T R q T R (σ) ≜ Q • s -1 (σ) (3.33)
and then, to compute the time-rescaled state trajectory x T R and corresponding time-rescaled input trajectory p T R from it

x = x T R (σ) ≜ FK • q T R (σ) p = p T R ≜ ID • q T R (σ) (3.34)
readily obtained using Forward Kinematics (FK) (see e.g. [Carpentier2019] for implementation) and Inverse Dynamics (ID). The kinematic reference trajectories required for the other tasks are readily computed from the time-rescaled trajectory q T R using FK. The derivation of the controller that we propose in this Chapter is now complete. It is schematically depicted on Figure 3.11.

Trajectory smoothing

Task changes produce discontinuities of target trajectories (q t , qt ). These occur either when lifting the foot off the ground (switching from SoT 1 to SoT 2 ), or when landing the foot (switching from SoT 2 to SoT 1 ).

STABILIZING THE CENTER OF MASS WITH THE SUPPORT LEG 75

Swing leg ctrl.

Support leg ctrl.

p T R (σ) ξ T R (σ) q T R (σ)

76CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING

The discontinuities at foot lift are smoothed when switching from SoT 1 to SoT 2 by interpolating between the nominal trajectory T sw and the last swing joint target position q sw t . This interpolation is performed using well-designed 5 th order Splines over the remaining duration of the step. The smoothing is performed once per step when the swing foot lifts off, hence the modified swing leg trajectory is entirely known at the beginning of the swing phase, and so the swing leg path P can be computed as explained in Section 3.1.2.

Besides, there are also target discontinuities when switching from SoT 2 to SoT 1 as, in the general case, the foot landing position in Cartesian space is not exactly on the trajectory due to CoM tracking errors. When the stabilization effect of admittance control is efficient, the CoM stays close to the nominal one. In this case, when tasks change occurs at impact, the swing foot lands almost at the expected place, and the targets are sufficiently close to make the discontinuity imperceptible. When the stabilization effect of admittance control is less efficient, the CoM can get noticeably far from the nominal one, which results in larger discontinuities due to the tasks change and the lack of target reconciliation strategy. We tackle the loss of balance from which they result in Chapter 4.

Impact of the modification of the SoT

We now expose the rationale of our modification of admittance control. For this purpose, we propose a formal analysis of the changes. We consider the exoskeleton in single support. It constitutes an open kinematic chain. We split this kinematic chain into a support sub-chain (the support leg), in contact with the ground, on which the admittance task bears, and a swing sub-chain (the swing leg), with no contact with the ground. We only study the impact of the SoT modification on the target CoM acceleration, assuming perfect tracking of the target joint acceleration q t . To lighten the notation, we do not distinguish target and actual positions for the rest of this section, i.e. q = q t . We advocate that, provided a sufficient number of DoFs are included in the support sub-chain, the admittance task (2.19) can be satisfied with performance similar to the original SoT 1 , exposed in Section 2.3. This is suggested by computations that we now detail.

Consider the schematic representation of the exoskeleton depicted on Figure 3.12 during the single support phase. The kinematic chain is opened, and it consists of 13 bodies B j linked together by 12 revolute articulations J i . Each body B j has an associated frame The subscript •| k is used for quantities defined w.r.t. the reference frame F k , expressed in the inertial frame W. Ω j is the angular velocity vector of frame F j w.r.t. the inertial frame W, and a k the acceleration of O k . The acceleration due to the swing leg motion is readily computed as m sp m sp + m sw csp (q sp , q sp ) + m sw m sp + m sw (ä 6 (q sp , q sp ) + Ω6 (q sp , qsp , q sp ) ∧ O 6 O c sp (q sp , q sw ))

csw = ä6 + Ω6 ∧ O 6 O c sw + csw | 6 + Ω 6 ∧ (Ω 6 ∧ O 6 O c sw ) + 2Ω 6 ∧ ċsw (3.

+

m sw m sp + m sw b(q sw , qsp , qsw , q sp , q sw ) ≜ f (q sp , qsp , q sp , q sw ) + m sw m sp + m sw b(q sw , qsp , qsw , q sp , q sw ) (3.37) with b(q sw , qsp , qsw , q sp , q sw ) = csw

| 6 + Ω 6 ∧ (Ω 6 ∧ O 6 O c sw ) + 2Ω 6 ∧ ċsw (3.38)
This term b, entirely determined by the known acceleration qsw of the swing sub-chain, can be seen as a bias in the equation (3.37) governing c.

As appears above, the first and second rows in equation (3.37), gathered into a single function f (q sp , qsp , q sp , q sw ), are controlled by the acceleration qsp of the 6 DoF of the support sub-chain. Provided a sufficient number of DoFs are used, the support sub-chain should be able to compensate for the bias b created by the swing sub-chain and generate the desired total CoM acceleration. Mathematically, if b belongs to the image of the mapping f , then it can be compensated for. The support leg possesses 6 DoF, while the CoM acceleration is a 3D quantity, therefore it seems possible to realize any given CoM acceleration.

To further validate this observation, the proposed split of DoF, and subsequent changes in the SoT, we perform comparative simulations at nominal velocity. The nominal curvilinear abscissa is our control variable for the swing leg σ(t) = s(t) (3.39) which yields readily computed using Forward Kinematics (FK) and Inverse Dynamics (ID), as in (3.34). Figure 3.13 depicts the error of the CoM acceleration task during one simulated step both for the original algorithm, using SoT 1 during the complete step (equivalent to the passive walking controller of Atalante described in Chapter 2), and for the proposed algorithm, using SoT 2 during single support phases and SoT 1 during double support phases. The error of the CoM acceleration task is the difference between the desired CoM acceleration cd , as defined by (2.19), and the target CoM acceleration, as defined by the target joint acceleration qt solution of the IK solution to SoT 1 or SoT 2 . As shown on Figure 3.13, the error magnitude of the CoM acceleration task during the single support phase (from 1.1 s to 2.05 s) is identical for the original algorithm and the proposed algorithm. In addition, Figure 3.13 includes leading and tailing double support phases (from 1.0 s to 1.1 s and from 2.05 s to 2.15 s), s.t. one can check that the error magnitude is not increased due to the sudden tasks changes at the single-double support phase changes. Thus, the comparison confirms that the proposed split has a sufficient number of actuation DoF to track the desired acceleration cd as well as the original admittance control scheme during single support phases.

q sw t = P • s(t) = T sw (t) qsw t = T • s(t)ṡ(t) = Ṫ sw (t) (

80CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING

Impact on safety of patient empowerment

In this section, we first study the stabilization performance of the proposed rehabilitation controller, both in simulation and during experiments. All the simulation results presented in this section have been obtained using the open-source simulator Jiminy [Duburcq2019]. Then, we provide a numerical illustration of the limits of the proposed synchronization strategy, which motivates the next Chapter.

Stability analysis

The safety of exoskeletons is directly linked to the stability of the walking exercises. Hence, the stability must be rigorously evaluated, as recommended in [Pinto-Fernandez2020]. To do so, we evaluate in this chapter the stability of the walking exercises enabled by this controller over 10-second simulations and experiments. We use a dummy instead of a patient during all experiments, as shown on Figure 3.14, because able-bodied users tend to involuntarily stabilize the exoskeleton using their upper bodies and/or legs (while disabled patients might not be able to act similarly).

In order to ensure the repeatability of the experiments, one must generate virtual patient efforts τ sw u and τ sp u , or the consequences of these efforts expressed through the variations of the patient velocity σu . Because Atalante does not possess any force sensor at the interfaces with the patient nor at the joint level, we do not have any experimental data to reproduce the patient efforts (τ sw u ,τ sp u ). Therefore, we generate virtual control variable velocity profiles σu to simulate the consequences of variations of the patient efforts, and not virtual efforts directly. We generate these profiles σu for arbitrarily tailored functions ηu as follows

σ u (t) = s • η u (t) (3.42)
We design the ηu functions with one non-nominal phase variable velocity change. These functions are parametrized with three variables: the nonnominal phase variable velocity η, the starting time of the velocity change t η , and the duration of the modification ∆τ . To generate a periodic signal, we make these functions start and end at nominal velocity. The η u phase variable function is created from a piecewise affine function, with a discontinuous phase variable velocity change, smoothed by a Gaussian filter. An example of an η u function is pictured on Figure 3.15. Note that these functions are not, a priori, exactly representative of patient efforts. However, they allow Figure 3.14: Dummy inside Atalante during passive walking (ensures that no user stabilizes the system).

for an intuitive and systematic exploration of a class of velocity profiles with various levels of difficulty. Note that, as we constrain the η u function to end at a specific final value η f , only a limited set of combinations of non-nominal velocity, duration and start time can be explored.

Simulation results

The stability performance evaluation of the controller proposed in this chapter consists of 10-second rigid-body walking simulations for each set of parameters. A simulation is deemed successful if the robot has not fallen at the end of the simulation. The aggregated results are shown on Figure 3.16 where, for each set of parameters (custom velocity duration, custom velocity), the rates of success are lumped in three bins: above 90%, between 35% and 90%, and below 35%. The custom velocity parameter η is varying along the Y axis of Figure 3.16, while the custom velocity duration parameter ∆τ is varying along the X axis. The heatmap on Figure 3.16 shows the stabilization algorithm has good performance only around the nominal velocity ηn ≜ 1.0, for custom velocities η between 0.8 and 1.5. With lower or higher velocities custom velocities η, the algorithm fails to stabilize the exoskeleton, especially when the velocity change is too long.

Surprisingly, the stabilization algorithm manages to stabilize some timerescaled trajectories which do not respect the constraints acting on the system. Indeed, while the nominal trajectory respects both the direct (the torque limits) and indirect (the unilateral force constraint) actuation limits, as detailed in Section 2.2, the time-rescaled nominal trajectory, defined by (3.33), violates the indirect actuation limit even for small variations around the nominal velocity. This is illustrated on Figure 3.17, where the rate of time-rescaled trajectories respecting the unilateral force constraint for each set of parameters (custom velocity duration ∆τ , custom velocity η) is reported. These rates, aggregated along the starting time of the velocity change t η dimension, are lumped in three bins: above 90%, between 35% and 90%, and below 35%. Comparing Figure 3.17 with Figure 3.16 shows that the stabilization algorithm manages to stabilize even some time-rescaled Figure 3.16: Heatmap of 10 s walking simulation success (good stabilization performance around the base scenario, poor ones in the low-velocity range). Success rate computed for each (custom velocity duration, custom velocity) couple, aggregating the (custom velocity start) dimension. trajectories which do not respect the unilateral force constraint. While this performance is surprising at first, they only reflect that the stabilization algorithm manages to dynamically find trajectories, close to the time-rescaled trajectories, which respect all the actuation constraints. In other words, these results reveal that there exist trajectories respecting all the actuation constraints in a neighborhood of trajectories that violate constraints, and the stabilization algorithm acts as a dynamic generation algorithm of these trajectories. The distance of these trajectories to the tracked trajectory, which does not respect the constraints, is reflected in the tracking error of the stabilization algorithm.

Experimental observations

The custom functions depicted in Section 3.3.1 are used again for experimental stability evaluation. As opposed to the simulations, the initial state of the exoskeleton-dummy system is standing still. Hence, a starting step, different from the next cyclic steps, is required. In order to get closer to the simulation setup, we manually stabilize this starting step and, if required, the two first cyclic steps by gently holding the exoskeleton side handles.

A reduced number of parameter combinations is tested to keep the total duration of the experiments at a reasonable level (i.e. being able to complete the experiments in one day). The custom velocity (t η , δτ, η) triplets are chosen to span as much as possible the simulations test range, but they are also adapted during the test sessions w.r.t. the first results. For instance, for each such triplet, we stop decreasing the custom velocity as soon as we find an unstable velocity profile. This decision is guided by the results depicted on Figure 3.16. For this reason, we only test velocities higher than ηmin = 0.6.

The stabilization results are shown on Figure 3.18. We could not try velocities higher than ηmax = 1.5 due to the electrical power limits of the exoskeleton. We see that the walking experiment is successfully stabilized for all the velocities in the [1.0, 1.5] velocity range, and for all velocity change durations. However, walking with our rehabilitation controller at velocities lower than the nominal one often leads to falling, even for short-duration changes. These results are mostly in accordance with the simulation results of Figure 3.16, but with lower performance in the low-velocity range.

The observed performance reduction from simulation to experiments can be explained by the model discrepancies with reality. Indeed, there are multiple effects unaccounted for in simulations, among which the most important one is the presence of flexibilities which result in several centimeters endeffector position errors [Vigne2021;Vigne2022]. Other factors which could also explain the performance reduction are, but not limited to, the unknown interaction of the dummy with the exoskeleton and model uncertainty (its total mass is precisely known, as with patients, but the limb masses and inertias are approximated), or sensor errors.

These experimental performances are not sufficient for the rehabilitation of deficient patients yet, as they will most likely have issues providing enough effort to go as fast as the nominal velocity. Reference quantities computed using the TR strategy and a simulated patient velocity σu as low as 60% (t η = 0.2 s, ∆η = 0.9 s and η = 0.6). Black horizontal lines represent the limits of the support polygon SP on the Y axis.

Performance analysis during a single step simulation

The shortcomings of our rehabilitation controller in terms of balance, especially in the low-velocity range, come from the choice of the Time Rescaling (TR) strategy described in Section 3.2.2. Indeed, this TR strategy does not take the unilateral contact constraints into account, as already exposed in Section 3.3.2. As a result, the CoP trajectory p T R is not confined to the support polygon SP, i.e. the state trajectory x T R does not respect the input-constrained dynamics, and the stabilization algorithm is unable to generate a CoP command which stabilizes the trajectory. This is illustrated on Figure 3.19 which reports the result of a simulation with a custom velocity η = 0.6 for δτ = 0.9 s. The reference CoP, in green, is not contained in the support polygon SP when modulating the trajectory at 60% of the nominal velocity. The final state (the endpoint of the red line) is different from the nominal final state x f ≜ FK • T (T f ) (the endpoint of the blue line). The 8 cm resulting error is sufficient to make the robot fall at the end of the step.

Chapter 4

Monitoring safety Chapitre 4 -Gestion de la sécurité: Dans ce chapitre, nous proposons une nouvelle stratégie de synchronisation des jambes et détaillons son implémentation dans le contrôleur développé au chapitre précédent. L'objectif de cette modification est d'augmenter l'équilibre lors des exercices de rééducations actifs avec Atalante. Cette nouvelle stratégie s'appuie sur un problème de contrôle optimal en temps final libre pour un modèle réduit du système qui est résolu en ligne à haute fréquence. Parce que résoudre un tel problème à une fréquence suffisamment élevée pour le confort du patient est difficile, nous avons développé une méthode de résolution originale réalisant une dichotomie sur la fonction de faisabilité d'un problème quadratique. La solution de ce problème fournit une nouvelle trajectoire pour le centre de masse du système, d'une durée cohérente avec la vitesse de la jambe de vol du patient, ou aussi proche que possible de cette vitesse, tout en respectant les contraintes dynamiques du système. Les modifications apportées au contrôleur sont donc doubles : la trajectoire de référence est remplacée par cette trajectoire optimale, et, lorsque nécessaire, la vitesse du patient est limitée pour préserver l'équilibre. Les performances de stabilisation de ce contrôleur sont évaluées en simulation et à l'aide d'expériences avec un mannequin. Des expériences illustrant le comportement général de ce contrôleur avec un patient sont également fournies.

The experiments of Section 3.3 have shown that the Time Rescaling (TR) strategy fails to generate trajectories that can be stabilized by admittance. In this chapter, we derive a novel Optimal Planning (OP) strategy and modifications to the controller from the previous chapter to solve this issue and increase performance.

The OP strategy leverages an Optimal Control Problem (OCP) to generate a CoM trajectory, meant to be tracked by admittance control, while taking the unilateral force constraint into account. The literature about generating CoM trajectories online is broad and dates back to the early 2000s, first as intricate preview regulators [Kagami2002; Kajita2003; Kagami2011], which require a CoP reference planner to provide the reference CoP to track, then as plain OCPs [Wieber2006b], which only require the contact sequence, and can even optimize their location [Herdt2010]. Since then, the literature about online walking trajectory generation has flourished, leveraging reduced [Feng2013;Caron2016] or complete [Erez2013; Koenemann2015; Dantec2022] dynamical models of the robots, to optimize first a reduced set of parameters (e.g. the horizontal motion of the CoM only [Wieber2006b]), and then gradually increase the number of free parameters (e.g. optimizing a combination the 2D CoM motion [Khadiv2016; Caron2020; Smaldone2021], CoM height [Caron2020], stepping location [Khadiv2016; Smaldone2021], and switching time [Khadiv2016;Smaldone2021]).

The OCP we formulate is specific for two reasons:

1. Our primary objective is to maximize patient expressiveness. Hence, our OCP has a free final time.

2. The optimal solution acts as a safety filter on the patient schedule.

Hence, it has to respect hard real-time constraints (typically running at least at 1kHz).

Among all the online CoM trajectory generation algorithms for biped robots published over the past 2 decades, we found only 5 explicit occurrences of algorithms optimizing the switching time (note that including quadrupeds adds only 2 occurrences). Namely, [Khadiv2016; Griffin2017; Smaldone2021; Ahn2021; Katayama2022], among which two demonstrated their algorithm in simulation only. In view of application, we did not consider these 2 algorithms with no experimental evidence of their performance. Among the three remaining algorithms, only one [Griffin2017] satisfies the hard 1kHz real-time constraint. However, this algorithm relies on the definition of the CoP trajectory as a constant piecewise function with specified durations of each piece, therefore the duration adjustment is not a free variable in their optimization, but rather comes as an input from a stabilizing feedback on the measured CoM state. Thus, this algorithm is not suitable for our use case, in which we require the OCP to maximize a cost formulated on the final time of the trajectory.

In details, we formulate an OCP as follows. We model the system as subjected to the LIP dynamics with the unilateral force constraint being the only actuation constraint. In addition, we limit the horizon to the end of the current step, and we fix the footprints to be the same footprints as in the nominal trajectory. The final time of the trajectory is free, and we set a penalty on its difference with some patient-dependent target time. This target time is computed by assuming the patient will produce efforts that keep the instantaneous swing leg velocity constant up until the end of the current step. To reflect that patient expressiveness has the highest priority over all other objectives, we formulate the OCP as a bi-level optimization problem, with the cost on the final time being the sole cost in the upper-level problem. The sublevel cost is only a control input regularizing cost. The importance of the sublevel comes from its constraints: the input-constrained LIP dynamics, plus some well-chosen final state equality constraint ensuring the recursive feasibility of the steps. Because of the sublevel constraints, the final time penalty can not always be nullified. Thus, our OP strategy reveals when the patient-driven swing leg velocity cannot safely be satisfied. In this case, the final time of the optimal solution provides us with the optimal safe velocity of the swing leg (i.e. which maximizes the patient desire by finding the closest final time to the patient target time).

The controller from Chapter 4 is adapted to this online CoM trajectory CHAPTER 4. MONITORING SAFETY generation strategy in two different ways. First, the admittance reference trajectory, used to control the support leg, is replaced by the solution to the OCP. Then, the patient swing leg velocity is replaced by the optimal velocity from the OCP. When the OCP manages to nullify the final time penalty, the optimal velocity is identical to the patient velocity and this change has no impact. When the final time penalty is not nullified, the patient velocity is filtered. Hence, the OCP also acts as a safety filter of the patient velocity. Therefore, it has to run at the same frequency as the VG controller. It is known from expert knowledge at Wandercraft that increasing the VG update frequency (3.14) directly increases the user's ability to drive the motion along the swing leg path P . It has been found that running VG at 1kHz yields a good trade-off between computation load and user efforts. Solving an OCP with free final time at such a high-frequency replanning is challenging, especially with the computational power limitations of the Atalante onboard computer. As we found no readily applicable solver in the literature, we have developed an original solving method of our bi-level OCP based on bisection on the feasibility function of Quadratic Program (QP)s. We show that this approach is optimal in almost all situations 1 encountered for the LIP dynamics in one dimension. We provide numerical hints that this is probably always true in cases of practical interest. Further, even in cases where our method does not find the optimal trajectory, i.e. when the upper-level cost is not minimized, the constraints of the sublevel are always satisfied. Therefore, even suboptimal trajectories are safe (i.e. respecting the dynamics and boundary conditions).

The chapter is organized as follows. First, we motivate and derive our OP strategy, along with the modifications it induces in the rehabilitation controller from Chapter 3. Then, we study the optimality of the bisection method used to determine the solution to the OCP. Finally, we quantify the stabilization performance of this controller with extensive simulations and during experiments with an able-bodied user.

Online Planning and velocity filter

In this section, we first describe the low-level fixed final time OCP with the high-frequency requirement in mind. Then, we compute a target final time from the patient-driven swing leg velocity, and we formulate the complete bi-level OCP with a cost on the final time at the upper level. Finally, we 1 i.e. for a large class of initial conditions and target final times

ONLINE PLANNING AND VELOCITY FILTER

93 describe how we modify the controller from the previous chapter for its use with the bi-level OCP.

Trajectory generation with fixed final time

The target frequency of the OCP solving is the same as the VG controller frequency, i.e. 1 kHz. For this reason, we use a reduced model of the system dynamics, namely the LIP dynamical model (2.17). Recall that it writes, for each dimension i ∈ {x, y} of the horizontal plane,

χi = Aχ i + Bu i (4.1)
with the one dimensional state χ i ≜ (c i , ċi ), and the CoP command u i = p i . This yields, in two dimensions,

χ = Aχ + Bu (4.2) with χ ≜ (χ x , χ y ), A = A O 2x2 O 2x2 A , B = B O 2x1 O 2x1 B , and u ≜ u x u y ∈ SP.
For some duration T ∈ R + , we define the set of admissible controls U ad (T )

U ad (T ) ≜ {u s.t. ∀t ∈ [0, T ], u(t) ∈ SP(t)} (4.3)
where SP is the support polygon and, from some initial state χ 0 and some input u ∈ U ad (T ), we denote χ u the forward integration of (4.2) from χ 0 , s.t.

∀t ∈ [0, T ], χ u (t) = e At χ 0 + t 0 e A(t-τ ) Bu(τ )dτ (4.4)
To get to a simple yet insightful OCP, we only optimize for the state trajectory and control input up until the end of the current step. For this reason, the support polygon SP, which varies with each contact phase, is constant over the optimization domain considered.

The VG controller restricts the freedom of the patient in such a way that the path of the patient-exoskeleton system's swing foot, in the back referential, is identical to the nominal path. Hence, with this framework, the swing foot location in the world frame W could theoretically be left free for the optimization to choose, by changing the CoM final location 2 . However, the landing location of the swing foot is constrained by the kinematic limits of the exoskeleton joints, but this constraint can not be incorporated in LIPbased optimization and would require a non-trivial approximate surrogate. For this reason, we do not include the landing location of the swing foot in the optimization variables. Instead, we consider it fixed, at the location defined by the nominal trajectory T , which respects the kinematic limits of the exoskeleton according to Section 2.2.

It is critical, for our application, to guarantee the feasibility of the next step i.e. the recursive feasibility of the OCP. Indeed, if not enforced, it would jeopardize the safety of the walking exercise, hence the safety of the patient. Classically, the recursive feasibility of LIP-based OCPs is either guaranteed with a terminal constraint on the state, or only softly guaranteed with an integral cost on some diverging quantities3 , guiding the planner toward recursively feasible trajectories [Wieber2016]. The soft guarantee is not enough in our case for safety reasons, hence we follow the terminal constraint approach. One type of such constraint is capturability constraints, i.e. a constraint ensuring that the robot will be driven to a complete rest in N steps. There is strong numerical evidence [Ciocca2017] that including a capturability constraint does not jeopardize the recursive feasibility, even when this constraint is recursively postponed into the future because the objective is to keep walking, and not to stop. However, this is not guaranteed, and we prefer not to take any risk for the safety of the patient. Thus, we include a different constraint, ensuring strict recursive feasibility of our OCP.

Designing a constraint that ensures the recursive feasibility of LIP-based OCPs requires some knowledge or assumptions about the desired motion after the optimization horizon, as studied in [Scianca2020]. In their work, they do not possess any offline trajectory to inform them about the desired plan after the horizon, hence they study several ways to estimate this plan. While they show that they manage to ensure the recursive feasibility of their OCP, their solution entails the design of a sophisticated terminal constraint. Luckily, we have an offline trajectory at our disposal to guide us, enabling us to design a simpler terminal constraint.

Assuming the desired plan is to keep walking, we design our terminal constraint as an equality constraint on the final state of the trajectory, with the nominal final state as value, i.e. the final state of the nominal trajectory. Indeed, the nominal trajectory satisfies the input-constrained full dynamics from its initial state to its final state in a certain time T nom . If there exist also a trajectory satisfying the input-constrained LIP dynamics joining the same initial and final conditions over the same duration T nom , then we can use the nominal final state as a constraint for our OCP. This LIP feasibility check is performed numerically by solving a QP problem at the beginning of each exercise. The nominal walking trajectory being cyclic, it is sufficient to check it over a single step.

Finally, our terminal constraint writes

χ u (T ) = χ f ≜ FK • T (T nom ) (4.5)
with T nom ∈ R + the duration of the nominal trajectory T . We denote Ω the set of feasible commands respecting the initial and final conditions (χ 0 , χ f )

Ω(χ 0 , χ f , T ) ≜ u ∈ U ad (T ), χ u (0) = χ 0 , χ u (T ) = χ f (4.6)
and the fixed final time OCP writes4

Problem 4.1. Given χ 0 , χ f and T , find u opt s.t.

u opt = arg min u∈Ω(χ 0 ,χ f ,T ) T 0 u 2 dt
This OCP constitutes the lower level of our bi-level strategy.

Bi-level trajectory generation with free final time

The patient's efforts directly affect the instantaneous swing leg velocity through (3.14). At each time, the swing leg velocity equivalently defines a forecasted step duration T as follows

T = L max -σ σ (4.7)
Imposing (4.7) as the final time of Problem 4.1 would synchronize the support and swing leg trajectories. However, Problem 4.1 may be unfeasible for some values of T . Therefore, we search for the closest step duration for which Problem 4.1 has a solution.

The patient desired velocity is thus only a target velocity from now on and denoted σt . It is directly computed from (3.14), saturated to be strictly positive, s.t.

σt = max   ϵ, T (σ) ⊤ K sw p (q sw -P (σ)) + K sw d qsw T (σ) T K sw d T (σ)   (4.8)
with ϵ ∈ R + * a small parameter. The target velocity is converted into a target time T t until the end of the current step as follows

T t = L max -σ σt (4.9)
Problem 4.1 has solutions only for durations in T(χ 0 , χ f ) ⊂ R + , the set of feasible times for which Ω is not empty

T(χ 0 , χ f ) ≜ {T > 0, Ω(χ 0 , χ f , T ) ̸ = ∅} (4.10)
However, there is no guarantee that a given T t is in T(χ 0 , χ f ). For this reason, we include the target time T t in the upper-level cost of a bi-level OCP as follows Problem 4.2. Given χ 0 , χ f and T t , find u opt and T opt as

T opt = arg min T ∈T(χ 0 ,χ f ) ||T -T t || 2 s.t. u opt = arg min u∈Ω(χ 0 ,χ f ,T ) T 0 u 2 dt
Problem 4.2 is a minimal time problem for input-constrained linear dynamics of dimension 4. The solutions T opt and u opt are determined using a bisection method on the final time T , granting high-numerical efficiency. We implement this algorithm at 1 kHz.

Bisection resolution on the feasibility: In details, at time t k , given an initial condition χ 0 , a target time T t k > 0, and an initial guess T 0 ∈ T(χ 0 , χ f ), we solve Problem 4.2 using bisection between T t k and T 0 on the feasibility function5 of Problem 4.1 (the lower-level of Problem 4.2). Classically, the search interval is reduced by a factor 2 N using bisection, where N is the maximum number of iterations (typically 10). We obtain the solution u opt k , defined over T opt k ∈ T(χ 0 , χ f ). Then, we update the initial condition χ 0 and initial guess T 0 using the optimal solution (u opt k , T opt k ) as follows

χ 0 = χ u opt k (dt) , T 0 = T opt k -dt (4.11)
where dt = t k+1 -t k > 0 is the duration between two control loop iterations. We repeat this procedure until T 0 ≤ 0.

ONLINE PLANNING AND VELOCITY FILTER

97

For a given T ∈ T(χ 0 , χ f ), determining u opt solution of Problem 4.1 is a fixed horizon input constrained LTI problem, which can be readily solved numerically because it is convex. A more challenging point is the description of the set T constraining Problem 4.2. It is the subject of Section 4.2, where we leverage our main result to characterize the nature of T and conclude on the optimality of the solution found using bisection.

Before detailing in Section 4.2 how we solve Problem 4.2, we describe how we use T opt , to filter the patient target velocity σt , and u opt , to define the support leg targets.

Modification of the rehabilitation controller

Two modifications of the rehabilitation described in Chapter 3 are necessary for using the solution of Problem 4.2.

First, the reference trajectory (ξ * , p * ) tracked by the admittance-based stabilizer, previously computed by the time rescaling of the nominal trajectory T , must be replaced by the optimal reference trajectory (ξ opt , p OP ) s.t.

ξ * = ξ opt , p * = p OP (4.12)
with p OP = u opt and ξ opt readily computed from x OP = χ u opt . Then, the swing leg schedule must be replaced by the optimal schedule σ opt computed from the optimal time T opt as follows σopt ≜ L max -σ opt (τ ) T opt (τ ) (4.13) with σ opt recursively computed from the optimal velocity σopt by numerically integration. The derivation of the controller that we propose in this Chapter is now complete. It is schematically depicted on Figure 4.1, with the Optimal Planning (OP) block detailed in Algorithm 1. This OP strategy induces two different changes. First, the reference trajectory (ξ * , p * ) is guaranteed to respect the input-constrained LIP dynamics. This is illustrated on Figure 4.2a, where, as opposed to Figure 3.19, the CoP reference trajectory p OP y is entirely contained in the support polygon (black horizontal lines), and the CoM trajectory c opt y is offsetted toward the y = 0.0 line. As a result, the endpoint of the forward integration of the full-state dynamics (2.10) (the endpoint of the red line) is close to the nominal final state (the blue cross): the state is successfully driven to the final state using the stabilization controlled Section 3.2. This recursively ensures the success of the walk. Second, the simulated target velocity σt is filtered by the OP Algorithm 1 Optimal Planning (OP)

Require: χ 0 Require: T 0 ∈ T(χ 0 , χ f ) while T 0 > 0 do

T t k ← Lmax-σ opt k-1 σt k ▷ Convert target velocity to target time T opt k = arg min T ∈T(χ 0 ,χ f ) ||T -T t k || 2 ▷ Bisection on QP feasibility function s.t. u opt k = arg min u∈Ω(χ 0 ,χ f ,T ) T 0 u 2 dt ▷ Solved as a QP σopt k ← Lmax-σ opt k-1
T opt k ▷ Convert optimal time to optimal velocity

σ opt k ← σ opt k-1 + σopt k dt ▷ Integrate optimal velocity over dt χ 0 ← χ u opt k (dt) ▷ Update current state T 0 ← T opt k -dt ▷ Update current step duration output (σ opt k , u opt k (dt), χ 0 )
▷ Output toward VG and admittance end while strategy. This is illustrated on Figure 4.2b, where the optimal velocity σopt is identical to the target velocity σt until 1.1 s, i.e. the target time T t is in the set of feasible times T(χ 0 , χ f ). After 1.1 s, the optimal velocity σopt is different from σt , and the optimization finds the duration T opt ∈ T(χ 0 , χ f ) closest to T t . In the next section, we show our bisection almost always finds the optimal feasible final time.

Bisection on trajectory final time

The support polygon SP is constant over the optimization horizon and equal to the footprint of the support foot, which is assumed to be rectangular, s.t.

u ∈ SP ≜ [u mx , u M x ] × [u my , u M y ] (4.14)
with u mx < u M x and u my < u M y . Since the local frame is defined by the support foot, this yields that the x and y dynamics of the LIP model (2.17) are decoupled. Therefore, with the trajectory duration T fixed, the set of feasible commands is also decoupled, i.e. 4.2: Reference quantities computed using the OP strategy and a simulated target velocity. The support polygon limits are not violated using the OP strategy, the CoM reference is successfully tracked, and the simulated target velocity is followed as much as possible.

Ω(χ 0 , χ f , T ) = Ω x (χ 0 x , χ f x , T ) × Ω y (χ 0 y , χ f y , T ) (4

BISECTION ON TRAJECTORY FINAL TIME

101

with

Ω i (χ 0 i , χ f i , T ) ≜ u i ∈ U adi (T ), χ u i i (0) = χ 0 i , χ u i i (T ) = χ f i (4.16) and U adi (T ) ≜ {u i s.t. ∀t ∈ [0, T ], u i (t) ∈ [u mi , u M i ]} (4.17) for i ∈ {x, y}.
Hence, the sole link between the sets of feasible commands in x and y dimensions is the trajectory duration T , which has to belong to both sets of feasible times for trajectories to exist with this duration, i.e.

T(χ 0 , χ f ) = T x (χ 0 x , χ f x ) ∩ T y (χ 0 y , χ f y ) (4.18)
In the next section, we characterize the set T.

Set of admissible final times

In this section, we characterize the sets T x and T y . Since the analysis is independent of the direction, we use the notation χ (resp. u ) to denote χ x or χ y (resp. u x or u y ).

Preliminary results on minimum and maximum time solutions Lemma 1. For all (χ 0 , χ f ) ∈ R 4 and T > 0, if there exists a solution u ∈ Ω(χ 0 , χ f , T ), then a minimum time solution (noted u min ) always exists and, when the set T is upper-bounded, a maximum time solution (noted u max ) exists. They are global optima.

Proof. Equation (A.2) is linear, and U ad is compact and convex, hence, when a solution u ∈ Ω(χ 0 , χ f , T ) exists, a minimum time solution u min exists from [Liberzon2012, Theorem 4.3].

When the set T is upper-bounded, we note T its supremum. Given a sequence (T k , u k ) s.t. lim k→∞ T k = T , consider the sequence (T , ũ k ) of prolonged u k on [T k , T ] by the null function, then the proof provided in [Liberzon2012, Theorem 4.3] is straightforwardly extended to the ũ k sequence, yielding the existence of u max . Hence, T is maximum.

We denote

(T min , u min ) ≜ arg min T, u ∈ Ω(χ 0 , χ f , T > 0) (T max , u max ) ≜ arg min -T u ∈ Ω(χ 0 , χ f , T > 0) (4.19)
The following Theorem 1 fully characterizes the sets T x,y and is instrumental to efficiently solve the bi-level OCP. Proof. The proof of Theorem 1 is provided in Appendix A.

Theorem 1 reveals that T is composed of one or two intervals in each of x and y dimensions. This result also holds for the intersection of the two sets, hence it holds for T = T x ∩ T y . This is instrumental for numerically solving Problem 4.2. Knowing that T is composed of one or two intervals, the solution is simply the projection of T t onto them. As detailed in Section 4.2.2, numerical evidence hints toward T being composed of a single interval in our practical case, therefore the projection is readily obtained by a bisection method applied to the feasibility function of a quadratic program.

Numerical method

Theorem 1 describes T. We notice, numerically, that T is most often a single interval for the sets of initial conditions and target time typically encountered in walking rehabilitation, which enables us to use bisection to solve Problem 4.2. To compute T, we use the nominal terminal condition χ f , which is the actual terminal condition set in Problem 4.2, and a time-varying initial condition χ t , computed from the nominal trajectory T as follows

χ t = FK • T (t) (4.20)
In order to ease the visualization of T, we rather compute the set of feasible velocities V defined according to (4.9) as follows 4.4b shows that V y (χ t y , χ f y ) is an interval for all initial conditions χ t (varying with time). However, Figure 4.4b reveals that V x (χ t x , χ f x ) is the union of two intervals for initial conditions χ t evaluated at times near the end of the nominal trajectory duration. This numerically illustrates our main Theorem 1, as it states that T can be the union of two intervals.

V(χ t , χ f ) = L max -s(t) T , T ∈ T(χ t , χ f ) (4.21)
Nonetheless, Figure 4.5, which depicts the intersection of V x and V y , i.e. V(χ t , χ 0 ) according to (4.18), reveals that V is actually always a single interval along the nominal CoM trajectory. This does not guarantee that V is always a single interval along all possible CoM trajectories, but of all the numerical exhaustive searches that we have performed along various CoM trajectories (recursively defined through the optimal solutions of Problem 4.2 with various target velocity profiles ηt ), none revealed V being the union of two intervals. However, if, for any set of initial and final condition (χ 0 , χ f ), the set of feasible times T becomes the union of two intervals, then only the optimality of the upper-level of Problem 4.2 would be lost (i.e. the bisection method would return a feasible duration T ̸ = T opt ). The solution u opt , by-product of the bisection method, would still minimize Problem 4.1, and therefore, respect the input-constrained LIP dynamics and the two boundary conditions (χ 0 , χ f ). Thus, the trajectory would still be safe for the robot to stabilize around, and only the optimality of the rehabilitation objective would be lost.

QP resolution and feasibility check

Problem 4.1 defines a fixed-time OCP that can be addressed using a direct numerical method. Conveniently, the input signal is represented by a piecewise C 1 function in between non-uniform nodes. The dynamics and the value of the integral cost are exactly represented using the first-order hold quadrature. This allows the expression of boundary conditions and input constraints under an affine form in a finite number of variables, and the cost as a quadratic function of these variables. The same discretization procedure is employed in the x -y directions, resulting in a QP with 4P + 4 variables and 4P + 4 affine constraints. The outcome of the QP resolution is a feasibility boolean, and, when it is feasible, a solution to Problem 4.1.

Towards real-time implementation

Typical numerical setups consider 2P + 2 = 10 variables, with N = 10 maximum number of iterations. The employed software is a streamlined implementation of the positive definite QP dual algorithm from [Goldfarb1983] specifically coded in C++ for this application to minimize any overheads. The problem is treated as dense. Mean CPU times on a laptop with turboboost disabled (s.t. the CPU frequency is similar to the embedded computer of Atalante) is 0.2 ms, varying between 0.07 ms and 0.75 ms, which are lower than the 10 ms reported in [Caron2017] and the 100 ms reported in [Pon-ton2021] where similar online planning problems are addressed. They are also lower or equal to those reported in [Fernbach2020;Caron2016] where fixed-time online planning problems are solved. These computation times are in line with the target 1 kHz frequency of the OCP solving.

Simulation results

In this section, we first illustrate the performance of our bisection search of the optimal trajectory duration with a highly varying simulated patient input. Then, we study the increase of stability offered by the use of the OP strategy instead of the TR strategy using extensive full-body simulations of the patient-exoskeleton system under a variety of simulated patient inputs.

Results for highly varying target velocity profile

We first perform a simulation with a strongly varying target velocity profile σt . Similarly to Section 3.3.1, we generate this target velocity profile σt from an arbitrarily tailored function η t as follows

σ t (t) = s • η t (t) (4.22)
A strongly oscillating velocity ηt along the geometric path is considered to simulate the behavior of a highly demanding patient. The nominal velocity is 1 and the variations are ±50%. This defines a signal t → T t (t) from (4.9) and the curvilinear abscissa s. For reference, the set of feasible velocities V(χ u opt , χ f ) is estimated by exhaustive numerical search at each instant. Figure 4.6 depicts the corresponding target velocity ηt and optimal velocity ηopt ≜ d dt s -1 • σ opt . As is visible in Figure 4.6, the resolution of Problem 4.2 satisfies the simulated patient's intent at the beginning of the simulation. Gradually the feasible set gets more stringent and at some point, near t = 0.4 s, the proposed algorithm has to intervene. The desired time T t is no longer feasible on many occasions. The situation worsens until the end of the simulation. Notably, in the end, the walking velocity has to be sped up significantly. This comes from the numerical discretization of the horizon of QP problem solving Problem 4.1.

Results on full-body simulations

Extensive closed-loop rigid-body simulations of the patient-exoskeleton system have been performed to evaluate the stability increase of the walking exercises offered by the use of the OP strategy (in particular w.r.t. the TR strategy). To simulate the behavior of the patient, we consider the same smoothed piecewise velocity signals as the ones described in Section 3.3.1. Figure 4.7 reports the results (for each duration ∆τ and magnitude ηt , a vast list of possible starting times t η for the square disturbance is considered, and we report the success rate). Figure 4.7 (left) reports simulation results obtained with the controller derived in Chapter 3, hence using the TR strategy, while Figure 4.7 (right) reports the results obtained using the same controller with the modifications presented in this chapter, hence using the OP strategy. Again, a simulation is considered stable if the simulated patient-exoskeleton system walks for at least 10 s without falling.

These results show the substantial improvement in the stability of the walking exercises provided by the use of OP strategy over the TR one. In details, a total of 2917 simulations have been conducted. In total, less than 8% of considered velocity profiles are failing with the OP strategy, while more than 30% were with the TR strategy. Therefore, safety is ensured for almost all considered cases, except for some very low-velocity cases with long durations. A careful examination of simulation results reveals that fall occurs mostly when the slow-down takes place at the late stages of the step.

Experimental results

Two types of experiments are conducted. First, we compare the stability performance of the two strategies and conclude on the vast superiority of the OP strategy. Then, we explore the performance of the OP strategy in terms of compliance with the patient schedule.

Stability comparison

We perform extensive stability comparisons as in Section 3.3.3 by replacing the patient with a dummy and simulating the behavior of the patient as described in Section 3.3.2. An experiment consists of walking 10 steps in a straight line. A practitioner keeps hold of the two lateral exoskeleton handles and is allowed to create an effort with one finger on each hand only. This creates a very low upper-bound on the external forces 8 .

The reported results on Figure 4.8 show a great safety improvement offered by the OP strategy in the low-velocity range, below 90% of nominal velocity, compared to the TR strategy. They stress that, using the OP strategy, the proposed controller is completely preserving the balance of the system for velocities about as low as 70% of the nominal velocity, and velocities as low as 50% of the nominal velocity provided the change duration is lower or equal to 300 ms. These results are in complete alignment with the stability results obtained in simulation and reported in Section 4.2.3.

In addition, the mean computation time during these experiments is 0.039 ms, and varies between 0.22 ms and 0.11 ms, which validates the efficiency of our bisection-based solver (running on a i7-1185G7E CPU at fixed 1.8 GHz frequency).

8 Stabilization of patients walking with Atalante is still requiring some very low external efforts, probably due to the inaccuracy of the patient model and/or their loose fastening to the exoskeleton for safety reasons. Able-bodied users doing walking rehabilitation can walk fully autonomously by moving their upper body to help maintain balance. The proposed setup somewhat emulates this. 9 The white spaces in this figure corresponds to unfeasible values of the parameters violating the constraint σ * ≤ Lmax.

Rehabilitation: experiments with an able-bodied user

We report below the results of two walking experiments with an able-bodied user using the proposed controller with the OP strategy:

• a 10-steps walking experiment, during which the user drives the walking velocity between 25% and 150% of the nominal walking velocity;

• a single-step walking experiment, during which the user brings the walking velocity to 0 for 2 seconds, then resumes walking.

10-steps walking experiment

We first report a single step velocity (on Figure 4.9) and CoM trajectory (on Figure 4.10). Figure 4.9 reveals how the OP strategy accounts for the unilateral contact constraint by projecting the target velocity on the set of feasible velocities. In detail, during the first 450ms of the single support phase, the OP strategy leaves the target velocity unchanged because the solution of Problem 4.2 is T opt = T t ( σopt in green completely overlaps σt in blue). Hence, the patient's desire is fully respected. During the remaining 500 ms, the OP strategy starts filtering the target velocity to preserve the balance of the system, σ opt ̸ = σt . Gradually, σ opt is constrained around 71% normalized velocity.

The OP strategy also wisely chooses the CoM reference trajectory and satisfies the final state constraint. Figure 4.10 shows the final state constraint is satisfied by both the replanned (OP) and time-rescaled (TR) CoM trajectories' endpoints. However, the CoM trajectory computed with the OP strategy is very different from the one with TR strategy, and, in particular, exhibits a minimum on the Y axis 2.5 cm closer to the support foot (centered at 0.0 cm) than the nominal trajectory: the exoskeleton sways its hip toward the support foot to accommodate for the user's low-velocity desire while remaining stable with the OP strategy, while the TR strategy only accommodates for the user desire without taking care of the balance of the system. The OP strategy's behavior is consistent with human behavior.

Finally, Figure 4.11 shows σopt over the whole experiment, with double support areas in grey. During this experiment, the user varies the level of effort produced by their legs during the single support phases. Note that (3.14) is only used during these phases while the user's desire is ignored during double support phases10 .

Mid-step complete stop experiment

We report a single step velocity (on Figure 4.12) and CoM trajectory (on Figure 4.13). Figure 4.12 reveals how the OP strategy accounts for the unilateral contact constraint by letting the user drive the walking velocity to almost 0% of the nominal walking velocity for approximately 2 seconds, when possible, then projects the target velocity on the set of feasible velocities by limiting the target slow-down during the late part of the step to roughly 40% of the nominal walking velocity.

The CoM trajectory computed with OP strategy, reported on Figure 4.13, differs even more from the one with TR strategy than reported previously. In particular, it exhibits a minimum on the Y axis 5 cm closer to the support foot (centered at 0.0 cm) than the nominal trajectory: the exoskeleton sways its hip toward the support foot twice more than during the previous experiment to accommodate for the user desire to stop complete walking in the middle of the step. This illustrative experiment reveals the great adaptability to the user desire offered by the OP strategy. The proposed control methodology updates the existing controller of Atalante in two ways: on one hand, the walking velocity is relaxed according to the patient's willingness and capacity to participate in the motion, enabling walking rehabilitation with Atalante. On the other hand, the patient schedule is monitored and followed at best while respecting the dynamics of the system, guaranteeing the safety of the walking exercises. The proposed controller is the first safe walking rehabilitation controller designed for self-balanced exoskeletons and the main contribution of the thesis.

Thesis outcome

Using the presented empower and monitor control methodology, we have been able to modulate by ±50% the temporal execution of two straight-line walking trajectories of practical interest: the simple flat-foot trajectory, at an average translational velocity of 15 cm/s, and the more anthropomorphic foot-rolling trajectory, at an average translational velocity of 40 cm/s. Videos obtained with an able-bodied user are given below: 116 CHAPTER 5. DISCUSSION AND PERSPECTIVES

• Experiments with the flat-foot trajectory: https://youtu.be/_1A-2 nLy5ZE.

• Experiments with the foot-rolling trajectory: https://youtu.be/ZKO ouUXfTvw.

Limitations and possible extensions

We have explored the possibility of offering temporal freedom to the user with the self-balanced exoskeleton Atalante. The obtained results are encouraging. However, the clinical evaluation of the proposed control methodology remains to be done. Among the various metrics of importance, one should consider, without loss of generality, the level of restored functional ambulation ability (for instance evaluated with the Functional Ambulation Categories test), the fatigue level of physiotherapists, and the sense of safety of patients.

On the methodological side, many ways of improvement shall be considered in the future. It is possible to combine temporal freedom with spatial freedom, which has already been considered in the literature. For this purpose, one could introduce constraint relaxation in the spirit of [Gurriet2020], and/or formulate a more general OCP allowing one to include additional variables such as the location of the swing-foot landing. Additionally, a mathematical formulation could also consider the reverse execution of the trajectories as required in some repetitive physical training, or reproduce the behavioral description of the two-thirds power law [Ivanenko2002] to better fit human intent.

Besides, for improved stabilization properties, it could be of interest to consider other tasks than the considered admittance task, like foot damping control (see [Caron2019]), or reduction of the disturbances coming from the flexibilities, see [Vigne2021].

Finally, the class of trajectories under consideration could be extended. Two extensions of dominant interest are turn-in-place and slope/stairway trajectories. Interestingly, these two use cases challenge the representativeness of the LIP model for mobile/changing reference frames along trajectories and some adjustments are necessary. Two distinct paths can be followed to improve the dynamical model in Problem 4.2 and extend the proposed methodology: leverage the recent progress of online whole-body MPC (see [Dantec2022;Mastalli2022]), or compute an efficient oracle offline using machine learning (see [Castillo2021;Siekmann2021;Singh2022]). defining a control u as detailed above, over τ ∈ [0, t f ] we note the solution x seq ≜ x u of the differential equation ẋ = Ax + Bu which is ϕ(x 0 , seq, τ ) ≜ x seq (τ ) = e Aτ x 0 + τ 0 e A(τ -s) Bu(s)ds By extension, we define ϕ(x f , seq, -τ ) ≜ e -Aτ x fτ 0 e -A(τ -s) Bu(T -s)ds

We define several subsets of R 2 as follows

D ≜ {(x 1 , x 2 ), x 2 = -x 1 } D + ≜ {(x 1 , x 2 ), x 2 > -x 1 } , D -≜ {(x 1 , x 2 ), x 2 < -x 1 } and U - m ≜ {(x 1 , x 2 ), x 1 < u m } , U + m ≜ {(x 1 , x 2 ), x 1 > u m } U - M ≜ {(x 1 , x 2 ), x 1 < u M } , U + M ≜ {(x 1 , x 2 ), x 1 > u M } Finally,

we define two open double cones

C M ≜ {D + ∩ U - M } ∪ {D -∩ U + M } , C m ≜ {D -∩ U + m } ∪ {D + ∩ U - m }
Next, we study the solutions of minimal and maximal time OCPs (4.19). This study stresses the role of several regions in the phase plane being key in the reachability of a target x f from an initial condition x 0 .

A.2 Regions of interest in the phase portrait Lemma 2. The solution u min is bang-bang, i.e. takes only values in {u m , u M }, with a maximum number of one switch. The same property holds for u max when it exists.

• x 0 ∈ R 258 • x 0 ∈ R 2 \R 258 : Bd x 0 ∈ R 5 : ∞ x 0 ∈ R 2 (x 0 ∈ R 8 ) x f ∈ R 123 : Bd (x f ∈ R 789 )
x f ∈ R 456 : ∞ Proof. We split the proof according to the location of x 0 in the phase plane and, when necessary, the location of x f . Only cases corresponding to nonempty T are considered.

For x 0 ∈ R 2 \{R 258 }, the argument stems from the monotonic divergence of x 1 . For instance x 0 ∈ R 147 , there exists ϵ > 0, s.t. x u 1 (0) ≤ u m -ϵ. Then, using (A.4), one easily shows that ∀t > t 0 , ẋ1 = ω(x 1 -u) ≤ -ωϵ. Therefore, the final time is upper bounded by x 1 (t 0 )-x f 1 ωϵ . A similar inequality is obtained for x u 1 (t 0 ) ≥ u M + ϵ to cover R 369 . Hence, T is upper bounded. For x 0 ∈ R 5 \D, which is entirely covered by C m ∪ C M , and is stable by symmetry w.r.t. D. Lemma 3 permits to build a sequence that periodically returns to x 0 , prolonging infinitely any solution from x 0 . Hence, T is not upper-bounded.

For x 0 ∈ R 5 ∩ D, for all possible values of u, the tangent vector field at x 0 is orthogonal to D and does not vanish. For any short time, the preceding rationale applies. For x 0 ∈ R 2 and x f ∈ R 123 , one has ẋ2 < 0, therefore x 2 is decreasing, hence ẋ2 ≤ -ω(x f 2 +u m ). Therefore, T is upper-bounded by x 0 2 -x f 2 ω(x f 2 +um) < ∞. For x 0 ∈ R 2 and x f ∈ R 456 , there exists a sequence from any x 0 s.t., for some t > 0, x w ≜ ϕ(x 0 , seq, t) ∈ R 5 . In addition, any x f can be accessed from this waypoint x w through a sequence (u m , a, b) or (u M , a, b), with a, b > 0. Therefore, a transient from x 0 to x f passing through x w can be arbitrarily prolonged with sequences periodically returning to x w . Hence, T is not upper-bounded.

The case x 0 ∈ R 8 the analysis is similar to x 0 ∈ R 2 . This completes the proof. Proof. Lemma 4 shows that for T to be bounded either x 0 ∈ R 2 \R 258 , or (x 0 , x f ) ∈ R 2 × R 123 , or (x 0 , x f ) ∈ R 8 × R 789 .

As we only consider the case of bounded T in this section, Lemma 1 shows the existence of solutions of (4.19). In general, there are at most two bang-bang sequences with one switch between x 0 and x f which are noted seq m ≜ (u m , a m , b M ) and seq M ≜ (u M , a M , b m ). Further, according to Lemma 2, the two controls u min and u max are bang-bang with at most one switch. Hence, either u min = seq m and u max = seq M , or the other way around.

By definition, T ⊂ [T min , T max ]. When T min = T max , T is a singleton, hence is trivially convex. We now assume T min < T max . The rest of the proof depends on the location of (x 0 , x f ) relative to D. (x 0 , x f ) ∈ D -2 (on the same side) The case is identical to the previous case, x 0 belonging to R 9 ∩ D -and x f belonging to R 9 ∩ D + ∪ R 6 . This completes the proof.

A.5 Convexity of unbounded T cases

Let us define J + a subset of R 4 as follows

J + ≜                    x 0 ∈ R 25 ∩ D + , x f ∈ R 56 ∩ D + s.t. ∃x D , x d ∈ D ∩ R 5 , t D > 0, t d > 0,       
x d1 > x D1

x f = ϕ(x D , u m , t D )

x 0 = ϕ(x d , u M , -t d )
and J -by symmetry of J + w.r.t. D. The set J + is partially pictured in The torque vector τ for all the joints of the exoskeleton is now computed as a high-gain PD controller plus gravity compensation τ (σ, σ) = τ P D (σ, σ, q j , qj ) + τ GC (q) (B.3) with τ P D (σ, σ, q j , qj ) = K p (P (σ) -q) + K d (T (σ) σqj ) (B.4)

where σ ∈ [0, L max ] is still the control variable and q j ∈ R 12 is the vector of joint positions. Following the VG approach, the control variable velocity reads σ(σ, q j , qj ) ≜ T (σ) ⊤ K p (q j -P (σ)) + K d qj T (σ) T K d T (σ) (B.5)

The control variable σ defines the patient schedule by the relation

η : t ∈ [0, t f ] → s -1 • σ(t) ∈ [0, T f ] (B.

B.2 Monitor next step feasibility: a conservative approach

Motivated by the good results of the OP strategy, we ground our patient freedom saturation on the same reachability problem for the input-constrained LIP dynamics (2.17). The fixed-horizon OCP writes the same way as Problem 4.1, but the set of feasible commands Ω is adapted to include the CoM path constraint (B.10), yielding

Ω(χ 0 , χ f , T ) ≜        u ∈ U ad (T ),        χ u (0) = χ 0 χ u (T ) = χ f Cχ u (t) ∈ γ, ∀t              
with C = 1 0 0 0 0 0 1 0 and γ the manifold defined by the nominal CoM

path γ = x ∈ R 2 , ∃t ∈ [0, T f ] s.t. x = c nom (t) (B.12)
In general, γ is a non-trivial manifold. Hence, the constraint Cχ u (t) ∈ γ, ∀t can not be written as an affine constraint, which makes Problem 4.1 impossible to transcribe as a QP problem. This discards the solving methodology of Problem 4.2 proposed in Chapter 4.

Luckily, further investigations of the nominal CoM trajectory during double support 1 , depicted on This yields a few important remarks

• the X and Y CoM dynamics are linearly dependent, hence we only need to solve Problem 4.2 along one axis;

• the CoM velocity and acceleration are tangential to the CoM db 0 CoM db f vector;

• the CoP location is constrained to lie on the segment defined by the intersection of, on one hand, the (infinite) line passing through CoM db 0 and CoM db f , and the double support polygon, on the other hand (by (2.17), the remark above and the unilateral contact constraints). • the saturation computed using dimensions X and Y do not exactly overlap, which is a consequence of the real final nominal CoM velocity is not exactly tangential to the CoM db 0 CoM db f vector.

This alternative control law for double support phases has only been tested in simulation with arbitrarily shaped patient schedules η to simulate the impact of a patient. It acts as expected and stabilizes simulations with simulated slowdowns during the double support phase which makes the system fall without the proposed saturation. Experiments supporting these simulation observations remain to be done.

  Exoskeletons for walking rehabilitation . . . . . . . . . 18 1.1.1 A concise view of walking rehabilitation . . . . . . 18 1.1.2 Using exoskeletons during training . . . . . . . . . 20 1.2 Passive walking with a self-balanced exoskeleton . . . . 24 1.2.1 Description of the Atalante platform . . . . . . . . 24 1.2.2 Passive walking control architecture: plan then execute . . . . . . . . . . . . . . . . . . . . . . . . 25 1.3 Active gait training with self-balanced exoskeletons . . . 28 1.3.1 Spatial and temporal freedom for active gait training 28 1.3.2 Exploiting spatial freedom with self-balanced exoskeletons . . . . . . . . . . . . . . . . . . . . . . . 29 1.3.3 Contribution of the thesis: exploit temporal freedom and maintain balance . . . . . . . . . . . . . . 30

Figure 1 . 1 :

 11 Figure 1.1: Phases of gait cycle (from [Kour2019]).

20

  Figure 1.2: Walking rehabilitation examples (from Helen Hayes Hospital https://youtu.be/g__BYaS9viw).

  Atalante Figure1.3d), designed to balance themselves during walking and not requiring crutches or external help.

  Figure 1.3: Exoskeleton examples: (a) platform-based exoskeleton, (b) crutched exoskeleton, (c) self-balanced exoskeleton, (d) self-balanced exoskeleton.

Figure 1 . 5 :

 15 Figure 1.5: Elements of Atalante's passive walking controller.

  Active gait training with self-balanced exoskeletons 1.3.1 Spatial and temporal freedom for active gait training A variety of active walking controllers have been proposed for active gait training with platform-based or crutched exoskeletons. They allow the patient to participate in the motion by relaxing the tracking of a reference trajectory 15 .

  Figure 1.6: Comparison between spatial and temporal freedom on one joint trajectory. The larger envelope on the right makes temporal freedom easier to perceive for patients.

Figure 1

 1 Figure 1.7: Proposed control architecture for active gait training: splitleg controllers carefully synchronized with patient schedule through online replanning.

Figure 2 . 1 :

 21 Figure 2.1: Definition of the sagittal, frontal and transverse planes (adapted from [Winter2009].

Figure 2 . 2 :

 22 Figure 2.2: Atalante picture (Left). Atalante schematic hardware (Right).

Figure 2 . 3 :

 23 Figure 2.3: Transitions graph of flat-foot trajectory.

Figure 2 . 4 :

 24 Figure 2.4: CoM, CoP and DCM of the trajectory generated by the flat-foot OCP (no double support phase between single support ones).

  4 s for the (STARTD) and (STOPD) phases. The step length is fixed approximately to 40 cm. The CoM, CoP and DCM of the trajectory solution of the foot-rolling OCP are displayed on Figure 2.7, and the swing foot clearance during the (CS) phase is illustrated on Figure 2.6 (Right). The comparison of the swing foot clearance given on Figure 2.6 shows the improvements offered by the foot-rolling OCP over the flat-foot OCP, namely: longer step-length, higher foot clearance and more anthropomorphic toe-off/heel-strike phases. The average translational velocity is approximately 40 cm/s.

Figure 2 . 6 :

 26 Figure 2.6: Swing foot clearance comparison between the trajectory generated by the flat-foot (Left) and the foot-rolling (Right) OCP. Lines representing the sole of the swing foot at various instants of the (CS) phase, for the flat-foot OCP, and the lumped (CD)-(CS)-(CD) phase, for the foot-rolling OCP in the XZ-plane. The foot-rolling steps are longer and more anthropomorphic.

  Figure 2.7: CoM, CoP and DCM of the trajectory generated by the footrolling OCP (with double support phase between single support ones).

SoT 1 Level 0 :

 10 tracking of swing and support feet pose by (2.40); Level 1: tracking of CoM acceleration by (2.42), realizing (2.19); Level 2: tracking of back orientation by (2.43); Level 3: tracking of joint by (2.45).

  Figure 3.1: Curvilinear abscissa s (Left) and curvilinear velocity ṡ (Right) over one step, computed from the swing leg trajectory T sw obtained from the flat-foot OCP (notice that s is strictly monotonous).

  Figure 3.3: Curvilinear abscissa convergence may vary with initialization.

Figure 3 . 4 :

 34 Figure 3.4: Normalized f (defined by (3.20)) with qsw = 0: existence of local minima depends on the fictitious time t f i .

  Figure 3.5: Normalized f and Euclidean distance with qsw ̸ = 0.

68CHAPTER 3 .

 3 Figure 3.6: Normalized f with qsw = 0 and a f i ̸ = 0: orthogonal components create and move local minima.

  Figure 3.7: Torque norm and Euclidean distance (null measured velocity and a f i = 0) with anisotropic gain matrices.

  Figure 3.8: Simulated joint position and velocity (Black: nominal state (T sw , Ṫ sw ), blue: time rescaled nominal state (P (σ), T (σ) σ), orange: fictitious state (q sw f i , qsw f i )).

Figure 3 . 9 :

 39 Figure 3.9: Experimental joint position and velocity (Black: nominal state (T sw , Ṫ sw ), blue: time rescaled nominal state (P (σ), T (σ) σ), orange: measured state (q sw m , qsw m )).

Stack of Tasks SoT 2

 2 Level 0: swing leg joint position and support foot pose tracking; Level 1: CoM acceleration tracking (2.19); Level 2: back orientation tracking; Level 3: regularizing joint tracking.

  32) 74CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING with K sw p , K sw d ≻ 0 constant diagonal gain matrices. As the swing and support legs are only defined during single support phases, we keep the SoT 1 during double support phases 3 .

Figure 3

 3 Figure 3.11: Split-leg rehabilitation controller using the naive TR synchronization strategy.

  Figure 3.12: Atalante's kinematic chain during single support phases.

  36) where all the quantities in the second line are independent of the acceleration qsp of the 6 DoF of the support sub-chain. To highlight the dependence of 78CHAPTER 3. EMPOWERING THE PATIENT FOR GAIT TRAINING the various factors w.r.t. the variables of both legs, we express c as c =

  Figure 3.13: CoM acceleration IK task error in simulation (same error magnitude of the original and proposed algorithm). Original algorithm: SoT 1 . Proposed algorithm: SoT 2 . Left: forward direction. Right: lateral direction.

82CHAPTER 3 .

 3 Figure 3.15: η u (Left) and ηu (Right) functions, along with the three parameters (emulation of patient-induced slow-downs).

84CHAPTER 3 .

 3 Figure 3.17: Heatmap of time-rescaled nominal trajectories violating the unilateral force constraint. Rate of time-rescaled trajectories respecting the constraint computed for each (custom velocity duration, custom velocity) couple, aggregating the (custom velocity start) dimension.

Figure 3 .

 3 Figure 3.18: Heatmap of 10 s walking experiments success. Success rate computed for each (custom velocity duration, custom velocity) couple, aggregating the (custom velocity start) dimension. The results are in accordance with simulation, with lower performance in the low-velocity range.

86CHAPTER 3 .

 3 Figure 3.19: Reference CoP u T R y (violating the support polygon limits in black), CoM c T R y and measured CoM c y (diverging toward the end of the step) positions along the Y axis of the inertial frame W. Reference quantities computed using the TR strategy and a simulated patient velocity σu as low as 60% (t η = 0.2 s, ∆η = 0.9 s and η = 0.6). Black horizontal lines represent the limits of the support polygon SP on the Y axis.

Figure 4

 4 Figure 4.1: Split-leg rehabilitation controller using the careful OP synchronization strategy.

  .15) (a) Reference CoP u OP y , CoM c OP y and measured CoM c y positions along the Y axis of the inertial frame W. Black horizontal lines represent the support polygon limits on the Y axis.(b) Simulated target velocity σt (with t η = 0.2 s, ∆η = 0.9 s and η = 0.6) and optimal velocity σopt .

Figure

  Figure 4.2: Reference quantities computed using the OP strategy and a simulated target velocity. The support polygon limits are not violated using the OP strategy, the CoM reference is successfully tracked, and the simulated target velocity is followed as much as possible.

Figure 4 . 3 :

 43 Figure 4.3: Illustration of the main result, Theorem 1.

4. 2 .

 2 BISECTION ON TRAJECTORY FINAL TIME 103 (a) x dimension (b) y dimension

Figure 4 . 4 :

 44 Figure 4.4: Feasible velocities set V(χ t , χ f ) along the nominal CoM trajectory (computed from T ). A hole appears at the end of the step for the x dimension.

Figure 4 .

 4 Figure 4.4 displays the set of feasible velocities V for each dimension x and y, over time, estimated by exhaustive numerical search.Figure4.4b shows that V y (χ t y , χ f y ) is an interval for all initial conditions χ t (varying with time). However, Figure4.4b reveals that V x (χ t x , χ f x ) is the union of two intervals for initial conditions χ t evaluated at times near the end of the nominal trajectory duration. This numerically illustrates our main Theorem 1, as it states that T can be the union of two intervals.

Figure

  Figure 4.4 displays the set of feasible velocities V for each dimension x and y, over time, estimated by exhaustive numerical search.Figure4.4b shows that V y (χ t y , χ f y ) is an interval for all initial conditions χ t (varying with time). However, Figure4.4b reveals that V x (χ t x , χ f x ) is the union of two intervals for initial conditions χ t evaluated at times near the end of the nominal trajectory duration. This numerically illustrates our main Theorem 1, as it states that T can be the union of two intervals.

Figure 4 .

 4 Figure 4.5: Feasible velocities set V(χ t , χ f ) along the nominal CoM trajectory (computed from T ). The hole of the x dimension is filtered out by the intersection with the y dimension.

Figure 4 . 6 :

 46 Figure 4.6: Velocity of the trajectory (1 is the nominal velocity). (Dotted blue): target velocity ηt . (Green dots): set T determined by an exhaustive search, for reference. (Solid black): optimal velocity ηopt from the OP strategy.

Figure 4 . 7 :

 47 Figure 4.7: Comparison of rate of success heatmaps for velocity variations having various durations and magnitudes 7 . The OP strategy increases the stabilization performance in the low-velocity range (w.r.t. the TR strategy). Left: TR strategy. Right: OP strategy.

Figure 4 . 8 :

 48 Figure 4.8: Comparison of experimental stability for velocity variations having various durations and magnitudes 9 . Experimental results are in accordance with the simulations, with lower performance in the low-velocity range. Green: stable without external help. Orange: stable with practitioner help. Red: unstable. Left: TR strategy. Right: OP strategy.

Figure 4 . 9 :

 49 Figure 4.9: Effective velocity σopt and target velocity σt over a step with an able-bodied user (step 6 of Figure 4.11). Black curves: lower and upper limits of the set of feasible velocities.

Figure 4 .

 4 Figure 4.10: CoM from TR and OP strategies over a step with an able-bodied user (step 6 of Figure 4.11). Grey areas: double support phases. White areas: single support phases. Left: X axis. Right: Y axis.

Figure 4 .

 4 Figure 4.11: Effective velocity σopt over a 10-steps experiment with an ablebodied user. Grey areas: double support phases. White areas: single support phases.

Figure 4 .

 4 Figure 4.12: Effective velocity σopt and target velocity σt over a step comprising a full stop, with an able-bodied user. Black curves: lower and upper limits of the set of feasible velocities.

Figure 4 .

 4 Figure 4.13: CoM from TR and OP strategies over a step comprising a full stop, with an able-bodied user. Grey areas: double support phases. White areas: single support phases. Left: X axis. Right: Y axis.

Proof.

  Consider the HamiltonianH(t, λ 0 , λ, x, u) = µ + λ(t)(Ax + Bu)Using the PMP, the adjoint equation and solution writeλ = -∂H ∂x = -λ(t)A λ(t) = λ 0 e -At

Figure A. 1 :

 1 Figure A.1: Phase diagram for (A.2) with u m = -1 (red) and u M = 2 (blue).C m covers R 1 ∩ D + , R 5 ∩ D -, R 9 ∩ D -and R 8 . C M covers R 1 ∩ D + , R 5 ∩ D + , R 9 ∩ D -and R 2 .

Figure A. 2 :

 2 Figure A.2: Graph of all possible cases of non-empty T.

A. 4 .

 4 CONVEXITY OF BOUNDED T CASES 123 A.4 Convexity of bounded T cases Lemma 5. When T is bounded, T = [T min , T max ].

Figure A. 3

 3 Figure A.3 shows the trajectory steering the system from x 0 to x f using seq m and seq M . The trajectories do not cross each other. The trajectory corresponding to seq m and the reverse trajectory corresponding to seq M define a positively oriented closed curve. The region R 147 ∩ D -being positively invariant, the curve is strictly included in D -2 .For all x ∈ D -,

  Figure A.1 (all possible values of x 0 are colored in red when x f varies along the dashed line). It plays a particular role in Lemma 6 as it is the only one where boundary conditions yield a non-convex set T. Lemma 6. When T is unbounded, if (x 0 , x f ) / ∈ J + ∪J -, then T = [T min , +∞[, otherwise, there exists A < B s.t. T = [T min , A] ∪ [B, +∞[. Proof. Following Lemma 4, a careful investigation of the graph in Figure A.2 reveals that for T to be unbounded we have (x 0 , x f ) ∈ R 258 × R 456 . By symmetry of the vector field (rotation of π about ( um+u M 2 , -um+u M 2 ) T ), we now only consider a pair (x 0 , x f ) ∈ R 25 ∩ D + × R 456 . (x 0 , x f ) / ∈ J + In all such cases, there exists a sequence (u m , a, b, c, d), with a ≥ 0, b > 0, c > 0, d ≥ 0 steering x 0 to x f with a single intersection with D ∩ R 5 . This sequence can be easily extended in the vicinity of D ∩ R 5 (which excludes equilibria) to increase the transient time by any desired arbitrarily small increment ϵ > 0. Iteratively, this construction allows to infinitely increase the transient time by a continuous constructive process. Also, the same type of sequence with other values for a, b, c, d can generate a smooth collection of trajectories approaching the minimum time trajectory. The proof of Appendix A.4 yields the conclusion with the continuous mapping T : (a, b, c, d) → a + b + c + d.

Figure A. 3 :

 3 Figure A.3: Phase diagram for (A.2) with u = u m (dotted red), u = u M (dotted blue), minimum time trajectory x min (solid blue line), maximum time trajectory x max (solid red line), and x seq v for some a ∈ [0, a m ] (dashed black line).

  Figure B.1: Nominal CoM trajectory during one step (double and single support) and footprint of the support polygon during double support. The nominal CoM path during double support is almost a segment (the dotted orange segment).

  Figure B.2: Backward trajectories from x f (final state of flat-foot trajectory with double support). Nominal trajectory in green.

136

  Figure B.3: Lower and upper saturation bounds computed using backward CoM trajectories either along X or Y axis. Bounds differ with the dimension used for their computation because the nominal CoM path during the double support phase is not exactly a segment.

  SoT used by the original passive walking controller of Atalante SoT 2 SoT used by the proposed active walking controller during single support phase Q nominal trajectory mapping time in R + to the space of generalized positions in R 18 reference CoP command in R 2 [m] p d desired CoP position in R 2 [m] p T R time-rescaled CoP trajectory mapping curvilinear abscissa σ ∈ R + to CoP in R 2 [m] p OP optimal CoP trajectory mapping curvilinear abscissa σ ∈ R + to CoP in R 2 [m] x T R time-rescaled CoM state trajectory mapping curvilinear abscissa σ ∈ R + to the CoM state in R 4 [m, m.s -1 ] x OP optimal CoM state trajectory mapping curvilinear abscissa σ ∈ R +

			Glossary
		ċ	CoM velocity the inertial frame, in R 3 [m]
		ċ *	horizontal reference CoM velocity in R 2 [m.s -1 ]
		c *	horizontal reference CoM acceleration in R 2 [m.s -1 ]
		cd	desired CoM acceleration in R 3 [m.s -2 ]
	ω Glossary LIP natural frequency in R + [s -1] ξ DCM in R 2 [m]
		ξ *	reference DCM in R 2 [m]
		p
	CBF Control Barrier Functions CoM Center of Mass CoP Center of Pressure DCM Divergent Component of Motion DoF Degrees of Freedom FK Forward Kinematics HZD Hybrid Zero Dynamics ID Inverse Dynamics IK Inverse Kinematics IMU Inertial Measurement Unit LIP Linear Inverted Pendulum MPC Model Predictive Control NLP Non Linear Program OCP Optimal Control Problem OP Optimal Planning PD Proportional Derivative QP Quadratic Program SoT 1 η u simulated user phase variable in R [s]
	SoT ηu	Stack-of-Tasks simulated user phase variable velocity in R
	TR m	Time Rescaling Total mass of the patient-exoskeleton system [kg]
	VG L	Virtual Guides angular momentum about the CoM in R 3 [rad.s -1 ]
	c	CoM position w.r.t. the inertial frame, in R 3 [m]
	c *	horizontal reference CoM position in R 2 [m]
			13

T nominal joint trajectory mapping time in R + to joint space in R 12 [rad] T sw nominal swing leg joint trajectory mapping time in R + to joint space in R 6 [rad] P parametric curve mapping curvilinear abscissa in R + to joint space in R 6 [rad] T Frenet-Serret unit tangent vector to the curve P mapping curvilinear abscissa in R + to joint space in R 6 [rad] s curvilinear abscissa mapping time in R + to arc length, in joint space, in R + [rad] σ control variable in R + [rad] σ control variable derivative in R [rad.s -1 ] σu simulated user velocity in R [rad.s -1 ] *

  Passive walking is the mode in which exoskeletons fully assign the patients' leg motion, irrespective of the patient's efforts. This mode allows passive gait training (i.e. rehabilitation with no participation required from the patient,

see

[Quintero2012; Bortole2015]

). By contrast, active walking is the mode in which exoskeletons allow patients to significantly modify the motion of the patient-exoskeleton system through their efforts (see Section 1.3.1 for examples). Such active walking modes are used for active gait training (i.e. rehabilitation with participation required from the patient). The control methodology of Atalante described in Section 1.2.2 allows passive gait training.

  The nominal SoT 1 of Atalante . . . . . . . . . . . 53

	34	CHAPTER 2.	ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER
	Contents
		2.1 Atalante: description and dynamical models . . . . . .	35
			2.1.1 Definitions of the frames of reference . . . . . . . . 35
			2.1.2 Description of Atalante . . . . . . . . . . . . . . . 35
	Chapter 2	2.1.3 Kinematic model . . . . . . . . . . . . . . . . . . . 36 2.1.4 Lagrangian dynamics . . . . . . . . . . . . . . . . 39
			2.1.5 Newton-Euler dynamics . . . . . . . . . . . . . . . 40
		2.2 Offline trajectory generation . . . . . . . . . . . . . .	40
			2.2.1 Generic walking motion optimal control problem . 40
	2.2.2 Flat-foot trajectory optimal control problem . . . 43 Atalante: models and passive 2.2.3 Foot-rolling trajectory optimal control problem . . 44
	2.3 Center of Mass stabilization and admittance controller . walking controller 2.3.1 DCM control using admittance . . . . . . . . . . . 49 47
			2.3.2 Hierarchical Inverse Kinematic problems: the
			Stack-of-Tasks framework . . . . . . . . . . . . . . 50
			2.3.3
				Marine Pétriaux, Florent Di Meglio, and Nicolas
		Petit. "Fast Replanning of a Lower-Limb Exoskeleton Trajectories
		for Rehabilitation". In: 2022 IEEE 61st Conference on Decision and
		Control (CDC). Dec. 2022
		• Maxime Brunet, Marine Pétriaux, Florent Di Meglio, and Nicolas Petit.
	Enabling Safe Walking Rehabilitation on the Exoskeleton Atalante: Chapitre 2 -Atalante: modèles et contrôleur de marche passif: Dans ce
	Experimental Results. Apr. 17, 2023. preprint to appear in 2023 IEEE chapitre, nous décrivons le robot Atalante et les principaux modèles utilisés
	pour décrire le système patient-exosquelette. Ensuite, nous décrivons le International Conference on Robotics and Automation (ICRA) contrôleur nominal d'Atalante qui permet de réaliser des exercices passifs de
	rééducation de la marche, composé de deux éléments principaux. Le premier A video describing the main experiments done in this thesis is available est un problème de commande optimale résolu hors-ligne pour obtenir des at https://youtu.be/_1A-2nLy5ZE. trajectoires du système. Nous décrivons en particulier deux trajectoires de
	marches qui seront utilisées comme cas d'usage du travail présenté dans ce
	manuscrit : une trajectoire de marche simple, sans déroulé du pied, et une
	trajectoire de marche plus complexe, avec un déroulé du pied. Le deuxième
	élément composant ce contrôleur est une loi de contrôle utilisée en ligne pour
	stabiliser ces trajectoires lors de leur exécution.

  Then, because there is no physical actuator directly controlling the CoP p, admittance control[Caron2019] is used to drive it to the desired value p d (2.18). The admittance control methodology encompasses a general class the measured CoP, c * the horizontal state trajectory acceleration, and A ≻ 0 a diagonal gain matrix to be tuned. This admittance task is included in a hierarchical Inverse Kinematics (IK) problem called a Stack-of-Tasks (SoT)[Mansard2009] framework. The hierarchy yields that this equation holds at almost all times. We describe this problem and SoT framework in the next section.

	CHAPTER 2.	ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER
	of control laws having a force regulation objective to be satisfied through
	kinematic actuators. This methodology considers one or several kinematic
	feedback laws, called admittance tasks, to regulate the force errors to zero.
	It is used, for instance, in [Yokoi2003; Englsberger2012; Li2012].
	On Atalante, there is only one admittance task, a horizontal CoM ac-
	celeration task, which defines the desired horizontal CoM acceleration cd as
	follows
		cd = c * + A(p -p d )	(2.19)
	with p

x,y , ċx,y ) (see (2.8)), ξ * the DCM reference trajectory, computed from the horizontal state trajectory (c * , ċ * ), p * the associated CoP trajectory, and k p , k d , k i ≻ 0 three diagonal matrices. Note that the actual CoP being bounded, regulating the DCM instead of the CoM gives the "Best CoM-ZMP 9 regulator" in a sense defined in

[Sugihara2009]

. Formula (2.18) encodes a feedback which, for any k p > 0, k i = 0, k d = 0, drives any initial condition of the LIP in the controllable region to 0. [Caron2019] extends the formula to k i > 0, k d > 0.

  .8. Dans ce chapitre, nous proposons un nouveau contrôleur pour l'exosquelette Atalante qui donne au patient la possibilité de participer au mouvement. Plus précisément, nous laissons le patient modifier la vitesse de la jambe de vol en temps réel tout en stabilisant, grâce à la jambe de support, le centre de masse du système autour de la trajectoire de référence synchronisée sur la jambe de vol. Nous utilisons le principe des guides virtuels sur la jambe de vol pour libérer le mouvement du patient le long de la trajectoire de référence. Nous modifions la technique de stabilisation utilisée par le contrôleur passif d'Atalante pour n'utiliser que la jambe de support. Une simple dilatation temporelle de la trajectoire nominale est utilisée pour synchroniser la trajectoire de référence du centre de masse avec la jambe de vol. Les performances de stabilisation de ce contrôleur sont évaluées en simulation et à l'aide d'expériences avec un mannequin. Choice of the Virtual Guides definition space . . . 60 3.1.2 Construction of a parametric curve of reference . . 60 3.1.3 Virtual Guides controller . . . . . . . . . . . . . . 61 3.1.4 Behavior of Virtual Guides . . . . . . . . . . . . . 63 3.2 Stabilizing the Center of Mass with the support leg . . 72 3.2.1 Modification of admittance's SoT . . . . . . . . . . 73 3.2.2 Synchronization of legs by Time Rescaling (TR) . 74 3.2.3 Trajectory smoothing . . . . . . . . . . . . . . . . 74 3.2.4 Impact of the modification of the SoT . . . . . . . 76 3.3 Impact on safety of patient empowerment . . . . . . .

	CHAPTER 2.	ATALANTE: MODELS AND PASSIVE WALKING CONTROLLER
	Contents			Patient
	3.1 Satisfying the patient intent with the swing leg . . . . .	59
	Chapter 3	3.1.1		Exoskeleton
				)
			, τ sp
			(τ sw
	PD ctrl. (2.46) Empowering the patient for gait , q sp (q sw training
				q t
			q *	CoP ctrl. (2.19) cd SoT 1 IK c *	Stabilization controller
				p d
				p *	DCM ctrl. (2.18)
				ξ *
			Offline	Trajectory	Q

) p ξ Figure 2.8: Atalante's passive walking controller: offline trajectories stabilized with admittance control. Chapitre 3 -Donner de la liberté au patient: 80 3.3.1 Stability analysis . . . . . . . . . . . . . . . . . . . 80 3.3.2 Simulation results . . . . . . . . . . . . . . . . . . 81 3.3.3 Experimental observations . . . . . . . . . . . . . . 83 3.3.4 Performance analysis during a single step simulation 86

  6) with t f the final time defined by the first time the control variable σ reaches the end of the curvilinear trajectory s i.e. σ(t f ) = L max .The CoM trajectory c is defined by the parametric curve P and the control variable σ according toc(t) = FK • P • σ(t) = FK • T • η(t)The position c of the CoM trajectory follows the nominal CoM path c nom . Consequently, the final CoM position c(t f ) is located at the final nominal CoM location c nom (T f ), i.e. c(t f ) = c nom (T f ) (B.10) • The velocity ċ(t) of the CoM trajectory is tangential to the nominal CoM velocity dc nom dη , time-rescaled by the patient schedule η. Consequently, the final CoM velocity ċ(t f ) is also tangential to the final nominal CoM velocity dc nom dη (T f ), i.e.

	B.2.	MONITOR NEXT STEP FEASIBILITY: A CON-SERVATIVE APPROACH	133
	• ċ(t f ) = η(t f )	dc nom dη	(T f )	(B.11)
						(B.7)
		By definition of T ,
			c(t) = c nom • η(t)	(B.8)
		with c nom ≜ FK • T the nominal CoM trajectory. Hence, the CoM velocity
		ċ reads	ċ(t) = η(t)	dc nom dη	• η(t)	(B.9)
		Inspecting (B.8) and (B.9), one can easily notice that

Neurological disorders encompass numerous disorders related to nervous system disorders such as spinal cord injuries, cerebral palsy, or stroke.

Ambulation is the act of moving around by walking. There is no rehabilitation objective during ambulation, only a mobility one.

We refer to realism of the walking for the ratio of patient weight supported by their legs.

We refer to safety of the walking for the probability of falling during walking.

We refer to intensity of the walking exercise for the duration of the exercise and the level of effort produced by the patient's legs.

Roughly speaking there is a 1-1 ratio between the mass of the human and the mass of the exoskeleton.

The benefit of a third family of control methodologies [Castillo2021; Siekmann2021; Singh2022], leveraging the recent progress of machine learning, is still under investigation for self-balanced exoskeletons [Duburcq2020; Duburcq2022b] and does not constitute the primary control paradigm of Atalante.

The same control methodology can be applied to realize a variety of motions like walking in a straight line, side-ways, turning in place, or climbing stairs (see[Huynh2021] and accompanying video https://youtu.be/4bmuNlnREaA).

The reduced-state trajectory being generated using an approximate model of the dynamics, some conservatism during planning must be considered, which limits the performance. Some of the most advanced works [Mansard2018; Dantec2022] on humanoid robots control exhibit online motion generation and tracking considering the full dynamics of the robots to alleviate this limitation. However, even by today's computational power standards, these algorithms are still too heavy a computational burden to run online at a high update rate on the embedded computers available in robots and have to run on remote computers. Other recent approaches attempt to ensure a consensus between reduced and complete dynamics, as proposed in[Budhiraja2019].

 14 An alternative is to compute the torques directly by Inverse Dynamics (ID) [Her-zog2016;Kuindersma2016; Koolen2016]

Some control laws do not use any reference trajectory. This is the case, for instance, in[Ekkelenkamp2005; Vallery2008; Murray2015]. However, these paradigms deviate further from Atalante's passive walking paradigm, therefore we do not rely on them for inspiration.

A PD controller with deadband is employed to help the patient to stay in the vicinity of the reference trajectory while offering spatial freedom around the reference trajectory.

It is updated during the walk.

This Henke axis is contained in the sagittal plane, and is equivalent to the X axis depicted on Figure2.1, but tilted by

degrees upward around the Y axis.

The velocity and acceleration of each joint are obtained through numerical differentiation of their position.

The patient efforts are internal to the patient-exoskeleton system. Hence, their only effect appear in the joint position and can be expressed as additional joint torques.

Paraplegic patients cannot voluntarily produce efforts with their legs. Patients suffering from spinal cord injuries may produce erratic involuntary efforts (this syndrome is called spasticity). The efforts are unpredictable by nature and are therefore unmodeled.

The duration of each step being free to choose by the optimizer, the waypoints are defined w.r.t. the total duration of each step, and not at absolute times.

The stabilization algorithm also helps to mitigate all the uncertainties like the unknown patient model, the unmodeled effects like the flexibilities of the structure[Vigne2021], the actuator dynamics and the dynamics of their low-level current control laws.

The CoP is only indirectly commanded by the joint actuators through the (unknown) contact equation, hence we call it a virtual actuator.

On planar ground, the Zero Moment Point (ZMP) and the CoP coincide at all times. Hence, we do not distinguish them in this thesis.

We explore the possibility to modify the controllers during the double support phases in Appendix B.

Keeping during the double support phases is the safest approach, using all DoF to stabilize the CoM dynamics. An opposite approach, using all DoF to follow the patient schedule, is sketched in Appendix B.

Choosing the CoM location is not equivalent to choosing the position of the back, which defines the landing location of the swing leg. However, the latter cannot directly be included in a LIP-based OCP. Choosing the CoM position at the end of the step is a surrogate, which is reasonable as the back frame is roughly at the center of the mass distribution of the patient-exoskeleton system.

According to the LIP model, minimizing the integral of any derivative of the CoM over a sufficiently long enough horizon is enough[Wieber2019]

The presented quadratic cost function can easily be changed to incorporate extra tuning parameters to affect performance, without loss of generality.

We treat the feasibility function as a boolean.

The white spaces in this figure correspond to unfeasible values of the parameters violating the constraint σ < Lmax.

More precisely, the reference trajectory used during double support is computed once, at the beginning of the step. For this, we use the OP strategy and the mean velocity of the previous step, for the sake of continuity.

The OCP used to generate the flat-foot trajectory, described on Section

2.2 in the YX-plane, has been modified to include a double support phase for the purpose of the development presented in this appendix.

Remerciements

This simulation reveals that the choice of the reference trajectory according to the TR strategy is too naive to ensure the safety of the walking exercise, and calls for a more careful choice of legs synchronization strategy. This is the topic of the next chapter.

As detailed in Section 4.2, the x and y dynamics of the LIP model (2.17) are totally decoupled, and the Theorem 1 is about the nature of T(χ 0 , χ f ) in one dimension. Therefore, we conduct the proof of the theorem in one dimension and χ notations in place of x or y dimension will be omitted in this appendix.

The proof of Theorem 1 is organized as follows. First, after some preliminary change of coordinates and definitions in Appendix A.1, we exhibit in Appendix A.2 particular regions of the phase portrait which serve to decompose the proof. Then, we study the boundedness of T in Appendix A.3. Finally, we assess its convexity properties in Appendix A.4 and Appendix A.5 by differentiating according to the boundedness of T.

A.1 Preliminaries

A.1.1 Decoupling dimensions

To ease the analysis of the set of feasible times T(χ 0 , χ f ), we perform the following change of coordinates, which diagonalizes the dynamics (2.17) with

Then, (2.17) takes the diagonal form ẋ ≜ Ax + Bu with

The solution of (A.2) with input u, from the initial condition x 0 ∈ R 2 is denoted x u . For all vector variables, a subscript 1 or 2 indicates the first or second coordinate. Zero-order hold of u for a duration d = t 2 -t 1 yields the solution

The change of coordinate between χ and x is trivial, s.t. the set of feasible times T(χ 0 , χ f ) can be equivalently denoted T(x 0 , x f ). For brevity purposes, we omit the dependence of the set of feasible times T on (χ 0 , χ f ) when convenient.

A.1.2 Definitions

A piecewise constant control input u having N steps over an interval [0, t f ] is defined using a finite (irreducible) partition

with u taking values only in {u m , u M }. For convenience, it is described by its first value and the durations, e.g. for 3 steps of respective durations d 1 , d 2 , d 3 starting with u m , a sequence

For any initial condition x 0 , and any seq

A.2. REGIONS OF INTEREST IN THE PHASE PORTRAIT

and the switching function is

If λ 0 1 λ 0 2 < 0, then there exists a unique switching time 1 2ω log(-

) for which Γ changes sign. Otherwise, Γ has a constant sign. This concludes the proof for u min . The proof regarding u max is identical.

Lemma 2 highlights the importance of the phase portrait in Figure A.1 corresponding to constant control values u m and u M . It is split into nine open regions, some of them being open semi-infinite strips, whose boundaries are the trajectories passing through the equilibrium points for u m and u M . We denote each region R i , i = 1, ..., 9. Also, we will note R ijk... ≜ R i ∪R j ∪R k ∪... for any number of indexes. Notice two interesting properties: i) the locus of intersecting parallel arcs is D and ii) the subsets R 147 and R 369 are positively invariant under the controlled flow.

Next, the following result states that in the cone C m (resp. C M ), the flow corresponding to u m (resp. u M ) reaches a point symmetric to the initial condition w.r.t. the line D. This property is instrumental in the proof.

Lemma 3. For all x in the double cones C m ∪ C M , we have

Proof. In the double cones C m and C M , f is well-defined as a function of its arguments. A direct calculus with (A.4) yields the conclusion.

A.3 Boundedness of the set of feasible times T

Depending on the values of x 0 and x f , the set T can be empty (∅), bounded (Bd), or unbounded (∞). which maps the duration a of the first arc of x seq v to the total duration a + b + c. Define g as

The intermediate point The proof is identical, replacing D -by D + , the trajectory corresponding to seq m and the reverse trajectory corresponding to seq M defining a negatively oriented closed curve.

If the (Euclidean) distance d(x f , D) between x f and D is strictly lower than the distance d(x 0 , D) between x 0 and D, then Lemma 3 states that ϕ(x f , u M , -t 0 ) = Sx f = ϕ(x f , u m , -t 1 ), for some t 0 , t 1 ≥ 0. We use the same constructive proof between Sx f and x f with the sequence

, for some t 0 , t 1 ≥ 0. We use the same constructive proof between Sx 0 and x 0 with seq B .

If d(x f , D) = d(x 0 , D), the proof directly follows from x 0 = Sx f to x f , this situation is illustrated in

Now, consider a trajectory from x 0 to x f leaving Γ 1 ∪ Γ 2 . A detailed investigation of the phase portrait gives that this trajectory leaves

From all solutions, starting from x ii , reaching x f and leaving Γ 1 ∪ Γ 2 at x ii , the one passing through Sx ii displayed on Figure A.5 is the only bang-bang solution. Hence, it is the fastest according to Lemma 2. We denote Γ ii the union of this solution and the fastest solution red from x 0 to x ii and construct Γ i similarly by passing through

We now show that Γ i is faster than Γ ii . As they share the arc from Sx i to x f , we only need to use that the one-arc solution from x 0 to Sx i is faster than the 3-arc solution passing through x ii and Sx ii . This holds because the only bang-bang solution from x 0 to Sx i is the minimum time solution (according to Lemma 2). This holds for all x ii , hence all solution leaving Γ 1 ∪ Γ 2 is slower than Γ i . As Γ i is crossing D, it can be arbitrarily extended. This completes the proof. In this appendix, we sketch an alternative control methodology for active gait training during double support phases which aims at enabling the patient to participate in the motion during these phases. It could be beneficial to patients' rehabilitation because transferring the weight from one leg to the other, which occurs during double support phases, is also important to regain ambulation capabilities. The sketched control law offers temporal freedom around the offline trajectory using the same VG methodology used during single support phases (see in Section 3.1), with a substantial difference: all the joints of the exoskeleton are constrained to the parametric curve P (computed from the offline trajectory T instead of T sw ). Hence, no DoF is left free for the stabilization of the CoM dynamics. This is the opposite choice to the one described in Chapter 3, where all the joints are used for the stabilization of the CoM dynamics during double support phases.

In addition, we propose a saturation of the patient schedule, during double support, which ensures the feasibility of the next single support phase. This ensures the OP strategy can still be used during single support phases.

MOTS CLÉS

Rééducation robotique, Exosquelette, Stabilisation, Marche bipède, Contrôle optimal RÉSUMÉ Les exosquelettes auto-équilibrés, grâce à leurs algorithmes de contrôle qui permettent une marche autonome et stable, offrent aux personnes handicapées des capacités de marche retrouvées. Cette thèse vise à permettre aux patients d'utiliser ces systèmes à des fins de rééducation selon le concept de rééducation active de la marche. Selon notre interprétation, ce concept implique de redonner une certaine liberté de mouvement aux patients. Cet objectif est en conflit avec la garantie de la sécurité, qui implique de prendre le contrôle du mouvement. Le principal sujet abordé dans cette thèse est le défi que représente la conciliation de ces deux objectifs. Notre solution est double : donner de la liberté au patient et en garantissant l'équilibre. Tout d'abord, nous introduisons une division dans l'architecture de contrôle existante, permettant aux patients de participer activement au mouvement de la jambe de vol. Parallèlement, l'algorithme de contrôle de la jambe de support régule les forces de contact. Deuxièmement, nous résolvons de manière embarquée des problèmes de contrôle optimal en temps final libre sur un modèle réduit afin de générer des trajectoires stabilisables. La durée de ces trajectoires sert à déterminer le moment où le système doit reprendre le contrôle de l'exécution de la marche. Le manuscrit présente les mises à jour nécessaires de l'architecture de contrôle ainsi que le développement de l'algorithme de gestion de la sécurité. Des résultats expérimentaux sont fournis pour évaluer l'efficacité de notre approche lors d'exercices de rééducation de la marche en utilisant une trajectoire sans et une trajectoire avec déroulé du pied.

ABSTRACT

Self-balanced exoskeletons offer individuals with disabilities recovered ambulation capabilities thanks to their embedded control algorithms which realize autonomous and stable walking. This thesis is concerned with enabling patients to use these systems for rehabilitation purposes according to the concept of active gait rehabilitation. In our interpretation, this concept entails giving back some freedom of motion to the patients. This objective is conflicting with ensuring safety, which implies taking control over the motion. The main topic addressed in the thesis is the challenge of reconciling these two objectives. Our solution is twofold: empower the patient and monitor safety. Firstly, we introduce a split in the control architecture, allowing patients to actively participate in the swing leg motion. Concurrently, the control algorithm of the support leg regulates the unilateral contact forces. Secondly, we solve optimal control problems in free final time on a reduced model to generate stabilizable trajectories. The duration of these trajectories serves as a determinant for when the control system should take back control over the execution of the walking pattern. The manuscript presents the necessary control architecture updates along with the development of the safety management algorithm. Experimental results are provided to assess the effectiveness of our approach in both flat-foot and foot-rolling walking exercises.

KEYWORDS

Rehabilitation Robotics, Exoskeleton, Stabilization, Humanoid and Bipedal Locomotion, Optimal control