Iterative learning control of FES applied to the upper extremity for rehabilitation

ARTICLE IN PRESS Control Engineering Practice 17 (2009) 368–381

Contents lists available at ScienceDirect

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

Iterative learning control of FES applied to the upper extremity for rehabilitation C.T. Freeman a,, A.-M. Hughes b, J.H. Burridge b, P.H. Chappell a, P.L. Lewin a, E. Rogers a a b

School of Electronics and Computer Science, University of Southampton, Southampton SO17 1BJ, UK School of Health Sciences, University of Southampton, Southampton SO17 1BJ, UK

a r t i c l e in f o

a b s t r a c t

Article history: Received 1 August 2007 Accepted 29 August 2008 Available online 21 October 2008

Iterative learning control schemes are used to apply functional electrical stimulation to the triceps of unimpaired subjects in order to perform trajectory tracking tasks. The subjects supply no voluntary effort and a robotic workstation is used to constrain their movement, impose known dynamics at the point of interaction with the robot, and provide assistive torque about the shoulder. Results from 18 subjects are presented and show that a high level of performance can be achieved using the proposed method. In addition to illustrating how stimulation and robotics can be successfully combined in order to perform reaching tasks, the results provide justification for the system to be subsequently used by stroke patients for rehabilitation. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Iterative learning control Functional electrical stimulation Robotics Tracking control Muscle models

1. Introduction Strokes affect between 174 and 216 people per 100,000 population in the UK each year (Mant, Wade, & Winner, 2004). Since around two-thirds of patients in England survive their stroke, there are more than 900,000 people who have had a stroke living in England of which 50% will be left disabled and dependent (Hudson, Scharaschkin, Wilkins, & Taylor, 2005). Half of all acute stroke patients starting rehabilitation will have a marked impairment of function in one arm of whom only about 14% will regain useful sensory-motor function (Royal College of Physicians, 1989). The stroke patient is unable to practice arm movement because of impaired motor control and the delay in recovery may lead to a decreased likelihood of recovery occurring at all (Castro-Alamancos, Garcia-Segura, & Borrell, 1992). Functional electrical stimulation (FES) can provide the experience of moving for the patient, and consequently may limit the problem of inability to practice movement because of impaired motor control. It has been used successfully to improve recovery of upper limb motor control (Burridge & Ladouceur, 2001) and recent studies have shown that when stimulation is associated with a voluntary attempt to move the limb, improvement is enhanced (Burridge & Ladouceur, 2001). A hypothesis has been proposed to explain this enhanced recovery observed when stimulation is applied to correspond with a subject’s intended

 Corresponding author.

E-mail address: [email protected] (C.T. Freeman). 0967-0661/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.conengprac.2008.08.003

movement (Rushton, 2003), which suggests that the degree of functional recovery is closely related to the accuracy with which the stimulation assists the subject’s completion of a repeated task. In this paper, the task will be to repeatedly track a trajectory which is predefined so that the subject’s voluntary intention is exactly known. This then allows FES controllers to precisely assist their desired movement in order to promote the therapeutic effect gained over repeated performance of the tracking task. The trajectories used will be elliptical and will approximate reaching movements commonly encountered in everyday life. During treatment they will be tracked using the impaired arm, and are designed to target the difficulty most stroke patients have in performing full elbow extension. The exact sequence is that the patient has to try to complete the movement, their arm is then returned to the starting position, and the task repeated. During each repetition of the tracking task, the performance can be measured and a natural approach is to use this information to adjust the stimulation to produce improved performance on the next (and successive) attempts. This scenario fits exactly within the iterative learning control (ILC) setting for processes which are required to continually repeat the same task over a finite interval with resetting between trials. The salient feature of ILC is to use information from previous trials to update the current trial input so that performance is successively improved, making it a natural candidate for the stroke rehabilitation task considered here. ILC has been able to provide high levels of accuracy in the presence of significant plant uncertainty in a wide variety of applications (see Ahn, Chen, & Moore, 2007, for a comprehensive

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

review). Prior to the application of ILC, a feedback controller is typically first implemented to act as a prestabilizer. Using this arrangement, even simple ILC structures have been shown to provide significant improvement in tracking performance compared with the use of feedback controllers alone (Ratcliffe, 2005). Using data from tests conducted on stroke patients, it has been shown that the presence of voluntary effort occurring during the application of electrical stimulation to the arm can be accounted for by an additive disturbance acting on the stimulation input signal (see Freeman et al., 2008a, for details). Furthermore, if the stimulation does not change significantly over repeated tracking tasks (as can be ensured using ILC), this disturbance has been seen to be reasonably small and to change only slowly between trials. Therefore a simple approach to incorporating voluntary effort is to design algorithms which exhibit sufficient robustness to operate in the presence of this form of disturbance appearing on the plant input. It has been shown that ILC provides such a framework (see, for example, Moore, 1992). The initial step is therefore to confirm the performance of controllers in the absence of voluntary effort (the results gained in this process also providing evidence that the controller is also appropriate for use by paralysed subjects in order to perform accurate reaching movements necessary for daily living). Although ILC is an obvious choice for the current application, it is possible that a suitable method exists among those already applied in this area. Unfortunately, open-loop methods for the control of FES for upper limb movements (see for example Davoodi & Andrews, 2004; Popovic & Popovic, 1998) have not been able to provide the high level of performance which is necessary to fully promote the required association between the subject’s intended movement, and the action of the applied electrical stimulation in realizing it. A wide variety of modelbased controllers for electrical stimulation have been implemented which may be able to produce such accuracy. These include the use of H1 (Hunt, Jaime, & Gollee, 2001), optimal control (Hunt, Munih, & Donaldson, 1997) and fuzzy logic control (Davoodi & Andrews, 1998) for paraplegic standing, and the use of neural networks (Hatwell, Oderkerk, Sacher, & Inbar, 1991) and data-driven control (Previdi, Ferrarin, Savaresi, & Bittanti, 2005) for control of the paraplegic knee joint. Few such model-based schemes have been applied to upper limb movement, exceptions principally comprising the use of neural nets for paraplegic arm movements (see for example Lan, Feng, & Crago, 1994; Tresadern, Thies, Kenney, Howard, & Goulermas, 2006). Unfortunately these require extensive training and have unresolved stability issues due to their black-box structure. Another factor is that any suitable control method selected for the stimulation must also operate in the presence of voluntary effort supplied by the patient. A simple method of achieving this is for the controller to directly use electromyographic (EMG) or myoelectric activity of the muscle being stimulated (see for example Thorsen, Spadone, & Ferrarin, 2001). However, model-based control methods have not yet incorporated this information since it does not directly relate to the force or torque generated by the muscle, and because the signal is often either weak and unreliable or that the artefact produced by the stimulation signal corrupts the natural EMG signal (although in this case blanking techniques may be applied). ILC is one of the few advanced controllers that has previously been applied to control the upper limb, although a high level of performance has not been achieved in practice (Dou, Tan, Lee, & Zhou, 1999). The motivation for this paper is therefore the lack of controllers for upper limb electrical stimulation that are capable of providing a high level of performance, but which do not require extensive tuning. A future concern is then the need for the controller to operate whilst voluntary effort is simultaneously supplied by the subject.

369

The need for a controlled environment in which to apply stimulation has inspired the design and construction of an experimental test facility incorporating a five-link planar robotic arm and an overhead trajectory projection system (see Freeman et al., 2008b, for further details). The subject is seated with their arm strapped to the robot, and the trajectory is projected onto a target mounted above their hand. Their task is to repeatedly track it using a combination of voluntary control and surface FES applied to muscles in their impaired shoulder and arm. During the tests carried out, the robotic arm provides control over the dynamics experienced by the patient, and produces an assistive torque when necessary. Whilst there exist several robotic devices for the application of robotic therapy to stroke patients through purely mechanical manipulation of their arms (e.g. the MIT Manus system, Krebs, Hogan, Aisen, & Volpe, 1998), this form of treatment has hitherto not been combined with the application of FES. In this paper a dual control scheme is formulated to govern both the application of stimulation and the assistance provided by the robotic arm. Tests are conducted on 18 unimpaired subjects who do not provide voluntary effort. The results illustrate how stimulation and robotics can be successfully combined in order to perform reaching tasks, and confirm the efficacy of the system prior to its use by stroke patients. This study has received ethical approval (reference no. 505-12/1).

2. Workstation description The robotic workstation is shown in Fig. 1 and consists of a five-link planar robotic arm rigidly connected to an overhead projection system. The links of the robot are labelled 1–5, and the upper arm and forearm are labelled u and f, respectively. The vectors x0 , y0 and z0 are components of the robotic base coordinate frame, and x1 , y1 and z1 are those of the human arm base coordinate frame, the two systems being related by a translation. Two coaxially mounted DC brushless motors actuate links 1 and 2, and a 4000-line encoder is mounted on each motor shaft to record the angle of these links. The subject is strapped to the extreme fifth link, and grips a cushioned handle which is rigidly connected to a 6 axis force/torque sensor which records the force they apply to the robotic end effector. Forces can be measured of up to 200 N applied in the horizontal plane with a resolution of 0.0122 N. The fifth link also contains a 4000-line encoder to measures its angle, and LEDs to provide visual feedback of the tracking performance. The robotic arm is used

Fig. 1. Unimpaired subject using the robotic workstation.

ARTICLE IN PRESS 370

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

Fig. 2. Geometry of combined human arm and robotic arm system.

to constrain the subject’s arm, to impose forces on the endeffector that make the task feel ‘natural’ to the subject, to apply assistance during the performance of tracking tasks, and to move the patient’s arm when necessary. The subject’s task during the treatment is to track a range of trajectories that are projected onto a target positioned above their hand. All tests are controlled using a graphical user interface running on the host PC. Further details of the design and functionality of the workstation, and the instrumentation it contains, can be found in Freeman et al. (2008b). A description of the identification procedures necessary to produce an accurate model of the robot for use in the controller is also provided. Fig. 2 shows the geometry of the dual human and robotic system. The robot joint angle vector is qr ¼ ½W1 W2 T , where W1 and W2 are the joint angles corresponding to links 1 and 2, respectively. The axes of all the robot joints are parallel to z0 , and the links are labelled to correspond with Fig. 1. The torque supplied by the motors is given by sr ¼ ½t1 t2 T, where t1 and t2 are applied to actuate joints 1 and 2, respectively. It is assumed that the subject interacts with the robot by applying a vector of forces and torques at the point Q, which has a z0 component of lz . This corresponds to the vertical distance of the subject’s hand above the force/torque sensor. The vector of the components of these forces applied in the x0 and y0 directions is given by hr . The subject’s arm is strapped to the fifth link of the robot and the human arm model therefore includes the properties of this link, which will be neglected from the model of the robot. To ensure the robot’s safe interaction with an unknown environment, a form of impedance control (see Colgate & Hogan, 1998, for details) is used to govern the torque demand supplied to the motors. This is a strategy that has previously been used to control robotic therapy devices (see Krebs et al., 1998), although alternative approaches could be used (such as force or compliant control Schutter & Brussel, 1988; Siciliano & Villani, 1999). The controller compensates for the inertial and damping properties of the robotic arm, and produces the relationship hr ¼ K K x x~ r þ K Bx x_~ r þ K Mx x€~ r

(1)

at the point Q. Here x^ r is the reference position, x~ r ¼ x^ r  xr , xr ¼ kr ðqr Þ, x_ r ¼ J r ðqr Þq_ r and x€ r ¼ J r ðqr Þq€ r þ J_ r ðqr ; q_ r Þq_ r . The direct kinematics equation for the robotic system is xr ¼ kr ðqr Þ, and J Tr ðqr Þ is the Jacobian of that system. When the robot is moved freely by the subject in the absence of assistance, the gain matrices are set as K K x ¼ 0, K Bx ¼ K Bx I and K Mx ¼ K Mx I, where I is the identity matrix. The values of K Bx and K Mx assume positive values and are tuned to create a ‘natural’ feel, and an additional requirement is that x^ r equals a constant value. When the robot is required to move the subject’s arm along predefined trajectories, it is necessary to set K K x ¼ K K x I with K K x 40. The three gains are then tuned to produce the required tracking performance. The

form of the gain matrices for the case where the robot applies assistance during tracking tasks is described in Section 4. Further details of the robotic controller are given in Freeman et al. (2008b).

3. Human arm model In this section a mathematical description of the subject’s arm is derived, consisting of a model of the passive dynamical system to which the torque generating properties of the stimulated muscle is then added. 3.1. Passive system Fig. 3 provides a more detailed description of the geometry of the constrained human arm model that appeared in Fig. 2. The first link represents the upper arm, from the shoulder joint to the elbow, with length ðlu1 þ lu2 Þ. The second link represents the forearm, from the elbow to the thumb web, with length ðlf 1 þ lf 2 Þ. The constraint means that the forearm must lie in the horizontal plane, and rotation is possible about the axis along the upper arm. The point, Q, denotes where the subject’s hand grasps the robot, and only forces and torques along unconstrained directions are shown (only rotation about the axis parallel to z1 is unconstrained at this point, rotation not being possible about the two orthogonal axes). The triceps has been selected for stimulation since stroke patients typically experience problems with shoulder and elbow extension during reaching tasks (Lum, Burgar, & Shor, 2004). Actuation of the triceps has been modelled as a torque, T b X0, acting about an axis orthogonal to both the upper arm and forearm. Components of the forces in the x1 and y1 directions applied by the subject’s hand at the point of interaction with the robot are denoted by F x1 and F y1 , respectively. A point, xr , in the robot coordinate system is expressed in the human arm system by a point, xa , such that xr ¼ xa þ ½lx ly lz T

(2)

where lx and ly represent the x0 and y0 components of the shoulder joint position in the robot coordinate system. If V u and V f are vectors aligned along the upper arm and forearm, respectively, the horizontal constraint requires that V f has a z1 component of zero. In this case the unitary axis about which T b is applied is given by Ve ¼

Vu  Vf 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½suf sg cuf sg  sf cg T jV u  V f j 1  c2 c2 f

(3)

g

where cf and cg denote cosðWf Þ and cosðgÞ, respectively, and similarly sf and sg denote sinðWf Þ and sinðgÞ. In addition cu and su are used to denote cosðWu Þ and sinðWu Þ, respectively, and suf and cuf

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

371

passive properties of the joint. Since g is invariant, (6) means F mp ðb; b_ Þ is accounted for when using the most general form of F a ðqa ; q_ a Þ that has been considered. It has therefore been decided not to add additional forms of friction to the existing model. Details of the procedures used to establish the parameters appearing in the model, and also additional information relating to the subjects whose experimental results are presented in Section 6, are given in Freeman et al. (2008a).

4. Robotic assistance and trajectory choice The action of the robotic arm necessary to make the task a feasible yet productive one is described in this section. The following points therefore concern the choice of trajectory and role of the robot during the performing of the task:

Fig. 3. Geometry of constrained human arm.

denote sinðWu þ Wf Þ and cosðWu þ Wf Þ, respectively. For arbitrary g; V e has a z1 component of zero when Wf ¼ 0 and is equal to V^ e ¼ ½su cu 0T

(4)

As Wf increases, the angle that V e rotates about the upper arm, starting from V^ e , is then given by 0 1 ! ^ c s Ve  Ve f g B C ¼ arccos@qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA (5) aðWf Þ ¼ arccos jV e jjV^ e j 1  c2f c2g The elbow angle is equal to     Vu  Vf ¼ arccos cf cg bðWf Þ ¼ arccos jV u jjV f j

(6)

(i) The trajectories will be elliptical reaching tasks for each subject’s dominant arm, and should be achievable given their identified arm model. (ii) The triceps muscle will provide the sole actuating torque about the elbow, and the robotic arm will use the control scheme given by (1) to make the dynamics about the elbow feel ‘natural’ to the subject. (iii) The robotic arm will provide a torque acting about the subject’s shoulder in order to track the reference in a manner which is entirely governed by the angle of the forearm. This then makes the task feasible without lessening the role played by the triceps. The robotic control scheme needed to achieve these goals is derived in the next subsection.

The dynamic model of the constrained arm can now be expressed in the form

4.1. Robotic control scheme

(7) Ba ðqa Þq€ a þ C a ðqa ; q_ a Þq_ a þ F a ðqa ; q_ a Þ ¼ sa  J Ta ðqa Þha h qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiiT where qa ¼ ½Wu Wf T , sa ¼ 0  T b sf cg = 1  c2f c2g , ha ¼ ½F x1 F y1 T

Firstly the human arm model, described by (7), and the endeffector dynamics, given by (1), are combined to give the system Ba ðqa Þq€ a þ C a ðqa ; q_ a Þq_ a þ F a ðqa ; q_ a Þ ¼ s þ J T ðq ÞðK x~ þ K x~_ þ K

and the remaining terms are given in Appendix A. The relationship ha ¼ hr is a consequence of the connection between the robotic and human arm systems. Various forms of the friction term F a ðqa ; q_ a Þ have been considered in order to arrive at a compromise between repeatability and the accuracy of the overall model. The most general form considered is

in which x^ a ¼ x^ r  ½lx ly lz T has been set to produce x~ a ¼ x~ r . To separate the dynamics of the end effector in the directions corresponding to the human arm joint angles, it is then necessary to set

F a ðqa ; q_ a Þ ¼ ½F a1 ðWu ; W_ u Þ F a2 ðWf ; W_ f ÞT

K K x x~ a þ K Bx x_~ a þ K Mx x€~ a

(8)

a

Kx a

a

a

Bx a

€

~ aÞ Mx x

_~ €~ ~ ¼ J T a ðqa ÞðK K q qa þ K Bq qa þ K M q qa Þ

in which F a1 ðÞ and F a2 ðÞ are piecewise linear functions.

where q~ a ¼ q^ a  qa and q^ a ¼ 3.2. Muscle model

1 ka ðx^ a Þ.

(10)

(11) This produces

Ba ðqa Þq€ a þ C a ðqa ; q_ a Þq_ a þ F a ðqa ; q_ a Þ ¼ s þ K q~ þ K q_~ þ K q€~

(12)

To account for the action of the triceps, an established model of the torque, T b , generated by electrically stimulated muscle acting about a single joint is given by

To meet the requirements of points (2) and (3) above, it is also required that

T b ðb; b_ ; u; tÞ ¼ gðu; tÞ  F ma ðb; b_ Þ þ F mp ðb; b_ Þ

K K q ¼ diagfK K 1 ; 0g;

(9)

where u denotes the stimulation pulsewidth applied, and b is the joint angle (see, for example Shue, Crago, & Chizeck, 1995). A Hammerstein structure incorporating a static non-linearity, hIRC ðuÞ, representing the isometric recruitment curve, cascaded with linear activation dynamics, hLAD ðtÞ, produces the first term, gðu; tÞ. The activation dynamics can be adequately captured using a critically damped second order system (Baratta & Solomonow, 1990). The term F ma ðb; b_ Þ models the multiplicative effect of the joint angle and joint angular velocity on the active torque developed by the muscle. The term F mp ðb; b_ Þ accounts for the

a

Kq

Bq

a

Mq

a

a

K Bq ¼ diagfK B1 ; K B2 g

K Mq ¼ diagfK M1 ; K M2 g

(13)

and q^ a ¼ ½W^ u cT where c is a constant and K K 1 ; K B1 ; K B2 ; K M1 ; K M2 X0. This allows a choice of arbitrary second order dynamics to be imposed about the shoulder and the damping and inertia about the elbow to be prescribed. This produces the expression 2 3 _ € K K 1 W~ u þ K B1 W~ u þ K M1 W~ u 4 5 (14) sa þ K W_  K W€ B2

f

M2

f

ARTICLE IN PRESS 372

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

for the right-hand side of (12), and provides the required dynamic relationship for both components of the torque. However, to achieve this with the proposed control scheme has involved the _ € assumption W^ f ¼ c, so that W^ f ¼ W^ f ¼ 0. This is unsuitable for the desired references, and so the control scheme is not appropriate _ for the intended tracking task. If instead it is assumed that W^ u and €^ Wu are sufficiently small, then the controller can still provide asymptotic stability to externally applied forces by using 2 3 K K 1 W~ u  K B1 W_ u  K M1 W€ u 5 (15) sa þ 4 K B2 W_ f  K M2 W€ f h iT in place of (14). This can be produced by choosing q^ ¼ W^ u W^ f

a

and by exchanging (1) for hr ¼ K K x x~ r  K Bx x_ r  K Mx x€ r

(16)

The requirement given by (11) is now replaced by K K x x~ a  K Bx x_ a  K Mx x€ a ~ _ € ¼ J T a ðqa ÞðK K q qa  K Bq qa  K Mq qa Þ

(17)

To arrive at the required values of W^ u and W^ f , components of (17) are compared to give " # K K 1 ðW^ u  Wu Þ ðq Þ K K x ðx^ a  xa Þ ¼ J T a a 0 " # ^ K K 1 ðWu  Wu Þ cuf ¼ (18) ðlu1 þ lu2 Þcg sf suf This leads to a solution K Kx ¼

K K 1 ðW^ u  Wu Þ I jx^ a  xa jðlu1 þ lu2 Þcg sf

(19)

and "

cuf x^ a ¼ xa þ jx^ a  xa j suf

# (20)

so that x^ a is not uniquely defined, but can be any point lying on a line extending along the forearm and passing through xa . To achieve the tracking task it must therefore be set equal to the point of intersection with the trajectory. This is shown in Fig. 4, in which xa ðtÞ ¼ ka ðqa ðtÞÞ, with the trajectory defined by 



qa ðtÞ ¼ ½Wu ðtÞ Wf ðtÞT ;

t 2 ½0 T

(21)

Eliminating t from the components provides the relationship Wu ¼ CðWf Þ. The reference point is then defined formally as 02 31 B6 C ^ 7 B6 7C x^ a ¼ Xðxa ; CðÞÞ:¼ka B6 CðWf Þ 7C @4 5A W^ f  02 31 " #  cuf CðW^ f Þ  ka @4 5A ¼ xa þ l  suf  W^

Fig. 4. Trajectory geometry.

(ii) With the use of robotic assistance, it is shown in Appendix C that the behaviour of the electrically stimulated forearm can be approximated by the system represented schematically in Fig. 6. The approximations are based on experimentally confirmed properties of the human arm model, and the ability of the robotic control system to provide accurate tracking of W^ u by Wu.  Given T b X0, point (1) necessitates that Wf ðtÞ be monotone decreasing. With reference to the geometry of the task shown in Fig. 4, this requires that the line joining xa ðTÞ and xa ð0Þ pass through the origin of the arm system, since at both these points the direction of the elliptical trajectory is orthogonal to the major axis.  (iii) Under the assumption that Wf ðtÞ is tracked perfectly, the torque that must be applied to the arm system shown in Fig. 6 is given by 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1  1  c2f  ðtÞ c2g  AfK B W_  ðtÞ T b ðtÞ  @ f 2 sf  ðtÞ cg 

(22)

where l is a scalar. The complete control system is shown in Fig. 5, and the necessary controller gains are provided in Appendix B. 4.2. Trajectory selection Having established the robotic control system, certain points will now be made concerning the nature of the trajectories used and the selection of K B2 and K M2 that govern how the task feels to the subject: 

(i) For CðWf Þ to be a one–one continuous function, both Wu ðtÞ and Wf ðtÞ must be monotone.



(23)

and the stimulation required in order to produce it, u ðtÞ, must satisfy Z t T b ðtÞ hIRC ðu ðtÞÞhLAD ðt  tÞ dt ¼ (24)   0 F ma ðb ðtÞ; b_ ðtÞÞ 

f



þ ðK M2 þ ba3 ÞW€ f ðtÞ þ F a2 ðWf ðtÞ; W_ f ðtÞÞg

where b ðtÞ corresponds to the reference trajectory, so that b ðtÞ ¼ bðWf ðtÞ; tÞ. This expression clearly limits the magnitude and rate of change of any achievable torque trajectory. In practice the existence of a solution to (24) has been ensured by selecting slow trajectories that comprise half ellipse segments that are comfortably within both the robot’s and the subject’s workspace. The start and end points are chosen so that they can be reached by a smooth extension of the elbow, and are individually calculated for each subject depending on their maximum reach capability. The gains, K B2 and K M2 , are selected in order to mimic a realistic activity, to provide a high level of stability in response to sudden stimulation inputs, and to necessitate a moderate level of work to be generated by the subject’s muscles in order to track the trajectory. Their effect on the overall control scheme is discussed in Section 5.1.

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

373

Fig. 5. Human arm system with robotic assistance.

Fig. 6. Approximate model of stimulated human arm.

controller is thus given by

5. Control schemes

0 Two strategies have been selected to control the stimulation applied to the subject. The first consists of a linearizing controller in a simple feedback arrangement, the second involves the addition to this of an ILC feedforward controller, two types of which are considered. The controllers have been tested on 18 unimpaired subjects aged between 50 and 65 years (mean age 57 years, 10 months with standard deviation 5 years, 4 months). These subjects are age matched with the stroke patients on whom the controllers will later be applied. Representative results from a single subject are presented to enable detailed examination of performance. In addition, summary results of the performance of all subjects are included. 5.1. Linearizing controller

1

u ¼ hIRC

@ðF a2 ðWf ; W_ f Þ þ wÞ

to address the multiplicative effect of the first two non-linear terms, and then F a2 ðWf ; W_ f Þ is added to the input in order to cancel the additive effect of the third non-linear term. The linearizing

(26)

and the validity of this approach will now be investigated. For conciseness, attention is restricted to the case where F a2 ðÞ and F ma ðÞ are functions of their first argument only. The resulting system is shown in Fig. 7 in which gðWf Þ ¼

sf cg F ma ðbðWf ÞÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  c2f c2g

(27)

and P denotes multiplication of inputs. In Appendix D an equivalent state-space expression for this system is derived, and then linearized at time t¯ (see (D.4)). When the first derivatives of F a2 ðWf Þ and gðWf Þ are zero for all t, this linearized system corresponds to the desired relationship

1

The first component of the linearizing controller is hIRC ðÞ, which is the inverse of the isometric recruitment curve that has been identified for each subject. The controller’s subsequent action is thenq motivated ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi by the form of the remaining non-linear terms sf cg = 1  c2f c2g , F ma ðb; b_ Þ and F a2 ðWf ; W_ f Þ appearing in Fig. 6. The value of all three functions can be shown to vary only slowly when the trajectories considered in this paper are followed perfectly (see Freeman et al., 2008a for details). This can be seen by examining graphs of each function using the substitutions   Wf ¼ Wf ðtÞ, W_ f ¼ W_ f ðtÞ, b ¼ b ðtÞ and b_ ¼ b_ ðtÞ. The control action taken is to attempt to remove the effect of these functions in order to produce a system which approximates the linear activation dynamics in series with the linear arm dynamics (the transfer functions appearing in the left and right sub-systems, respectively, of Fig. 6). To achieve this, the controller next applies the gain term qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  1  c2f c2g (25) ¼ ðsf cg = 1  c2f c2g Þ1 ðF ma ðb; b_ ÞÞ1 sf cg F ma ðb; b_ Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  1  c2f c2g A sf cg F ma ðb; b_ Þ

Wf w

ðsÞ ¼

s2

o2n 1 ¼ GðsÞ  þ 2son þ o2n sððba3 þ K M2 Þs þ K B2 Þ

(28)

So that (28) adequately approximates the system behaviour when the derivatives are non-zero, trajectories must be chosen over which g 1 ðtÞ, g 2 ðtÞ and F 0a2 ðÞ are small. Since the terms in (D.4) involving g 1 ðt¯ Þ; g 2 ðt¯ Þ also include multiplication by C p Ap, their effect can also be reduced by ensuring that the dynamics of the arm system are slower than the activation dynamics. This also reduces the effect of the terms involving F 0a2 ðWf ð¯t ÞÞ, but this can be more clearly seen by examining Bode plots of the system given by (D.4) with g 1 ð¯t Þ; g 2 ð¯t Þ ¼ 0 using various values of F 0a2 ðWf ð¯t ÞÞ. To examine whether the choice of trajectories and arm dynamics that are used in Section 6 correspond to a model that is well approximated by (28), simulations have been conducted in which its output is compared with that of the linearizing controller applied to the full model shown in Fig. 6. The applied input, w, has been chosen to result in approximate tracking of the more rapid of the demands used. Fig. 8 shows results using the identified model parameters of one of the subjects tested. The model outputs are in close agreement which supports the use of GðsÞ to approximate the combined linearizing controller and arm system in the remainder of this paper. Care must be taken,

ARTICLE IN PRESS 374

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

Fig. 7. Stimulated arm system and linearizing controller.

2.2

Angle (rad)

1.8 1.4 1 ϑ∗f (t)

0.6

output using full model output using G(s)

0.2 0

2

4

6

8

10

Time (s) Fig. 8. Comparison of GðsÞ with full arm model.

however, to ensure that the control schemes subsequently considered are robust to the modelling error that is present in the system. Linearization along a trajectory is a well established approach which allows the linearized dynamics to be used to infer properties of the non-linear system when the state variables and input are close to those of the linearized system. In particular, if the resulting time-varying system is stable then the non-linear system also has that property in the neighbourhood of the trajectory (Shamma, 1996). This technique can be applied to assess the local stability and robustness of the proposed control system by using trajectories comprising experimental test data or those resulting from simulations using the system model. Stability of the linearized system can then be assessed using well known methods for linear time-varying systems (Kamen, 1996). This can then be repeated whilst varying the model parameters in order to gain a broader picture of the system robustness and performance properties. This approach can then be extended to provide local convergence properties for the ILC algorithms discussed in Sections 5.3 and 5.4. Since the monotonic convergence criteria given by (34) and (36) both hold for the case of time-varying plant dynamics, they can be inspected using the time-varying system resulting from linearizing the plant along a suitable trajectory. This trajectory can again be formed by simulating the system, or by direct application of experimental data. As before, model parameters can then be varied in order to gain a broader picture of ILC convergence and robustness properties. 5.2. Feedback controller The next stage is to choose a feedback controller to supply the torque demand, w, necessary for the system to track the intended references. A proportional plus derivative (PD) controller has been selected and its transfer function approximated using CðsÞ ¼

Kds þ 1

K d s þ 1

(29)

1 with 16 XX 20 . The resulting system is shown schematically in Fig. 9. The level of stimulation that first produces a response from the triceps, um , is used to supply an offset so that the feedback system operates within the torque generating capabilities of the muscle. The feedback controller is then tuned for each subject with an emphasis on robustness, since stability is of greater concern than accurate tracking, and therefore a conservative bandwidth and high gain and phase margins are desired. To illustrate the role of the end-effector dynamics on the system bandwidth, Bode plots of the closed-loop system are shown in Fig. 10 for the ratios ðba3 þ K M2 Þ=K B2 ¼ 0:03, 0:07, 0:4, 0:8 and 1:5, respectively. These have been created using an experimentally identified value of on ¼ 0:85p, and, for ease of comparison, they have all been tuned using the standard Zeigler–Nichols rules. The corresponding closed-loop systems are denoted by P a ; P b ; . . . P e, and their respective bandwidths are ob ¼ 0:53, 0.45, 0.42, 0.23 and 0.16 Hz. Alteration of ðba3 þ K M2 Þ=K B2 cannot produce bandwidths much in excess of these for the given controller and tuning method due to the limiting factor of the muscle dynamics. It is desirable to select a ratio that makes the task feel natural to the patient, but does not lead to too narrow a system bandwidth which would necessitate an excessive controller effort and correspondingly large levels of muscle torque in order to accomplish the task. The next subsections examine the effect of the end-effector dynamics on the ILC controller performance.

5.3. Phase-lead ILC Since the control task is composed of tracking the same trajectory over the course of a number of trials, between which the human arm is returned to the start position by the robot, it is an appropriate one on which to apply ILC. The phase-lead ILC algorithm has been shown to provide excellent results despite its simplicity and limited parameter set (see, for example, Freeman, Lewin, & Rogers, 2007; Wang & Longman, 1996, for further details). It is given in discrete time by vkþ1 ðzÞ ¼ vk ðzÞ þ Lzl ek ðzÞ

(30)

where L is a scalar gain, l is the phase-lead in samples, k is the trial number, and vk and ek are the input and error signals for the kth trial. The controller structure is shown schematically in Fig. 11,  in which the error is calculated as eðtÞ ¼ Wf ðtÞ  Wf ðtÞ. ILC is commonly used in conjunction with a feedback controller, in the manner shown, in order to increase closed-loop stability and noise rejection. From (30) it is then possible to obtain the relationship ekþ1 ðzÞ ¼ ek ðzÞ  Lzl PðzÞek ðzÞ

(31)

with PðzÞ ¼

CðzÞGðzÞ 1  CðzÞGðzÞ

(32)

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

375

Magnitude (dB)

Fig. 9. Block diagram of feedback control system.

0 −50 −100 −150 10−1

100 Frequency (rads)

101

100

101

Phase (deg)

0 −100

ωb = 0.53 Hz ωb = 0.45 Hz ωb = 0.42 Hz ωb = 0.234 Hz ωb = 0.1603 Hz

−200 −300 10−1

Fig. 10. Bode plots of linearized feedback systems for different end-effector dynamics.

where GðzÞ and CðzÞ are the discretized representations of (28) and (29), respectively. The relationship ekþ1 ðzÞ ¼ ek ðzÞð1  LPðzÞzl Þ

(33)

yields the monotonic convergence criterion joT s

j1  LPðe

jolT s

Þe

jo1

order to produce an extremely robust system at the expense of convergence speed. Reducing L can be shown to increase robustness to model uncertainty over all frequencies, so this is a sensible approach to take. It can also be shown that choosing the phase-lead to maximize the convergence at a given frequency also achieves a maximum robustness to gain and phase uncertainty at that frequency (Ha¨to¨nen et al., 2004), which means that there is no compromise between robustness and convergence speed with respect to choosing this parameter. Therefore the phase-lead will be selected to maximize the convergence over a suitable frequency range, and a zero-phase filter implemented to ensure stability at all higher frequencies (this being advisable even if the convergence criterion is satisfied). For general application, it is natural for this frequency range to correspond to the bandwidth of the system. The system with ob ¼ 0:45 Hz can be seen to have maximum convergence over its bandwidth using a phase-lead of 30 samples (although a phase-lead of 20 samples produces greater stability over higher frequencies, which is an important factor if a filter is not used). Similarly the system with ob ¼ 0:53 Hz has a maximum convergence over its bandwidth of 20 samples. When using these phase-leads the system with the higher bandwidth has the property of monotonic convergence at greater frequencies than the system with lower bandwidth. For further discussion and guidelines on the choice of phase-lead and filter design, see Freeman, Lewin, and Rogers (2005). 5.4. Adjoint ILC

(34)

for o up to the Nyquist frequency. Satisfying this at a given frequency is a sufficient condition for monotonic convergence at that frequency (this is a well established convergence criterion, see for example, Elci, Longman, Phan, Juang, & Ugoletti, 1994). Furthermore, the convergence speed is dictated by the magnitude of the left-hand side; if it is close to zero, convergence will occur in a single trial, while if it is greater than one, divergence is likely to occur at that frequency. Therefore L and l are chosen so that the left-hand side is minimized to provide the greatest convergence over those frequencies present in the demand. Frequencies above this, or those at which plant uncertainty may cause the criterion to be violated, are removed through the use of a non-causal zerophase filter applied to the error. A more in-depth discussion relating to the choice of parameters is given in Freeman et al. 1 (2007). The sample time, T s ¼ 40 , will be used since this corresponds to the frequency at which stimulation pulses are applied to the patient. This frequency is synchronized with the robotic control system and each pulse is produced with a delay of less than 10 ms. The effect of varying the phase-lead is shown for two of the systems that were considered in the previous section. In Fig. 12 the monotonic convergence criterion (34) is shown for the system with bandwidth ob ¼ 0:45 Hz using a variety of phase-leads, l. In Fig. 13 the criterion is shown for the system with bandwidth ob ¼ 0:53 Hz. A value of L ¼ 0:2 has been chosen in both cases in

The adjoint ILC algorithm is given in discrete form by vkþ1 ðzÞ ¼ vk ðzÞ þ bP ðzÞek ðzÞ

(35)

where P ðzÞ is the adjoint of the plant model used. An attractive feature of the method is that, with a sufficiently small scalar multiplier, b, it is guaranteed to satisfy the condition for monotonic convergence over all frequencies, and hence ensure a satisfactory transient response (Chen & Longman, 2002). The monotonic convergence criterion in this case is j1  bPðejoT s ÞP  ðejoT s Þjo1

(36)

which leads to 0objPðejoT s Þj2 o2

(37)

for o up to the Nyquist frequency. The algorithm has been found to provide a high level of robustness to plant uncertainty (Freeman et al., 2005). To examine the effect of the choice of gain, plots of the monotonic convergence criterion (36) can be drawn for the two systems examined in the last subsection for a variety of b. In this way it can be shown that both systems satisfy the criterion for all values of b between 0 and 1. For frequencies where the closedloop system has a gain close to unity, approximate convergence will occur in a single trial for values of b close to one. In practice, however, b is chosen to be significantly lower in order to increase the system robustness at the expense of the convergence rate.

ARTICLE IN PRESS 376

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

Fig. 11. Block diagram of the ILC system.

|1 − L P(e jωTs) e jω λ|

1.1 λ=0 λ = 10 λ = 20 λ = 30 λ = 40 λ = 50

1 0.9 0.8 0.7 10−1

100 101 Frequency (rad/s)

102

Fig. 12. Monotonic convergence criterion for the system with ob ¼ 0:45 Hz for L ¼ 0:2 and various l.

|1 − LP(e jωTs) e jωλ|

1.1 λ=0 λ = 10 λ = 20

1

λ = 30 λ = 40 λ = 50

0.9 0.8 0.7 10−1

100 101 Frequency (rad/s)

102

Fig. 13. Monotonic convergence criterion for the system with ob ¼ 0:53 Hz for L ¼ 0:2 and various l.

As with the phase-lead algorithm, the convergence rate that can be expected at a given frequency is closely connected to the system bandwidth. The convergence rate can be increased within a given range by convolving P  ðzÞ with a zero-phase band-pass filter and substituting in (35) (see Freeman et al., 2007, for details). The same effect could be achieved by replacing P ðzÞ with the plant inverse P 1 ðzÞ, and then using a low-pass filter. However, the increase in convergence speed is gained at the price of reduced robustness, and so a zero-phase band-pass filter will not be used unless the convergence rate is found to be unacceptably low for the desired application.

6. Experimental results The arm model of each of the 18 unimpaired subjects tested was first identified using tests and procedures described in

Freeman et al. (2008a). These involved taking measurements of the arm, establishing its maximum range of movement, and then fitting a circle to the trajectory traced out by the elbow, in order to provide g and the position of the shoulder joint. The arm was then held stationary while FES was applied to the triceps using a ramp signal, in order to produce the functions, hIRC ðuÞ and hLAD ðtÞ, in the muscle model, using deconvolution and an non-linear optimization procedure. Stimulation sequences and kinematic trajectories, imposed on the arm by the robot, were then applied, and an LMS optimization was subsequently used to arrive at the remaining model parameters. Fig. 14 shows an example of the shape of the reference used to produce the results in this paper in relation to the position of the subject’s shoulder joint. The reference was set at an angle of 20 from the y-axis but was individually calculated for each subject so that it extended their arm from 55% to 95% of their total arm length over the course of the movement. Two trajectories have been created by moving along this reference at two different speeds. Each trajectory started with a waiting period when it was set equal to the starting point of the reference. The ‘slow’ trajectory lasted for 12.5 s in total (a 5 s waiting period and a 7.5 s movement along the reference), and the ‘fast’ trajectory lasted for 10 s (a 5 s waiting period followed by a 5 s movement). The waiting period was included to allow the ILC update to begin before the arm was required to move. Before each trial began, the subject’s arm was moved to the initial position by the robot and then released only when the trajectory started. The subjects were not shown the trajectory before or during the test. The values of K B2 and K M2 that dictate the end-effector characteristic were set at 5:78 N m=rad s1 and 0:29 N m=rad s2 , respectively, since these created a natural feel to the system and allowed the chosen trajectories to be accomplished with moderate effort whilst not limiting the bandwidth excessively. For a typical identified value of ba3 these gains produce ðba3 þ K M2 Þ=K B2  0:1, which leads to a closed-loop system bandwidth of approximately 0:45 Hz, when calculated using the same muscle dynamics and controller tuning procedure as adopted in Section 5.1. The values of on ¼ 0:85p and ba3 ¼ 0:271 were identified for the subject whose experimental results are presented in Sections 6.2 and 6.3. 6.1. Linearizing PD controller The feedback controller gains were tuned for each subject in the manner described in Section 5.1. Two or three repetitions of the slow trajectory were used to fine-tune the gains given by the standard Zeigler–Nichols rules, and these values were used for the duration of the tests. In the remainder of this section ILC results are presented for one of those subjects tested, the PD controller gains used in this case were K p ¼ 10 and K d ¼ 2. The effect of the linearizing feedback controller in the absence of ILC can be seen by inspection of the first trial of the ILC results given here.

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

−0.6

λ = 15 λ = 20 λ = 25 λ = 30 λ = 35 λ = 40

RMS Error (m)

0.03

−0.7 y position (m)

377

Reference Experimental data Shoulder joint position Initial arm position Final arm position Elbow position

−0.8 −0.9

0.02

0.01

0

−1

1

2

3

4

−1.1

5 6 Trial No.

7

8

9

10

Fig. 15. Single subject phase-lead ILC results for the slow trajectory using L ¼ 0:2 and various l.

−1.2 0.2

0.3

0.4 0.5 x position (m)

0.6

0.7

0.8

Fig. 14. Reference trajectory and position of subject tested.

6.2. Phase-lead ILC Fig. 15 shows results using the phase-lead algorithm for various values of l. The slow trajectory is used together with a value of L ¼ 0:2. The learning gain has been chosen conservatively so that the speed of convergence has been sacrificed for greater robustness. The zero-phase filter mentioned in Section 5.3 was not implemented since no signs of instability were observed over the course of any of the experimental tests conducted. It can be seen that the best performing phase-lead appears to be 35 samples, which is similar to that found by inspection of the monotonic convergence criterion for P b, a system with similar bandwidth, discussed in Section 5 and shown in Fig. 12. The root mean square (RMS) error corresponding to this phase-lead converges to approximately 5 mm and has a minimum value of 3.2 mm. Fig. 16 shows the results obtained using L ¼ 0:2 and l ¼ 35 in greater detail. It can be seen that during the first five trials the error reduces monotonically and the reference is tracked extremely well. This reflects the improvement in tracking accuracy that simple ILC schemes can provide compared with the use of feedback controllers alone. In order to examine the results further, Fig. 17 shows the associated stimulation input and ILC update. It can be seen from (a) that the stimulation applied during k ¼ 1 saturates at 300 ms, but that the effect of further trials removes this saturation and produces lower levels of stimulation over the course of the trial. ILC also leads to the application of stimulation during the initial 5 s waiting period before movement is required. Fig. 17(b) shows that the updated reference, vk , gets progressively closer to a fixed trajectory over repeated trials, indicating convergence.

6.3. Adjoint ILC Fig. 18 shows results using the adjoint algorithm. The slow trajectory is again used and b ¼ 0.2. The plant model is labelled P, and to examine the robustness of the algorithm, the models P a , Pc , P d and P e , described in Section 5.1, have also been used in its place. The model Pb has not been used due to its close similarity with P. It is clear that the algorithm is capable of exhibiting robustness to significant model uncertainty and produces minimum error and convergence rates comparable to the performance of the phaselead algorithm. Use of the actual plant model results in convergence to approximately 5 mm in four trials and this level

Angle (rad)

0.1

ϑu∗(t) k=1 k=2 k=3 k=4 k=5

1

0.5

0 0

2

4

6 Time (s)

8

10

12

4

6 Time s

8

10

12

2 Angle (rad)

0

1.5

ϑu∗(t) k=1 k=2 k=3 k=4 k=5

1 0.5 0 0

2





Fig. 16. Single subject experimental tracking of (a) Wu ðtÞ and (b) Wf ðtÞ for the slow trajectory using phase-lead ILC with L ¼ 0:2 and l ¼ 35.

of error is maintained over the remaining trials. Tracking results using P are very similar to those obtained using the phase-lead update, both in terms of the monotonicity of convergence and the level of stimulation applied. 6.4. Experimental comparison The performance of the phase-lead and adjoint algorithms has been compared against norm-optimal ILC (NOILC), an approach which has received considerable attention in the ILC literature due to its mature theoretic basis (Amann, Owens, & Rogers, 1996). The implementation of this technique requires extensive computation, but it is a natural choice to assess the effectiveness of the proposed methods since it has been found to outperform a variety of other approaches in experimental tests (Ratcliffe, 2005). In NOILC the input on the ðk þ 1Þth trial is chosen to minimize Jkþ1 ðvkþ1 Þ ¼ kekþ1 k2Y þ kvkþ1  vk k2U

(38)

where U and Y denote the input and output function spaces, respectively. The appearance of vkþ1  vk in (38) allows control

ARTICLE IN PRESS 378

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

0.03

200

RMS Error (m)

Stimulation, u (μs)

300

k=1 k=5 k = 10

100 0 0

2

4

6 Time (s)

8

10

Phase−lead Adjoint NOILC, Q = 0.5 NOILC, Q = 1

0.02

NOILC, Q = 2

0.01

12 0 1

2

3

4

5 6 Trial No.

7

8

9

10

1

2

3

4

5 6 Trial No.

7

8

9

10

0

0.04 k=1 k=2 k=3 k=4 k=5

−0.1 −0.2 0

2

4

6 Time (s)

8

10

12

0.03 0.02 0.01

Fig. 17. Single subject experimental (a) stimulation and (b) updated input results for the slow trajectory using phase-lead ILC with L ¼ 0:2 and l ¼ 35.

∗

0

Fig. 19. Single subject comparison of phase-lead and adjoint algorithms with NOILC for (a) the slow trajectory, and (b) the fast trajectory.

Pa

0.03 RMS Error (m)

RMS Error (m)

ILC input, vk

0.1

P Pc∗

Pd∗

0.02

6.5. Results from multiple subjects

Pe∗

0.01

0 1

2

3

4

5 6 Trial No.

7

8

9

10

Fig. 18. Single subject adjoint ILC results using for the slow trajectory with b ¼ 0:2 and various plant models.

over the possibility of excessive control input signals. The noncausal solution vkþ1 ðzÞ ¼ vk ðzÞ þ P  ðzÞekþ1 ðzÞ

Given the similarity in performance of the phase-lead and adjoint ILC algorithms, group results are presented for the adjoint controller implementation only. Fig. 20 shows results from the first four subjects tested. These indicate that convergence can be achieved within five trials for both trajectories, and that an RMS tracking error of less than 10 mm over the course of the movement (and in some instances less than 5 mm) is possible in all cases. Table 1 shows the mean and standard deviation of the RMS error obtained during the last trial, for all subjects tested. Trajectory 1 refers to the reference previously described, whilst trajectory 2 is a reference whose inclination from the y-axis is increased to 40 . The results confirm that the control methods that have been implemented are capable of producing high levels of tracking accuracy.

(39) 7. Discussion

is transformed into a causal implementation incorporating both feedback and feedforward actions. To apply NOILC (38) is expressed as a linear quadratic cost in ekþ1 and vkþ1  vk using weighting values of Q and R, respectively (see Amann et al., 1996, for information concerning the practical application of NOILC). A Kalman estimator has been used to provide the states of the linear plant approximation (D.4) that are required in the implementation, derived using output and state covariance weights of 1 and 10, respectively. Fig. 19 shows representative results using a variety of Q values (it is the ratio of the weights that influences the optimal solution so that R has been set at unity) compared with the use of the phase-lead and adjoint algorithms. It can be seen that the performance of NOILC does not exceed that of the simpler ILC methods, despite its far more complex structure.

The experimental results show that superior tracking performance has been achieved compared with alternative control methods that have been applied to the upper limb. These include the previous application of ILC (Dou et al., 1999), the use of openloop controllers (Davoodi & Andrews, 2004; Popovic & Popovic, 1998), and the limited number of closed-loop controllers that have been experimentally applied in this area (Lan et al., 1994; Tresadern et al., 2006). The major advantage of the proposed control method, however, is the simplicity of tuning, and the absence of any training experiments. The simple ILC laws that have been examined also permit the degree of assistance to easily be changed via the use of a forgetting factor (see Freeman et al., 2007) which can also be used to promote voluntary effort when

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

379

8. Conclusions and further work

RMS Error (m)

0.03

0.02

0.01

0

RMS Error (m)

1

2

3

4

5 6 Trial No.

7

8

9

10

A study using 18 unimpaired subjects has been conducted to investigate the feasibility of tracking reaching trajectories using FES applied to the triceps muscle in combination with robotic assistance. Results have shown that a high level of tracking performance was achieved by all subjects. The two ILC schemes implemented produce similar results and significantly improve upon the performance achieved when solely using a linearizing PD controller. The results confirm the efficacy of a combined robotic assistance and stimulation system that has been designed for the purpose of stroke rehabilitation. Tests have recently begun using the system with five stroke patients in treatment sessions comprising up to 25 visits of 1 h duration.

0.04

Acknowledgements

0.03

This work is supported by the Engineering and Physical Sciences Research Council (EPSRC). Grant no. EP/C51873X/1.

0.02 Appendix A. Biomechanical arm model

0.01 0 1

2

3

4

5 6 Trial No.

7

8

9

10

Fig. 20. Adjoint ILC results for (a) the slow trajectory, and (b) the fast trajectory, using b ¼ 0:2 for four subjects.

Table 1 Mean of last trial RMS error for all 18 subjects (standard deviation in brackets)

The additional terms appearing in the biomechanical model of the human arm given by (7) are 2 3 " # 2ca1 W_ f ca1 W_ f ba1 ba2 5, Ba ðqa Þ ¼ ; C a ðqa ; q_ a Þ ¼ 4 ba2 ba3 ca1 W_ u ca2 W_ f J Ta ðqa Þ "

ðlu1 þ lu2 Þcg su  ðlf 1 þ lf 2 Þsuf

ðlu1 þ lu2 Þcg cu þ ðlf 1 þ lf 2 Þcuf

ðlf 1 þ lf 2 Þsuf

ðlf 1 þ lf 2 Þcuf

¼

#

(A.1) Type

Slow Fast

Mean of last trial RMS error (mm)

with

Trajectory 1

Trajectory 2

9.41 (5.67) 12.18 (6.94)

10.22 (6.07) 12.93 (6.87)

2

ba1 ¼ mu ðlu1 cg Þ2 þ Iu þ mf ðlf 1 þ ððlu1 þ lu2 Þcg Þ2 þ 2ðlu1 þ lu2 Þcg lf 1 cf Þ þ If 2

ba2 ¼ mf ðlf 1 þ ðlu1 þ lu2 Þcg lf 1 cf Þ þ If applied to stroke patients. Another advantage is that these laws are easy to tune which is an important concern when conducting tests with patients. Although ILC algorithms have moved beyond these relatively simple structure types and now encompass a wide range of plant models and control law structures, the approach taken here was to apply ILC laws with the simplest structure which could meet the necessary performance requirements. This was reflected in Section 6.4 in which results were compared with those using a far more complex ILC scheme and similar levels of performance were observed. Future extensions to the experimental set-up will involve removal of the constraints to the patient’s arm to allow it to be lifted against gravity. To guarantee the subject’s safety, however, it is necessary to first stabilize their arm using a sling or passive mechanical support that gently returns their arm to a neutral position in the task workspace. Alternatively, this could be substituted by a FES feedback controller applied to the anti-gravity muscles of the arm. If more complex three-dimensional tasks are used and it is not feasible to clearly display a fixed trajectory, it may be necessary to replace ILC by controllers which use muscle EMG signals to control the applied FES in accordance with the subject’s voluntary intention.

2

ba3 ¼ mf lf 1 þ If þ Ie

sg 1  c2f c2g

!2

ca1 ¼ mf ðlu1 þ lu2 Þcg lf 1 sf !

s2g c2g cf sf

ca2 ¼ 2Ie

(A.2)

ð1  c2f c2g Þ3

Appendix B. Robotic assistance controller gains Suitable choices for the robotic controller matrices used in the control scheme presented in Section 4.1 satisfy 8 1 < K Bx ¼ J T a ðqa ÞK Bq J a ðqa Þ (B.1) _ : K M ðJ ðq Þq€ þ J ðq ; q_ Þq_ Þ ¼ J T ðq ÞK M x

a

a

a

a

a

a

a

a

a

q

with the gains given by (13). Using (22) the remaining controller gain (19) can then be written explicitly as  K K 1 CðW^ f Þ  Wu I (B.2) K K x ðxa ; CðÞÞ ¼ lðlu1 þ lu2 Þcg sf

ARTICLE IN PRESS 380

C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

Appendix C. Approximate robot assisted arm model

arm dynamics be represented by the state-space system ½Ap ; Bp ; C p  so that the relationship between y and Wf is given by

The robot assisted arm model shown in Fig. 5 can be simplified by substituting equations describing the action of the robotic assistance into the human arm model. The bottom row of (10), with the expressions given by (16) and (17), the control gains (13), and using the general form of friction function, is 2 3 2 3 W_ u W€ u ½ba2 ba3 4 € 5 þ ½ca1 W_ u ca2 W_ u 4 _ 5 þ F a2 ðWf ; W_ f Þ 0

Wf

1

Wf

B sf cg C ¼ T b @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA  K B2 W_ f  K M2 W€ f 1  c2f c2g

(C.1)

þ

þ ca2 W_ u W_ f þ F a2 ðWf ; W_ f ÞÞ

(D.4)

(D.5)

in which g 2 ðtÞ ¼ g ðWf ðtÞÞ=gðWf ðtÞÞ. (C.2)

CðWf Þ, the system is then independent of Wu , and can be written solely in terms of Wf , giving 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1(    1  c2f c2g A K B W_ f þ K M þ ba3 þ ba2 d CðWf Þ W€ f Tb ¼ @ 2 2 dWf sf c g !   2 2 2 d d d CðWf Þ þ ca2 CðWf Þ W_ f þ ba2 2 CðWf Þ þ ca1 d Wf dWf dWf ) þF ðW ; W_ Þ (C.3) f

In practice ca1 , ca2 5ba2 ; ba3 (see Freeman et al., 2008a, for further details) and it will be assumed that the trajectories are chosen so that there are no sudden changes in Wu , as Wf is varied, so that the derivatives of CðÞ are small. The latter assumption has already been used when deriving the robotic assistance control scheme, and is easy to satisfy in practice. The terms involving Ie are also small, a typical value being 0:005 kg m2 (see Freeman et al., 2008a, for details). Therefore (C.3) can be approximated as 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1  1  c2f c2g AfK B W_ f þ ðK M þ ba3 ÞW€ f þ F a2 ðWf ; W_ f Þg Tb  @ (C.4) 2 2 sf c g

Appendix D. Linearized arm system and linearizing controller In order to linearize the combined arm system and linearizing controller shown in Fig. 7, the linear activation dynamics can be represented by the state-space system ½Am ; Bm ; C m , and the relationship between it w and y is then given by   zðtÞ x_^ m ðtÞ ¼ Am x^ m ðtÞ þ Bm gðWf ðtÞÞ yðtÞ ¼ C m x^ m ðtÞgðWf ðtÞÞ (D.1) where z ¼ w þ F a2 ðWf Þ. The transformation xm ðtÞ ¼ x^ m ðtÞgðWf Þ allows the system to absorb the effect of gðWf Þ to become x_ m ðtÞ ¼ ðAm þ Ig 1 ðtÞÞxm ðtÞ þ Bm zðtÞ yðtÞ ¼ C m xm ðtÞ

then linearizing the complete system at t ¼ t¯ yields 2 3 " # Am  IC p Ap xp ðt¯ Þg 1 ðt¯ Þ W x_~ m ðtÞ 4 5¼ Bp C m Ap  Bp F 0a2 ðWf ðt¯ ÞÞC p x_~ p ðtÞ # " # " " # x~ m ðtÞ x~ m ðtÞ 0 ~ ~  w; Wf ¼ ½0 C m  þ x~ p ðtÞ x~ p ðtÞ Bm

00

2 2 2 case W_ u ¼ W_ f ðd=dWf ÞCðWf Þ and W€ u ¼ W€ f ðd=dWf ÞCðWf Þ þ W_ f ðd =dWf Þ

f

(D.3)

W ¼ IC p Ap ðxm ðt¯ Þg 1 ðt¯ Þ þ xp ðt¯ ÞC p xm ðt¯ Þðg 2 ðt¯ Þ  g 1 ðt¯ Þ2 ÞÞ þ Bm F 0a2 ðWf ð¯t ÞÞC p

It will now be assumed that the robotic assistance system provides perfect tracking of W^ u by Wu, so that Wu ¼ CðWf Þ. In this

a2

Wf ðtÞ ¼ C p xp ðtÞ

where

so that 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1  1  c2f c2g AðK B W_ f þ ðK M þ ba3 ÞW€ f þ ba2 W€ u Tb ¼ @ 2 2 sf c g 2 ca1 W_ u

x_ p ðtÞ ¼ Ap xp ðtÞ þ Bp ðyðtÞ  F a2 ðWf ÞÞ

(D.2)

_ Wf ðtÞÞ=gðWf ðtÞÞ. The effect of gðÞ is therefore to in which g 1 ðtÞ ¼ gð directly shift the eigenvalues of its state transition matrix. Let the

References Ahn, H. S., Chen, Y., & Moore, K. L. (2007). Iterative learning control: Brief survey and categorization. IEEE Transactions on Systems, Man, and Cybernetics Part C: Applications and Reviews, 37(6), 1099–1121. Amann, N., Owens, D. H., & Rogers, E. (1996). Iterative learning control using optimal feedback and feedforward actions. International Journal of Control, 65, 277–293. Baratta, R., & Solomonow, M. (1990). The dynamic response model of nine different skeletal muscles. IEEE Transactions on Biomedical Engineering, 37(3), 243–251. Burridge, J. H., & Ladouceur, M. (2001). Clinical and therapeutic applications of neuromuscular stimulation: A review of current use and speculation into future developments. Neuromodulation, 4(4), 147–154. Castro-Alamancos, M., Garcia-Segura, L., & Borrell, J. (1992). Transfer of function to a specific area of the cortex after induced recovery from brain damage. European Journal of Neuroscience, 4(9), 853–863. Chen, K., & Longman, R. W. (2002). Stability issues using FIR filtering in repetitive control. Advances in the Astronautical Sciences, 206, 1321–1336. Colgate, J. E., & Hogan, N. (1988). Robust control of dynamically interacting systems. International Journal of Control, 48(1), 65–88. Davoodi, R., & Andrews, B. J. (1988). Computer simulation of FES standing up in paraplegia: A self-adaptive fuzzy controller with reinforcement learning. IEEE Transactions on Rehabilitation Engineering, 6(2), 151–161. Davoodi, R., & Andrews, B. J. (2004). Fuzzy logic control of FES rowing exercise in paraplegia. IEEE Transactions on Biomedical Engineering, 51(3), 541–543. Dou, H., Tan, K. K., Lee, T. H., & Zhou, Z. (1999). Iterative learning control of human limbs via functional electrical stimulation. Control Engineering Practice, 7, 315–325. Elci, H., Longman, R. W., Phan, M., Juang, J. N., & Ugoletti, R. (1994). Discrete frequency based learning control for precision motion control. In Proceedings of the 1994 IEEE international conference on systems, man, and cybernetics (Vol. 3, pp. 2767–2773), San Antonio, TX. Freeman, C. T., Hughes, A. M., Burridge, J. H., Chappell, P. H., Lewin, P. L., & Rogers, E. (2008a). A model of the upper extremity using surface FES for stroke rehabilitation. ASME Journal of Biomechanical Engineering, in press. Freeman, C. T., Hughes, A. M., Burridge, J. H., Chappell, P. H., Lewin, P. L., & Rogers, E. (2008b). A robotic workstation for stroke rehabilitation of the upper extremity using FES. Medical Engineering and Physics, in press, doi:10.1016/j.medengphy.2008.05.008. Freeman, C. T., Lewin, P. L., & Rogers, E. (2005). Experimental evaluation of iterative learning control algorithms for non-minimum phase plants. International Journal of Control, 78(11), 826–846. Freeman, C. T., Lewin, P. L., & Rogers, E. (2007). Further results on the experimental evaluation of iterative learning control algorithms for non-minimum phase plants. International Journal of Control, 80(4), 569–582. Hatwell, M. S., Oderkerk, B. J., Sacher, C. A, & Inbar, G. F (1991). The development of a model reference adaptive controller to control the knee joint of paraplegics. IEEE Transactions on Automatic Control, 36(6), 683–691. Ha¨to¨nen, J. J., Harte, T., Owens, D. H., Ratcliffe, J. D., Lewin, P. L., & Rogers, E. (2004). Iterative learning control—What is it all about? In Proceedings of IFAC workshop on adaption and learning in control and signal processing and the IFAC workshop on periodic control systems (pp. 547–553), Yokohama, Japan. Hudson, J., Scharaschkin, A., Wilkins, M., & Taylor, K. (2005). Reducing brain damage: Faster access to better stroke care. London: National Audit Office. Hunt, K. J, Jaime, R. P., & Gollee, H. (2001). Robust control of electrically-stimulated muscle using polynomial H1 design. Control Engineering Practice, 9, 313–328. Hunt, K. J., Munih, M., & Donaldson, N. de N. (1997). Feedback control of unsupported standing in paraplegia—part I: Optimal control approach. IEEE Transactions on Rehabilitation Engineering, 5(4), 331–340.

ARTICLE IN PRESS C.T. Freeman et al. / Control Engineering Practice 17 (2009) 368–381

Kamen, E. W. (1996). Fundamentals of linear time-varying systems. CRC controls handbook (pp. 451–468). Boca Raton, FL: CRC Press. Krebs, H. I., Hogan, N., Aisen, M. L., & Volpe, B. T. (1998). Robot-aided neurorehabilitation. IEEE Transactions on Rehabilitation Engineering, 6, 75–87. Lan, N., Feng, H. Q., & Crago, P. E. (1994). Neural network generation of muscle stimulation patterns for control of arm movements. IEEE Transactions on Rehabilitation Engineering, 2(4), 213–224. Lum, P. S., Burgar, C. G., & Shor, P. C. (2004). Evidence for improved muscle activation patterns after retraining of reaching movements with the MIME robotic system in subjects with post-stroke hemiparesis. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 12(2), 186–194. Mant, J., Wade, D., & Winner, S. (2004). Health care needs assessment: Stroke. In A. Stevens, J. Raftery, J. Mant, & S. Simpson (Eds.), Health care needs assessment: The epidemiologically based needs assessment reviews. Oxford: Radcliffe Medical Press. Moore, K. L. (1992). Iterative learning control for deterministic systems. Berlin: Springer. Popovic, D., & Popovic, M. (1998). Tuning of a nonanalytical hierarchical control system for reaching with FES. IEEE Transactions on Biomedical Engineering, 45(2), 203–212. 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. Ratcliffe, J. (2005). Iterative learning control of multi-axis systems. Ph.D. thesis, School of Electronics and Computer Science, University of Southampton, UK.

381

Royal College of Physicians (1989). Stroke: Towards better management. London: Royal College of Physicians. Rushton, D. N. (2003). Functional electrical stimulation and rehabilitation—an hypothesis. Medical Engineering and Physics, 25(1), 75–78. Schutter, J. D., & Brussel, H. V. (1988). Compliant robot motion ii. A control approach based on external control loops. International Journal of Robotics Research, 7(4), 18–33. Shamma, J. S. (1996). Linearization and gain-scheduling. CRC controls handbook (pp. 388–397). Boca Raton, FL: CRC Press. Shue, G., Crago, P. E., & Chizeck, H. J. (1995). Muscle-joint models incorporating activation dynamics, moment-angle, and moment-velocity properties. IEEE Transactions on Biomedical Engineering, 42(2), 213–223. Siciliano, B., & Villani, L. (1999). Robot force control. Boston, MA: Kluwer Academic Publishers. Thorsen, R., Spadone, R., & Ferrarin, M. (2001). A pilot study of myoelectrically controlled FES of upper extremity. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 9(2), 161–168. Tresadern, P., Thies, S., Kenney, L., Howard, D., & Goulermas, J. Y. (2006). Artificial neural network prediction using accelerometers to control upper limb FES during reaching and grasping following stroke. In Proceedings of the 28th annual international conference of the IEEE engineering in medicine and biology society (pp. 2916–2919), New York, USA. Wang, Y., & Longman, R. W. (1996). Use of non-causal digital signal processing in learning and repetitive control. Advances in the Astronautical Sciences, 90, 649–668.