Théa Lancien
email: thea.lancien@centralesupelec.fr

Kevin Prieur

Daniel Durox

Sébastien Candel

Ronan Vicquelin

Leading point behavior during the ignition of an annular combustor with liquid n-heptane injectors

Keywords: Two-phase flow, Large-Eddy Simulation, Ignition, Eulerian-Eulerian approach, Annular combustor

Experimental and numerical investigations of ignition in combustors with multiple burners have recently emerged and have provided new insights on the last phase of ignition in gas turbine-like annular geometries where the flame propagates from burner to burner. Previous comparisons between calculations and experiments of light-round in a laboratory scale annular com-

Leading point behavior during the ignition of an annular combustor with liquid n-heptane injectors Paper length: 8 pages, Equivalent word count: 6200 words Supplemental material: A list of supplemental material is included on the last page

All figures in the present article can be printed in black and white, however, if color printing is judged necessary, the author agree to pay color reproduction charges if applicable.

bustion chamber have demonstrated the ability of large-eddy simulation to predict such processes for perfectly premixed conditions and, more recently, for n-heptane spray injection. The present analysis focuses on two additional operating points with liquid n-heptane sprays and the turbulent flame propagation in the two-phase mixture is examined through the behavior of its leading points. The validation of the light-round process is characterized in terms of ignition delays. The detailed analysis of the propagation through the definition of a leading point enables to highlight some key phenomena responsible for the flame behavior, such as the influence of the liquid droplet

Introduction

One of the critical design issues of any aero-engine combustor is that of a smooth, safe and reliable ignition over a wide operability range. In the case of a multiple injector combustor, a successful ignition can typically be decomposed into three main stages [START_REF] Lefebvre | Gas Turbine Combustion[END_REF]: The first two stages correspond to the flame kernel initiation and to the kernel's expansion and propagation to stabilize a flame on a single burner. The last stage of the ignition process involves the flame propagation from each ignited burner to the next until the full system is ignited. In the case of an annular chamber, this phase is called the light-round. Compared to single burner studies [6][7][8][9][10], this last stage has been less well documented. The interest in the need to cope with burner to burner initiation has recently led to a number of studies of ignition on multiple injector systems.

Results of the first large-eddy simulation (LES) of the ignition of a helicopter combustion chamber [11] have shown the strong impact of burnt gases volumetric expansion on the flame propagation velocity, but comparisons with experiments were not reported. Flame propagation has been studied experimentally in a linear array of injectors with varying inter-injector spacing with gaseous premixed injection [12] and liquid n-heptane injection [10]. Interinjector distance was found to modify the pattern of flame motion giving rise to two propagation modes. This feature has been retrieved in LES [12] and also observed in an annular combustor [13,14]. Nonetheless, such linear configurations do not account for specific phenomena arising in annular geometries. In 2013, two experimental studies have been reported on successful light-round sequences in idealized annular systems, with premixed propane-air injection [15] and methane-air injection [16]. Both studies highlighted the influence of the mixture equivalence ratio and bulk velocity on the light-round duration. Later on, investigations have reported the effects of non-premixed injection [13], inter-injector spacing [13,14], and liquid fuel injection [17], bringing new insights on the ignition of realistic gas turbines.

The first comparisons of light-round simulations with experiments [18,19] have demonstrated the ability of LES to successfully predict the duration and flame propagation behavior in an annular fully premixed propane-air combustor and for several operating points. A first analysis of the flame propagation pattern has distinguished several stages in the light-round, thus refining the classical three phase decomposition of the process [START_REF] Lefebvre | Gas Turbine Combustion[END_REF]. Accounting for the spray of droplets presents an additional complexity, but is essential when dealing with realistic aeronautical configurations. Prieur et al. [17] carried out the first experimental visualization and characterization of ignition in a laboratory scale annular combustor fed with liquid spray. A large-eddy simulation of the corresponding light-round sequence with liquid n-heptane [20] has recently enabled the first comparison of ignition in an annular combustors involving a two-phase flow. Results have shown that the global behavior of the flame and the light-round durations are accurately predicted, validating, for this case, the employed modeling for the dispersed spray with a monodisperse Eulerian approach. Identical phases as those identified in gaseous conditions in [18] have been observed. As noted in previous studies, the effect of burnt gases volumetric expansion yields an enhanced propagating velocity due to an induced flow in the fresh stream.

The nature of the mechanisms controlling the flame propagation velocity during the light-round remains an open question. One possible explanation of the turbulent flame burning velocity might rely on the leading points of the advancing flame fronts [21][22][23][24]. Based on observations of particular flow dynamics, a leading point scenario is here considered to identify some of the main mechanisms involved in the flame propagation velocity. Inhibiting or promoting effects linked to the local curvature (and related to the Markstein length) are not considered because their order of magnitude is lower than that pertaining to the flow dynamics.

In the present work, the behavior of these leading points is examined during the ignition of an annular combustor equipped with liquid spray injectors.

The experimental and numerical configurations are described in Sec. 2. Numerical results are analyzed for three operating points: the case considered in [20] and two additional operating conditions (Sec. 3). Simulations are first compared with experimental data in terms of light-round duration and sector ignition delays. The trajectory and dynamics of the leading point of each flame are then considered. Finally, the observed interactions between the propagating flame front, its leading point and the liquid spray provides new insights on ignition in annular combustors with liquid spray injectors.

Experimental and numerical configurations

Experimental set-up

The model scale annular combustor MICCA-Spray (Fig. 1) is described in [17]. It comprises sixteen swirl spray injectors equally distributed in a circular pattern. A two-phase mixture of air and liquid n-heptane is injected in the chamber through sixteen swirler units (see the sketch in Fig. 1). The other end of the chamber is open and the exhaust gases are released into the atmosphere. As indicated in [17] and [20], the evolution of the flame in MICCA-Spray is This is achieved by keeping a constant fuel mass flow rate while varying the air flow rate.

Numerical configuration

All simulations are carried out with the LES solver AVBP [25], that solves the three dimensional compressible Navier-Stokes equations for reactive flows.

A two-step Taylor-Galerkin weighted residual central distribution scheme (TTGC), third order in time and space [26] is used for both gaseous and Eulerian liquid phases. The numerical set-up is identical to that described in [20]. The liquid phase is modeled using an Eulerian mono-disperse approach and the droplets are assumed to have a uniform temperature and their evaporation follows the Abramzon-Sirignano model [27]. The description of the liquid phase was chosen as a compromise between computational cost and accuracy and the parameters of the model, in particular the injected mono-disperse diameter, were carefully evaluated and presented in [20,[START_REF] Lancien | Int. Conf. on Multiphase Flow[END_REF].

In the case described in [20], the calculation of local physical quantities used in the evaporation law, such as viscosity, did not depend on the local gaseous composition. This has been corrected for the three cases, including case φ = 0.89, which has been simulated anew, with an improvement of the results. Combustion of the gaseous n-heptane is represented by a two-step reduced chemical mechanism obtained in [START_REF] Paulhiac | Modélisation de la combustion d'un spray dans un bruleur aéronautique[END_REF][START_REF] Shum-Kivan | [END_REF]. The unresolved flame structures and the interactions between flame and turbulence are described using TFLES model [31][32][START_REF] Legier | Proc. Sum. Prog[END_REF] and the efficiency function derived in [START_REF] Charlette | [END_REF]. The boundary conditions are imposed using characteristic boundary conditions [35,36]. In the present study, relight conditions of the combustor are considered: in the experiment, the combustor walls are hot and have reached thermal equilibrium after a period of operation of the system. Under these conditions, one may assume that heat losses from burnt gases are signifi-cantly decreased and that adiabatic boundary conditions are applicable at the chamber walls. The parameters for the liquid injection such as droplet diameter and liquid dispersed phase density have been determined through a separate study [20,[START_REF] Lancien | Int. Conf. on Multiphase Flow[END_REF]. Particular care is taken in obtaining the initial solutions for the three operating points [20]. The computation requires a mesh size of 288 million elements and the simulation of one light-round sequence, including cold flow convergence, takes 3.5 million CPU hours on 6144 processors of the French national cluster CURIE. The ignition is triggered experimentally by a spark plug, located close to the injector in sector S 0 , that generates a propagating kernel. This complex phenomenon is simulated by defining a sphere of hot burnt products that represents the initial hot spot. the predicted duration falls in the experimental scatter of points, and φ = 0.8, where the duration is overestimated by 5%, which remains acceptable.

However, the relative error of 30% in the φ = 1.0 simulation indicates that this calculation is less adequate, probably due to the locally rich conditions at the leading point which are discussed later.

Another assessment of the simulation consists in comparing durations between the ignition of two consecutive burners. The PM measurements give access to the experimental evolution of the OH* light emissions in each sector on the H-side of the chamber. The instant the emission is maximum in one sector is taken as its ignition time, allowing the calculation of the delays between consecutive sectors. For each operating point, two or three measurements are plotted in Fig. 4 as symbols, each symbol standing for one set of measurements, highlighting some variability in the experimental results. In the simulations, the integrated heat release inside each sector is calculated and, similarly to the experiment, the time when a maximum is reached is taken as the sector ignition time. Values are reported in Fig. 4, for both H+ (thin solid lines) and H-(thick solid lines). As was done for the global light-round duration, the numerical and experimental results are synchronized with respect to the ignition of the first sector (S 1 or S -1 ), that sets the origin of times. Cases φ = 0.8 (Fig. 4 left) and φ = 0.89 (Fig. 4 center)

show an excellent agreement between experiment and simulation. Both H+ and H-evolutions are close to that of the experimental sequence, indicating that the accuracy observed in the global light-round duration is linked to the ability of the simulation to capture the flame propagation velocity. On the H+ side of case φ = 0.8, the plot shows that the flame front slows down between sectors S 6 and S 7 , indicating that the 5% difference observed in the light-round duration is due to a deceleration of one of the flame fronts towards the end of the sequence, when the two fronts collide head on and the flow is notably modified by the flame. Finally, the plot on the right in Fig. 4 allows to identify the source of error for the case φ = 1.0 for both sides of the chamber. The simulation systematically overestimates the ignition delay between sectors, leading to a gradually increasing error in ignition time, which explains the larger error in the light-round duration. Contrarily to the previous case, the error does not arise towards the end of the sequence, but is formed by an accumulation throughout the simulation, indicating that the propagation mechanisms are less well captured under these conditions.

Leading points dynamics

Similarly to the observations detailed in [20], the phases identified in [18,19] can be retrieved in the three cases simulated in this study (not shown). Once the kernel has propagated to the vicinity of the injector (phase I) in S 0 , an arch-like flame is formed that propagated towards sectors S 1 and S -1 and the chamber exit, corresponding to phase II. Upon reaching the exit, the flame separates into two fronts, thus going into phase III, where the propagation is driven from the bottom of the chamber, until the fronts merge in the other side of the annulus (phase IV). These phases are observed in both experiment and simulation, for all three cases, showing that the numerical set-up is able to retrieve the behavior of the flame, even without the correct propagation speed in case φ = 1.0.

For a closer study of the flame and an understanding of its propagation mechanisms, it is interesting to track its leading point (LP). For each half chamber, the leading point is defined as the most forward point in the azimuthal direction on the c = 0.9 iso-surface, c being the progress variable defined as c = (Y CO 2 + Y CO )/(Y eq CO 2 + Y eq CO ), which ensures its monotony through the flame front. The value of c = 0.9 was chosen because it is close to the peak of heat release in the flame in the Arrhenius law, but the same post-processing has been performed for the leading point defined by tracking c = 0.1 and c = 0.5 iso-surfaces. This yields qualitatively similar results leading to the same conclusions. For example, the effect of the progress variable threshold on the leading point trajectory is given in supplemental material. The shape of the flame fronts is strongly influenced by the local flow field, so there can be a competition between several flame elements. Due to its definition, the LP's trajectory features some discontinuities in its evolution and characteristics. Instants where the LP crosses each injector azimuthal position are plotted in Fig. 4 as dashed lines (thin dashed lines for H+ and thick dashed lines for H-), synchronized at sector S 1 or S -1 as with the other data sets. The consistency between the data extracted from the LP and from the sector-integrated heat release, which is a more global measure, confirms that one possible scenario for the flame propagation involves the leading-point entraining the turbulent flame brush. Some insight can be gained by observing the trajectory of the LP and examining the interaction of the flame tip with the aerodynamic flow. Figure S2 shows the axial position of the LP for both halves of the chamber and for all three cases. For each case and during the whole propagation, the LP is located below 50 mm, which corresponds to the bottom 25% of the chamber.

The LP in H-(dashed lines) exhibits a relatively smoother behavior, with fewer discontinuities. On the H+ side however, some common features can be identified (less pronounced for φ = 0.89). During the first third of the angular propagation, which corresponds to phase II, discontinuities indicate that the LP is lifted up due to the arch-like flame. The distance of the LP with respect to the chamber backplane gradually increases until it is suddenly brought back. This trend is repeated several times, until some point after the half of the trajectory, which varies depending on the case, where the LP definitely drops down near the backplane. This behavior suggests that favorable conditions exist near the chamber inlet, probably due to the presence of the swirling flow generated by the injectors in combination with richer pockets where droplets are trapped by the outer recirculation zones. This can be evidenced by looking at the value of the local gaseous equivalence ratio at the leading point, plotted in Fig. 6, for the three cases. The equivalence ratio is computed based on the gaseous species atoms balance. A horizontal black dashed line highlights the injected global equivalence ratio φ. In all cases and during the whole propagation, the LP is systematically in a mixture that is richer than φ. Fewer discontinuities are observed here than in the LP position, indicating that the equivalence ratio does not vary much in the LP propagation region and therefore that the competing flame tips burn at similar φ. When the flame fronts are located in the last third of the chamber, some large swings occur in the local equivalence ratio evolution that can be surprising. The tangential positions where these events take place correspond to the positions of the injectors, suggesting that in the presence of the highly turbulent sprays, the flame fronts encounter large φ variations, due in particular to the presence of pockets of liquid fuel. The presence of the series of swirling jets can also be an explanation for this concentration of fuel vapor at the bottom of the chamber. This is further discussed in the next section.

Flame-spray interactions

As noted in other similar studies ([11, 15, 18, 19]), the volumetric expansion of the burnt gases in the semi-confined chamber has a notable impact on the fresh stream before the flame. It is therefore interesting to analyze its interaction with the liquid spray. This can be observed by considering a cylindrical cut inside the chamber and unfolding this surface. In this representation the spark plug corresponds to the left and right boundaries of the unfolded cylindrical surface and the flames propagate towards the center of the image. Figure 7 shows the gaseous tangential velocity u θ field (top image) and the liquid mass density α l ρ l (bottom image), for the case φ = 0.8 at an instant chosen to observe phase III, where two propagating fronts are present. Similar behaviors are observed in the other operating cases. It appears clearly that the propagation of the flames generates an azimuthal flow that pushes the fresh gases several sectors ahead of the flame. During the light-round, the leading point evolves in this region of large azimuthal velocity which dominates the LP velocity (see Fig. 8). The bottom image in Fig. 7 indicates that the flow induced ahead of the flame perturbs the liquid phase distribution, evidenced by the sprays injected by the un-ignited injectors, which are inclined away from the flame, a feature that is observed experimentally. Some pockets of droplets are also created and pushed away from the flame fronts, forming patterns similar to those found with jets in cross-flows (see insert on the right in Fig. 7). This phenomenon is due to the strong azimuthal velocity which forms wakes on the downstream side of the swirling streams established by the different injectors.

The induced flow impact on the liquid phase can be further evidenced by observing the temporal evolution of some quantities at selected points in the chamber. Figure 9 displays the evolution of α l ρ l and the local gaseous equivalence ratio φ g , for the three cases, at the point [x = 2.5 mm, r = 0.175 m, θ = 9π/8], which corresponds to a point located near the chamber back-plane at the frontier between S 4 and S 5 , on the H+ side. This point is chosen to be near the leading point trajectory. To highlight the crossing of the flame, the progress variable c is also plotted (dotted line). In all cases, the same trend is observed; although a steady state regime had been reached before the ignition, the local liquid mass density begins to increase approximately 10 ms before the crossing of the flame (evidenced by the increase in the progress variable). This shows that an accumulation of liquid is created by the flow generated ahead of the flame fronts, so that the local conditions seen by the flame differ from those corresponding to the steady-state cold two-phase flow field. As a consequence, the local gaseous equivalence ratio is increased close to the flame front, when this accumulation of droplets evaporates in the vicinity of the flame. This increase is expected to be even larger for fuels that are less volatile than n-heptane. It appears that in the three cases under consideration, the modification of the initial heterogeneous field of equivalence ratio by the flame propagation is not negligible and systematically yields an equivalence ratio at the LP that is larger than the global value.

The larger error of the simulated case φ = 1.0 shown in Fig. 4 could be due to a more pronounced heterogeneity with rich conditions prevailing in the leading point region, a situation that is less well tackled by the combustion model (see Fig. 6).

Conclusion

Large eddy simulations of the light-round sequence in a laboratory scale annular combustor have been carried out with liquid spray injection using an Eulerian framework and a mono-disperse representation of the spray. Three different cases with varying global equivalence ratios are compared to experimental data. Cases φ = 0.8 and φ = 0.89 indicate that the simulation is able to accurately predict the light-round duration as well as the ignition delay between the injectors, while the case φ = 1.0 is less accurately described with an underestimation of the absolute flame velocity. It is found that the volumetric expansion resulting from chemical conversion in the flame induces an azimuthal velocity in the fresh stream at a distance of several sectors ahead of the flame. This induced flow modifies conditions in the fresh stream. A detailed analysis of local dynamics is carried out through the definition of a leading point for each flame front. The three cases investigated present similar preferential leading point trajectories in the combustor. Analysis of local values at the leading point indicates that the flow induced by the flame interacts with the liquid phase, modifying the liquid volume fraction in the wakes of the swirling jets established by the various injectors. This gives rise to a heterogeneity in composition and higher equivalence ratios in the fresh stream in front of the flames. These two features influence the leading point motion and trajectory in the chamber and affect the absolute flame velocity. Future investigations will consider the impact of more accurate spray modeling taking into account the polydispersity of the liquid phase. It is also intended to change the turbulent combustion model to improve predictions under heterogeneous conditions and high equivalence ratio values. Further analysis will be necessary to assess the leading-point and mean flame brush scenarios and their probable interplay in the flame propagation process.

T. Lancien a, 1

 1 , K. Prieur a,b , D. Durox a , S. Candel a , R. Vicquelin a a Laboratoire EM2C, CNRS, CentraleSupélec, Université Paris-Saclay, 3 rue Joliot Curie 91192 Gif Sur Yvette cedex, France b Safran Tech, E&P, Chateaufort 78772 Magny-Les-Hameaux, France Colloquium: Gas Turbine and Rocket Engine Combustion Word Count: Method 2, counting from a two-column formatted version of the paper.

Figure 1 :

 1 Figure 1: Direct view of the MICCA-Spray combustion chamber [17]. The sketch at the bottom right represents a cut in one swirler unit showing the gaseous (G arrows) and liquid (L arrow) injections.

Figure 2 :

 2 Figure 2: Schematic top view of the MICCA-Spray backplane showing the swirlers positions, photo-multipliers (PM), spark plug and swirl direction.

  recorded with a high-speed camera equipped with a CH* filter featuring a bandwidth of 420 -440 nm and photo-multipliers (PM) that are sensitive to the light emission of the OH* chemical radical. In the present study, three operating points are considered, for three different global equivalence ratios at the same nominal power P = 79.3 kW: φ = 0.8, φ = 0.89 and φ = 1.0.

3 .

 3 Results and discussions3.1. Light-round durations and sector ignition delaysCalculated light-round durations are compared with those of the experimental sequences. The duration is defined as the delay between the ignition of the S 1 sector and the merging of the fronts in S 8 . Such a definition allows to focus on the light-round sequence, without taking into account the first instants of the simulation where the physics of the kernel induces some variability in its initial growth. Experimentally, these instants are determined using the recorded images of the flame evolution. Some uncertainty of a few frames arises during the determination of sector S 1 ignition time which, through the camera frame rate, leads to uncertainties on the light-round duration.The experimental and numerical values for the light-round durations are shown in Fig.3. Very good results are achieved for cases φ = 0.89, for which

Figure 3 :

 3 Figure 3: Light-round durations as a function of the global equivalence ratio. Diamond symbols stand for experimental results. The durations predicted by the present simulations are represented by red circles. The green square indicates the light-round from [20] before fixing the evaporation model. Estimated measurement uncertainties are also displayed in this graph.

Figure 4 :

 4 Figure 4: Ignition delay times for each sector in the experiments (in symbols, each set stands for a different experimental sequence) and in the simulations (in the form of lines). Data in the left, middle and right plots respectively correspond to φ = 0.8, are φ = 0.89 and φ = 1.0. Numerical data extracted from the heat release in H+ (thin solid lines) and in H-(thick solid lines), extracted from the leading point position in H+: (thin dashed lines) and in H-(thick dashed lines).

Figure 5 :

 5 Figure 5: Leading points height against their azimuthal position for H+ (solid lines) and H-(dashed lines) and all three cases: φ = 0.8 at the top, φ = 0.89 in the middle and φ = 1.0 at the bottom.

Figure 6 :

 6 Figure 6: Time evolution of the equivalence ratio at the leading point in H+ (solid lines) and H-(dashed lines) and for φ = 0.8 (top), φ = 0.89 (middle) and φ = 1.0 (bottom).The horizontal black dashed lines mark the global equivalence ratio for each case.

Figure 7 :

 7 Figure 7: Top: visualization of the tangential gaseous velocity on a cylinder of radius r = 0.175 m unfolded on a plane surface for the case φ = 0.80 at instant t = 17 ms.Bottom: visualization of the liquid mas density α l ρ l on the same unfolded plane. The radius was chosen so that the cut would cross the center of each injector. The lateral sides of the unfolded cylinder correspond to the location of the first ignited injector. Tangential velocities are counted positive from left to right. Two iso-lines of the progress variable c = 0.1 and c = 0.9 (in black) highlight the position of the flame fronts. The insert on the right represents the impact of the flame progress (red iso-surface) on the liquid spray, evidenced by a blue iso-surface at α l ρ l = 0.01 kg.m -3 .

Figure 8 :

 8 Figure 8: Zoom on the flame front leading point (LP) area on an unfolded cylindrical cut at r = r LP , colored by the tangential velocity. The flame is evidenced by white iso-lines of the progress variable c = 0.1 and c = 0.9. The LP is highlighted by the black circle.

Figure 9 :

 9 Figure 9: Temporal evolution of the liquid mass density (solid lines, right axis) and the gaseous equivalence ratio (dashed lines, left axis) for the three cases φ = 0.8 (top), φ = 0.89 (middle) and φ = 1.0 (bottom) at point [x = 2.5 mm, r = 0.175 m, θ = 9π/16]. The evolution of the progress variable c is added (dotted lines, left axis) to highlight the crossing of the flame.

Acknowledgments

This work was granted access to the HPC resources of TGCC and IDRIS under allocations GENCI 2017-A0022B10118 and PRACE n14-2016153551. This work was supported by grant ANR14-CE23-0009-01 of the French Agence Nationale de la Recherche in the TIMBER project.

List of supplemental material

Figure S2: Leading points trajectories projected on the chamber backplane. Top: the blue, red and green lines respectively stand for φ = 0.8, φ = 0.89 and φ = 1.0. The LP is defined on the c = 0.9 iso-surfaces. Bottom: Trajectories for case Φ = 0.8 for the LP defined on three different iso-surfaces, c = 0.1 (blue), c = 0.5 (red) and c = 0.9 (green).


Ines Jammeli
email: jammeli.ines@gmail.com

Senior Member, IEEE Ahmed Chemori
email: ahmed.chemori@lirmm.fr

Huiseok Moon

Salwa Elloumi
email: salwa.elloumi@gmail.com

Senior Member, IEEE Samer Mohammed
email: samer.mohammed@u-pec.fr

An Assistive Explicit Model Predictive Control Framework for a Knee Rehabilitation Exoskeleton

Keywords: Assistive Robotics, Knee joint powered orthosis, Model predictive control, Rehabilitation

This study focuses on the control of an actuated knee joint orthosis. The proposed solution is a novel model predictive control framework dedicated to assistive and rehabilitation purposes. This framework includes (i) an exact inputto-state feedback linearization, (ii) a model predictive controller (MPC or EMPC), considering input/output constraints, (iii) a least-squares dynamic parameters identification, (iv) a nonlinear disturbance observer for the estimation of the wearer's torque, (v) a Lyapunov-based stability analysis of the resulting closedloop system, and (vi) a reference trajectory generator. The proposed framework has been validated via real-time experiments performed on three healthy subjects wearing the knee joint orthosis. Various experimental scenarios have been considered, including assistive and resistive rehabilitation tasks in a sitting position and walking with normal/abnormal gait patterns. The obtained results indicate the efficiency of the proposed predictive controllers with respect to a conventional PID controller in terms of tracking performance, required torque, and wearer comfort.

I. INTRODUCTION

W ITH recent medical advances, life expectancy increases steadily. According to the World Health Organization, the world population aged 60 years and older is expected to reach a total of 2 billion by 2050; however, it was approximately 900 million in 2015 [START_REF]Ageing and health[END_REF]. This population is exposed to health risks caused by weakened muscle strength, which hinders their ability to walk as frequently as normal and adversely affects their walking stability, thereby making them dependent on others. Consequently, the aging of the population and the physical deterioration of the elderly have become a global socioeconomic problem [START_REF] Chen | Recent developments and challenges of lower extremity exoskeletons[END_REF]. This issue calls for considerable attention on how to assist this population, as well as people with lower and/or upper-limb pathologies in their daily life, especially regarding mobility and autonomy [START_REF] Mohammed | Towards Intelligent Lower Limb Wearable Robots: Challenges and Perspectives -State of the Art[END_REF].

One of the best remedies for reduced mobility is rehabilitation. Conventional intensive therapies are usually adopted in clinical centers to help people recover their voluntary movements. However, the problem with these therapies is that they are only effective when they are intensive [START_REF] Sherwani | RISE-based adaptive control for EICoSI exoskeleton to assist knee joint mobility[END_REF]. Furthermore, repetition is a key element in this case, which allows the brain to reprogram the motion sequence. Therefore, this process is time and capital intensive, and requires the strength of both the patient and therapist. However, intensive long-term rehabilitation therapy is not always an option owing to its expensive cost, and the insufficiency of qualified staff. To reduce the burden on care services, several initiatives have been set up to promote assistive technologies, such as the promising technologies based on wearable robots. These devices are mechatronic systems, equipped with sensors and actuators, and embodied by the human upper and/or lowerlimbs, which provide the following functions: (i) augmenting physical human capabilities at upper/lower-limbs, (ii) assisting people with reduced mobility for achieving daily living activities, and (iii) automating the rehabilitation of human joints and muscles to recover and improve the control of the wearer's limbs [START_REF] Mohammed | Lower-Limb Movement Assistance through Wearable Robots: State of the Art and Challenges[END_REF]. This considerable focus on wearable robots can be explained by their ability to reproduce repetitive tasks that require strength and robustness. They can autonomously perform these tasks faster than therapists, with a better level of accuracy, without getting tired, and without requiring a third party. Hence, they may promote a reduction in patient fatigue. They enable long training sessions with optimal consistency, as well as measurements for the user to track the desired gait patterns, which may help to accelerate the rehabilitation process [START_REF] Chen | Knee exoskeletons for gait rehabilitation and human performance augmentation: A state-of-theart[END_REF]. Sitting and standing up movements are essential for most human activities. Given their importance, several wearable robots have been developed to help weakened people perform these movements. Shepherd et al. [START_REF] Shepherd | Design and Validation of a Torque-Controllable Knee Exoskeleton for Sit-to-Stand Assistance[END_REF] designed a knee exoskeleton using a torque-controllable series elastic actuator capable of providing the torques and speeds required for sitto-stand. In [START_REF] Huo | Active Impedance Control of a Lower Limb Exoskeleton to Assist Sit-to-Stand Movement[END_REF], the authors proposed an intention-based active impedance control (AIC) strategy applied to a lower-limb exoskeleton for wearers with a partial loss of their lower-limb muscular strength. In [START_REF] Junius | Mechatronic Design of a Sit-to-Stance Exoskeleton[END_REF], the authors designed an exoskeleton to assist the wearer during sit to stand tasks by delivering assistance-as-needed to patients or elderly people with muscle weakness. In [START_REF] Kamalia | Trajectory generation and control of a knee exoskeleton based on dynamic movement primitives for sit-to-stand assistance[END_REF], a trajectory generator and an impedance controller were designed to calculate the required torques for a knee exoskeleton to assist the wearer in moving along a predicted trajectory while allowing the wearer to partly control the speed of the sit-to-stand movement. If sitting and standing up movements are essential for most human activities, walking is even more essential, and it plays a key role in everyday life. Losing the ability to walk presents a real challenge, which makes the patient dependent on others. To address this challenge and enable people with reduced mobility to recover their independence, significant attention has been given to exoskeletons. To assist people while walking, Wang et al. [START_REF] Wang | Comfort-Centered Design of a Lightweight and Backdrivable Knee Exoskeleton[END_REF] developed a highly back-drivable, partial-assist knee orthosis designed to improve excessive shear forces, reduce joint misalignment, and minimize the distal weight. In [START_REF] Zhu | Design and Validation of a Partial-Assist Knee Orthosis with Compact, Backdrivable Actuation[END_REF], the authors presented the design of a partial-assist knee orthosis by improving its back-drivability, weight, and size. They reduced the motor winding temperature and enabled higher output torques by adopting an encapsulation technology. In [START_REF] Lv | On the Design and Control of Highly Backdrivable Lower-Limb Exoskeletons[END_REF], the authors presented a task-invariant and usercooperative approach, known as energy shaping, to provide assistance by augmenting body energetics. Although some of the aforementioned controllers provide fairly satisfactory performances, none of them can provide optimal control while considering the system constraints. In this study, we focus on the control of exoskeleton intelligently communicating and sensitive to intention (EICOSI) system, an active lower-limb orthosis acting on the knee joint level, developed for assistive and rehabilitation purposes. This system has a simple structure, is easy to don and doff, and practical for people suffering from knee joint impairments. Various control approaches have been proposed for EICOSI such as PID control, adaptive control [START_REF] Mohammed | Adaptive Control of a Human-Driven Knee Joint Orthosis[END_REF] [START_REF] Rifai | Augmented L1 Adaptive Control of an Actuated Knee Joint Exoskeleton: From Design to Real-Time Experiments[END_REF], and high-order sliding mode control [START_REF] Mohammed | Towards Intelligent Lower Limb Wearable Robots: Challenges and Perspectives -State of the Art[END_REF].

In this study, a novel model predictive control framework (illustrated in Fig. 1), dedicated to assistive and rehabilitation purposes, is proposed. This framework includes (i) an exact input-to-state feedback linearization, (ii) a model predictive controller (MPC) considering input/output constraints, (iii) a least-squares dynamic parameters identification, (iv) a nonlinear disturbance observer (NDO) for the estimation of the wearer's torque, (v) a Lyapunov-based stability analysis of the resulting closed-loop system, and (vi) a reference trajectory generator. Although MPC has several advantages in comparison with conventional controllers, to the best of our knowledge, the proposed controller has not been used in an assistive context with actuated orthosis. This control allows desired performances to be expressed in a single cost function that can be easily adjusted according to each user. Moreover, the tracking error is guaranteed along with the smoothness of the orthosis behavior ensured by the minimization of the cost function in the MPC control formulation. This smoothness could have a considerable and favorable impact on the wearer's comfort. Furthermore, the tracking error shows the ability of the controller to follow a desired trajectory while taking into account external perturbations and uncertainties. Providing satisfactory tracking performances as prescribed by therapist doctors within a rehabilitation process will ensure that the rehabilitation purposes in terms of muscle reinforcement and joint range of movement enhancement are met. The MPC formulation is also characterized by its ability to implicitly consider constraints on the control inputs (torques), as well as on the states and outputs (position, velocity, and acceleration), which may ensure the safety of the wearer, and the stability of the overall control strategy. In addition, MPC has an anticipation feature, which enables it to anticipate a forthcoming change in the desired behavior of the system to be controlled. MPC can also be considered a robust control that can compensate for parametric uncertainties resulting from a poor identification process, or for external disturbances such as the exerted human torques or disturbances originating from the interaction with the ground. The proposed framework has been validated via real-time experiments performed on healthy subjects (a) in a sitting position while asked to exert an assistive/resistive effort, and (b) during a walking activity while considering the interaction with the ground, and following a reference walking pattern generated by the proposed trajectory generator. The exerted human torque and the effect of the ground interaction during the stance phase were estimated using the proposed NDO. The rest of the paper is organized as follows: A description of the EICOSI active lower-limb orthosis is presented in section II. The proposed MPC approaches are discussed in section III. Section IV is dedicated to the proposed NDO, as well as the proposed walking trajectory generator. The resulting closedloop system stability analysis is discussed in section V. Section VI gathers all the common methods between the conducted experiments and illustrates the experimental protocol. The obtained experimental results are presented and discussed in section VII. Finally, section VIII provides some concluding remarks, and projections for future works.

II. EICOSI ACTIVE ORTHOSIS A. Description and modelling of EICOSI orthosis

EICOSI active orthosis is a single DoF prototype adopted in this study to assist the wearer at the knee joint level for flexion/extension movements in a sitting position and during a walking activity. Fig. 2 shows the subject (a) in a sitting rehabilitation position and (b) walking on a treadmill while wearing the actuated orthosis. EICOSI comprises two segments attached separately to the wearer's thigh and shank using appropriate braces. The shank-foot of the wearer is considered as the single rigid segment, and it freely moves around the knee joint. The orthosis is driven using an EC-4pole Maxon 200W motor controlled by an ESCON 50/8 driver. A real-time communication with FPGA card is established. To guarantee an efficient and portable system, as well as relatively high output torque, a compact transmission system is designed using a gear motor, ball screw, transmission belt, and cable drive. The reduction ratio from the motor to the joint side is 264:1. The actuation system can deliver a torque up to 18 N.m, which is estimated using a current sensor. The motor is equipped with an incremental encoder to measure its angular position, which is used to compute the knee joint angle. The angular velocity and acceleration are obtained by numerical derivatives (first and second-order respectively) of the angular position. The thigh angle is measured using an IMU (inertial measurement unit: MTw Awinda by Xsens © ) attached to the thigh brace, rather than the user's shank, to obtain less noise from ground interaction, and facilitate the measurements at the sitting position. The experiments were performed using the LabView software. This active orthosis enables its user to perform flexion/extension movements with a relative angle ranging from 0 • to 135 • , where the latter corresponds to the maximum knee joint-orthosis flexion and the former to the full knee joint-orthosis extension. For dynamic modelling, the application of Newton's second law gives:

J θ = i τ i = τ e + τ h -τ g -τ f -τ l ,
where τ e , τ h , τ g , and τ f represent the exoskeleton torque, human torque, gravitational toque, and friction torque, respectively. The latter primarily comprises static and viscous frictions, expressed by Asign( θ(t)) + B( θ(t)), where A and B respectively represent the coefficients of static and viscous friction torques of the entire system, including the human leg and the exoskeleton. It is noteworthy that the adopted friction model is advantageous because it is linear in its parameters A and B. τ l represents the lumped torque triggered by the other body segments and the eventual external forces. J is the inertia of the human lower limb wearing the active orthosis. The wearer's shank and the embodied actuated orthosis system can be dynamically modeled by the following second-order equation:

J θ = -T g sin(θ -θ t ) -Asign( θ) -B( θ) + τ e + τ h -τ l (1)
where θ, θ, and θ represent the angular position, velocity, and acceleration of the entire lower part system (including the wearer's shank and the active orthosis lower part), respectively. θ t , and θt are the angular position and velocity of the upper part (including the wearer's thigh and the upper part of the orthosis [START_REF] Mohammed | Nested saturation based control of an actuated knee joint orthosis[END_REF]), respectively. In this work, we are interested in the application of the proposed MPC-based framework on EICOSI orthosis, both in a sitting rehabilitation position and for a walking activity. The human gait cycle can be categorized into two different phases as illustrated in Fig. 3. The stance phase, accounting for 60% of the cycle, during which the foot is placed on the ground, and the swing phase, accounting for 40% of the cycle, during which the foot is no longer in contact with the ground. In the swing phase, the lumped torque τ l can be expressed by the following equation [START_REF] Huo | Impedance Reduction Control of a Knee Joint Human-Exoskeleton System[END_REF]:

τ l = m s d s l t ( θt (t)cos(θ(t)) -θ2 t (t)sin(θ(t))) (2)
where m s , l t , and d s represent the lumped mass of the human shank and the corresponding orthosis segment, human thigh length, and length of the human shank part from the knee joint center to the center of mass of the shank-foot, respectively. In addition, θt and θt are the velocity and acceleration of the upper part, respectively. For the stance phase, the contact forces with the ground are considered external perturbations.

Because the studied active orthosis has only one degree of freedom, and does not include an ankle actuator, the ground does not interact directly with the exoskeleton; hence, the reaction forces could be considered as external disturbances to the system. In this case, the effect of the ground contact can be considered within the torque estimated by the proposed NDO.

B. Dynamic identification of parameters

The proposed control framework (cf. Fig. 1) includes a nonlinear state feedback and an MPC, both of which are modelbased control strategies. Consequently, the dynamic model of the system is significantly beneficial in the performance of the proposed control scheme. To minimize modeling errors, an identification process is performed as follows. We ask the subject to wear the exoskeleton, stay in a sitting position, and be passive (i.e. not generating any voluntary movement) during the identification process. Because the human-orthosis is considered to be a single system, we assume that τ h , τ l , θt , and θt are zero and that θ t = π 2 . Hence, the inverse dynamic model can be represented by the following equation:

τ e = J θ -T g cos(θ) + Asign( θ) + B θ (3)
which can be reformulated as affine in the parameters:

τe(t) = -cos(θ(t)) sign( θ(t)) θ(t) θ(t)    Tg A B J    = γ T (t)P (4)
with γ T (t) = -cos(θ(t)) sign( θ(t)) θ(t) θ(t) and P = [T g A B J] T . The identification algorithm is based on the adoption of the inverse model ( 4), which is linear in the dynamic parameters of the system. It enables the estimation of these parameters, providing either a measurement or an estimation of the joint torque τ e , as well as the joint position θ, velocity θ, and acceleration θ. Let us define T = [τ e (1) τ e (2) . . . τ e (n)] T , and Υ = [γ T (1) γ T (2) . . . γ T (n)] T , with n as the number of samples. The parameters of our system, i.e T g , A, B, and J are computed by minimizing the error:

ε(t) = τ e (t) -γ T (t)P (5)
via the optimization of the following cost function:

ψ(P ) = 1 n n t=1 ε 2 (t). (6)
The optimal parameters can be obtained by zeroing the term

dψ(P )
dP . Computing P = (Υ T Υ) -1 Υ T T , facilitates the identification of the studied system parameters [START_REF] Krneta | Recursive Least Squares Method in Parameters Identification of DC Motors Models[END_REF]. Different series of various tests were performed using a basic proportional controller while applying chirp signals as inputs. For each test, a different amplitude was used. This enables us to obtain the system response to several amplitudes and frequencies. Each parameter has been computed as the average of the obtained values from different tests, and are summarized in TABLE I.

C. Input-to-state feedback linearization

The proposed MPC approach is a discrete-time controller based on a linear model of the system. However, the studied active orthosis has a continuous nonlinear model; therefore, it should be linearized and discretized to integrate it into the MPC controller. The approximated linearization of the EICOSI system (1) is not feasible owing to the existence of the term sign( θ), which gives it a highly nonlinear feature. Consequently, we propose the use of an exact input-to-state linearization approach, which can transform the nonlinear model into an equivalent linear alternative via a change of variables and a suitable nonlinear feedback. Let us consider the nonlinear state feedback:

τ e = α(X) + β(X)v (7)
By selecting the adequate control input τ e , this approach facilitates the emergence of a linear input-output map between the linearized input v and the output y. For the sitting position, we set:

   α(X) = -T g cos(θ) + Asign( θ) + B( θ) -τ h β(X) = J X = [x 1 x 2 ] T = [θ θ] T
The nonlinear state feedback τ e can then be expressed as:

τ e (t) = -T g cos(θ) + Asign( θ) + B( θ) -τ h + Jv. ( 8
)
For the standing position, we set:

   α(X) = T g sin(θ -θ t ) + Asign( θ) + B( θ) -τ h + τ l β(X) = J X = [x 1 x 2 ] T = [θ θ] T
which leads to the following linearizing control input:

τ e (t) = T g sin(θ -θ t ) + Asign( θ) + B( θ) -τ h + τ l + Jv (9)
Introducing the control input (8) (for sitting rehabilitation, with θ t = π 2 ) or ( 9) (for walking activity) into the dynamics (1) results in the following closed-loop (inner loop in the proposed framework of Fig. 1) linearized model:

   ẋ1 = x 2 ẋ2 = v y = x 1 (10)
Because the expression of the position reference trajectory θ ref is known, then the velocity and acceleration reference trajectories θref and θref , respectively can be easily computed using successive time derivatives. For the aim of stability analysis [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] of the resulting closed-loop dynamics introduced in section V hereafter, the definition of the following state errors are required:

X = X -X ref = x1 x2 = θ -θ ref θ -θref (11)
ẋ2 = θ -θref = v -v ref = ṽ with v ref = θref (12)
The obtained linear system [START_REF] Kamalia | Trajectory generation and control of a knee exoskeleton based on dynamic movement primitives for sit-to-stand assistance[END_REF] can now be discretized using Euler's method and a sampling period T s , which leads to the following discrete state-space representation:

X(k + 1) = A X(k) + Bṽ(k) ỹ(k) = C x(k) (13)
with A = 1 T s 0 1 , B = 0 T s , and C = 1 0 .

Owing to the exact feedback linearization and the system discretization, the formulation of the MPC problem can now be easily defined, as introduced in the following section.

III. PROPOSED MPC APPROACHES

A. A brief background on MPC

MPC is a strategy widely adopted in the process industry and successfully implemented in various applications. It is an advanced control scheme based on the process model, which is used to predict the system future behavior over a certain prediction horizon, such that the best control action is selected [START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF]. A major advantage of MPC is its ability to implicitly consider constraints imposed on the system inputs and/or outputs by including them in the design process [START_REF] Kohler | A nonlinear model predictive control framework using reference generic terminal ingredients[END_REF]- [START_REF] Camacho | Model Predictive Control[END_REF]. Its simple and intuitive design principle makes it a very popular choice in both the industry and in academia. The methodology of MPC is characterized by the following strategy: At each time instant k, the output of the system is measured. Based on this measurement, an optimal control sequence is computed. This sequence minimizes the cost function, which differs from one application to another, according to the user's expectations, over a finite horizon. Generally, it is assumed that the control signal varies over a control horizon with length N c (less than or equal to the length of the prediction horizon N , and remains constant over the interval [N c , N ], with N as the length of the prediction horizon. Only the first sample of this sequence is applied, and the same steps are repeated at time k + 1 by shifting the prediction horizon. This update allows the consideration of uncertainties and/or disturbances that may occur while controlling the system [START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF]. The main limitation of this scheme is that, by considering the constraints imposed by the system, the MPC becomes relatively greedy in terms of computing time and memory space [START_REF] Li | Real-time Optimal Energy Management Strategy for a Dual-Mode Power-Split Hybrid Electric Vehicle Based on an Explicit Model Predictive Control Algorithm[END_REF].

B. Application of MPC to EICOSI orthosis

Here, we propose the control of EICOSI orthosis using an MPC (block (ii) in the proposed framework of Fig. 1), to track the reference trajectory while smoothing the orthosis behavior. After identifying the system parameters, linearizing, and discretizing its dynamics, we can apply the proposed MPC approach to the studied lower-limb active orthosis as explained hereafter. The associated optimization problem can be expressed:

V (k) = min( 1 2 ṽT →k H ṽ→k + X(k) T F ṽ→k ) + 1 2 X(k) T Y X(k)                s.t θmin ≤ θ(k + i) ≤ θmax for i = 1..Nc θmin ≤ θ(k + i) ≤ θmax for i = 1..Nc ṽmin ≤ ṽ(k + i) ≤ ṽmax for i = 1..Nc X(0) = X(k) X(k + i + 1) = A X(k + i) + Bṽ(k + i) for k ≥ 0 ỹ(k + i) = C X(k + i) for k ≥ 0 (14)
Iterating the state-space (13) over the prediction horizon gives:

     X(k+1) X(k+2)
. . .

X(k+N )      =      A A 2 . . . A N      X(k)+       B 0 • • • 0 AB B . . . . . . . . . . . . . . . 0 A N-1 B • • • AB B           ṽ(k) ṽ(k+1) . . . ṽ(k+N -1)     (15)
which can be rewritten in a more compact form as:

X→k+1 = P X X(k) + H X ṽ→k (16)
The first constraint in ( 14) is imposed by the mechanical stops of the exoskeleton. The second constraint is imposed by the physical limit of the actuator in terms of maximum speed. The angular speed of the knee joint is related to the motor angular speed. Accordingly, knowing the range of the motor, the equivalent range of the knee joint speed can be computed. Knowing that the reference trajectory is priori known and bounded in terms of θ ref , θref , and θref , the third constraint is discussed hereafter. Because the proposed control framework includes two control loops (an inner loop based on a nonlinear state feedback control law and an outer loop based on the proposed MPC/EMPC controller), we designed our EMPC controller considering the constraints on the control input τ e . In this study, the constraints on the linearized input have been computed based on those imposed on τ e . Using the relationship between these two variables expressed in [START_REF] Shepherd | Design and Validation of a Torque-Controllable Knee Exoskeleton for Sit-to-Stand Assistance[END_REF], to maintain τ e within its admissible values, we have to ensure that the right-hand side of this equation satisfies these bounds. By imposing constraints on the angular position and velocity, and because τ h and τ l are finite-valued torques, the constraints on the linearized input, which forces τ e not to exceed its admissible values, can be deduced. Regarding the effect of the choice of these constraints on the cost function, it is noteworthy that altering each of these limits will change the size of the search space of the optimization problem. Accordingly, small limits will reduce the search space and may lead to suboptimal solutions, or even in extreme cases, to an unfeasible problem. However, dealing with an important number of varying parameters may cause discomfort to the wearer or a physical saturation. In summary, tuning is performed based on a compromise between the closed-loop desired performance and the physical limits of the system (in terms of joint position, speed, acceleration, and motor torque). In the proposed cost function V (k) [START_REF] Huo | Impedance Reduction Control of a Knee Joint Human-Exoskeleton System[END_REF], all the target performances are gathered and expressed mathematically. Specifically in this study, we attempt to reduce the tracking position error, as well as smoothen the EICOSI response by reducing the chattering phenomenon and the tracking velocity error. This approach was also selected because it helps to ensure guaranteeing an optimal trajectory tracking, even in the presence of parametric uncertainties originating from the identification process, and in the presence of external disturbances emerging from the interaction between the wearer and/or the ground. The weighting matrices Q and R in [START_REF] Huo | Impedance Reduction Control of a Knee Joint Human-Exoskeleton System[END_REF], are related to the state tracking errors and the acceleration, respectively while Q f is related to the terminal constraint. The cost function [START_REF] Huo | Impedance Reduction Control of a Knee Joint Human-Exoskeleton System[END_REF] was selected to meet our expectations in terms of tracking performance, as well as realize the closed-loop stability of the studied system. Initially, we set some conditions on the cost function tuning parameters (Please refer to Section V for more information), then we selected those that allow the desired performance to be achieved while satisfying these conditions:

V (k) = N -1 i=0 ( X(k + i) 2 Q + ṽ(k + i) 2 R ) + X(k + N ) T Q f (N ) X(k + N ) (17)
This cost function can be represented in a compact form as:

V (k) = N -1 i=0 l( X(k+i), ṽ(k+i)) + F ( X(k+N )) (18)
By substituting X(k+i) = A i X(k) + i-1 j=0 A j Bṽ k+i-1-j into (18) and ( 14), we obtain the following equation :

         V (k) = min( 1 2 ṽT →k H ṽ→k + X(k) T F ṽ→k ) + 1 2 X(k) T Y X(k) s.t Gṽ →k ≤ W +E X(k) (19)
where H, F, G, E, Y and W are easily obtained from Q, R, Q f and ( 18) and ( 14). The constraints imposed by the system are omitted while calculating the MPC solution. This facilitates a significant reduction in computation time, and thus enables its real-time validation even with a small sampling period. In this case, the optimal control sequence is obtained by zeroing grad ṽ→k V (k).

C. Background on Explicit Model Predictive Control (EMPC)

In the second part of this study, the constraints imposed by the system are considered while maintaining the same small sampling period. Even if efficient QP solvers based on interior point and active set methods exist, the computation of the control action ṽ(t) requires a significant online effort that makes MPC not applicable to fast systems with small sampling time, which is similar to the case of our study. Knowing that the sampling period is too small to enable the application of constrained MPC in real-time, we propose computing the possible control values offline using an EMPC method. This method is based on multi-parametric quadratic programming, which pre-solves the QP offline for the entire set of states. By considering X(t) as the vector of parameters, the EMPC moves the computation time offline by solving the optimization problem for all possible values of X(t) of interest, and making this dependence explicit. It converts the MPC into a continuous and piecewise-affine function of the parameter vector; thus enabling its application in realtime for fast systems. Generally, instead of computing the control online by minimizing the cost function, the problem is converted to a lookup table containing the optimal values of the gains that allow the MPC control law to be obtained explicitly as a function of the states, which considerably reduces the computation time [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF].

D. Application of EMPC to EICOSI orthosis

Here, we consider the entire optimization problem as expressed in [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] (while considering the constraints imposed by the system). Let us consider X = X X ref τ e T as the parameter vector of the problem in [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], which will be considered as a multi-parametric quadratic programming problem. Based on the Karush-Kuhn-Tucker (KKT) conditions for optimality, the authors of [START_REF] Bardaro | MPC-based control architecture of an autonomous wheelchair for indoor environments[END_REF] proved that the optimizer function ṽ * →k (x) is piecewise affine and continuous over the set of feasible parameters. They also verified that the corresponding optimal value V (k) * is continuous, convex, and piecewise quadratic. The algorithm of EMPC begins by selecting an arbitrary starting parameter vector X 0 such as the origin. It then solves the QP problem to obtain the optimal solution ṽ(X 0 ) * and identify the subset of all active constraints represented by Gṽ( X) = S X + W , as well as the subset of inactive constraints represented by Ĝṽ( X) ≤ Ŝ X + Ŵ .

IV. NDO AND REFERENCE TRAJECTORY GENERATOR A. Application of NDO on the studied system

In typical rehabilitation scenarios, when the wearer is asked to develop a muscle activation, the human torque τ h is considered an external disturbance. In this case, we propose an NDO, the block (iv) in the framework of Fig. 1, to estimate its values [START_REF] Mohammed | Nonlinear Disturbance Observer Based Sliding Mode Control of a Human-Driven Knee Joint Orthosis[END_REF]. Accordingly, the dynamics (1) can be rewritten as follows:

Ẋ = F 1 (X) + G 1 (X)τ e + G 2 (X)d (20)
where

G 1 (X) = G 2 (X) = 0 1 J
During the sitting rehabilitation position, the human torque τ h is considered an external disturbance and is estimated using the proposed NDO: d = τh and F 1 (X) is expressed as:

F 1 (X) = θ -1 J (-T g cos(θ) + Asign( θ) + B( θ)) (21)
while for the standing position, we have two different cases related to the gait cycle phases. The torque τ l differs between these two cases, as mentioned in Section II. During the swing phase, the lumped torque τ l is represented by (2); whereas, in the stance phase, it is considered an external disturbance, and can be estimated with the human torque using the NDO.

During the swing phase, the human torque τ h is considered as an external disturbance and is estimated using the NDO: d = τh and

F 1 (X) = θ -1 J (Tg sin(θ -θt)+Asign( θ)+B( θ) + τ l ) (22)
while during the stance phase, the human τ h and lumped τ l torques are both estimated using the NDO: d = τh + τl and

F 1 (X) = θ -1 J (Tg sin(θ -θt)+Asign( θ)+B( θ)) (23)
The proposed nonlinear observer is designed:

d = z + p(X) ż = L(-F 1 (X) -G 1 (X)τ e -G 2 (X)(z + p(X))) (24)
where

d = τh , p(X) = k 1 θ + k 2 θ, L = ∂p(X) ∂X = (k 1 k 2 ),
where k 1 and k 2 are two positive constants that allows us to deduce the unknown human/human-lumped torque values (considered here as an external disturbance) without using an additional sensor, thus reducing both the cost and dimensions of the studied active orthosis.

B. Reference trajectories generator

In this study, we adopt a time-based algorithm that facilitates the generations of a real-time knee joint reference trajectory by interpolating a number of samples matched with each phase sample speed (e.g., before stance/swing length) using a cubic spline function. The pre-allocated samples are extracted from the gait profile analyses of 20 healthy subjects [START_REF] Hutin | Walking velocity and lower limb coordination in hemiparesis[END_REF]. The resulting generated trajectory is illustrated in Fig. 3.

V. CLOSED-LOOP STABILITY ANALYSIS

To validate the asymptotic stability of the resulting closedloop system, we need to demonstrate, via the direct method presented in [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], that the cost function ( 18) is decreasing. Verifying that the objective function is a Lyapunov candidate ensures that this function decreases; hence, the controlled system is stable:

V (k + 1) -V (k) = ∆V (k) < 0 for ( X, ṽ) = (0, 0) ∆V (k) = 0 if ( X, ṽ) = (0, 0) (25)
In [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], the closed-loop asymptotic stability of an MPC system is validated by satisfying:

V (k + 1) -V (k) + l( X(k), ṽ(k)) ≤ 0 ( 26
)
and the following assumptions A1-A4: A1 :

X f ⊂ X, X f closed, 0 ∈ X f state constraint satisfied in X f . A2 : K f (x) ∈ U, ∀x ∈ X f control constraint satisfied in X f . A3 : f (x, K f (x)) ∈ X f , ∀x ∈ X f X f is a positively invariant under K f (.) . A4 : [F + ](x, K f (x)) ≤ 0, ∀x ∈ X f F (.) is a local
Lyapunov function . By designing the cost function as expressed in ( 17)-( 18), we can easily verify that:

l( X(k+i), ṽ(k+i)) = X(k+i) 2 Q + ṽ(k+i) 2 R ( 27
)
is a Lyapunov function candidate because

l( X(k + i), ṽ(k + i)) > 0 f or ( X, ṽ) = (0, 0) l(0, 0) = 0 (28)
In this case, inequalities ( 25) and ( 26) become equivalent because l( X(k), ṽ(k)) is a Lyapunov function. The cost function at the time step k + 1 is equal to:

V (k+1) = N-1 i=0 ( X(k+i+1) 2 Q + ṽ(k+i+1) 2 R ) + X(k+ N +1) T Q f (N +1) X(k+ N +1) (29)
A simple mathematical reformulation gives the expression:

V (k+1) = N i=1 ( X(k+i) 2 Q + ṽ(k+i) 2 R ) + X(k+ N +1) T Q f (N +1) X(k+ N +1) (30)
In the following, we express V (k + 1) as a function of V (k) to simplify the computation of ∆V (k):

V (k+1) = V (k)+ X(k+N ) T Q(N ) X(k+N )-l( X(k), ṽ(k)) +R(N )ṽ(k+ N ) 2 -X(k+ N ) T Q f (N ) X(k+ N ) + X(k+ N +1) T Q f (N +1) X(k+ N +1). (31)
Because X(k+ N +1) = A X(k+ N ) + Bṽ(k+ N ), then:

V (k+1)-V (k)+l( X(k), ṽ(k)) = X(k+ N ) T (Q(N ) -Q f (N ) +A T Q f (N +1)A) X(k+ N ) +2 X(k+ N ) T A T Q f (N +1)Bṽ(k+ N ) +(R(N )+B T Q f (N +1)B)ṽ(k+ N ) 2 (32)
By choosing l( X, ṽ) as a Lyapunov function l(0, 0) = 0 and l( X, ṽ) > 0 , ∀ ( X, ṽ) = (0, 0) , if we compute (32) and demonstrate that its right hand side is negative, we can then prove that the cost function is a Lyapunov candidate and consequently deduce the asymptotic stability of the closedloop system. Accordingly:

V (k+1) -V (k) + l( X(k), ṽ(k)) ≤ 0 ⇒ ∆V (k+1) ≤ -l( X(k), ṽ(k)) < 0 (33)
Equation ( 33) is verified for all ( X, ṽ) = (0, 0) because l( X, Ṽ ) is a Lyapunov function. This can be satisfied if the tuning parameters Q, Q f , R, N , and N c are carefully selected.

To ensure that, we impose a few conditions and summarize them in (34). If these conditions are satisfied, the negativity of the right hand side of (32) is guaranteed; thus the asymptotic stability of the controlled system is proved. This is expressed in the following equation as:

       Q(N ) -Q f (N ) + A T Q f (N + 1)A ≤ 0 R(N ) + B T Q f (N + 1)B ≤ 0 X(k+ N ) T A T Q f (N + 1)Bṽ(k+ N ) ≤ 0 V (k) ≥ 0, ∀k ≥ 0 (34)
We choose the weighting matrices Q and Q f as follows:

Q =                q 2 0 0 m 2 0 0 0 0 • • • 0 0 0 0 0 0 0 0 q 3 0 0 m 3 . . . . . . . . . . . . . . . . . . . . . q N 0 0 m N 0 0 0 0 0 0 0 0 • • • 0 0 0 0 q 0 0 m               
Therefore, the weighting matrix applied to the N th term of the cost function is equal to Q(N ) = q 0 0 m .

In addition, we can select the terminal-weighting matrix

as Q f (N + i) = -(N + i)(q + 1) 0 0 -N m(i + 1) N .
In this case, we obtain:

M (N ) = Q(N ) -Q f (N ) + A T Q f (N + 1)A (35)
which can be represented as

M (N ) = a b c d
, where a = -1 , b = c = -T s (N + 1)(q + 1) , d = (N + 1)(m -T 2 s (q + 1)) -2 N N m. To ensure that M (N ) is a negative semidefinite matrix, its principle minors must satisfy the conditions ∆ 1 = det(a) ≤ 0 and ∆ 2 = det(M ) ≥ 0. Their computation leads to ∆ 1 = -1, then the first condition is always satisfied and

∆ 2 = -m(N+1)+T 2 (N+1)(q+1)+2 N N m-T 2 (N+1) 2 (q+1) 2 .
To satisfy the latter condition, we choose m as follows:

m ≥ T 2 (N + 1)(q + 1)(N q + N + q) N (2 N -1) -1 (36)
To ensure that R(N ) + B T Q f (N + 1)B ≤ 0, the following condition should be satisfied:

R(N ) ≤ 2 N N T 2 m (37)
Let us now consider the third condition in (34). The state-space equation at time step k + N is equal to:

X(k+ N ) = A X(k+ N -1) + Bṽ(k+ N -1) (38)
Computing (38) with A and B as the respective state and input matrices, defined in [START_REF] Lv | On the Design and Control of Highly Backdrivable Lower-Limb Exoskeletons[END_REF], leads to:

x2 (k+ N ) = x2 (k+ N -1) + T s ṽ(k+ N -1) (39)
For a control horizon satisfying N c ≤ N -1, we can write:

x2 (k+ N ) = x2 (k+ N -1) + T s ṽ(k+ N ) (40) because ṽ(k + i) = ṽ(N c ) , ∀i ≥ N c , then ṽ(k+ N ) = ṽ(k+ N -1)).
Using the Z-transform, we obtain:

x2 (z) = T s z z -1 ṽ(z) (41)
which is equivalent to:

x2 (k+ N ) = T s Γ(k+ N )ṽ(k+ N ) (42)
with Γ(k +N ) being the unit step function. This proves that x(k+N ) and ṽ(k+N ) have the same sign; hence, the following inequality is always satisfied:

X(k+ N ) T A T Q f (N + 1)Bṽ(k+ N ) = -2 N T s N mx 2 (k+ N )ṽ(k+ N ) ≤ 0 (43)
With Q and Q f equal to those chosen in this study, and the tuning parameters q, m, N , and R(N ) satisfying the following conditions:

       m ≥ T 2 s (N +1)(q+1)(N q+N +q) N (2 N -1)-1 R(N ) ≤ 2 N N T 2 s m R(N ), q, m, N ≥ 0 0 ≤ N c ≤ N -1 (44)
this proves that (33) is satisfied for all k ≥ 0 and ( X, ṽ) = (0, 0). For ( X, ṽ) = (0, 0), the cost function is equal to:

V (k) = N -1 i=0 l(0, 0) + X(k+ N ) T Q f (N ) X(k+ N ) = 0 (45)
Accordingly, l(0, 0) = 0 and X(k+ N ) = 0 because ( X, ṽ) = (0, 0). This proves that the cost function V (k) is a Lyapunov candidate because the following conditions are satisfied for all k ≥ 0: ∆V (k) < 0 ∀( X, ṽ) = (0, 0)

V (k) = 0 for ( X, ṽ) = (0, 0) (46)
This concludes the asymptotic stability of the resulting closedloop system. The proposed EMPC controller inherits both the stability and performance properties of the MPC [START_REF] Bardaro | MPC-based control architecture of an autonomous wheelchair for indoor environments[END_REF].

VI. EXPERIMENTAL METHODS

In this section, we give a brief overview of the experimental operating conditions, the performed experimental scenarios, and the proposed metrics used to assess the performances of the proposed controllers. Our experiments were conducted with three healthy subjects (gender: 3 males; age: 30 ± 3 years old; height: 1.72 ± 0.04 m; weight: 67 ± 8 kg) while wearing the exoskeleton (cf. TABLE I). The wearers, who were explicitly informed about the experimental protocol and the expectations of this study, gave their consent before participating in the experiments. All subjects can ensure a complete flexion and extension of the knee joint without spasticity or contracture. Constraints on the knee joint position and velocity of the orthosis between the full knee extension (0 • ) and flexion (135 • ), are guaranteed by adopting mechanical stops. All precautions were taken to protect the health of the subjects who served as research subjects. Precautions were also taken to protect the privacy of the research subjects and their confidentiality. In section VII, we attempt to validate the proposed MPC-based control framework on the subjects wearing the active orthosis in three rehabilitation scenarios (i) in a sitting position where the hip angle is constant without ground interaction, and while walking on a treadmill with (ii) normal, and (iii) abnormal, walking gaits, respectively. In the last two scenarios, the interaction with the ground is considered as an external perturbation, and the hip angle is no longer constant. The proposed NDO is used to estimate the exerted human/human-lumped torque. To assess the performances of the proposed controllers, the following metrics have been analyzed: (i) the EMG activation of the muscles during walking and (ii) the exoskeleton torque required by the AAFO to provide functional assistance and satisfactory tracking of the desired trajectory. Indeed, exoskeleton torque, human torque and muscle activation, are good indicators of the subject involvement in the rehabilitation process. The EMG signals were sampled at 100 Hz. The raw EMG signals were filtered using a high-pass fourth order Butterworth filter with a cutoff frequency at 30 Hz, full wave rectified, and then filtered using a low-pass fourth order Butterworth filter with a cutoff frequency at 2 Hz [28] Figs. 4, 5, and 7 illustrate the evolution of the angular position and velocity, estimated human torque, and measured exoskeleton torque. Processed muscular activities of both Vastus Medialis (VM) and Biceps Femoris (BF) are also illustrated. The reference trajectories representing the walking profile of a healthy subject are depicted in red while the measured ones are in black for S1 and S2 and in green for S3. Each scenario was repeated three times with the concerned subject, and each time a different controller is applied (PID / MPC and EMPC). It should be mentioned that the wearers were only informed about the experimental protocol and the expectations of the study, without any indication on which controller was running in each session to not alter their performances and induce any bias in their movements. For each subject, the controller order was chosen randomly and only the operator was aware about which controller was running. The participants had no previous experiences with the exoskeleton, they were simply asked to perform some training movements with the robot, in both sitting and walking conditions prior to the experiment sessions, to get used to the robot and ensure they can perform movements comfortably. Each experimental session lasted for 60-s per controller/subject of a treadmill walking (60-s for the PID, 60-s for MPC, 60-s for the EMPC) while wearing the exoskeleton. After each session, participants were asked to sit on a chair during 3 minutes for rest. To compare the performance of each controller at the same scale and for the sake of clarity, we have extracted some of the obtained results with each controller (PID,MPC and EMPC) and have concatenated them together in Figs. [START_REF] Sherwani | RISE-based adaptive control for EICoSI exoskeleton to assist knee joint mobility[END_REF], 5 and 7. During the swing phase, the lumped torque is computed using (2), while it is considered as an external disturbance and is estimated with the human torque via the NDO during the stance phase.

VII. RESULTS AND DISCUSSIONS

In this section, we present and discuss the obtained results of the conducted experimental study to evaluate the performances of the proposed control framework with MPC and EMPC in three different scenarios. The performances of the proposed controllers are compared with a tuned PID controller, based on trial-and-error, until it reaches its best performances in terms of settling time, precision tracking, and stability. Today, PID controllers are used in most automatic process control applications in the industry to regulate flow, temperature, pressure, position, speed, level, and several other industrial process variables. In wearable robotics, beyond academia, PID-based controllers have also been adopted in several commercialized exoskeletons.

A. Scenario 1: Sitting position with two rehabilitation modes

In this first case, we ask subject S1 to wear the exoskeleton and to perform a movement of flexion/extension while sitting, without any interaction with the ground. This 20-s rehabilitation exercise is divided into two different parts. During the first 10-s, the wearer is asked to provide an assistive force by following the reference trajectory displayed on the laptop screen in real-time. During the last 10-s, he is asked to resist towards the desired trajectory by providing a resistive effort in the opposite direction. All the generated desired trajectories are presented in real-time on a laptop screen in front of the wearer while the change between the different phases is indicated by a beep. Only one gait cycle data per assistive/resistive phase and per controller (PID/MPC/EMPC) was cropped and all of them were concatenated together in Fig 4. This Figure clearly demonstrates that both MPC and EMPC exhibit better performances than PID controller in terms of position and velocity tracking. It can be noted that during the assistive scenario (first 10-s), the three controllers succeeded in tracking the desired trajectory. When comparing the results of the three controllers with respect to the resistive scenario, it can be demonstrated that both MPC and EMPC exhibited better tracking performances than PID while preserving a relatively smooth motion. The PID controller was unable to reach the peak values of the reference trajectory and exhibited a poor tracking performance.

B. Scenario 2: Healthy walking pattern

In this scenario, subject S1 is asked to walk on a treadmill for 30-s while wearing the exoskeleton. The reference trajectory is generated via a fixed-time based method, as described in Section IV-B. Fig. 5 gathers the obtained results with the three used controllers, where only 10-s of each controller performance are cropped and concatenated together. This figure indicates that MPC and EMPC ensure better position and velocity tracking performances than the PID controller. Because the NDO-estimated torque includes τ h and τ l together, it would be more reliable to use the EMG signals to interpret the results of the walking experiments. The EMG curves of both BF and VM muscles exhibit higher values when PID controller is used, which indicates the superiority of the MPC based controllers in providing better assistance than the PID controller, during the walking activity. MPC [START_REF] Kamalia | Trajectory generation and control of a knee exoskeleton based on dynamic movement primitives for sit-to-stand assistance[END_REF][START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF] and EMPC [START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF]30] controllers.

C. Scenario 3: Abnormal walking pattern

Because the walking gaits of people with reduced mobility differ from those of healthy people (as shown in Fig. 6), and to validate the robustness of the proposed approaches, two healthy subjects (S2 and S3) were asked to mimic the behavior of people suffering from lower-limb muscular weaknesses. These subjects did not participate in the previous set of experiments with a normal gait. S2 and S3 were trained to mimic abnormal gaits while wearing the active orthosis and walking on a treadmill, with a pre-selected speed of 2 km/h for 60-s, under three different assistive scenarios. Each experimental test had a duration of 60 s, where the useful data were extracted from the last 30 s, while the first 30 s were used to enable the wearer to adapt his walking pattern to the abnormal gait. Indeed, 10 s of each test (PID, MPC, and EMPC) were cropped and concatenated together in Fig. 7 for the two subjects to enable the comparison of the controllers. It can be deduced from Figs. 7 and8 that in terms of position and velocity knee joint tracking, both MPC and EMPC exhibit better results than the conventional PID controller. In fact, a relatively important tracking error can be observed in PID controller case, in particular, during the stance phase, owing to the ground reaction force acting as external perturbation. This error was reduced when using the MPC. In the case of the EMPC scheme, the NDO compensates for this disturbance in the same time step. Similarly, the error is reduced during the swing phase while using both MPC and EMPC. Furthermore, the MPC exhibits less control effort and better tracking performances than the PID. In addition, the EMPC required further control torque to ensure better tracking performances than both the PID and MPC controllers. Fig. 7: S2/S3 -Walking with an abnormal pattern while using PID [0,10], MPC [START_REF] Kamalia | Trajectory generation and control of a knee exoskeleton based on dynamic movement primitives for sit-to-stand assistance[END_REF][START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF] and EMPC [START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF]30] controllers.

It is worth noting that the EMG recordings exhibit less activation in the MPC case than the PID and EMPC controllers. In fact, the VM muscular group exhibits an important activation in the PID and EMPC cases to compensate the body weight during the heel strike of the stance phase. Furthermore, during late swing phase, the BF muscular group exhibits an important activation to ensure sufficient clearance of the swinging leg following the knee joint flexion. For both subjects, the EMG results for MPC exhibit relatively reduced peak activation, as illustrated in Fig 7 . In summary, according to the proposed metrics and the obtained results, in general, the proposed predictive controllers (MPC and EMPC) outperform the PID controller. Furthermore, although EMPC requires a larger control input than the conventional MPC, it exhibits the best performances among the controllers, in terms of trajectory tracking for all subjects. The obtained experimental results of scenario 3 are illustrated in the video available at: https : //youtu.be/tfvGgg4aFm0.

VIII. CONCLUSION AND FUTURE WORK

In this study, we propose a novel MPC framework dedicated to assistive and rehabilitation purposes, by adopting an actuated knee joint orthosis. The proposed framework includes (i) an exact input-to-state feedback linearization, (ii) a dedicated MPC or EMPC, (iii) a least-squares dynamic parameters identification, (iv) an NDO, (v) a Lyapunov-based stability analysis, and (vi) a fixed time based trajectory generator. Because considering the imposed constraints makes the MPC controller consume more computing time and memory space, these constraints are omitted in the first proposed controller (MPC). To address this limitation, an explicit MPC (EMPC) scheme is proposed, which can consider the system constraints. The proposed framework was validated via realtime experiments performed on three healthy subjects wearing the active orthosis while sitting, and while walking on a treadmill with normal and abnormal gait walking patterns. The performances of the proposed controllers are compared with those of a conventional PID controller. The obtained results clearly indicate that the proposed predictive controllers outperform the PID controller in terms of tracking performances, orthosis generated torque, and EMG activities of the involved muscles. Our future work will focus on realtime MPC control implementation with larger population of subjects including real patients, while considering the effects of (i) model uncertainties, (ii) external disturbances, and (iii) noises on the control performances and stability, and reducing the conservative behavior of MPC.

Fig. 1 :

 1 Fig. 1: Block diagram of the proposed control framework.

Fig. 2 :

 2 Fig. 2: Images of the subject wearing the EICOSI exoskeleton (a) in a sitting rehabilitation position, and (b) walking on a treadmill.

Fig. 3 :

 3 Fig.3: Illustration of the walking cycle of the lower limb while wearing the exoskeleton, and with the generated walking trajectory (for both reference and abnormal gaits).

Fig. 4 :

 4 Fig.4: S1 -Using PID, MPC, and EMPC controllers in a sitting position while exerting an assistive effort during the first 10-s and a resistive effort during the last 10-s.

Fig. 5 :

 5 Fig.5: S1 -Walking with a normal pattern while using PID [0,10], MPC[START_REF] Kamalia | Trajectory generation and control of a knee exoskeleton based on dynamic movement primitives for sit-to-stand assistance[END_REF][START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF] and EMPC[START_REF] Yue | An Efficient Model Predictive Control for Trajectory Tracking of Wheeled Inverted Pendulum Vehicles with Various Physical Constraints[END_REF]30] controllers.

Fig. 6 :

 6 Fig. 6: Average knee joint angle profiles in saggital plane, EMG profiles of the vastus medialis (VM) and biceps femoris (BF) along 100 % gait cycles for different conditions. (a): the subject walks with abnormal gait without wearing the knee joint orthosis, (b): the subject walks with abnormal gait while wearing the orthosis but without torque assistance, (c): the subject walks with abnormal gait with the torque assistance from the orthosis using proposed control method, and (d): the subject walks with normal gait while wearing the orthosis but without torque assistance. The right bar plot shows the VM and BF EMG activities in the different cases by using trapezoidal numerical integration with respect to 100 cycle points.

Fig. 8 :

 8 Fig. 8: S2/S3 -Root mean square (RMS) errors for angular position, torque, VMEMG and BFEMG, versus controllers.

TABLE I :

 I Identified parameters of EICOSI orthosis for the three participating subjects.

	Parameter	Symbol	S1	S2	S3
	Static friction coefficient (N.m)	A	2.0082	1.2067	0.7613
	Viscous friction coefficient N.m.s.rad -1	B	1.713	3.238	2.4539
	Inertia Kg.m 2	J	0.4325	0.2594	0.2525
	Gravity torque (N.m)	Tg	9.4199 10.4741 3.4379