Feedback Linearization for a Generalized Multivariable T-S Model

Abstract

This study presents a novel optimal fuzzy logic control (FLC) strategy based on feedback linearization for the regulation of multivariable nonlinear systems. Building upon an enhanced Takagi–Sugeno (T-S) model previously developed by the authors, the proposed method incorporates a refined parameter-weighting scheme to optimize both local and global approximations within the T-S framework. This approach enables improved selection and minimization of the performance index. The effectiveness of the control strategy is validated through its application to a two-link serial robotic manipulator. The results demonstrate that the proposed FLC achieves robust performance, maintaining system stability and high accuracy even under the influence of noise and load disturbances, with well-damped system behavior and negligible steady-state error.

1. Introduction

Controlling complex systems remains one of the most significant challenges and achievements in modern control theory. Particularly, nonlinear systems, multiple-input multiple-output (MIMO) systems, and systems with significant uncertainties are common in real-world applications yet notoriously difficult to control using conventional techniques. To address these challenges, ongoing research continues to explore and refine analysis tools and control strategies that offer systematic solutions for such complex problems.

This paper introduces a novel controller designed to manage systems exhibiting all three of these challenging characteristics: nonlinearity, multiple variables, and uncertainty. The proposed approach integrates feedback linearization with fuzzy modeling, offering a robust and optimal control solution capable of handling a wide range of operational conditions.

Feedback linearization has emerged as a powerful nonlinear control strategy, gaining widespread adoption over the past few decades. Its key advantage lies in transforming a nonlinear open-loop system into a linearized closed-loop system, typically with a broader effective operating region compared to classical Jacobian-based linearization around an equilibrium point [1]. This method has demonstrated success across diverse domains: from classic pendulum systems [2] and power electronics [3] to highly nonlinear and uncertain platforms such as space payloads [4] and robotic manipulators [5].

Though not universally applicable, feedback linearization provides a flexible framework that accommodates advanced control techniques, including sliding mode control (SMC) [6], and supports integration with learning-based methods, such as neural network-enhanced linearization for systems with dynamic uncertainties [7].

However, even the most sophisticated control strategies rely heavily on accurate system modeling. In this context, fuzzy logic has proven especially effective for modeling complex, nonlinear systems. By encoding expert knowledge or empirical data into rule-based representations, fuzzy modeling enables the derivation of structured system models suitable for control design [8].

Among the various fuzzy modeling approaches, Takagi–Sugeno (T-S) models [9] have become particularly prominent (e.g., see [10,11,12,13]). These models offer powerful approximation capabilities and support rigorous stability analysis and control design, even for MIMO systems [14,15,16]. Their architecture, based on locally affine models blended through fuzzy membership functions (MFs), facilitates the use of data-driven identification techniques, enhancing their adaptability to real-world systems.

Recent research has continued to enhance the accuracy and computational efficiency of T-S models. Techniques aimed at reducing identification complexity, improving convergence speed, and ensuring robustness have been explored in multiple studies [17,18]. Other contributions have focused on theoretical aspects, such as the pros and cons of different kinds of models [19]; stability guarantees; and convergence criteria for fuzzy systems [20,21]. These developments help extend the applicability of T-S models to higher-order systems with more inputs and outputs [22].

A particularly relevant application domain is that of robotic systems, which are classic examples of nonlinear systems with parameter uncertainties and time-varying external disturbances. T-S fuzzy models have proven highly suitable for modeling and controlling these systems [23,24]. Their robustness under challenging operating conditions, such as sudden load changes or sensor noise, has been demonstrated in a variety of practical studies [25,26].

Furthermore, fuzzy FLC is particularly advantageous for enabling interaction with the robot’s environment, especially in unstructured settings where the nature of obstacles is dynamic and often unknown. Notable applications include goods transportation in storage facilities [27], as well as scenarios involving physical interaction with the environment, such as during collisions or when connecting to infrastructure [28].

This article contributes to an ongoing line of research by the authors, building on prior work ([29,30,31]) that proposed improved methodologies for estimating local and global T-S models. In particular, these efforts addressed limitations of traditional identification methods, such as restrictions on non-overlapping membership functions, through a novel weighting parameter approach that reduces the computational cost without sacrificing accuracy.

The current study extends these contributions by enhancing the Fuzzy Logic Controller with Linear Quadratic Regulator (FLC-LQR) architecture introduced in [29]. The new controller includes a time-based control action for trajectory convergence and introduces mechanisms to maintain performance even in the presence of modeling errors and external disturbances.

The integration of Fuzzy Logic Controllers (FLCs) with feedback linearization techniques has also gained attention in recent studies. Several studies have proposed hybrid methods to leverage the strengths of both paradigms (e.g., [32,33,34]). For instance, in [35] the authors present a fuzzy feedback linearization control for nonlinear T-S systems. The main difference from the proposed work is that in [35], the designed algorithm was applied to a single-input–single-output system, while in this work, the proposed controller is applied to multivariable systems. Another difference is that the robustness of the proposed controller is shown by subjecting it to disturbance and noise effects.

The main contributions of this paper are as follows:

The development of an optimal feedback linearization-based FLC, grounded in the authors’ improved T-S modeling framework and integrated with Linear Quadratic Regulator (LQR) optimal control.

Validation of the proposed controller through implementation on a two-link robotic manipulator, demonstrating robust performance, zero steady-state error, and effective disturbance rejection.

The structure of this paper is as follows: Section 2 reviews the fundamentals of feedback linearization. Section 3 introduces a control algorithm, ensuring zero steady-state error. Section 4 details the proposed optimal FLC-LQR methodology. Section 5 discusses fuzzy model identification and estimation, including the use of parameter weighting and constraints of the T-S framework. Section 6 outlines the implementation procedure. Section 7 presents the experimental validation on a two-link robot and compares the proposed method with conventional T-S-based approaches. Section 8 concludes the paper with final remarks and perspectives for future work.

2. Feedback Linearization

Feedback linearization is a widely adopted methodology in nonlinear control system design that has seen significant development over the past few decades. Its core principle lies in transforming a nonlinear system, either partially or fully, into a linear one through a carefully designed feedback mechanism. Unlike conventional Jacobian linearization, which locally approximates nonlinear behavior around an equilibrium point, feedback linearization achieves global or semi-global linearization by explicitly canceling out nonlinear dynamics through state and input transformations. This process enables the direct application of linear control techniques to the resulting linearized closed-loop system.

In this section, we introduce the foundational concepts and formal procedures involved in feedback linearization. In the following section, this methodology will be further extended to include an integral action, aiming to eliminate steady-state errors in the controlled system. However, it is worth noting that practical implementations of feedback linearization often require supplementary mechanisms (such as anti-windup strategies and signal filtering) to account for real-world issues not captured by the theoretical model, including actuator saturation and measurement noise.

We start by modeling the system to be controlled as an affine system in its canonical form. It is modeled as follows:

(1)$\begin{array}{ccc}\hfill \dot{x}\left(t\right)& =& A_0\left(x\left(t\right)\right)+A\left(x\left(t\right)\right)x\left(t\right)+B\left(x\left(t\right)\right)u\left(t\right)\hfill \\ \hfill y\left(t\right)& =& Cx\left(t\right)\hfill \end{array}$

where $\dot{x}\left(t\right)$ denotes the time derivative of the state vector $x\left(t\right)$; $u\left(t\right)$ is the system input; $y\left(t\right)$ is the output; $A_0\left(x\left(t\right)\right)$ represents the drift dynamics or the nonlinear bias term; $A\left(x\right(t\left)\right)$ and $B\left(x\right(t\left)\right)$ are state-dependant matrices defining the linear and input coupling terms, respectively; and C is the output matrix, assumed to be constant in this formulation. The second equation dictates the output of the system, which in this particular model is only dependent on the state vector based on the relation matrix C. The dimensions of the variables and functions involved are specified as

$\begin{array}{ccc}\hfill x& :& \Re ⟶{\Re }^n\hfill \\ \hfill y& :& \Re ⟶{\Re }^q\hfill \\ \hfill u& :& \Re ⟶{\Re }^m\hfill \\ \hfill A_0& :& {\Re }^n⟶{\Re }^n\hfill \\ \hfill A& :& {\Re }^n⟶{\Re }^{n\times n}\hfill \\ \hfill B& :& {\Re }^n⟶{\Re }^{n\times m}\hfill \\ \hfill C& :& {\Re }^{q\times n}\hfill \end{array}$

We assume that the system satisfies the condition $n\ge m$, ensuring that the number of state variables is at least as large as the number of inputs. This forms the basis upon which the feedback linearization procedure is constructed in the subsequent sections.

The canonical form of the system’s matrices has the following structure:

$A_0=\left[\begin{array}{c}{A}_{01}\\ A_{02}\\ ⋮\\ A_{0m}\end{array}\right],A_{0i}=\left[\begin{array}{c}0\\ ⋮\\ 0\\ a_{0i}\end{array}\right]\in {\Re }^{n_i}$

$A=\left[\begin{array}{cccc}{A}_{11}& A_{12}& \dots & A_{1m}\\ A_{21}& A_{22}& \dots & A_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ A_{m1}& A_{m2}& \dots & A_{mm}\end{array}\right]$

$A_{ii}=\left[\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 1\\ a_{ii1}& a_{ii2}& a_{ii3}& \dots & a_{ii{n}_i}\end{array}\right]\in {\Re }^{n_i\times n_i}$

$A_{ij}=\left[\begin{array}{ccccc}0& 0& 0& \dots & 0\\ 0& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 0\\ a_{ij1}& a_{ij2}& a_{ij3}& \dots & a_{ij{n}_j}\end{array}\right]\in {\Re }^{n_i\times n_j}$

$B=\left[\begin{array}{c}{B}_1\\ B_2\\ ⋮\\ B_m\end{array}\right],B_i=\left[\begin{array}{cccc}0& 0& \dots & 0\\ ⋮& ⋮& \dots & ⋮\\ 0& 0& \dots & 0\\ b_{i1}& b_{i2}& \dots & b_{im}\end{array}\right]\in {\Re }^{n_i\times m}$

where $\left\{n_1,n_2,\cdots ,n_m\right\}$ denote the dimensions of the controllable subspaces associated with each input.

Our primary objective is to reach the desired final zero state, regardless of the initial state. To this end, we define the following matrix:

$S=\left[\begin{array}{cccc}{S}_{11}& S_{12}& \dots & S_{1m}\\ S_{21}& S_{22}& \dots & S_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ S_{m1}& S_{m2}& \dots & S_{mm}\end{array}\right]\in {\Re }^{m\times n}$

$S_{ii}=\left[\begin{array}{ccccc}0& 0& \dots & 0& 1_{n_i}\end{array}\right]\in {\Re }^{1\times n_i}$

$S_{ij}=\left[\begin{array}{ccccc}0& 0& \dots & 0& 0\end{array}\right]\in {\Re }^{1\times n_j}$

If the feedback system is required to have a system matrix $A_f$ with the same structure as $A\left(x\right(t\left)\right)$, then the control action must be

(2)$u\left(t\right)=-[SB{\left(x\left(t\right)\right]}^{-1}S[A_0\left(x\left(t\right)\right)+(A\left(x\left(t\right)\right)-A_f)x\left(t\right)]$

assuming that the matrix $SB\left(x\right(t\left)\right)$ is non-singular throughout the entire operational workspace. Under this assumption, the feedback system becomes

(3)$\dot{x}\left(t\right)=A_0\left(x\left(t\right)\right)+A\left(x\left(t\right)\right)x\left(t\right)-B\left(x\left(t\right)\right){\left[SB\left(x\left(t\right)\right)\right]}^{-1}S[A_0\left(x\left(t\right)\right)+(A\left(x\left(t\right)\right)-A_f)x\left(t\right)]$

Note that

$SB=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1m}\\ b_{21}& b_{22}& \dots & b_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ b_{m1}& b_{m2}& \dots & b_{mm}\end{array}\right]\in {\Re }^{m\times m}$

and therefore,

$B{\left[SB\right]}^{-1}=\left[\begin{array}{c}{B}_1^{*}\\ B_2^{*}\\ ⋮\\ B_m^{*}\end{array}\right]\in {\Re }^{n\times m}$

$\begin{array}{cc}{B}_i^{*}& =B_i{\left[SB\right]}^{-1}=\left[\begin{array}{cccc}0& 0& \dots & 0\\ ⋮& ⋮& \dots & ⋮\\ 0& 0& \dots & 0\\ b_{i1}& b_{i2}& \dots & b_{im}\end{array}\right]{\left[\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1m}\\ b_{21}& b_{22}& \dots & b_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ b_{m1}& b_{m2}& \dots & b_{mm}\end{array}\right]}^{-1}\\ & =\left[\begin{array}{cccccc}0& \dots & 0& 0& \dots & 0\\ ⋮& \dots & ⋮& ⋮& \dots & ⋮\\ 0& \dots & 0& 0& \dots & 0\\ 0& \dots & 1_{n_{i}i}& 0& \dots & 0\end{array}\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}$

This means that

$B{\left[SB\right]}^{-1}=S^T$

fulfilling

$B{\left[SB\right]}^{-1}S=S^{T}S\in {\Re }^{n\times n}$

$S^{T}S=\left[\begin{array}{cccc}{\tilde{S}}_{11}& {\tilde{S}}_{12}& \dots & {\tilde{S}}_{1m}\\ {\tilde{S}}_{21}& {\tilde{S}}_{22}& \dots & {\tilde{S}}_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{S}}_{m1}& {\tilde{S}}_{m2}& \dots & {\tilde{S}}_{mm}\end{array}\right]$

$\begin{array}{c}\hfill {\tilde{S}}_{ii}=\left[\begin{array}{cccc}0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮\\ 0& \dots & 0& 0\\ 0& \dots & 0& 1_{n_i{n}_i}\end{array}\right]\in {\Re }^{n_i\times n_i},S_{ij}=0\in {\Re }^{n_i\times n_j}\end{array}$

Therefore, it will satisfy the following:

$B{\left[SB\right]}^{-1}S{A}_0=A_0$

$A-B{\left[SB\right]}^{-1}S[A-A_f]=A_f$

Consequently, the system adopts the following structure:

(4)$\dot{x}\left(t\right)=A_{f}x\left(t\right)$

3. Steady State Error

We now turn our attention to the problem of ensuring that, at the steady state, the system output variables $y_p$ reach their set points. A conventional approach to achieve zero steady-state error is to incorporate integral action, which is equivalent to introducing a set of q additional states defined as follows:

$\dot{w}\left(t\right)=y_p-y\left(t\right)=y_p-Cx\left(t\right)w:\Re \to {\Re }^q$

We define the following extended state vector:

(5)$\begin{array}{c}\hfill x_e\left(t\right)=\left[\begin{array}{c}w\left(t\right)\\ x\left(t\right)\end{array}\right]\end{array}$

Consequently, $w\left(t\right)$ represents the integral of the position error, and the complete system is given by

(6)$\begin{array}{ccc}\hfill {\dot{x}}_e\left(t\right)& =& \left[\begin{array}{c}0\\ A_0\left(x\left(t\right)\right)\end{array}\right]+\left[\begin{array}{cc}0& -C\\ 0& A\left(x\right(t\left)\right)\end{array}\right]x_e\left(t\right)+\left[\begin{array}{c}0\\ B\left(x\right(t\left)\right)\end{array}\right]u\left(t\right)+\left[\begin{array}{c}I\\ 0\end{array}\right]y_p\hfill \\ \hfill y\left(t\right)& =& \left[\begin{array}{cc}0& C\end{array}\right]x_e\left(t\right)\hfill \end{array}$

For this augmented system, the following control action is proposed:

(7)$u\left(t\right)=-[SB{\left(x\left(t\right)\right]}^{-1}S[A_0\left(x\left(t\right)\right)+(A\left(x\left(t\right)\right)-A_f)x\left(t\right)-Kw\left(t\right)]$

where

$K=\left[\begin{array}{c}{K}_1\\ K_2\\ ⋮\\ K_m\end{array}\right],K_i=\left[\begin{array}{ccc}0& 0& 0\\ ⋮& ⋮& ⋮\\ 0& 0& 0\\ K_{i1}& K_{i2}& K_{iq}\end{array}\right]\in R^{n_i\times q}$

It can be verified that

(8)$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =A_0\left(x\left(t\right)\right)+A\left(x\left(t\right)\right)x\left(t\right)-B\left(x\left(t\right)\right){\left[SB\left(x\left(t\right)\right)\right]}^{-1}S[A_0\left(x\left(t\right)\right)+(A\left(x\left(t\right)\right)-A_f)x\left(t\right)]\hfill \\ & +B\left(x\left(t\right)\right){\left[SB\left(x\left(t\right)\right)\right]}^{-1}SKw\left(t\right)\hfill \\ \hfill \dot{x}\left(t\right)& =A_{f}x\left(t\right)+Kw\left(t\right)\hfill \end{array}$

and therefore, the feedback system takes the following form:

(9)$\begin{array}{c}{\dot{x}}_e\left(t\right)=\left[\begin{array}{cc}0& -C\\ K& A_f\end{array}\right]x_e\left(t\right)+\left[\begin{array}{c}I\\ 0\end{array}\right]y_p\hfill \\ y\left(t\right)=\left[\begin{array}{cc}0& C\end{array}\right]x_e\left(t\right)\hfill \end{array}$

The dynamics of this system are governed by the new system matrix, and if it stable, the system will reach a steady state:

(10)$\begin{array}{cc}\hfill x_s=\underset{t\to \infty }{lim}{x}_e\left(t\right)& =\left[\begin{array}{c}{w}_s\\ x_s\end{array}\right]\Rightarrow x_e^{\prime }=0\\ \hfill y_s=C{x}_s\end{array}$

It can be verified that

(11)$\begin{array}{c}\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{cc}0& -C\\ K& A_f\end{array}\right]\left[\begin{array}{c}{w}_s\\ x_s\end{array}\right]+\left[\begin{array}{c}I\\ 0\end{array}\right]y_p\Rightarrow \begin{array}{c}0=-C{x}_s+y_p\hfill \\ y_s=\left[\begin{array}{cc}0& C\end{array}\right]x_e\left(t\right)\Rightarrow y_s=C{x}_s=y_p\hfill \end{array}\hfill \end{array}$

This implies that the control objective of achieving zero steady-state error has been met.

4. Optimal Controller Design

To select the appropriate control matrices $A_f$ and K, this work proposes designing these matrices based on a representative central operating point of the system:

(12)$\begin{array}{ccc}\hfill \dot{x}\left(t\right)& =& A_{c}x\left(t\right)+B_{c}u\left(t\right)\hfill \\ \hfill y\left(t\right)& =& Cx\left(t\right)\hfill \end{array}$

Therefore, the extended system incorporating integral action is given by

(13)$\begin{array}{ccc}\hfill {\dot{x}}_e\left(t\right)& =& \left[\begin{array}{cc}0& -C\\ 0& A_c\end{array}\right]{\dot{x}}_e\left(t\right)+\left[\begin{array}{c}0\\ B_c\end{array}\right]u\left(t\right)+\left[\begin{array}{c}I\\ 0\end{array}\right]y_p\hfill \\ \hfill y\left(t\right)& =& \left[\begin{array}{cc}0& C\end{array}\right]{\dot{x}}_e\left(t\right)\hfill \end{array}$

Assuming that $y_p=0$, it follows that

(14)$\begin{array}{c}{\dot{x}}_e\left(t\right)=\left[\begin{array}{cc}0& -C\\ 0& A_c\end{array}\right]{\dot{x}}_e\left(t\right)+\left[\begin{array}{c}0\\ B_c\end{array}\right]u\left(t\right)\hfill \end{array}$

We now seek the control action $u\left(t\right)$ that drives the state from an arbitrary initial condition $x_e\left(t_0\right)$ to the desired final zero state. To achieve this, we employ an LQR to minimize the quadratic performance criterion:

(15)$J={\int }_{t_0}^{\infty }(x_e^t\left(t\right)Q{x}_e\left(t\right)+u^{t}Ru)dt$

where $Q\in {\Re }^{(n+q)\times (n+q)}$ and $R\in {\Re }^{m\times m}$ are symmetric positive definite matrices. The optimal control law is determined by the following equation:

(16)$\begin{array}{c}\hfill u\left(t\right)=-\left[\begin{array}{cc}{K}_w& K_x\end{array}\right]x_e\left(t\right)\end{array}$

(17)$\begin{array}{c}\hfill \left[\begin{array}{cc}{K}_w& K_x\end{array}\right]=R^{-1}{B}_c^{t}L\end{array}$

where $L\in {\Re }^{n\times n}$ is the solution to the following Riccati equation:

(18)$0=-Q+L{B}_c{R}^{-1}{B}_c^{t}L-L{A}_c-A_c^{t}L$

Thus, the closed-loop systems will behave as

(19)$\begin{array}{cc}\hfill {\dot{x}}_e\left(t\right)& =\left[\begin{array}{cc}0& -C\\ 0& A_c\end{array}\right]x_e\left(t\right)+\left[\begin{array}{c}0\\ B_c\end{array}\right]u\left(t\right)\hfill \\ & =\left[\begin{array}{cc}0& -C\\ K& A_c\end{array}\right]x_e\left(t\right)-\left[\begin{array}{c}0\\ B_c\end{array}\right]\left[\begin{array}{cc}{K}_w& K_x\end{array}\right]x_e\left(t\right)\hfill \end{array}$

(20)$\begin{array}{ccc}\hfill {\dot{x}}_e\left(t\right)& =& \left[\begin{array}{cc}0& -C\\ -B_c{K}_w& A_c-B_c{K}_x\end{array}\right]x_e\left(t\right)\hfill \end{array}$

Therefore, the goal is for the feedback linearized system to emulate the optimal system described in (20). In other words, the desired feedback system matrix is proposed as

(21)$K=-B_c{K}_w$

(22)$A_f=A_c-B_c{K}_x$

where K is the feedback matrix associated with the integral position error $w\left(t\right)$ and $A_f$ is the desired feedback system matrix.

This approach ensures optimal dynamics at a representative central operating point.

5. Fuzzy T-S Model Parameter Estimation

To obtain the system matrices given in (1), we propose using the T-S identification method [9], which estimates system parameters by minimizing a performance index. An enhancement to the traditional T-S modeling approach is presented by the authors in [29].

The proposed method employs functions such as the following for identification:

$f:{\Re }^n⟶\Re$

$y=f(x_1,x_1,\dots ,x_n)$

These functions can be transformed into if–then rules for an $n^{th}$-order system:

$\begin{array}{c}\hfill S^{(i_1\dots i_m)}:If{z}_{1}is{M}_1^{i_1}and \dots z_{m}is{M}_m^{i_m}then\widehat{y}=p_0^{(i_1\dots i_m)}+p_1^{(i_1\dots i_m)}{x}_1+\dots +p_n^{(i_1\dots i_m)}{x}_n\end{array}$

where $\left\{y, x_1, x_2, \dots , x_n\right\}$ are measurable variables, and $\left\{z_1, z_2, \dots , z_m\right\}$ is a subset of $\left\{x_1, x_2, \dots , x_n\right\}$.

The fuzzy estimation of the output is given by

(23)$\begin{array}{c}\hfill \widehat{y}=\sum _{i_1=1}^{r_1}\dots \sum _{i_m=1}^{r_m}{\beta }^{(i_1\dots i_m)}\left(z\right)\left[p_0^{(i_1\dots i_m)}+p_1^{(i_1\dots i_m)}{x}_1+\dots +p_n^{(i_1\dots i_m)}{x}_n\right]\end{array}$

where

(24)$\begin{array}{c}\hfill {\beta }_k^{(i_1\dots i_m)}\left(z\right)={\displaystyle \frac{{\mu }_{1i_1}\left(z_1\right){\mu }_{2i_2}\left(z_2\right)\dots {\mu }_{m{i}_m}\left(z_m\right)}{{\sum }_{j_1=1}^{r_1}\dots {\sum }_{j_m=1}^{r_m}{\mu }_{1j_1}\left(z_1\right){\mu }_{2j_2}\left(z_2\right)\dots {\mu }_{m{j}_m}\left(z_m\right)}}\end{array}$

and ${\mu }_{j{i}_j}\left(z_j\right)$ denotes the membership function corresponding to the set $M_j^{i_j}$.

Fuzzy system parameters can be obtained by minimizing the following quadratic performance index, where the series $\left\{x_{1k},x_{2k},\dots ,x_{nk},y_k\right\}$ represents a set of N input/output samples:

(25)$J=\sum _{k=1}^N{(y_k-{\widehat{y}}_k)}^2={∥Y-XP∥}^2$

with

$\begin{array}{ccc}\hfill Y& =& {\left[\begin{array}{cccc}{y}_1\hfill & y_2\hfill & \dots \hfill & y_N\hfill \end{array}\right]}^T\hfill \\ \hfill P& =& {\left[\begin{array}{ccc}{p}_0^{(1\dots 1)}{p}_1^{(1\dots 1)}{p}_2^{(1\dots 1)}\hfill & \dots \hfill & p_n^{(1\dots 1)}\dots p_0^{(r_1\dots r_m)}\dots p_n^{(r_1\dots r_m)}\hfill \end{array}\right]}^T\hfill \\ \hfill X& =& \left[\begin{array}{ccc}{\beta }_1^{(1\dots 1)}{\beta }_1^{(1\dots 1)}{x}_{11}\dots {\beta }_1^{(1\dots 1)}{x}_{n1}\hfill & \dots \hfill & {\beta }_1^{(r_1\dots r_m)}\dots {\beta }_1^{(r_1\dots r_m)}{x}_{n1}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}⋮\hfill \\ {\beta }_N^{(1\dots 1)}{\beta }_N^{(1\dots 1)}{x}_{1N}\dots {\beta }_N^{(1\dots 1)}{x}_{nN}\hfill & \dots \hfill & {\beta }_N^{(r_1\dots r_N)}\dots {\beta }_N^{(r_1\dots r_m)}{x}_{nN}\hfill \end{array}\right]\hfill \end{array}$

If X is a full-rank matrix, a solution can be obtained directly as

(26)$\begin{array}{ccc}\hfill J& =& {∥Y-XP∥}^2={(Y-XP)}^T(Y-XP)\hfill \\ \hfill \nabla J& =& X^T(Y-XP)=X^{T}Y-X^{T}XP=0\hfill \\ \hfill P& =& {\left(X^{T}X\right)}^{-1}{X}^{T}Y\hfill \end{array}$

A limitation of the methodology presented in [9] is that it cannot be applied to the most common type of MFs, such as those illustrated in Figure 1.

The membership functions

${\mu }_{i1}\left(v_i\right)={\displaystyle \frac{b_i-v_i}{b_i-a_i}}\mathrm{and}{\mu }_{i2}\left(v_i\right)={\displaystyle \frac{v_i-a_i}{b_i-a_i}}$

are defined over the interval $[a_i,b_i]$, which must satisfy the following condition:

$\begin{array}{cc}{\mu }_{i1}\left(a_i\right)=1\hfill & {\mu }_{i1}\left(b_i\right)=0\hfill \\ {\mu }_{i2}\left(a_i\right)=0\hfill & {\mu }_{i2}\left(b_i\right)=1\hfill \\ {\mu }_{i1}\left(v_i\right)+{\mu }_{i2}\left(v_i\right)=1\hfill \end{array}$

It has been demonstrated [29,30,31] that, in this common case, the matrix X is rank-deficient. Consequently, $X^{T}X$ is not invertible, preventing the direct application of the aforementioned method. A detailed proof can be found in [31].

When the matrix X is not of full rank, an efficient approach with low computational complexity, based on the well-known parameter-weighting method, is presented in [31].

Another advantage of this method is its capability to tune the T-S model parameters using experimental data. To achieve this, we begin with an initial parameter estimation as follows:

(27)$p_0={[p_0^0{p}_1^0{p}_2^0\dots p_n^0]}^T$

which may be obtained, for example, through a preliminary least squares analysis of the system under known operating conditions. This initial estimate serves as a reference in the calculation of the final system parameters. The fuzzy model parameters are then determined by minimizing the following function:

(28)$\begin{array}{cc}\hfill J& =\sum _{k=1}^N{(y_k-{\widehat{y}}_k)}^2+{\gamma }^2\sum _{i_1=1}^{r_1}\dots \sum _{i_m=1}^{r_m}\sum _{j=0}^m{(p_j^0-p_j^{(i_1\dots i_m)})}^2\hfill \\ & ={∥Y-XP∥}^2+{\gamma }^2{∥P_0-P∥}^2\hfill \\ & ={∥\left[\begin{array}{c}Y\hfill \\ \gamma P_0\hfill \end{array}\right]-\left[\begin{array}{c}X\hfill \\ \gamma I\hfill \end{array}\right]P∥}^2\hfill \\ & ={∥Y_a-X_{a}P∥}^2\hfill \end{array}$

where

$P_0=\underset{r_1.r_2\dots r_m}{\underbrace{{[p_0{p}_0\dots p_0]}^T}}$

where $\gamma$ is a parameter representing the confidence level of the initial estimation [29,31].

When applied to a multivariable system, it is composed of a vector-valued function of the following form:

$f:{\Re }^n⟶{\Re }^q$

$y=f(x_1,x_2,\dots ,x_n)y\in {\Re }^q$

or, equivalently, as a set of q scalar functions

$y_j=f_j(x_1,x_2,\dots ,x_n)y\in \Re ,j=1\dots q$

For each of these functions, a MF can be modeled as

$\begin{array}{c}\hfill S_j^{(i_1\dots i_m)}:If{z}_{1}is{M}_{j1}^{i_1}and{z}_{2}is{M}_{j2}^{i_2}and\dots z_{m}is{M}_{jm}^{i_m}then\\ \hfill {\widehat{y}}_j=p_{j0}^{(i_1\dots i_m)}+p_{j1}^{(i_1\dots i_m)}{x}_1+p_{j2}^{(i_1\dots i_m)}{x}_2+\dots +p_{jn}^{(i_1\dots i_m)}{x}_n\end{array}$

When applied to the entire model, the basic system takes the following form:

$\begin{array}{ccc}\hfill Y_N& =& {\left[\begin{array}{cccc}{y}_{j1}\hfill & y_{j2}\hfill & \dots \hfill & y_{jN}\hfill \end{array}\right]}^T\hfill \\ \hfill P_j& =& {\left[\begin{array}{ccc}{p}_{j0}^{(1\dots 1)}{p}_{j1}^{(1\dots 1)}{p}_{j2}^{(1\dots 1)}\hfill & \dots \hfill & p_{jn}^{(1\dots 1)}\dots p_{j0}^{(r_1\dots r_m)}\dots p_{jn}^{(r_1\dots r_m)}\hfill \end{array}\right]}^T\hfill \\ \hfill X_j& =& \left[\begin{array}{ccc}{\beta }_{j1}^{(1\dots 1)}{\beta }_{j1}^{(1\dots 1)}{x}_{11}\dots {\beta }_{j1}^{(1\dots 1)}{x}_{n1}\hfill & \dots \hfill & {\beta }_{j1}^{(r_1\dots r_m)}\dots {\beta }_{j1}^{(r_1\dots r_m)}{x}_{n1}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}⋮\hfill \\ {\beta }_{jN}^{1\dots 1}{\beta }_{jN}^{(1\dots 1)}{x}_{1N}\dots {\beta }_{jN}^{(1\dots 1)}{x}_{nN}\hfill & \dots \hfill & {\beta }_{jN}^{(r_1\dots r_m)}\dots {\beta }_{jN}^{(r_1\dots r_m)}{x}_{nN}\hfill \end{array}\right]\hfill \end{array}$

accompanied by the tuning equation:

(29)$\begin{array}{ccc}\hfill J_j& =& {∥\left[\begin{array}{c}{Y}_j\hfill \\ {\mathrm{\Gamma }}_j{P}_{j0}\hfill \end{array}\right]-\left[\begin{array}{c}{X}_j\hfill \\ {\mathrm{\Gamma }}_j\hfill \end{array}\right]P_j∥}^2\hfill \end{array}$

Furthermore, it is common for all premises of the system to be identical, implying that in such cases,

$X_1=X_2=\dots =X_q=X$

If, in addition, the weighting factors ${\gamma }_j^{(i_1\dots i_m)}$ for all parameters of $\Gamma$ are equal, then

$J={∥\left[\begin{array}{cccc}{Y}_1\hfill & Y_2\hfill & \dots \hfill & Y_q\hfill \\ \mathrm{\Gamma }{P}_{10}\hfill & \mathrm{\Gamma }{P}_{20}\hfill & \dots \hfill & \mathrm{\Gamma }{P}_{q0}\hfill \end{array}\right]-\left[\begin{array}{c}X\hfill \\ \mathrm{\Gamma }\hfill \end{array}\right]\left[\begin{array}{cccc}{P}_1\hfill & P_2\hfill & \dots \hfill & P_m\hfill \end{array}\right]∥}^2$

$J={∥Y_a-X_{a}P∥}^2$

with the solution remaining as

(30)$P={\left(X_a^t{X}_a\right)}^{-1}{X}_a^t{Y}_a$

A detailed proof is provided in [29].

6. Implementation of the Proposed Method

The proposed control method is summarized in the following steps, which are organized into offline and online tasks.

Offline tasks:

Obtain the T-S model through experimentation: $A_0\left(x\left(t\right)\right)$, $A\left(x\right(t\left)\right)$, $B\left(x\right(t\left)\right)$, and C.

Determine a central working point for the system; this may correspond to the central rule of the T-S model or the $p_0$ parameters from Section 5: $A_c$ and $B_c$.

Compute the desired system matrix and the steady-state error cancellation matrix using LQR: $A_f$ and K.

Determine the matrix S.

Online tasks:

Update the system matrices at the current instant: $A_0\left(x\left(t\right)\right)$, $A\left(x\right(t\left)\right)$, and $B\left(x\right(t\left)\right)$ at the current time.

Update the control action $u\left(t\right)$ at the current time using Equation (7).

7. Example: Two-Link Robot

To illustrate the application of the proposed control methodology, consider a simple two-link robot, as shown in Figure 2. Each articulated joint is driven by a motor equipped with an encoder and a tachometer, providing torque as well as measurements of position and velocity. The system input consists of desired joint positions ${\theta }_1$ and ${\theta }_2$ or joint position trajectories ${\theta }_{d1}\left(t\right)$ and ${\theta }_{d2}\left(t\right)$, which are specified by the robot’s control unit. These conditions are typical in industrial robots tasked with following prescribed paths.

The Lagrangian dynamical model of this type of robot is given by

(31)$\begin{array}{c}\hfill \left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{c}{\ddot{\theta }}_1\\ {\ddot{\theta }}_2\end{array}\right]+\left[\begin{array}{cc}-h{\dot{\theta }}_2+{\beta }_1& -h{\dot{\theta }}_1-h{\dot{\theta }}_2\\ h{\dot{\theta }}_1& {\beta }_2\end{array}\right]\left[\begin{array}{c}{\dot{\theta }}_1\\ {\dot{\theta }}_2\end{array}\right]+\left[\begin{array}{c}{g}_1\\ g_2\end{array}\right]=\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\end{array}\right]\end{array}$

where $\theta ={\left[{\theta }_1{\theta }_2\right]}^T$ denotes the joint angles in radians, and $\tau ={\left[{\tau }_1{\tau }_2\right]}^T$ represents the input torque at each joint in Newton-meters. The terms $H_{ij}$ correspond to the elements of the inertia matrix, defined as

$\begin{array}{ccc}\hfill H_{11}& =& m_1{l}_{c1}^2+I_1+m_2(l_1^2+l_{c2}^2+2l_1{l}_{c2}cos{\theta }_2)+I_2\hfill \\ \hfill H_{12}& =& m_2{l}_1{l}_{c2}cos{\theta }_2+m_2{l}_{c2}^2+I_2\hfill \\ \hfill H_{21}& =& H_{12}\hfill \\ \hfill H_{22}& =& m_2{l}_{c2}^2+I_2\hfill \end{array}$

Here, $m_i$ is the mass of the $i^{th}$ link in kilograms; $l_{ci}$ is the distance from the ${(i-1)}^{th}$ joint to the center of mass of the $i^{th}$ link (whose total length is $l_i$), measured in meters; and $I_i$ denotes the moment of inertia of the $i^{th}$ link, expressed in kilogram–meter squared. The term h is defined as

$\begin{array}{ccc}\hfill h& =& m_2{l}_1{l}_{c2}sin{\theta }_2\hfill \end{array}$

The viscous friction coefficients ${\beta }_1$ and ${\beta }_2$, given in kilogram–meter squared per second, form part of the centripetal and Coriolis force matrix. Finally, $g_1$ and $g_2$ represent the gravity compensation terms for each link, with g = 9.81 m/s² denoting the gravitational acceleration, where

$\begin{array}{ccc}\hfill g_1& =& m_1{l}_{c1}gsin{\theta }_1+m_{2}g(l_{c2}sin({\theta }_1+{\theta }_2)+l_{1}sin{\theta }_1)\hfill \\ \hfill g_2& =& m_2{l}_{c2}gsin({\theta }_1+{\theta }_2)\hfill \end{array}$

In our specific example, we consider a robot with the following characteristics:

$m_1=2\mathrm{kg},l_1=1\mathrm{m},l_{c1}=0.5\mathrm{m},I_1=0.25\mathrm{kg}{\mathrm{m}}^2,{\beta }_1=0.05{\displaystyle \frac{\mathrm{kg}{\mathrm{m}}^2}{s}}$

$m_2=1.5\mathrm{kg},l_2=1\mathrm{m},l_{c2}=0.75\mathrm{m},I_2=0.12\mathrm{kg}{\mathrm{m}}^2,{\beta }_2=0.05{\displaystyle \frac{\mathrm{kg}{\mathrm{m}}^2}{s}}$

For

${\theta }_1={\dot{\theta }}_1={\theta }_2={\dot{\theta }}_2=0$

it can be verified that

${\tau }_{10}=g_{10}=0,{\tau }_{20}=g_{20}=0$

7.1. T-S Model

The states of the model are defined as

$x={\left[\begin{array}{cccc}{\theta }_1& {\dot{\theta }}_1& {\theta }_2& {\dot{\theta }}_2\end{array}\right]}^T$

In this nonlinear system, the input torque pair $[{\tau }_1,{\tau }_2]$ is linearly dependent, implying that the T-S rules do not depend on the control inputs. Suppose we assign three fuzzy sets to each state variable; given that there are four state variables, the complete rule base will consist of $3^4=81$ rules of the following form:

$\begin{array}{c}{S}^{(i_1\dots i_n)}:\mathit{If}{\mathit{x}}_{\mathit{1}}\mathit{is}{\mathit{M}}_{\mathit{1}}^{{\mathit{i}}_{\mathit{1}}}\mathit{and}\dots {\mathit{x}}_{\mathit{4}}\mathit{is}{\mathit{M}}_{\mathit{4}}^{{\mathit{i}}_{\mathit{4}}}\mathit{then}\\ \\ \dot{x}=\left[\begin{array}{c}0\\ a_{20}^{\left(i_1{i}_2{i}_3{i}_4\right)}\\ 0\\ a_{40}^{\left(i_1{i}_2{i}_3{i}_4\right)}\end{array}\right]+\left[\begin{array}{cccc}0& 1& 0& 0\\ a_{21}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{22}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{23}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{24}^{\left(i_1{i}_2{i}_3{i}_4\right)}\\ 0& 0& 0& 1\\ a_{41}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{42}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{43}^{\left(i_1{i}_2{i}_3{i}_4\right)}& a_{44}^{\left(i_1{i}_2{i}_3{i}_4\right)}\end{array}\right]x\\ +\left[\begin{array}{cc}0& 0\\ b_{21}^{\left(i_1{i}_2{i}_3{i}_4\right)}& b_{22}^{\left(i_1{i}_2{i}_3{i}_4\right)}\\ 0& 0\\ b_{41}^{\left(i_1{i}_2{i}_3{i}_4\right)}& b_{42}^{\left(i_1{i}_2{i}_3{i}_4\right)}\end{array}\right]\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\end{array}\right]\end{array}$

The workspace is defined as follows: ${\theta }_1\in [0,\pi ],{\dot{\theta }}_1\in [-1,1],{\theta }_2\in [-{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}}],\mathrm{and}{\dot{\theta }}_2\in [-1,1]$. For this example, we employ five membership functions for the position variables and three for the velocity variables. Figure 3, Figure 4, Figure 5 and Figure 6 display the membership functions for $\theta ,{\dot{\theta }}_1,{\theta }_2\mathtt{and}{\dot{\theta }}_2$, respectively.

By applying the methodology described in [29,31], we obtain the following results:

$\begin{array}{c}{S}^{1111}:If{x}_{1}is{M}_1^{i_1}and\dots x_{4}is{M}_4^{i_1}then\\ {\dot{x}}_2=-2.2816-14.8609x_1+1.2924x_2-2.1967x_3+0.1681x_4+0.3355u_1-0.3346u_2\\ {\dot{x}}_4=8.4921+24.2250x_1-3.1540x_2+0.7544x_3+0.4721x_4-0.2624u_1+1.1547u_2\\ ⋮\\ S^{4444}:If{x}_{1}is{M}_1^{i_4}and\dots x_{4}is{M}_4^{i_4}then\\ {\dot{x}}_2=-1.5761+1.7165x_1-1.7015x_2-1.6274x_3+1.1794x_4+0.2618u_1-0.3398u_2\\ {\dot{x}}_4=6.0187-3.8358x_1+4.1864x_2+8.5766x_3-2.1167x_4-0.1041u_1+1.1526u_2\end{array}$

Using this approximation, the mean square error for the first joint is 0.0016, with a maximum absolute error of 1.1055. For the second joint, the mean square error is 0.0040, with the same maximum absolute error of 1.1055.

The parameter-weighting method ensures that the obtained results represent the best possible approximation of the T-S model as the interval coverage approaches completeness. However, this theoretical limit cannot be fully achieved in practice, as it renders the computation intractable. Nevertheless, this method consistently yields more accurate results than the original T-S modeling approach.

The outputs of the system are the joint angles of the first and second links, denoted by ${\theta }_1$ and ${\theta }_2$, respectively, where

(32)$\begin{array}{ccc}\hfill C& =& \left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]\hfill \end{array}$

7.2. Proposed Optimal Controller

The matrix S is defined as follows:

(33)$\begin{array}{ccc}\hfill S& =& \left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]\hfill \end{array}$

First, the feedback system matrix must be computed. The linearized system is obtained by applying the minimum least squares method around the equilibrium point defined by

${\left[\begin{array}{cccc}{\displaystyle \frac{\pi }{2}}& 0& 0& 0\end{array}\right]}^T$

$\left[\begin{array}{c}{\tau }_{10}\\ {\tau }_{20}\end{array}\right]=\left[\begin{array}{c}35.525\\ 11.025\end{array}\right]$

The resulting linear model of the system is given by

(34)$\begin{array}{ccc}\hfill \dot{x}& =& \left[\begin{array}{cccc}0& 1& 0& 0\\ -5.8137& 0.0072& 0.4226& -0.1265\\ 0& 0& 0& 1.0000\\ 7.9809& -0.0089& -3.4335& 0.4205\end{array}\right]x\left(t\right)+\left[\begin{array}{cc}0& 0\\ 0.7778& -1.5877\\ 0& 0\\ -1.5968& 4.1704\end{array}\right]u\left(t\right)\hfill \end{array}$

The selection of the weighting matrices Q and R in the cost function was performed using an iterative trial-and-error approach. This involved systematically varying the elements of the matrices to evaluate their influence on the system’s performance. A commonly adopted strategy in the literature involves choosing Q as a diagonal matrix with relatively large values on the diagonal (e.g., $Q=10000I$) to emphasize state penalties while setting R as the identity matrix, $R=I$, to impose a standard control effort penalty. In the context of this study, the weighting matrices were specified as follows:

$\begin{array}{ccc}\hfill Q& =& \left[\begin{array}{cccccc}1000& 0& 0& 0& 0& 0\\ 0& 1000& 0& 0& 0& 0\\ 0& 0& 1000& 0& 0& 0\\ 0& 0& 0& 100& 0& 0\\ 0& 0& 0& 0& 1000& 0\\ 0& 0& 0& 0& 0& 100\end{array}\right],\hfill \\ \hfill R& =& \left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\hfill \end{array}$

Accordingly, the following optimal controller matrices are obtained:

$\begin{array}{ccc}\hfill K& =& \left[\begin{array}{cc}0& 0\\ 23.6613& -50.6547\\ 0& 0\\ -48.0394& 132.7923\end{array}\right]\hfill \end{array}$

$\begin{array}{ccc}\hfill A_f& =& \left[\begin{array}{cccc}0& 1.0000& 0& 0\\ -32.2124& -9.3863& 65.1860& 16.7827\\ 0& 0& 0& 1.0000\\ 61.5397& 14.7298& -172.6900& -45.9968\end{array}\right]\hfill \end{array}$

7.3. Case Examples

Case 1:

The setpoint is set to the state ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{4}}& 0& 0& 0\end{array}\right]}^T$ starting from the zero initial condition. The transient response of each robot joint is shown in Figure 7.

Case 2:

The setpoint is set to the state ${\left[\begin{array}{cccc}0& 0& -{\displaystyle \frac{\pi }{4}}& 0\end{array}\right]}^T$ starting from the zero initial condition. The transient response of each robot joint is illustrated in Figure 8.

Case 3:

The setpoint is initially set to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{2}}& 0& 0& 0\end{array}\right]}^T$ at 1 s, and then changed to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{2}}& 0& -{\displaystyle \frac{\pi }{4}}& 0\end{array}\right]}^T$ at 10 s, starting from zero initial conditions. The transient response of each robot joint is shown in Figure 9.

Disturbance effect:

The robot is assumed to be subjected to a load of mass p at a certain instant. Accordingly, Equation (31) can be rewritten as follows:

(35)$\begin{array}{cc}\hfill & \left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\ddot{\theta }+\left[\begin{array}{cc}-h{\dot{\theta }}_2& -h{\dot{\theta }}_1-h{\dot{\theta }}_2\\ h{\dot{\theta }}_1& 0\end{array}\right]\dot{\theta }+\left[\begin{array}{cc}\beta 1& 0\\ 0& \beta 2\end{array}\right]\dot{\theta }+\left[\begin{array}{c}{g}_1\\ g_2\end{array}\right]=\tau \hfill \end{array}$

where $\theta ={\left[{\theta }_1{\theta }_2\right]}^T$ denotes the joint angles, $\tau ={\left[{\tau }_1{\tau }_2\right]}^T$ represents the input torques at each joint, and

$\begin{array}{ccc}\hfill H_{11}& =& m_1{l}_{c1}^2+I_1+m_2(l_1^2+l_{c2}^2+2l_1{l}_{c2}cos{\theta }_2)+I_2+p(l_1^2+l_2^2+2l_1{l}_{2}cos{\theta }_2)\hfill \\ \hfill H_{12}& =& m_2{l}_1{l}_{c2}cos{\theta }_2+m_2{l}_{c2}^2+I_2+p{l}_2^2\hfill \\ \hfill H_{21}& =& H_{12}\hfill \\ \hfill H_{22}& =& m_2{l}_{c2}^2+I_2+p{l}_2^2\hfill \\ \hfill h& =& m_2{l}_1{l}_{c2}sin{\theta }_2+p{l}_1{l}_{2}sin{\theta }_2\hfill \\ \hfill g_1& =& m_1{l}_{c1}gsin{\theta }_1+m_{2}g(l_{c2}sin({\theta }_1+{\theta }_2)+l_{1}sin{\theta }_1)+pg(l_{2}sin({\theta }_1+{\theta }_2)+l_{1}sin{\theta }_1)\hfill \\ \hfill g_2& =& m_2{l}_{c2}gsin({\theta }_1+{\theta }_2)+p{l}_{2}gsin({\theta }_1+{\theta }_2)\hfill \end{array}$

Equation (31) can be expressed in the following compact form:

(36)$\begin{array}{c}\hfill H\left(\theta \right)\ddot{\theta }+C(\theta ,\dot{\theta })\dot{\theta }+B\dot{\theta }+g\left(\theta \right)=\tau \end{array}$

Case 4:

The setpoint is set to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{4}}& 0& 0& 0\end{array}\right]}^T$ at 1 s, and then changed to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{4}}& 0& -{\displaystyle \frac{\pi }{4}}& 0\end{array}\right]}^T$ at 10 s, starting from zero initial conditions. At 20 s, the robot is subjected to an external load of $p=0.5$ kg. The transient response of each robot joint is shown in Figure 10.

Case 5:

The setpoint is set to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{4}}& 0& 0& 0\end{array}\right]}^T$ at 1 s, and then changed to ${\left[\begin{array}{cccc}{\displaystyle \frac{3\pi }{2}}& 0& -{\displaystyle \frac{\pi }{4}}& 0\end{array}\right]}^T$ at 10 s, starting from zero initial conditions. At 20 s, the robot is subjected to a load of $p=0.5$ kg and measurement noise in the angle with standard deviation $\sigma =1^{\circ }$. Figure 11 depicts the transient response of the two-link robot.

Case 6:

To further analyze the robustness of the system, the setup from Case 5 is repeated with the measurement noise increased to $\sigma =2^{\circ }$. The transient response, shown in Figure 12, demonstrates that the system’s behavior remains largely unaffected despite the increased disturbances and noise.

Case 7:

Finally, Figure 13 shows the transient response of the two-link robot starting from non-zero initial conditions, specifically ${45}^{\circ }$ in the first joint and ${10}^{\circ }$ in the second.

The example cases clearly demonstrate that the proposed controller achieves robust and stable performance even in the presence of disturbances, load variations, and measurement noise. It produces a smooth, well-damped transient response while ensuring zero steady-state error.

8. Conclusions

This paper has presented a feedback linearization scheme integrated with Takagi–Sugeno (T-S) fuzzy systems to design a robust controller for nonlinear multivariable systems. The improved T-S modeling approach and the proposed controller construction enable the development of a linearized closed-loop system that is stable from any initial condition and demonstrates strong performance even in the presence of noise and load variations.

A two-link robot was used as a case study to illustrate the effectiveness, accuracy, and robustness of the closed-loop system. Simulation results confirm the smoothness and robustness of the controller, demonstrating that it maintains performance despite disturbances, load changes, and measurement noise. Overall, the controller exhibits well-damped behavior with zero steady-state error.

Author Contributions

Conceptualization, B.M.A.-H. and A.J.; methodology, B.M.A.-H.; software, A.J.; validation, B.M.A.-H. and A.J.; formal analysis, B.M.A.-H.; investigation, B.M.A.-H.; resources, A.J.; data curation, B.M.A.-H.; original draft preparation, A.J.; review and editing, B.M.A.-H. and J.B.R.; visualization, J.B.R.; supervision, B.M.A.-H.; project administration, B.M.A.-H.; funding acquisition, A.J. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 Fuzzy system membership functions.

Figure 2 Two-link robot.

Figure 3 Membership functions of ${\theta }_1$.

Figure 4 Membership functions of ${\dot{\theta }}_1$.

Figure 5 Membership functions of ${\theta }_2$.

Figure 6 Membership functions of ${\dot{\theta }}_2$.

Figure 7 The two-link robot transient response and the final value of ${\displaystyle \frac{3\pi }{4}}$ of the first link.

Figure 8 The two-link robot transient response and the final value of ${\displaystyle \frac{-\pi }{4}}$ of the second link.

Figure 9 The two-link robot transient response and the final values of ${\displaystyle \frac{3\pi }{4}}$ and ${\displaystyle \frac{-\pi }{4}}$ of the first and second links, respectively.

Figure 10 The two-link robot transient response and the load effect.

Figure 11 The two-link robot transient response and the load effect and measurement noise of $\sigma =1^{\circ }$.

Figure 12 The two-link robot transient response and the load effect and measurement noise of $\sigma =2^{\circ }$.

Figure 13 The two-link robot transient response with initial conditions (${45}^{\circ },{10}^{\circ }$).

References

1. Floares, A. Feedback linearization using neural networks applied to advanced pharmacodynamic and pharmacogenomic systems. Proceedings of the International Joint Conference on Neural Networks 2005; Montreal, QC, Canada, 31 July–4 August 2005; Volume 1, pp. 173-178.

2. Herníndez-Santiago, J.; Garrido-Meléndez, J.; Escobedo-Trujillo, B.A.; Alcalá-Peréa, G. Exact Feedback Linearization Technique Applied to an Inverted Pendulum Coupled to a DC Motor. Proceedings of the 2023 IEEE International Conference on Engineering Veracruz (ICEV), Boca del Río; Veracruz, Mexico, 23–26 October 2023; IEEE: Piscataway, NJ, USA, 2023.

3. Swain, N.; Malik, S.; Pati, N. Design and Analysis of Step up Regulator using Exact Feedback Linearization by State Feedback Approach. Proceedings of the 2021 19th OITS International Conference on Information Technology (OCIT); Bhubaneswar, India, 16–18 December 2021; IEEE: Piscataway, NJ, USA, 2021.

4. Fang, Z.; Yu, Z.; Huang, Q.; Wang, Y.; Gu, X. Research on design and control method of active vibration isolation system based on piezoelectric Stewart platform. Sci. Rep.; 2025; 15, 944. [DOI: https://dx.doi.org/10.1038/s41598-024-84980-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39762358]

5. Mihaly, V.; Şuşcă, M.; Birlescu, I.; Sim, S.; Pisla, D.; Dobra, P. Robust Feedback Linearization for Serial Robots. Proceedings of the 2024 IEEE International Conference on Automation, Quality and Testing, Robotics (AQTR); Cluj-Napoca, Romania, 16–18 May 2024; IEEE: Piscataway, NJ, USA, 2024.

6. Gao, Q.; Liu, L.; Feng, G.; Wang, Y.; Qiu, J. Universal Fuzzy Integral Sliding-Mode Controllers Based on T–S Fuzzy Models. IEEE Trans. Fuzzy Syst.; 2013; 22, pp. 350-362. [DOI: https://dx.doi.org/10.1109/TFUZZ.2013.2254717]

7. Sarbaz, M.; Manthouri, M.; Zamani, I. Rough Neural Network and Adaptive Feedback Linearization Control Based on Lyapunov Function. Proceedings of the 2021 7th International Conference on Control, Instrumentation and Automation (ICCIA); Tabriz, Iran, 23–24 February 2021; IEEE: Piscataway, NJ, USA, 2021.

8. Park, C.-W.; Kang, H.-J.; Yee, Y.-H.; Park, M. Numerical robust stability analysis of fuzzy feedback linearisation regulator based on linear matrix inequality approach. IEE Proc. Control. Theory Appl.; 2002; 149, pp. 82-88. [DOI: https://dx.doi.org/10.1049/ip-cta:20020255]

9. Takagi, T.; Sugeno, M. Fuzzy Identification of Systems and Its Applications to Modeling and Control. IEEE Trans. Syst. Man Cybern.; 1985; 15, pp. 116-132. [DOI: https://dx.doi.org/10.1109/TSMC.1985.6313399]

10. Zeng, K.; Zhang, N.Y.; Xu, W.L. A comparative study on sufficient conditions for Takagi-Sugeno fuzzy systems as universal approximators. IEEE Trans. Fuzzy Syst.; 2000; 8, pp. 773-780.

11. Al-Hadithi, B.M.; Jiménez, A.; Matía, F. Variable Structure Control with Chattering Reduction of a Generalized T-S Model. Asian J. Control.; 2012; 15, pp. 155-168. [DOI: https://dx.doi.org/10.1002/asjc.518]

12. Gang, F. A survey on analysis and design of model-based fuzzy control systems. IEEE Trans. Fuzzy Syst.; 2006; 14, pp. 676-697. [DOI: https://dx.doi.org/10.1109/tfuzz.2006.883415]

13. El-Sousy, F.F.M.; Elmorshedy, M.F.; Amin, M.M.; Mohammed, O.A. Robust Adaptive Feedback Linearization Control Using Online Neural-Network Estimators for Uncertain Linear Induction Motor Drive System. Proceedings of the 2021 13th International Symposium on Linear Drives for Industry Applications (LDIA); Wuhan, China, 1–3 July 2021; IEEE: Piscataway, NJ, USA, 2021.

14. Huang, C.-J.; Li, T.-H.S.; Chen, C.-C. Fuzzy Feedback Linearization Control for MIMO Nonlinear System and Its Application to Full-Vehicle Suspension System. Circuits Syst. Signal Process.; 2009; 28, pp. 959-991. [DOI: https://dx.doi.org/10.1007/s00034-009-9126-3]

15. Chung-Cheng, C.; Chao-Hsing, H.; Ying-Jen, C.; Yen-Feng, L. Disturbance attenuation of nonlinear control systems using an observer-based fuzzy feedback linearization control. Chaos Solitons Fractals; 2007; 33, pp. 885-900.

16. THSLi, CJHuang, CCChen, Novel fuzzy feedback linearization strategy for control via differential geometry approach. ISA Trans.; 2010; 49, pp. 348-357. [DOI: https://dx.doi.org/10.1016/j.isatra.2010.02.004]

17. Xiaoshen, L.; Xuehai, Y.; Mingzuo, J.; Chunling, Z. Fuzzy inference modeling method of time-varying system based on TS fuzzy system. Proceedings of the 2015 12th International Conference on Fuzzy Systems and Knowledge Discovery (FSKD); Zhangjiajie, China, 15–17 August 2015; IEEE: Piscataway, NJ, USA, 2015.

18. Salgado, C.M.; Viegas, J.L.; Azevedo, C.S.; Ferreira, M.C.; Vieira, S.M.; Sousa, J.M.C. Takagi–Sugeno Fuzzy Modeling Using Mixed Fuzzy Clustering. IEEE Trans. Fuzzy Syst.; 2016; 25, pp. 1417-1429. [DOI: https://dx.doi.org/10.1109/TFUZZ.2016.2639565]

19. Zeng, X.J. A comparison between TS fuzzy systems and affine TS fuzzy systems as nonlinear control system models. Proceedings of the 2014 IEEE International Conference on Fuzzy Systems (FUZZ-IEEE); Beijing, China, 6–11 July 2014; IEEE: Piscataway, NJ, USA, 2014.

20. Naizheng, S.; Junmin, L.; Lunle, Y. Adaptive fuzzy control of nonlinear systems based on multiple inputs TS fuzzy bilinear model. Proceedings of the 31st Chinese Control Conference; Hefei, China, 25–27 July 2012; IEEE: Piscataway, NJ, USA, 2012.

21. Chang-Woo, P. Robust Stable Fuzzy Control via Fuzzy Modeling and Feedback Linearization with Its Applications to Controlling Uncertain Single-link Flexible Joint Manipulators. J. Intell. Robot. Syst.; 2004; 39, 131147. [DOI: https://dx.doi.org/10.1023/b:jint.0000015344.84152.dd]

22. Kim, S.W.; Lee, J.J. Fuzzy control of a class of multivariable nonlinear systems subject to parameter uncertainties: Model reference approach. Int. J. Approx. Reason.; 2001; 26, pp. 129-144. [DOI: https://dx.doi.org/10.1016/s0888-613x(00)00068-2]

23. Purwar, S.; Kar, I.N.; Jha, A.N. Adaptive control of robot manipulators using fuzzy logic systems under actuator constraints. Fuzzy Sets Syst.; 2005; 152, pp. 651-664. [DOI: https://dx.doi.org/10.1016/j.fss.2004.11.012]

24. Labiod, S.; Boucherit, M.S.; Guerra, T.M. Adaptive fuzzy control of a class of MIMO nonlinear systems. Fuzzy Sets Syst.; 2005; 15, pp. 59-77. [DOI: https://dx.doi.org/10.1016/j.fss.2004.10.009]

25. Song, Z.; Yi, J.; Zhao, D.; Li, X. A computed torque controller for uncertain robotic manipulator systems: Fuzzy approach. Fuzzy Sets Syst.; 2006; 154, pp. 208-226. [DOI: https://dx.doi.org/10.1016/j.fss.2005.03.007]

26. Tsai, C.H.; Wang, C.H.; Lin, W.S. Robust fuzzy model following control of robot manipulators. IEEE Trans. Fuzzy Syst.; 2000; 8, pp. 462-469. [DOI: https://dx.doi.org/10.1109/91.868952]

27. Xu, Y.; Bao, R.; Zhang, L.; Wang, J.; Wang, S. Embodied intelligence in RO/RO logistic terminal: Autonomous intelligenttransportation robot architecture. Sci. China Inf. Sci.; 2025; 68, pp. 1-17. [DOI: https://dx.doi.org/10.1109/91.868952]

28. Chen, Z.; Zhan, G.; Jiang, Z.; Zhang, W.; Rao, Z.; Wang, H.; Li, J. Adaptive impedance control for docking robot via Stewart parallel mechanism. ISA Trans.; 2024; 155, pp. 361-372. [DOI: https://dx.doi.org/10.1109/91.868952]

29. Al-Hadithi, B.M.; Jiménez, A.; Matía, F. A new approach to fuzzy estimation of Takagi-Sugeno model and its applications to optimal control for nonlinear systems. Appl. Soft Comput.; 2012; 12, pp. 280-290. [DOI: https://dx.doi.org/10.1016/j.asoc.2011.08.044]

30. Jiménez, A.; Al-Hadithi, B.M.; Matía, F. An Optimal T-S Model for the Estimation and Identification of Nonlinear Functions. WSEAS Trans. Syst. Control.; 2008; 3, pp. 897-906.

31. Jiménez, A.; Al-Hadithi, B.M.; Matía, F. Improvement of Takagi-Sugeno Fuzzy Model for the Estimation of Nonlinear Functions. Asian J. Control.; 2012; 14, pp. 320-334. [DOI: https://dx.doi.org/10.1002/asjc.310]

32. Louda, S.; Karkar, N.; Seghir, F.; Boutalbi, O. Fuzzy Dynamic Feedback Linearization for Efficient Mobile Robot Trajectory Tracking and Obstacle Avoidance in Autonomous Navigation. Int. J. Robot. Control. Syst.; 2025; 5, pp. 881-901. [DOI: https://dx.doi.org/10.1002/asjc.310]

33. Yesil, B.; Sahin, S. Real-Time Implementation of a Microcontroller-Based Coupled-TankWater Level Control System with Feedback Linearization and Fuzzy Logic Controller Algorithms. Sensors; 2025; 25, 1279. [DOI: https://dx.doi.org/10.3390/s25051279]

34. Al-Ani, F.R.; Lutfy, O.F.; Al-Khazraji, H. Optimal Synergetic and Feedback Linearization Controllers Design for Magnetic Levitation Systems: A Comparative Study. J. Robot Control (JRC); 2024; 6, pp. 22-30. [DOI: https://dx.doi.org/10.18196/jrc.v6i1.24452]

35. Al-Hadithi, B.M.; Jiménez, A.; Pérez-Oria, J.; Alonso, L. Optimal Control Using Feedback Linearization for a Generalized T-S Model. IFIP International Federation for Information Processing 2014, AIAI 2014, Hangzhou, China, 17–20 October 2014; IFIP AICT 436 Springer: Berlin/Heidelberg, Germany, 2014; pp. 466-475.