Rehabilitation through functional electrical stimulation (FES) is a treatment for people with neurological conditions (NCs), such as stroke and spinal cord injury [1] and [2]. FES induces involuntary muscle contractions to perform a functional movement by applying an electric potential across the motor neurons of a muscle. To improve motor function and overall quality of life, multiple efforts in the rehabilitation field use FES with rehabilitation robots to facilitate human-robot therapy [3]. Stationary FES cycling is a common human-robot rehabilitative therapy for people with movement impairments resulting from NCs [4]. FES cycling has the benefits of both FES and rehabilitation robotics; it is a preferred therapy because there is minimal risk of a fall, and the repetition of coordinated limb movements improves motor skills and nervous system reorganization [5].
Optimal controllers can be established by assigning a userdefined cost to the states and control inputs, which penalizes the state and the magnitude of the control input. Through the cost function, a balance can be obtained between the accuracy of the limb motion versus the level of control effort, allowing potential tradeoffs between comfort, performance, duration of exercise, and muscle fatigue. The only results that apply optimal control methods to FES applications are Wanjiku A. Makumi, Max L. Greene, Kimberly J. Stubbs, and Warren E. Dixon are with the Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL, USA. Email: {makumiw, maxgreene12, kimberlyjstubbs, wdixon}@ufl.edu.
This research is supported in part by NSF Award number 1762829, Office of Naval Research Grant N00014-13-1-0151, AFOSR award number FA9550-18-1-0109, and AFOSR award number FA9550-19-1-0169. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the sponsoring agency. [6] and [7]. These results use extremum seeking, a modelfree online optimization tool, to adjust a closed-loop PID controller to minimize the cost function for upper limb electrical stimulation.
Optimal control problems can be solved via the Hamilton-Jacobi-Bellman (HJB) equation [8]. By solving the HJB equation to determine the optimal value function, an optimal control policy can be developed [8]. Generally, the HJB equation does not have a closed-form analytic solution for nonlinear systems. Motivated by the challenges of solving the HJB, especially in real-time, approximate dynamic programming (ADP) has emerged as a method to yield an approximate solution. Specficially, ADP uses a reinforcement learning (RL)-based actor-critic framework to approximate the value function in real-time [9]. Neural networks (NNs) are generally used within ADP to approximate the unknown optimal value function, but other function approximation methods could also be used [10].
In traditional adaptive control, the uncertain parameter estimates are updated using an error feedback as a performance metric; in ADP, the Bellman error (BE) is used as feedback on the level of suboptimality. Specifically, the BE is used to update the NN parameters to improve the value function approximation online. BE extrapolation yields faster policy learning over a domain by evaluating the BE over user-defined, off-trajectory regions of the state space [11]. Sufficient off-trajectory data must be selected to achieve adequate exploration. The value function approximation is updated according to the on-and off-trajectory BE.
Due to the potential benefits of using an optimal controller, it is advantageous to apply ADP to the cycling system. However, the system switches between two actuation methods: the rider's muscles and the cycle's electric motor. Therefore, FES-cycling is a switched (also called hybrid) system, which requires switched (hybrid) system analysis and design methods [12]. Until recently, switching has not been investigated in the context of ADP. The result in [13] develops a framework to estimate the optimal feedback control policy online while switching between multiple dynamic system models. When analyzing switched systems, a common problem is the growth and discontinuity of Lyapunov functions at switching instances [14]. This growth and discontinuity problem is overcome in [13] in which a dwell-time analysis is developed to determine the minimum time necessary before the system can switch to a different subsystem (i.e., a minimum dwelltime). This provides a framework to switch between the two different modes of the FES-cycling controller and show stability of the overall switched system.
Motivated by our previous results in [13] and [14], this paper implements a continuous-time ADP-based tracking controller that allows for switching between mutliple cycle actuation methods to track a desired cadence. Uniformly ultimately bounded (UUB) stability of the overall switched system is proven. Moreover, the developed controller is also proven to converge to a neighborhood of the optimal controller.
For notational brevity, time-dependence is omitted while denoting trajectories of the dynamic systems. For example, the trajectory x (t), where x : R ≥0 → R n , is denoted as x ∈ R n and referred to as x instead of x (t). For example, the equation f + h (y, t) = g (x) should be interpreted as f (t) + h ((y, t) , t) = g (x (t)). The gradient denote the modulo operator where, generally, m is the dividend, p is the divisor, and r is the remainder. In this paper, the quantity or function belonging to the k th mode of the switched system is denoted with the subscript k.
Following the development in [15], the dynamics of the combined one-legged cycle and rider system are τ = M (q) q + V c (q, q) q + G (q) + P (q, q) + b c q, (1) where q, q, and q ∈ R denote the angle, angular velocity, and angular acceleration, of the crank arm respectively. M : R → R >0 denotes the inertia matrix, V c : R × R → R denotes the centripetal-Coriolis matrix, G : R → R denotes the gravitational effects, P : R × R → R denotes the passive viscoelastic tissue forces, b c ∈ R >0 denotes the cycle's viscous damping effect, and τ denotes the torque applied by the quadriceps muscle and the cycle motor, which is subsequently defined.
The torque is applied by two different actuators, corresponding to either the torque due to the FES-induced muscle contractions or the torque due to the electric motor. Given the need to use the different actuators at different times, we define two sets: Q, when the crank angle is in the kinematically effective quadricep region, and Q c , when the crank angle is in the region of poor kinematic efficiency [16]. Let Q ⊂ [0 • , 360 • ) denote where electrical stimulation is active and Q c denote the complement of Q, where the electric motor is active.
The torque τ : R × R → R in ( 1) is defined as
where b 1 : R × R → R >0 is the assumed known muscle control effectiveness, u 1 ∈ R is the muscle control input, b 2 ∈ R is the known motor control constant, and u 2 ∈ R is the motor control input. From (2) the dynamics for each mode are [15] b k u k = M (q) q + V c (q, q) q + G (q) + P (q, q) + b c q, (3)
where k represents the active switched subsystem. Let k ∈ S, where S ≜ {1, 2} is the switching index set.
Following the development in [17], the dynamics in (3) can be rewritten in the control-affine formfoot_0
where x ≜ [q, q] T , and a subsequently defined control input u k ∈ R represents the control input for the k th system. The drift dynamics f : R 2 → R 2 are defined as f (x) ≜ q M (q)
-1 (-V (q, q) q -G (q) -P (q, q) -b c q) , and the control effectiveness g k : R 2 → R 2 is defined as
The control objective is to track a time-varying continuously differentiable signal x d ∈ R 2 . To quantify the tracking objective, the tracking error e ∈ R 2 is defined as e ≜ x -x d . Using the technique in [17], the control affine dynamics in (4) can be expressed as
where
The following properties and assumptions facilitate the development of the desired approximate optimal tracking controller.
Property 1. The drift dynamics f are continuously differentiable [15], which, using [18,Lemma 3.2], means that f is a locally Lipschitz function and f (0) = 0.
Property 2. The control effectiveness matrix g k is continuously differentiable [15] and therefore a locally Lipschitz function [18,Lemma 3.2]. The matrix g k is bounded such that 0 < ∥g k (x)∥ ≤ g k ∀x ∈ R n , where g k ∈ R >0 is the supremum over all x of the maximum singular value of g k (x), respectively, for all k. It follows that ∥G
Assumption 1. The desired trajectory is upper-bounded by a known positive constant
Based on the above assumptions, the trajectory tracking component of the controller
The control objective is to solve the infinite-horizon optimal tracking problem i.e. to find a control policy µ k that minimizes the cost function
where
and Q k is PD so that the cost in (6) does not depend on the desired trajectory.
Property 3. The state cost matrix Q k satisfies q k ≤ Q k ≤ q k where q k , q k ∈ R >0 are the minimum and maximum eigenvalues of Q k , respectively.
The infinite horizon value function (i.e. the cost-to-go) for the k th mode V * k : R 4 → R ≥0 is defined as
where U ⊂ R is the action space for µ k .
Assumption 3. The optimal value function V * k is continuously differentiable for all k ∈ S [17].
The optimal transient control policy µ * k : R 4 → R is defined as
Each k th optimal value function and optimal control policy satisfy the HJB equation
which has the boundary condition V * k (0) = 0.
The optimal value function V * k is unknown for general nonlinear systems. Let Ω ⊂ R 4 be a compact set such that ζ ∈ Ω. The value function can be approximated with a NN in Ω by invoking the Stone-Weierstrass Theorem to obtain
where W k ∈ R L is a vector of unknown weights, ϕ : R 4 → R L is a user-defined vector of basis functions, and ϵ k : R 4 → R is the bounded function reconstruction error. 2 Substituting (10) into (8) yields a NN representation of the optimal control policy The ideal weights W k are unknown a priori; hence, an approximation of W k is desired. The critic weight estimate vector Ŵc,k ∈ R L is substituted into (10) to obtain the optimal value function estimate Vk :
The actor weight estimate vector Ŵa,k ∈ R L is substituted into (11) to obtain the optimal transient control policy estimate μk : R 4 × R L → R, defined as
The overall controller u k ∈ R is defined as
To calculate the BE δ k : R 4 × R L × R L → R, the optimal value function V * k (ζ) and the optimal control policy µ * k (ζ) in ( 9) are replaced by the approximate optimal value function Vk ζ, Ŵc,k and the approximate optimal control policy μk ζ, Ŵa,k , respectively, where
The value of the BE indicates how close the actor and critic weight estimates are to their respective ideal weight values. By subtracting ( 9) from ( 14), substituting (10)- (13), and denoting the difference between the actual and ideal weight values by Wc,k ≜ W k -Ŵc,k and Wa,k ≜ W k -Ŵa,k the analytical form of the BE in ( 14) is
where
Although they are equivalent, ( 14) is used in implementation and ( 15) is used in the stability analysis.
Using the control policy given in ( 13), the current system state, the critic weight estimate, and the actor weight estimate, the estimated BE in ( 14) can be evaluated to calculate the instantaneous BE denoted by δ k ζ, Ŵc,k , Ŵa,k at each time instance t ∈ R ≥0 . The exploration versus exploitation problem is well-known for learning-based control methods. In results such as [21], an exploration signal is required to successfully explore the operating domain. Results such as [11] use BE extrapolation, which simultaneously evaluates the BE along the system trajectory and at user-defined points in the state space. The BE extrapolation technique eliminates the need for the signal by providing simulation of experience, thus yielding a better value function mation [11].
The BE is extrapolated from the user-defined off-trajectory points {ζ i :
i=1 set by the user, where N k ∈ N denotes a user-specified number of overall extrapolation trajectories in the compact set Ω. The tuple (Σ c,k , Σ a,k , Σ Γ,k ) represents the data stacks defined as Σ c,k ≜ 1
is a user-defined gain, and Γ k : R L×L is a time-varying least-squares gain matrix. Each subsystem has its own distinct set of data, gain values, and update laws. Assumption 5. On the compact set, Ω, a finite set of offtrajectory points {ζ i :
where c k is a constant scalar lower bound of the value of each input-output data pair's minimum eigenvalues for the k th subsystem [11].
The critic and actor weights are updated according to the subsequent laws while each mode is active. In the weight update laws, η c1,k , η c2,k , η a1,k , η a2,k , λ k ∈ R are positive constant learning gains, and Γ k , Γ k ∈ R >0 are the upper and lower projection operator bounds for Γ k . The critic update 3
where
The actor update law for the k th mode Ẇa,k ∈ R L is defined as
where
Ŵc,k +η c2,k Σ a,k Ŵc,k . The operator proj {•} denotes the smooth projection operator defined in [22, Appendix E, Eq. E.4] and is designed such that Ŵc,k ∈
where 1 {•} denotes the indicator function. 4 While the k th mode is inactive Ẇc,k = 0 L×1 , Γk = 0 L×L , and Ẇa,k = 0 L×1 . 5Remark 2. Under Assumptions 1-3, the PD solution of the HJB equation is the optimal value function for each system. The approximation of the PD solution to the HJB equation is guaranteed by the appropriate selection of Lyapunov-based update laws and initial weight estimates [23].
It is possible for a switched system to become unstable, even if the individual subsystems of a switched system are stable [12,Ch. 3]. Hence, the stability of each subsystem must be investigated along with the switching between the systems. In the subsequent development, k subsystems, each with a class of dynamics in (4), are analyzed with the control policy in (13) and update laws outlined in ( 16)- (18).
Since the state penalty matrix Q k is PSD, the optimal value function V * k is PSD and is not a valid Lyapunov function. However, a nonautonomous form of the optimal value function denoted as
, and is PD and decrescent [17]. To facilitate the stability analysis, let z k ∈ R 2+2L be a concatenated state vector defined as
According to [17] and [18,Lemma 4.3], (19) can generally be bounded as
To facilitate the subsequent dwell-time analysis, the following more restrictive assumption is required.
Assumption 6. The optimal value function V * t,k (e, t) can be bounded by the square of the norm of its argument times a positive constant, i.e., 6
for all k ∈ S, where
Using Assumption 6, ( 19) can be bounded as
where
The normalized regressors ω k ρ k and
, and G ϕ,k is bounded as
Remark 3. Using the projection operator from the critic update law in (16) and [22,
Using the projection operator from the actor update in (17) and [22,Lemma E.1], -W T a,k Ẇa,k is bounded from above as
To facilitate the subsequent analysis, let R ∈ R >0 be the radius of a compact ball B R ⊂ R 2 × R L × R L centered at the origin.
Theorem 1. While each subsystem is active, if Assumptions 1-6 hold, the control policy in (13) and the weight update laws in ( 16)-( 18) are implemented, and the conditions converges to a neighborhood of the optimal control policy µ * k . Proof: Using (5) and the fact that V * t,k (e, t) = V * k (ζ) , ∀e ∈ R 2 , t ∈ R ≥0 and taking the time derivative of the candidate Lyapunov function in (19) yields
where the fact that d dt Γ -1 k = Γ -1 k Γk Γ -1 k is used. Under the sufficient gain conditions in (21) and (22), and using ( 9), (15), and the update laws in ( 16)-( 18), the expression in (24) can be bounded as
for all k ∈ S and t ∈ R ≥0 , where
, and l k is a positive constant that depends on the control gains and NN bounding constants in Assumption 4. Using the bounds in (20), the time derivative in (24), Λ k , and ( 23), [18,Thm. 4.18] can be invoked to prove that z k is UUB such that lim sup t→∞ ∥z k ∥ ≤ 2α 2,k l k α 1,k Λ k , and the transient tracking control policy μk converges to a neighborhood of the optimal control policy µ * k . Since z k ∈ L ∞ , it follows that e, Wc,k , Wa,k ∈ L ∞ , and since μk ∈ L ∞ and x d ≤ x d , it follows that u k ∈ L ∞ . Furthermore, every trajectory z k that is initialized in the ball B R is bounded such that z k ∈ B R , ∀t ∈ R ≥0 , ∀k ∈ S. Since z k ∈ B R , it follows that the individual elements of z k lie in a compact set, i.e. e, Wc,k , and Wa,k lie in a compact set. Additionally, since x d ≤ x d , then the concatenated state ζ ∈ Ω, ∀t ∈ R ≥0 , ∀k ∈ S, which facilitates value function approximation. Remark 4. See [11] for insight into satisfying the gain conditions in (21) and (22).
Let t ON k ∈ [0, t] denote a time instant when the k th subsystem of the switching sequence is activated. Let t OF F k ∈ [0, t] denote a time instant when the k th subsystem in the switching sequence is deactivated. The dwell-time in any active mode of a subsystem denoted by τ k ∈ R >0 is defined as τ k = t OF F k -t ON k and represents the amount of time a subsystem must be active before switching to the next. The minimum dwell time for any active mode of a system is denoted by τ * ∈ R >0 . There are a finite number of switches, and N σ ∈ N <∞ denotes the number of switching events. The sequence of time instants at which a switching event occurs is defined as t ON Nσ , such that 0
The stability analysis proves that each subsystem is UUB while active. However, the Lyapunov function for the overall switching system may instantaneously increase due to the change in the optimal value function and set of new weight parameters. The value function corresponding to mode k +1,
, may be larger than the value function corresponding to mode k, V * t,k (ζ). Similarly, the magnitude of the actor and critic weight errors could be larger in mode k + 1 than in mode k. Therefore, a dwell time condition must be designed to account for switching between the subsystems which ensures that the switched system is stable [12,Ch. 3].
Theorem 2. The system consisting of a family of subsystems with the dynamics in (4) with a properly designed minimum dwell-time, τ * ∈ R >0 ensures that the tracking error, critic estimate errors, and actor estimate errors will converge to a neighborhood of the origin in the sense that ∥z
is the maximum ultimate bound for all subsystems, and T ∈ R ≥0 is the time required to reach the ultimate bound.
The proof follows that of [13,Theorem 2] and is available upon request.
Remark 5. From Section II, the system switches modes based on if the crank angle q belongs to Q or Q c , i.e., the switching is state-based; furthermore the switches occur more frequently at higher desired cadence values. The user cannot directly control the time of the switching instances. For the system to be stable using the previous analysis, the dwell-time must be significantly smaller than the time required to travel through the regions Q and Q c . The dwelltime τ * is composed of many user-selected parameters. Notably, τ * can be decreased by increasing the decay rate γ 0 , i.e. stronger convergence parameters result in a shorter dwelltime. The dwell time τ * is inversely proportional to the decay rate γ 0 , and γ 0 is proportional to Λ. Hence, maximizing Λ will decrease the dwell time. This is achieved by maximizing the term min λ min (Q k ) , 1 6 η c2,k c k , 1 8 (η a1,k + η a2,k ) . While this maximization decreases the dwell-time so that it is significantly smaller than the time dictated by the desired trajectory, there are some practical drawbacks. A larger state cost matrix Q k will increase the penalty on the error; paired with larger actor and critic learning gains η a1,k , η a2,k , and η c2,k , this may lead to a more aggressive controller, which may cause rider discomfort. Furthermore, increasing c k relies on using more BE extrapolation data pairs, which may become computationally intensive. Motivated by these practical considerations, additional analysis methods that can potentially eliminate the need for a minimum dwell-time are motivated.
This paper develops an ADP-based controller for switched cycle dynamics while achieving a time-varying tracking objective. The stability of each individual subsystem is proven by a Lyapunov-based analysis, and the stability of the overall switched system is proven via a dwell-time analysis. The entire switched system is proven to be UUB such that the control policy is proven to converge to a neighborhood of the optimal policy and to track the cadence within a neighborhood of its desired value. Future work will investigate the application of the developed controller to an FES-cycling testbed and the development of analysis methods free of minimum dwell-time requirements.
The cycle-rider dynamics do not differ between modes. The only difference between the switching modes is the actuation methods.
To focus the scope of this manuscript, each switched system will use the same dimension vector of basis functions ϕ (ζ) i.e., L 1 = L 2 = ... = L k
Using (18) ensures that each Γ k ≤ ∥Γ k ∥ ≤ Γ k for all t ∈ R >0 .
The update laws will not update a subsystem k's weight estimates or least-squares matrix unless subsystem k is active.
Authorized licensed limited to: of Florida. Downloaded on September 12,2022 at 15:09:26 UTC from IEEE Xplore. Restrictions apply.
Authorized licensed use to: University of Florida. Downloaded on September 12,2022 at 15:09:26 UTC from IEEE Xplore. Restrictions apply.