A Robust Nonlinear Control Strategy for Unsupported Paraplegic Standing Using Functional Electrical Stimulation: Controller Synthesis and Simulation

Background: Functional electrical stimulation (FES) applies electrical pulses to paralyzed muscles to restore their function. Closed-loop control of an FES system, incorporating the control strategies, promises to improve the performance of FES systems. Therefore, the purpose of this paper is to design a new control strategy applicable to restoring the upright standing in paraplegic patients through FES. The control strategies proposed in the previous works based on controlling the angular joint position and none of them focused on controlling the center of pressure (CoP) dynamics directly. Since the CoP is representative of posture balance dynamics, in this study, the adopted FES based control strategy designed to control the CoP dynamics directly. Methods: In the proposed control strategy, two controllers determine the stimulation intensity of ankle muscles in a manner to restrict the CoP to a specific zone. According to the proposed strategy, until the CoP confined to the stable zone, an adaptive controller is active. When the CoP goes out of the stable zone, the adaptive controller is deactivated. Then, a sliding mode controller is activated instead of the adaptive controller. In this manner, not only the posture balance can be guaranteed, but also the emerged balance dynamics can be similar to the elicited balance dynamics in the healthy subjects. Results: In this study, extended evaluations carried out through the simulation studies. According to the achieved results, the proposed control strategy is not only robust against the external disturbances, but also insensitive to the initial postural conditions. Conclusion: The achieved results prove the acceptable performance of the proposed control strategy.


Introduction
The functional electrical stimulation (FES) applies electrical stimulation pulses to muscle nerves in order to produce muscle contractions. 1 Also, FES is a technique that widely used to restore lost motor function. 2,3 Unsupported standing using FES can enable the paraplegic patients to perform manual functions at self-selected leaning positions. 4 Also, standing brings other benefits to spinal cord injury patients such as decreasing bone osteoporosis, prevention of pressure sores, prevention of orthostatic hypotension, and improvement of the digestive system. [5][6][7] During recent years, different control approaches, applicable to restoring the arm-free standing in paraplegia via FES, have been proposed. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] Masani et al 8,9 introduced a feedback proportional-derivative controller to stabilize the body during quiet standing. The designed controller could facilitate stable standing even with closed-loop time delays longer than 75 ms. 8,9 Masani et al 10 also used a proportional integral derivative (PID) controller for balance control during standing due to FES. The used PID controllers determined the stimulation intensity of the ankle plantar flexor and dorsiflexor muscle groups. 10 Soetanto et al 11 presented an innovative control strategy for stabilizing the standing posture of paraplegic patients. The applied control strategy based on a combination of a proportional-derivative controller and a predictive mechanism. Huryn et al 12 assessed the viability and performance of an intermittent control strategy in human balance control. They showed that stability can improve by incorporating the prediction mechanism. 12 Hunt et al 13- 15 have also studied the performance of a control system with a nested-loop structure in which the outer loop was associated with regulation of the inclination angle, while the inner loop was associated with regulation of the total ankle torque. They used linear quadratic Gaussian (LQG), controllers. They also used the pole assignment method instead of LQG to achieve more extended periods of unsupported standing. 15 Matjacić and Bajd 16,17 proposed a new control scheme designed based upon a linear quadratic regulator optimal controller. They integrated the voluntary control of the unaffected upper limbs with the artificial control of the ankle joint movement. 16,17 The Mihelj and Munih 18 also addressed the integration of the voluntary trunk activity with artificial controlling the ankle torque. Audu et al 19 also proposed and designed a PID feedback controller for upright posture control using FES. Kobravi and Erfanian 20 proposed a new robust control framework based on the synergistic combination of an adaptive nonlinear compensator with sliding mode control. Adopting an agonist-antagonist muscle coactivation mechanism was one of the main contributions of their work. 20 Later, they utilized the proposed control approach 20 for unsupported paraplegic standing using FES. 21 In most of the mentioned works, the linear control structures associated with linear time-invariant systems have used, [8][9][10][11][12][13][14][15][16][17][18][19] while the musculoskeletal system has nonlinear and time-varying properties. 24 In most of them, the necessity of agonist-antagonist muscle co-activation was not regarded. [13][14][15] Also, no previous works focused on controlling the center of pressure (CoP) dynamics directly, while in the healthy subjects, the central nervous system controls the CoP directly to achieve a stable balance. Thus, this study aimed to design and evaluate a control strategy for controlling the unsupported paraplegic standing through controlling the CoP dynamics. 25 In other words, the goal of the designed control mechanism was realizing a dynamic balance rather than pushing the inclination angle toward a desired fixed position. In this manner, dynamic stability, instead of static stability, will be obtained. 26

Materials and Methods
The Model of Virtual Patient In this study, a musculoskeletal model used as a virtual patient. This virtual patient model includes the skeletal part and the muscular parts. Each part of the model will elaborate in the following subsections.

a) Skeletal Part
The skeletal part is an inverted pendulum that rotates around the ankle. This part has shown in Figure 1. 27 The dynamical equation for the inverted pendulum is given by (1), Where θ is the inclination angle, M is the body segment mass, I 0 is the moment of inertia, h is the distance of CoM from the ankle position, g is the acceleration of gravity, u(τ) is ankle torque, and ε(τ) is external disturbance torque. The set of parameters (i.e. M, I 0 , h and g) for the skeletal model have taken from 27 . In an inverted pendulum model of standing, the difference between the CoP and CoM is proportional to the horizontal acceleration of CoM in both the sagittal (anterior/posterior direction, A/P) and frontal (medial/ lateral direction, M/L) planes. Accordingly, if the horizontal acceleration of CoM is small, the CoP and CoM can be comparable and approximately equal. So, it can be supposed that horizontal displacement of CoM (y g ) and horizontal displacement of CoP (y p ) are equal. 28 It means that if ) sin(θ h y g = , y p can be defined as Accordingly, CoP y and the CȯP y will be defined by (2) and (3), respectively.
Where h is the position of CoM , and θ is the deviation angle. In reality, the CoP can estimate using force sensors data located on the shoe's insole.
by defining x 1 = Ɵ , x 2 =Ɵֹ , x 3= Ӫ using Euler scheme and substituting into (1), yields: Single joint inverted pendulum model in which h is distance of CoM from the ankle position, Ɵ is sway angle, u is total ankle torque, M is the body segment mass, g is acceleration of gravity, y g horizontal displacement of CoM, and y p horizontal displacement of CoP. 27 Int Clin Neurosci J. Vol During the quiet standing, the ankle strategy has applied to preserve the balance against the small external perturbation. Since this study focused on the preliminary evaluation of the proposed control strategy, the simulation studies and analyses carried out under the condition that the ankle strategy alone has applied. In this context, the hip joints and knee joints fixed in the fullextended position, and the balance preserved through activation and deactivation of the ankle extensor/flexor muscles. Accordingly, adopting the inverted pendulum rotating around the ankle joint, as the skeletal part of the musculoskeletal model, can be accountable.

b) Muscular Part
The muscular part of the virtual model consists of two muscle models with a different set of parameters. Figure  2 shows a schematic representation of the muscle model used in this study. Each muscle model includes muscle activation, muscle contraction, and segmental dynamics. 29 The recruitment curve modeled using a piecewise linear recruitment function. 29 A Gaussian function models the moment-angle relation, and the moment-angular velocity relation derived from a linear approximation of the forcevelocity relation described in. 29 The output of the muscle contraction block is the active muscle moment (M act ) computed by equation (5).

M
Act M M act m a mv = . . (5) Where M ma is the output of the moment-angle block, and M mv is the output of the moment-velocity block. At last, the total external moment in equation (6)  In order to have an uncertain model, we added a random term with uniform distribution to the intended parameters (some parameters of the skeletal part and some parameters of the muscular parts), which has the variance equal to the half of the real value of the parameter.

The Proposed Cooperative Control Strategy
The proposed strategy is a cooperative based control approach. The CoP lonely can consider as stability criteria. [27][28][29][30][31] For the sake of obtaining dynamic stability, a stable zone has determined so that confiding the CoP signal to this zone means eliciting a desired dynamic postural control. In other words, the proposed cooperative  29 The model describes the muscle activation properties using the muscle recruitment curve. The muscle contraction properties (moment-angle and moment-velocity), joint mechanical properties (joint-elasticity and joint-viscosity), and the leg dynamics have also described using this model. Unsupported Standing Control Using FES journals.sbmu.ac.ir/Neuroscience http based control strategy seeks to bring back the CoP to a specific zone and then seeks to prevent the CoP from going outside the zone. In other words, convergence toward the zero moment point has not desired. Figure 3 shows the structure of the proposed control strategy.
The decentralized control strategy has applied to deal with muscle co-contraction. Two different closed loop control systems envisioned for adjusting the stimulation pattern of each ankle muscle group (Figure 3). In the designed control strategy, whenever the COP signal and • COP signal are out of a specific elliptical boundary, two sliding mode controllers are used to adjust the stimulation pattern of ankle plantar flexor and dorsiflexor muscle groups, and whenever the CoP signal and • COP signal are within a specific elliptical boundary, two stable adaptive controllers are applied to adjust the stimulation pattern of ankle plantar flexor and ankle dorsiflexor muscle groups. In this manner, adaptive controllers, along with two sliding mode controllers assigned for adjusting the stimulation pattern of each ankle muscle group. It is worth noting that adaptation laws of the utilized adaptive controllers derived in a manner that can only guarantee the stability of the control system not asymptotically stability of control system. Therefore, moving toward the zeromoment point, which means obtaining static stability, is not the control goal. Such behavior is in conformity with the behavior emerging during normal upright standing. 30 As mentioned previously, an elliptical boundary has determine as a decision boundary concerning the instantaneous values of CoP and  Where a and b are horizontal and vertical diagonals of the ellipse. It has reported that during the quiet standing, the angular excursions of involved joints are less than degrees. 27 By assigning the mentioned range for the inclination angle and using equations (2) and (3) the stable intervals relating to COP and • COP will be as follows : The equation (10) shows the decision making mechanism related to activating two designed controllers and adjusting the stimulation intensity of the ankle muscle groups. are adaptive controllers' outputs. PW1 is the pulse width value of the delivered electrical pulse to the first muscle group (ankle dorsiflexor muscle), and PW2 is the pulse width value of the delivered electrical pulse to the second muscle group (ankle plantar flexor muscle).
The designed controllers will further elaborate in the next subsections.

Discrete Sliding Mode Control
Sliding mode control (SMC) has efficiently applied for controlling the various linear and nonlinear dynamical systems. 32 Since the musculoskeletal system has nonlinear dynamics, a nonlinear controller has to be utilized. 33 The sliding mode controllers are less sensitive Figure 3. The Schematic of the Structure of the Proposed Control Strategy. According to a decision-making mechanism, at each time instance, one of the adopted controllers (SMC or adaptive controller) is selected intermittently for the determination of muscle stimulation pulse width. The adaptive controller selected whenever the CoP signal and its time derivative value lie within a specific zone. Otherwise, the SMC will be activated. journals.sbmu.ac.ir/Neuroscience http concerning external disturbances and parametric uncertainties; therefore, they are appropriate for many real applications. 34 The design of SMC has two parts. The first part is the design of a sliding surface as a different guide the trajectories toward the origin. The second part is designing an appropriate control law. 35 The main feature of SMC is pushing the states of the system from the initial states toward a predefined sliding surface using a switching control law. 36 The musculoskeletal model used in this study consists of two muscle groups and a skeletal part. If each muscle-joint system has considered a subsystem, the musculoskeletal model is a coupled system. In this study, a unique sliding mode prediction model (SMPM) for a discrete-time nonlinear uncertain coupled system has presented. The approach introduced in, 37 has been utilized here. Consider the nonlinear system containing two coupled subsystems as follows 37 : is the state vector, u 1 (k) and u 2 (k) are the control signals deliver to two subsystems, x 1 (k) and x 2 (k) are the state variables of the first subsystems, x 3 (k) and x 4 (k) are the state variables of the second subsystems, f 1 (k), f 2 (k), g 1 (k) and g 2 (k) are nonlinear discrete functions, and d 1 (k), d 2 (k) represent dynamic coupling, parameter uncertainty, unmolded dynamics, and external disturbances.
To implement the control strategy, at first, a special SMPM is created as follows 37 : Where the s i (k) is sliding variable, x 2i (k) and x 2i-1 (k) are state variables of subsystems, and ϭ i is a constant. In the next step by applying feedback correction and receding horizon optimization, sliding mode control law is obtained as follows 37 : Where λ i is a tunable parameter, and h i (k) is a nonlinear discrete function which is derived as follows, g i (x(k)) is a nonlinear discrete function describing the uncertain dynamic system 37 : Where f i (x(k)) is a nonlinear discrete function describing the uncertain dynamic system, is computed using the available system information at each time instance. Also, ϭ i is a constant which its value can significantly influence on the stability and dynamic performance of the designed sliding mode controller. Also, γ i is designable parameter is the coefficient of the feedback correction term . More details about the utilized discrete sliding mode control can be found in. 37,38 Discrete Adaptive Control The design process of an adaptive controller usually contains three steps: the first step is choosing a control law, the second step is choosing an adaptation law for adjusting control parameters, and the third step is analyzing the output convergence properties of the resulting control system. 39 For the system defined by (1), the adaptive feedback control law and the update law are given by (15) and (16), respectively 40 .
Where K i (k) is control gain, which has been defined in appendix A in more details, F i (x(k)) is a nonlinear function, q is a positive constant, x i (k) is state variable, B 0 is a symmetric and sign definite matrix, Y i is a positive identity matrix, and P i satisfies the condition More details about A o and ∧ A , can be found in Fu et al 40 and Rugh, 41 respectively.
The stability of the used adaptive controller is guaranteed using the adaptation law (equation 16) based on Lyapunov stability theory 40 that has been explained in appendix B.

Results
Extended simulation studies were carried out on a musculoskeletal model as a virtual patient. The achieved results will be elaborated in details through the subsequent sections.

Effects of Initial Conditions
Two different initial conditions disturbing the balance of the musculoskeletal model ( ( ) The positive initial condition is equivalent to forward leaning and the negative initial condition is equivalent to backward leaning.
The results of using the sliding mode controllers have shown in Figure 4 and Figure 5. The results of using the adaptive controllers have shown in Figure 6 and Figure  7. The results of using the cooperative control strategy have shown in Figure 8 and Figure 9. The parameters of two sliding mode controllers selected as Table 1. 37,38 The parameters of the adaptive controllers have chosen heuristically, and according to some principles presented in the literature. [40][41][42] while ( ) , the parameters of two adaptive controllers were selected as follows: 90 u , a 2 10 = , , the parameters of two adaptive controllers were selected as follows:   , both flexor and extensor muscles contracted. However, the contraction of the extensor muscle happened sooner. Figures 5, 7, and 9 show that when the initial condition ( ) , both flexor and extensor muscles contracted. Nevertheless, the contraction of the flexor muscle happened sooner.
It could see that whenever the CoP and • COP are out of the desired ranges (defined by the equations (8) and (9)), the sliding mode controllers generate the control signal to contract the muscle (For example the time between 0.05 s to 0.1 s in Figures 8 and 9). Also, whenever the CoP and     The intriguing point that can address concerning to the applied adaptive controller is the dependence of the controller parameters on the initial condition of the body posture. Though the parameters of the adaptive controller depend on the initial inclination angle, such dependence on the initial body position did not influence the performance of the controller. As Figures 6 and 7 show, the CoP confined to the interval (-0.0758 < Cop ≤ 0.0758), which regarded as the stable interval.
In this study, because of the perturbed standing situations, CoP variations may seem fast. Clearly, the stable boundary of CoP ( ) 0.0758 CoP ≤ is less than the peak of CoP and the stable boundary of CoP velocity . 0 CoP is near to the peak of CoP velocity which have reported in. 43 The frequency of variation of the torque which has shown in Figures 8 and 9 is near to what has reported by Callahan et al, 44 some differences between real muscle and muscle model can be attributed to the inevitable unmolded dynamics of the chosen muscle model as a virtual muscle.

Effects of Pulse Form External Disturbances
In this step, two distinct pulse form external disturbances (one positive, one negative) exerted to the musculoskeletal model at the time 40 seconds after starting the control process. The time duration of each disturbance pulse was 1 second and its amplitude was 20 N.m (which was approximately 30% of peak-to-peak generated torque during the disturbance-free trial).
The results of using the sliding mode controllers have been shown in Figure 10 and Figure 11. The results of using the adaptive controllers have been shown in Figure  12 and Figure 13. The results of using the cooperative control strategy have been shown in Figure 14 and Figure  15.
The parameters of two sliding mode controllers were selected as Table 1.
While applying the positive external pulse, the parameters of two adaptive controllers were selected as follows: While applying the negative external pulse, the parameters of two adaptive controllers were selected as follows:       Figures 10, 12, and 14 show that while applying the positive pulse, both flexor and extensor muscles contracted. However, contraction of the extensor muscle happens sooner, or the value of the extensor torque was more. Also, Figures 11, 13, and 15 show that while applying the negative pulse, both flexor and extensor muscles contracted. Nevertheless, the contraction of the flexor muscle happens sooner, and the value of the flexor torque is more. It could see that whenever the CoP and • COP are out of the desired ranges (defined by the equations 8 and 9) sliding mode controllers generated control signal to contract the muscle (For example the time between 40.5 s to 40.6 s). Also, whenever the CoP and • COP are in the desired ranges, the adaptive controller generated the control signal to contract the muscle (For example: the time 41.12 s to 41.16 s).
Also, it could see that whenever the CoP and • COP are out of the desired ranges (defined by the equations 8 and 9), sliding mode controller activated again and generated the control signal to contract the muscle (For example the time between 40.61 s to 40.65 s in Figure 15). The performance of the cooperative control strategy evaluated, when two pulse form external disturbances (one positive, one negative) exerted at the time 40 s. The time duration of disturbance pulses was 4 s and their amplitude was 20 N.m. Figures 16 and 17 show the results. The frequency of variation of the torque which has seen in Figures 14-17 is near to what has been reported by Callahan et al, 44 some differences between real muscle and muscle model can attribute to the inevitable unmolded dynamics of the chosen muscle model as the virtual muscle.

Effects of Pulse Trains Form External Disturbances
In this step, two distinct pulse trains form external disturbances (one positive, one negative). It has exerted to the musculoskeletal model. The time duration of each disturbance pulse was 10 s, and its amplitude was 20 N.m. The results of using the sliding mode controllers have shown in Figure 18 and Figure 19. The results of using the adaptive controllers have shown in Figure 20 and Figure      Figure 22 and Figure 23.
The parameters of two sliding mode controllers selected as Table 1.
While applying the positive external pulse train, the parameters of two adaptive controllers were selected as follows:  While applying the negative external pulse train, the parameters of two adaptive controllers were selected as follows:

Quantitative Evaluations
At last, the performance of the three used control schemes compared. For this aim, three quantitative measures, including steady state error, settling time, and the number of times that the CoP and • COP entered and exited from the stable zone, were computed and analyzed.
When the predictive SMC used as a lonely control strategy, the control goal was making the output error converged to zero. Therefore, the steady state error was    Table 2 shows the computed values of steady state errors achieved when the SMC used as a lonely control strategy. The average value is less than 0.2°. Since the range of motion of the inclination angle can be between -90˚ to 90 ˚, such steady state error can be acceptable.
In the second step, the settling time was computed. The settling time is the time required for the system output (inclination angle) to reach within a boundary around the zero. The computed times have been reported in Table 3. Obviously, using the cooperative control strategy conduced to observe the shortest settling time (Average value: 3.7 s).
In the third stage, the number of times that the COP and • COP entered the stable zone and the number of times that the COP and • COP got out the stable zone were counted. Counting results have been shown in Table  4. Clearly, the most number of times which the COP and • COP were confined to the desired zone was observed when the adaptive controller was adopted. But, in terms of counting results, the performance of the adaptive controller and combined control method are comparable.

Discussion
In the previous subsections, the evaluation results of the proposed control strategy discussed from a technical point of view. Scrutinizing, in the achieved results, reveals some different aspects demonstrated in the following subparts.

Robustness Concerning Time-Varying Properties
As mentioned previously, the random numbers with uniform distribution added to some parameters of the musculoskeletal model. In this manner, the plant became a non-deterministic and time-varying system. The acceptable performance of the controller in such context proves the robustness of the control system concerning the time-varying properties of the musculoskeletal system. One of the notorious events that give rise to time-varying behavior is muscle fatigue. FES expedites the occurrence of muscle fatigue. Although the used muscle model does not simulate the occurrence of muscle fatigue, according to the results, it can expect that this controller can cope with muscle fatigue.

Imitating the Normal Standing Dynamics
In all the previous similar works, the control strategies mostly designed in a manner that pushes the inclination angle toward a set point to preserve the posture balance during the quiet standing. However, in this manner, the controller imposes a closed-loop dynamic that mimics  Unsupported Standing Control Using FES journals.sbmu.ac.ir/Neuroscience http the behavior which differs from the normal behavior of the musculoskeletal system observed during upright posture standing. In practice, imposing such closed-loop dynamics may at least result in a discomfortable standing bothering the patient. As emphasized previously, the main idea underlying the proposed control strategy is preserving the CoP and CoP velocity in a stable zone while their trajectories do not move towards the origin. In this manner, the balance dynamics can nearly imitate the balance dynamics of a healthy subject. Consequently, it can expect that utilizing such a control strategy not only can preserve the balance but also results in a comfortable standing of the patients.

Comparison to the Similar Works
As mentioned previously, controlling the CoP signal has not addressed directly in previous research works. Therefore, quantitative comparison to the similar presented works is cumbersome. However, the qualitative comparison seems rational. According to the results, the CoP dynamics controlled in a manner that it behaves in compliance with the CoP dynamics elicited in healthy subjects. However, using posture control strategies designed based on controlling the angle positions, as was done in the previous works, such dynamics conformity cannot be guaranteed. In fact, according to the achieved results, not only the joint angles were confined to the anatomical range of motion, but also the CoP dynamics has controlled directly.

Conclusion
In this work, a novel cooperative control strategy utilized for controlling the upright standing in paraplegia using FES. The cooperative control strategy designed based on a decision mechanism determined when the SMC or the adaptive controller has to be activated or deactivated. The results of the extended simulation studies on a musculoskeletal model, as a virtual patient, proves the acceptable and promising performance of the proposed strategy. Overall, it can conclude that the proposed cooperative control strategy is preferable in comparison to the single SMC and single adaptive controller.
Since the implementation of such experimental setup (implementation of electrical stimulators, proper hardware for bilateral communication with PC, related software, etc) was very time consuming, we would instead to present first the results of the comprehensive evaluation of the proposed control strategy on a musculoskeletal model. The next step of this research is providing a proper experimental setup for evaluation on the patients.