Computers, Materials & Continua DOI:10.32604/cmc.2022.022451
Article

Takagi–Sugeno Fuzzy Modeling and Control for Effective Robotic Manipulator Motion

Izzat Al-Darraji1,2, Ayad A. Kakei2, Ayad Ghany Ismaeel3, Georgios Tsaramirsis4, Fazal Qudus Khan5, Princy Randhawa6, Muath Alrammal4 and Sadeeq Jan7,*

1Department of Automated Manufacturing, University of Baghdad, Baghdad, 10001, Iraq
2Department of Mechanical Engineering, University of Kirkuk, Kirkuk, 36001, Iraq
3Department of Computer Technical Engineering, Al-kitab University, Kirkuk, 36001, Iraq
4Higher Colleges of Technology, Abu Dhabi Women’s College, Abu Dhabi, 41012, UAE
5Department of Computer Science, University of Swat, Shangla Campus, Alpurai, 19100, Shangla, Pakistan
6Department of Mechatronics Engineering, Manipal University Jaipur, Jaipur, 302004, India
7Department of Computer Science & IT, University of Engineering & Technology Peshawar, Peshawar, 25000, Pakistan
*Corresponding Author: Sadeeq Jan. Email: sadeeqjan@uetpeshawar.edu.pk
Received: 06 August 2021; Accepted: 07 September 2021

Abstract: Robotic manipulators are widely used in applications that require fast and precise motion. Such devices, however, are prompt to nonlinear control issues due to the flexibility in joints and the friction in the motors within the dynamics of their rigid part. To address these issues, the Linear Matrix Inequalities (LMIs) and Parallel Distributed Compensation (PDC) approaches are implemented in the Takagy–Sugeno Fuzzy Model (T-SFM). We propose the following methodology; initially, the state space equations of the nonlinear manipulator model are derived. Next, a Takagy–Sugeno Fuzzy Model (T-SFM) technique is used for linearizing the state space equations of the nonlinear manipulator. The T-SFM controller is developed using the Parallel Distributed Compensation (PDC) method. The prime concept of the designed controller is to compensate for all the fuzzy rules. Furthermore, the Linear Matrix Inequalities (LMIs) are applied to generate adequate cases to ensure stability and control. Convex programming methods are applied to solve the developed LMIs problems. Simulations developed for the proposed model show that the proposed controller stabilized the system with zero tracking error in less than 1.5 s.

Keywords: Nonlinear robot manipulator; precise fast robot motion; flexible joints; motor friction; Takagy–Sugeno fuzzy control; modeling nonlinear flexible robot system

1  Introduction

Recently, robots have been applied widely in situations that require precise movement at high speeds. In this type of robotic application, the flexibilities in joints [1] and friction [2] are critical factors of the modeling and control processes. To provide precise tracking despite of the existing high nonlinearity effects of joint flexibilities and motor friction, advanced control techniques are essential to be taken in the control design stage. Generally, nonlinearity is an essential issue that was and is still the focus of researchers [3–7]. Techniques such as backstepping control [8], impedance control [9], and sliding mode control [10] are some of the methods used in the control stage. The main issue of these techniques is that the high order nonlinearity for fast motion robotics cannot be solved. In contrast, Takagy–Sugeno Fuzzy Model (T-SFM) is an effective method for representing the high order nonlinear systems in terms of a combination of state equations [11–13]. In this study, initially, the nonlinear state space representation of the robot system is derived. It is assumed that the high nonlinearities of 3rd order in the equation of spring torque of the flexible joint and friction model of Striebeck effect. These nonlinear sources in the state space equations are substituted by rules of T-SFM. The T-SFM is implemented to linearize the derived state space equation. A controller is designed from the developed T-SFM using the method of Parallel Distributed Compensation (PDC). The prime concept of the designed controller is to deduce all the fuzzy rules to compensate for all rules of the fuzzy model. Furthermore, the linear matrix inequalities (LMIs) are applied to generate adequate cases for approving the stability and control purpose issues. Convex programming methods are applied to solve the developed LMIs problems.

This paper is organized as follows. In Section 2, the related work of this study is presented. In Section 3, the state and output equations of the robot arm are derived. In Section 4, the nonlinear robot arm model is represented by the Takagi–sugeno model. Section 5 describes the designe of the controller for the nonlinear robot system. In Section 6, the results of this paper are verified through simulation tests. Finally, the conclusions are written in Section 7.

2  Related Work

Nonlinearity is an essential issue that has been considered by various researchers in the control of robot manipulators. In [14], T-SFM is implemented with a sliding mode controller to solve the issues of the systems having nonlinear dynamic behavior. This proposed controller showed the high efficiency of avoiding chattering within very accurate tracking. In [15], the nonlinear parts in the dynamic system equation are identified by using T-SFM model. Next, using the obtained model from the T-SFM technique, a new controller is designed depending on estimation principle, i.e., not whole measurements. It is shown that nonlinear partial differential equation systems can be controlled using an incomplete number of actuators and sensors. In [16], T-SFM with decomposition controller approach are applied to control the desired trajectory of aircraft with taking into consideration existing both the inaccuracies in the aircraft model and disturbances. The introduced control technique stabilized the closed loop system and tracked in an asymptotically, for a reference step input signal, the pitch angle of the aircraft. In [17], T-SFM approach is applied to represent a partly-active chair suspension system of electrorheological damper. The application of the T-SFM approach has simplified the design process of the H∞ controller. Compared to the existing control technique, the introduced Takagy–Sugeno control technique improved the performance of the the electrorheological damper partly active chair suspension system. In [18], T-SFMs are applied in the robust control of non-constant speed wind turbines which utilizes a generator of twice-fed induction type. The suggested Takagy–Sugeno control method provides the optimum power under considering non-constant wind speed. In [19], an adaptive T-SFM is proposed for piezoelectric actuators to overcome the issue of nonlinearity behavior due to the hysteresis features which degrades the tracking performance. The presented T-SFM showed its efficiency in controlling the piezoelectric actuators without needing the mathematical representation of the hysteresis model. Furthermore, the values of the T-SFM are tuned online to handle the errors of tracking. In [20], an adaptive T-SFM is implemented in the design of a permanent magnetic generator that uses a turbine system. The adaptive T-SFM is assured of the ability to supply electricity in a robust and reliable way. In [21], T-SFM is applied to control an omnidirectional ball robot manipulator of motors which is equipped on two orthogonal planes. The proposed T-SFM control method satisfied high performance results. In which, the model rules corresponding to zero and five degrees limited the omnidirectionally ball robot manipulator performance by narrowing the control range. In [22], TS-FM technique is applied with observers in commercial vehicles to accomplish the appropriate sensors that are necessary to realize precise torque sensing. The rule of the TS-FM was to treat the increasing of nonlinearities in the driving as opposed to load when the speed of the vehicle increased. Based on the above features of applying T-SFM in various applications, this study is focused on the issues of nonlinearities due to joint flexibilities and motor friction consideration for fast and precise motions robot manipulator by T-SFM control.

3  Development of State Space Model

In this section, the state and output equations of the robot arm model are developed. Selection of the appropriate components [23] and modeling is the important step that should be implemented before any development [24,25]. Firstly, equations of motion of the robot arm shown in Fig. 1 are derived. The robot arm model has three main parts: motor, gear, and robot arm [26]. The inertia of the motor, gear and robot arm are denoted by Jmotor, Jgear, and Jarm, respectively. The robot arm system is actuated by the input torque τin of the motor. Besides this, the friction effect is assumed to act on the motor within a nonlinear torsion coefficient Kf,motor. Due to the flexibility effect, the input torque rotates the motor, the gear and robot arm within an unequal angular position θmotor, θgear, and θarm, respectively. The flexibility of the gearbox is modeled by a nonlinear torsion stiffness Kgb. The gear ratio of the gearbox is assumed one, i.e., ngear=1. Since this study focused on the nonlinearity of the joints, the flexibility in the robot arm is modeled by a linear spring of torsion stiffness Karm. Besides the effect of flexibility in the gearbox and the robot arm, the damping coefficients in the gearbox and the robot arm are considered dgb and darm, respectively.

Figure 1: The model of the nonlinear robot arm

The nonlinear torque friction part is considered as the effect of coulomb friction and Striebeck effect. The coulomb friction is assumed as [27].

KC,f,motor=μNθ˙motor,(1)

The total nonlinear torque friction effect is considered as Tustin friction model [28]:

Kf,motor=(KC,f,motor+(KS,f,motor−KC,f,motor)e−θ˙motor10)sign(θ˙motor)+Kv,f,motorθ˙motor(2)

where KS,f,motor, vs, and Kv,f,motor denote static friction, stribeck velocity, and viscous friction, respectively. On the other hand, the nonlinear torsion torque of gear model is assumed as

τgb=kgb,1(θ˙motor−θ˙gear)+kgb,2(θ˙motor−θ˙gear)3(3)

Applying newton second law, the equation of motion for the motor inertia, gear inertia, and arm inertia is

Jmotorθ¨motor=−Kgb(θmotor−θgear)−dgb(θ˙motor−θ˙gear)−Kf,motorθ˙motor+τin,(4)

Jgearθ¨gear=Kgb(θmotor−θgear)+dgb(θ˙motor−θ˙gear)−Karm(θgear−θarm)−darm(θ˙gear−θ˙arm),(5)

Jarmθ¨arm=Karm(θgear−θarm)+darm(θ˙gear−θ˙arm),(6)

Respectively. To get the state space model, assume the following state variables: x1=θmotor, x2=θgear, x3=θarm, x4=θ˙motor, x5=θ˙gear, x6=θ˙arm, and the output is the angular velocity of the motor that is assumed y=x4. thus, the differential states are,

x˙1=x4,(7)

x˙2=x5,(8)

x˙3=x6,(9)

x˙4=1Jmotor(−Kgb(x1−x2)−dgb(x4−x5)−Kf,motorx4+τin),(10)

x˙5=1Jgear(Kgb(x1−x2)+dgb(x4−x5)−Karm(x2−x3)−darm(x5−x6)),(11)

x˙6=1Jarm(Karm(x2−x3)+darm(x5−x6)),(12)

The differential states from Eqs. (7)–(12) and the output formula, i.e., the angular velocity of the motor, can be arranged in state space model format as

[x˙1x˙2x˙3x˙4x˙5x˙6]=A(t)[x1x2x3x4x5x6]+B(t),(13)

y=C(t),

where

A(t)=[000100000010000001−KgbJmotorKgbJmotor0−dgb−Kf,motorJmotordgbJmotor0KgbJgear−Kgb−KarmJgearKarmJgeardgbJgear−dgb−darmJgeardarmJgear0KarmJarm−KarmJarm0darmJarm−darmJarm],

B(t)=[000τin00],

C(t)=x4.

The obtained state model in Eq. (13) will be implemented in the next section to be represented by the T-SFM approach. The numerical results of this study are based on the physical parameters of the robot model presented in Tab. 1.

After deriving the state equations of our proposed nonlinear robot model, T-SFM can be implemented to linearize the robot system as explained in the next section.

4  Implementation of T-SFM

T-SFM is implemented to represent the nonlinear state space model of the robot arm. In addition, the T-SFM has linearized the nonlinear terms in state space equations. The nonlinear robot arm model will be denoted by T-SFM as

x˙(t)=∑j=1rωj(z(t))(Ajx(t)+Bju(t)+aj),(14)

y(t)=∑j=1rωj(z(t))(Cjx(t)+ci),(15)

where r,Z(t), ωj(z(t)), and “aj & cj” denotes the number of local models, scheduling variables vector, normalized membership function, and the biases of the jth local model, respectively. Within T-SFM, the nonlinear robot arm model is exemplified in a close combination of state variables. Consider the range of the state variables as follow: x1(t)∈[0,2π], x2(t)∈[0,2π], x3(t)∈[0,2π], x4(t)∈[0,10], x5(t)∈[0,10], and x6(t)∈[0,10]. In Eq. (13), the nonlinear elements due to nonlinear properties of Kf,motor and Kgb are assumed: z1(t)=−KgbJmotor, z2(t)=KgbJmotor, z3(t)=−dgb−Kf,motorJmotor, and z4(t)=−Kgb−KarmJgear. Inserting these formulas in the value of matrix A of Eq. (13), results:

A(t)=[000100000010000001z1(t)z2(t)0z3(t)dgbJmotor0z2(t)z4(t)KarmJgeardgbJgear−dgb−darmJgeardarmJgear0KarmJarm−KarmJarm0darmJarm−darmJarm],(16)

Under our considered range values, the minimum and maximum values of z1(t), z2(t), z3(t), z4(t) are calculated considering the values of the parameters which are involved in their formulas. In terms of Eq. (14), the scheduling variables can be represented now as

z1(t)=L1(z1(t)).Max z1(t)+L2(z1(t)).Min z1(t),(17)

z2(t)=M1(z2(t)).Max z2(t)+M2(z1(t)).Min z2(t),(18)

z3(t)=N1(z3(t)).Max z3(t)+N2(z3(t)).Min z3(t),(19)

z4(t)=T1(z4(t)).Max z4(t)+T2(z4(t)).Min z4(t)(20)

where L1(z1(t))+L2(z1(t))=1, M1(z2(t))+M2(z1(t))=1, N1(z3(t))+N2(z3(t))=1, and T1(z4(t))+T2(z4(t))=1. Consequently, the membership functions are obtained as:

L1=z1(t)−Min z1(t)z1(t).Max z1(t)−z1(t).Min z1(t),(21)

L2=z1(t)−Max z1(t)z1(t).Min z1(t)−z1(t)Max z1(t),(22)

M1=z2(t)−Min z2(t)z2(t).Max z2(t)−z2(t).Min z2(t),(23)

M2=z2(t)−Max z2(t)z2(t).Min z2(t)−z2(t)Max z2(t),(24)

N1=z3(t)−Min z3(t)z3(t).Max z3(t)−z3(t).Min z3(t),(25)

N2=z3(t)−Max z3(t)z3(t).Min z3(t)−z3(t)Max z3(t),(26)

T1=z4(t)−Min z4(t)z4(t).Max z4(t)−z4(t).Min z4(t),(27)

T2=z4(t)−Max z4(t)z4(t).Min z4(t)−z4(t)Max z4(t),(28)

these membership functions L1, L2, M1, M2, N1, N2, T1, and T2 are named as Big_1, Small_1, Big_2, Small_2, Big_3, Small_3, Big_4, Small_4, respectively.

5  Controller Design

The presented design of the controller is closely related to the feature of the derived T-SFM in the previous section. This feature is that the terms of state equations z1(t), z2(t), z3(t), and z4(t) are not constant. These terms are varying corresponding to the torsion torque Kgb and motor friction Kf,motor depends on the velocities of motor and gear of the state variables. This property was the source of the high nonlinearity in the robot arm model. Consequently, for linearization purpose, the nonlinear robot arm system is expressed by the bellow fuzzy model:

Model rule 1

If z1(t) is Small_1 and z2(t) is Small_2 and z3(t) is Small_3 and z4(t) is Small_4 then

x˙(t)=A1x(t)+B1u(t)

Model rule 2

If z1(t) is Small_1 and z2(t) is Small_2 and z3(t) is Small_3 and z4(t) is Big_4 then

x˙(t)=A2x(t)+B2u(t)

Model rule 3

If z1(t) is Small_1 and z2(t) is Small_2 and z3(t) is Big_3 and z4(t) is Small_4 then

x˙(t)=A3x(t)+B3u(t)

Model rule 4

If z1(t) is Small_1 and z2(t) is Small_2 and z3(t) is Big_3 and z4(t) is Big_4 then

x˙(t)=A4x(t)+B4u(t)

Model rule 5

If z1(t) is Small_1 and z2(t) is Big_2 and z3(t) is Small_3 and z4(t) is Small_4 then

x˙(t)=A5x(t)+B5u(t)

Model rule 6

If z1(t) is Small_1 and z2(t) is Big_2 and z3(t) is Small_3 and z4(t) is Big_4 then

x˙(t)=A6x(t)+B6u(t)

Model rule 7

If z1(t) is Small_1 and z2(t) is Big_2 and z3(t) is Big_3 and z4(t) is Small_4 then

x˙(t)=A7x(t)+B7u(t)

Model rule 8

If z1(t) is Small_1 and z2(t) is Big_2 and z3(t) is Big_3 and z4(t) is Big_4 then

x˙(t)=A8x(t)+B8u(t)

Model rule 9

If z1(t) is Big_1 and z2(t) is Small_2 and z3(t) is Small_3 and z4(t) is Small_4 then

x˙(t)=A9x(t)+B9u(t)

Model rule 10

If z1(t) is Big_1 and z2(t) is Small_2 and z3(t) is Small_3 and z4(t) is Big_4 then

x˙(t)=A10x(t)+B10u(t)

Model rule 11

If z1(t) is Big_1 and z2(t) is Small_2 and z3(t) is Big_3 and z4(t) is Small_4 then

x˙(t)=A11x(t)+B11u(t)

Model rule 12

If z1(t) is Big_1 and z2(t) is Small_2 and z3(t) is Big_3 and z4(t) is Big_4 then

x˙(t)=A12x(t)+B12u(t)

Model rule 13

If z1(t) is Big_1 and z2(t) is Big_2 and z3(t) is Small_3 and z4(t) is Small_4 then

x˙(t)=A13x(t)+B13u(t)

Model rule 14

If z1(t) is Big_1 and z2(t) is Big_2 and z3(t) is Small_3 and z4(t) is Big_4 then

x˙(t)=A14x(t)+B14u(t)

Model rule 15

If z1(t) is Big_1 and z2(t) is Big_2 and z3(t) is Big_3 and z4(t) is Small_4 then

x˙(t)=A15x(t)+B15u(t)

Model rule 16

If z1(t) is Big_1 and z2(t) is Big_2 and z3(t) is Big_3 and z4(t) is Small_4 then

x˙(t)=A16x(t)+B16u(t)

In this way, the nonlinear robot model can be described in the following general formula

Modelrulej:Ifz1(t)isM1jandz2(t)isM2jandz3(t)isM3jandz4(t)isM4jThenx˙(t)=Ajx(t)+Bju(t),j=1,…,16(29)

Applying the defuzzification principle of fuzzy method, the output of the system is obtained as:

x˙(t)=∑j=1rωj(z(t))(Ajx(t)+Bju(t)+aj)∑j=1rωj(z(t))(30)

y(t)=∑j=1rωj(z(t))(Cjx(t)+ci)∑j=1rωj(z(t))(31)

where ωj(z(t)) denotes the membership function for the model rule j. Hence

ωj(z(t))=∏v=1rMvj(xv(t))(32)

For the obtained T-S fuzzy model rule, the technique of state feedback is applied to develop the following control rules

Modelrulej:If z1(t)isM1jandz2(t)isM2jandz3(t)isM3jandz4(t)isM4jThenu(t)=Kjx(t),j=1,…,16(33)

The PDC technique is implemented in this study to find the solution of the system using the obtained T-SFM. Considering Eqs. (29) and (33), PDC technique is applied to design the controller based on T-S fuzzy approach as follow

u(t)=∑j=116ωj(z(t))Kjx(t)∑j=116ωj(z(t))(34)

Our developed mode rules in Eq. (29) can be asymptotically stable by a potential positive matrix Q when satisfying the following conditions:

QAjT+AjQ+VjTBjT+BjVj<0,j=1,2,…,16(35)

QAjT+AjQ+QAkT+AkQ+VkTBjT+BjVk+VjTBkT+BkVj<0,j<k≤16(36)

Q=P−1>0(37)

where Vj=KjQ. Consequently, Kj in Eq. (34) can be calculated by implementing the above LMI conditions.

6  A Numerical Example

In this section, simulation tests with MatLab R2017b are implemented to verify the performance of the designed controller for the robot arm model of the parameters that are listed in Tab. 1 with initial arm position and angular velocity 2∘ and 0rad/s, respectively. The control issue of the robot arm system shown in Fig. 1 is assumed to reach the robot arm a desired position by applying a motor torque input. Furthermore, it is assumed that the robot arm is stabilized at a specific angle. The robot arm model is specified in Eq. (13). By using Eq. (35), for 16 rules, we get 16 LMI of the robot arm as follows:

QA1T+A1Q+V1TB1T+B1V1<0,(38)

QA2T+A2Q+V2TB2T+B2V2<0,(39)

QA3T+A3Q+V3TB3T+B3V3<0,(40)

QA4T+A4Q+V4TB4T+B4V4<0,(41)

QA5T+A5Q+V5TB5T+B5V5<0,(42)

QA6T+A6Q+V6TB6T+B6V6<0,(43)

QA7T+A7Q+V7TB7T+B7V7<0,(44)

QA8T+A8Q+V8TB8T+B8V8<0,(45)

QA9T+A9Q+V9TB9T+B9V9<0,(46)

QA10T+A10Q+V10TB10T+B10V10<0,(47)

QA11T+A11Q+V11TB11T+B11V11<0,(48)

QA12T+A12Q+V12TB12T+B12V12<0,(49)

QA13T+A13Q+V13TB13T+B13V13<0,(50)

QA14T+A14Q+V14TB14T+B14V14<0,(51)

QA15T+A15Q+V15TB15T+B15V15<0,(52)

QA16T+A16Q+V16TB16T+B16V16<0,(53)

where Q is obtained as explained in Eq. (37). For j<k≤16, from Eq. (36), we have QAjT+AjQ+QAkT+AkQ+VkTBjT+BjVk+VjTBkT+BkVj<0. Hence, there are 240 LMI can be designed~ as:

•   j=1, k=2, j=1, k=3, j=1, k=4, j=1, k=5, j=1, k=6, j=1, k=7, j=1, k=8, j=1, k=9, j=1, k=10, j=1, k=11, j=1, k=12, j=1, k=13, j=1, k=14, j=1, k=15, j=1, k=16, 

•   j=2, k=3, j=2, k=4, j=2, k=5, j=2, k=6, j=2, k=7, j=2, k=8, j=2, k=9, j=2, k=10, j=2, k=11, j=2, k=12, j=2, k=13, j=2, k=14, j=2, k=15, j=2, k=16, 

•   j=3, k=4, j=3, k=5, j=3, k=6, j=3, k=7, j=3, k=8, j=3, k=9, j=3, k=10, j=3, k=11, j=3, k=12, j=3, k=13, j=3, k=14, j=3, k=15, j=3, k=16,

•   j=4, k=5, j=4, k=6, j=4, k=7, j=4, k=8, j=4, k=9, j=4, k=10, j=4, k=11, j=4, k=12, j=4, k=13, j=4, k=14, j=4, k=15, j=4, k=16,

•   j=5, k=6, j=5, k=7, j=5, k=8, j=5, k=9, j=5, k=10, j=5, k=11, j=5, k=12, j=5, k=13, j=5, k=14, j=5, k=15, j=5, k=16,

•   j=6, k=7, j=6, k=8, j=6, k=9, j=6, k=10, j=6, k=11, j=6, k=12, j=6, k=13, j=6, k=14, j=6, k=15, j=6, k=16, 

•   j=7, k=8, j=7, k=9, j=7, k=10, j=7, k=11, j=7, k=12, j=7, k=13, j=7, k=14, j=7, k=15, j=7, k=16,

•   j=8, k=9, j=8, k=10, j=8, k=11, j=8, k=12, j=8, k=13, j=8, k=14, j=8, k=15, j=8, k=16,

•   j=9, k=10, j=9, k=11, j=9, k=12, j=9, k=13, j=9, k=14, j=9, k=15, j=9, k=16,

•   j=10, k=11, j=10, k=12, j=10, k=13, j=10, k=14, j=10, k=15, j=10, k=16,

•   j=11, k=12, j=11, k=13, j=11, k=14, j=11, k=15, j=11, k=16,

•   j=12, k=13, j=12, k=14, j=12, k=15, j=12, k=16, 

•   j=13, k=14, j=13, k=15, j=13, k=16,

•   j=14, k=15, j=14, k=16,

•   j=15, k=16.

One of the essential considerations of LMI design is the interconnection of membership j and k. For instance, assume j=11 and k=12. Thus, the interaction between the 11th fuzzy rule and the 12th fuzzy rule should be considered. The obtained LMI in this case is

QA11T+A11Q+QA12T+A12Q+V12TB11T+B11V12+V11TB12T+B12V11<0,(54)

Referring to Eq. (41), T-SFM controller has been designed by applying PDC technique

u=w1(z1(t))K1x(t)+w2(z2(t))K2x(t)+w3(z3(t))K3x(t)+w4(z4(t))K4x(t)+w5(z5(t))K5x(t)+w6(z6(t))K6x(t)+w7(z7(t))K7x(t)+w8(z8(t))K8x(t)+w9(z9(t))K9x(t)+w10(z10(t))K10x(t)+w11(z11(t))K11x(t)+w12(z12(t))K12x(t)+w13(z13(t))K13x(t)+w14(z14(t))K14x(t)+w15(z15(t))K15x(t)+w16(z16(t))K16x(t)(55)

where Ki=ViP, i=1,2,…,16. Utilizing both of LMI and YALMIP toolboxes, the values of Q and Vi are obtained. Consequently, Ki is obtained as:

K1=[4263390868115229],K2=[1298908474829502],K3=[6015891876714980],K4=[3173442379126245],K5=[8903956238914291],K6=[4189692298265229],K7=[5103798467123108],K8=[6186961288665439],K9=[−7139−4700−6559−5999],K10=[−5654−6357−7383−9072],K11=[−8049−4568−9641−4689],K12=[−7678−6473−7213−6541],K13=[−6984−8867−6293−4796],K14=[−5180−8744−1746−6244],K15=[−3523−7650−8102−9719],K16=[−5780−7217−6834−7422].