# Robust T-S Fuzzy Control of Electrostimulation for Paraplegic Patients considering Norm-Bounded Uncertainties

Authors: Márcio Roberto Covacic; Marcelo Carvalho Minhoto Teixeira; Aparecido Augusto de Carvalho; Rodrigo Cardim; Edvaldo Assunção; Marcelo Augusto Assunção Sanches; Henrique Shuiti Fujimoto; Maxwell Simões Mineo; Anderson Ross Biazeto; Ruberlei Gaino
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 from0° to 30° 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 B (viscous coefficient), J (inertial moment), τ (time constant), and G (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 f˜21x1 are computed considering the angle variation from 0° to 60°, that is, −30°≤x1≤30°. The range of values of B, J, τ, and G 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 f˜21x1 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 f˜21x1.In [9], the nonlinear term f˜21x1 is considered as the uncertain parameter, and the parameters B, J, τ, and G 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 d.Figure 1
Fiber recruitment curve (the black circles are activated fibers) [7].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.Fori=1,2…,q, the i-th rule of the continuous-time T-S fuzzy model is described as(1)rulei:ifz1tisℳ1iand…zptisℳpi,thenx˙t=Aixt+Biut,yt=Cxt.In the fuzzy model (1), ℳji,j=1,2,…,p, is the fuzzy set j of rule i, and z1t,…,zpt are the premise variables. Let μjizjt be the membership function of the fuzzy set ℳji, and define(2)wizt=∏j=1pμjizjt,zt=z1,…,zpT.Sinceμjizjt≤0, one has, for i=1,2,…,q,(3)wizt≥0,i=1,…,q.The resulting fuzzy model is the weighted mean of the local models, given by(4)x˙t=∑i=1qwiztAixt+Biut∑i=1qwizt=∑i=1qαiztAixt+Biut,=Aαzt+Bαutyt=Cxt,where(5)α=α1,…,αqT,αizt=wizt∑i=1qwizt,αizt≥0,i=1,…,q,∑i=1qαizt=1.


### 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(6)rulei:ifz1tisℳ1iand…zptisℳpi,thenut=−Kixt.The fuzzy global regulator is given by(7)ut=−∑i=1qwiztKixt∑i=1qwizt=−∑i=1qαiztKixt=−Kαxt,where α=α1,…,αqT.So, the closed-loop system is composed as follows:(8)x˙t=∑i=1q∑j=1qαiztαjztAi−BiKjxt.


## 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(6)rulei:ifz1tisℳ1iand…zptisℳpi,thenut=−Kixt.The fuzzy global regulator is given by(7)ut=−∑i=1qwiztKixt∑i=1qwizt=−∑i=1qαiztKixt=−Kαxt,where α=α1,…,αqT.So, the closed-loop system is composed as follows:(8)x˙t=∑i=1q∑j=1qαiztαjztAi−BiKjxt.


## 3. Norm-Bounded Uncertainties

According to [21], a control system with norm-bounded uncertainties is described by(9)x˙t=An+δAxt+Bn+δBut,yt=Cxt,with xt∈ℝn, ut∈ℝm, and yt∈ℝp. Matrices An, Bn, and Cn are composed by the nominal values of the uncertain parameters.The uncertain matricesδA and δB can be represented as(10)δA=LΔRA,δB=LΔRB,where(11)L=ℓ11ℓ12⋯ℓ1νℓ21ℓ22⋯ℓ2ν⋮⋮⋱⋮ℓi1ℓi2⋯ℓiν⋮⋮⋱⋮ℓn1ℓn2⋯ℓnν,Δ=diagδ1,δ2,…,δν,Δ=maxδi≤1,RA=ra11ra12⋯ra1j⋯ra1nra21ra22⋯ra2j⋯ra2n⋮⋮⋱⋮⋱⋮raν1raν2⋯raνj⋯raνn,RB=rb11rb12⋯rb1j⋯rb1mrb21rb22⋯rb2j⋯rb2m⋮⋮⋱⋮⋱⋮rbν1rbν2⋯rbνj⋯rbνm.So, matrixδA is given by(12)δA=LΔRA=ϕ11ϕ12⋯ϕ1nϕ21ϕ22⋯ϕ2n⋮⋮⋱⋮ϕn1ϕn2⋯ϕnn,with(13)ϕij=∑k=1νℓikδkrakj,for i,j=1,2,…,n. Then, considering (13), one obtains n2ν equations that are given by(14)ℓikrakj=∂ϕij∂δk,for i,j=1,2,…,n and k=1,2,…,ν.Similarly, matrixδB is given by(15)δB=LΔRB=β11β12⋯β1mβ21β22⋯β2m⋮⋮⋱⋮βn1βn2⋯βnm,where(16)βij=∑k=1νℓikδkrbkj,for i=1,2,…,n and j=1,2,…,m. Then, considering (16), one obtains nmν equations that are given by(17)ℓikrbkj=∂βij∂δk,for i=1,2,…,n, j=1,2,…,m, and k=1,2,…,ν.


## 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 is(18)Jθ¨v=−mglsinθv−Ms−Bθ˙+Ma,whereJ is the inertial moment of the shin-foot complexθv is the knee angle between the shin and the thigh on the sagittal planeθ=θv+90°θ. is the knee angular velocitym is the mass of the shin-foot complexg is the gravitational accelerationl is the distance between the knee and the shin-foot complex mass centerB is the viscous friction coefficientMs is the torque due to the stiffness componentMa is the knee active torque generated by electrical stimulationFigure 2
Shin-foot complex and its parametersθ, θv, and Ma.The stiffness moment is defined as(19)Ms=λe−Eθθ−ω,where λ and E 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 (Ma) and the electrical stimulation pulse width (Pu) can be adequately related by the transfer function described in (20), where G 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]:(20)Hs=MasPus=G1+sτ.So, from [1], the knee joint mathematical model is given by(21)x˙t=Aαxt+Bαut,yt=Cxt,with(22)xt=ΔθvtΔθ˙vtΔMat=θvt−θv0θ˙vt−θ˙v0Mat−Ma0,ut=PNt,yt=Δθvt,where θv0, θ.v0, and Ma0 are the knee angle, the knee angular velocity (θ˙v0=0), and the active torque at the operation point. The input PNt is the pulse width of the electrical stimulation.(23)Aα=010f˜21x1−BJ1J00−1τ,Bα=00Gτ,C=100.The nonlinear termf˜21x1 is given by(24)f˜21x1=1Jx1mglsinx1t+θv0−1Jx1×λexp−Ex1t+θv0+π2x1t+θv0+π2−w+Ma0,where(25)Ma0=mglsinθv0+λexp−Eθv0+π2θv0+π2−w.


### 4.1. Local Fuzzy Models of the Knee Joint considering Norm-Bounded Uncertainties

At the system described by (21), the uncertain parameter f˜21x1 given in (24) belongs to the interval f˜21min≤f˜21x1≤f˜21max. Then, the number of rules is q=2, that is, the T-S fuzzy system has two local models. So, appropriate functions α1 and α2 are defined as(26)α1=f˜21x1−f˜21maxf˜21min−f˜21max,α2=f˜21x1−f˜21minf˜21max−f˜21min=1−α1.The nonlinear functionf˜21x1 is exactly represented by the convex combination of these two local models: f˜21x1=α1f˜21min+α2f˜21max [2, 3, 5, 11, 13]. Note that, from (26), α1≥0, α2≥0, and α1+α2=1.So, considering two local models, matricesAα and Bα are(27)Aα=α1A1n+δA1+α2A2n+δA2,Bα=α1B1n+δB1+α2B2n+δB2,α1≥0,α2≥0,α1+α2=1,A1n+δA1=010f˜21min−BJ1J00−1τ,A2n+δA2=010f˜21max−BJ1J00−1τ,B1n+δB1=B2n+δB2=00Gτ,C=100,with f˜21x1 described in (24) and (25), which belongs to the interval(28)f˜21min≤f˜21≤f˜21max.Functionsα1 and α2 are given by (17), and the matrices that describe the local models are(29)A1n+δA1=010f˜21mina122a12300a133,A2n+δA2=010f˜21maxa222a22300a233,B1n+δB1=00b131,B2n+δB2=00b231,C=100,where(30)a122=a222=−BJ,a123=a223=1J,a133=a233=−1τ,b131=b231=Gτ,with(31)0<Bmin≤B≤Bmax,0<Jmin≤J≤Jmax,0<τmin≤τ≤τmax,0<Gmin≤G≤Gmax.So,(32)A1n=010f˜21mina122na123n00a133n,A2n=010f˜21maxa222na223n00a233n,B1n=00b131n,B2n=00b231n,δA1=0000δa122δa12300δa133,δA2=0000δa222δa22300δa233,δB1=00δb131,δB2=00δb231,with(33)a122n=a122max+a122min2,a123n=a123max+a123min2,a133n=a133max+a133min2,b131n=b131max+b131min2,a222n=a222max+a222min2,a223n=a223max+a223maxmin2,a233n=a233max+a233min2,b231n=b231max+b231min2,a122min−a122max2≤δa122≤a122max−a122min2,a123min−a123max2≤δa123≤a123max−a123min2,a133min−a133max2≤δa133≤a133max−a133min2,b131min−b131max2≤δb131≤b131max−b131min2,a222min−a222max2≤δa222≤a222max−a222min2,a223min−a223max2≤δa223≤a123max−a223min2,a233min−a233max2≤δa233≤a233max−a233min2,b231min−b231max2≤δb231≤b231max−b231min2,knowing that(34)a122min≤a122≤a122max,a123min≤a123≤a123max,a133min≤a133≤a133max,b131min≤b131≤b131max,a222min≤a222≤a222max,a223min≤a223≤a223max,a233min≤a233≤a233max,b231min≤b231≤b231max.So,(35)δa122=a122max−a122min2δ11,−1≤δ11≤1,(36)δa123=a123max−a123min2δ12,−1≤δ12≤1,(37)δa133=a133max−a133min2δ13,−1≤δ13≤1,(38)δb131=b131max−b131min2δ14,−1≤δ14≤1,(39)δa222=a222max−a222min2δ21,−1≤δ21≤1,(40)δa223=a223max−a223min2δ22,−1≤δ22≤1,(41)δa233=a233max−a233min2δ23,−1≤δ23≤1,(42)δb231=b231max−b231min2δ24,−1≤δ24≤1.Replacing (35)–(38) in δA1 and δB1, one has(43)δA1=0000ρ122δ11ρ123δ1200ρ133δ13,(44)δB1=00ρ134δ14,where(45)ρ122=a122max−a122min2,(46)ρ123=a123max−a123min2,(47)ρ133=a133max−a133min2,(48)ρ134=b131max−b131min2.So, matricesδA1 and δB1 are decomposed according to (10), where(49)L1=ℓ111ℓ112ℓ113ℓ114ℓ121ℓ122ℓ123ℓ124ℓ131ℓ132ℓ133ℓ134,Δ1=δ110000δ120000δ130000δ14,RA1=ra111ra112ra113ra121ra122ra123ra131ra132ra133ra141ra142ra143,RB1=rb111rb121rb131rb141.Following (13), the elements of δA1 are functions of δ11, δ12, δ13, and δ14, given by(50)ϕ1ij=ℓ1i1δ11ra11j+ℓ1i2δ12ra12j+ℓ1i3δ13r1a3j+ℓ1i4δ14ra14j,for i,j=1,2,3, where ϕ1ij are the elements of matrix δA1 in (43).From (16), the elements of δB1 are also the function of δ11, δ12, δ13, and δ14, given by(51)β1i1=ℓ1i1δ11rb111+ℓ1i2δ12rb121+ℓ1i3δ13rb131+ℓ1i4δ14rb141,for i=1,2,3 (j=1), where β1i1 are the elements of matrix δB1 in (44).According to (14) and (17), the elements of L1, RA1, and RB1 are solutions of the following equations:(52)ℓ111ra111=0,ℓ112ra121=0,ℓ113ra131=0,ℓ114ra141=0,ℓ111ra112=0,ℓ112ra122=0,ℓ113ra132=0,ℓ114ra142=0,ℓ111ra113=0,ℓ112ra123=0,ℓ113ra133=0,ℓ114ra143=0,ℓ111rb111=0,ℓ112rb121=0,ℓ113rb131=0,ℓ114rb141=0,ℓ121ra111=0,ℓ122ra121=0,ℓ123ra131=0,ℓ124ra141=0,ℓ121ra112=ρ122,ℓ122ra122=0,ℓ123ra132=0,ℓ124ra142=0,ℓ121ra113=0,ℓ122ra123=ρ123,ℓ123ra133=0,ℓ124ra143=0,ℓ121rb111=0,ℓ122rb121=0,ℓ123rb131=0,ℓ124rb141=0,ℓ131ra111=0,ℓ132ra121=0,ℓ133ra131=0,ℓ134ra141=0,ℓ131ra112=0,ℓ132ra122=0,ℓ133ra132=0,ℓ134ra142=0,ℓ131ra113=0,ℓ132ra123=0,ℓ133ra133=ρ133,ℓ134ra143=0,ℓ131rb111=0,ℓ132rb121=0,ℓ133rb131=0,ℓ134rb141=ρ134,where ρ122, ρ123, ρ133, and ρ134 are given in (45)–(48).So, the solution of the set of equations (52) is(53)L1=0000ρ122/ra112ρ123/ra1230000ρ133/ra133ρ134/ra131,RA1=0ra112000ra12300ra133000,RB1=000ra131,For simplicity, one can choosera112=ra123=ra133=rb131=1. Therefore, matrices L1, RA1 and RB1 are(54)L1=0000ρ122ρ1230000ρ133ρ134,RA1=010001001000,RB1=0001,with ρ122,ρ123,ρ133, and ρ134 given in (45)–(48).Similarly, replacing (39)–(42) in δA2 and δB2, one has(55)δA2=0000ρ222δ21ρ223δ2200ρ233δ23,(56)δB2=00ρ234δ24,where(57)ρ222=α222max−α222min2,(58)ρ223=α223max−α223min2,(59)ρ233=α233max−α233min2,(60)ρ234=α231max−α231min2.So, doing the same previous calculations, matricesδA2 and δB2 are decomposed according to (10), where(61)L2=0000ρ222/ra212ρ223/ra2230000ρ233/ra233ρ234/rb231,Δ2=δ210000δ220000δ230000δ24,RA2=0ra212000ra22300ra233000,RB2=000rb231.For simplicity, one can choosera212=ra223=ra233=rb231=1. Therefore, matrices L2,RA2 and RB2 are(62)L2=0000ρ222ρ2230000ρ233ρ234,RA2=010001001000,RB2=0001,with ρ222, ρ223, ρ233, and ρ234 given in (57)–(60).


## 4.1. Local Fuzzy Models of the Knee Joint considering Norm-Bounded Uncertainties

At the system described by (21), the uncertain parameter f˜21x1 given in (24) belongs to the interval f˜21min≤f˜21x1≤f˜21max. Then, the number of rules is q=2, that is, the T-S fuzzy system has two local models. So, appropriate functions α1 and α2 are defined as(26)α1=f˜21x1−f˜21maxf˜21min−f˜21max,α2=f˜21x1−f˜21minf˜21max−f˜21min=1−α1.The nonlinear functionf˜21x1 is exactly represented by the convex combination of these two local models: f˜21x1=α1f˜21min+α2f˜21max [2, 3, 5, 11, 13]. Note that, from (26), α1≥0, α2≥0, and α1+α2=1.So, considering two local models, matricesAα and Bα are(27)Aα=α1A1n+δA1+α2A2n+δA2,Bα=α1B1n+δB1+α2B2n+δB2,α1≥0,α2≥0,α1+α2=1,A1n+δA1=010f˜21min−BJ1J00−1τ,A2n+δA2=010f˜21max−BJ1J00−1τ,B1n+δB1=B2n+δB2=00Gτ,C=100,with f˜21x1 described in (24) and (25), which belongs to the interval(28)f˜21min≤f˜21≤f˜21max.Functionsα1 and α2 are given by (17), and the matrices that describe the local models are(29)A1n+δA1=010f˜21mina122a12300a133,A2n+δA2=010f˜21maxa222a22300a233,B1n+δB1=00b131,B2n+δB2=00b231,C=100,where(30)a122=a222=−BJ,a123=a223=1J,a133=a233=−1τ,b131=b231=Gτ,with(31)0<Bmin≤B≤Bmax,0<Jmin≤J≤Jmax,0<τmin≤τ≤τmax,0<Gmin≤G≤Gmax.So,(32)A1n=010f˜21mina122na123n00a133n,A2n=010f˜21maxa222na223n00a233n,B1n=00b131n,B2n=00b231n,δA1=0000δa122δa12300δa133,δA2=0000δa222δa22300δa233,δB1=00δb131,δB2=00δb231,with(33)a122n=a122max+a122min2,a123n=a123max+a123min2,a133n=a133max+a133min2,b131n=b131max+b131min2,a222n=a222max+a222min2,a223n=a223max+a223maxmin2,a233n=a233max+a233min2,b231n=b231max+b231min2,a122min−a122max2≤δa122≤a122max−a122min2,a123min−a123max2≤δa123≤a123max−a123min2,a133min−a133max2≤δa133≤a133max−a133min2,b131min−b131max2≤δb131≤b131max−b131min2,a222min−a222max2≤δa222≤a222max−a222min2,a223min−a223max2≤δa223≤a123max−a223min2,a233min−a233max2≤δa233≤a233max−a233min2,b231min−b231max2≤δb231≤b231max−b231min2,knowing that(34)a122min≤a122≤a122max,a123min≤a123≤a123max,a133min≤a133≤a133max,b131min≤b131≤b131max,a222min≤a222≤a222max,a223min≤a223≤a223max,a233min≤a233≤a233max,b231min≤b231≤b231max.So,(35)δa122=a122max−a122min2δ11,−1≤δ11≤1,(36)δa123=a123max−a123min2δ12,−1≤δ12≤1,(37)δa133=a133max−a133min2δ13,−1≤δ13≤1,(38)δb131=b131max−b131min2δ14,−1≤δ14≤1,(39)δa222=a222max−a222min2δ21,−1≤δ21≤1,(40)δa223=a223max−a223min2δ22,−1≤δ22≤1,(41)δa233=a233max−a233min2δ23,−1≤δ23≤1,(42)δb231=b231max−b231min2δ24,−1≤δ24≤1.Replacing (35)–(38) in δA1 and δB1, one has(43)δA1=0000ρ122δ11ρ123δ1200ρ133δ13,(44)δB1=00ρ134δ14,where(45)ρ122=a122max−a122min2,(46)ρ123=a123max−a123min2,(47)ρ133=a133max−a133min2,(48)ρ134=b131max−b131min2.So, matricesδA1 and δB1 are decomposed according to (10), where(49)L1=ℓ111ℓ112ℓ113ℓ114ℓ121ℓ122ℓ123ℓ124ℓ131ℓ132ℓ133ℓ134,Δ1=δ110000δ120000δ130000δ14,RA1=ra111ra112ra113ra121ra122ra123ra131ra132ra133ra141ra142ra143,RB1=rb111rb121rb131rb141.Following (13), the elements of δA1 are functions of δ11, δ12, δ13, and δ14, given by(50)ϕ1ij=ℓ1i1δ11ra11j+ℓ1i2δ12ra12j+ℓ1i3δ13r1a3j+ℓ1i4δ14ra14j,for i,j=1,2,3, where ϕ1ij are the elements of matrix δA1 in (43).From (16), the elements of δB1 are also the function of δ11, δ12, δ13, and δ14, given by(51)β1i1=ℓ1i1δ11rb111+ℓ1i2δ12rb121+ℓ1i3δ13rb131+ℓ1i4δ14rb141,for i=1,2,3 (j=1), where β1i1 are the elements of matrix δB1 in (44).According to (14) and (17), the elements of L1, RA1, and RB1 are solutions of the following equations:(52)ℓ111ra111=0,ℓ112ra121=0,ℓ113ra131=0,ℓ114ra141=0,ℓ111ra112=0,ℓ112ra122=0,ℓ113ra132=0,ℓ114ra142=0,ℓ111ra113=0,ℓ112ra123=0,ℓ113ra133=0,ℓ114ra143=0,ℓ111rb111=0,ℓ112rb121=0,ℓ113rb131=0,ℓ114rb141=0,ℓ121ra111=0,ℓ122ra121=0,ℓ123ra131=0,ℓ124ra141=0,ℓ121ra112=ρ122,ℓ122ra122=0,ℓ123ra132=0,ℓ124ra142=0,ℓ121ra113=0,ℓ122ra123=ρ123,ℓ123ra133=0,ℓ124ra143=0,ℓ121rb111=0,ℓ122rb121=0,ℓ123rb131=0,ℓ124rb141=0,ℓ131ra111=0,ℓ132ra121=0,ℓ133ra131=0,ℓ134ra141=0,ℓ131ra112=0,ℓ132ra122=0,ℓ133ra132=0,ℓ134ra142=0,ℓ131ra113=0,ℓ132ra123=0,ℓ133ra133=ρ133,ℓ134ra143=0,ℓ131rb111=0,ℓ132rb121=0,ℓ133rb131=0,ℓ134rb141=ρ134,where ρ122, ρ123, ρ133, and ρ134 are given in (45)–(48).So, the solution of the set of equations (52) is(53)L1=0000ρ122/ra112ρ123/ra1230000ρ133/ra133ρ134/ra131,RA1=0ra112000ra12300ra133000,RB1=000ra131,For simplicity, one can choosera112=ra123=ra133=rb131=1. Therefore, matrices L1, RA1 and RB1 are(54)L1=0000ρ122ρ1230000ρ133ρ134,RA1=010001001000,RB1=0001,with ρ122,ρ123,ρ133, and ρ134 given in (45)–(48).Similarly, replacing (39)–(42) in δA2 and δB2, one has(55)δA2=0000ρ222δ21ρ223δ2200ρ233δ23,(56)δB2=00ρ234δ24,where(57)ρ222=α222max−α222min2,(58)ρ223=α223max−α223min2,(59)ρ233=α233max−α233min2,(60)ρ234=α231max−α231min2.So, doing the same previous calculations, matricesδA2 and δB2 are decomposed according to (10), where(61)L2=0000ρ222/ra212ρ223/ra2230000ρ233/ra233ρ234/rb231,Δ2=δ210000δ220000δ230000δ24,RA2=0ra212000ra22300ra233000,RB2=000rb231.For simplicity, one can choosera212=ra223=ra233=rb231=1. Therefore, matrices L2,RA2 and RB2 are(62)L2=0000ρ222ρ2230000ρ233ρ234,RA2=010001001000,RB2=0001,with ρ222, ρ223, ρ233, and ρ234 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 Section4. 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.Table 1
Parameter values obtained experimentally in [2] for a 45-year-old paraplegic patient (last column of Table 2 at page 127 of [2])
ParameterValueUnitJInertial moment of the shank-foot complex0.2881kgmMass of the shank-foot complex3.391kglDistance between the knee and the center of the mass of the shank-foot complex0.245mBViscous coefficient0.2756N·m·s/radλCoefficient of the exponential term1.3193N·m/radECoefficient of the exponential term0.0394ωResting elastic knee angle2.3185radτTransfer function coefficient (time constant of the pole)0.11796sGTransfer function coefficient (static gain)22856N·m/sThe 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.Table 2
Eigenvalues ofAni−BniKj,i,j=1,2, when only stability is specified.
Eigenvalue ofAn1−Bn1K1Eigenvalue ofAn2−Bn2K2Eigenvalue ofAn1−Bn1K2Eigenvalue ofAn2−Bn2K1−5.863299313142837×105−5.679155054760271×105−5.679155054760272×105−5.863299313142835×105−2.5533482793+j5.1449448541−2.5512811186+j4.6837870864−2.5512811186+j5.1448427531−2.5533482793+j4.6838993430−2.533482793−j5.1449448541−2.5512811186−j4.6837870864−2.5512811186−j5.1448427531−2.5533482793−j4.6838993430The nonlinear termf˜21x1 is described in (24) and (25). Considering the parameters of Table 1, knowing that g = 9.8 m/s2 and taking the operation point angle as θv0=π/6 rad, these values are replaced in (25), obtaining Ma0=3.8028 N·m. For this operation point, the term f˜21x1 belongs to the interval f˜21min≤f˜21x1≤f˜21max, where(63)f˜21min=−29.0297,f˜21max=−24.4981.The new input of the system,PN, is defined from the system input, Pu, and is known as the unreferenced pulse width [11, 22–24]. It is given by(64)PN=Pu−Ma0G.Since the input is the pulse width which is applied on the skin of the patient, its value must be positive, that is,PN>0. So,(65)PN>−Ma0G.


### 5.1. Norm-Bounded Uncertainties

From the parameters presented in Table1, a ±20% variation is considered on the values of J, B, τ, and G, that is,(66)Jmin=0.2305,Jmax=0.3457,Bmin=0.2205,Bmax=0.3307,τmin=0.0944,τmax=0.1416,Gmin=18285,Gmax=27427.Given the minimum and the maximum values off˜21x1, one has the system with two local models, each one described by (8), where(67)An1=010−29.0297−1.03633.615600−8.8307,An2=010−24.4981−1.03633.615600−8.8307,Bn1=00209910,Bn1=00209910,C1=C2=100,and matrices δA1 and δB1 are described by (10), where(68)L1=0000−0.39860.72310000−1.766180734,Δ1=diagδ11,δ12,δ13,δ14,RA1=010001001000,RB1=0001,with −1≤δ1i≤1, for i=1,2,3,4, and matrices δA2 and δB2 are described by (10), where(69)L2=0000−0.39860.72310000−1.766180734,Δ2=diagδ21,δ22,δ23,δ24,RA2=010001001000,RB2=0001,with −1≤δ1i≤1, for i=1,2,3,4.


## 5.1. Norm-Bounded Uncertainties

From the parameters presented in Table1, a ±20% variation is considered on the values of J, B, τ, and G, that is,(66)Jmin=0.2305,Jmax=0.3457,Bmin=0.2205,Bmax=0.3307,τmin=0.0944,τmax=0.1416,Gmin=18285,Gmax=27427.Given the minimum and the maximum values off˜21x1, one has the system with two local models, each one described by (8), where(67)An1=010−29.0297−1.03633.615600−8.8307,An2=010−24.4981−1.03633.615600−8.8307,Bn1=00209910,Bn1=00209910,C1=C2=100,and matrices δA1 and δB1 are described by (10), where(68)L1=0000−0.39860.72310000−1.766180734,Δ1=diagδ11,δ12,δ13,δ14,RA1=010001001000,RB1=0001,with −1≤δ1i≤1, for i=1,2,3,4, and matrices δA2 and δB2 are described by (10), where(69)L2=0000−0.39860.72310000−1.766180734,Δ2=diagδ21,δ22,δ23,δ24,RA2=010001001000,RB2=0001,with −1≤δ1i≤1, for i=1,2,3,4.


## 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.Theorem1, 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, Ai and Bi, i=1,…,q, is decomposed as Ai=Ani+δAi and Bi=Bni+δBi, respectively, with δAi=LiΔiRAi,δBi=LiΔiRBi, and Δi≤1.
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 εiji,j=1,…,q such that the following LMIs are satisfied:(70)Ψii∗∗RAiQ−RBiMi−εiiI∗LiT0−εii−1I<0,1≤i≤q,(71)ϒij∗∗∗∗RAiQ−RBiMj−εijI∗∗∗RAiQ−RBjMi0−εijI∗∗LiT00−εii−1I∗LjT000−εii−1I<0,1≤i≤j≤q,where(72)Ψii=QAniT+AniQ−MiTBniT−BniMi,ϒij=QAniT+AniQ+QAnjT+AnjQ−MjTBniT−BniMj−MiTBnjT−BnjMi,Q=P−1, and Mi=KiP−1, where ∗ denotes the transposed elements in the symmetric positions. From the solution of the aforementioned LMIs, output feedback matrices Ki are obtained from Ki=MiQ−1.
Multiplying (70), at left and at right, by diagI,I,εiiI, the following LMI is obtained:(73)Ψii∗∗RAiQ−RBiMi−εiiI∗εiiLiT0−εiiI<0,1≤i≤q,where Q>0 and(74)Ψii=QAniT+AniQ−MiTBniT−BniMi.
Now, multiplying (71), at left and at right, by I,I,I,εijI,εijI, the following LMI is obtained:(75)ϒij∗∗∗∗RAiQ−RBiMj−εijI∗∗∗RAiQ−RBjMi0−εijI∗∗εijLiT00−εijI∗εijLjT000−εijI<0,1≤i≤j≤q,where Q>0 and(76)ϒij=QAniT+AniQ+QAnjT+AnjQ−MjTBniT−BniMj−MiTBnjT−BnjMi.
So, it is observed from (74), ϒij<0, 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 Q=QT,Mi and scalars εij, for i,j=1,...,q, such that LMIs (73), (75) and Q>0 hold.
For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections5 and 6, the number of rules is q=2. So, to achieve the fuzzy control law (7) that stabilizes system (4), the problem consists of finding matrices Q=QT,M1,andM2 and scalars ε11,ε22,andε12 such that the following LMIs hold:(77)Ψ11∗∗RA1Q−RB1M1−ε11I∗ε11L1T0−ε11I<0,(78)Ψ22∗∗RA2Q−RB12M2−ε22I∗ε22L2T0−ε22I<0,(79)ϒ12∗∗∗∗RA1Q−RB1M2−ε12I∗∗∗RA2Q−RB2M10−ε12I∗∗ε12LiT00−ε12I∗ε12LjT000−ε12I<0,Q>0,where(80)Ψ11=QAn1T+An1Q−M1TBn1T−Bn1M1,Ψ22=QAn2T+An2Q−M2TBn2T−Bn2M2,ϒ12=QAn1T+An1Q+QAn2T+An2Q−M2TBn1T−Bn1M2−M1TBn2T−Bn2M1.
From the solution of the aforementioned LMIs, matricesK1 and K2 are given by(81)K1=M1Q−1,K2=M2Q−1.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,x˙=Aox,x∈ℝn,Ao∈ℝn×n, according to [26], the decay rate is defined as the maximum value of the real constant γ>0 such that(82)limt⟶∞eγtxt=0,∀x0∈ℝn,holds for all trajectories xt.The conditionVx≤−2γVx, in a system x˙=Aox, for all trajectories of xt, for Vx=xTPx,P=PT>0, is equivalent to(83)AoTP+PAo+2γP=Ao+γITP+PAo+γI≤0.Then, replacing matrixAnibyAni+γI 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 γ are(84)Ψii∗∗RAiQ−RBiMi−εiiI∗εiiLiT0−εiiI<0,1≤i≤q,(85)ϒij∗∗∗∗RAiQ−RBiMj−εijI∗∗∗RAiQ−RBjMi0−εijI∗∗εijLiT00−εijI∗εijLjT000−εijI<0,1≤i≤j≤q,where Q>0, and(86)Ψii=QAniT+AniQ−MiTBniT−BniMi+2γQ,ϒij=QAniT+AniQ+QAnjT+AnjQ−MjTBniT−BniMj−MiTBnjT−BnjMi+4γQ,Q=P−1,andMi=KiP−1, where ∗ denotes the transposed elements in the symmetric positions. From the solution of the aforementioned LMIs, output feedback matrices Ki are obtained from Ki=MiQ−1.For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections5 and 6, the number of rules is q=2. 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 Q=QT,M1,andM2 and scalars ε11,ε22,andε12 such that the following LMIs hold:(87)Ψ11∗∗RA1Q−RB1M1−ε11I∗ε11L1T0−ε11I<0,(88)Ψ22∗∗RA2Q−RB12M2−ε22I∗ε22L2T0−ε22I<0,(89)ϒ12∗∗∗∗RA1Q−RB1M2−ε12I∗∗∗RA2Q−RB2M10−ε12I∗∗ε12LiT00−ε12I∗ε12LjT000−ε12I<0,Q>0,where(90)Ψ11=QAn1T+An1Q−M1TBn1T−Bn1M1+2γQ,Ψ22=QAn2T+An2Q−M2TBn2T−Bn2M2+2γQ,ϒ12=QAn1T+An1Q+QAn2T+An2Q−M2TBn1T−Bn1M2−M1TBn2T−Bn2M1+4γQ.From the solution of the aforementioned LMIs, matricesK1andK2 are given by(91)K1=M1Q−1,K2=M2Q−1.


## 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,x˙=Aox,x∈ℝn,Ao∈ℝn×n, according to [26], the decay rate is defined as the maximum value of the real constant γ>0 such that(82)limt⟶∞eγtxt=0,∀x0∈ℝn,holds for all trajectories xt.The conditionVx≤−2γVx, in a system x˙=Aox, for all trajectories of xt, for Vx=xTPx,P=PT>0, is equivalent to(83)AoTP+PAo+2γP=Ao+γITP+PAo+γI≤0.Then, replacing matrixAnibyAni+γI 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 γ are(84)Ψii∗∗RAiQ−RBiMi−εiiI∗εiiLiT0−εiiI<0,1≤i≤q,(85)ϒij∗∗∗∗RAiQ−RBiMj−εijI∗∗∗RAiQ−RBjMi0−εijI∗∗εijLiT00−εijI∗εijLjT000−εijI<0,1≤i≤j≤q,where Q>0, and(86)Ψii=QAniT+AniQ−MiTBniT−BniMi+2γQ,ϒij=QAniT+AniQ+QAnjT+AnjQ−MjTBniT−BniMj−MiTBnjT−BnjMi+4γQ,Q=P−1,andMi=KiP−1, where ∗ denotes the transposed elements in the symmetric positions. From the solution of the aforementioned LMIs, output feedback matrices Ki are obtained from Ki=MiQ−1.For the fuzzy system of the electrostimulation for paraplegic patients, described in Sections5 and 6, the number of rules is q=2. 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 Q=QT,M1,andM2 and scalars ε11,ε22,andε12 such that the following LMIs hold:(87)Ψ11∗∗RA1Q−RB1M1−ε11I∗ε11L1T0−ε11I<0,(88)Ψ22∗∗RA2Q−RB12M2−ε22I∗ε22L2T0−ε22I<0,(89)ϒ12∗∗∗∗RA1Q−RB1M2−ε12I∗∗∗RA2Q−RB2M10−ε12I∗∗ε12LiT00−ε12I∗ε12LjT000−ε12I<0,Q>0,where(90)Ψ11=QAn1T+An1Q−M1TBn1T−Bn1M1+2γQ,Ψ22=QAn2T+An2Q−M2TBn2T−Bn2M2+2γQ,ϒ12=QAn1T+An1Q+QAn2T+An2Q−M2TBn1T−Bn1M2−M1TBn2T−Bn2M1+4γQ.From the solution of the aforementioned LMIs, matricesK1andK2 are given by(91)K1=M1Q−1,K2=M2Q−1.


## 7. System Control with Input Constraint

Consider the input controlut given in (4) and an initial state x0. 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 ut in [24]:(92)maxut2=maxuTtut≤μ0,for all t≥0.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:(93)1x0Tx0Q≥0,QMiTMiμ02I≥0,for an initial state x0andi=1,2,…,q. 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 Sections5 and 6, the number of rules is q=2. So, the LMIs that guarantee the input constraint are(94)1x0Tx0Q≥0,(95)QM2TM1μ02I≥0,(96)QM2TM2μ02I≥0,for an initial state x0. 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 by(97)Pu=PN−Ma0G,whose value must, naturally, be positive. So, the input needs to follow the condition(98)PN>−Ma0G.To follow this condition, input constraint (92) is added for the input, where(99)ut=PN,μ0<Ma0G.Forθv0=π/6, one has μ0<1.6638×10−4. The initial state is(100)x0=θv00−Ma0T=−π60−3.8028T.


## 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 Sections4 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(101)Q=0.000000033491903−0.0000000329330510.000000000438928−0.0000000329330510.000000901932404−0.0000009786037260.000000000438928−0.0000009786037260.003521496953692,M1=0.0000000001326330.0000000019478500.009833371724050,M2=−0.000000000060590−0.0000000006168270.009524542815504,ϵ11=5.47606985268709×10−7,ϵ22=5.28067137087779×10−7,ϵ12=5.72890931362815×10−7.From the solution above, matricesK1andK2 of the control law (7) are(102)K1=3.0594607331076503.1445786641197532.793257471726990,K2=2.9546876169511703.0427255360966762.705530990497405.Table2 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, 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 B,J,τ,andG, from the initial state x0 given in (100), with α described in (17).Figure 3
Simulation result of the nominal fuzzy system, considering only stability.
(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 Figure4, shows this fast convergence.Figure 4
Time zoom onut, for the nominal fuzzy system, considering only stability.In Figures5–7, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, 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 B,J,τ,andG, as described in Table 3. The results show that the controlled system is stable for all possible values of these parameters.Figure 5
Simulation result of the uncertain fuzzy system forx1t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 6
Simulation result of the uncertain fuzzy system forx2t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 7
Simulation result of the uncertain fuzzy system forx3t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Table 3
Values of the parametersB,J,τ,andG for each vertex.
VertexBJτG10.22050.23050.09441828520.22050.23050.09442742730.22050.23050.14161828540.22050.23050.14162742750.22050.34570.09441828560.22050.34570.09442742770.22050.34570.14161828580.22050.34570.14162742790.33070.23050.094418285100.33070.23050.094427427110.33070.23050.141618285120.33070.23050.141627427130.33070.34570.094418285140.33070.34570.094427427150.33070.34570.141618285160.33070.34570.141627427


### 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(103)Q=0.000000006293797−0.000000049443200−0.000000000516894−0.0000000494432000.000000530987030−0.000003234652073−0.000000000516894−0.0000032346520730.011980248303538M1=0.0000000000712810.0000000022129720.031657939760077M2=0.0000000000069860.0000000010802840.031467438088653ϵ11=1.76838667145756×10−6,ϵ22=1.75650519970404×10−6,ϵ12=1.83398991648706×10−6.From the aforementioned solution, matricesK1andK2 of the control law (7) are(104)K1=4.7488325886532150.6042030740259820.026588450497982×102,K2=4.7192562261429570.6004528724567690.026428423547721×102.Table4 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, related to the nominal system, when decay rate γ=4 is specified. In Figure 8, the simulation result is shown for the system with the nominal values of B,J,τ,andG, from the initial state given in (100), with α described in (17).Table 4
Eigenvalues ofAni−BniKj,i,j=1,2, when decay rate γ=4 is specified.
Eigenvalue ofAn1−Bn1K1Eigenvalue ofAn2−Bn2K2Eigenvalue ofAn1−Bn1K2Eigenvalue ofAn2−Bn2K1−5.580376787035519×105−5.546786114873048×105−5.546786114873052×105−5.580376787035501×105−9.1076841424−9.0380952953−9.1077701925−9.0380260712−74.1012735586−74.1553615536−74.0856866552−74.1709316311Figure 8
Simulation result of the nominal fuzzy system, considering decay rateγ=4.
(a)(b)(c)(d)The convergence of the inputut 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.Figure 9
Time zoom onut, for the nominal fuzzy system, considering decay rate γ=4.In Figures10–12, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, from the initial state x0 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 B,J,τ,andG, as described in Table 3. The results show that the controlled system holds the decay rate specification for all possible values of these parameters.Figure 10
Simulation result of the uncertain fuzzy system forx1t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 11
Simulation result of the uncertain fuzzy system forx2t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 12
Simulation result of the uncertain fuzzy system forx3t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(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 ut=PNandμ0=Ma0/G=1.6638×10−4, 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 μ0. 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 μ0=200, the solution of these LMIs is(105)Q=0.000028233983376−0.000007455366912−0.000000768978327−0.0000074553669120.000760237097218−0.000470832034624−0.000000768978327−0.0004708320346248.320872522948498×104,M1=0.0000002385985840.0000020050600065.711656533062334×104,M2=0.000000279742955−0.0000022438071935.684676514801396×104,ε11=2.640736437667783,ε22=2.571381549244325,ε12=2.599534617684041.From the aforementioned solution, matricesK1 and K2 of the control law (7) are(106)K1=0.1029770780903820.4287809456504960.686449446488605,K2=0.1025153633904570.4211791572092130.683206565769427.Table5 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, related to the nominal system, when input constraint (92), with μ0=200, is specified. In Figure 13, the simulation result is shown for the system with the nominal values of B,J,τ,andG, from the initial state given in (100), with x=0 described in (17).Table 5
Eigenvalues ofAni−BniKj,i,j=1,2, when input constraint (92), with μ0=200, is specified.
Eigenvalue ofAn1−Bn1K1Eigenvalue ofAn2−Bn2K2Eigenvalue ofAn1−Bn1K2Eigenvalue ofAn2−Bn2K1−1.440973312768043×105−1.434166564077318×105−1.434166564077322×105−1.440973312768037×105−1.6473497848+j5.1825399965−1.6325949039+j4.7302916980−1.6325949036+j5.1872191983−1.6473497851+j4.7251599715−1.6473497848−j5.1825399965−1.6325949039−j4.7302916980−1.6325949036−j5.1872191983−1.6473497851−j4.7251599715Figure 13
Simulation result of the nominal fuzzy system, considering input constraintμ0=200.
(a)(b)(c)(d)The convergence of the inputut 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.Figure 14
Time zoom onut, for the nominal fuzzy system, considering input constraint μ0=200.In Figures15–17, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, from the initial state x=0 given in (100), with α described in (17). Each figure contains 16 curves, each one representing one combination of the extreme (minimum or maximum) values B,J,τ,andG, as described in Table 3.Figure 15
Simulation result of the uncertain fuzzy system forx1t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 16
Simulation result of the uncertain fuzzy system forx2t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 17
Simulation result of the uncertain fuzzy system forx3t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(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 Figure3, 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 Figure8, 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 Figures8 and 9, the control law assumes too high values, so this control signal is impracticable for the real system.To reduce the magnitude ofu(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.


## 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(101)Q=0.000000033491903−0.0000000329330510.000000000438928−0.0000000329330510.000000901932404−0.0000009786037260.000000000438928−0.0000009786037260.003521496953692,M1=0.0000000001326330.0000000019478500.009833371724050,M2=−0.000000000060590−0.0000000006168270.009524542815504,ϵ11=5.47606985268709×10−7,ϵ22=5.28067137087779×10−7,ϵ12=5.72890931362815×10−7.From the solution above, matricesK1andK2 of the control law (7) are(102)K1=3.0594607331076503.1445786641197532.793257471726990,K2=2.9546876169511703.0427255360966762.705530990497405.Table2 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, 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 B,J,τ,andG, from the initial state x0 given in (100), with α described in (17).Figure 3
Simulation result of the nominal fuzzy system, considering only stability.
(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 Figure4, shows this fast convergence.Figure 4
Time zoom onut, for the nominal fuzzy system, considering only stability.In Figures5–7, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, 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 B,J,τ,andG, as described in Table 3. The results show that the controlled system is stable for all possible values of these parameters.Figure 5
Simulation result of the uncertain fuzzy system forx1t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 6
Simulation result of the uncertain fuzzy system forx2t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 7
Simulation result of the uncertain fuzzy system forx3t, considering only stability. Each curve represents one combination of the extreme (minimum or maximum) values of B, J, τ, and G, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Table 3
Values of the parametersB,J,τ,andG for each vertex.
VertexBJτG10.22050.23050.09441828520.22050.23050.09442742730.22050.23050.14161828540.22050.23050.14162742750.22050.34570.09441828560.22050.34570.09442742770.22050.34570.14161828580.22050.34570.14162742790.33070.23050.094418285100.33070.23050.094427427110.33070.23050.141618285120.33070.23050.141627427130.33070.34570.094418285140.33070.34570.094427427150.33070.34570.141618285160.33070.34570.141627427


## 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(103)Q=0.000000006293797−0.000000049443200−0.000000000516894−0.0000000494432000.000000530987030−0.000003234652073−0.000000000516894−0.0000032346520730.011980248303538M1=0.0000000000712810.0000000022129720.031657939760077M2=0.0000000000069860.0000000010802840.031467438088653ϵ11=1.76838667145756×10−6,ϵ22=1.75650519970404×10−6,ϵ12=1.83398991648706×10−6.From the aforementioned solution, matricesK1andK2 of the control law (7) are(104)K1=4.7488325886532150.6042030740259820.026588450497982×102,K2=4.7192562261429570.6004528724567690.026428423547721×102.Table4 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, related to the nominal system, when decay rate γ=4 is specified. In Figure 8, the simulation result is shown for the system with the nominal values of B,J,τ,andG, from the initial state given in (100), with α described in (17).Table 4
Eigenvalues ofAni−BniKj,i,j=1,2, when decay rate γ=4 is specified.
Eigenvalue ofAn1−Bn1K1Eigenvalue ofAn2−Bn2K2Eigenvalue ofAn1−Bn1K2Eigenvalue ofAn2−Bn2K1−5.580376787035519×105−5.546786114873048×105−5.546786114873052×105−5.580376787035501×105−9.1076841424−9.0380952953−9.1077701925−9.0380260712−74.1012735586−74.1553615536−74.0856866552−74.1709316311Figure 8
Simulation result of the nominal fuzzy system, considering decay rateγ=4.
(a)(b)(c)(d)The convergence of the inputut 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.Figure 9
Time zoom onut, for the nominal fuzzy system, considering decay rate γ=4.In Figures10–12, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, from the initial state x0 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 B,J,τ,andG, as described in Table 3. The results show that the controlled system holds the decay rate specification for all possible values of these parameters.Figure 10
Simulation result of the uncertain fuzzy system forx1t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 11
Simulation result of the uncertain fuzzy system forx2t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 12
Simulation result of the uncertain fuzzy system forx3t, considering decay rate γ=4. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(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 ut=PNandμ0=Ma0/G=1.6638×10−4, 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 μ0. 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 μ0=200, the solution of these LMIs is(105)Q=0.000028233983376−0.000007455366912−0.000000768978327−0.0000074553669120.000760237097218−0.000470832034624−0.000000768978327−0.0004708320346248.320872522948498×104,M1=0.0000002385985840.0000020050600065.711656533062334×104,M2=0.000000279742955−0.0000022438071935.684676514801396×104,ε11=2.640736437667783,ε22=2.571381549244325,ε12=2.599534617684041.From the aforementioned solution, matricesK1 and K2 of the control law (7) are(106)K1=0.1029770780903820.4287809456504960.686449446488605,K2=0.1025153633904570.4211791572092130.683206565769427.Table5 shows the eigenvalues of Ani−BniKj,i,j=1,2, with the aforementioned matrices K1andK2, related to the nominal system, when input constraint (92), with μ0=200, is specified. In Figure 13, the simulation result is shown for the system with the nominal values of B,J,τ,andG, from the initial state given in (100), with x=0 described in (17).Table 5
Eigenvalues ofAni−BniKj,i,j=1,2, when input constraint (92), with μ0=200, is specified.
Eigenvalue ofAn1−Bn1K1Eigenvalue ofAn2−Bn2K2Eigenvalue ofAn1−Bn1K2Eigenvalue ofAn2−Bn2K1−1.440973312768043×105−1.434166564077318×105−1.434166564077322×105−1.440973312768037×105−1.6473497848+j5.1825399965−1.6325949039+j4.7302916980−1.6325949036+j5.1872191983−1.6473497851+j4.7251599715−1.6473497848−j5.1825399965−1.6325949039−j4.7302916980−1.6325949036−j5.1872191983−1.6473497851−j4.7251599715Figure 13
Simulation result of the nominal fuzzy system, considering input constraintμ0=200.
(a)(b)(c)(d)The convergence of the inputut 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.Figure 14
Time zoom onut, for the nominal fuzzy system, considering input constraint μ0=200.In Figures15–17, the simulation result is presented for x1t,x2t,andx3t, respectively, considering all combinations of the minimum or maximum values of the parameters B,J,τ,andG, from the initial state x=0 given in (100), with α described in (17). Each figure contains 16 curves, each one representing one combination of the extreme (minimum or maximum) values B,J,τ,andG, as described in Table 3.Figure 15
Simulation result of the uncertain fuzzy system forx1t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 16
Simulation result of the uncertain fuzzy system forx2t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Figure 17
Simulation result of the uncertain fuzzy system forx3t, considering input constraint μ0=200. Each curve represents one combination of the extreme (minimum or maximum) values of B,J,τ,andG, as described in Table 3.
(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 Figure3, 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 Figure8, 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 Figures8 and 9, the control law assumes too high values, so this control signal is impracticable for the real system.To reduce the magnitude ofu(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.



---

*Source: 4624657-2020-08-12.xml*