<figure><img src="image_1.png" /></figure>

HAL
open science

# Joint state and parameter estimation for discrete-time Takagi-Sugeno model

Krishnan Srinivasarengan, José Ragot, Didier Maquin, Christophe Aubrun

<figure><img src="image_2.png" /></figure>

## To cite this version:

Krishnan Srinivasarengan, José Ragot, Didier Maquin, Christophe Aubrun. Joint state and parameter estimation for discrete-time Takagi-Sugeno model. 13th European Workshop on Advanced Control and Diagnosis, ACD 2016, Nov 2016, Lille, France. hal-01399791

HAL Id: hal-01399791

https://hal.science/hal-01399791v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Joint State and Parameter Estimation for Discrete-Time Takagi-Sugeno Model

Krishnan Srinivasarengan, José Ragot, Didier Maquin, Christophe Aubrun

CRAN UMR 7039, Université de Lorraine, Vandoeuvre-lès-Nancy, Cedex, France
CRAN UMR 7039, CNRS, France

E-mail: krishnan.srinivasarengan@univ-lorraine.fr, jose.ragot@univ-lorraine.fr, didier.maquin@univ-lorraine.fr, christophe.aubrun@univ-lorraine.fr

**Abstract.** Takagi-Sugeno (T-S) models have been popular in analyzing nonlinear systems. Observer designs for T-S models have focused on both cases with measured and unmeasured premise variables. However, the unmeasured premise variables have to be one of the states of the system. If one of the inputs is a premise variable, these approaches are not applicable. A recent literature proposed a joint state and parameter estimation for time-varying unknown parameters in systems of a special type [1]. A discrete-time version of the results is developed in this paper. The modeling approach is described and the stability analysis of the estimation is given using the Lyapunov method with the $L_2$ approach handling the uncertainties.

## 1. Introduction

Takagi-Sugeno (T-S) type of systems is a class of quasi-LPV (Linear Parameter Varying) systems whose model structure could depend on inputs, outputs, and states of the system. Observer design for T-S models has been popular for a long time. The sector nonlinearity (SNL) approach is typically used to obtain the T-S equivalent model of a nonlinear system that are an exact representation within a sector. However, the use of this approach gives rise to unmeasured states as premise variables. Since the work of [2], there has been significant interest in developing observers for systems with unmeasured premise variable. These works have been extended to cascaded [3], distributed [4] and other such scenarios. However, there are limited works when it comes to dealing with the case when premise variables contain an unknown input. In this circumstance, the unknown input observer formulations such as [5] wouldn't be applicable.

In [6], an interesting approach was proposed to represent time varying unknown parameters using the sector nonlinearity approach. The estimation of the unknown parameters followed, in spirit, the approaches proposed for observer design with unmeasured premise variables. In [1], this was combined with a T-S system with unmeasured premise variables to develop a joint state and parameter estimation observer. The unknown parameters are represented using the sector nonlinearity approach, assuming that their bounds are known. This leads to a set of weighting functions that depend on estimated parameters and forms a product of estimated weighting functions. Various approaches assume the knowledge of Lipschitz constant that bound the difference between the estimated and actual weighting function. However, the convex sum property

<figure><img src="image_4.png" /></figure>
along with matrix theory results are used to bound the uncertain components. Further, an L₂ approach is used to guarantee a bounded error on the estimations.

In this work, the discrete-time version of this result is solved. On the outset, the derivation follows similar steps, however, the discrete-time version of some preliminary results (especially the Bounded Real Lemma) makes the task of adapting the existing results to discrete-time a significant challenge. The key contributions are in formulating the problem so as to avoid non-linear matrix inequalities, resolving nonlinear matrix inequalities by simplifying them to linear matrix inequalities (LMIs), developing bounds on the time varying uncertainties and providing the necessary LMI conditions for the convergence of the estimation. The organization is as follows: The Sec.2 provides the preliminaries necessary for the results developed in this work including the definition of model structures. The main results are presented in Sec.3 and an example to illustrate the results are given in 4. The Sec.5 provides some concluding remarks with comments on the way forward.

## 2. Preliminaries and Problem Formulation

Takagi-Sugeno models are of the form,

$$ 
\begin{aligned}
x(k + 1) &= \sum_{i=1}^{r} \mu_i(z(k)) [A_i x(k) + B_i u(k)] \\
y(k) &= C x(k)
\end{aligned}
\quad (1) $$

In this work, the output equation is considered to be linear. However, the results can be readily extended to cases with the output equation as,

$$ y(k) = \sum_{i=1}^{r} \mu_i(z(k)) [C_i x(k) + D_i u(k)] \quad (2) $$

Here, $r = 2^{n_p}$, where $n_p$ is the number of premise variables. Further,

$$ x \in \mathbb{R}^n, \quad u \in \mathbb{R}^{n_u}, \quad z \in \mathbb{R}^{n_p}, \quad y \in \mathbb{R}^{n_y} $$

and

$$ A_i \in \mathbb{R}^{n_x \times n_x}, \quad B_i \in \mathbb{R}^{n_x \times n_u}, \quad C \in \mathbb{R}^{n_x \times n_y}, \forall i \quad (3) $$

$\mu_i(z(k))$ corresponds to the weighting function obtained from the product of corresponding membership functions of the premise variables (see for e.g., [7] for more details). The weighting functions satisfy the convex sum property, such that,

$$ \sum_{i} \mu_{i}(.) = 1 \quad \text{and} \quad 0 \le \mu_{i}(.) \le 1, \forall i \quad (4) $$

### 2.1. Preliminary Results

The following results would be used during the proof of the results in this work,

**Lemma 2.1** [8] The following LMIs, for appropriately suitable dimensions, are equivalent:

$$ A^T P A - Q < 0, \quad P > 0 \quad (5) $$

$$ \begin{bmatrix} -Q & A^T P \\ PA & -P \end{bmatrix} < 0, \quad P > 0 \quad (6) $$
**Lemma 2.2** [9] Consider two matrices X and Y with appropriate dimensions, a time varying matrix $\Delta(k)$ and a positive scalar $\lambda$. The following property is verified:

$$X^T \Delta^T(k) Y + Y^T \Delta(k) X \leq \lambda X^T X + \lambda^{-1} Y^T Y \quad (7)$$

for $\Delta^T(k)\Delta(k) \leq I$

**Lemma 2.3** (Discrete-time Bounded Real Lemma [10]) For a discrete-time system, of the form,

$$\begin{aligned} x(k + 1) &= Ax(k) + Bu(k) \\ y(k) &= Cx(k) + Du(k) \end{aligned} \quad (8)$$

the Bounded Real Lemma equivalent LMI condition for stability with an $\mathcal{L}_2$ gain $\Gamma$ is,

$$P > 0, \quad \begin{bmatrix} A^T P A - P + C^T C & A^T P B + C^T D \\ B^T P A + D^T C & D^T D + B^T P B - \Gamma \end{bmatrix} \leq 0 \quad (9)$$

**Lemma 2.4** (Schur's Complement) For a symmetric matrix M, given by,

$$M = \begin{bmatrix} A & B \\ B^T & C \end{bmatrix}$$

if C is invertible, then the following properties hold:

(i) $M \succ 0$ iff $C \succ 0$ and $A - BC^{-1}B^T \succ 0$

(ii) if $C \succ 0$, then $M \succeq 0$ iff $A - BC^{-1}B^T \succeq 0$.

## 2.2. System Model Structure

Consider discrete-time Takagi Sugeno type systems of the form,

$$\begin{aligned} x(k + 1) &= \sum_{i=1}^{r} \mu_i(x) [A_i(\theta(k))x(k) + B_i(\theta(k))u(k)] \\ y(k) &= Cx(k) \end{aligned} \quad (10)$$

where, $\theta \in \mathbb{R}^{n_{\theta}}$. In this paper, it is considered that $\theta$ is in a bilinear form with the premise variables, which could be the states, outputs, and the inputs of the system, that is:

$$\begin{aligned} A_i(\theta(k)) &= A_i + \sum_{j=1}^{n_{\theta}} \bar{A}_{ij} \theta_j(k) \\ B_i(\theta(k)) &= B_i + \sum_{j=1}^{n_{\theta}} \bar{B}_{ij} \theta_j(k) \end{aligned} \quad (11)$$

These matrices could be rewritten by applying sector nonlinearity on the unknown parameter, assuming that the extremum values are known. With this assumption, the time varying matrix is then rewritten as,

$$\begin{aligned} A_i(\theta(k)) &= A_i + \sum_{j=1}^{n_{\theta}} \sum_{m=1}^{2} \tilde{\mu}_j^m(\theta_j(k)) \bar{A}_{ij} \theta_j^m \\ B_i(\theta(k)) &= B_i + \sum_{j=1}^{n_{\theta}} \sum_{m=1}^{2} \tilde{\mu}_j^m(\theta_j(k)) \bar{B}_{ij} \theta_j^m \end{aligned} \quad (12)$$
where, $\tilde{\mu}_j^m(\theta_j(k))$ refers to the weighting function corresponding to mth submodel of $\theta_j$. Referring to further steps in [1], the overall system could be rewritten as,

$$ 
\begin{aligned}
x(k + 1) &= \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(x) \tilde{\mu}_j(\theta) [A_{ij}x(k) + B_{ij}u(k)] \\
y(k) &= Cx(k)
\end{aligned} 
\qquad (13) 
$$

where $s \triangleq 2^{n_{\theta}}$,

$$ \forall i, j, \quad A_{ij} = A_i + \sum_{j=1}^{n_{\theta}} \theta_j^m \bar{A}_j \qquad B_{ij} = B_i + \sum_{j=1}^{n_{\theta}} \theta_j^m \bar{B}_j \qquad (14) $$

where, $\theta_j^k$ is the minimum or maximum value of $\theta_j$ depending upon the submodel under consideration ($m = 1, 2$).

### 2.3. Observer Model Structure

The structure of the observer envisaged to design has a Luenberger like structure for the state estimation part and a first order form for the parameter estimation part, it is given by,

$$ \hat{x}(k+1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [A_{ij}\hat{x}(k) + B_{ij}u(k) + L_{ij}(y(k) - \hat{y}(k))] \qquad (15) $$

$$ \hat{\theta}(k+1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [K_{ij}(y(k) - \hat{y}(k)) - \eta\theta(k)] \qquad (16) $$

The gains $L_{ij} \in \mathbb{R}^{n_x \times n_y}$ and $K_{ij} \in \mathbb{R}^{n_{\theta} \times n_y}$ are to be estimated with the gain $\eta \in \mathbb{R}^{n_{\theta} \times n_{\theta}}$ is chosen. The choice of $\eta$ shall typically be in the form of a diagonal matrix. The reason for this choice is to avoid the presence of nonlinear matrix inequality terms that arise. In the initial work [1], this was introduced to avoid a marginal stability condition for the error dynamics. As discussed in [11], choosing this reduces the number of variables in the final LMI to solve and hence allows for a computational tractable problem.

## 3. Results

The main result in this work is stated as the theorem below.

**Theorem 3.1** Given the system model of the form (13), there exists an observer of the form (15)-(16), if there exists $P_0, P_1, R_{ij}, F_{ij}, \lambda_1, \lambda_2, \lambda_3, \lambda_4, \Gamma_2^k$ ($\forall i, \forall j, \forall k \in (0, 1, 2, 3)$), such that,

$$ P_0 > 0 \qquad (17) $$

$$ P_1 > 0 \qquad (18) $$

$$ \left( 
\begin{array}{cccccccc}
-P+Q_e & Q_A & \Phi_{ij}^T P & 0 & P_A A & P_A B & 0 & 0 \\
Q_A^T & T_{22} & 0 & Q_B^T & 0 & 0 & 0 & 0 \\
P \Phi_{ij} & 0 & -P & 0 & 0 & 0 & 0 & 0 \\
0 & Q_B^T & 0 & -P & 0 & 0 & P_0 A & P_0 B \\
A^T P_A & 0 & 0 & 0 & -\lambda_1 I & 0 & 0 & 0 \\
B^T P_A & 0 & 0 & 0 & 0 & -\lambda_3 I & 0 & 0 \\
0 & 0 & 0 & A^T P_0 & 0 & 0 & -\lambda_2 I & 0 \\
0 & 0 & 0 & B^T P_0 & 0 & 0 & 0 & -\lambda_4 I
\end{array}
\right) < 0 \qquad (19) $$
Here, $Q_e \triangleq \text{diag}(Q_{e0}, Q_{e1})$ in line with $P \triangleq \text{diag}(P_0, P_1)$ and

$$P_A \triangleq A_{ij}^T P_0 - C^T R_{ij}^T \qquad (20)$$

$$T_{22} = \begin{bmatrix} -\Gamma_2^x + (\lambda_1 + \lambda_2)E_A^T E_A & 0 & 0 & 0 \\ 0 & -\Gamma_2^x + (\lambda_3 + \lambda_4)E_B^T E_B & 0 & 0 \\ 0 & 0 & -\Gamma_2^\theta & 0 \\ 0 & 0 & 0 & -\Gamma_2^{\Delta\theta} \end{bmatrix} \qquad (21)$$

$$Q_A = \begin{bmatrix} 0 & 0 & -C^T F_{ij}^T (I + \eta) & -C^T F_{ij}^T \\ 0 & 0 & -\eta^T P_1 (I + \eta) & -\eta^T P_1 \end{bmatrix} \qquad Q_B = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & (I + \eta)^T P_1 \\ 0 & P_1 \end{bmatrix}$$

$$\Phi_{ij}^T P = \begin{bmatrix} P_A & -C^T F_{ij}^T \\ 0 & -\eta^T P_1 \end{bmatrix}$$

and,

$$A = [A_{11} \quad \dots \quad A_{rs}], \quad B = [B_{11} \quad \dots \quad B_{rs}],$$

$$E_A = [I_{n_x} \quad \dots \quad I_{n_x}], \quad E_B = [I_{n_u} \quad \dots \quad I_{n_u}]$$

$$\Sigma_A(k) = \text{diag}(\delta_{11}(k)I_{n_x}, \dots, \delta_{rs}(k)I_{n_x}), \quad \Sigma_B(k) = \text{diag}(\delta_{11}(k)I_{n_u}, \dots, \delta_{rs}(k)I_{n_u}),$$

$$\delta_{ij}(k) = \mu_i(x)\tilde{\mu}_j(\theta) - \mu_i(\hat{x})\tilde{\mu}_j(\hat{\theta}) \qquad (22)$$

The observer gains are given by,

$$K_{ij} = P_1^{-1} F_{ij}, \quad L_{ij} = P_0^{-1} R_{ij} \qquad (23)$$

**Proof** Given that the system model and observer are not in same form, comparing the error dynamics would be difficult. Hence the original system (13) is represented in an uncertain-like form,

$$x(k + 1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [(A_{ij} + \Delta A(k))x(k) + (B_{ij} + \Delta B(k))u(k)] \\ y(k) = Cx(k) \qquad (24)$$

with,

$$\Delta A(k) = A \Sigma_B(k) E_A, \quad \Delta B(k) = B \Sigma_B(k) E_B \qquad (25)$$

where, the components are as given in (22). Due to the convex sum property of the weighting functions, we also have,

$$-1 \le \delta_{ij}(k) \le 1, \quad \Sigma_A^T(k)\Sigma_A(k) \le I \quad \Sigma_B^T(k)\Sigma_B(k) \le I \qquad (26)$$

Defining $e_\theta(k) \triangleq \theta(k) - \hat{\theta}(k)$ and $\Delta\theta(k) \triangleq \theta(k + 1) - \theta(k)$, the estimation errors are given as below:

$$e_x(k + 1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [(A_{ij} - L_{ij}C)e_x(k) + \Delta A(k)x(k) + \Delta B(k)u(k)] \qquad (27)$$

$$e_{\theta}(k + 1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [\theta(k + 1) - K_{ij} C e_x(k) + \eta \hat{\theta}] \\ = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [\Delta\theta(k) + (I + \eta)\theta(k) - K_{ij} C e_x(k) - \eta e_\theta(k)] \qquad (28)$$
The augmented error equation is given by,

$$e_a(k + 1) = \sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) [\Phi_{ij} e_a(k) + \Psi_{ij}(k) \tilde{u}] \quad (29)$$

where,

$$\Phi_{ij} = \begin{bmatrix} A_{ij} - L_{ij}C & 0 \\ -K_{ij}C & -\eta \end{bmatrix} \quad \text{and} \quad \Psi_{ij}(k) = \begin{bmatrix} \Delta A(k) & \Delta B(k) & 0 & 0 \\ 0 & 0 & I + \eta & I \end{bmatrix} \quad (30)$$

with $e_a \triangleq [e_x \ e_\theta]^T$ and $\tilde{u}(k) \triangleq [x(k) \ u(k) \ \theta(k) \ \Delta\theta(k)]^T$.

Consider a quadratic Lyapunov candidate $V(k) = e_a^T(k)Pe_a(k)$ and an extended stability condition such that,

$$V(k + 1) - V(k) + e_a^T(k)Q_e e_a(k) - \tilde{u}^T(k)\Gamma_2\tilde{u}(k) < 0 \quad (31)$$

where $Q_e$ is a scaling factor for the errors that is chosen and $\Gamma_2 = \text{block}(\Gamma_2^x, \Gamma_2^u, \Gamma_2^\theta, \Gamma_2^{\Delta\theta})$ is the $L_2$ gain value under which the uncertainties $(\tilde{u}(t))$ are expected to be bounded. By applying the discrete-time Bounded Real Lemma (BRL) (Lemma 2.3) we get,

$$\begin{bmatrix} \Phi_{ij}^T P \Phi_{ij} - P + Q_e & \Phi_{ij}^T P \Psi_{ij}(k) \\ (\Phi_{ij}^T P \Psi_{ij}(k))^T & \Psi_{ij}^T(k) P \Psi_{ij}(k) - \Gamma_2 \end{bmatrix} < 0 \quad (32)$$

There are nonlinear terms in this matrix condition. Specifically, the term $\Phi_{ij}P\Phi_{ij}$ is not quadratic as it appears, but it includes cross terms where the index for the first $\Phi$ is different from that of the second $\Phi$. However, based on the illustration in Theorem 17 in [12] on non-symmetric terms in the discrete-time models, it could be shown that ignoring the cross terms would not affect the negativity of the inequality condition.

Given the structure of $\Phi$ and $\Psi(k)$ in (30), there are likely nonlinear terms arising due to the quadratic structure in (32). At first however, applying the Lemma 2.1 on the (1, 1) and (2, 2) blocks of (32) leads to,

$$\begin{align} \Phi_{ij}^T P \Phi_{ij} - P + Q_e &\equiv \begin{bmatrix} -P + Q_e & \Phi_{ij}^T P \\ P \Phi_{ij} & -P \end{bmatrix} < 0 \nonumber \\ \Psi_{ij}^T(k) P \Psi_{ij} - \Gamma_2 &\equiv \begin{bmatrix} -\Gamma_2 & \Psi_{ij}^T(k) P \\ P \Psi_{ij}(k) & -P \end{bmatrix} < 0 \tag{33} \end{align}$$

Inserting these into (32), we get,

$$\begin{bmatrix} -P + Q_e & \Phi_{ij}^T P \Psi_{ij}(k) & \Phi_{ij}^T P & 0 \\ (\Phi_{ij}^T P \Psi_{ij}(k))^T & -\Gamma_2 & 0 & \Psi_{ij}^T(k) P \\ P \Phi_{ij} & 0 & -P & 0 \\ 0 & P \Psi_{ij}(k) & 0 & -P \end{bmatrix} < 0 \quad (34)$$

In what follows, the aim is to reduce the matrix inequality in (34), which is both time varying and possibly nonlinear, to an LMI. The key steps are as follows:

* Handling the nonlinearity in the matrix entries

- Enforce a structure on the positive definite matrix $P$ such that, $P = \text{diag}(P_0, P_1)$.
– Rewrite quadratic terms in variables into new variables.

• Handling the time varying terms in the matrix entries

– Partition the time varying terms so that the Lemma 2.2 could be applied

The diagonal structure enforcing of the Lyapunov function is followed by the following change of variables:

$$R_{ij} = P_0 L_{ij}, \quad F_{ij} = P_1 K_{ij} \qquad (35)$$

This leads to the terms in (1,2), (1,3), and (2,4) blocks in (34) to become,

$$\Phi_{ij}^T P \Psi_{ij}(k) = \begin{bmatrix} (A_{ij}^T P_0 - C^T R_{ij}^T) \Delta A(k) & (A_{ij}^T P_0 - C^T R_{ij}^T) \Delta B(k) & -(C^T F_{ij}^T)(I+\eta) & -C^T F_{ij}^T \\ 0 & 0 & -\eta^T P_1(I+\eta) & -\eta^T P_1 \end{bmatrix} \quad (36)$$

and,

$$\Phi_{ij}^T P = \begin{bmatrix} A_{ij}^T P_0 - C^T R_{ij}^T & -C^T F_{ij}^T \\ 0 & -\eta^T P_1 \end{bmatrix}$$

$$\Psi_{ij}^T(k) P = \begin{bmatrix} \Delta A^T(k) P_0 & 0 \\ \Delta B^T(k) P_0 & 0 \\ 0 & (I + \eta)^T P_1 \\ 0 & P_1 \end{bmatrix} \quad (37)$$

The terms $\Psi_{ij}(k)$ and $\Phi_{ij}^T P \Psi_{ij}(k)$ in (34) have time varying terms which need to be bounded. To do it, the original error dynamics equation in (29) is split into terms with and without time varying components, that is,

$$\sum_{i=1}^{r} \sum_{j=1}^{s} \mu_i(\hat{x}) \tilde{\mu}_j(\hat{\theta}) (Q_{ij} + \mathcal{L}_{ij}(k) + \mathcal{L}_{ij}^T(k)) < 0 \quad (38)$$

where,

$$Q_{ij} = \begin{bmatrix} -P + Q_e & Q_A & \Phi_{ij}^T P & 0 \\ Q_A^T & -\Gamma_2 & 0 & Q_B \\ P \Phi_{ij} & 0 & -P & 0 \\ 0 & Q_B^T & 0 & -P \end{bmatrix} \quad (39)$$

with $Q_A$ and $Q_B$ given by,

$$Q_A = \begin{bmatrix} 0 & 0 & -C^T F_{ij}^T (I + \eta) & -C^T F_{ij}^T \\ 0 & 0 & -\eta^T P_1 (I + \eta) & -\eta^T P_1 \end{bmatrix}, \quad Q_B = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & (I + \eta)^T P_1 & 0 \\ 0 & 0 & P_1 & 0 \end{bmatrix} \quad (40)$$

And the time varying terms are gathered as,

$$\begin{align}
\mathcal{L}_{ij}(k) &= \begin{bmatrix} 0 & \mathcal{L}_A(k) & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \mathcal{L}_B(k) & 0 & 0 \end{bmatrix} \nonumber \\
\mathcal{L}_A(k) &= \begin{bmatrix} (A_{ij}^T P_0 - C^T R_{ij}^T) \Delta A(k) & (A_{ij}^T P_0 - C^T R_{ij}^T) \Delta B(k) & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \nonumber \\
\mathcal{L}_B(k) &= \begin{bmatrix} P_0 \Delta A(k) & P_0 \Delta B(k) & 0 & 0 \end{bmatrix} \tag{41}
\end{align}$$
and their corresponding transposes. To bound the time varying matrices, the Lemma 2.2 is used. Each of the four time varying components (two in $\mathcal{L}_{ij}(k)$ and two in its transpose) would be bounded by applying the Lemma individually.

We define, for simplified notation purposes, $P_A \triangleq A_{ij}^T P_0 - C^T R_{ij}^T$. We can write each of the summation element in an expanded form by splitting $\mathcal{L}_{ij}(k)$ as follows:

$$ \mathcal{L}_{ij}(k) = \begin{bmatrix} 0 & [P_A \Delta A(k) & 0] & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & P_0 \Delta A(k) & 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & P_0 \Delta B(k) & 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} \quad (42) \\ = \mathcal{L}_{A1}(k) + \mathcal{L}_{A2}(k) + \mathcal{L}_{B1}(k) + \mathcal{L}_{B2}(k) \quad (43) $$

The zeros are not of the same size, but have appropriate dimensions. The $\Delta A(k)$ and $\Delta B(k)$ contain the time varying terms $\Sigma_A(k)$ and $\Sigma_B(k)$ as given in (26). To do this, the matrices need to be split making use of the components as in (25), so that we get,

$$ \mathcal{L}_{A1} = \begin{pmatrix} P_A & \mathcal{A} \\ 0 & \end{pmatrix} \Sigma_A(k) \begin{pmatrix} 0 & E_A & 0 \end{pmatrix} \qquad \mathcal{L}_{B1} = \begin{pmatrix} P_A & \mathcal{B} \\ 0 & \end{pmatrix} \Sigma_B(k) \begin{pmatrix} 0 & E_B & 0 \end{pmatrix} \quad (44) $$

$$ \mathcal{L}_{A2} = \begin{pmatrix} 0 \\ P_0 \mathcal{A} \\ 0 \end{pmatrix} \Sigma_A(k) \begin{pmatrix} 0 & E_A & 0 \end{pmatrix} \qquad \mathcal{L}_{B2} = \begin{pmatrix} 0 \\ \mathcal{B}^T P_0 \\ 0 \end{pmatrix} \Sigma_B(k) \begin{pmatrix} 0 & E_B & 0 \end{pmatrix} \quad (45) $$

And this could be further extended to their corresponding transposes. For each of the matrix and transpose pair, we can apply the Lemma 2.2, so that,

$$ \mathcal{L}_{A1}(k) + \mathcal{L}_{A1}^T(k) \le \lambda_1^{-1} \begin{pmatrix} P_A & \mathcal{A} \\ 0 & \end{pmatrix} (\mathcal{A}^T P_A \quad 0) + \lambda_1 \begin{pmatrix} 0 \\ E_A^T \\ 0 \end{pmatrix} \begin{pmatrix} 0 & E_A & 0 \end{pmatrix} \\ \le \begin{bmatrix} \lambda_1^{-1} P_A \mathcal{A} \mathcal{A}^T P_A & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & \lambda_1 E_A^T E_A & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad (46) $$

$$ \mathcal{L}_{B1}(k) + \mathcal{L}_{B1}^T(k) \le \begin{bmatrix} \lambda_3^{-1} P_A \mathcal{B} \mathcal{B}^T P_A & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \lambda_3 E_B^T E_B & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad (47) $$
Similarly,

$$ \mathcal{L}_{A2}(k) + \mathcal{L}_{A2}^{T}(k) \le \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \begin{bmatrix} \lambda_2^{-1} P_0 A A^T P_0 & 0 \\ 0 & 0 \end{bmatrix} \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & \begin{bmatrix} \lambda_2 E_A^T E_A & 0 \\ 0 & 0 \end{bmatrix} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \quad (48) $$

$$ \mathcal{L}_{B2}(k) + \mathcal{L}_{B2}^{T}(k) \le \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \begin{bmatrix} \lambda_4^{-1} P_0 B B^T P_0 & 0 \\ 0 & 0 \end{bmatrix} \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & \lambda_4 E_B^T E_B & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \quad (49) $$

The zeros inside the matrices are of different sizes and so the position of the nonzero elements are not identical as visible. The overall bound time varying would be,

$$ \mathcal{L}_{ij}(k) + \mathcal{L}_{ij}^{T}(k) \le \begin{pmatrix} \begin{bmatrix} \mathcal{L}^{1} & 0 \\ 0 & 0 \end{bmatrix} & 0 & 0 & 0 \\ 0 & \mathcal{L}^{2} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \begin{bmatrix} \mathcal{L}^{3} & 0 \\ 0 & 0 \end{bmatrix} \end{pmatrix} \quad (50) $$

with,

$$ \mathcal{L}^{1} = \lambda_{1}^{-1} P_{A} A A^{T} P_{A} + \lambda_{3}^{-1} P_{A} B B^{T} P_{A}, \quad (51) $$

$$ \mathcal{L}^2 = \begin{bmatrix} (\lambda_1 + \lambda_2)E_A^T E_A & 0 & 0 & 0 \\ 0 & (\lambda_3 + \lambda_4)E_B^T E_B & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad (52) $$

$$ \mathcal{L}^3 = \lambda_2^{-1} P_0 A A^T P_0 + \lambda_4^{-1} P_0 B B^T P_0 \quad (53) $$

This would mean, we could rewrite the matrices of (38) with bounds as ($\forall i, j$),

$$ Q_{ij} + \mathcal{L}_{ij}(k) + \mathcal{L}_{ij}^{T}(k) \le Q_{m} \quad (54) $$

$$ Q_m = \begin{bmatrix} \begin{pmatrix} -P_0 + Q_{e0} + \mathcal{L}^1 & 0 \\ 0 & -P_1 + Q_{e1} \end{pmatrix} & Q_A & \Phi_{ij}^T P & 0 \\ Q_A^T & T_{22} & 0 & Q_B^T \\ P \Phi_{ij} & 0 & \begin{pmatrix} -P_0 & 0 \\ 0 & -P_1 \end{pmatrix} & 0 \\ 0 & Q_B^T & 0 & \begin{pmatrix} -P_0 + \mathcal{L}^3 & 0 \\ 0 & -P_1 \end{pmatrix} \end{bmatrix} \quad (55) $$

Here, $Q_e \triangleq \text{diag}(Q_{e0}, Q_{e1})$ and

$$ T_{22} = \begin{bmatrix} -\Gamma_2^x + (\lambda_1 + \lambda_2)E_A^T E_A & 0 & 0 & 0 \\ 0 & -\Gamma_2^u + (\lambda_3 + \lambda_4)E_B^T E_B & 0 & 0 \\ 0 & 0 & -\Gamma_2^\theta & 0 \\ 0 & 0 & 0 & -\Gamma_2^{\Delta\theta} \end{bmatrix} \quad (56) $$
The nonlinear terms of $\mathcal{L}^1$ and $\mathcal{L}^3$ could be converted to linear terms using Schur's complement (Lemma 2.4). The two terms in each of these two entries would contribute to four matrices, that is,

$$ \lambda_1^{-1} P_A \mathcal{A} \mathcal{A}^T P_A < 0 \Rightarrow 0 - (P_A \mathcal{A}) \lambda_1^{-1} (\mathcal{A}^T P_A) < 0 \Rightarrow \begin{bmatrix} 0 & P_A \mathcal{A} \\ \mathcal{A}^T P_A & -\lambda_1 I \end{bmatrix} \quad (57) $$

The same steps would lead to the other three nonlinear terms to give an overall matrix condition such that,

$$ \left[ \begin{array}{c} \begin{pmatrix} -P_0 + Q_{e0} & 0 \\ 0 & -P_1 + Q_{e1} \end{pmatrix} \begin{array}{c} Q_A \\ Q_A^T \end{array} \begin{pmatrix} \Phi_{ij}^T P \\ 0 \end{pmatrix} \begin{array}{c} 0 \\ Q_B^T \end{array} \end{pmatrix} \\ \begin{pmatrix} P \Phi_{ij} & 0 \\ 0 & 0 \end{pmatrix} \begin{array}{c} T_{22} \\ 0 \end{array} \begin{pmatrix} -P_0 & 0 \\ 0 & -P_1 \end{pmatrix} \begin{array}{c} 0 \\ 0 \end{array} \end{pmatrix} \\ \begin{pmatrix} A^T P_A & 0 & 0 & 0 \\ B^T P_A & 0 & 0 & 0 \\ 0 & 0 & 0 & A^T P_0 \\ 0 & 0 & 0 & B^T P_0 \end{pmatrix} \end{array} \begin{bmatrix} P_A A & P_A B & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & P_0 A & P_0 B \end{bmatrix} \right] \quad (58) $$

or in a simplified form, it leads to (19). Hence the proof. ■

## 4. Simulation Example

Consider a simple second order nonlinear system,

$$ \begin{aligned} x_1(k+1) &= -0.5x_1^2(k) + 0.5x_2(k) + x_1(k)\theta(k) \\ x_2(k+1) &= 0.8x_1(k) + (1 - 0.5\theta)u(k) \\ y(k) &= x_1(k) \end{aligned} \quad (59) $$

Considering that the state $x_1(k) \in [0, 1]$ and the $\theta(k) \in [-1, 1]$ and consider $z \triangleq x_1$ as the premise variables, the system matrices for (13) are given by,

$$ A_{11} = \begin{bmatrix} -1 & 0.5 \\ 0.8 & 0 \end{bmatrix}, \quad A_{12} = \begin{bmatrix} 1 & 0.5 \\ 0.8 & 0 \end{bmatrix}, \quad A_{21} = \begin{bmatrix} -1.5 & 0.5 \\ 0.8 & 0 \end{bmatrix}, \quad A_{11} = \begin{bmatrix} 0.5 & 0.5 \\ 0.8 & 0 \end{bmatrix} \quad (60) $$

and with

$$ B_{11} = B_{13} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \quad B_{21} = B_{22} = \begin{bmatrix} 0 \\ 0.5 \end{bmatrix} \quad (61) $$

The weighting functions are given by,

$$ \mu_1(\hat{x}) = 1 - \hat{x}, \quad \mu_2(\hat{x}) = \hat{x}, \quad \tilde{\mu}_1(\hat{\theta}) = \frac{1 - \hat{\theta}}{2}, \quad \tilde{\mu}_2(\hat{\theta}) = \frac{\hat{\theta} + 1}{2} \quad (62) $$

For developing the observer, the aspects were considered:

* The value of $\eta$ was fixed at -0.995. Apart from avoiding the nonlinear matrix term, this would also avoid $\Phi_{ij}$ to become marginal stable in (30).

* A condition to ensure that the value of $K_{ij}$ are sufficiently larger than that of $\eta$ as discussed in [13], the following LMI condition was introduced

$$ F_{ij} > P_1 \eta $$
<figure><img src="image_5.png" /><figcaption>Figure 1. States and the parameter of the system and their estimates</figcaption></figure>

* The chosen values for $Q_e = diag(Q_{e0}Q_{e1})$ as in (34) are given by

$$Q_{e0} = \begin{bmatrix} 0.001 & 0 \\ 0 & 0.001 \end{bmatrix}, Q_{e1} = 0.001$$

The resultant observer gains were as follows:

$$L_{11} = \begin{bmatrix} -1 \\ 0.8 \end{bmatrix}, L_{12} = \begin{bmatrix} 1 \\ 0.8 \end{bmatrix}, L_{21} = \begin{bmatrix} -1.5 \\ 0.8 \end{bmatrix}, L_{22} = \begin{bmatrix} 0.5 \\ 0.8 \end{bmatrix} \quad \text{and} \quad K_{ij} = 1.014 \ \forall i, j \quad (63)$$

The state and the parameter estimation results are shown in the Fig. 1. As could be seen, when there is an error in the state estimates, the parameter estimation component does not take corrective actions. Only when the state estimates settle down does the parameter estimation component kicks in, which is visible in the form of the delay in the parameter estimation. The input used for the simulation as well as the variation of the four weighting functions during the simulation are shown in the Fig.2.

## 5. Concluding Remarks

This paper provided the results for a joint state and parameter estimation observer for a discrete time Takagi Sugeno models. The techniques illustrated to handle the nonlinear matrix inequalities and time varying terms could be useful in similar discrete-time version implementation of other results. A proof for the results were given as well and a simple academic example illustrated the results. In future, one exploration could be to use piece-wise Lyapunov
<figure><img src="image_6.png" /><figcaption>**Figure 2.** States and the parameter of the system and their estimates</figcaption></figure>

function to avoid the problems with the nonlinearity of matrix inequalities. Further, these results could be applied to some system models of practical systems.

## References

[1] Bezzaoucha S, Marx B, Maquin D and Ragot J 2013 State and parameter estimation for nonlinear systems: a Takagi-Sugeno approach *American Control Conference, ACC 2013, Washington, DC, USA, June 17-19*

[2] Bergsten P, Palm R and Driankov D 2001 Fuzzy observers. *FUZZ-IEEE* pp 700–703

[3] Lendek Z, Babuška R, De Schutter B et al. 2009 *IEEE Transactions on Fuzzy Systems* **17** 641–653

[4] Lendek Z, Babuška R and De Schutter B 2011 *Fuzzy Sets and Systems* **174** 1–30

[5] Ichalal D, Marx B, Ragot J and Maquin D 2009 State and unknown input estimation for nonlinear systems described by Takagi-Sugeno models with unmeasurable premise variables 17th *Mediterranean Conference on Control and Automation MED* (IEEE) pp 217–222

[6] Bezzaoucha S, Marx B, Maquin D and Ragot J 2013 State and parameter estimation for time-varying systems: a Takagi-Sugeno approach 5th *Symposium on System Structure and Control, IFAC Joint Conference 2013 SSSC, TDS, FDA*

[7] Tanaka K and Wang H O 2004 *Fuzzy control systems design and analysis: a linear matrix inequality approach* (John Wiley & Sons)

[8] Kruszewski A 2006 *Lois de commande pour une classe de modèles non linéaires sous la forme Takagi-Sugeno: Mise sous forme LMI* Ph.D. thesis Université de Valenciennes et du Hainaut-Cambresis

[9] Zhou K and Khargonekar P P 1988 *Systems & Control Letters* **10** 17–20

[10] de Souza C E and Xie L 1992 *Systems & Control Letters* **18** 61–71

[11] Srinivasarengan K, Ragot J, Maquin D and Aubrun C 2016 Nonlinear joint state-parameter observer for vav damper position estimation 3rd *Conference on Control and Fault-Tolerant Systems, SysTol 2016*

[12] Blanco Y 2001 *STABILISATION DES MODELES TAKAGI-SUGENO ET LEUR USAGE POUR LA COMMANDE DES SYSTEMES NON LINEAIRES* Ph.D. thesis Université des Sciences et Technologies de Lille

[13] Srinivasarengan K, Ragot J, Maquin D and Aubrun C 2016 Takagi-sugeno model based nonlinear parameter estimation in air handling units 4th *IFAC International Conference on Intelligent Control and Automation Sciences, ICONS 2016*