Shared Control of a Powered Exoskeleton and Functional Electrical Stimulation Using Iterative Learning

A hybrid exoskeleton comprising a powered exoskeleton and functional electrical stimulation (FES) is a promising technology for restoration of standing and walking functions after a neurological injury. Its shared control remains challenging due to the need to optimally distribute joint torques among FES and the powered exoskeleton while compensating for the FES-induced muscle fatigue and ensuring performance despite highly nonlinear and uncertain skeletal muscle behavior. This study develops a bi-level hierarchical control design for shared control of a powered exoskeleton and FES to overcome these challenges. A higher-level neural network–based iterative learning controller (NNILC) is derived to generate torques needed to drive the hybrid system. Then, a low-level model predictive control (MPC)-based allocation strategy optimally distributes the torque contributions between FES and the exoskeleton’s knee motors based on the muscle fatigue and recovery characteristics of a participant’s quadriceps muscles. A Lyapunov-like stability analysis proves global asymptotic tracking of state-dependent desired joint trajectories. The experimental results on four non-disabled participants validate the effectiveness of the proposed NNILC-MPC framework. The root mean square error (RMSE) of the knee joint and the hip joint was reduced by 71.96 and 74.57%, respectively, in the fourth iteration compared to the RMSE in the 1st sit-to-stand iteration.

be driven to zero. Thus, the output, y, is defined as 13 y θ−h.
where L 2 f y is the 2 nd -order Lie derivative of y, L g M L f y and L g F L f y are the decoupling matrices, and 15 d ∈ R N is the disturbance of the system output. 16 We consider thatȳ 2 =ẏ (i) , where superscript (i) = 1, 2, ..., N shows i th row of a vector.

17
As the ultimate motivation is to develop a novel Iterative learning control (ILC) approach, (23) is expressed 18 for a k th iteration as where σf 1,k + f 2,k (i) is equal to i th element of L 2 f y and is expressed as a sum of structured and 20 unstructured uncertain nonlinear terms. Specifically, (σf 1,k ) (i) is the linearly parameterizable (structured is the motor control constant gain for o th electric motor input. u M,o and is assumed to be known. B F,j,k is 25 an unknown control gain associated with j th normalized FES input, u F,j,k . Normalized FES input, u F,j,k is 26 bounded as follow N F is the total number of electrical stimulator and N m is the total number of motors. For simplicity 28 purposes, the superscript i (i = 1, 2, ..., N ) will be dropped hereafter in all the corresponding notations.

29
We assume that the terms f 2,k and B F,k can be represented using two ideal NNs as follows where X k ∈ R 2N +1 is the augmented input vector for the aforementioned two NNs and is defined as and P ∈ R 2N +1×N in are the ideal weight matrices for term f 2,k 33 and Q ∈ R N Ω ×N is the ideal weight vector for B F,k . The input layer is made of 2N + 1 neurons. N 34 is the output layer neurons number, and the hidden layer numbers of neurons in the two NNs are N in 35 and N Ω . Λ k : R N in → R N 2 +1 is the first NN activation function in (26) that maps the input layer to 36 the hidden layer, where P ∈ R (2N +1)×N in is the weight matrix corresponding to the augmented input. 37 φ j,k : R 2N +1 → R N Ω is the activation function in (27) that maps the input layer to the output layer.

38
ε 1,k ∈ R N and ε 2,k ∈ R N are the unknown functional reconstruction errors for the two NNs, which all elements of them are bounded with ε 1,j,k ≤ε 1 and ε 2,j,k ≤ε 2 , respectively, whereε 1 ,ε 2 ∈ R + (Chen et al. (2010)). By the NN universal approximation property the bounds on the functional reconstruction 41 errors can be made smaller by choosing a larger number of neurons in the layers (Lewis et al. (2002)).

42
Based on the subsequent stability analysis, the estimates of ideal weight matrices in the NNS,Ŵ k j ,P k j , 43 andQ j,k are updated using the following gradient method where ρ 1 ∈ R + and ρ 2 ∈ R + are user defined constants, j represents j th element of a matrix and E f 2,k is 47 defined as

49
Using (17), by adding and subtracting N f j=1 ψ j,k u F,j,k to (24) results in where Q j,k = Q j −Q j,k is the weight estimation error, and β ε,j = ε 2,j,k − ( j ) and it is bounded by 51β ε,j ∈ R + . (32) describes the closed loop dynamic of the output. It shows how the stimulation amplitude 52 level and motor torques are contributing to the output dynamics. Additionally, it shows how the estimation 53 errors for control gain matrix and the system dynamics can have influence on the behavior of the output of 54 the closed loop system. Therefore, by substituting (7) in (4) and considering (11) , the following equation 55 is resulted based on (32)

Frontiers
Using (33) and (3), the dynamics of the sliding surface can be obtained that Based on (34), the sliding surface's dynamic equation can be rewritten as ..,N, the subsystem in (24) can reach the equilibrium point (s i , e 0,i ) = (0, 0) 60 asymptotically, if the control inputs are selected as (4) and (16).

61
The following Lyapunov-like functional is considered tr Q T j,kQj,k . t 0 is the start time of iterations, t is the elapsed time after the start of an iteration and 64 γ, q c , ξ ∈ R + are constants. It is considered that at the beginning of each iteration, the exoskeleton 65 wearer has same sitting initial position. Subsequent mathematical procedure proves that the tracking error 66 and estimation error converge to zero based on the energy difference between two successive iterations.

67
Accordingly, for V (1) k , the energy difference between each two iterations can be written as 68

∆V
(1) Then the following formula can be obtained The following equation is derived by substituting (8) to equation (38)

∆V
(1) This is a provisional file, not the final typeset article

74
The third energy function difference between two successive iterations can be obtained as Based on the theorem in Chen et al. (2012), the following equation can be derived By substituting (5) to (43) and substituting the results to (42), the following expression can be obtained The fourth energy function difference between two iterations is given as Frontiers By considering the update laws in (28), (29), and (31), the fourth energy function difference will be given 80 as follows The fifth energy function difference between two iterations can be written as

83
By considering Q j,k (t 0 ) = 0, ∆V k can be written as

87
To show the convergence of the output tracking error and parameters estimation error, based on the 88 fact that s k ∈ R, define γ = 1, and α 1 = then combine the aforementioned five energy difference terms for the two successive iterations, finally, the following ultimate inequality can be resulted By selecting β 2 > 2 and α 2 > 2 λ 2 , (51) can be more simplified as follow where α 3 ∈ R + and α 4 ∈ R + . Therefore, it can be concluded that 93 ∆V k ≤ −κ(s k , s k−1 , I k , I k−1 ,f 2,k ,f 2,k−1 ,σ k ,σ k−1 , Q j,k−1 ).
Because κ is a class K function, V k is monotonically decreasing. In order to prove the bound of V k , by 94 using (8), (35), and (19), then the time derivative of V 0 is derived as below Based on (31) and (5),V 0 can be further simplified as It follows from (52) that V k is bounded. Accordingly, I k , Q j,k , σ k , f 2,k , and s k are also bounded.

98
Furthermore, based on the (52), it can be written that

Frontiers
Because V k is monotonically decreasing but it is lower bounded by zero, then the following conclusion is Based on (35), by considering (25) and the bounds of I k , Q j,k , σ k , f 2,k , and s k ,ṡ k is also bounded.

102
Accordingly, by Barbalat-like lemma presented in Xu and Yan (2004); Sun (2009), lim k→∞ s k = 0 103 uniformly on [t 0 , t]. Given that the sliding surface dynamics is Hurwitz, therefore, after s k converges to 104 zero asymptotically, the output error converges to zero exponentially.

APPENDIX C. VIRTUAL CONSTRAINT DESIGN
This section elaborates on desired trajectory generation using virtual constraints.
106 h(θ) in (22) is designed using the Bezier polynomials as (Westervelt et al. (2007)) In (59) Z is an integer, showing the number of Bezier polynomial terms, k is the parameter that is going 109 to be optimized, and w is calculated according to the following equation where Θ + and Θ − are maximum value and minimum value of the Θ(θ). Θ(θ) is defined as 111 Θ(θ) = ζ 1 θ 1 + ζ 2 θ 2 + ... + ζ n θ N where ζ i ∈ R should be chosen such that Θ(θ) is monotonically increasing. θ 1 , θ 2 , ...θ N represents for is depicted in Fig. 1. In the first step, the optimization is performed offline based on the system model. The

116
computed virtual constraint is embedded in the controller as the desired profile.

117
In this study, a sit-to-stand task was chosen for investigating the performance of the controller. For this 118 task, the knee joints' trajectories change monotonically, either decreasing or increasing. Therefore, both the 119 left or right knee joint angle can be selected as a phase variable. Each leg cannot have an independent phase 120 variable because it causes misalignment and miscoordination between the left and right legs. Therefore we 121 only selected the right knee joint angle as the phase variable to coordinate both legs' joints in this study.

122
Additionally, because all joints are actuated, one of the joints (here the right knee) follows a predesigned Figure 1. The optimization algorithm to determine the optimal Bezier polynomial coefficients time-dependent trajectory. All other joints use the state-dependent virtual constraint that is a function of 124 the actual right knee angular position (phase variable). It should also be noted that the right knee desired 125 velocity is also a function of the right knee phase variable.

APPENDIX D. MPC ALLOCATION ALGORITHM
The steps of the model predictive allocation algorithm can be found in Table 1.  Initialization: r = 0 (1a) The convergence tolerance is set to ε j . (1b) θ(t r ),θ(t r ) are measured. (1c) Feedback controller and virtual constraint are used to get h k , and total torque demand, where τ k ∈ [t r , t r + t p ]. (1d) An initial control trajectory is chosenū F,j,k (τ k ) ∈ U [tr,tr+tp] , where τ k ∈ [t r , t r + t p ].

2
Optimal Solution Searching: (2a) For solving the costates, integration backward in time is done for l    (ii) if r has exceeded the max iteration limit, N t , quit. (iii) otherwise r = r + 1 and reiterate gradient step from (1a).