Abstract
This manuscript presents a Takagi–Sugeno fuzzy control for a mathematical model of the knee position of paraplegic patients using functional electrical stimulation (FES). Each local model of the fuzzy system is represented considering norm-bounded uncertainties. After obtaining the model of FES with norm-bounded uncertainties, the fuzzy control strategy is designed through the solution of linear matrix inequalities (LMIs) using the conditions available in the literature, which consider these norm-bounded uncertainties. The strategy considers decay rate and constraints on the input signal. The model is simulated in the Matlab environment using the numerical parameters measured by experimental tests from a paraplegic patient.
1. Introduction
Several researchers have used functional electrical stimulation (FES) to restore some motion activities of people with injured spinal cord [1]. However, FES is not yet a regular clinical method because the amount of effort involved in using actual stimulation systems still outweighs the functional benefits they provide. One serious problem of using FES is that artificially activated muscles fatigue at a faster rate than those activated by the natural physiological processes. Due to this problem, a considerable effort has been directed toward developing FES systems based on closed-loop control. The movement is measured in real time with several types of sensors, and the stimulation pattern is modulated accordingly [1]. The dynamics of the lower limb is represented by a nonlinear second-order model, which considers the gravitational and inertial characteristics of the anatomical segment as well as the damping and stiffness properties of the knee joint.
In this paper, we present a Takagi–Sugeno nonlinear system with the aim of controlling the position of the leg of a paraplegic patient. The controller was designed in order to change the angle of the knee joint from to when electrical stimulation is applied in the quadriceps muscle.
The authors considered the leg mathematical model proposed by Ferrarin and Pedotti [1], with the parameter values given in [2, 3]. The parameters (viscous coefficient), (inertial moment), (time constant), and (static gain) of the shank-foot complex model have the nominal values given in [2, 3], but with a 20% tolerance range around these nominal values, that is, these values are in the range between 80% and 120% of their nominal values. The minimum and the maximum values of the nonlinear term are computed considering the angle variation from to , that is, . The range of values of , , , and are considered as norm-bounded uncertainties, whose analysis requires a lower number of linear matrix inequalities (LMIs), compared to polytopic uncertainty analysis, obtaining a lower computational cost. For the case studied in this manuscript, with two local models and four uncertain parameters, the control design methods that consider polytopic uncertainty analysis require the solution of a set of 49 LMIs, while, for the norm-bounded uncertainty analysis, only 4 LMIs are required. In this paper, the proposal for the knee position control design of paraplegic patients with functional electrical stimulation (FES) considers that the parameters of the mathematical model of the system are uncertain, whose uncertainties are bounded in norm. To the authors’ knowledge, the Takagi–Sugeno (T-S) fuzzy control considering norm-bounded uncertainties, applied to the knee joint movement of the paraplegic patient, was not published yet.
The simplest design technique to obtain a design model for nonlinear plants is its linearization at an interest point. However, this linearized design described is not adequate when the system operates far from the operation point. A possible solution for this problem is the nonlinear plant representation by T-S fuzzy models, whose idea consists on the description of the nonlinear system as a combination of a certain number of local linear models. So, the global model is obtained by the fuzzy combination of these local linear models.
The T-S fuzzy methodology can also be applied if the nonlinear model contains polytopic uncertainties, as in [4]. The representation of the uncertainties considers that the uncertain system structure varies according to a convex combination of some vertices that limit the polytope, where each vertex is described as a fuzzy combination of local linear models.
Fuzzy control theory is useful because the fuzzy systems can approximately represent real systems with a precision that can be specified by the designer. Furthermore, there are several types of models, suitable for different applications, since linguistic models for modeling a given system, even the T-S models, whose structure is suitable for control applications. The systematic procedure to design fuzzy control systems involves the fuzzy model construction for nonlinear systems [5].
The parallel distributed compensation (PDC) [6] in fuzzy regulator design can be used to stabilize nonlinear systems described by fuzzy models. The idea is to design a compensator for each fuzzy rule. For each rule, there exists an associated controller. The resulting global fuzzy regulator, which is nonlinear in general, is a fuzzy combination of each individual linear regulator. The PDC offers a procedure to design a regulator for each T-S fuzzy model, where each control rule is designed from the correspondent plant T-S model rule.
An important fact that motivates the use of the representation of a broad class of nonlinear plants by Takagi–Sugeno fuzzy models for designing suitable controllers is that they usually allow a design procedure, with PDC approach, based on LMIs. LMI-based designs can offer conditions for the stability of the equilibrium point, and also, it is possible to specify other performance indexes, such as decay rate, constraints of the input and/or output, and the minimization of the cost, even for plants with uncertain parameters. Furthermore, the procedure for finding a feasible solution, when there exists one, becomes a convex optimization problem, and there exist easy methods such as polynomial time algorithms for solving this class of problem. With the aforementioned relevant facts, many researchers are using Takagi–Sugeno fuzzy models for obtaining adequate controllers for solving complex nonlinear control problems.
Muscle is a highly complex nonlinear system [7], capable of producing the same output for a variety of inputs. A property exploited by the physiologically activated muscle is its effort to minimize fatigue [8]. Considering that when the quadriceps is electrically stimulated, its response is nonlinear, we used T-S fuzzy models in order to design a controller for the knee angle variation.
Ferrarin and Pedotti [1] showed that, for the conditions considered in their experiments, a simple one-pole transfer function was able to model the relationship between stimulus pulse width and active muscle torque. The nonlinear term is analyzed in a T-S fuzzy representation. So, the mathematical model of the functional electrical stimulation of the knee angle of the paraplegic patient is represented as a T-S fuzzy combination of two local models, considering the minimum and the maximum values of the nonlinear term .
In [9], the nonlinear term is considered as the uncertain parameter, and the parameters , , , and of the shank-foot complex model have definite values. However, these parameters assume different values each day, depending on the health conditions of the patient at each moment. Fatigue can also change these values. So, these parameters have unknown values, but these can be considered in determinate ranges, bounded by the minimum and the maximum values of each parameter [10].
According to Santos et al. and Gaino [9, 11], the application of FES on the quadriceps muscle, more particularly on the motor neurons of a person, causes an involuntary contraction of the muscle of the leg, that is, causes an action potential (AP). FES is also known by neuromuscular electrical stimulation (NMES) [11].
In order to obtain the muscle contraction, the amplitude (or intensity) and duration of the electrical stimulus must be inside specific bounds. Then, the AP is generated and propagates in both directions of the nerve fiber [9]. Complex mechanisms of electrochemical stimuli occur in the neuromuscular structure causing the process of excitation-contraction coupling responsible for the movement of the leg [11]. The modulation of the force, by the number of muscle fibers recruited, and the speed of fiber recruitment depend on several parameters. Some of these parameters include the proximity of the nerve fiber and the electrode, the electrode diameter, and the variation of the number of active states of the fibers by the variation of the amplitude or pulse duration [9, 11]. As can be seen in Figure 1, the degree of muscle activation () is a nonlinear function that depends on the duration of the stimulus .
By the knowledge of the authors, scientific studies on the application of T-S controllers to control the leg position of paraplegic patients are interesting, relevant, and challenging research topics. Recently, few articles have been published in this area, as, for instance, [12, 13]. In [11] and subsequent articles [2, 3, 14–17], published by our research group, some theories have been employed to control knee joint movement using neuromuscular electrical stimulation. In [18], an article was published describing control application in paraplegic patients.
The T-S fuzzy strategy for control of the knee joint angle of a paraplegic position using LMIs and considering polytopic uncertainties was studied in [19], but other papers involving T-S fuzzy control with polytopic uncertainties were not found in the literature. The T-S fuzzy control using LMIs and considering norm-bounded uncertainties was not found in the literature. In [20], LMI stabilization conditions considering polytopic and norm-bounded uncertainties were presented for fractional-order systems, but the method does not consider fuzzy models.
The list of the variables used in this manuscript is presented in Appendix.
2. General Takagi–Sugeno Fuzzy Representation
Certain classes of nonlinear systems can be exactly represented with T-S fuzzy models, using the method described in [5]. According to this construction method, the local models are obtained in function of the operation region.
For , the -th rule of the continuous-time T-S fuzzy model is described as
In the fuzzy model (1), , is the fuzzy set of rule , and are the premise variables. Let be the membership function of the fuzzy set , and define
Since , one has, for ,
The resulting fuzzy model is the weighted mean of the local models, given bywhere
2.1. Regulators with Takagi–Sugeno Fuzzy Models
The PDC [6] offers a procedure to design a regulator for the T-S fuzzy model, where each control rule is designed from the correspondent plant T-S model rule. The designed fuzzy regulator shares the local controllers, each one given by
The fuzzy global regulator is given bywhere .
So, the closed-loop system is composed as follows:
3. Norm-Bounded Uncertainties
According to [21], a control system with norm-bounded uncertainties is described bywith , , and . Matrices , , and are composed by the nominal values of the uncertain parameters.
The uncertain matrices and can be represented aswhere
So, matrix is given bywithfor . Then, considering (13), one obtains equations that are given byfor and .
Similarly, matrix is given bywherefor and . Then, considering (16), one obtains equations that are given byfor , , and .
4. Knee Joint Model
The mathematical model of the leg employed in this manuscript was proposed in [1]. This model relates the applied pulse width to the torque generated on the knee joint. In the modeling [1], the leg was considered as an open kinematic system composed of two stiff segments: the thigh and the shin-foot complex, as shown in Figure 2. From [1], the equilibrium equation around the knee joint iswhere is the inertial moment of the shin-foot complex is the knee angle between the shin and the thigh on the sagittal plane is the knee angular velocity is the mass of the shin-foot complex is the gravitational acceleration is the distance between the knee and the shin-foot complex mass center is the viscous friction coefficient is the torque due to the stiffness component is the knee active torque generated by electrical stimulation
The stiffness moment is defined aswhere and are the coefficients of the exponential terms and is the elastic rest angle of the knee. In [1], it was observed that the torque applied to the muscle () and the electrical stimulation pulse width () can be adequately related by the transfer function described in (20), where and are positive constants obtained from parametric identification, when a stimulus is applied on the quadriceps and angular variation is read, as mentioned in [1]:
So, from [1], the knee joint mathematical model is given bywithwhere , , and are the knee angle, the knee angular velocity (), and the active torque at the operation point. The input is the pulse width of the electrical stimulation.
The nonlinear term is given bywhere
4.1. Local Fuzzy Models of the Knee Joint considering Norm-Bounded Uncertainties
At the system described by (21), the uncertain parameter given in (24) belongs to the interval . Then, the number of rules is , that is, the T-S fuzzy system has two local models. So, appropriate functions and are defined as
The nonlinear function is exactly represented by the convex combination of these two local models: [2, 3, 5, 11, 13]. Note that, from (26), , , and .
So, considering two local models, matrices and arewith described in (24) and (25), which belongs to the interval
Functions and are given by (17), and the matrices that describe the local models arewherewith
So,withknowing that
So,
Replacing (35)–(38) in and , one haswhere
So, matrices and are decomposed according to (10), where
Following (13), the elements of are functions of , , , and , given byfor , where are the elements of matrix in (43).
From (16), the elements of are also the function of , , , and , given byfor (), where are the elements of matrix in (44).
According to (14) and (17), the elements of , , and are solutions of the following equations:where , , , and are given in (45)–(48).
So, the solution of the set of equations (52) is
For simplicity, one can choose . Therefore, matrices , and arewith , and given in (45)–(48).
Similarly, replacing (39)–(42) in and , one haswhere
So, doing the same previous calculations, matrices and are decomposed according to (10), where
For simplicity, one can choose . Therefore, matrices and arewith , , , and given in (57)–(60).
5. Numerical Model of the Paraplegic Patient
The mathematical model of the system, which relates the pulse width applied for the muscle to the torque generated at the knee joint, is presented in Section 4. The parameters of the system, presented in [2, 3], were obtained experimentally. The parameters related to the shank-foot complex and their numerical values are given in Table 1.
The values of these parameters were experimentally measured, for a 45-year-old paraplegic patient [2, 3]. However, several factors, such as changes of temperature, fatigue, and spasm, cause physiological changes on musculature that must be considered on process control. Furthermore, for other people, the physiological characteristics may be completely different since these characteristics depend on several factors, for instance, the age, weight, physical activities, and health conditions. So, calibration is needed before the beginning of the tests. The adjustments to perform a desired movement are made after the identification of each patient at a specific current and frequency.
The nonlinear term is described in (24) and (25). Considering the parameters of Table 1, knowing that = 9.8 m/s2 and taking the operation point angle as rad, these values are replaced in (25), obtaining N·m. For this operation point, the term belongs to the interval , where
The new input of the system, , is defined from the system input, , and is known as the unreferenced pulse width [11, 22–24]. It is given by
Since the input is the pulse width which is applied on the skin of the patient, its value must be positive, that is, . So,
5.1. Norm-Bounded Uncertainties
From the parameters presented in Table 1, a ±20% variation is considered on the values of J, B, τ, and G, that is,
Given the minimum and the maximum values of , one has the system with two local models, each one described by (8), whereand matrices and are described by (10), wherewith , for , and matrices and are described by (10), wherewith , for .
6. System Control with Norm-Bounded Uncertainties
For the T-S fuzzy system given in (4), where each local model is described in (8), the objective is to obtain a fuzzy control law, given in (7), such that the controlled system is stable.
Theorem 1, given in [25], gives a sufficient condition for the stabilization of system (4) by the fuzzy control law (7), with the uncertain matrices described in (6).
Theorem 1 (see [25]). Consider continuous-time T-S fuzzy system (4), with q local models, where each local model is described as in (9) and (10), that is, each matrix, and , , is decomposed as and , respectively, with , and .
So, continuous-time T-S fuzzy system (4) is asymptotically stabilizable via the T-S fuzzy model-based state-feedback controller (7) if there exist a symmetric positive definite matrix, some matrices, and some scalars such that the following LMIs are satisfied:where, and , where denotes the transposed elements in the symmetric positions. From the solution of the aforementioned LMIs, output feedback matrices Ki are obtained from .
Multiplying (70), at left and at right, by , the following LMI is obtained:where andNow, multiplying (71), at left and at right, by , the following LMI is obtained:where andSo, it is observed from (74), that T-S fuzzy system (4), with the local models described in (9) and (10) and the control law (7), is asymptotically stable if there exist matrices and scalars , for such that LMIs (73), (75) and hold.
For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections 5 and 6, the number of rules is . So, to achieve the fuzzy control law (7) that stabilizes system (4), the problem consists of finding matrices and scalars such that the following LMIs hold:whereFrom the solution of the aforementioned LMIs, matrices and are given by
Remark 1. For the case studied in this manuscript, with two local models and four uncertain parameters, the control design methods that consider polytopic uncertainty analysis required the solution of a set of 49 LMIs, while for the norm-bounded uncertainty analysis, only 4 LMIs are sufficient. A smaller number of LMIs requires a lower computational cost so that the LMI set is less conservative, and the solution is obtained more fastly.
6.1. Decay Rate with Norm-Bounded Uncertainties
Sometimes, only stability is not sufficient to get a suitable performance for a control system. Frequently, the transient response must also be specified. Given a linear time-invariant system, , according to [26], the decay rate is defined as the maximum value of the real constant such thatholds for all trajectories .
The condition in a system for all trajectories of for is equivalent to
Then, replacing matrix in (73) and (75), the LMIs that guarantee for T-S fuzzy system (4) with the local models described in (9) and (10), the fuzzy control law (7), and a decay rate greater than arewhere , and, where denotes the transposed elements in the symmetric positions. From the solution of the aforementioned LMIs, output feedback matrices are obtained from .
For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections 5 and 6, the number of rules is . So, to achieve a control law (7) that guarantees for system (4), where each local model is described in (9) and (10) and decay rate is greater than γ, the goal is to find matrices and scalars such that the following LMIs hold:where
From the solution of the aforementioned LMIs, matrices are given by
7. System Control with Input Constraint
Consider the input control given in (4) and an initial state . In various situations, it is necessary to specify a bound for the input control to avoid a too large magnitude of this input. Consider the following restriction on the input in [24]:for all
This condition is guaranteed by adding new LMIs in the LMI set that determines stability or decay rate. Input constraint (92) is guaranteed by the LMIs:for an initial state The LMIs that guarantee the input constraint need to be added to the LMI set that determines stability or decay rate.
For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections 5 and 6, the number of rules is . So, the LMIs that guarantee the input constraint arefor an initial state . LMIs (94)–(96) need to be added to the LMI set that determines stability or decay rate, given in (77)–(80) or (87)–(89).
For the FES, the control signal is the electric pulse width applied to the muscle that is given bywhose value must, naturally, be positive. So, the input needs to follow the condition
To follow this condition, input constraint (92) is added for the input, where
For one has The initial state is
8. Numerical Results and Simulation
After obtaining the mathematical model of FES for the knee joint of a paraplegic patient, considering norm-bounded uncertainties, as described in Sections 4 and 5, the LMIs presented in Section 6 were solved using the Control Toolbox of MATLAB, version 2012a [27].
8.1. Stability
Although the plant is already stable, the first action was to solve LMIs (77)–(80) to obtain a control law (7) that stabilizes system (4). The solution of these LMIs is
From the solution above, matrices of the control law (7) are
Table 2 shows the eigenvalues of , with the aforementioned matrices , related to the nominal system, when only stability is specified. In Figure 3, the simulation result is shown for the system with the nominal values of , from the initial state given in (100), with α described in (17).
(a)
(b)
(c)
(d)
The convergence of the input is very fast, reaching zero at a short time. Although some values seem to be negative, the real value is zero. A zoom on the time, from 0 to 20 s, in Figure 4, shows this fast convergence.
In Figures 5–7, the simulation result is presented for , respectively, considering all combinations of the minimum or maximum values of the parameters , from the initial state given in (100), with described in (17). Each figure contains 16 curves, each one representing one combination of the extreme (minimum or maximum) values of , as described in Table 3. The results show that the controlled system is stable for all possible values of these parameters.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
8.2. Decay Rate
To obtain a faster transient response, decay rate was specified. For this purpose, LMIs (87)–(89) were solved. The solution of these LMIs is
From the aforementioned solution, matrices of the control law (7) are
Table 4 shows the eigenvalues of , with the aforementioned matrices , related to the nominal system, when decay rate is specified. In Figure 8, the simulation result is shown for the system with the nominal values of , from the initial state given in (100), with α described in (17).
(a)
(b)
(c)
(d)
The convergence of the input is very fast, reaching zero at a short time. Although some values seem to be negative, the real value is zero. A zoom on the time, from 0 to 20s, in Figure 9, shows this fast convergence.
In Figures 10–12, the simulation result is presented for , respectively, considering all combinations of the minimum or maximum values of the parameters , from the initial state given in (100), with described in (17). Each figure contains 16 curves, each one representing one combination of the extreme (minimum or maximum) values of , as described in Table 3. The results show that the controlled system holds the decay rate specification for all possible values of these parameters.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
8.3. Input Constraint
Since the pulse width (63) must be positive, input constraint (92), with , was considered in the problem. Then, LMIs (94)–(96) were added to LMIs (77)–(80).
When stability and input constraint were specified, the Matlab Control Toolbox did not find the solution for LMIs (77)–(80) and (94)–(96), with small values of . Specification of a very small bound on the control input can make the problem too conservative, and it is not easy to find a solution. For , the solution of these LMIs is
From the aforementioned solution, matrices K1 and K2 of the control law (7) are
Table 5 shows the eigenvalues of , with the aforementioned matrices , related to the nominal system, when input constraint (92), with , is specified. In Figure 13, the simulation result is shown for the system with the nominal values of , from the initial state given in (100), with described in (17).
(a)
(b)
(c)
(d)
The convergence of the input is very fast, reaching zero at a short time. Although some values seem to be negative, the real value is zero. A zoom on time, from 0 to 60 s, in Figure 14, shows this fast convergence.
In Figures 15–17, the simulation result is presented for , respectively, considering all combinations of the minimum or maximum values of the parameters , from the initial state given in (100), with described in (17). Each figure contains 16 curves, each one representing one combination of the extreme (minimum or maximum) values , as described in Table 3.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
When decay rate and an input constraint were specified, no solution for LMIs (87)–(89) and (94)–(96) was found. This fact may have occurred because the specification of decay rate and an input constraint has increased the conservativeness of the problem.
8.4. Discussion of the Results
Although the plant is already stable, the first purpose is to get a control law that guarantees the stability of the system, for any exact value of the uncertain parameters, within the considered set of values. As seen in Figure 3, for the nominal plant, and in Table 2 and Figures 5–7, considering the minimum or maximum value of each uncertain parameter, the controlled system is stable, and its state variables converge fastly to the equilibrium point. The control law converges to zero at a short time, as seen in Figures 3 and 4.
To obtain a faster transient response, decay rate was determined. As seen in Figure 8, for the nominal plant, and in Table 4 and Figures 10–12, considering the minimum or maximum value of each uncertain parameter, the convergence of the state variables is faster compared with the first case. In Tables 2 and 4, one can note that, for the second case, the closed-loop poles of the system are further from the imaginary axis than for the first case.
However, as observed in Figures 8 and 9, the control law assumes too high values, so this control signal is impracticable for the real system.
To reduce the magnitude of u(t), a constraint on the input, described in (92), was specified. However, no solution was found for small values of μ0. Forμ0 = 200, the control law stabilizes the system, but the convergence of the state variables is slower than the previous cases, as seen in Figure 13, for the nominal plant, and in Table 5 and Figures 15–17, considering the minimum or maximum value of each uncertain parameter. The control law is also too high, but it converges fastly to zero, as observed in Figures 13 and 14. No solution was found when decay rate and an input constraint were specified.
9. Conclusion
In this manuscript, a mathematical model for the knee position control for a paraplegic patient, following the description in [1], was obtained using the parameters given in [2, 3] and considering norm-bounded uncertainties. The uncertain model is a T-S fuzzy combination of linear local models. After finding the T-S fuzzy model of the uncertain nonlinear system, a T-S fuzzy state-feedback control has been designed for this system using LMIs. The initial objective, although the plant is already stable, is to obtain a stable system. Then, to obtain a faster transient response, decay rate was specified. Finally, an input constraint was specified to guarantee that the pulse width is always positive. The simulation results show the efficiency of the control.
Appendix
The variables used in this manuscript are described in Tables 6–9. The variables used in the T-S fuzzy representation are presented in Table 6. Table 7 lists the variables used to describe the norm-bounded uncertainties. Note that, on T-S fuzzy description, a new index can be added to these variables, corresponding to the fuzzy rule that is being used. The knee joint parameters are given in Table 8, and the variables used in the LMI control strategy are presented in Table 9.
Data Availability
The system parameter data used to support the findings of this study are published in [2, 3] and are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the Coordenação de Aperfeiçoamentode Pessoal de Nível Superior-Brasil (CAPES)- Finance Code 001, Conselho Nacional de Desenvolvimento Científico e Tecnológico -(CNPq), and Pró-Reitoria de Pesquisa e Pós-Graduação da Universidade Estadual de Londrina (PROPPG-UEL), from Brazil.