A Data-Driven Approach of Takagi-Sugeno Fuzzy Control of Unknown Nonlinear Systems

: A novel approach to build a Takagi-Sugeno (T-S) fuzzy model of an unknown nonlinear system from experimental data is presented in the paper. The neuro-fuzzy models or, more speciﬁcally, fuzzy basis function networks (FBFNs) are trained from input–output data to approximate the nonlinear systems for which analytical mathematical models are not available. Then, the T-S fuzzy models are derived from the direct linearization of the neuro-fuzzy models. The operating points for linearization are chosen using the evolutionary strategy to minimize the global approximation error so that the T-S fuzzy models can closely approximate the original unknown nonlinear system with a reduced number of linearizations. Based on T-S fuzzy models, optimal controllers are designed and implemented for a nonlinear two-link ﬂexible joint robot, which demonstrates the possibility of implementing the well-established model-based optimal control method onto unknown nonlinear dynamic systems.


Introduction
The Takagi-Sugeno (T-S) fuzzy model is a powerful and practical engineering tool for modeling and control of complex nonlinear systems. It proves to be a universal function approximator that can approximate any smooth nonlinear functions to any degree of accuracy [1,2] and is less sensitive to the curse of dimensionality than other fuzzy models [3]. The concept of T-S fuzzy model is similar to the piecewise linear approximation approaches in nonlinear control, which linearizes a system at a set of selected operating points and designs a local linear feedback controller for each linear model. However, since the overall control action is switching among the local linear controllers according to system states and thus there is only one local controller active at a certain time in such approaches, it can only ensure the stability and performance of the control system at the neighborhood of selected operating points [4], In contrast, the T-S fuzzy model approximates the entire nonlinear system by fuzzy inference among local linear models so that the overall control action can be generated by aggregation of local linear control laws [5]. Therefore, it empowers a paradigm of designing controllers for local linear models while analyzing stability for the global nonlinear system [6]. The T-S fuzzy-model-based control that blends feedback controllers from local models is referred to as Parallel Distributed Compensation (PDC) scheme, in which the stability of the overall control system is assessed through Lyapunov stability analysis, especially by the Linear Matrix Inequality (LMI) technique [7][8][9][10][11][12].
To take advantage of T-S fuzzy-model-based control, the identification of T-S fuzzy model has attracted great research interest. There are two kinds of methods for establishing T-S fuzzy models. One is linearizing the original system at a series of operating points when an analytical model of the system is available. The other is the consecutive structure and parameter identification from the data generated by the unknown system [6], which is more of interest to us. The structure identification refers to the selection of locations of fuzzy where R k denotes the kth fuzzy rule, k ∈ {1, 2, · · · , p}, p is the number of fuzzy rules, F j k (j = 1, 2, · · · , v) are the input fuzzy sets, x(t) ∈ R n is the state variable vector, u(t) ∈ R m is the input variable vector, z(t): = [z 1 , z 2 , · · · , z v ] are a subset of measurable or observable variables in the state and input vectors that are used for fuzzification, and (A k , B k , d k ) are the matrices of the kth local model [6]. If the constant bias vector is not null, d k = 0, for some fuzzy inference rules, the corresponding local models are affine models instead of linear models. By using the product fuzzy inference, singleton output membership functions, and centroid defuzzifier, a T-S fuzzy model in continuous-time state-space form can be organized as: where A(µ) = ∑ p k=1 µ k A k , B(µ) = ∑ p k=1 µ k B k , d(µ) = ∑ p k=1 µ k d k and µ k denotes the normalized fuzzy membership function represents the i-th membership function in fuzzy set F i k of the k-th fuzzy rule. Due to the fact that the membership functions are nonlinear (e.g., triangular or Gaussian), the model in Equation (1), as an aggregation of local linear or affine models, is also a nonlinear model in nature.

Identification of Takagi-Sugeno Fuzzy Model
The T-S fuzzy model is expected to attain a good approximation of not only the local dynamics of the underlying system to facilitate local compensator design, but also the global dynamics to guarantee the overall control performance. However, considering them together is difficult because building constituent local linear models from input-output data is not always straightforward and the tradeoff between the local and global accuracy should be carefully addressed [2].
In this paper, a novel approach that utilizes the neuro-fuzzy model to obtain a T-S fuzzy model is proposed. The neuro-fuzzy models are based on the fusion of fuzzy inference systems and neural networks. Among diverse neuro-fuzzy models, the one proposed by Wang and Mendel [20], which is also known as fuzzy basis function network (FBFN), has gained much attention as it has a similar structure to radial basis function networks (RBFNs) and thus can adopt the training methods that are already established for the RBFN. It has been proved that the FBFNs can uniformly approximate any real continuous nonlinear functions on a compact set to arbitrary accuracy [21].
First, from the input-output data, a set of multi-input-single-output FBFNs with product fuzzy inference, Gaussian membership functions (MFs), and centroid defuzzifier are trained by the least square algorithm [22] to approximate the nonlinear system whose analytical mathematical model is not available. Each FBFN can be represented as a fuzzy system in layered network form (see Figure 1), and is used to approximate the dynamics of one state variable. For example, the FBFN for the p-th state variable is constructed using fuzzy rules as: where l ∈ {1, 2, · · · , M p }, and Mp is the number of fuzzy rules for the p-th FBFN. Through fuzzy inference and defuzzification, the FBFN based on the above fuzzy rules can be written as: .
where x: = [x 1 , x 2 , · · · , x n ] is the state variable vector, u: = [u 1 , u 2 , · · · , u m ] is the input vector, and w p l is the weighting factor. The x ip l and σ ip l are the centers and widths of in fuzzy set A ip l , and the u ip l and σ uip l are the centers and widths of Gaussian MFs e − 1 2 ( where H.O.T denotes the higher order terms in the local model, which w in the following analysis, Similarly, Next, the T-S fuzzy model is derived from the neuro-fuzzy model by linearizing Equation (2). Linearization about an operating point (x k, u k ) results in where H.O.T denotes the higher order terms in the local model, which will be neglected in the following analysis, f(x k , u k ) = [ f 1 (x k , u k ), . . . , f n (x k , u k )] T , and By taking the partial derivative of Equation (2) with respect to x q , one can get where Similarly, where ∂ f p (x, u)/∂u q has the same expression as Equation (5) except that: For simplicity, and also by neglecting the higher-order terms, Equation (3) can be rewritten as where The affine term d k is non-null, i.e., d k = 0 in general even if the operating point is an equilibrium point other than the origin [23]. Including the affine terms will offer a more accurate local approximation of the system's dynamics around the operating points.
Then, the operating points used for linearization need to be chosen carefully to achieve a good global approximation of the original system. Meanwhile, it is desirable to minimize the number of linearization so that the number of fuzzy rules in the T-S fuzzy model can be reduced and the fuzzy controller synthesized in the following section will be more computationally efficient. To reduce the number of linearization points, it is necessary to know which variables among z = [z 1 , z 2 , · · · , z v ] in fuzzy inference are the major sources of nonlinearity. These variables should be assigned more positions while the variables that are minor sources can be assigned fewer positions. A measure to roughly recognize the source of nonlinearity by inspecting the variation of linearized A-matrix is presented in Section 5. If the i-th variable z i is to be assigned p i positions {z i,1 , . . . , z i,pi }, then by combination of the ∑ v i=1 p i positions, there will be p = ∏ v i=1 p i operating points in total, i.e., (x k , u k ) (k = 1, 2, · · · , p).
After the number of design positions p i for each z i is selected, their optimal positions can be searched by the Evolutionary Strategy (ES) [24] inside the operating range to produce the p optimal operating points to build the T-S fuzzy model: The non-dimensional error index (NDEI) between the responses predicted by the T-S fuzzy model .
x i,TS [k] and the measurement data .
x i,M [k] is a reliable criterion of the quality of global approximation and can be used as the objective function to be minimized. The input-output data pairs used for validation of the neuro-fuzzy model can be reused here so that no extra experiments need to be conducted. N is the number of pairs in the validation data set, n is the number of state variables and .
x i,M is the mean value of .
The ES optimization is controlled by three parameters: the maximum number of generations t, the parent population µ, and the offspring population λ. The optimization starts with initializing a parent pool with µ individuals. Then, in each generation, a pair of parents is randomly selected to produce an offspring via recombination and mutation, as explained in [24], until λ offspring have been generated. Each offspring is used to build a T-S model at the operating point associated with it and the NDEI of this model is recorded. The best µ offspring will form the parent pool for the next generation. If the optimized set of operating points after t generations doesn't yield a T-S fuzzy model with satisfactory global approximation, the number of operating points p will be increased, and the optimization will be repeated, as illustrated in Figure 2.
The ES optimization is controlled by three parameters: the maximum numbe erations t, the parent population µ, and the offspring population λ. The optimizati with initializing a parent pool with µ individuals. Then, in each generation, a pai ents is randomly selected to produce an offspring via recombination and mutatio plained in [24], until λ offspring have been generated. Each offspring is used to b S model at the operating point associated with it and the NDEI of this model is r The best µ offspring will form the parent pool for the next generation. If the optim of operating points after t generations doesn't yield a T-S fuzzy model with sat global approximation, the number of operating points p will be increased, and the zation will be repeated, as illustrated in Figure 2.

Controller Synthesis for the T-S Fuzzy Model
In order to design a global stabilizing controller for the nonlinear system u identified T-S fuzzy model, the parallel distributed compensation (PDC) frame adopted [25]. Local feedback rule is designed as a compensator for each local mo a global fuzzy controller is constructed by the aggregation of local compensato the same fuzzy inference system in the T-S fuzzy model [6]: The fuzzy controller is aggregated as: where µk is the normalized membership function same as in Equation (1). Assume that the fuzzy controller is designed to minimize the performance i

Controller Synthesis for the T-S Fuzzy Model
In order to design a global stabilizing controller for the nonlinear system using the identified T-S fuzzy model, the parallel distributed compensation (PDC) framework is adopted [25]. Local feedback rule is designed as a compensator for each local model and a global fuzzy controller is constructed by the aggregation of local compensators using the same fuzzy inference system in the T-S fuzzy model [6]: The fuzzy controller is aggregated as: where µ k is the normalized membership function same as in Equation (1). Assume that the fuzzy controller is designed to minimize the performance index: where r is the command input vector to be tracked, x is the state vector, u is the input vector, t 0 is the initial time, t f is the final time and Q and R are symmetric positive definite matrices to be determined by the designer, then for each local controller, the optimal control action u k can be derived as [26]: where K k is given by Appl. Sci. 2021, 11, 62 7 of 15 P k is found by solving the Continuous Algebra Riccati Equation: The stability condition of an affine T-S fuzzy control system based on a quadratic Lyapunov function is given in [27]. The equilibrium point (x = x k , u = u k ) of the control system is asymptotically stable in large if there exists a common positive definite matrix P = P T > 0 and scalars τ ijq ≥ 0 such that: where the G i,j , η i,j , T ijq , u ijq and v ijq are defined as: Equation (16) belongs to bilinear matrix inequalities (BMIs) and can be solved in an iterative LMI manner. The details of the iterative linear matrix inequality (ILMI) algorithm can be found in [27]. The stability condition in Equation (16) can be integrated into the ES optimization procedure shown in Figure 2 as a constraint so that the fuzzy control system designed based on the optimized operating point set is guaranteed to be stable.
The fuzzy controller is implemented in a full state-feedback manner. If some state variables are not measurable during operation, then an observer needs to be designed for each local model and the fuzzy observer is constructed by aggregation of local observers with a fuzzy inference system. In [28], it has been proved that the separation principle for linear systems also holds for T-S fuzzy systems, and thus the fuzzy controller and fuzzy observer can be designed independently.

Application Example
In this section, the data-driven T-S fuzzy model identification and control are implemented on a flexible two-link joint robot to demonstrate the effectiveness of the proposed approach.
Flexible robot manipulators possess various advantages over the rigid ones: they require less material, allow higher manipulation speed while consume less power, and are safer to operate due to reduced inertia. However, controlling flexible robot manipulators for precise positioning could a challenging task because of the high precision required for positioning, oscillation due to flexibility, highly nonlinear and distributed dynamics of the system, as well as the difficulty in establishing an accurate model [29]. The picture and the schematic of the two-link flexible-joint robot manipulator to be dealt with in this paper are shown in Figure 3.
The robot is described by 8 state variables: θ 1 , angle of the 1st link; θ 2 , angle of the 2nd link; θ 3 , angle of the 1st motor; θ 4 , angle of the 2nd motor; and the four angular velocities. There are two input variables: T 1 , the torque of the 1st motor and T 2 , the torque of the 2nd motor. The state vector and input vector of this system are defined as: quire less material, allow higher manipulation speed while consume less power, and are safer to operate due to reduced inertia. However, controlling flexible robot manipulators for precise positioning could a challenging task because of the high precision required for positioning, oscillation due to flexibility, highly nonlinear and distributed dynamics of the system, as well as the difficulty in establishing an accurate model [29]. The picture and the schematic of the two-link flexible-joint robot manipulator to be dealt with in this paper are shown in Figure 3. The robot is described by 8 state variables: θ1, angle of the 1st link; θ2, angle of the 2nd link; θ3, angle of the 1st motor; θ4, angle of the 2nd motor; and the four angular velocities. There are two input variables: T1, the torque of the 1st motor and T2, the torque of the 2nd motor. The state vector and input vector of this system are defined as: The nonlinear equation of motion of this robot can be expressed as [30]: where   M θ is the inertia matrix,   , V θ θ  is the vector of Coriolis and centrifugal functions, C is the viscous damping coefficient matrix,   D θ  is the Coulombic friction vector, K is the stiffness coefficient matrix, and T is the input torque vector. The inertia matrix is given by: The nonlinear equation of motion of this robot can be expressed as [30]: where M(θ) is the inertia matrix, V θ, θ is the Coulombic friction vector, K is the stiffness coefficient matrix, and T is the input torque vector. The inertia matrix is given by: r 2 M 1 (θ) = p 1 + m 2 a 2 2 + J 2 + 2(l 1 m 2 a 2 ) cos(θ 2 ) m 2 a 2 2 + J 2 + (l 1 m 2 a 2 ) cos(θ 2 ) m 2 a 2 2 + J 2 + (l 1 m 2 a 2 ) cos(θ 2 ) m 2 a 2 2 + J 2 (19) where m i is the lumped mass of components, J i is the moment of inertia of components, l i is the length of links, and a 1 and a 2 denote the offset from the center of gravity of the first and second link to the first and second joint, respectively. In addition, b 1 is the distance between the second motor and the first joint, and r is the gear ratio of the chain drives. The vector of the Coriolis and centrifugal functions is and the viscous damping matrix can be written as where c i is the viscous friction coefficients. The vector of Coulombic friction is given by where d i is the friction torque at each joint. The matrix of stiffness coefficients is given by where k i denotes the coefficients of the torsional springs at the flexible joints. The torque vector is Table 1 lists the estimated values of the robot's physical parameters. However, deriving the mathematical model and obtaining an accurate estimation of each parameter is quite difficult and time-consuming. The purpose of this paper is to develop a methodology of T-S model construction and controller design when the analytical mathematical model detailed above is not available. Hence, the Equations (18)-(24) and the parameters in Table 1 are only used to validate the fuzzy controller via simulation. In experiments, four encoders are mounted to measure the four arm and motor angles, and the angular velocities and accelerations are obtained from the angular position signals by the central finite difference method: where θ i [k], k = 1, 2, · · · are the values of discrete angular position measurements and ∆T = 0.001 s is the sampling time. To reduce noises in the velocity and acceleration signals that mainly originate from the quantization of the position signal, a zero-phase lowpass Butterworth filter with 150 Hz cutoff frequency was applied to the position signals.
To avoid the transient effect from filtering in both directions when initializing the filter states, experiments were started with a 0.5 s rest period. This period and the last 0.5 s of the experiment were removed from the data set.
To obtain the neuro-fuzzy model of the robot, the motors were excited by sine sweep torques. A combination of a 5-s 1 Hz sine wave and a subsequent 5-s sine sweep signal with an initial frequency of 1 Hz and a final frequency of 5 Hz was used to excite the shoulder motor and elbow motor, as shown in Figure 4. Then, 2000 training data are evenly drawn from the collected input-output pairs and used to train the neuro-fuzzy model. 500 testing data that are different from the training data are used to test the accuracy of the neuro-fuzzy model and the subsequent T-S fuzzy models. The T-S fuzzy models in this paper are all constructed with triangular membership functions due to their simplicity.
FOR PEER REVIEW 11 of 16 shoulder motor and elbow motor, as shown in Figure 4. Then, 2000 training data are evenly drawn from the collected input-output pairs and used to train the neuro-fuzzy model. 500 testing data that are different from the training data are used to test the accuracy of the neuro-fuzzy model and the subsequent T-S fuzzy models. The T-S fuzzy models in this paper are all constructed with triangular membership functions due to their simplicity. To determine the proper number of operating points for linearization, the grade of the nonlinearity of each state variable is evaluated as follows: first, the system is linearized at the equilibrium point x = 0, where the result is Ax=0. Then by randomly changing one state variable xi in the operating range and fix all the other state variables, a series of matrices Axi are obtained. They are compared to Ax=0 in terms of the 2-norm of difference. The results are summarized in Table 2. The idea is that if the state variable xi is linear or affine in the system, then identical A matrices should be obtained at different positions of xi. Otherwise, xi should be regarded as nonlinear. From Table 2, it can be seen that the variables 2  , 1   and 2   are sources of nonlinearity around the origin, among which θ2 is the most significant one. The other variables can be treated as linear.   Table 2.
In Figure 5, it can be seen that the outputs predicted by the T-S fuzzy models are very close to the testing experimental data. The NDEIs of the neuro-fuzzy model and the two T-S fuzzy models, which are used as the criteria of global approximation, are listed in Table 3. It can be seen that the T-S fuzzy models derived from the neuro-fuzzy model can provide an accurate approximation of the original nonlinear system. The 12 points T-S fuzzy model with optimized operating points gives an even better approximation than the 27 points one with uniformly distributed operating points.  To determine the proper number of operating points for linearization, the grade of the nonlinearity of each state variable is evaluated as follows: first, the system is linearized at the equilibrium point x = 0, where the result is A x=0 . Then by randomly changing one state variable x i in the operating range and fix all the other state variables, a series of matrices A xi are obtained. They are compared to A x=0 in terms of the 2-norm of difference. The results are summarized in Table 2. The idea is that if the state variable x i is linear or affine in the system, then identical A matrices should be obtained at different positions of x i . Otherwise, x i should be regarded as nonlinear. From Table 2, it can be seen that the variables θ 2 , . θ 1 and . θ 2 are sources of nonlinearity around the origin, among which θ 2 is the most significant one. The other variables can be treated as linear. Table 2. Inspection of the source of nonlinearity.  Table 2.
In Figure 5, it can be seen that the outputs predicted by the T-S fuzzy models are very close to the testing experimental data. The NDEIs of the neuro-fuzzy model and the two T-S fuzzy models, which are used as the criteria of global approximation, are listed in Table 3. It can be seen that the T-S fuzzy models derived from the neuro-fuzzy model can provide an accurate approximation of the original nonlinear system. The 12 points T-S fuzzy model with optimized operating points gives an even better approximation than the 27 points one with uniformly distributed operating points.

Neuro-fuzzy Model
Because the gear ratio is 5, the command inputs for θ3 and θ4 are given by The command inputs for the angular velocities are   The performance index is defined as: Figure 5. Experimental data and predicted outputs from the T-S fuzzy models. To test the controller performance, the reference commands to be tracked are given by Because the gear ratio is 5, the command inputs for θ 3 and θ 4 are given by The command inputs for the angular velocities are The performance index is defined as: The fuzzy optimal control systems derived from the T-S fuzzy models with 27 and 12 points were simulated using the mathematical model in MATLAB with initial condition x = 0. To obtain an assessment of the proposed control algorithm, it is necessary to compare this paradigm to alternative control techniques. Therefore, a PID controller was also designed: where e 1 and e 2 are tracking error of θ 1 and θ 2 . The PID gains were based on [32], in which a PID controller was designed for this robot manipulator. Figure 6 shows the simulated responses of the flexible robot system with the fuzzy controllers and PID controller, while the mean square errors (MSEs) of command tracking are compared in Table 4. The performances of the fuzzy controllers are evidently better than that of the PID controller. The PID responses exhibit larger overshoot and longer settling time, which leads to larger MSE. The results of the 12 points T-S fuzzy model are almost the same as the results of the 27 points one, as the non-significant difference of modeling errors may only yield trivial effect on the control system. Since the number of linearization is reduced, the 12 points T-S fuzzy model is more computationally efficient.
where e1 and e2 are tracking error of θ1 and θ2. The PID gains were based on [32], in w a PID controller was designed for this robot manipulator. Figure 6 shows the simul responses of the flexible robot system with the fuzzy controllers and PID controller, w the mean square errors (MSEs) of command tracking are compared in Table 4. The formances of the fuzzy controllers are evidently better than that of the PID controller. PID responses exhibit larger overshoot and longer settling time, which leads to la MSE. The results of the 12 points T-S fuzzy model are almost the same as the results o 27 points one, as the non-significant difference of modeling errors may only yield tr effect on the control system. Since the number of linearization is reduced, the 12 poin S fuzzy model is more computationally efficient.   The fuzzy optimal controllers were then implemented on the robot. The experimental results are shown in Figure 7. The achieved tracking performances are good and close to the simulation results, except for some lags and oscillations along the trajectories.  The fuzzy optimal controllers were then implemented on the robot. The experime results are shown in Figure 7. The achieved tracking performances are good and clo the simulation results, except for some lags and oscillations along the trajectories. T oscillations are caused by the flexible joints of the robot, which are very difficult to c pletely eliminate.

Conclusions
An explicit procedure of establishing the T-S fuzzy models of unknown nonli systems from experimental data is presented: (1) the local constituent models are obta by direct linearization of neuro-fuzzy models trained from data so that they are close proximations to the local linearization of the original nonlinear systems; (2) The opera points for linearization are optimized using the evolutionary strategy to achieve g global approximation with a reduced number of linearization. The controller design b on T-S fuzzy model has also been discussed. The derived fuzzy optimal controller applied to a nonlinear flexible-joint robot system and compared with the alternative trol technique, which demonstrated the effectiveness of the proposed method for con ling nonlinear dynamic systems whose analytical models are not available. Comp with the existing T-S model methods, the proposed method could more effectively dress the tradeoff between local and global approximation for complex systems. Mo ver, since the computational methods of the proposed method are based on a set of w developed tools (e.g., FBFN training, ES optimization), it can be readily used for d driven system identification and controller design in realistic applications. The fu work might include studying the effect of uncertainty in data and the modeling erro the neuro-fuzzy and T-S fuzzy models, such that the robustness of the T-S fuzzy con scheme can be more explicitly addressed.

Conclusions
An explicit procedure of establishing the T-S fuzzy models of unknown nonlinear systems from experimental data is presented: (1) the local constituent models are obtained by direct linearization of neuro-fuzzy models trained from data so that they are close approximations to the local linearization of the original nonlinear systems; (2) The operating points for linearization are optimized using the evolutionary strategy to achieve good global approximation with a reduced number of linearization. The controller design based on T-S fuzzy model has also been discussed. The derived fuzzy optimal controller was applied to a nonlinear flexible-joint robot system and compared with the alternative control technique, which demonstrated the effectiveness of the proposed method for controlling nonlinear dynamic systems whose analytical models are not available. Compared with the existing T-S model methods, the proposed method could more effectively address the tradeoff between local and global approximation for complex systems. Moreover, since the computational methods of the proposed method are based on a set of welldeveloped tools (e.g., FBFN training, ES optimization), it can be readily used for datadriven system identification and controller design in realistic applications. The future work might include studying the effect of uncertainty in data and the modeling errors of the neuro-fuzzy and T-S fuzzy models, such that the robustness of the T-S fuzzy control scheme can be more explicitly addressed.