Toward lower limbs movement restoration with input–output feedback linearization and model predictive control through functional electrical stimulation

Control Engineering Practice 20 (2012) 182–195

Contents lists available at SciVerse ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Toward lower limbs movement restoration with input–output feedback linearization and model predictive control through functional electrical stimulation S. Mohammed a,n, P. Poignet b, P. Fraisse b, D. Guiraud b a b

LISSI, University Paris Est Cre´teil—(UPEC), 120–122 rue Paul Armangot, 94400 Vitry Sur Seine, France LIRMM, DEMAR team, University of Montpellier II, CNRS/INRIA, 161 rue Ada, 34392 Montpellier, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 September 2008 Accepted 19 October 2011 Available online 10 November 2011

Functional electrical stimulation (FES) is used to excite paralyzed muscles that are no longer controlled by patients with spinal cord injuries (SCI). Appropriate stimulation patterns are chosen to stimulate intact muscles, in order to extend their overall performance, postponing thus the muscular fatigue during the daily activities such as standing, standing up, sitting down and walking. This paper presents the modeling and the control of a knee joint actuated by the quadriceps muscles. Appropriate stimulation patterns are computed as a function of the desired lower-limb knee joint movements. Parameters of the biomechanical model are identified based on experimental kinematic data. Model predictive control (MPC) is applied to the input–output feedback linearized (IOFL) system. IOFL allows linearization by inverting the system dynamics through a nonlinear feedback transformation. The described control approach is validated through different simulation scenarios for knee flexion– extension. Internal dynamics stability is mathematically proved and performances are compared to those produced by classical pole placements method. The controller has shown satisfactory results in terms of regulation, stability and robustness with respect to external disturbances. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Muscle modeling Rehabilitation engineering Functional electrical stimulation Input–output feedback linearization Model predictive control

1. Introduction Skeletal muscle is activated via a complex mechanism, which involves the brain, the spinal cord and the peripheral nervous system (Salmons, 1994). When muscles are no more controlled by the central nervous system, such as for spinal cord injured (SCI) patients, functional electrical stimulation (FES) can be used. Even though open-loop control strategies do not account for any changes in the muscle performances induced by fatigue or load changes, they are widely used in clinics due to their relatively simple implementation (Bajd et al., 1981; Guiraud, Stieglitz, Koch, Divoux, & Rabischong, 2006; Kobetic & Marsolais, 1994). On the contrary, closed-loop control strategies are rarely adopted in clinical use due to: (i) the need of sensors that have to be placed on the patient, (ii) safety issues, (iii) embedded computational cost. In both cases, in order to fulfill fair tracking performances, controllers should be based on realistic models which may be difficult to estimate accurately (Previdi, Ferrarin, Savaresi, & Bittanti, 2005). Thus, a trade-off between complexity and accuracy n

Corresponding author. Tel.: þ33 1 41807318; fax: þ33 1 41807376. E-mail addresses: [email protected] (S. Mohammed), [email protected] (P. Poignet), [email protected] (P. Fraisse), [email protected] (D. Guiraud). 0967-0661/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.conengprac.2011.10.010

must be found. Regarding closed loop control, classical linear controllers fail to provide satisfactory tracking performances because of the inherent nonlinearities in the musculoskeletal system (Riener & Fuhr, 1998). Some authors have used a simple PID controller (Wood et al., 1998), a knee extension controller (Poboroniuc, Wood, Donaldson, Fuhr, & Riener, 2003), a combination of feedback and feedforward controllers or an adaptive approach (Ferrarin, Palazzo, & Riener, 2001), a combination of an adaptive nonlinear compensator with a sliding mode control (Kobravi & Erfanian, 2009; Shen, 2009). Others have used a first or second order switching curve in the state space to control patient movements: the on/off controller (Mulder, Veltink, & Boom, 1992) and the ONZOFF controller (Poboroniuc, Fuhr, Wood, Riener, & Donaldson, 2002) in the so-called ‘‘controller-centered’’ strategies. The so-called ‘‘subject centered’’ strategies, PDMR: patientdriven motion reinforcement (Riener & Fuhr, 1998) and CHRELMS: Control by Handle REactions of Leg Muscle Stimulation (Donaldson & Yu, 1996), introduce voluntary contributions of the patient’s upper body as a part of the control scheme. Many other movements through FES have been achieved (Alon, Levitt, & McCarthy, 2007; Braz, Russold, & Davis, 2009; Chen, He, & Hsueh, 2009; Freeman, Davies, Lewin, & Rogers, 2008; Freeman et al., 2009; Vette, Masani, Kim, & Popovic, 2009). Grasping and releasing controls of the forearm, controls of the wrist, elbow and shoulder

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

Nomenclature bl F F_ Fq g Np Nu I J Ks L0 lc

shape factor normalized force–length relation derivative of the normalized force–length relation extension force generated by the quadriceps gravitational acceleration prediction horizon control horizon inertia about the knee joint cost function stiffness of the muscle serial element thigh length distance from the knee joint center to the shank center of mass

have been also investigated (Bhadra & Peckham, 1997). Within the literature on the muscle modeling for FES, some studies have treated the dynamic contraction of the human muscle (Bestel & Sorine, 2000; Zahalak, 1981). Others have represented the muscle as a nonlinear function of recruitment with dynamic activation, angle and angular velocity dependencies (Riener & Fuhr, 1998; Veltink, Chizeck, Crago, & El-Bialy, 1992). The aforementioned works have used either static representations of the muscle model in a closed loop control scheme or have represented, independently, detailed models of the muscle contractions/relaxations and activation dynamics at the muscular fiber level. A Hammerstein structure has been used in Le, Markovsky, Freeman, and Rogers (2010) to represent the electrically stimulated muscle models of subjects who have incomplete paralysis. Iterative learning techniques have been used to control the system and to identify the model parameters (Freeman et al., 2009; Le et al., 2010). In this paper, an input–output feedback linearization (IOFL) cascaded with MPC controller is proposed and evaluated, in simulations, to highlight its performance in terms of its capability of tracking a pre-defined reference trajectory and its robustness with respect to external disturbances. Controller performances are compared with those of a classical pole placement controller. In Section 2, the model of the knee joint actuated by the quadriceps muscles is presented. Its state space formulation and its parameters identification based on experimental setup are also shown. Section 3 details the IOFL of the biomechanical model, its internal dynamics stability and the corresponding MPC controller. Comparison study of the controllers’ performances is addressed in Section 4. Discussion and conclusions are stated in Sections 5 and 6.

2. Biomechanical model 2.1. System modeling The nonlinear muscle model (El-Makssoud, Guiraud, & Poignet, 2004) used in this study forms an extended physio-mathematical formulation of the macroscopic Hill and microscopic Huxley concepts. It reflects the dynamic phenomena occurring during muscle’s contraction and relaxation. Functional electrical stimulation signal is characterized by three parameters, as they are delivered by the stimulator: the pulse width PW, the signal frequency freq and the current amplitude Ia. Frequency of the stimulation signal is kept constant while PW or Ia are computed based on the use of a control strategy. The number of recruited motor units increases as a function of both PW and I of the stimulus and is represented by an activation model that expresses the ratio

Lc Lc0 Liq Lq Ls Ls0 m Q, R rp uch

l

m

183

muscle contractile element length resting muscle contractile element length distance between the knee joint rotation center and the insertion point of the quadriceps on the shank length of the quadriceps muscle serial element length resting muscle serial element length mass of the shank weighting matrices for the cost function pulley radius chemical control input to the muscle model knee joint elasticity factor knee joint damping factor

of recruited fibers. The contraction dynamics are expressed by a set of nonlinear differential equations representing the mechanical model. Parameters of the knee joint model, actuated by the quadriceps muscle, are identified on paraplegic patients and healthy subjects using various algorithms (DeLeva, 1996; Gautier & Poignet, 2002; Gill & Murray, 1978). Kinematic data such as lower limb positions, velocities and accelerations are computed using a video-based motion capture system (Vicon 360—Oxford Metrics, Oxford, United Kingdom). 2.2. Knee joint model Closed-loop control of the muscles actuating the knee joint of a paraplegic patient constitutes a prerequisite step toward any upward mobility such as standing up, standing, walking, climbing stairs, etc. The biomechanical model chosen in this study consists of two segments representing the shank and the thigh connected to each other by one degree of freedom revolute joint (knee joint). The thigh is assumed to be fixed with respect to the patient while the shank is free to move around the knee joint (Fig. 1). The quadriceps muscle acts as an agonist extensor muscle (Fq causes the knee extension) while the knee flexion is induced by gravity (y ¼ 01 corresponds to full extension of the knee and y ¼ 901 corresponds to the resting position). Fq, which is related to the stimulation signal parameters through the muscle modeling formulation (Section 2.3), represents the input of the biomechanical model while y is the corresponding output. This model has been chosen for many reasons: it is easy to implement by reducing the number of degrees of freedom of the multilink human musculo-skeletal system, it allows the reproduction of knee movements (flexion/extension) and it allows the control, in future studies, of the antagonist muscles (quadriceps/hamstrings) during co-contraction. In Fig. 1, the geometric equations allow the evaluation of the quadriceps length Lq as a function of the knee angle y: qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi Lq ðyÞ ¼ L20 r 2p þr p y þ L2iq r 2p ð1Þ From the above equation, one can deduce the relative elongation of the quadriceps: qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi L20 r 2p þ r p y þ L2iq r 2p Lq0 Lq Lq0 Eq ¼ ¼ ð2Þ Lq0 Lq0 Lq0 corresponds to the initial quadriceps length. The moment arm of the quadriceps is assumed to be constant and equal to the pulley radius rp (Fig. 1). The equation of motion is a nonlinear second order equation. It is a function of the knee joint inertia I,

184

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

Fig. 1. Biomechanical model of the knee joint actuated by the quadriceps agonist muscles.

Fig. 2. Mechanism of force generation from FES signal: DEMAR muscle model.

gravity, joint damping factor m and elasticity factor l: Iy€ ¼ mglc cosðyÞmy_ ly þ r p F q

ð3Þ

lc corresponds to the distance from the knee joint center to the shank’s center of mass. Identification of the above parameters is realized based on an experimental setup and is presented in Section 2.6. 2.3. Muscle model 2.3.1. Force generation by FES Generally, open-loop stimulation patterns are empirically tuned upon doctor’s estimation. This has proved that fast muscular fatigue and thus loss of natural desired movements are rapidly reached. In order to synthesize optimal stimulation patterns and to carry out model-based closed loop control, a better understanding of the muscle modeling is needed. The muscle model used in this study is based on internal physiological characteristics gathering two levels: a microscopic level involving the sliding actin–myosin filaments and a macroscopic level represented by contractile and elastic elements (El-Makssoud, 2005; El-Makssoud et al., 2004). This model describes the internal physiological mechanisms resulting from the application of FES. As shown in Fig. 2, the stimulation signal induces a force applied at the skeletal model through the activation and the mechanical model of the muscle. The activation model expresses the ratio of the recruited fibers a through the muscle fiber recruitment model and the chemical control input variable uch through the chemical control generator, while the mechanical model (detailed below) represents the force generation as a function of a, uch and the force length relationship (see Section 2.6.2). uch represents the

Fig. 3. Muscle model with the contractile element Ec, serial element Es and parallel element Ep.

chemical control input that is linked to the muscle contraction and relaxation through the calcium level release (Bestel & Sorine, 2000). Activation model translates the stimulation signal parameters into two-control variables a and uch. Formulations of these two functions are detailed in El-Makssoud et al. (2004). a and uch give, through the mechanical model, the corresponding force generated by the muscle. The resulted muscle force leads to the knee joint position by means of the knee joint dynamic model. The model shown in Fig. 2 is adopted in this study to develop strategies for simulation, motion synthesis and motor control during the clinical restoration of movement.

2.3.2. Mathematical formulation Fig. 3 illustrates the muscle model having a parallel element Ep and two elements in series: the serial element Es and the contractile element Ec. The contractile element Ec of length Lc is freely extensible and is capable of shortening when activated. The elastic element Es, of length Ls, models the isometric deformation of the muscular fiber. Since Ec and Es are in series, Fc is equal to Fs.

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

quadriceps muscle model. The state model of the knee-quadriceps can be expressed as 8   s0 aF m x1 su x2 x1 sv ax1 r p x4 > > _ > uch su x1 uch  ¼ s a K þ s q x m v 1 0 > > 1þ px1 sv qx2 L0 ð1 þ px1 sv qx2 Þ > >   > > > < x_ ¼ s0 aF m su x2 u þ bx1 r p x4 sv ax2 r p x4 2 1þ px1 sv qx2 ch L0 ð1þ px1 sv qx2 Þ > > > _ x ¼ x4 > > > 3 > > 1 > > x : _ 4 ¼ ½x2 r p lx3 mx4 mglc cos x3  I

The parallel element Ep of length L (L ¼ Lc þLs ) is added to account for elasticity of the muscle. This parallel element represents the passive resistivity of the muscle and has very limited influence on the overall generated force. In the following, the effect of this element has been neglected since quadriceps hyperextension is not considered during flexion–extension of the knee joint. The mechanical model has been described by two sets of differential equations (4) where the outputs are Kc and Fc representing, respectively, the stiffness and the force generated by the contractile element. Km and Fm correspond to the maximum values of Kc and Fc, respectively. These equations can be expressed as follows (El-Makssoud, 2005):   8 s0 aF m K c su F c K c sv aK c > _ > uch  E_ > < K c ¼ s0 aK m su K c þ sv q 1 þ pK sv qF 1 þpK c sv qF c c c     s0 aF m su F c bK c sv aF c > > _ > : F c ¼ 1þ pK s qF uch þ 1þ pK s qF E_ v v c c c c

ð6Þ 2.5. A simplified nonlinear biomechanical model The system (4) coupled with the knee joint dynamics (3) is a nonlinear multivariable system (6). The goal is to obtain a feasible solution before applying the controller to the original nonlinear plant. Some assumptions are justified as follows:

ð4Þ s0 ¼

1þ su , 2

L L Ec ¼ c c0 , Lc0

a¼

L0 , Lc0

b ¼ L0 ,

L L Es ¼ s s0 , Ls0

p¼

1 , Ks

LL0 E¼ , L0

q¼

1 Lc0 K s

 The chemical control uch is considered as a positive constant

L ¼ Lc þ Ls

Lc0, Ls0 and L0 represent respectively the length at rest position of the contractile element, the serial element and the parallel element. Ks represents the serial element stiffness. Ec , Es and E represent respectively the relative elongation of the contractile, serial and parallel elements, su ¼ signðuch Þ and sv ¼ signðE_ c Þ are sign functions related respectively to the control and the velocity of the contractile element. For a complete mathematical description of the muscle model, the reader could refer to the literature (El-Makssoud et al., 2004).



ð5Þ

ð7Þ

Taking into account the above assumptions, the muscle-knee model (6) will have a simplified nonlinear form, where u ¼ a is the control input (uch is considered constant). The new plant model is detailed in Mohammed (2006) and can be expressed by

x¼½x1 . . . x4 T ¼ ½K c F c y y_ T is the state vector and u ¼½uch aT the control vector. The variable y represents the knee joint angle and the variables K c , F c , uch , a represent the state variables of the

400

90 Simplified model Original model

350

80

300 Knee position θ (°)

Extension force Fc (N)

L0 1 _ E_  Fc Lc0 K s Lc0

Derivation of (7) is presented in Appendix A. In the above conditions and regarding the movement range of the knee joint, the term ð1=Lc0 K s ÞF_ c is less than 103 ðL0 =Lc0 ÞE_ and consequently can be neglected.

The model of the quadriceps muscle actuating the knee joint is formulated as a nonlinear state model: x_ ¼ f ðx,t,uÞ

(su ¼ signðuch Þ ¼ 1Þ and 9uch 9 ¼ uch ) corresponding to a muscular fiber fusion, this is validated by the fact that the clinical stimulation frequency is much greater than the muscular fiber fusion threshold (El-Makssoud, 2005). Consequently, during stimulation, no relaxation but only contraction of the muscular fibers occurs which leads to fibers shortening (sv ¼ 1 and E_c 9 ¼ E_c Þ. Since only dynamic movements of the knee joint are considered and isometric stimulations are excluded, the stiffness Ks of the serial element, representing the tendon, is much greater than the stiffness Kc of the contractile element:

E_c ¼

2.4. Muscle-knee: state space model

250 200 150 100

70 60 50 40 30

50 0

185

Simplified model Original model

0

0.5

1 Time (s)

1.5

2

20

0

0.5

1 Time (s)

Fig. 4. Nonlinear knee–muscle models: original and simplified.

1.5

2

186

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

the following set of differential equations: 8 x_ ¼ a1 x1 þ a2 x1 x4 þ a3 Fðx3 Þu > > > 1 > < x_ 2 ¼ a4 x1 x4 þ a1 x2 þ a2 x2 x4 þ a5 Fðx3 Þu x_ 3 ¼ x4 > > > > : x_ 4 ¼ a6 x2 þ a7 cosðx3 Þ þ a9 x3 þ a8 x4

ð8Þ

with a1 ¼ uch , rp a6 ¼ , I

a2 ¼

rp , Lc0

mglc a7 ¼ , I

a3 ¼ uch K m , m a8 ¼ , I

a4 ¼ r p ,

l a9 ¼ , I

a5 ¼ uch F m 2

!2 3 l1 5 Fðx3 Þ ¼ exp4 bl

The above hypotheses are validated through simulations as shown in Fig. 4. It presents the responses of both the original nonlinear model (6) and the simplified one (8). Fig. 4(a) shows the force generated by both models. The original nonlinear model presents some oscillations reflecting the contraction–relaxation cycle. Fig. 4(b) presents the knee position output of both models. It is observed that the simplified model fits well with the original model in terms of force generation as well as angular knee position. 2.6. Parameters identification The parameters of the biomechanical knee–muscle system (8) are identified using different protocols. The geometric parameters such as the insertion point of the quadriceps muscle on the shank Liq, the muscle length L0 and the moment arm rp are identified based on the Hawkins model (Hawkins & Hull, 1990) and using the Levenberg–Marquardt algorithm (Gill & Murray, 1978). The knee joint dynamic parameters, such as joint elasticity l and damping factors m, are identified through a linear least square algorithm (Gautier & Poignet, 2002). Muscle parameters such as the maximal muscle force Fm and the force–length relationship Fðx3 Þ are identified using nonlinear interpolation. Inertial parameters (shank mass and inertia) are estimated based on DeLeva

(1996) regression equations. Muscle parameters, not yet identified on humans, such as the muscle stiffness and the contractile– elastic muscle length distribution, are taken from the literature (El-Makssoud et al., 2004).

2.6.1. Identification of the knee joint parameters Kinematic data of the knee joint are measured through the Vicon& motion analysis system, using the passive pendulum test. Details of the experimental procedure are shown below. This test is carried out by recording the knee joint angle y and velocity y_ during passive movements. Table 1 summarizes the identified knee joint parameters for a given subject, which correspond to the thigh length L0, the quadriceps moment arm rp, the quadriceps muscle insertion point Liq and their standard deviations. These parameters are satisfactorily identified and closely correspond to those found in the literature (Kromer, 1993). The standard deviations are less than 4%. As shown in Fig. 5, the muscle length computed by both the model described above (Fig. 1) and the Hawkins model match closely. The dynamic parameters of the knee joint such as knee elasticity and damping factors are identified based on the following equation, using the linear least square method. During passive movement, the active moment produced by the muscle is equal to zero. This moment is then omitted in (3) as follows: Iy€ ¼ mglc cos ylymy_

Three subjects have participated in this study: two paraplegics and one healthy subject as shown in Table 2. Three exclusion criteria have been applied to choose appropriate participants:

 complete flexion and extension of the knee joint,  no spasticity or contracture during flexion and extension of the knee joint,

 any pathology of the knee joint. Table 2 Data related to people who participated in the experiments.

Table 1 Knee joint parameter identification. Parameter Value (m) Standard deviation (%)

L0 0.3726 0.335

rp 0.0397 3.230

Liq 0.0401 3.138

Subject

Age (year)

Height (m)

Weight (kg)

SCI level

Accident date (months)

Healthy Paraplegic 1 Paraplegic 2

28 29 29

1.82 1.80 1.88

86 80 100

– T5 T12

– 5 12

120

Reference Estimation

120

110

110

100 Knee angle (°)

Knee angle (°)

130

100 90

80 70

70

60

0

1

2

3 4 Time (s)

5

6

Reference Estimation

90

80

60

ð9Þ

50

0

1

2

3 4 Time (s)

5

6

Fig. 5. Parameters identification of the biomechanical model (the reference model is the Hawkins model).

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

187

Fig. 6. (a) Vicon analysis system. (b) Healthy subject. (c) Paraplegic patient.

Table 3 Estimation of the anthropometric parameters (inertial moment, mass and length).

Table 4 Dynamic parameter identification.

Subject

I (kg m2)

m (kg)

l (cm)

Subject

m (N m s/rad)

sd (%)

l (N m/rad)

sd (%)

Healthy Paraplegic 1 Paraplegic 2

0.1682 0.1536 0.2092

4.8852 4.5830 5.6162

19.46 19.2 20.26

Healthy Paraplegic 1 Paraplegic 2

0.0838 0.0659 0.0897

6.4845 13.3589 10.6973

0.17 0.1095 0.0529

2.266 6.0050 10.1739

The identification process, applied to both healthy and paralyzed subjects, is performed with the subject lying semi-supine and the lower leg hanging over the edge of a chair (Fig. 6). Vicon 360 is an optical motion-capture system. The kinematic data are collected through two 50-Hz cameras while a third camera is placed in front of the subject and aligned along the axis of the lower limb movements. Three reflective passive markers are placed on the subject’s joints (hip, knee and ankle) in order to calculate the knee joint movement assuming that the movement remains in the sagittal plane. This motion analysis system has the advantage of not overwhelming the subject by an external sensor that could affect the accuracy of the identification such as body attached goniometers and accelerometers. Two to three tests are performed on each subject in the presence of an assistant who remains beside the patient throughout all the experiments. The pendulum test is carried out: the operator raises the shank of the subject to a given angle (approximately 451) then allows the shank to swing freely until it reaches the resting position (901), while measuring the EMG signals of the muscles actuating the knee joint (quadriceps and hamstrings). The EMG signal analysis serves only to identify and reject any undesired voluntary/reflex muscle contractions when the EMG signal level exceeds a given threshold value. The anthropometric parameters, such as the shank mass and inertia, are estimated by using the regression equations proposed by DeLeva (1996), and by measuring the weight and height of each subject. The anthropometric parameters of the subjects who participated in the identification process are shown in Table 3. The dynamic parameters of the knee joint such as elasticity and damping factors are identified by means of linear least square methods (Gautier & Poignet, 2002). The knee angle position are extracted from the kinematic data, while the velocity and acceleration are computed by numeric derivation using a Butterworth low-pass filter. This filter has been adopted because it provides maximum passband flatness and its cut-off frequency is chosen empirically for each participant upon the knee movement amplitude and the subject shank inertia. The identified parameters as well as their standard deviations (sd) are shown in Table 4. These results are close to those found in the literature (Ferrarin et al., 2001), computed in the same context.

2.6.2. Identification of the muscle parameters The force–length relationship (Riener & Fuhr, 1998) expressed by (10) and the maximal isometric force generated by the quadriceps muscle are identified using a special experimental platform (Fig. 8) (Mohammed, 2006). The platform is equipped, as shown in Fig. 9, by a force sensor, position sensor, mechanical shank and a foot blocking system, allowing force and position measurements during isometric contraction tests. The isometric force generated by the quadriceps muscle has been measured at different knee joint angles ranging from 01 to 1001 as shown in Fig. 7. Maximal stimulation current is set to 150 mA and maximal pulse width is set to 500 ms: 2 !2 3 l1 5 ð10Þ FðyÞ ¼ exp4 bl With l ¼ L=Lq0 , L represents the current length of the quadriceps muscle and Lq0 represents its length at the resting position. From (10), it can be noted that bl can be easily identified by expressing bl as a function of the measured force and the measured joint angle y. The remaining parameters, not yet experimentally identified, particularly the muscle stiffness and the contractile– elastic muscle length distribution, are taken from El-Makssoud et al. (2004) and are reported in Table 5.

3. Controlling the knee joint using input–output feedback linearization and model predictive control Unlike a previous study (Mohammed, 2006) where the control of the knee joint was achieved using an approximate local linearization, this study computes the exact IOFL of the knee muscle model. A locally linearized system is sensitive to external disturbances and allows only a satisfactory control around the local linearization neighborhood. In the case of input–output feedback linearization (IOFL), the system dynamics are decoupled into external and internal dynamics. While the external dynamics are easily controlled, stability of the internal dynamics must be ensured in order to guarantee the stability of the overall system (Isidori, 1990). This paper adopts a MPC due to the relative slow internal dynamics of the biomechanical model, the constraints on

188

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

Fig. 7. Isometric quadriceps force measurement at different knee joint angles.

Table 5 Parameters of the quadriceps muscle. Muscle parameter

Quadriceps

Unit

Ks Lc0 Ls0

1  104 41  10  2 8  10  2

N/m m m

the input stimulation current and the output knee joint position. The MPC formulation, also known as receding horizon control (Allgower, Badgwell, Qin, Rawling, & Wright, 1999, Chap. 12), is usually stated as an optimization problem subject to physical coherent constraints and is solved with classical optimization algorithms. The MPC is widely used in different applications due to its interesting properties (Camacho & Bordons, 1995). These controllers have been very successfully used in the industry, especially in refining and chemical processes (Allgower et al., 1999; Nevistic & Morari, 1996). They have also been used in robot-based applications (Azevedo, Poignet, & Espiau, 2002; Richalet et al., 1996). In this paper, the MPC has been applied to a continuous time physiological musculo-skeletal model as opposed to the black box models used in the past (Schauer & Hunt, 2000). In this case, the nonlinearities are compensated by the IOFL inner-loop while the MPC outer-loop is used to guarantee stability and to reject undesired disturbances. 3.1. Feedback linearization Since a simplified nonlinear model of the biomechanical knee– muscle model is used, additional approximate linearization (Appendix B) would lead to an inaccurate model that does not reflect the physiological modeling cited above. That is why, an exact IOFL of the system (8) is used instead of using its approximation. IOFL is a systematic method that allows linearization by inverting the knee–muscle nonlinear plant dynamics (Fig. 2) through a nonlinear feedback transformation (Fig. 10). Subsequently, the simplified nonlinear model can be formulated as ( x_ ¼ f ðxÞ þ gðxÞu ð11Þ y ¼ hðxÞ with 2

f 1 ðxÞ

3

2

a1 x1 þa2 x1 x4

3

6 7 6 7 a4 x1 x4 þ a1 x2 þ a2 x2 x4 6 f 2 ðxÞ 7 6 7 7 6 7 f ðxÞ ¼ 6 6 f ðxÞ 7 ¼ 6 7 x 4 3 4 5 4 5 f 4 ðxÞ a6 x2 þ a7 cos x3 þa9 x3 þ a8 x4

ð12Þ

gðxÞ ¼ ½g 1 ðxÞ,g 2 ðxÞ,g 3 ðxÞ,g 4 ðxÞT ¼ Fðx3 Þ½a3 ,a5 ,0; 0ÞT and hðxÞ ¼ x3 . The relative degree of the nonlinear system (8) is computed using the Lie derivative (L) mathematic tool and is equal to 3 (Appendix C.1). The system (11) can be formulated as (Isidori, 1995; Slotine & Weiping, 1991) ( ðrÞ y ¼ v, yð0Þ ¼ y0 ð13Þ Z_ ¼ cðy y_ . . . yðr1Þ Þ, Zð0Þ ¼ Z0

Z, of dimension ðnrÞ, is chosen such that the nonlinear transforT mation z ¼ ½z ZT T is a diffeomorphism where z ¼ ½y y_ . . . yðr1Þ T and c is a nonlinear function. The new form (13) is called the normal form of Byrnes–Isidori. The second equation of (13) represents the internal dynamics of the nonlinear system (11). As a result, the IOFL allows the dissection of the nonlinear dynamics behavior (Fig. 11); therefore, the system description (13) can be expressed in the new framework as follows (Guemghar, 2005): 8 _ > i A f1; 2, . . . ,r1g > < zi ¼ zi þ 1 , ð14Þ z_r ¼ v > > : Z_ ¼ cðZ, zÞ where z1 ð0Þ ¼ y0 and Zð0Þ ¼ Z0 . System (11) is globally linearizable and a linear controller will be used to control the system’s output. The equation Z_ ¼ cðZ, zÞ represents the internal dynamics and corresponds to the last ðnrÞ equations. If uncontrolled, the internal dynamics can lead to system instability. When dealing with IOFL, the nonlinear system dynamics are decoupled into an external part represented by a linear relation between input and output, and an internal unobservable part (Slotine & Weiping, 1991). When designing a control strategy, the whole system’s dynamics should be taken into account. While the external part is easily controlled by a linear controller so that the output behaves as desired, the internal dynamics should stay bounded and their stability should be addressed carefully. Consequently, the internal dynamics stability cannot be neglected and are characterized by studying their zero dynamics. By definition, zero dynamics are computed using the normal form, by maintaining the system’s output y to zero or to a constant value, in which case all its derivatives with respect to time will vanish. The zero dynamics describe a movement restricted to a continuous surface of dimension ðnrÞ, defined by z ¼ 0. In order to operate in zero dynamics zone, the input control u is computed such that the output y remains equal to zero, that is to say yðrÞ ¼ v ¼ 0. Where v, as shown in Fig. 10, represents the controller output applied to the linearized system. The system dynamics can then be written in the following form: ( z_ ¼ 0 ð15Þ Z_ ¼ cðZ,0Þ By definition, system (15) represents the zero dynamics of the initial system (11) (Slotine & Weiping, 1991). Derivation of the internal dynamics and their stability are developed in Appendix C.2. 3.2. Pole placement controller Suppose that k1, k2 and km are the controller parameters. By using the first equation of (13), which represents the linearized system or the external dynamics, a classical pole placement control is applied. The controller output v, that controls the FES signal, is applied to the linearized system and can be expressed as follows: _ 2 e€ v ¼ yð3Þ km ek1 ek d

ð16Þ

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

189

Fig. 8. Experimental platform and muscle parameter identification.

Fig. 9. Experimental platform.

Fig. 10. Input–output linearization and model predictive controller.

where e ¼ yyd represents the tracking error between the reference trajectory yd and the current trajectory y. The tracking error converges exponentially toward zero with respect to the following differential equation: eð3Þ þ k2 e€ þ k1 e_ þkm e ¼ 0

ð17Þ

The input control of the nonlinear system (Fig. 10) u is derived by Isidori (1995) and takes the following form: u¼

Fig. 11. Input–output feedback linearization (Guemghar, 2005).

1 Lg L2f y

...

€ y€ d Þk1 ðy _ y_ d Þkm ðyyd Þ ½L3f y þ yd k2 ðy

ð18Þ

As a result, the input–output behavior between the external input v and the output y is fully linearized. In Fig. 12, the controller robustness with respect to variation in initial knee angle position

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

100

90

90

80

80

70

Knee position (°)

Knee position (°)

190

70 60 50 40 30 20

60 50 40 30 20 10

0

0.5

1 Time (s)

1.5

0

2

0

0.5

1 Time (s)

1.5

2

Fig. 12. Controller robustness with respect to variation in initial position.

Fig. 13. Details of the MPC application.

has been shown. The controller succeeds to track the desired position with a shift of 101 from the rest position (901) for both step and continuous desired positions. 3.3. Model predictive controller The proposed control strategy shown in Fig. 10 is based on the use of IOFL cascaded with model predictive control (MPC). The biomechanical knee–muscle model is input–output feedback linearized through decoupling of the model internal dynamics from the input–output model behavior. Consequently, the cascade control scheme consists of two loops: an inner loop representing the IOFL and an outer loop representing the model predictive control of the internal dynamics. The IOFL could be viewed as a pre-compensator before applying MPC while predictive control could be seen as a systematic way of controlling internal dynamics (Guemghar, 2005). The MPC is usually formulated as a constrained optimization problem: minJðkÞ ¼ N

uk p

NX p 1

Jyk þ i9k ydk þ i9k J2Q þ

min

k

JDuk þ i9k J2R

ð19Þ

i¼0

i¼0

subject to 8 > < umin r uk r umax ymin ryk rymax > : Du r Du r Du

NX u 1

ð20Þ max

where y, yd and Du, represent respectively the current and the reference system’s output and the variation in control vector. In Eq. (19), Jyk þ i9k ydk þ i9k J2Q is equal to ðyyd ÞTk þ i9k Q ðyyd Þk þ i9k and JDuk þ i9k J2R is equal to DuTk þ i9k RDuk þ i9k . R and Q are diagonal

weighting matrices of the cost function chosen such as Q Z0 and R4 0 in order to guarantee the system’s stability (Akesson, 2003). In the following simulations, R and Q are tuned empirically and set respectively to 0.5 and 1. For each sample, a finite horizon optimal control problem is solved over a fixed interval of time called the prediction horizon Np. Predictive control (Fig. 13) consists of computing the control vector uk at instant k for consecutive inputs over the control horizon Nu by optimizing the objective function J(k) (19) subject to constraints (20). After Nu samples, the control signal is assumed to be fixed over a horizon of dimension ðNp Nu Þ. When the solution of the optimal control problem is obtained, the value of the first control variable in the optimal trajectory is applied to the process. The rest of the predicted control variable trajectory is discarded and the entire procedure is repeated at the next sampling interval (Akesson, 2003).

4. Results The simulation results presented in this study are performed in a Matlab-Simulink environment using the ‘‘ode45’’ integration algorithm with variable step size. The MPC optimization problem is solved using the mpcsim function of the Matlab Model Predictive Control Toolbox,1 which simulates closed-loop systems with saturation constraints on the manipulated variables. According to the values used in practical experiments, sampling period and stimulation frequency are set respectively to 0.01 s and 25 Hz. Simulation results are shown in Fig. 14 where initial 1

http://www.mathworks.com/products/mpc

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

191

100

Step PP MPC

Knee position (°)

90 80 70 60 50 40

0

1

2

3

4

5 Time (s)

6

7

8

9

Normalized input control

1

10

PP MPC

0.8 0.6 0.4 0.2 0

0

1

2

3

4

5 Time (s)

6

7

8

9

10

Fig. 14. MPC and PP controllers: desired position corresponds to a 451 knee extension (MPC is the model predictive control and PP is the pole placement controller).

Knee position (°)

100

MPC PP Disturbance Step

80 60 40 20 0

0

1

2

3

4

5 Time (s)

6

7

8

9

10

Normalized input control

1 MPC PP

0.8 0.6 0.4 0.2 0

0

1

2

3

4

5

6

7

8

9

10

Time (s) Fig. 15. Controllers’ behaviors with respect to external disturbance of 101.

conditions correspond to the rest position (y ¼ 901). The desired position is lowered to y ¼ 451 at t ¼2 s, which corresponds to half extension of the knee joint induced by activating the quadriceps muscle. A comparison has been drawn in terms of input control and output regulation between the performances of the MPC and those of a classical linear controller based on pole placement (PP).

The PP controller serves as a reference. Both controllers converge to the desired position in a finite time with respect to the lower/ upper bounds of the control input. These controllers have been simulated under the same initial conditions. Unlike the PP controller which takes a longer time to converge to the desired position and presents, at the same time, an input saturation

192

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

Reference

Mass

Gravity

Inertia

Viscosity

Knee joint position (°)

95 90 80 70 60 50 40

0

1

2

3

4

5

6

Time (s) Fig. 16. Robustness with respect to parameter uncertainties.

100

91.5

Knee joint angle (°)

Knee joint angle (°)

90 80 70 60

91

90.5

50 90 40

0

2

4 6 Time (s)

8

10

0

0.2

0.4 0.6 Time (s)

0.8

1

Fig. 17. Guaranteed constraints while tracking a continuous desired trajectory.

during the transient period, the MPC based controller converges to the desired position in a relatively minimal time. In the MPC case, the input control shows a tolerable overshoot. The optimization time required to reach the desired position is around 20 ms which is encouraging for a real time implementation. In order to study the robustness of these controllers, a position disturbance that corresponds to a prompt and limited knee flexion is applied. In terms of position regulation, Fig. 15 shows that both controllers converged toward the desired position. In terms of input control, the PP controller is very sensitive to this disturbance while the MPC controller is much less sensitive and performs better with respect to external disturbances. Robustness of MPC with respect to parameter uncertainties has been also investigated. Fig. 16 shows the controller performances under parameter uncertainties such as inertial parameters, gravity center position and damping factor coefficient. It can be noticed that the MPC controller coupled with the IOFL succeeds to stabilize the knee joint output position to the desired value and shows consequently robustness with respect to parameter uncertainties (Fig. 17). The robustness of the original nonlinear system with respect to parameter uncertainties has not been proved yet theoretically.

5. Discussion This study presents a new approach of modeling the knee joint actuated by the quadriceps muscles through FES. The parameters

of the biomechanical knee–muscle system have been identified using different protocols and algorithms (Levenberg–Marquardt, linear least square, nonlinear interpolation, regression equations). Identification processes are conducted on both paraplegic and healthy subjects. The model parameters used in the illustrative simulations are taken for the healthy subject. It results in relatively high damping values of the knee joint, which are needed to avoid undesired oscillation in the movement. It should be noted that very high damping values may induce an alteration of the knee joint flexion–extension movements and will require consequently an increasing FES stimulation signal power. Performances of the MPC defined by Eqs. (19) and (20) and cascaded with the exact IOFL biomechanical model are presented in Section 4. The prediction and control horizon Np and Nu respectively are computed as a function of the system time constant: Np ¼30 and Nu ¼10. These parameters should not be assigned too large in order to avoid unnecessary solutions to the optimization problem for each sample. Alternatively, too small horizons may cause system instability. The constrained input u ¼ a is the recruitment variable since the recruitment curve is a static function and the pulse width as well as the stimulation current can be computed by simply inverting the recruitment curve. The recruitment variable is constrained between 0 and 1, representing respectively the cases of no fiber recruitment and maximal fiber recruitment. The controlled variable which is the controlled output is constrained to remain between y ¼ 01 (hyperextension) and y ¼ 901 (rest position). Only the knee position is recalled to update the

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

control input that corresponds to the desired trajectory. Performances of the MPC have been investigated by comparison to those of a classical pole placement controller, particularly, in terms of stability, regulation, finite time convergence as well as robustness with respect to external disturbances. The pole placement PP controller parameters k1, k2 and km are computed in order to ensure an exponential convergence toward zero of the tracking error as well as to have a better performance with respect to external perturbations. Performances of the PP controller would be better than the MPC in terms of overshoot and settling time if perturbations are not taken into consideration. Furthermore, the advantage of using MPC over PP controller is the capability of the MPC to add explicitly, in its formulation, constraints on states, inputs and outputs while these constraints could not be controlled if applying the PP controller to the biomechanical model.

6. Conclusions In this paper, a cascaded control strategy based on IOFL and MPC has been presented with reference to an existing experimental setup. The MPC offers the possibility to integrate explicitly into its formulation, constraints on input, output and measured states. IOFL ensures decoupling of the system dynamics into internal and external dynamics. While the external dynamics are simply controlled, stability of the system’s internal dynamics is mathematically proved. This controller is applied to a physiomathematical nonlinear muscle model actuating the knee joint for flexion/extension movements. Dynamic and geometric parameters of the biomechanical model were identified based on experimental data recorded using a video based motion analysis and a special experimental platform. Satisfactory stability and tracking error of the proposed strategy were achieved after a finite time delay relatively to a classical pole placement controller. Its performances have been assessed in the presence of external disturbances. The results show that the constraints on input and output were respected. The knee joint was chosen in this study because it is a key joint to be controlled in all the main movements of the lower limbs. Further experiments will be conducted to assess the knee joint control through three important movements: (i) standing up where the main issue is to have a controller that acts in synergy with the voluntary movement of the upper body, the knee joint provides an important amount of the torque (ii) quiet standing where the knee joint must be kept locked with the minimum stimulations levels to limit fatigue, (iii) assisted gait with a walker where the knee joint has to follow a desired trajectory.

Appendix A. Muscle contractile element stiffness The muscle contractile element stiffness as formulated in Eq. (7) is presented in the following: 8 L ¼ Ls þ Lc > > > > LL0 > > L ¼ L0 ð1 þ EÞ ) E ¼ > > L0 > < Ls Ls0 ðA:1Þ Ls ¼ Ls0 ð1 þ Es Þ ) Es ¼ > > Ls0 > > > > > Lc Lc0 > > : Lc ¼ Lc0 ð1þ Ec Þ ) Ec ¼ Lc0 Equations referred in (A.1) lead to

EL0 ¼ Ec Lc0 þ Es Ls0 ) EL0 ¼ Ec Lc0 þðLs Ls0 Þ

ðA:2Þ

193

Since serial element Es is in serial with the contractile element Ec, the force Fs developed by Es is equal to Fc. On the other hand, the force Fs is linked linearly to the serial stiffness Ks through the following formulation (El-Makssoud, 2005): F c ¼ F s ¼ K s ðLs Ls0 Þ ) Ls Ls0 ¼

Fc Ks

ðA:3Þ

Eqs. (A.2) and (A.3) lead to

EL0 ¼ Ec Lc0 þ

Fc L0 1 L0 1 _ Fc ) Ec ¼ E F c ) E_ c ¼ E_  Lc0 K s Lc0 K s Ks Lc0 Lc0

ðA:4Þ

Appendix B. Linearization The plant model in its simplified form can be expressed by the following set of differential equations: 8 x_ 1 ¼ a1 x1 þ a2 x1 x4 þ a3 Fðx3 Þu > > > > < x_ 2 ¼ a4 x1 x4 þ a1 x2 þ a2 x2 x4 þ a5 Fðx3 Þu ðB:1Þ x_ 3 ¼ x4 > > > > : x_ 4 ¼ a6 x2 þ a7 cosðx3 Þ þ a9 x3 þ a8 x4 A linearization of this nonlinear system around an operating point ðx,uÞ has been computed using the Jacobian of the system x_ ¼ f ðx,uÞ as well as its constraints. The linear model can be formulated as follows: ( x_ ¼ A:x þB:u y ¼ C:xþ D:u with 2 6 6 A¼6 6 4 2

a1 þ a2 x 4

0

a4 x 4

a1 þ a2 x 4

a3 u F_ ðx 3 Þ a5 u F_ ðx 3 Þ

0

0

0

0

a6

a9 a7 sinðx 3 Þ

a3 uFðx 3 Þ

a2 x 1

3

7 a4 x 1 þa2 x 2 7 7 7 1 5 a8

3

6 a uFðx Þ 7 3 7 6 5 B¼6 7, 4 5 0

C ¼ ½0 0 1 0,

D ¼ ½0

0

Appendix C. Input–output feedback linearization derivation C.1. The relative degree The relative degree of the nonlinear system (8) is computed using the Lie derivative (L) mathematic tool: 8 < Lg Li1 i A f1; 2, . . . ,r1g f hðxÞ ¼ 0, ðC:1Þ : Lg Lr1 hðxÞ a0 f with the relative degree r A N and 1 r r r n, 8x A Rn . The Lie derivatives of hðxÞ over f(x) and g(x) are defined as follows: Lf hðxÞ ¼

@h f ðxÞ @x

et

Lg hðxÞ ¼

and 8 > Lg Lif hðxÞ ¼ Lg ½Lif hðxÞ > > < Lif hðxÞ ¼ Lf ½Li1 f hðxÞ, > > > : L0 hðxÞ ¼ hðxÞ f

@h gðxÞ @x

i A f1; 2, . . . ,r1g

ðC:2Þ

194

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

The Lie derivatives of hðxÞ with respect to f ðxÞ and gðxÞ are given below: 8 > Lg L0f hðxÞ ¼ 0 > > > > > > L0f hðxÞ ¼ x3 > > > > > < L1f hðxÞ ¼ x4 ðC:3Þ > Lg L1f hðxÞ ¼ 0 > > > > 2 > > > Lf hðxÞ ¼ a6 x2 þ a7 cos x3 þa9 x3 þa8 x4 > > > > : Lg L2f hðxÞ ¼ a5 a6 Fðx3 Þ Lg L2f hðxÞ ¼ Lg Lr1 hðxÞ is a0 since a5 a 0, a6 a0 and Fðx3 Þ a 0 f (conditions always true in the knee range of motion). Consequently, r is equal to 3. C.2. Internal dynamics derivation and stability proof In the current study, the normal form (13) developed in Section 3.1 takes the following form: 8 z1 ¼ hðxÞ ¼ x3 > > < z2 ¼ Lf hðxÞ ¼ x4 ðC:4Þ > > : z3 ¼ L2f hðxÞ ¼ a6 x2 þ a7 cos x3 þ a9 x3 þ a8 x4 Since the system dimension n is equal to 4, the remaining equation of the transformation (C.4) corresponds to the system internal dynamics of dimension (nr) that is equal to 1. The internal dynamics, which are represented by the normal form of Eq. (C.4), must satisfy the following condition (Isidori, 1995): T   @Z @Z @Z @Z @Z @Z Lg Z ¼ gðxÞ ¼ 0 0 a3 Fðx3 Þ a5 Fðx3 Þ @x @x1 @x2 @x3 @x4 @x3 ¼

@Z @Z a3 Fðx3 Þ þ a5 Fðx3 Þ ¼ 0 @x1 @x2

ðC:5Þ

A solution to this equation is as follows:

ZðxÞ ¼

a5 x1 x2 a3

ðC:6Þ

It can easily be verified that this solution satisfies Eq. (C.5). Considering the state vector z ¼ ½z1 z2 z3 ZT , its Jacobian takes the form 3 2 3 2 0 0 1 0 z1 6 7 60 0 0 17 z 7 6 7 @ 6 7 6 27¼6 6 6 7 @x 4 z3 5 4 0 a6 a9 a7 sin x3 a8 7 5 a5 1 0 0 Z a3 This matrix is not singular for every x, its rank is equal to 4 (a3 ,a5 ,a7 and a8 are not null and x3 4arcsinða9 =a7 Þ). The inverse transformation of (C.4) is given by 8 a3 > > > x1 ¼ a a ½a6 Z þ z3 a7 cos z1 a9 z1 a8 z2  ¼ dðzÞ > > 5 6 > > < 1 x2 ¼ ½z3 a7 cos z1 a9 z1 a8 z2  ¼ oðzÞ ðC:7Þ a6 > > > > x3 ¼ z1 > > > :x ¼z 4

2

Consequently, this transformation is available and the system dynamics can be represented in the normal form where the system output and all its derivatives with respect to time will be constrained to zero. The internal dynamics are then represented by   a Z_ ¼ a1 Z þ 5 a2 a4 z2 dðzÞ ðC:8Þ a3 The zero dynamics z1 ¼z2 ¼z3 ¼0 are described by

Z_ ¼ a1 Z

ðC:9Þ

Since a1 ¼ uch is always negative, the global system presents stable zero dynamics that converge exponentially toward zero and consequently, toward stable internal dynamics.

References Akesson, J. (2003). Operator interaction and optimization in control systems. Ph.D. Thesis. Department of Automatic Control Lund Institute of Technology Lund. Allgower, F., Badgwell, T. A., Qin, S. J., Rawling, J. B., & Wright, S. J. (1999). Nonlinear predictive control and moving horizon estimation—An introductory overview. in: P. Frank (Ed.), Advances in control: Highlights of ECC 99 (pp. 391–449). Springer-Verlag. Alon, G., Levitt, A. F., & McCarthy, P. A. (2007). Functional electrical stimulation enhancement of upper extremity functional recovery during stroke rehabilitation: A pilot study. Neurorehabilitation and Neural Repair, 21, 207. Azevedo, C., Poignet, P., & Espiau, B. (2002). On line optimal control for biped robots. In International federation of automatic control (IFAC) (pp. 749–757), Barcelona, Spain. Bajd, T., Kralj, A., Sega, J., Turk, R., Benko, H., & Strojnik, P. (1981). Use of a twochannel functional electrical stimulator to stand paraplegics. Physical Therapy, 61, 526–527. Bestel, J., & Sorine, M. (2000). A differential model of muscle contraction and applications. In Schloessmann seminar on mathematical models in biology, chemistry and physics. Bad Lausick, Germany: Max Plank Society, May 19–23. Bhadra, N., & Peckham, P. H. (1997). Peripheral nerve stimulation for restoration of motor function. Journal of Clinical Neurophysiology, 14, 378–393. Braz, G. P., Russold, M., & Davis, G. M. (2009). Functional electrical stimulation control of standing and stepping after spinal cord injury: A review of technical characteristics. Neuromodulation, 12, 180–190. Camacho, E., & Bordons, C. (1995). Model predictive control in the process industry. New York: Springer. Chen, C. C., He, Z. C., & Hsueh, Y. H. (2009). An EMG feedback control functional electrical stimulation cycling system. Journal of Signal Processing Systems, 64(2), 195–203. DeLeva, P. (1996). Adjustments to Zatsiorsky–Seluyanov’s segment inertia parameters. Journal of Biomechanics, 29, 1223–1230. Donaldson, N., & Yu, C. (1996). FES standing control by handle reactions of leg muscle stimulation (CHRELMS). IEEE Transactions on Rehabilitation Engineering, 4, 280–284. El-Makssoud, H. E. (2005). Modelling and identification of the skeletal muscle under functional electrical stimulation. Ph.D. Thesis. LIRMM-University of Montpellier II. El-Makssoud, H. E., Guiraud, D., & Poignet, P. (2004). Mathematical muscle model for electrical stimulation control strategies. In IEEE international conference on robotics and automation ICRA’04 (pp. 1282–1287), New Orleans, LA, USA. Ferrarin, M., Palazzo, F., & Riener, R. (2001). Model-based control of FES-induced single joint movements. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 9, 245–257. Freeman, C., Davies, I., Lewin, P., & Rogers, E. (2008). Iterative learning control of upper limb reaching using functional electrical stimulation. In The 17th IFAC World congress (pp. 13444–13449), Seoul, Korea, July 6–11. Freeman, C., Hughes, A. M., Burridge, J., Chappell, P., Lewin, P., & Rogers, E. (2009). Iterative learning control of FES applied to the upper extremity for rehabilitation. Control Engineering Practice, 17, 368–381. Gautier, M., & Poignet, P. (2002). Closed loop identification by inverse model of physical parameters of mechatronic systems. Journal Europe´en des Syste mes Automatise´s, JESA, 36, 465–480. Gill, P. E., & Murray, W. (1978). Algorithms for the solution of the nonlinear leastsquares problem. SIAM Journal on Numerical Analysis, 15(5), 977–992. Guemghar, K. (2005). On the use of input–output feedback linearization techniques for the control of nonminimum-phase systems. Ph.D. Thesis. Switzerland: Ecole Polytechnique Federale de Lausanne, EPFL. Guiraud, D., Stieglitz, T., Koch, K. P., Divoux, J. L., & Rabischong, P. (2006). An implantable neuroprosthesis for standing and walking in paraplegia: 5-year patient follow-up. Journal of Neural Engineering, 3, 268–275. Hawkins, D., & Hull, M. (1990). A method for determining lower extremity muscletendon lengths during flexion/extension movements. Journal of Biomechanics, 23, 487–494. Isidori, A. (1990). Nonlinear control systems ((2nd ed.). Springer-Verlag. Isidori, A. (1995). Nonlinear control systems ((3rd ed.). Springer-Verlag. Kobetic, R., & Marsolais, B. (1994). Synthesis of paraplegic gait with multichannel functional neuromuscular stimulation. IEEE Transactions on Rehabilitation Engineering, 2(2), 66–79. Kobravi, H. R., & Erfanian, A. (2009). Decentralized adaptive robust control based on sliding mode and nonlinear compensator for the control of ankle movement using functional electrical stimulation of agonist antagonist muscles. Journal of Neural Engineering, 6, 046007. Kromer, V. (1993). Muscular forces analysis during locomotion—Rigid body approach and simulation in flexible plans mechanics by finite element. Ph.D. Thesis. France: Institut National Polytechnique de Lorraine. Le, F., Markovsky, I., Freeman, C., & Rogers, E. (2010). Identification of electrically stimulated muscle models of stroke patients. Control Engineering Practice, 18, 396–407.

S. Mohammed et al. / Control Engineering Practice 20 (2012) 182–195

Mohammed, S. (2006). Contribution to movement synthesis and control of skeletal muscle under functional electrical stimulation. Ph.D. Thesis. LIRMM-University of Montpellier II. Mulder, A., Veltink, P., & Boom, H. (1992). On/off control in FES-induced standing up: A model study and experiments. Medical and Biological Engineering and Computing, 30, 205–212. Nevistic, V., & Morari, M. (1996). Robustness of MPC-based schemes for constrained control of nonlinear systems. In International federation of automatic control (IFAC) congress (pp. 25–30), San Francisco, USA. Poboroniuc, M., Wood, D., Donaldson, N., Fuhr, T., & Riener, R. (2003). Closed-loop control for FES-based restoration of standing in paraplegia. In 2nd World congress of the international society of physical and rehabilitation medicine—ISPRM (pp. 201–204), Prague, Czech Republic. Poboroniuc, M. S., Fuhr, T., Wood, D., Riener, R., & Donaldson, N. (2002). FESinduced standing-up and sitting down control strategies in Paraplegia. In FESnet conference (pp. 1–3), Glasgow, UK. Previdi, F., Ferrarin, M., Savaresi, S. M., & Bittanti, S. (2005). Closed-loop control of FES supported standing up and sitting down using virtual reference feedback tuning. Control Engineering Practice, 13, 1173–1182. Richalet, J., Ata-Doss, A., Arber, C., Kuntze, H. B., Jacubasch, A., & Schill, W. (1996). Predictive functional control: Application to fast and accurate robots. In International federation of automatic control (IFAC) congress (Vol. 4, pp. 251–258), Munich, Germany.

195

Riener, R., & Fuhr, T. (1998). Patient-driven control of FES-supported standing up: A simulation study. IEEE Transactions on Rehabilitation Engineering, 6, 113–123. Salmons, S. (1994). Exercise, stimulation and type transformation of skeletal muscle. International Journal of Sports Medicine, 15, 136–141. Schauer, T., & Hunt, K. J., 2000. Nonlinear predictive control of knee-joint angle using FES. In 5th annual conference of the international functional electrical stimulation society IFESS (pp. 425–428), Aalborg, Denmark. Shen, X. (2009). Nonlinear model-based control of pneumatic artificial muscle servo systems. Control Engineering Practice, 18, 311–317. Slotine, J., & Weiping, L. (1991). Applied nonlinear control. Prentice Hall. Veltink, P., Chizeck, H., Crago, P., & El-Bialy, A. (1992). Nonlinear joint angle control for artificially stimulated muscle. IEEE Transactions on Biomedical Engineering, 39, 368–380. Vette, A. H., Masani, K., Kim, J. Y., & Popovic, M. R. (2009). Closed-loop control of functional electrical stimulation-assisted arm-free standing in individuals with spinal cord injury: Feasibility study. Neuromodulation: Technology at the Neural Interface, 12, 22–32. Wood, D., Harper, V., Barr, F., Taylor, P., Phillips, G., & Ewins, D. (1998). Experience in using knee angles as part of a closed-loop algorithm to control FES-assisted paraplegic standing. in: 6th international workshop on functional electrical stimulation (pp. 137–140). Zahalak, G. I. (1981). A distribution-moment approximation for kinetic theories of muscular contraction. Mathematical Biosciences, 55, 89–114.