Analyzing Uncertain Dynamical Systems after State-Space Transformations into Cooperative Form: Verification of Control and Fault Diagnosis

Abstract

When modeling real-life applications, uncertainty appears in the form of, for example, modeling approximations, measurement errors, or simply physical restrictions. Those uncertainties can either be treated as stochastic or as bounded, with known limits in the form of intervals. The latter is considered in this paper for a real-life application in the form of an electrical circuit. This is reasonable because the electrical circuit is subject to uncertainties, mainly due to circuit element tolerances and variable load conditions. Since conservative worst-case limits for such parameters are commonly known, interval methods can be applied. The aim of this paper is to demonstrate a possible overall handling of the given uncertain system. Firstly, this includes a control and a reliable computation of the states’ interval enclosures. On the one hand, this can be used to predict the system’s behavior, and on the other hand to verify the control numerically. Here, the implemented feedback control is based on linear matrix inequalities (LMIs) and the states are predicted using an interval enclosure technique based on cooperativity. Since the original system is not cooperative, a transformation is performed. Finally, an observer is implemented as a diagnosis tool regarding faulty measurements or component failures. Since adding a state-of-the-art observer would destroy this structure, a cooperativity-preserving method is applied. Overall, this paper combines methods from robust control design and interval-based evaluations, and presents a suitable observer technique to show the applicability of the presented methods.

1. Introduction

There are different reasons for the occurrence of uncertainty. It may appear due to model simplifications, the approximation of nonlinearities, imprecise parameter knowledge and/or order reduction, as well as physical and numerical restrictions of the system itself. Uncertainty caused by measurement noise and sensor inaccuracies are further examples. In any case, uncertainties can be treated either stochastically or as bounded quantities in terms of worst-case scenarios, where the lower and upper bounds are summarized in scalar or multi-dimensional intervals. The considered application scenario in this paper is an uncertain passive, second-order electric network as a first glance for these types of systems. As the reader will see, the presented theoretical aspects seem to be difficult to apply in real life. However, this paper aims to show their applicability and discuss their limitations. An extension to the presented methods by introducing strategies for the multi-sectioning of uncertain parameter domains as well as an application to a step-down converter can be found in [1]. In the presented electrical system, in preparation for the extension to a step-down converter, the load resistance can vary in a wide range and is treated as the main source of uncertainty. However, since the bounds are known, interval arithmetic is used (see [2]). We will treat the fundamental electrical circuit with a full control-oriented approach, which generally includes a controller and an observer. Since our system is uncertain, those steps are not as straightforward as for exactly known systems. Hence, the first step is to find a suitable controller gain for the uncertain system with an LMI-based method. To verify the control (i.e., to verify its resulting stable behavior),we want to predict the resulting state reliably despite the given uncertainties. Normally, this would be done by a Picard iteration with a subsequent tightening step evaluating a temporal Taylor series expansion of the initial value problem (IVP) [3,4]. However, this tends to lead to overestimation [2] due to the so-called wrapping effect. An alternative solution avoiding this could be to use cooperativity, which has already been investigated in several papers [5,6,7]. Its advantage is the simplification of several tasks, that is, the computation of guaranteed state enclosures as well as forecasting worst-case bounds for selected system outputs in feedback and predictive control [8], the identification of unknown parameters [9,10], and state prediction with the aim of fault diagnosis [11]. If a system is cooperative, the worst-case bounds of the trajectories of the uncertain system can be computed following the element-wise inequalities (1) as two separate systems with point-valued parameters and states, while assuring that all possible states lie within their solutions

$${\mathbf{f}}_v\left(\mathbf{v}\right)=\dot{\mathbf{v}}\left(t\right)\le \dot{\mathbf{x}}\left(t\right)\le \dot{\mathbf{w}}\left(t\right)={\mathbf{f}}_w\left(\mathbf{w}\right)\mathrm{with}\mathbf{x}\in \left[\mathbf{v};\mathbf{w}\right].$$

(1)

A system is cooperative (as a sufficient criterion) if for an autonomous dynamic system

$$\dot{\mathbf{x}}\left(t\right)=\mathbf{f}\left(\mathbf{x}\left(t\right)\right),\mathbf{x}\in {\mathbb{R}}^n,$$

(2)

all off-diagonal elements $J_{i,j}$, $i,j\in \{1,\dots ,n\}$, $i\ne j$, of the corresponding Jacobian

$$\mathbf{J}=\frac{\partial \mathbf{f}\left(\mathbf{x}\right)}{\partial \mathbf{x}}$$

(3)

are non-negative according to

$$J_{i,j}\ge 0,i,j\in \{1,\dots ,n\},i\ne j.$$

(4)

Matrices with this structure are also called Metzler. For those, state trajectories $\mathbf{x}\left(t\right)$ starting in the positive orthant

$${\mathbb{R}}_{+}^n=\{\mathbf{x}\in {\mathbb{R}}^n\mid x_i\ge 0, \forall i\in \{1,\dots ,n\}\}$$

(5)

are guaranteed to stay in this positive orthant for all $t\ge 0$ because

$${\dot{x}}_i\left(t\right)=f_i\left(x_1,\dots ,x_{i-1},0,x_{i+1},\dots x_n\right)\ge 0$$

(6)

holds for all components $i\in \{1,\dots ,n\}$ of the state vector as soon as the state variable $x_i$ reaches the value $x_i=0$. This property is often referred to as positivity of the system model (2) [12]. There are system models that are naturally cooperative, like in the fields of biological, chemical, and medical applications (e.g., compartmental models in epidemiology). However, the presented application scenario does not show this property when derived by first-principles techniques. To use the advantages of the decoupled bounding systems in (1) regardless, a transformation into cooperative form is employed based on the findings in [6]. Note that the transformation approach used is only suitable for systems with purely real eigenvalues, which has to be considered especially when deriving the control. A combination of the methods with the sequence of application and their relations is shown in Figure 1.

Firstly, a control design is applied to the uncertain system. The resulting model is then transformed into a cooperative state-space representation with which verified interval enclosures can be computed. These can be used to tune the control design, further optimizing the robust control approach. Finally, a cooperativity-ensuring observer design can be applied as a form of fault diagnosis tool, further securing the overall system dynamics to follow the requirements set by the user. Note that throughout this paper, the general method is applied step by step (following Figure 1) to the electrical circuit for a better understanding. For this purpose, Section 2 derives the model of the fundamental electrical network by first-principles techniques. Section 3 first presents a general control design including uncertainties and using LMIs and then shows the results of applying this design to the electrical circuit. A transformation into cooperative form follows in Section 4, with the same order of first giving the theory and presenting numerical results as a conclusion. The simulation results of the presented methods can be found in Section 5. An observer is added as a tool for, for example, fault diagnosis of the measurements used in the control. However, a cooperativity-preserving approach is applied in Section 6 to maintain the structure from before. Note that a further paper [13] presents the possibility of including a cooperativity-preserving controller, which could be applied to the given example in rearranging the order of implementation into transformation, control, and observer. Finally, Section 7 gives a conclusion and an outlook on future work.

Note that throughout this paper, intervals are denoted by square brackets, such as

$$\mathbf{x}\in \left[\mathbf{x}\right]=\left[\underline{\mathbf{x}};\overline{\mathbf{x}}\right],$$

(7)

where, for example, $\mathbf{x}$ is a point-valued state vector and $\left[\mathbf{x}\right]$ its interval-valued enclosure with $\underline{\mathbf{x}}$ and $\overline{\mathbf{x}}$ as the lower and upper bounds, respectively. Matrices and scalars in interval notation are given accordingly.

2. Modeling of an Electrical Circuit

To model the system’s dynamical behavior mathematically, the simplified RLC model shown in Figure 2 is used, including the parasitic inner resistances $R_{\mathrm{L}}$ and $R_{\mathrm{C}}$ for the inductivity and the capacity, respectively.

Here, ${\check{R}}_{\mathrm{L}}=R_0+R_{\mathrm{L}}$ combines a limiting resistance $R_0$ and the inner resistance $R_{\mathrm{L}}$ of the inductivity. Furthermore,$R_{\mathrm{S}}={\tilde{R}}_{\mathrm{S}}+\Delta R_{\mathrm{S}}$ is a summed-up variable load, where $\Delta R_{\mathrm{S}}$ is implemented as a series connection of various resistances that can be activated and deactivated by semiconductive switches while ${\tilde{R}}_{\mathrm{S}}$ is always present to avoid a short circuit at the system’s output terminals. Two voltage loops are described by the mesh equations

$$u_{\mathrm{E}}=u_{{\check{R}}_{\mathrm{L}}}+u_{\mathrm{L}}+u_{R_{\mathrm{C}}}+u_{\mathrm{C}}$$

(8)

and

$$u_{\mathrm{C}}+u_{R_{\mathrm{C}}}=u_{R_{\mathrm{S}}},$$

(9)

and the Kirchhoff’s node equation gives

$$i_{R_{\mathrm{S}}}+i_{\mathrm{C}}=i_{\mathrm{L}}.$$

(10)

All ohmic resistances are governed by the component equations

$$u_{R_i}=R_i·i_{R_i}$$

(11)

with $i\in \{\mathrm{L},\mathrm{S},\mathrm{C}\}$, while the inductivity and capacity are represented by

$$u_{\mathrm{L}}=L·{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_{\mathrm{L}}$$

(12)

and

$$i_{\mathrm{C}}=C·{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{u}_{\mathrm{C}},$$

(13)

respectively. Variations of the magnetic field energy are characterized by the first-order ordinary differential equation

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_{\mathrm{L}}={\displaystyle \frac{1}{L}}\left[-\left({\check{R}}_{\mathrm{L}}+{\displaystyle \frac{R_{\mathrm{S}}{R}_{\mathrm{C}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}\right)·i_{\mathrm{L}}-\left(1-{\displaystyle \frac{R_{\mathrm{C}}}{\left(R_{\mathrm{S}}+R_{\mathrm{C}}\right)}}\right)·u_{\mathrm{C}}+u_{\mathrm{E}}\right]$$

(14)

and changes of the electric field energy by

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{u}_{\mathrm{C}}={\displaystyle \frac{R_{\mathrm{S}}}{C\left(R_{\mathrm{S}}+R_{\mathrm{C}}\right)}}·i_{\mathrm{L}}-{\displaystyle \frac{1}{C\left(R_{\mathrm{S}}+R_{\mathrm{C}}\right)}}·u_{\mathrm{C}}.$$

(15)

Since the applied theories in this paper are based on the state-space representation, the modeling is finalized by deriving it. Both types of energy in the derivation of (14) and (15) represent the physical storage expressed by the state variables $x_1=i_{\mathrm{L}}$ and $x_2=u_{\mathrm{C}}$. The resulting state-space representation becomes

$$\dot{\mathbf{x}}=\left[\begin{array}{cc}-{\displaystyle \frac{1}{L}}\left({\check{R}}_{\mathrm{L}}+{\displaystyle \frac{R_{\mathrm{S}}{R}_{\mathrm{C}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}\right)& {\displaystyle \frac{1}{L}}\left({\displaystyle \frac{R_{\mathrm{C}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}-1\right)\\ {\displaystyle \frac{R_{\mathrm{S}}}{C\left(R_{\mathrm{S}}+R_{\mathrm{C}}\right)}}& -{\displaystyle \frac{1}{C\left(R_{\mathrm{S}}+R_{\mathrm{C}}\right)}}\end{array}\right]·\mathbf{x}+\left[\begin{array}{c}{\displaystyle \frac{1}{L}}\\ 0\end{array}\right]·u_{\mathrm{E}}$$

(16)

with the output equation

$$y=u_{R_{\mathrm{C}}}+u_{\mathrm{C}}=\left[\begin{array}{cc}{\displaystyle \frac{R_{\mathrm{S}}{R}_{\mathrm{C}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}& 1-{\displaystyle \frac{R_{\mathrm{C}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}\end{array}\right]·\mathbf{x}$$

(17)

according to Figure 2. This output is a reasonable way to determine the power transported to the consumer $R_{\mathrm{S}}$ and is given here to complete the state-space representation. However, it is not used explicitly in the applied theory because the state variables are assumed to be measurable. 

Table 1. Parameters of the low-power electrical circuit.

Variable	Unit/Value	Meaning
L	$1\mathrm{H}$	Inductance
C	$2\mathrm{mF}$	Capacity
$R_{\mathrm{S}}$	$\left[0.1;3\right]\Omega$	Ohmic resistance of the load
$R_{\mathrm{C}}$	$\left[0.1;0.6\right]\Omega$	Ohmic resistance of the capacity
${\check{R}}_{\mathrm{L}}$	$100\Omega$	Overall ohmic resistance of the inner resistance of the inductivity and a constant limiting resistance

outlines the chosen parameters, including the uncertain resistances of $R_{\mathrm{S}}$ and $R_{\mathrm{C}}$.

Here, the inductance is implemented with the help of a gyrator circuit that turns capacitors into virtual inductors with non-zero internal resistances [14], which is a common approach to apply such large inductance values in a low-power network. Note further that the ohmic resistances are given in interval representation because of their uncertainty, which is due to their variability. The resulting eigenvalues of the system matrix,which appears in Equation (16), are given in Figure 3 with respect to the variable load resistance. For the presented example, the system was evaluated for both $R_{\mathrm{C}}=0.1\Omega$ (blue line) and $R_{\mathrm{C}}=0.6\Omega$ (black line), resulting in purely real eigenvalues, which is to be kept for the controlled system in order to successfully apply the transformation into cooperativity later on. All values of $R_{\mathrm{C}}$ between the two limits would lie in this range, and are thus not visualized.

3. Robust State-Feedback Control

Assuming that all the state variables can be measured, a control is designed to stabilize the system with a vanishing reference signal. In general, a linear state feedback controller for a point-valued system is derived according to

$$\mathbf{u}=-\mathbf{Kx}$$

(18)

with a constant gain $\mathbf{K}$. Then, the closed loop system can be described by

$$\dot{\mathbf{x}}=\mathbf{Ax}+\mathbf{Bu}=\mathbf{Ax}-\mathbf{BKx}={\mathbf{A}}_{\mathrm{C}}\mathbf{x},$$

(19)

where

$${\mathbf{A}}_{\mathrm{C}}=\left(\mathbf{A}-\mathbf{BK}\right)$$

(20)

holds. For a robust stabilization of the system dynamics despite the given uncertainties, the given intervals have to be taken into account. For this, linear matrix inequality techniques are well suited. Here, we make use of a polytopic uncertainty representation [15], with a parameter-dependent system matrix $\mathbf{A}\left(\mathbf{p}\right)$ which overapproximates the system model (16) to achieve robustness of the control design. This results in the convex combination of extremal system realizations of

$$\begin{array}{cc}\hfill \mathcal{D}=\{& \left[\mathbf{A}\left(\mathit{\xi }\right),\mathbf{B}\left(\mathit{\xi }\right)\right]|\left[\mathbf{A}\left(\mathit{\xi }\right),\mathbf{B}\left(\mathit{\xi }\right)\right]\hfill \\ \hfill & =\sum _{\nu =1}^{n_{\nu }}{\xi }_v·\left[{\mathbf{A}}_{\nu },{\mathbf{B}}_{\nu }\right];\sum _{\nu =1}^{n_{\nu }}{\xi }_{\nu }=1;{\xi }_{\nu }\ge 0\},\hfill \end{array}$$

(21)

where the vertex matrices are denoted by ${\mathbf{A}}_{\nu }={\mathbf{A}}_{\nu }\left(\mathbf{p}\right)$ as well as ${\mathbf{B}}_{\nu }={\mathbf{B}}_{\nu }\left(\mathbf{p}\right)$, with each of them depending on the vector of independent parameters $\mathbf{p}\in {\mathbb{R}}^{n_{\mathrm{p}}}$ in an affine manner. Those independent parameters are contained in the interval box $\left[\mathbf{p}\right]=\left[\underline{\mathbf{p}};\overline{\mathbf{p}}\right]$ with the component-wise defined bounds ${\underline{p}}_i\le p_i\le {\overline{p}}_i$, $i\in \{1,\dots ,n_{\mathrm{p}}\}$. In general, the evaluation of $\mathbf{A}\left(\mathbf{p}\right)$ and $\mathbf{B}\left(\mathbf{p}\right)$ is performed for each of the vertices

$$\mathcal{P}=\left\{\left[\begin{array}{c}{\underline{p}}_1\\ {\underline{p}}_2\\ ⋮\\ {\underline{p}}_{n_{\mathrm{p}}}\end{array}\right],\left[\begin{array}{c}{\overline{p}}_1\\ {\underline{p}}_2\\ ⋮\\ {\underline{p}}_{n_{\mathrm{p}}}\end{array}\right],\dots ,\left[\begin{array}{c}{\overline{p}}_1\\ {\overline{p}}_2\\ ⋮\\ {\overline{p}}_{n_{\mathrm{p}}}\end{array}\right]\right\}=\left\{{\mathbf{p}}^{\langle 1\rangle },\dots ,{\mathbf{p}}^{\langle n_{\nu }\rangle }\right\}.$$

(22)

For the presented model with the two uncertain parameters $R_{\mathrm{S}}$ and $R_{\mathrm{C}}$, the domain $\mathcal{D}$ according to (21) can be parameterized by $n_{\nu }=4$ extremal realizations that need to be considered for the robust control design. Another advantage—in addition to the inclusion of the uncertainty of the LMI-based approach—is the inclusion of a so-called $\Gamma$-region of desired eigenvalue locations as well as optimality criteria such as robust $H_2$ and $H_{\infty }$ tasks [16]. In the presented case, we will make use of the former to enhance the closed-loop performance with feasible eigenvalue regions defined by domains of strict negative definiteness of

$${\mathbf{F}}_{\Gamma }\left(s\right)={\mathbf{D}}_0+s{\mathbf{D}}_1+\overline{s}{\mathbf{D}}_1^T\prec 0,$$

(23)

see [15]. Note that in general all LMIs have to be strictly positive or negative in order to compute them properly [17]. Here, the Laplace variable $s\in \mathbb{C}$ corresponds to the set of all eigenvalues if (23) is applied to the linear system model (16), and its conjugate complex $\overline{s}$. This must hold for all eigenvalues of the closed-loop system.

In (23), the real-valued parameter matrices ${\mathbf{D}}_0={\mathbf{D}}_0^T$ and ${\mathbf{D}}_1$ provide a certain flexibility to define $\Gamma$-stability regions such as ellipses, hyperbolas, parabolas, cones, and strips in the complex plane [18]. For a direct implementation into the computation of the controller gain, we envisage a computationally tractable solution. Therefore, the inequality ${\mathbf{F}}_{\Gamma }\prec 0$ is generally reformulated into an equivalent LMI

$${\mathbf{D}}_0\otimes \mathbf{P}+{\mathbf{D}}_1\otimes \left(\mathbf{AP}\right)+{\mathbf{D}}_1^T\otimes {\left(\mathbf{AP}\right)}^T\prec 0,$$

(24)

where ⊗ denotes the Kronecker product, applying a Lyapunov design according to [15] for any point-valued system matrix $\mathbf{A}$. If all eigenvalues of a real-valued system matrix $\mathbf{A}$ lie within the interior of the region (23), a positive definite matrix $\mathbf{P}={\mathbf{P}}^T\succ 0$ exists that fulfills the matrix inequality [15]. Here, $\mathbf{P}$ is a matrix that defines a Lyapunov function $V\left(\mathbf{x}\right)=\frac{1}{2}{\mathbf{x}}^T\mathbf{P}\mathbf{x}>0$ for $\mathbf{x}\ne {\mathbf{x}}_{\mathrm{s}}=\mathbf{0}$, where ${\mathbf{x}}_{\mathrm{s}}$ represents a stationary state, with which the stability of the dynamic system $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$ can be proven. Now, ${\mathbf{D}}_0={\mathbf{D}}_0^T$ and ${\mathbf{D}}_1$ can be used to specify the users’ needs concerning different stability regions as mentioned above. Pure Hurwitz stability can be included by setting $\gamma =0$ in $F_{\Gamma }=2\gamma +s+\overline{s}\prec 0$, which corresponds to real-valued scalars $D_0=2\gamma$ and $D_1=1$due to its simplicity. Note that an absolute stability margin can be guaranteed by choosing $\gamma >0$, which leads to $\Re \left\{s\right\}<-\gamma$ [16]. In the presented application, the eigenvalues shall be forced to stay purely real as a requirement for the transformation procedure. Hence, a sufficiently large damping ratio is introduced by setting

$${\mathbf{D}}_0=\mathbf{0}\mathrm{and}{\mathbf{D}}_1=\left[\begin{array}{cc}sin\left(\theta \right)& cos\left(\theta \right)\\ -cos\left(\theta \right)& sin\left(\theta \right)\end{array}\right]$$

(25)

with a successive, automatic minimization of $\theta$ up to $\theta ={\displaystyle \frac{\pi }{8}}$. In the sequel, it is shown that the resulting eigenvalues of the considered system are purely real. Now, a state-of-the-art Lyapunov approach is applied directly for the control design of the uncertain system in the form of a sufficient condition. Here, the system matrix in Equation (24) is replaced by the expression for a controlled system matrix from Equation (20), which under consideration of the convexity of the domain (21) leads to a reformulation of Equation (24) in the design LMIs

$${\mathbf{D}}_0\otimes \mathbf{Q}+{\mathbf{D}}_1\otimes (\mathbf{Q}{\mathbf{A}}_{\nu }^T-{\mathbf{Y}}^T{\mathbf{B}}_{\nu }^T)+{\mathbf{D}}_1^T\otimes ({\mathbf{A}}_{\nu }\mathbf{Q}-{\mathbf{B}}_{\nu }\mathbf{Y})\prec 0.$$

(26)

Note that this final control design formulation does include all extremal realizations with $\nu \in \{1,\dots ,2^{n_{\mathrm{p}}}\}$.This means that robust stability for the uncertainty representation (21) and (22) with eigenvalues that are compatible with the domain ${\mathbf{F}}_{\Gamma }\prec 0$ defined in (23) is achieved as soon as a joint solution $\mathbf{Q}\succ 0$, $\mathbf{Y}$ of the LMI (26) has been found that is valid for each of the vertices of index $\nu$. The linearizing change of variables is reverted by $\mathbf{P}={\mathbf{Q}}^{-1}$ and, generally, $\mathbf{K}=\mathbf{Y}\mathbf{P}$, which becomes ${\mathbf{k}}^T={\mathbf{y}}^T\mathbf{P}$ in the given single-input case. The controller gains

$${\mathbf{k}}^T=\left[\begin{array}{cc}1.211& 0.291\end{array}\right]\times {10}^4$$

(27)

are obtained by numerically solving the problem using well-known LMI toolboxes like the semidefinite programming solver SeDuMi [19] in combination with the modeling and optimization toolbox YALMIP [17] for MATLAB. The controlled system matrix is enclosed by

$${\left[\mathbf{A}\right]}_{\mathrm{C}}=\left[\begin{array}{cc}\left[-1.223;-1.221\right]& \left[-0.291;-0.290\right]\\ \left[0.001;0.751\right]& \left[-0.251;-0.013\right]\end{array}\right]\times {10}^4$$

(28)

with the resulting eigenvalue locations shown in Figure 4, which were calculated by gridding $R_S$ and $R_C$. Here, one can clearly see distinct eigenvalue regions for the two eigenvalues of the considered system as well as the fact that they stay purely real.

For a matrix to be Metzler, all off-diagonal entries need to be non-negative; thus, the presented system matrix is obviously not Metzler because it has a strictly negative off-diagonal entry. Hence, a transformation into cooperative form is performed in the following to validate the control design in the time domain.

4. Transformation into a Cooperative Form

A generalized method to transform time-varying systems and systems with uncertainties containing purely real eigenvalues into cooperative form was presented in [5]. This was redesigned into an optimization task to solve it computationally in [6]. As we ensured that our system would keep purely real eigenvalues, this transformation is used. A short summary of the method is given here, followed by the numerical results for our given system. Firstly, the requirements published in [5] are given. It is assumed that the uncertain system matrix can be expressed by the element-wise defined inequality

$${\mathbf{Z}}_{\mathrm{a}}-\Delta \le \mathbf{Z}:={\mathbf{A}}_{\mathrm{C}}\le {\mathbf{Z}}_{\mathrm{a}}+\Delta .$$

(29)

Here, $\Delta$ consists of the (symmetric) worst-case bounds of all entries in the interval-valued matrix ${\left[\mathbf{A}\right]}_{\mathrm{C}}$, which results in a symmetric midpoint matrix ${\mathbf{Z}}_{\mathrm{a}}={\mathbf{Z}}_{\mathrm{a}}^T$ in (29). Now, we need to search for a Metzler matrix $\mathbf{R}=\mu {\mathbf{E}}_n-\mathsf{\Gamma }$, which has the same eigenvalues as ${\mathbf{Z}}_{\mathrm{a}}$, with a constant $\mu \in \mathbb{R}$ and a diagonal matrix $\mathsf{\Gamma }\in {\mathbb{R}}^{n\times n}$; ${\mathbf{E}}_n\in {\mathbb{R}}^{n\times n}$ is a matrix with all elements equal to 1 and the identity matrix $\mathbf{I}$ is of order n. According to [5], if $\mathrm{eig}\left(\mathbf{R}\right)=\mathrm{eig}\left({\mathbf{Z}}_{\mathrm{a}}\right)$, there exists an orthogonal matrix $\mathbf{S}\in {\mathbb{R}}^{n\times n}$ such that ${\mathbf{S}}^T\mathbf{ZS}$ is Metzler provided that ${\mu >n\left|\right|\Delta \left|\right|}_{\max }$, where ${\left|\right|\Delta \left|\right|}_{\max }$ denotes the maximum, element-wise, absolute value of $\Delta$.

The next part sums up the findings of [6] and its general formulation of the computationally feasible optimization problem formulated with LMI constraints. To find ${\mathbf{Z}}_{\mathrm{a}}$ and $\Delta$, two cases are distinguished. If the system matrix is diagonally dominant, ${\mathbf{Z}}_{\mathrm{a}}$ is chosen to represent the diagonal entries of the original system matrix. On the other hand, if the system matrix is not diagonally dominant initially (as is the case for the presented application scenario), the element-wise defined interval midpoint matrix $\mathrm{mid}\left\{{\left[\mathbf{A}\right]}_{\mathrm{C}}\right\}$ is transformed into a diagonal structure (except for numerical round-off errors) by determining a new system matrix

$${\widehat{\mathbf{A}}}_{\mathrm{C}}={\mathbf{V}}^{-1}{\mathbf{A}}_{\mathrm{C}}\mathbf{V}.$$

(30)

Here, $\mathbf{V}$ is defined by the floating-point approximation of the n linearly independent real-valued eigenvectors of the matrix $\mathrm{mid}\left\{{\left[\mathbf{A}\right]}_{\mathrm{C}}\right\}$ and ${\mathbf{Z}}_{\mathrm{a}}$ is set to be a diagonal matrix with the asymptotically stable, real eigenvalues of $\mathrm{mid}\left\{{\left[\mathbf{A}\right]}_{\mathrm{C}}\right\}$. Note that interval arithmetic software libraries such as INTLAB are used to handle round-off errors in the matrix inversion in (30) [20]. The required worst-case bounds in $\Delta$ are chosen as

$$\delta =\max \left({|\left[\mathbf{A}\right]}_{\mathrm{C}}-{\mathbf{Z}}_{\mathrm{a}}|\right)\mathrm{or}\delta =\max \left(|{\left[\widehat{\mathbf{A}}\right]}_{\mathrm{C}}-{\mathbf{Z}}_{\mathrm{a}}|\right),$$

(31)

respectively, with $\Delta =\delta ·{\mathbf{E}}_n$ and a maximization that is carried out over all matrix entries after determining their absolute values in an element-by-element manner. From [5], it is known that ${\mu }^{☆}={n\left|\right|\Delta \left|\right|}_{\max }$ marks the lower bound for $\mu$. From [5] and the short summary above, it follows that $\mathbf{R}={\mathbf{S}}^T{\mathbf{Z}}_{\mathrm{a}}\mathbf{S}$ holds. Furthermore, ${\mathbf{S}}^T\mathbf{S}=\mathbf{I}$ needs to be fulfilled to ensure the orthogonality of the transformation matrix. To include both conditions into the optimization problem of finding a suitable matrix $\mathbf{S}$, LMIs [15,21] are introduced. For this, the requirements are reformulated into the positive definite matrix inequalities

$$-\mathbf{R}+{\mathbf{S}}^T{\mathbf{Z}}_{\mathrm{a}}\mathbf{S}\succ \mathbf{0}\mathrm{and}\mathbf{I}-{\mathbf{S}}^T\mathbf{S}\succ \mathbf{0},$$

(32)

whose iterative solution in combination with the following cost function (36) leads—in the limit case—to the same results as the direct solution of the equality problem. Due to opposite signs of the quadratic terms in $\mathbf{S}$, the norm of $\mathbf{S}$ is bounded both from below and above, resulting in the chosen signs of the inequalities. In the next step, these quadratic matrix inequalities are converted into linear ones by applying the Schur complement formula according to

$$\left[\begin{array}{cc}-\mathbf{R}& {\mathbf{S}}^T\\ \mathbf{S}& -{\mathbf{Z}}_{\mathrm{a}}^{-1}\end{array}\right]\succ \mathbf{0}\mathrm{and}\left[\begin{array}{cc}\mathbf{I}& {\mathbf{S}}^T\\ \mathbf{S}& \mathbf{I}\end{array}\right]\succ \mathbf{0},$$

(33)

respectively. Obviously, $\mathbf{R}$ is again defined as

$$\mathbf{R}=\overline{\mu }{\mathbf{E}}_n-\mathsf{\Gamma },\overline{\mu }>\mu ,$$

(34)

where the LMI constraints

$$\mathsf{\Gamma }\succ \mathbf{0}\mathrm{and}{\mathbf{R}}^T\mathbf{Q}+\mathbf{Q}\mathbf{R}\prec \mathbf{0}$$

(35)

with $\mathbf{Q}\succ \mathbf{0}$, here chosen as $\mathbf{Q}=\mathbf{I}$, represents the fact that the midpoint matrix to be transformed is assumed to be asymptotically stable. This is guaranteed by the control designed before the transformation. To find a unique solution for the transformation matrix $\mathbf{S}$ fulfilling all requirements, the LMIs (33)–(35) are solved not only for $\mathbf{S}$, but for the diagonal matrix $\mathsf{\Gamma }$ (which is not restricted to identical entries for all diagonal elements), as well as for the scalar $\overline{\mu }$ together with a minimization of the cost function

$$J=\mathrm{tr}\left(\mathsf{\Gamma }\right)+\mathrm{tr}({\mathbf{Z}}_{\mathrm{a}}\mathbf{S}-\breve{\mathbf{S}}\mathbf{R})-\kappa ·\mathrm{tr}({\breve{\mathbf{S}}}^T\mathbf{S}-\mathbf{I})$$

(36)

with the problem-dependent parameter $\kappa >0$. As this parameter serves as a tool to prevent the trivial solution $\mathbf{S}=\mathbf{I}$, it has to be chosen large enough to fulfill this, but small enough not to dominate (36). Therefore, we start with a large value, and reduce it if no solution can be found. The optimization task, which includes only point-valued matrices, is solved in an iterative manner, where $\breve{\mathbf{S}}$ denotes the solution of the last successful evaluation of the LMI-constrained optimization task to render this cost function linear in $\mathbf{S}$. Algorithm 1 gives an overview of the solution procedure.

Algorithm 1: LMI-based computation of the transformation matrix $\mathsf{\Theta }=\mathbf{VS}$.

To enhance the numerical convergence, it starts with some $\mu <{\mu }^{☆}$, which is gradually increased with a line-search rule ${\mu }_{+}=\mu +\Delta \mu$ with $\Delta \mu >0$ until $\mu$ becomes equal to (or larger than) the desired value ${\mu }^{☆}$. This is done because the solution for $\mu =0$ corresponds to the known starting point $\mathbf{S}=\mathbf{I}$, specifying the initialization $\breve{\mathbf{S}}=\mathbf{I}$. Note that these randomly chosen initialization values of $\Delta \mu$ in the line-search may lead to different results. Therefore, the algorithm is implemented in such a way that the calculation is repeated 10 times. If the newly calculated intervals are tighter than the ones before, they are added to the list as a new optimum, if not, the solution is discarded. In the given scenario this leads to only two entries in

Table 2. Hulls of interval enclosures over the whole time horizon of $3\mathrm{ms}$.

Step	${\underline{\mathit{i}}}_{\mathbf{L}}$ in A	${\overline{\mathit{i}}}_{\mathbf{L}}$ in A	${\underline{\mathit{u}}}_{\mathbf{C}}$ in V	${\overline{\mathit{u}}}_{\mathbf{C}}$ in V
1	$-\underline{2.879}897$	$\phantom{-}\underline{2.884}897$	$-\underline{0.0017}43$	$\phantom{-}\underline{5.0017}43$
2	$-2.879542$	$\phantom{-}2.884542$	$-0.001731$	$\phantom{-}5.001731$

showing the calculated hulls of the interval enclosures where the underlined digits from the first execution of the algorithm correspond to the digits that are identical to those of the tightest solution.

As mentioned, the original system matrix does not have a diagonally dominant form, which leads to a combined transformation matrix of

$$\mathsf{\Theta }=\mathbf{VS}=\left[\begin{array}{cc}\phantom{-}0.7768& -0.5841\\ -0.9680& -0.2714\end{array}\right],$$

(37)

resulting in a system matrix

$${\widehat{\mathbf{A}}}_{\mathrm{C}}=\left[\begin{array}{cc}\left[-0.428;-0.260\right]& \left[0.218;0.286\right]\\ \left[0.540;0.763\right]& \left[-1.015;-0.925\right]\end{array}\right]\times {10}^4,$$

(38)

which is Metzler. The simulations are done for a time horizon of $3\mathrm{ms}$ with the initial states ${\mathbf{x}}^T\left(0\right)\in {\left[\begin{array}{cc}\left[0;5\times {10}^{-3}\right]& \left[0;5\right]\end{array}\right]}^T.$

5. Simulation Results

In this section, a simulation of the transformed system, the original model of which is given in Section 2, is analyzed, including the approach presented in the section before. This is done to give a more insightful view on the numerical results of Section 3 and Section 4. A state-of-the-art Taylor series expansion [4,22], without preconditioning the state equations by a QR or similar method, is used to compare the solutions of the presented transformation method. Here, a Taylor series expansion of order 2 is applied. A further increase of the order did not yield better results. Note that preconditioning is omitted to show the raw result of the method. Possibilities to reduce the interval width for enclosures calculated by a Taylor series expansion are given in [22]. Figure 5 shows the predicted interval enclosure of the first state $i_{\mathrm{L}}$.

Although the lower and upper bounds reach the desired operating point in less than $20\mathrm{ms}$, an overapproximation in the starting phase occurs due to the transformation—after all, it is transformed in two directions (transforming the system into a cooperative form and then transferring the computed state enclosures back into the original form to show the physical states)—when applying the method based on cooperativity. The Taylor series expansion also suffers from overestimation right from the beginning but is not able to reduce this afterwards. This is because the overapproximation of one step is mapped onto the next and so on, increasing the overestimation with each step. For a better comparison, Figure 6 shows the varying interval width for the cooperativity-based and the Taylor-series-based methods in logarithmic scale for the y-axis. As a result, an exponential growth for the Taylor series expansion method becomes visible, while the method based on cooperativity results in an exponential decrease of the interval width.

Here, it becomes clear that the failure of the Taylor series expansion is far more difficult to handle, as the interval widths widen towards infinite diameters in a very short time if no countermeasures such as preconditioning of the state equations are performed (Figure 6). Since this method is based on a discrete second-order Taylor approximation procedure, discretization errors are captured by further additive interval bounds [23], underlying the fact of the growing overestimation as mentioned before. If at all, this method is only possible for the offline prediction of state intervals and fails completely for systems with fast dynamics, as is the case in the presented real-life application. Nevertheless, a possible improvement for the starting phase could be achieved by intersecting both solution techniques. However, concerning the second state $u_{\mathrm{C}}$ in Figure 7, this is not necessary, as the cooperativity-based method does not suffer from higher overestimation of the initial interval, and hence already presents the better of both solutions directly from the beginning of the simulation.

In general, the Taylor series approach without preconditioning is not suited for the given system, especially in later phases. However, the cooperativity-based method shows that it is not only possible to calculate the verified states but that they also show the stable behavior of the controlled system. To underline the cooperativity and positivity of the system, Figure 8 and Figure 9 are included to show the respective states in their transformed coordinates.

6. Observer

In the following, an observer is added to the system. This represents a kind of fault diagnosis instrument to evaluate the measurements of the states that were so far assumed to be forthcoming for the controller. In the case that the true system and the observer model coincide, the interval observer outputs will enclose the true system state. If the real system violates the observed bounds, a sensor fault is detected. However, additional applications may be the estimation of the load resistance as well as the detection of component failures. Hence, the controller acts on the measurements alone and the observer is purely instated as a (separate) monitoring system diagnosing the closed-loop system dynamics. However, adding a state-of-the-art observer would destroy the cooperative structure, so finding a cooperativity-preserving observer is necessary so that

$$\left({\underline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\mathbf{HC}}_{\mathrm{m}}\right)·\widehat{\mathbf{v}}\left(t\right)+\mathit{\upsilon }=\dot{\widehat{\mathbf{v}}}\left(t\right)\le \dot{\widehat{\mathbf{x}}}\left(t\right)\le \dot{\widehat{\mathbf{w}}}\left(t\right)=\left({\overline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\mathbf{HC}}_{\mathrm{m}}\right)·\widehat{\mathbf{w}}\left(t\right)+\mathit{\omega }$$

(39)

holds with

$$\begin{array}{cc}\hfill {\underline{\widehat{\mathbf{A}}}}_{\mathrm{C}}& =inf\left({\widehat{\mathbf{A}}}_{\mathrm{C}}\left(\left[\mathbf{p}\right]\right)\right)\hfill \\ \hfill {\overline{\widehat{\mathbf{A}}}}_{\mathrm{C}}& =sup\left({\widehat{\mathbf{A}}}_{\mathrm{C}}\left(\left[\mathbf{p}\right]\right)\right)\hfill \end{array}$$

(40)

and

$$\begin{array}{cc}\hfill \mathit{\upsilon }& =inf\left(\mathbf{H}\left[{\mathbf{y}}_{\mathrm{m}}\right]\right)\hfill \\ \hfill \mathit{\omega }& =sup\left(\mathbf{H}\left[{\mathbf{y}}_{\mathrm{m}}\right]\right)\hfill \end{array}$$

(41)

for including possible worst-case bounds of measurements which lead to the uncertain measurement vector $\left[{\mathbf{y}}_{\mathrm{m}}\right]={\mathbf{y}}_{\mathrm{m}}+\left[-\Delta {\mathbf{y}}_{\mathrm{m}};\Delta {\mathbf{y}}_{\mathrm{m}}\right]$, cf. [10]. For the considered application scenario, a method given in [24] was used. Firstly, the system matrix

$${\mathbf{A}}_{\mathrm{O}}\left(\mathbf{p}\right)={\widehat{\mathbf{A}}}_{\mathrm{C}}\left(\mathbf{p}\right)-{\mathbf{HC}}_{\mathrm{m}}={\widehat{\mathbf{A}}}_{\mathrm{C}}\left(\mathbf{p}\right)-\mathcal{H}$$

(42)

has to become asymptotically stable, which can be investigated by

$${\mathbf{A}}_{\mathrm{O}}\left(\mathbf{p}\right)·{\mathbf{P}}_{\mathrm{O}}+{\mathbf{P}}_{\mathrm{O}}·{\mathbf{A}}_{\mathrm{O}}^T\left(\mathbf{p}\right)\prec \mathbf{0},$$

(43)

where ${\mathbf{P}}_{\mathrm{O}}$ again defines a Lyapunov function like before. If (39) holds, an error vector

$$\mathbf{e}={\left[\begin{array}{cc}{(\widehat{\mathbf{v}}-\mathbf{v})}^T& {(\widehat{\mathbf{w}}-\mathbf{w})}^T\end{array}\right]}^T$$

(44)

considering the difference between estimated and true values for both lower and upper state bounds can be obtained. The resulting ODEs

$$\dot{\mathbf{e}}=\left[\begin{array}{cc}{\underline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\mathbf{HC}}_{\mathrm{m}}& \mathbf{0}\\ \mathbf{0}& {\overline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\mathbf{HC}}_{\mathrm{m}}\end{array}\right]\mathbf{e}+\left[\begin{array}{c}\mathbf{H}\\ \mathbf{H}\end{array}\right]\mathit{\zeta }$$

(45)

denote the estimation errors in dependence of a measurement tolerance vector $\mathit{\zeta }$. The augmented system output

$${\mathbf{y}}_{\infty }=\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ -\nu ·\mathbf{I}& \nu ·\mathbf{I}\end{array}\right]\mathbf{e}+\left[\begin{array}{c}\mathbf{I}\\ \mathbf{0}\end{array}\right]\mathit{\zeta }={\mathbf{C}}_{\infty }\mathbf{e}+{\mathbf{D}}_{\infty }\mathit{\zeta }$$

(46)

accounts for a comparison of the measurement errors $\mathit{\zeta }$ and the weighted ($\nu >0$) state diameter $(\widehat{\mathbf{w}}-\mathbf{w})-(\widehat{\mathbf{v}}-\mathbf{v})$ [10,13]. A respective H${}_{\infty }$-likeLMI optimization problem can be formulated as

$$\mathcal{L}\left(\mathsf{\Theta }\right):=\left[\begin{array}{ccc}\mathsf{\Theta }& \mathbf{H}& {\mathbf{P}}_{\mathrm{O}}{\mathbf{C}}_{\infty }^T\\ {\mathbf{H}}^T& -\mathbf{I}& {\mathbf{D}}_{\infty ,1}^T\\ {\mathbf{C}}_{\infty }{\mathbf{P}}_{\mathrm{O}}& {\mathbf{D}}_{\infty ,1}& -{\gamma }_{\infty }^2\mathbf{I}\end{array}\right]\prec 0$$

(47)

for both extremal systems $\mathsf{\Theta }\in \{\underline{\mathsf{\Theta }},\overline{\mathsf{\Theta }}\}$ with

$$\begin{array}{cc}\hfill \underline{\mathsf{\Theta }}& :={\underline{\mathbf{A}}}_{\mathrm{O}}·{\mathbf{P}}_{\mathrm{O}}+{\mathbf{P}}_{\mathrm{O}}·{\underline{\mathbf{A}}}_{\mathrm{O}}^T\hfill \\ \hfill \overline{\mathsf{\Theta }}& :={\overline{\mathbf{A}}}_{\mathrm{O}}·{\mathbf{P}}_{\mathrm{O}}+{\mathbf{P}}_{\mathrm{O}}·{\overline{\mathbf{A}}}_{\mathrm{O}}^T\hfill \end{array}$$

(48)

using the abbreviations

$$\begin{array}{cc}\hfill {\underline{\mathbf{A}}}_{\mathrm{O}}& =inf\left({\mathbf{A}}_{\mathrm{O}}\left(\left[\mathbf{p}\right]\right)\right)\hfill \\ \hfill {\overline{\mathbf{A}}}_{\mathrm{O}}& =sup\left({\mathbf{A}}_{\mathrm{O}}\left(\left[\mathbf{p}\right]\right)\right).\hfill \end{array}$$

(49)

Note that the suitable Lyapunov function candidate is denoted by a joint—valid for both $\underline{\mathsf{\Theta }}$ and $\overline{\mathsf{\Theta }}$—weighting matrix ${\mathbf{P}}_{\mathrm{O}}={\mathbf{P}}_{\mathrm{O}}^T\succ 0$. With a linearizing change of variables

$${\mathbf{Q}}_{\mathrm{O}}={\mathbf{Q}}_{\mathrm{O}}^T={\mathbf{P}}_{\mathrm{O}}^{-1}\succ 0\mathrm{with}{\mathbf{Y}}_{\mathrm{O}}^T={\mathbf{Q}}_{\mathrm{O}}\mathbf{H}={\mathbf{P}}_{\mathrm{O}}^{-1}\mathbf{H},$$

(50)

the observer gain $\mathbf{H}$ could be calculated. However, to preserve the cooperative structure, another requirement must be added. An obvious choice to keep all off-diagonal elements of ${\mathbf{A}}_{\mathrm{O}}$ non-negative for all possible parameter combinations is given by

$$\mathbf{H}={\left({\mathcal{K}\mathit{C}}_{\mathrm{m}}\right)}^T\mathrm{with}\mathcal{K}=\mathrm{diag}\{\mathit{\kappa }\},$$

(51)

where

$$\mathit{\kappa }=\left[\begin{array}{ccc}{\kappa }_1& \cdots & {\kappa }_m\end{array}\right],\mathrm{where}{\kappa }_i>0,i\in \{1,\cdots ,m\}$$

(52)

hold. Moreover, the matrix ${\mathbf{C}}_{\mathrm{m}}$ includes exactly one entry equal to 1 per row while the rest equals zero. For ${\mathbf{A}}_{\mathrm{O}}\in {\mathbb{R}}^{n\times n}$ this results in

$$\mathcal{H}={\mathbf{HC}}_{\mathrm{m}}={\mathbf{C}}_{\mathrm{m}}^T\mathcal{K}{\mathbf{C}}_{\mathrm{m}},$$

(53)

where

$$\mathcal{H}=\left\{\begin{array}{cc}{\kappa }_i\hfill & \mathrm{for}i=j\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(54)

with $i\in \{1,\cdots ,m\}$ and $j\in \{1,\cdots ,n\}$ representing the simplest case. Including these considerations, Equation (47) is rearranged into

$$\mathcal{N}\left(\mathbf{\Xi }\right):=\left[\begin{array}{ccc}\mathbf{\Xi }& {\check{\mathbf{Q}}}_{\mathrm{O}}·{\left(\mathcal{K}\mathit{C}\right)}^T& {\mathbf{C}}_{\infty }^T\\ \left(\mathcal{K}\mathit{C}\right)·{\check{\mathbf{Q}}}_{\mathrm{O}}& -\mathbf{I}& {\mathbf{D}}_{\infty ,1}^T\\ {\mathbf{C}}_{\infty }& {\mathbf{D}}_{\infty ,1}& -{\mu }_{\infty }\mathbf{I}\end{array}\right]\prec 0$$

(55)

with $\mathbf{\Xi }\in \left\{\underline{\mathbf{\Xi }},\overline{\mathbf{\Xi }}\right\}$ according to

$$\begin{array}{cc}\hfill \underline{\mathbf{\Xi }}& :={\mathbf{Q}}_{\mathrm{O}}{\underline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\check{\mathbf{Q}}}_{\mathrm{O}}{\mathbf{C}}^T\mathcal{K}\mathit{C}+{\underline{\widehat{\mathbf{A}}}}_{\mathrm{C}}^T{\mathbf{Q}}_{\mathrm{O}}-{\mathbf{C}}^T\mathcal{K}\mathit{C}{\check{\mathbf{Q}}}_{\mathrm{O}}\hfill \\ \hfill \overline{\mathbf{\Xi }}& :={\mathbf{Q}}_{\mathrm{O}}{\overline{\widehat{\mathbf{A}}}}_{\mathrm{C}}-{\check{\mathbf{Q}}}_{\mathrm{O}}{\mathbf{C}}^T\mathcal{K}\mathit{C}+{\overline{\widehat{\mathbf{A}}}}_{\mathrm{C}}^T{\mathbf{Q}}_{\mathrm{O}}-{\mathbf{C}}^T\mathcal{K}\mathit{C}{\check{\mathbf{Q}}}_{\mathrm{O}}.\hfill \end{array}$$

(56)

To ensure that (56) is linear despite a multiplicative coupling of ${\mathbf{Q}}_{\mathrm{O}}$ and $\mathcal{K}$, the solution is determined iteratively. Additionally, setting ${\kappa }_1=\cdots ={\kappa }_m>0$ yields another simplification. Applying this approach to the given application scenario results in

$$\mathbf{h}{\tilde{\mathbf{c}}}^T=\left[\begin{array}{cc}0.0389& 0.1372\\ 0.1372& 0.4849\end{array}\right],$$

(57)

which is a valid solution despite not complying with Equation (51). This obviously gives an overall system matrix ${\tilde{\mathbf{A}}}_{\mathrm{O}}={\widehat{\mathbf{A}}}_{\mathrm{C}}-\mathbf{h}{\tilde{\mathbf{c}}}^T$ which is still Metzler. Note that the output vector ${\mathbf{c}}_{\mathrm{m}}^T$ of the assumed measured output $u_{\mathrm{C}}$ also needs to be transformed into the new coordinates by ${\tilde{\mathbf{c}}}^T={\mathbf{c}}_{\mathrm{m}}^T{\left(\mathbf{VS}\right)}^{-1}$ with ${\mathbf{c}}_{\mathrm{m}}^T=\left[\begin{array}{cc}0& 1\end{array}\right]$, as opposed to Equation (17). A possible application of this observer is shown in Figure 10, where two possible usages came to mind.

On the one hand, a sensor fault could be detected by comparing the output of the observer $\left[\widehat{y}\right]$ with the measured output of the system $y_{\mathrm{m}}$. On the other hand, a system fault diagnosis could be realized by checking whether ${\mathbf{x}}_{\mathrm{C}}\in \left[{\widehat{\mathbf{x}}}_{\mathrm{C}}\right]$.

In [24], other methods to find cooperativity-preserving observers were presented. However, for some problems it is not possible to design a cooperativity-preserving observer. For those, an observer can be implemented using an LMI-based approach like the presented controller computation of this paper. However, this would raise the need to transform the resulting system into a cooperative form afterwards. A transformation of an observer into cooperative form can be found, for example, in [5,25]. In the latter, the system matrix is given point-wise in contrast to our presented paper. Note further that the presented methods can also be expanded to fractional-order systems, as mentioned in [26].

7. Conclusions and Future Work

An electrical circuit with uncertain parameters was modeled and a controller design was presented based on LMIs considering a polytopic overapproximation of the uncertain system that includes all possible parameter combinations. The system was then transformed into a cooperative form to compute interval enclosures of the predicted states. This transformation is also LMI-based and was derived from findings in [5] rewritten into a computationally feasible problem according to [6]. Simulation results show the successful robust control design of the presented method, which was compared with a Taylor series expansion design. It was found that the cooperativity-based approach is more suitable for applications with fast dynamics and wide intervals for uncertain parameters such as the given electrical network. However, crossovers with the Taylor series expansion, especially for the starting phase, might help in reducing the still present overestimation. Another improvement could be provided by applying the findings of [27]. Here, the author transforms a subsystem model into a cooperative form while the remaining dynamics are stated in a non-cooperative way. The non-cooperative part needs to be evaluated by classical interval tools such as the suggested Taylor method and acts on the cooperative components as an additive bounded disturbance. If this disturbance is sufficiently small, it becomes possible to reduce the overestimation in comparison to scenarios where the complete system model is transformed into a new frame of coordinates. The final part of our paper presents a cooperativity-preserving observer, which can be used to evaluate measured data as a form of fault diagnosis. Although a simulation with the controlled system output showed the same results for the state as in the simulation scenarios, the observer needs to be evaluated with measured data from a test rig in the future.