## Abstract

Functional electrical stimulation (FES) is prescribed as a treatment to restore motor function in individuals with neurological impairments. However, the rapid onset of FES-induced muscle fatigue significantly limits its duration of use and limb movement quality. In this paper, an electric motor-assist is proposed to alleviate the fatigue effects by sharing work load with FES. A model predictive control (MPC) method is used to allocate control inputs to FES and the electric motor. To reduce the computational load, the dynamics is feedback linearized so that the nominal model inside the MPC method becomes linear. The state variables: the angular position and the muscle fatigue are still preserved in the transformed state space to keep the optimization meaningful. Because after feedback linearization the original linear input constraints may become nonlinear and state-dependent, a barrier cost function is used to overcome this issue. The simulation results show a satisfactory control performance and a reduction in the computation due to the linearization.

## 1 Introduction

Functional electrical stimulation (FES) is a technique that can potentially restore function in impaired limbs by applying external currents to the nerves that innervate the muscles [1]. Tasks, such as standing and walking [2,3], or grasping and reaching [4,5], can be achieved by FES. Mainstream FES devices use transcutaneous electrodes to stimulate the nerves. The transcutaneous electrodes are applied to the surface of the skin [6]. FES via transcutaneous electrodes does not recruit muscles in a physiological manner, unlike voluntary muscle contractions that follow the Henneman's size principle. The principle states that slow-twitch and fatigue resistant muscles are recruited before the fast-twitch and less fatigue resistant muscles. During FES, this muscle recruitment order may be reversed, which is “nonphysiological” [7] and rapid onset of the muscle fatigue [8]. Feedback-based techniques usually overcome a fatigue-induced performance error by increasing stimulation amplitude or frequency, but this would only further aggravate the fatigue [9].

Recently, powered exoskeletons have been used to restore lower-limb function in individuals with spinal cord injury [10–13]. Powered exoskeletons can be reliably provide mobility for much longer duration compared to FES systems. For a longer usage, however, they may require larger batteries and heavier actuators to operate, and therefore may be cumbersome to wear [14].

To overcome the limitations of FES and the powered exoskeleton, a hybrid neuroprosthesis has also been proposed [15–23]. In the hybrid system, an electric motor assist can be incorporated to work along with FES so that the motor would help compensate for the muscle fatigue while the FES can reduce the battery size. However, how to optimally allocate control inputs between FES and the electric motor, based on a performance measure that minimizes muscle fatigue, limb angle error, and control effort, is an open research problem.

Model predictive control (MPC) is a widely used optimal control technique that can solve the aforementioned allocation problem [24–28]. MPC has a nominal model that enables the controller to predict the system behavior and accordingly plan an optimal control trajectory [28–30]. In our previous work [31], an MPC method based on a fast gradient-based search algorithm [30] was used to optimally determine the FES to achieve a leg extension motion with a hybrid neuroprosthesis. In Ref. [28], a motor was added to assist the FES. The prediction horizon was set to be 0.5 s, while the sampling frequency was 100 Hz. The limited prediction horizon length is due to the heavy computational load of the MPC. The computational load is mainly determined by the complexity of the Hamiltonian, which contains a running cost function and dynamical constraint. Therefore, a highly nonlinear and dynamical constraint, which represents the complex neuroprosthetic system, would cause a heavy computational load. Additionally, the nonlinear dynamical constraint may trap the optimization search at a local minimum [29,32]. This problem motivated us to linearize the system before implementing the MPC method.

Therefore, we are motivated to simplify the Hamiltonian to reduce the computational load. Musculoskeletal system can be linearized at the equilibrium point by using Jacobian linearization [24,33], but the linearized model can only be considered accurate near its linearization point. Feedback linearization method, however, applies exact linearization. Therefore, the model's validity is not restricted near a prescribed point [34]. When the system is feedback linearized the Hamiltonian is simplified, and due to convexity a global optimum to the optimization problem may be found.

The structure of the MPC method on controlling feedback linearized FES-Motor-driven neuroprosthetic system was studied in Ref. [26], where the prediction horizon length has been prolonged to over 1 s. However, when the inputs are linearly constrained the input constraints may become nonlinear and state-dependent after feedback linearization. In Ref. [26], the input constraints were not added in the prediction horizon. It means that if the search touches the boundary of input constraint, the finally solution may not be optimal. This was mentioned but not solved in Ref. [26], where only a saturation filter was added on the external loop.

This problem in the twisted constraint was first mentioned in Ref. [35]. Of course, if the original input is not constrained, the original optimization problem can be reduced to a linear quadratic regulation problem or can be solved as a convex optimization problem. For the hybrid neuroprosthetic system, however, because the muscle's response to the stimulation current is limited, a constraint on the FES must be imposed. Consequentially, due to feedback linearization, a nonlinear and state-dependent input constraint must be then considered. To deal with the constraint issue, some solutions were investigated in Refs. [35–37], but these methods may be inappropriate for an FES application because they may considerably reduce the size of the feasible region. To overcome this difficulty, a barrier function, which penalizes the input [38] if it is beyond a physical bound (by exponentially increasing the cost), is developed in this paper. Thus, the MPC algorithm can solve the optimization problem without imposing input constraints.

In this paper, we extend the previous work [26] by designing a barrier function to deal with the input constraints. The effect of the barrier function is examined, and the performance and computation time is compared with the original nonlinear system. Furthermore, despite feedback linearization, the state variables, include the angular position and the muscle fatigue, are still preserved in the state space. This consequentially keeps the optimization problem meaningful. Simulation results indicate that the gradient search-based MPC method with the barrier input function when applied on the feedback linearized human knee musculoskeletal system is able to achieve the control task satisfactorily. Further, it constrains the control signal and significantly reduces the computation time.

## 2 System Dynamics

where $J\u2208\mathbb{R}+$ and $m\u2208\mathbb{R}+$ denote the moment of inertia and the mass of human lower leg, respectively, $g\u2208\mathbb{R}+$ is the gravitational acceleration, $l\u2208\mathbb{R}+$ is the length from the knee to the center of the mass, $\varphi ,\varphi \u0307,\varphi \u0308\u2208\mathbb{R}$ denote the anatomical knee joint position, velocity, and acceleration, respectively, $\tau p\u2208\mathbb{R}$ is the damping and stiffness of the joint in human leg, i.e., passive musculoskeletal torque of the knee joint, $\tau ke\u2208\mathbb{R}$ is the unidirectional torque from human leg elicited by FES, and $\tau m\u2208\mathbb{R}$ is the electric motor torque.

*j*=

*1, 2, 3), and $\varphi 0\u2208\mathbb{R}$ are subject specific parameters,*

*u*∈ [0, 1] is the normalized electrical stimulation amplitude, and the fatigue variables, $\mu \u2208\mathbb{R}$, are modeled as [40]

_{ke}where $Tf\u2208\mathbb{R}+$ and $Tr\u2208\mathbb{R}+$ denote the muscle fatigue time constant and the recovery time constant, respectively. The fatigue variable *μ* ∈ [*μ*_{min,}*μ*_{max}], where *μ*_{max} = 1, which represents no fatigue, while $\mu min\u2208\mathbb{R}+$ means that the muscle has completely fatigued. Here, we also assume that 0 < *μ*_{min} < 1 [9].

where *α* = 1/*J*, *β* = *mglα*. In Eq. (4), and $u\u2208\mathbb{R}2\u2009u=[u1,u2]T=[\tau m,uke]T$.

## 3 Feedback Linearization

The conditions for performing feedback linearization are to check if the relative degree equals the order of the system. Relative degree, which is defined as the times of Lie derivation is taken on the system until the input explicitly occurs. For a multiple input–multiple output system to be input-state linearizable, the dimension of the input must be equal to the output in order to linearize this system [34,41,42].

It can be seen that in the original system the relation between the input *u* and the state *x* is nonlinear. After linearization, the input *v* and the state *z* show a linear relation. The controllability of the linearized system can be verified by inspecting the rank of the controllability matrix.

## 4 Model Predictive Control Development

### 4.1 Control Loop.

*z*and the transformed input

*v*are generated by MPC by minimizing a cost function of the linearized system in Eq. (9). The new pair of desired state

*z*

_{des}and input

*v*

_{des}are computed from the desired state

*x*

_{des}, desired input

*u*

_{des}, and the nonlinear system in Eq. (4)

The optimal control trajectory *v* is the input to the linearized system, which is the human knee system cascaded with the controller *u*(*x*, *v*). As this *u* is computed by considering a subsequently defined soft constraint, a saturation filter is added before it goes into the actual physical system.

### 4.2 Barrier Cost Function.

*v*becomes nonlinear and state-dependent. The motor torque is assumed to be unconstrained [43], but for FES, the input is normalized and confined within [0, 1]. This motivated us to add a barrier function. The barrier function is defined as

where *N _{i}*, for

*i*=

*1, 2, 3, 4 are user-defined positive integers. From Eq. (7), we can get Bar(*

*v*) by substituting

*u*by

*v*.

*Remark.* The parameters, *N _{i}*, are usually tuned to be large so that the penalty cost can prevent the control signal from deviating too much from the actuator's physical constraint. However, in a numerical search, large parameters may cause the cost to go to infinity. Therefore, there is a trade-off in the parameter selection. In practice, the parameter values can be selected from simulations.

### 4.3 Model Predictive Control Problem.

*k*∈ [0,

*T*], where [0,

_{N}*T*] denotes the period of the whole control process. The initial state of the nominal system $z\xaf(0)$ is set to be

_{N}*z*for the prediction horizon

_{k}*t*∈ [0,

*T*] starts at real time

*k*. The objective of the optimal control problem is to solve for the control input $v\xaf*(t)$ where

*t*∈ [0,

*T*], which minimizes the cost function $C(z\xaf,v\xaf)\u2208\mathbb{R}$ and the resulting optimal nominal state trajectory is defined as $z\xaf(t)*$ where

*t*∈ [0,

*T*]. Afterward, $vk=v\xafk*(\tau )$ will be applied to the real system. Here,

*τ*∈ [0,

*δt*] and

*δt*is the time that MPC takes to apply the control input. The integral cost function $l(z\xaf,v\xaf)$ and the terminal cost $V(z\xaf(T))$ are defined as

*v**(

*t*) for

*t*∈ [0,

*T*) solves the optimal control problem if it minimizes the Hamiltonian, where the Hamiltonian is defined as

where $\lambda (t)\u2208\mathbb{R}n$ is the costate vector.

To solve the optimal control problem (11), a gradient projection algorithm [30,44] is adopted. The detailed computation steps of the gradient projection algorithm can be found in Ref. [30], and summarized as in Table 1 for each time instant *k*.

One condition has to be satisfied is that the initial input must be in the physically feasible region of the actuators.

#### 4.3.1 Benefit of Using a Barrier Cost Function.

A contour of the cost function at time *k* is visualized in Fig. 3(a). In this figure, the state space is assumed to be one-dimensional, and the input is two-dimensional (2D), thus the cost function contour is three-dimensional. The dynamical constraint and the input constraints are 2D. The optimal solution search starts on a dynamical constraint plane, and at the optimal solution, the cost contour and the dynamical constraint plane are tangent to each other. If the solution is not found, the search proceeds to an input constraint plane. In Fig. 3(b), the three-dimensional contour is reduced to 2D. From this figure, we can notice that in the original optimization problem the dynamical constraint is nonlinear and the input constraint is linear. However, after the feedback linearization the dynamical constraint becomes linear and the input constraint becomes nonlinear (see Eq. (6)).

Therefore, if an optimum point is found on the boundary of the input constraint, the optimum might still be local due to the nonlinearity. Moreover, the computational burden may not be necessarily reduced if the search proceeds to the nonlinear input constraint. According to Eq. (6), the computation may be infeasible if the transformed input constraint becomes state-dependent. This is because during the search the input constraint can also change (see Fig. 3(c)) due to its state dependency. The barrier function can help by converting the input constraint to a penalty in the cost function (a soft constraint). Of course, when the barrier function is added, the cost contour near the barrier may also be state-dependent (see Fig. 3(d)). However, from a practical point of view the latter case is much easier to solve than to specify the state-dependent input constraint in the prediction horizon.

## 5 Results

The control loop was realized by using simulink of matlab while incorporating a gradient search-based MPC optimization algorithm composed in a *C* code. Discrete time-step was set to be 0.01, and the prediction horizon was set to be 1.2 s (the prediction was 0.5 s in our previous studies of nonlinear MPC [28,31], but it was prolonged in this research due to the linearization). In the cost function, the weight matrices were set to be *Q* = diag (680, 15, 4), *R* = diag (1, 80), and *P* was solved from the $PA+ATP\u2212PBR\u22121BTP+Q=0$. Solving *P* in this fashion is able to ensure the MPC stability when dealing with a regulation task [30,45]. In the barrier cost function, the parameters are set to be *N*_{1} = *N*_{3} = 8 and *N*_{2} = *N*_{4} = 2. The simulation was run to validate the method. In Fig. 4, the knee angle is regulated close to the desired position by the inputs computed from the MPC method. Figure 5 shows that the muscle fatigue index stays around 0.6–0.7 instead of dropping to 0. Figure 6 indicates that the deviation from the saturation bound is fairly small, but if the barrier function is turned off, the deviation from the physical bound could be huge. This validated the effectiveness of using the barrier function in the cost function. The motor torque is shown in Fig. 7. Moreover, in Fig. 8, the computation time of the linearized system is shown to be significantly reduced with a comparable control performance. (It is shown in Fig. 9, and details are summarized in Table 2, where the settling time is defined to be the time when the error is less than 5 deg while the overshoot is the maximum value above the reference. However, the nonlinear MPC used for comparison usually need over 0.02 s execution time at each time instance while the sampling time is 0.01 s, which means that it cannot be used although it shows a better performance than the linearized one.)

## 6 Conclusion

In this paper, an MPC method for a feedback-linearized hybrid neuroprosthetic system driven by FES and an electric motor-assist was designed and simulated. A barrier function is added in the cost function as a soft constraint. Simulations were run to validate the MPC method with the barrier function. As a result, the knee joint regulation was achieved, the onset of the FES-induced muscle fatigue was well controlled, and the FES was confined to the prescribed input constraint. Importantly, the MPC method needed a lower computational time when the feedback linearized model was used compared to the case when the original nonlinear model was used. The barrier cost function was shown to be able to prevent the control input from deviating away from the actuator's physical limitation. However, because the feedback linearized model is sensitive to the disturbance in the model parameters, it is less robust than the nonlinear model. This problem may be addressed by incorporating with a robust MPC method, which will be investigated in the future work. Overall, a potential implication is that it can facilitate a real-time implementation of this method in controlling a hybrid neuroprosthesis with a reduced computation load.

## Funding Data

National Science Foundation (Grant No. 1511139; Funder ID: 10.13039/501100008982).

National Institute of Health (Grant No. R03HD086529; Funder ID: 10.13039/100000002).