Transformation of Uncertain Linear Systems with Real Eigenvalues into Cooperative Form: The Case of Constant and Time-Varying Bounded Parameters

Abstract

Continuous-time linear systems with uncertain parameters are widely used for modeling real-life processes. The uncertain parameters, contained in the system and input matrices, can be constant or time-varying. In the latter case, they may represent state dependencies of these matrices. Assuming bounded uncertainties, interval methods become applicable for a verified reachability analysis, for feasibility analysis of feedback controllers, or for the design of robust set-valued state estimators. The evaluation of these system models becomes computationally efficient after a transformation into a cooperative state-space representation, where the dynamics satisfy certain monotonicity properties with respect to the initial conditions. To obtain such representations, similarity transformations are required which are not trivial to find for sufficiently wide a-priori bounds of the uncertain parameters. This paper deals with the derivation and algorithmic comparison of two different transformation techniques for which their applicability to processes with constant and time-varying parameters has to be distinguished. An interval-based reachability analysis of the states of a simple electric step-down converter concludes this paper.

1. Introduction

In previous work, Raïssi et al. [1,2] and Mazenc et al. [3] derived various techniques for a transformation of uncertain linear systems into cooperative state-space representations. An autonomous set of ordinary differential equations (ODEs)

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

(1)

is cooperative [4,5,6] if the property

$$x_i^{\langle 1\rangle }\left(t\right)\ge x_i^{\langle 2\rangle }\left(t\right)\mathrm{for}\mathrm{all}i\in \{1,\dots ,n\}$$

(2)

holds for two vectors ${\mathbf{x}}^{\langle 1\rangle }\left(0\right)$ and ${\mathbf{x}}^{\langle 2\rangle }\left(0\right)$ of initial conditions which satisfy the inequalities

$$x_i^{\langle 1\rangle }\left(0\right)\ge x_i^{\langle 2\rangle }\left(0\right)\mathrm{for}\mathrm{all}i\in \{1,\dots ,n\}.$$

(3)

Cooperativity can be checked by the sufficient sign conditions

$$J_{i,j}\left(\mathbf{x}\right)\ge 0,i,j\in \{1,\dots ,n\},i\ne j$$

(4)

for all off-diagonal elements of the Jacobian

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

(5)

of the right-hand side of the state equations evaluated for all reachable states $\mathbf{x}=\mathbf{x}\left(t\right)$. Matrices satisfying this non-negativity property for the off-diagonal elements are also denoted as Metzler matrices in the literature [1,5,7,8].

For such systems, general-purpose, interval-based or other set-valued solvers [9] (AWA [10,11], CAPD [12], Cosy-VI [13,14], DynIBEX [15,16], ValEncIA-IVP [17,18,19,20,21], VNODE-LP [22,23], VSPODE [24]) for initial value problems (IVPs) can be replaced by point-valued simulations of a finite number of extremal realizations that can be extracted from the interval box

$$\mathbf{x}\left(0\right)\in \left[{\mathbf{x}}_0\right]=\left[\mathbf{x}\right]\left(0\right):=\left[\underset{‾}{\mathbf{x}}\left(0\right);\overline{\mathbf{x}}\left(0\right)\right]$$

(6)

of initial conditions. These point-valued extremal system realizations are denoted as bounding systems throughout this paper.

A special case of these cooperative systems are so-called positive systems [25]. For initial conditions 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\}\},$$

(7)

their trajectories will remain in the positive orthant for all $t\ge 0$ due to

$${\dot{x}}_i\left(t\right)=\mathbf{f}\left(x_1,\dots ,x_{i-1},0,x_{i+1},\dots x_n\right)\ge 0\mathrm{for}\mathrm{all}i\in \{1,\dots ,n\}.$$

(8)

Then, worst-case enclosures in the form $v_i\left(t\right)\le x_i\left(t\right)\le w_i\left(t\right)$ are obtained from the decoupled IVPs [26,27]

$$\dot{\mathbf{v}}\left(t\right)=\mathbf{f}\left(\mathbf{v}\left(t\right)\right)$$

(9)

and

$$\dot{\mathbf{w}}\left(t\right)=\mathbf{f}\left(\mathbf{w}\left(t\right)\right)$$

(10)

with the initial conditions

$$\mathbf{v}\left(0\right)=\underset{‾}{\mathbf{x}}\left(0\right)\mathrm{as}\mathrm{well}\mathrm{as}\mathbf{w}\left(0\right)=\overline{\mathbf{x}}\left(0\right).$$

(11)

Transformations of linear state equations which enforce cooperativity can be classified into either point-valued approaches or into techniques that employ an interval-valued change of coordinates [28]. In general, point-valued transformations are applicable for parameter-dependent linear systems with purely real eigenvalues. However, finding such transformations becomes more complex for an increasing degree of uncertainty in the system matrix. A failure of the similarity transformation approach by Raïssi et al. [1,2] can be recognized by the alternative formulation proposed by the authors in [29]. There, the transformation was cast into an optimization problem constrained by linear matrix inequalities (LMIs). Those LMIs may become infeasible (as a verification that the solution procedure is not successful) or the step size of the iteration in [29] reduces below a certain threshold leading to an excessively slow progress (a non-strict indicator that a point-valued similarity transformation may not exist). Even if interval-valued Metzler matrices can be found for the transformation of a provable asymptotically stable state equation, the transformed dynamics matrix may consist of an unstable upper bound. This commonly results from the wrapping effect in interval computations [30,31,32]. Using the aforementioned approaches, it is impossible to systematically avoid this phenomenon. However, it can be countered by the subdivision procedure introduced in this paper.

In contrast to point-valued transformations, also interval-valued approaches exist in the literature. They were originally developed by Mazenc et al. [3] for systems with complex eigenvalues and can be used as a fallback procedure if the point-valued transformation approach is not successful [28]. As shown in Section 2.4 of this paper, they have to be treated with care if uncertain parameters are time-varying. Time-varying parameters occur typically if nonlinear state equations are cast into sets of quasi-linear models. Then, the dynamics’ matrices themselves become functions of the system states [33]. A drawback of these interval-valued similarity transformations is an inevitable introduction of the wrapping effect [10,30]. As demonstrated in this paper, this is often counterproductive with respect to the tightness of the obtained state enclosures.

In Section 2, fundamentals of cooperative dynamic systems and existing similarity transformation procedures are summarized. Section 3 deals with the novel subdivision-based approach before a simulation-based comparison is given in Section 4. This comparison does not only investigate point- and interval-valued similarity transformations but also performs a comparison with the Taylor model-based verified ODE solver verifyode [34]. The application scenario considered in Section 4 is an electric step-down converter with a predefined duty cycle. Finally, conclusions and an outlook on future work are given in Section 5.

2. Cooperativity-Enforcing Similarity Transformations

2.1. Special Case: Linear and Quasi-Linear Systems with Bounded Parameters

As a special case of the general system model introduced in (1), assume the parameter-dependent linear ODEs

$$\dot{\mathbf{x}}\left(t\right)=\mathbf{A}\left(\mathbf{p}\right)·\mathbf{x}\left(t\right).$$

(12)

Here, the parameter vector $\mathbf{p}$ can either be uncertain but constant or time-varying within its interval bounds. For a scenario where the latter case was accounted for during a robust control synthesis, the reader is referred to [33].

For a compact notation of the system model (12), which contains all possible parameter-dependent realizations of $\mathbf{A}\left(\mathbf{p}\right)$, we introduce an interval matrix $\left[\mathbf{A}\right]$ that is formed of the bounds $\left[A_{i,j}\right]$ for each entry $i,j\in \{1,\dots ,n\}$. Those bounds are obtained according to the generally not minimal interval extension [30,31]

$$\left\{\begin{array}{c}\mathbf{p}\in \left[\mathbf{p}\right]\hfill \\ {\mathcal{A}}_{i,j}=A_{i,j}\left(\mathbf{p}\right)\hfill \end{array}\right.⟹{\mathcal{A}}_{i,j}\in \left[A_{i,j}\right]=\left[{\underset{‾}{A}}_{i,j};{\overline{A}}_{i,j}\right]=\left[A_{i,j}\right]\left(\left[\mathbf{p}\right]\right).$$

(13)

In this paper, it is foremost desired to find a time-invariant change of coordinates for the case that the inequalities

$${\underset{‾}{A}}_{i,j}\ge 0$$

(14)

are not satisfied for at least one $i,j\in \{1,\dots ,n\}$, $i\ne j$, i.e., that the uncertain linear system model is not proven to be cooperative. The change of coordinates should then lead to a set of state equations

$$\dot{\mathbf{z}}\left(t\right)=\widehat{\mathbf{A}}·\mathbf{z}\left(t\right)\mathrm{with}\mathbf{z}\left(t\right):={\Theta }^{-1}·\mathbf{x}\left(t\right)$$

(15)

with an interval-valued system matrix $\widehat{\mathbf{A}}\in [\widehat{\mathbf{A}}]$ with non-negative off-diagonal elements ${\widehat{\underset{‾}{A}}}_{i,j}\ge 0$ for each $i\ne j\in \{1,\dots ,n\}$. Unfortunately, this change of coordinates may lead to the case that—despite the original system model (12) was asymptotically stable for all possible $\mathbf{p}\in \left[\mathbf{p}\right]$—the transformed matrix $[\widehat{\mathbf{A}}]$ contains unstable realizations due to the wrapping effect.

Especially in the area of control engineering, dynamic system models are often not purely autonomous as stated in (12). They typically contain additive input terms $\mathbf{B}·\mathbf{u}\left(t\right)$. Due to the fact that these terms can be included additively in Equation (12), they can also be included additively in the transformed model (15) by the term

$${\Theta }^{-1}·\mathbf{B}·\mathbf{u}\left(t\right),$$

(16)

where uncertainty needs to be taken into account by the same $inf\left(·\right)$ and $sup\left(·\right)$ operators that are included subsequently in Equation (17).

As described in the introduction, cooperativity (together with stability) allows for simplifying the simulation of the transformed system model (15). Guaranteed bounds $\mathbf{z}\left(t\right)\in \left[\mathbf{z}\right]\left(t\right)=\left[\mathbf{v}\left(t\right);\mathbf{w}\left(t\right)\right]$ of all reachable states in the new coordinate frame $\mathbf{z}$ are obtained by an evaluation of the following coupled set of state equations (which are a direct consequence of [35])

$$\begin{array}{cc}\hfill {\dot{v}}_i\left(t\right)& =inf\left([{\widehat{A}}_{i,i}]·v_i\left(t\right)\right)+inf\left(\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n[{\widehat{A}}_{i,j}]·\left[z_j\right]\left(t\right)\right)\hfill \\ \hfill {\dot{w}}_i\left(t\right)& =sup\left([{\widehat{A}}_{i,i}]·w_i\left(t\right)\right)+sup\left(\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n[{\widehat{A}}_{i,j}]·\left[z_j\right]\left(t\right)\right)\hfill \end{array}$$

(17)

with the resulting state bounds

$$\left[z_i\right]\left(t\right)=\left[v_i\left(t\right);w_i\left(t\right)\right],$$

(18)

where couplings between the vectors $\mathbf{v}\left(t\right)$ and $\mathbf{w}\left(t\right)$ can be ignored if the system is positive, i.e., $v_i\left(t\right)\ge 0$ holds for all $t\ge 0$ as discussed in (9)–(11).

Then, the bounding systems given in (17) simplify to

$${\dot{v}}_i\left(t\right)=\sum _{j=1}^{n}inf\left([{\widehat{A}}_{i,j}]\right)·v_j\left(t\right)\mathrm{and}{\dot{w}}_i\left(t\right)=\sum _{j=1}^{n}sup\left([{\widehat{A}}_{i,j}]\right)·w_j\left(t\right).$$

(19)

The expressions for the bounding systems (17) and (19) are also a direct consequence of the results given in [8] (Lemma 1) and [36] (Equation (15)).

2.2. Illustrating Example

As an illustrating example for the evaluation of uncertain system models according to Equations (17), consider a set of ODEs with the uncertain parameters and initial conditions

$$\widehat{\mathbf{A}}\in \left[\begin{array}{cc}\left[-2;1\right]& \left[0.5;1\right]\\ \left[0.5;1\right]& \left[-4;-2\right]\end{array}\right]\mathrm{and}\mathbf{z}\left(0\right)\in \left[\begin{array}{c}\left[-1;2\right]\\ \left[-1;2\right]\end{array}\right].$$

(20)

An equidistant gridding of all six independently uncertain intervals for the initial conditions and matrix entries into eight points for the initial states and four points for each entry in the matrix $\widehat{\mathbf{A}}$ (resulting in 16,384 simulations of systems of order $n=2$) yields the gray trajectories in Figure 1, while the black enclosures are obtained by a single simulation of the bounding system (17) of order $2n=4$.

The following two subsections make a distinction between changes of coordinates that are applicable if the original system (12) has only purely real eigenvalues or if it may also contain conjugate-complex ones. The respective two methods are only briefly presented in the following because they are a result of the previous publications [2,28,29].

2.3. Purely Real Eigenvalues

In this case, the transformation is based on the work of Raïssi et al. in [2]. There, a general method to transform uncertain systems with purely real eigenvalues into a cooperative form (15) was presented. To extend its applicability towards larger uncertainty bounds, an extension was derived in [29], where the original approach was redesigned into an LMI-constrained optimization task [37] to make it more generally applicable and to simplify its practical use.

We assume that the uncertain system matrix $[\mathbf{A}]$ can be enclosed by the element-wise defined inequality

$${\mathbf{Z}}_{\mathrm{a}}-\Delta \le \mathbf{Z}:={\mathbf{V}}^{-1}·\mathbf{A}·\mathbf{V}\le {\mathbf{Z}}_{\mathrm{a}}+\Delta$$

(21)

for all $\mathbf{A}\in [\mathbf{A}]$ defined in (13), where ${\mathbf{Z}}_{\mathrm{a}}={\mathbf{Z}}_{\mathrm{a}}^T$ in (21) is a point matrix. Typical choices for ${\mathbf{Z}}_{\mathrm{a}}$ are the midpoint of $[\mathbf{A}]$ if this matrix is diagonally dominant (with $\mathbf{V}=\mathbf{I}$) or the midpoint of a matrix obtained after a principal axis transformation using the invertible real-valued eigenvector matrix $\mathbf{V}$.

Next, a Metzler matrix $\mathbf{R}=\mu {\mathbf{E}}_n-\Gamma$ is searched for, which has the same eigenvalues as ${\mathbf{Z}}_{\mathrm{a}}$. For that, $\mu \in \mathbb{R}$ is a constant, ${\mathbf{E}}_n\in {\mathbb{R}}^{n\times n}$ is a matrix with all elements equal to 1, and $\Gamma \in {\mathbb{R}}^{n\times n}$ is a diagonal matrix. Following the procedure detailed in [2], there exists an orthogonal matrix $\mathbf{S}\in {\mathbb{R}}^{n\times n}$ such that ${\mathbf{S}}^T\mathbf{ZS}$ is Metzler with $\mathrm{eig}\left(\mathbf{R}\right)=\mathrm{eig}\left({\mathbf{Z}}_{\mathrm{a}}\right)$ and ${\mu >n\left|\right|\Delta \left|\right|}_{\max }$, where ${\left|\right|\Delta \left|\right|}_{\max }$ is the maximum absolute value of $\Delta$. The overall transformation according to (15) becomes $\Theta :=\mathbf{V}·\mathbf{S}$, which is point-valued and time-invariant. The LMI-based optimization in [29] automatizes the search for the matrix $\mathbf{S}$. It terminates successfully, if all realizations of the interval matrix $[\widehat{\mathbf{A}}]:={\Theta }^{-1}·\left[\mathbf{A}\right]·\Theta$ are Metzler.

Remark 1.

The LMI-based optimization may fail to find a solution because of two reasons. Either the LMI solver (e.g., SeDuMiin combination with Yalmip [38,39]) may produce the output of being infeasible. Then, it is ensured that the problem formulation does not have a suitable solution. However, previous work has shown that this is typically not the case. Mostly, the step size control procedure, which increases the parameter μ defined above up to the point where ${\mu >n\left|\right|\Delta \left|\right|}_{\max }$ is satisfied, progresses too slowly to find a solution in acceptable time. This is an indicator that the intervals are too wide to find a common point-valued transformation for the uncertain system model.

Instabilities and the failure of the LMI-based solution are countered by the subdivision in Section 3.

2.4. Mixed and Conjugate-Complex Eigenvalues

If the original system (12) contains conjugate-complex eigenvalues or if a solution with the help of the previous approach cannot be found for the complete uncertainty domain, the following procedure can be applied. The most important differences in comparison with the previous approach are:

• the transformation is usually time-varying [1,3];
• mapping the uncertainty directly into the location of the eigenvalues turns the transformation into an interval-valued expression.

Hence, evaluating the system matrix $\mathbf{A}\left(\mathbf{p}\right)$ in (12) for the whole range of parameters leads to a variability of the real and imaginary parts ${\sigma }_i$ and ${\omega }_i$ of the eigenvalues, where

$$\left[{\sigma }_i\right]=\left[{\underline{\sigma }}_i;{\overline{\sigma }}_i\right]\mathrm{and}\left[{\omega }_i\right]=\left[{\underline{\omega }}_i;{\overline{\omega }}_i\right]$$

(22)

are their respective interval bounds. Typically, it is necessary to determine these bounds with the help of interval techniques. Possible solution procedures (verifyeig) are included in the Matlab toolbox IntLab [40].

Assume that the system under investigation contains $\tilde{n}$ mutually disjoint conjugate-complex eigenvalue pairs with $n^{*}=n-2\tilde{n}\ge 0$ as the number of remaining real eigenvalues. To simplify the notation, assume further that all eigenvalues are sorted so that all complex pairs are listed first, cf. [28].

Now, an equivalent system representation in block diagonal structure with

$$\tilde{\mathbf{A}}=\mathrm{blkdiag}\left({\tilde{\mathbf{A}}}_1,\dots ,{\tilde{\mathbf{A}}}_{\tilde{n}+n^{*}}\right)$$

(23)

can be obtained, where pairs of conjugate-complex eigenvalues with

$${\omega }_i=-{\omega }_{i+1},i\in \{1,3,\dots ,2\tilde{n}-1\},$$

(24)

lead to the blocks

$${\tilde{\mathbf{A}}}_j\in \left[\begin{array}{cc}[{\sigma }_j]& [{\omega }_j]\\ -[{\omega }_j]& [{\sigma }_j]\end{array}\right]$$

(25)

with $j\in \{1,2,\dots ,\tilde{n}\}$, while an uncertain real eigenvalue is reflected by

$${\tilde{\mathbf{A}}}_j\in [{\sigma }_i],j=i-\tilde{n},i\in \{2\tilde{n}+1,2\tilde{n}+2,\dots ,n\}.$$

(26)

The respective transformation between the original system matrix $\mathbf{A}\left(\mathbf{p}\right)$ and the representation in (23) is given by the column-wise partitioned matrix

$$\tilde{\mathbf{T}}\in \left[\tilde{\mathbf{T}}\right]=\left[\left[{\tilde{\mathbf{T}}}_1\right],\dots ,\left[{\tilde{\mathbf{T}}}_{\tilde{n}+n^{*}}\right]\right]$$

(27)

with

$$\left[{\tilde{\mathbf{T}}}_j\right]=\left[\begin{array}{cc}\Re \left\{\left[{\mathbf{v}}_j\right]\right\}& \Im \left\{\left[{\mathbf{v}}_j\right]\right\}\end{array}\right],j\in \{1,\dots ,\tilde{n}\},$$

(28)

consisting of interval enclosures for the real and imaginary parts of the eigenvectors of the interval matrix $\left[\mathbf{A}\right]$ and the eigenvectors

$$\left[{\tilde{\mathbf{T}}}_j\right]=\left[{\mathbf{v}}_j\right],j=i-\tilde{n},$$

(29)

for each of the real eigenvalues with $i\in \{2\tilde{n}+1,2\tilde{n}+2,\dots ,n\}$.

In a second stage, a (generally) time-varying transformation is performed according to

$$\mathbf{z}={\mathbf{S}}^{-1}\left(t\right)·\tilde{\mathbf{z}},$$

(30)

where

$$\mathbf{S}\left(t\right)={\left({\mathbf{S}}^{-1}\left(t\right)\right)}^T=\mathrm{blkdiag}\left({\mathbf{S}}_1\left(t\right),\dots ,{\mathbf{S}}_{\tilde{n}+n^{*}}\left(t\right)\right)$$

(31)

consists of the orthogonal blocks

$${\mathbf{S}}_j\left(t\right)=\left[\begin{array}{cc}cos\left({\omega }_{j}t\right)& sin\left({\omega }_{j}t\right)\\ -sin\left({\omega }_{j}t\right)& cos\left({\omega }_{j}t\right)\end{array}\right],j\in \{1,2,\dots ,\tilde{n}\},$$

(32)

and

$${\mathbf{S}}_{\tilde{n}+1}=\dots ={\mathbf{S}}_{\tilde{n}+n^{*}}=1.$$

(33)

A symbolic simplification of the transformed state-space representation according to

$$\begin{array}{cc}\hfill \dot{\mathbf{z}}& ={\dot{\mathbf{S}}}^T\left(t\right)·\tilde{\mathbf{z}}+{\mathbf{S}}^T\left(t\right)·\dot{\tilde{\mathbf{z}}}\hfill \\ \hfill & =\left[\left[{\displaystyle \frac{\mathrm{d}{\mathbf{S}}^T\left(t\right)}{\mathrm{d}t}}+{\mathbf{S}}^T\left(t\right)\tilde{\mathbf{A}}\right]\mathbf{S}\left(t\right)\right]\mathbf{z}\hfill \\ \hfill & =\mathbf{N}·\mathbf{z}\hfill \end{array}$$

(34)

yields a purely diagonal system matrix

$$\mathbf{N}=\mathrm{blkdiag}\left({\sigma }_1·\mathbf{I},\dots ,{\sigma }_{\tilde{n}}·\mathbf{I},{\sigma }_{\tilde{n}+1},\dots ,{\sigma }_{\tilde{n}+n^{*}}\right)\mathrm{with}\mathbf{I}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

(35)

in the new coordinates. For simulation purposes, the system model (34) with (35) is evaluated for the interval matrix

$$\left[\mathbf{N}\right]=\mathrm{blkdiag}\left(\left[{\sigma }_1\right]·\mathbf{I},\dots ,\left[{\sigma }_{\tilde{n}}\right]·\mathbf{I},\left[{\sigma }_{\tilde{n}+1}\right],\dots ,\left[{\sigma }_{\tilde{n}+n^{*}}\right]\right),$$

(36)

where the change of coordinates (15) is based on the overall interval-valued transformation

$$\left[\Theta \right]\left(t\right)=[\tilde{\mathbf{T}}]·[\mathbf{S}]\left(t\right).$$

(37)

As shown before, the advantage of this transformation procedure is that it is equally applicable to systems with real and complex eigenvalues. Moreover, for the purely real case, it breaks again down to a time-invariant transformation because all matrices in (32) then become equal to the identity matrix. However, although the system dynamics according to (36) are fully decoupled in the new frame of coordinates, overestimation occurs due to the fact that an interval-valued transformation needs to be employed even if the initial conditions $\mathbf{x}\left(0\right)$ according to (6) were given by point-valued data.

Remark 2.

Because the approaches in Section 2.3 and Section 2.4 lead to time-invariant coordinate transformations with a single resulting set of bounding systems and because the time derivative in the square brackets of (34) cancels out for purely real eigenvalues, both procedures are equally applicable if the matrix $\mathbf{A}\left(\mathbf{p}\right)$ is influenced by either constant or time-varying bounded parameters $\mathbf{p}\in \left[\mathbf{p}\right]$. For complex eigenvalues, the transformation in Section 2.4 and its simplification according to (34) are only valid if all ${\omega }_j$ are constant for known time intervals on which the simulation is carried out. For fast, arbitrary variation rates of $\mathbf{p}$, implying fast changes of ${\omega }_j\left(t\right)$, this equation needs to be treated with special care.

3. Novel Subdivision-Based State-Space Transformation

As the bridging element between both transformations summarized in the previous Section 2.3 and Section 2.4, a novel subdivision-based procedure is introduced. Its basic components are summarized in the following Algorithms 1–4. In a first stage, the procedure given in Section 2.3 for systems with real eigenvalues is employed (line 1 in Algorithm 1). The novel components of this algorithm are activated in the case that—after a predefined wall-clock time—no suitable transformation of the uncertain system model into an interval-valued Metzler matrix has been found, or if the resulting transformed model contains unstable realizations despite provable asymptotic stability of the original state Equation (12).

To reduce the wrapping effect during the transformation of the state equations, this novel algorithm does not only try to find transformations for the complete parameter domain ${\mathcal{L}}^{\langle *\rangle }.\tilde{\mathbf{p}}$ (which denotes a storage element of the currently investigated interval box for the system parameter) as a whole. In addition, the list elements of yet undecided parameter subintervals (i.e., those intervals for which a solution of the desired transformation has not yet been found) also carry several ${\mathcal{L}}^{\langle *\rangle }.N$ subintervals, which are all stored in ${\mathcal{L}}^{\langle *\rangle }.{\mathbf{p}}_1$. In such a way, it becomes possible to individually investigate each of the corresponding subboxes, whether at least for some of them the algorithm of Section 2.3 produces stable realizations of the transformed system model.

According to Remark 1, a stopping criterion needs to be implemented which avoids an endless search for point-valued transformations. Currently, this is done by checking whether a timeout occurred (line 1 in Algorithm 1). Suitable alternatives would be progression rates for the parameter $\mu$ which fall below a certain threshold or a maximum number of non-successful trials. If this stopping criterion becomes active, the stable intervals—transformed into interval-valued Metzler matrices—are appended to the list ${\mathcal{L}}_{\mathrm{s}}$ (starting in line 1 of Algorithm 1, as well as Algorithm 2).

Algorithm 1: Splitting-based transformation

Algorithm 2: Stability analysis

For all undecided intervals (either without existing cooperativity-enforcing transformation or lacking the property of asymptotic stability), two different subdivision procedures are used in lines 1, 1, and 1 of Algorithm 1. Details about their distinction are summarized in Algorithms 3 and 4. Basically, they either break up the collection of $\mathcal{N}.N$ subintervals $\mathcal{N}.\tilde{\mathbf{p}}$ into $\mathcal{N}.N$ yet undecided individual boxes in Algorithm 3, or split the largest parameter box into smaller subintervals after appending the $\mathcal{N}.N$ individual entries to the undecided list in Algorithm 4.

Algorithm 3: Subdivision procedure A

Algorithm 4: Subdivision procedure B

In such a way, the list ${\mathcal{L}}_{\mathrm{s}}$, obtained after completion of the while-loop in Algorithm 1, contains individual state-space transformations into asymptotically stable subsystem models. The resulting transformed subsystems are always valid for uncertain but constant parameters. Temporal variations, however, are now no longer admissible within the complete initial parameter box. Instead, due to the subdivision into a collection of mutually independent subsystem models (corresponding to the length of ${\mathcal{L}}_{\mathrm{s}}$), temporal variations of parameters are only allowed within each individual box ${\mathcal{L}}_{\mathrm{s}}.\tilde{\mathbf{p}}$. Finally, in cases in which it is not a-priori known that the original system model has purely real eigenvalues and that the while-loop in line 1 of Algorithm 1—based on the procedure from Section 2.3—will, hence, terminate after a sufficiently large number of subdivisions, a fallback to the interval-based transformation of Section 2.4 should be included.

4. Application Scenario: Step-Down Converter

The algorithms for a transformation of dynamic system models into a cooperative state-space representation are now compared by means of numerical simulations for the state equations of a linear step-down converter circuit. All interval evaluations are performed in IntLab V. 12 [40]. For the numerical simulation of point-valued bounding systems, the floating-point solver ode23 with standard tolerance setting and the maximum step size ${10}^{-5}$ was used. This solver is sufficiently accurate and more effective than, for example, ode45 for the moderately stiff system model under consideration. For this application, the resulting approximation errors are also several orders of magnitude smaller than the absolute values of the corresponding states. This was checked by a computation of the corresponding solutions with the help of symbolic formula manipulation in Matlab. In addition to the presentation of the transformation approaches introduced in this paper, a comparison with the verified ODE solver verifyode is illustrated. This solver is also included in IntLab.

4.1. Modeling

According to Figure 2, a step-down converter is combined with a fuse for excessive current protection, a variable load $R_{\mathrm{S}}={\tilde{R}}_{\mathrm{S}}+\Delta R_{\mathrm{S}}$, and the variably connectable additional resistor ${\tilde{R}}_{\mathrm{C}}$. For simplicity of the following state equations, the diode forward bias is assumed to be zero. In future work, it can be included in a straightforward manner by means of an additional system parameter in the corresponding voltage loop equation.

Here, the equations (denoted for brevity without explicitly mentioning time arguments)

$$u_{\mathrm{E}}=u_{R_{\mathrm{L}}}+u_{\mathrm{L}}+u_{R_{\mathrm{C}}}+u_{\mathrm{C}}\mathrm{and}{u}_{\mathrm{C}}+u_{R_{\mathrm{C}}}=u_{R_{\mathrm{S}}},$$

(38)

with

$$u_{R_{\mathrm{S}}}=u_{{\tilde{R}}_{\mathrm{S}}}+\Delta u_{R_{\mathrm{S}}},$$

(39)

describe the two voltage loops, while

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

(40)

results from Kirchhoff’s node equation. The component equations of all Ohmic resistances are given by

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

(41)

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

$$u_{\mathrm{L}}=L·{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_{\mathrm{L}}\mathrm{and}{i}_{\mathrm{C}}=C·{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{u}_{\mathrm{C}}.$$

(42)

Here, the physical energy storage is expressed by the current $i_{\mathrm{L}}$ and the voltage $u_{\mathrm{C}}$, which results in the ODE

$$\begin{array}{cc}\hfill {\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}}}{R_{\mathrm{S}}+R_{\mathrm{C}}}}\right)·u_{\mathrm{C}}+u_{\mathrm{E}}\right)\hfill \end{array}$$

(43)

to describe variations of the magnetic field energy and in

$${\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}}$$

(44)

for changes of the electric field energy. The corresponding state-space representation with

$$\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}{i}_{\mathrm{L}}\\ u_{\mathrm{C}}\end{array}\right]$$

(45)

results in

$$\begin{array}{cc}\hfill \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}}\hfill \\ \hfill & =\mathbf{A}·\mathbf{x}+\mathbf{b}·u_{\mathrm{E}}.\hfill \end{array}$$

(46)

For the simulations in the following subsection with $\mathbf{x}\left(0\right)=\mathbf{0}$, the parameters in Equation (46) are set to $L=1\mathrm{H}$, ${\check{R}}_{\mathrm{L}}=100\Omega$ (together mimicking a real inductor with non-zero internal resistance; which can be implemented for a test rig with the help of an operational amplifier-based gyrator circuit), $C=2\mathrm{mF}$, $R_{\mathrm{C}}\in \left[0.1;0.6\right]\Omega$, and $R_{\mathrm{S}}\in \left[0.1;3\right]\Omega$.

In addition, the input voltage of the system is defined as a piecewise constant, periodic signal of duration $T=5\mathrm{ms}$ according to

$$u_{\mathrm{E}}=\left\{\begin{array}{cc}5\mathrm{V}\hfill & \mathrm{for}\left(i-1\right)·T\le t<\left(i-1\right)·T+\eta ·T\hfill \\ 0\mathrm{V}\hfill & \mathrm{for}\left(i-1\right)·T+\eta ·T\le t\le i·T,i\in \mathbb{N},\hfill \end{array}\right.$$

(47)

where the constant parameter $\eta =0.6$ specifies the length of the duty cycle.

4.2. Simulation-Based Comparison of Two Cooperativity-Enforcing Similarity Transformations

Considering the parameters listed in the previous subsection, it can be shown that there does not exist a common, point-valued transformation matrix $\Theta$ according to Section 2.3 that simultaneously leads to a Metzler matrix representation for the transformed system dynamics stated in Equation (15) and at the same time preserves the asymptotic stability of the original system (46).

Hence, a transformation of the system model into an interval-valued diagonal structure according to (34) and (36) has been performed first. After a subdivision of the parameter box

$$\left[\mathbf{p}\right]=\left[\begin{array}{cc}\left[R_{\mathrm{C}}\right]& \left[R_{\mathrm{S}}\right]\end{array}\right]$$

(48)

in a regular grid with $500\times 200$ subintervals, the following eigenvalue bounds were obtained by using IntLab’s verifyeig routine:

$$\left[{\sigma }_1\right]=\left[-115.15;-100.1\right]\mathrm{and}\left[{\sigma }_2\right]=\left[-2500;-124.2\right].$$

(49)

Due to the fact that these eigenvalues are purely real, the transformation matrix from (37) becomes time-invariant according to

$$\Theta \in \left[\begin{array}{cc}\left[0.0568;0.9948\right]& \left[0.0002;0.0352\right]\\ \left[0.1020;0.9984\right]& \left[0.9993;1.0000\right]\end{array}\right].$$

(50)

Here, the entries are displayed after outward rounding to four decimals. The simulation routine makes use of $\Theta$ to convert the results $\left[\mathbf{z}\right]\left(t\right)$ back into the original state variables $\mathbf{x}\left(t\right)$. The inverse

$${\Theta }^{-1}\in \left[\begin{array}{cc}\left[1.0052;45.9135\right]& \left[-1.6126;-0.0002\right]\\ \left[-45.8675;-0.1025\right]& \left[1.0000;2.6116\right]\end{array}\right]$$

(51)

of this matrix is required to express the influence of the control action according to (16). The interval enclosures for the temporal evolution of the state variables are depicted in Figure 3. These bounds are valid both, see Remark 2, for time-invariant uncertain parameters and for parameters that are changing their values arbitrarily within the complete uncertainty intervals.

In contrast, the novel solution procedure, with $N=10$ and $N=50$ in Algorithm 1 yields significantly tighter state enclosures as shown in Figure 4. Here, the solution is represented by the interval hull over all individual state enclosures, restricting the validity of the computed bounds to variations of the parameters within each of the boxes displayed in Figure 5. Please note that increasing the parameter N in the multi-sectioning strategy does not necessarily lead to tighter interval bounds. Hence, its systematic optimization in combination with a sensitivity-based selection of the splitting direction [41,42] will be especially helpful to make the procedure applicable to systems with a larger number of uncertain parameters in future work.

4.3. Comparison with a Taylor Model-Based Solution Approach

In general, there exist two different options to account for the uncertain parameter vector $\mathbf{p}\in \left[\mathbf{p}\right]$ defined in Equation (48) when using a general-purpose verified ODE solver. In this section, the solver verifyode [34,43] is employed to obtain a representative comparison.

When appending the uncertain parameters to the state vector $\mathbf{x}\left(t\right)$, introduced in (45), and accounting for the time invariance of these parameters by the so-called integrator disturbance model $\dot{\mathbf{p}}=\mathbf{0}$ in the state equations, the solver verifyode does not at all succeed in computing guaranteed state bounds. This results from a division by zero when substituting the Taylor model representations of both uncertain parameters into the system model (46). This problem is independent whether identical Taylor model orders are chosen for all four state variables or if the Taylor model order for the interval parameters is reduced to its minimum value 1 as described in [43].

A simulation of this augmented state-space representation becomes possible—even with the default setting of identical Taylor model orders for all state variables—if the uncertainty in both parameters $\mathbf{p}$ is reduced significantly. This, however, changes the system dynamics and can therefore not be used for a fair comparison with the proposed cooperativity-enforcing simulation approaches.

Successful simulations with this ODE solver become possible if the uncertain parameters are not appended to the state vector but if they are specified instead within the state equations as interval variables (datatype intval). Then, manually specifying the initial step size of the solver as $h_0={10}^{-4}$, the minimum step size as $h_{min}={10}^{-6}$, the Taylor model orders as either 10, 12 (the default setting), or 30, performing either a QR preconditioning or a curvilinear preconditioning with or without blunting (a strategy of widening the angles between the column vectors of an ill-conditioned preconditioning matrix to reduce overestimation [44]) and/or shrink wrapping (aiming at an elimination of additive error intervals and incorporating them in the Taylor model coefficients [45]), it is possible to simulate the system model.

For the comparison with the proposed solution approach described in Section 3, the outcome of verifyode is investigated for the following solver settings:

case 1

with default settings (order 12) for all options [43] except for the step size parameters and the choice of a QR preconditioning

verifyodeset(’h0’, 1e-4, ’h_min’, 1e-6, ’precondition’, 1);

case 2

with order 10, manually specified small tolerances, and QR preconditioning

verifyodeset(’order’, 10, ’shrinkwrap’, 0, ’precondition’, 1,  ’blunting’, 0, ...

 ’h0’, 1e-4, ’h_min’, 1e-6, ’loc_err_tol’, 1e-11, ’sparsity_tol’, 1e-20);

case 3

with order 30, manually specified small tolerances, and QR preconditioning

verifyodeset(’order’, 30, ’shrinkwrap’, 0, ’precondition’, 1, ’blunting’, 0, ...

 ’h0’, 1e-4, ’h_min’, 1e-6, ’loc_err_tol’, 1e-11, ’sparsity_tol’, 1e-20);

case 4

with order 10, manually specified small tolerances, and curvilinear preconditioning

verifyodeset(’order’, 10, ’shrinkwrap’, 0, ’precondition’, 3, ’blunting’, 0, ...

 ’h0’, 1e-4, ’h_min’, 1e-6, ’loc_err_tol’, 1e-11, ’sparsity_tol’, 1e-20);

case 5

with order 30, manually specified small tolerances, and curvilinear preconditioning

verifyodeset(’order’, 30, ’shrinkwrap’, 0, ’precondition’, 3, ’blunting’, 0, ...

 ’h0’, 1e-4, ’h_min’, 1e-6, ’loc_err_tol’, 1e-11, ’sparsity_tol’, 1e-20);

The simulations have shown that the cases 1–3 complete successfully up to the point $t=75\mathrm{ms}$ (the same final time instant as in Figure 3 and Figure 4, comprising 15 full duty cycles), while the state enclosures in the cases 4–5, in which the QR preconditioning was replaced by the curvilinear one blow up and break down shortly after $t=25\mathrm{ms}$ because the integration step size falls below the specified threshold. Setting the step size parameters to the default values as described in [43] does not resolve this issue. Activating the options for blunting and/or shrink wrapping did not influence the solutions for the considered application. Hence, they are not investigated further in this paper. Therefore, the following graphical comparison in Figure 6, Figure 7 and Figure 8 will be restricted to the cases 1–3.

When comparing Figure 6, Figure 7 and Figure 8, it can be seen clearly that the computed bounds for the electric current $i_{\mathrm{L}}$ obtained with the help of verifyode are tighter than those of the proposed approach but that the interval bounds for the capacitor’s voltage $u_{\mathrm{C}}$ are much wider if the Taylor model solver is employed. Especially the fact that this simulation is practically not able to predict the correct sign of the voltage $u_{\mathrm{C}}$ is counterproductive for practical applications where the simulation results may be employed to forecast the magnitude and direction of power flows towards a consumer (represented by the resistance $R_{\mathrm{S}}$ in the considered scenario). The new solution approach (without counting the offline parameter splitting) is faster by a factor of at least 100 than the alternative solver despite the use of multiple parameter intervals in the simulation if both are executed on the same computer using Matlab 2019b.

Remark 3.

Due to the linearity of the point-valued bounding systems produced by the new Algorithms 1–4, it is even possible (in the case when the dimension n is sufficiently small) to compute the worst-case state enclosures efficiently in closed form by using techniques for symbolic formula manipulation. The existence of such symbolic expressions can be exploited in future work to design fast model-predictive control strategies for uncertain dynamic systems.

Remark 4.

Due to cooperativity of the transformed system models—which are obtained with the approach presented in this paper—they are less affected by overestimation resulting from the wrapping effect, even if simulations were carried out for the corresponding interval boxes of the uncertain parameter and not for the decoupled bounding systems as suggested. Therefore, the obtained changes of coordinates could be employed also by general-purpose verified ODE solvers. This represents a further means to reduce overestimation in addition to the typically employed options such as QR preconditioning [11,34].

5. Conclusions and Future Work

In this paper, a novel subdivision-based transformation approach was presented that allows for computing the domains of reachable states of an uncertain dynamic system in terms of point-valued bounding models. This transformation significantly reduces computing times and overestimation in the computed results and outperforms general-purpose solvers such as verifyeig, as long as the assumption on limited parameter variabilities, detailed in this paper, are satisfied.

Methods, allowing the extension of the procedure further to mixed real and complex eigenvalues, resulting, for example from larger uncertainty in the resistances in the application considered in this paper will be investigated in future work. In addition, further research should be directed towards extensions for arbitrarily fast changes of parameters occurring in a time- or event-triggered framework as well as to interfacing the proposed methodology with techniques for gain scheduling control of uncertain dynamic systems [33]. Moreover, possible generalizations towards the simulation of uncertain fractional-order differential equations will be investigated [46,47,48].