A novel technique for stability analysis of linear time-invariant and nonlinear continuous-time control systems based on Markov parameters and Hankel matrices

Abstract

In this paper, we present an innovative approach for assessing the stability of linear time-invariant (LTI) systems. The key innovations in our methodology include the calculation of Markov parameters, the extension to nonlinear continuous-time systems, and significantly improved computational efficiency. Our novel approach leverages Markov parameters derived from the characteristic polynomial of the system. Unlike traditional methods, we employ the Leverrier-Takeno algorithm to efficiently determine the coefficients of the characteristic polynomial for state-space representation. The systematic calculation of Markov parameters using a specialized Maclaurin series expansion enhances our method’s precision. One of our major achievements is the rigorous verification of stability through Hankel matrices, ensuring positive definiteness using the Cholesky decomposition algorithm. In addition, we introduce a second stability criterion based on iterated polynomial long divisions, although we acknowledge its limitations in handling high-order systems due to computational inefficiency. Our pioneering contributions extend beyond LTI systems. We adapt our innovative techniques to analyze the local stability of nonlinear continuous-time systems by introducing the concept of Jacobian matrices and employing the indirect Lyapunov method. This expansion broadens the applicability of our approach to a wider range of real-world systems. Notably, our approach exhibits remarkable computational efficiency. According to CPU-time analysis, our first stability criterion outperforms the classical Routh–Hurwitz method by more than 12 times for dynamic systems with orders ranging from 2 to 180. This significant reduction in computation time positions our method as a valuable option for real-time stability analysis applications.

Keywords

Stability analysis method linear time-invariant systems Markov parameters Hankel matrices Cholesky decomposition Routh–Hurwitz method

Introduction

The stability of control systems is a crucial property that affects their performance, safety, and reliability. Therefore, a considerable amount of control engineering literature deals with the principles and techniques for assessing different stability definitions. The bounded-input-bounded-output (BIBO) stability is one of these concepts that plays a key role in the design and analysis of control systems. It ensures that the dynamic system will behave predictably and will not exhibit any unexpected behavior (Corriou, 2013). In fact, the nature of the roots of the characteristic polynomial associated with a system dictates its BIBO stability. If all the mentioned roots have negative real parts, the BIBO stability of the system is guaranteed. However, it is naïve and simplistic to solely rely on calculating the characteristic polynomial roots to determine a system’s stability because (a) there is no exact algebraic root-finding technique for polynomial equations of order five or higher, according to the Abel–Ruffini theorem, aka the “impossibility theorem,” (b) even using numerical root-finding algorithms, the results will sometimes suffer from large magnitudes of error if the polynomial is ill-conditioned (Cheney and Kincaid, 2009; Cohen, 1994).

Consequently, several stability assessment methods have been developed, which are based on identifying the regional location of the characteristic roots in the complex plane rather than their exact values. Examples of such methods are the Hurwitz criterion, the Routh–Hurwitz method, and Evans root locus method. Unfortunately, each of these methods has their own drawbacks. For instance, the Hurwitz and the Routh–Hurwitz methods are only applicable to real characteristic polynomials. They also cannot handle systems whose closed-loop transfer function includes non-algebraic terms, such as time delays. Furthermore, both methods are inefficient for higher-order systems from a computational viewpoint. The Evans root locus method is not applicable to multi-input-multi-output (MIMO) systems and cannot be applied to time-varying models. It can be time-consuming as a graphical method, especially for complicated higher-order systems and it does not account for time delays (Abolpour et al., 2023; Cogan et al., 2009; Fatoorehchi and Djilali, 2023; Kavitha and Ramakalyan, 2014; Shinskey, 1996). It is worthwhile to mention that Yang et al. (2022) have generalized the Routh–Hurwitz technique to fractional-order linear systems, and their method is called the fractional-order Routh–Hurwitz stability criterion.

The most recent BIBO stability criterion for LTI systems is an integrative method that is founded on an efficient approach for computing the Bézoutian matrix associated with the system under investigation, aiming subsequently to derive the system’s Hermite–Fujiwara matrix (Fatoorehchi and Ehrhardt, 2023).

Table 1 provides a brief overview of the advantages and drawbacks of the conventional stability criteria, facilitating quick comparison.

Table 1.

Comparison of stability analysis methods.

Stability analysis method	Advantages	Disadvantages
Hurwitz criterion	- Simple and systematic	- Cannot handle non-algebraic terms in the transfer function - Computation of n determinants is required - Cannot be used for nonlinear systems
Routh–Hurwitz criterion	- Improved computational efficiency compared with Hurwitz criterion	- Cannot handle non-algebraic terms in the transfer function - Computation of n determinants is required - Cannot be used for nonlinear systems
Evans root locus method	- Algorithmic nature of the method - Visualization of the roots’ movement	- Zeros and poles of the system must be calculated - The graphical nature of the method - Cannot be used for nonlinear systems
Bode’s stability criterion	- It is based on frequency response analysis - It requires simple mathematical operations	- The graphical nature of the method - Cannot be used for nonlinear systems
Fatoorehchi–Djilali method	- Applicable to state-space represented systems - It requires only simple matrix operations	- Cannot handle time delays in the system - It is not extended to account for nonlinear systems
Fatoorehchi–Ehrhardt method	- Applicable to state-space represented systems - It is extended to nonlinear time-continuous systems	- It requires advanced mathematics knowledge - Cannot handle time delays in the system
The proposed method	- Computational efficiency - Applicable to state-space represented systems - It is extended to nonlinear time-continuous systems	- Cannot handle time delays in the system - It is not a single formula and comprised several steps

It should be noted that fast or real-time stability analysis is particularly favored in practice for several industries. For these applications, the computationally efficient stability criteria are often preferred over more complex methods that require significant computational resources and time efforts. For instance, real-time stability assessment of power grid systems, with a typical frequency of 10 Hz, is essential to a transient-free power management. Furthermore, very fast stability analyses are required to control intensified unit operations of chemical engineering, such as spinning disk reactors or multiloop chemical reactors (Ghiasy et al., 2012; Sun et al., 2016; Yakhnin and Menzinger, 2004).

In this paper, we will develop a novel stability criterion for linear time-invariant (LTI) systems by exploiting the Maclaurin series expansion for evaluating the Markov parameters of rational functions. In the course of our algorithm, we will take benefit from the Leverrier-Takeno algorithm for calculating the characteristic polynomial, if the system is given in the state-space form, and the Cholesky decomposition for checking positive definiteness of the associated Hankel matrices. We will demonstrate that our approach is computationally efficient and valid for MIMO systems. A set of case studies involving distillation columns and a nonlinear boiler-turbine unit is given to demonstrate the methodology and application of our proposed method.

Mathematical background

In this section, we present the mathematical background that serves as the foundation for the methodology described in the subsequent section. Specifically, we introduce the necessary definitions and theorems. By laying out this mathematical background, we aim to provide a comprehensive and self-contained exposition of our approach, that is accessible for a broad audience.

Definition 2.1. Matrix $A \in C^{n \times n}$ is called Hermitian if it is equal to its own conjugate transpose.

Definition 2.2. A Hermitian matrix $A \in C^{n \times n}$ is called positive definite if $x^{*} A x > 0$ for all nonzero $x \in C^{n}$ . Note that asterisk denotes conjugate transpose. For real matrices, the definition is the same but conjugate transpose is replaced by transpose, that is the inequality changes to $x^{T} A x > 0$ .

Definition 2.3. A Hermitian matrix $A \in C^{n \times n}$ is called negative definite if $x^{*} A x < 0$ for all nonzero $x \in C^{n}$ . In case of real matrices, conjugate transpose should be changed to ordinary transpose.

Lemma 2.1. If we multiply a negative definite matrix by minus one, then the resulting matrix is positive definite.

Proof: The proof is simple. We just need to replace A by –A in the inequality of Definition 2.2, and then in view of Definition 2.3, conclude the statement.

Definition 2.4. The Cholesky decomposition of matrix $A \in C^{n \times n}$ is a factorization that allows for $A = R^{T} R$ , where $R = (r_{i, j}) \in C^{n \times n}$ is nonsingular, upper triangular with positive diagonal entries that are calculated by the following pseudocode

\begin{matrix} for j = 1 to n \\ for i = 1 to j - 1 \\ r_{i, j} = \frac{a_{i, j} - \sum_{k = 1}^{i - 1} r_{k, i} r_{k, j}}{r_{i, i}} \\ next \\ r_{j, j} = \sqrt{a_{j, j} - \sum_{k = 1}^{j - 1} r_{k, j}^{2}} \\ next \end{matrix}

Theorem 2.1. Matrix $A \in C^{n \times n}$ is positive definite if and only if it has a Cholesky decomposition (Roy and Banerjee, 2014).

Theorem 2.2. (Leverrier-Takeno’s method (Wayland, 1945)) Let $P (x) = x^{n} + \sum_{i = 1}^{n} c_{i} x^{n - i}$ be the characteristic polynomial associated with the matrix A.

Furthermore, let $σ_{k} = trace (A^{k}), k = 1, 2, \dots, n$ . The coefficients $c_{i}$ can be calculated using the following recursive relation

\begin{matrix} r a_{r} + σ_{1} a_{r - 1} + σ_{2} a_{r - 2} + \dots + σ_{r - 1} a_{1} + σ_{r} = 0, \\ r = 1, 2, \dots, n \end{matrix}

(1)

Methodology overview

Let the characteristic polynomial associated with the system matrix $A \in R^{n \times n}$ be given by

P (x) = x^{n} + a_{n - 1} x^{n - 1} + \dots + a_{1} x + a_{0}

(2)

We define the rational function $Q (x)$ as follows

\begin{matrix} Q (x) = \frac{P (- x)}{P (x)} = \frac{{(- 1)}^{n} x^{n} + {(- 1)}^{n - 1} a_{n - 1} x^{n - 1} + \dots - a_{1} x + a_{0}}{x^{n} + a_{n - 1} x^{n - 1} + \dots + a_{1} x + a_{0}} \end{matrix}

(3)

The Markov parameters associated with the rational function $Q (x)$ are denoted by $s_{- 1}, s_{0}, s_{1}, \dots$ and are defined by

Q (x) = s_{- 1} + \frac{s_{0}}{x} + \frac{s_{1}}{x^{2}} + \frac{s_{2}}{x^{3}} + \frac{s_{3}}{x^{4}} + \dots

(4)

where the convergence of the right-hand side of equation (4) does not play a role in defining $s_{i}$ .

For calculation of the Markov parameters, we will take advantage of the Maclaurin expansion. To do so, let us assume that the function $Q (1 / x)$ has a Maclaurin expansion and is represented by the subsequent power series

Q (1 / x) = \frac{P (- 1 / x)}{P (1 / x)} = b_{0} + b_{1} x + b_{2} x^{2} + b_{3} x^{3} + \dots

(5)

Next, by a simple reciprocal variable change, we can deduce that

Q (x) = b_{0} + \frac{b_{1}}{x} + \frac{b_{2}}{x^{2}} + \frac{b_{3}}{x^{3}} + \dots

(6)

If we compare equations (4) and (6), then we readily obtain the Markov parameters of the rational function $Q (x)$ as

\begin{matrix} s_{- 1} = b_{0}, \\ s_{0} = b_{1}, \\ s_{1} = b_{2}, \\ s_{2} = b_{3}, \\ ⋮ \end{matrix}

(7)

Subsequently, with the help of the Markov parameters, we can formally define the infinite Hankel matrix as

H = [\begin{matrix} s_{0} & s_{1} & s_{2} & s_{3} & \cdot \\ s_{1} & s_{2} & s_{3} & \cdot & \cdot \\ s_{2} & s_{3} & \cdot & \cdot & \cdot \\ s_{3} & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \end{matrix}]

(8)

It is clear from equation (8) that each Markov parameter fills an anti-diagonal of the infinite Hankel matrix.

Next, we define the submatrix $H_{(n, n)}$ as the one consisting of the first n rows and columns of the infinite Hankel matrix. In addition, we introduce the subsequent n-by-n lower triangular matrix L as

L = [\begin{matrix} s_{- 1} & 0 & 0 & \dots & 0 \\ s_{0} & - s_{- 1} & 0 & \dots & 0 \\ s_{1} & - s_{0} & s_{- 1} & \dots & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋱ & 0 \\ s_{n - 2} & - s_{n - 3} & s_{n - 4} & \dots & {(- 1)}^{n - 1} s_{- 1} \end{matrix}]

(9)

such that $L^{2} = I$ , where I denotes the identity matrix of dimension n.

Finally, we introduce the matrix ${\bar{H}}_{(n, n)}$ as

{\bar{H}}_{(n, n)} = L H_{(n, n)}

(10)

by which the stability of the system under study can be assessed according to the subsequent theorem.

Remark 3.1. The matrix ${\bar{H}}_{(n, n)}$ is symmetric.

Lemma 3.2. Multiplying all Markov parameters by a positive real number will not affect the stability status of the system.

Theorem 3.1. Assuming that the matrix $H_{(n, n)}$ is nonsingular, all the zeroes of the characteristic polynomial (2) have positive real parts if and only if the matrix ${\bar{H}}_{(n, n)}$ is positive definite. Furthermore, all the zeroes of equation (2) have negative real parts if and only if ${\bar{H}}_{(n, n)}$ is negative definite (Datta, 1979).

Corollary 3.1. (Stability criterion) The dynamic (control) system whose characteristic polynomial corresponds to equation (2) is BIBO stable if and only if matrix $- {\bar{H}}_{(n, n)}$ , as defined by equation (10), can be Cholesky-decomposed. Note that the minus sign is due to Lemma 2.1, Theorem 2.1, and Theorem 3.1.

Considering equations (8)–(10), it is important to underscore that our stability criterion relies solely on the first n Markov parameters, indicating that the subsequent Markov parameters (beyond the nth) hold no significance in the stability status of the studied system.

Another stability criterion can be formulated based on a distinct set of Markov parameters and an alternative definition of the Hankel matrix (see Theorem 3.2).

Theorem 3.2. Let $P (x)$ be real and denote the nth-order characteristic polynomial of an LTI system. Also, let $M (x)$ and $N (x)$ be of even degree and the odd degree parts of $P (x)$ , respectively. If n is even, then we expand the following rational function

\frac{N (x)}{M (x)} = \frac{δ_{1}}{x} - \frac{δ_{3}}{x^{3}} + \frac{δ_{5}}{x^{5}} - \dots - \frac{δ_{2 n - 1}}{x^{2 n - 1}}

(11)

and if n is odd, then we construct the following expansion

\frac{M (x)}{N (x)} = \frac{δ_{1}}{x} - \frac{δ_{3}}{x^{3}} + \frac{δ_{5}}{x^{5}} - \dots - \frac{δ_{2 n - 1}}{x^{2 n - 1}}

(12)

where $δ_{i}$ are the associated Markov parameters.

Now, the system characterized by $P (x)$ is stable if and only if the following particular Hankel matrix is positive definite (Katbab and Jury, 1990)

\hat{H} = [\begin{matrix} δ_{1} & δ_{2} & δ_{3} & \dots & δ_{n} \\ δ_{2} & δ_{3} & δ_{4} & \dots & δ_{n + 1} \\ δ_{3} & δ_{4} & δ_{5} & \dots & δ_{n + 2} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ δ_{n} & δ_{n + 1} & \dots & δ_{2 n - 2} & δ_{2 n - 1} \end{matrix}]

(13)

Corollary 3.2. The system with the characteristic polynomial, as mentioned in Theorem 3.2, is BIBO stable if and only if the matrix $\hat{H}$ , as in equation (13), is Cholesky-decomposable.

Note that the expansions in the right-hand sides of equations (11) and (12) should be derived from repeated polynomial long divisions, which is a time-intensive and tedious task especially for higher-order systems. Thus, we use Corollary 3.2 as our second stability criterion only for systems with a small number of state functions in the next section and rely on Corollary 3.1 as our main stability test.

In the end of this section, we emphasize that our stability criteria can be extended to nonlinear continuous-time systems for local stability analysis by replacing the state matrix A by the Jacobian matrix of the nonlinear system. This approach is known as the indirect Lyapunov criterion.

Theorem 3.3. Let $x_{e}$ be an equilibrium point of the following nonlinear autonomous system

\overset{\cdot}{x} = f (x)

(14)

such that

f (x_{e}) = 0

(15)

and also, let $J (x)$ denote the Jacobian matrix of $f$ in the point $x$ defined by

J (x) = [\begin{matrix} \frac{\partial f_{1}}{\partial x_{1}} & \dots & \frac{\partial f_{1}}{\partial x_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{n}}{\partial x_{1}} & \dots & \frac{\partial f_{n}}{\partial x_{n}} \end{matrix}]

(16)

It holds that if all eigenvalues of matrix $J (x_{e})$ have negative real parts, then system (14) is stable at $x_{e}$ (Adamy, 2022).

Remark 3.2. We can apply a simple variable change of type $x = y - h$ , where $h$ is a vector of appropriate constants, to convert system (14) to

\overset{\cdot}{y} = g (y)

(17)

whose equilibrium point is different from that of equation (14) and actually is transferred from $x_{e}$ to somewhere else. In this way, we can assess the local stability of any arbitrary point of system (14) with the help of Theorem 3.3.

For better illustration, the flowchart of our first and second stability criteria, the ones based on Corollary 3.1 and 3.2, is depicted in Figures 1 and 2.

Figure 1.

Flowchart of the first proposed stability criterion (Corollary 3.1).

Figure 2.

Flowchart of the second proposed stability criterion (Corollary 3.2).

Illustrative examples

Example 1. Binary distillation column with three trays

Rasmussen and Jørgensen (1999) modeled a binary distillation column comprising a total condenser, a reboiler, and three equilibrium (theoretical) trays. The dynamics of the column is governed by a system of five nonlinear ordinanry differential equations (ODEs) that are derived from appropriate mole balances

Reboiler : H_{r} \frac{d x_{1}}{dt} = (L + L_{f}) (x_{2} - x_{1}) + V (x_{1} - y_{1})

(18)

\begin{matrix} Theoretical Tray (No . 2) : H_{t} \frac{d x_{2}}{dt} = (L + L_{f}) (x_{3} - x_{2}) \\ + V (y_{1} - y_{2}) \end{matrix}

(19)

\begin{matrix} Theoretical Feed Tray (No . 3) : H_{t} \frac{d x_{3}}{dt} = L_{f} x_{f} + L x_{4} \\ - (L + L_{f}) x_{3} + V (y_{2} - y_{3}) \end{matrix}

(20)

\begin{matrix} Theoretical Tray (No . 4) : H_{t} \frac{d x_{4}}{dt} = L (x_{5} - x_{4}) + V (y_{3} - y_{4}) \end{matrix}

(21)

Condenser : H_{c} \frac{d x_{5}}{dt} = V (y_{4} - x_{5})

(22)

in which molar compositions of the liquid and vapor phases in the trays and the reboiler are subject to thermodynamic equilibrium simplistically modeled by

y_{i} = \frac{α x_{i}}{1 + (α - 1) x_{i}}

(23)

where $α$ is the relative volatility of the light-key and the heavy-key components averaged along the column. Actually, it is equation (23) that introduces the nonlinearity into the dynamic system. Constant molar flow rates were assumed for each phase. Furthermore, it was presumed that the liquid levels in the condenser and the reboiler are ideally controlled. For a specific set of conditions and input values, listed in Table 2, the aforementioned system has been linearized and the associated state matrix was given as

A = [\begin{matrix} - 3.9138 & 1.2459 & 0 & 0 & 0 \\ 5.6208 & - 4.5159 & 1.8688 & 0 & 0 \\ 0 & 2.6471 & - 2.8439 & 1.3688 & 0 \\ 0 & 0 & 0.9752 & - 1.8350 & 1.3688 \\ 0 & 0 & 0 & 0.3108 & - 1.0792 \end{matrix}]

(24)

Table 2.

Model parameters of the distillation column of Example 1.

$H_{c}$	Condenser holdup	30 mol
$H_{c}$	Reboiler holdup	30 mol
$H_{t}$	Tray holdup	20 mol
$α$	Relative volatility	5
$L_{f}$	Feed flow rate	10 mol/min
$x_{f}$	Molar composition of feed (light-keycomponent)	0.5
$L$	Liquid molar flow rate	27.3755 mol/min
$V$	Vapor molar flow rate	32.3755 mol/min

We examine the stability of the system at such operating conditions subsequently.

Using the Leverrier-Takeno method (Theorem 2.2), we derive the characteristic polynomial as

\begin{matrix} P (x) = x^{5} + 14.1878 x^{4} + 62.7713 x^{3} \\ + 102.3686 x^{2} + 51.8076 x + 1.7111 \end{matrix}

(25)

Expanding the function $Q (1 / x) = P (- 1 / x) / P (1 / x)$ using the Maclaurin theorem, we have

\begin{matrix} Q (1 / x) = \frac{P (- 1 / x)}{P (1 / x)} = - 1 + 28.3756 x - 402.5873 x^{2} \\ + 4135.3925 x^{3} - 36305.9620 x^{4} \\ + 295263.4155 x^{5} - 2312691.6661 x^{6} \\ + 17780972.9959 x^{7} - 135454080.0770 x^{8} \\ + 1027172849.3988 x^{9} + O (x^{10}) \end{matrix}

(26)

In view of equations (5), (7), and (26), we identify the Markov parameters of $Q (x)$ as listed below

\begin{matrix} \begin{matrix} s_{- 1} = - 1, \\ s_{0} = 28.3756, \\ s_{1} = - 402.5873, \\ s_{2} = 4135.3925, \\ s_{3} = - 36305.9620, \\ s_{4} = 295263.4155, \end{matrix} \begin{matrix} s_{5} = - 2312691.6661, \\ s_{6} = 17780972.9959, \\ s_{7} = - 135454080.0770, \\ s_{8} = 1027172849.3988, \\ ⋮ \end{matrix} \end{matrix}

(27)

Afterward, the infinite Hankel matrix is formally constructed as

H = [\begin{matrix} 28.3756 & - 402.5873 & 4135.3925 & - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 \\ - 402.5873 & 4135.3925 & - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot \\ 4135.3925 & - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot & \cdot \\ - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot & \cdot & \cdot \\ 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot & \cdot & \cdot & \cdot \\ - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 17780972.9959 & - 135454080.0770 & 1027172849.3988 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ - 135454080.0770 & 1027172849.3988 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{matrix}]

(28)

and from which the submatrix $H_{(5, 5)}$ is extracted

H_{(5, 5)} = [\begin{matrix} 28.3756 & - 402.5873 & 4135.3925 & - 36305.9620 & 295263.4155 \\ - 402.5873 & 4135.3925 & - 36305.9620 & 295263.4155 & - 2312691.6661 \\ 4135.3925 & - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 \\ - 36305.9620 & 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 \\ 295263.4155 & - 2312691.6661 & 17780972.9959 & - 135454080.0770 & 1027172849.3988 \end{matrix}]

(29)

Then, we define matrix L according to equation (9)

L = [\begin{matrix} - 1 & 0 & 0 & 0 & 0 \\ 28.3756 & 1 & 0 & 0 & 0 \\ - 402.5873 & - 28.3756 & - 1 & 0 & 0 \\ 4135.3925 & 402.5873 & 28.3756 & 1 & 0 \\ - 36305.9620 & - 4135.3925 & - 402.5873 & - 28.3756 & - 1 \end{matrix}]

(30)

Finally, by means of equation (10), we calculate matrix ${\bar{H}}_{(5, 5)}$

{\bar{H}}_{(5, 5)} = [\begin{matrix} - 28.3756 & 402.5873 & - 4135.3925 & 36305.9620 & - 295263.4155 \\ 402.5873 & - 7288.2647 & 81038.2814 & - 734940.0398 & 6065584.9067 \\ - 4135.3925 & 81038.2830 & - 929916.6163 & 8550735.6767 & - 71026271.7158 \\ 36305.9611 & - 734940.0522 & 8550735.5640 & - 79113132.3065 & 659061439.157 \\ - 295263.4037 & 6065585.0287 & - 71026270.6536 & 659061438.425 & - 5497911146.02 \end{matrix}]

(31)

Note that the calculated ${\bar{H}}_{(5, 5)}$ is approximately symmetric. The very minor discrepancies between some of symmetric pairs of the entries are due to round-off errors.

Subsequently, we try the Cholesky decomposition algorithm on the matrix $- {\bar{H}}_{(5, 5)}$ and obtain the following diagonal components of the matrix R, the one mentioned in Definition 2.4, as listed below

\begin{matrix} r_{1, 1} = 5.3268, \\ r_{2, 2} = 39.7044, \\ r_{3, 3} = 99.5377, \\ r_{4, 4} = 95.0096, \\ r_{5, 5} = 11.5668 . \end{matrix}

(32)

Since all diagonal entries $r_{j, j}$ are positive nonzero real numbers, we conclude that the Cholesky decomposition algorithm has run to completion and thus, the matrix $- {\bar{H}}_{(5, 5)}$ is positive definite. It immediately follows that the matrix ${\bar{H}}_{(5, 5)}$ is negative definite and hence, the system under study is deemed as stable in view of Corollary 3.1. This assessment is corroborated by knowing the eigenvalues of the matrix A are {–7.5492, –3.7373, –2.0185, –0.0355, –0.8473}, which obviously are left half-plane-located in the complex plane.

Now, we will examine the stability of the system by the second approach, that is, Corollary 3.2.

In view of Theorem 3.2 and equation (12), it is clear that we have to expand $M (x) / N (x)$ where

M (x) = 14.1878 x^{4} + 102.3686 x^{2} + 1.7111

(33)

and

N (x) = x^{5} + 62.7713 x^{3} + 51.8076 x

(34)

After a series of polynomial long divisions, we will derive the following nine Markov parameters according to equation (12)

\begin{matrix} \begin{matrix} δ_{1} = 14.1878, \\ δ_{2} = 0, \\ δ_{3} = 788.2180, \\ δ_{4} = 0, \\ δ_{5} = 48744.1469, \\ δ_{6} = 0, \end{matrix} \begin{matrix} δ_{7} = 3018897.7843, \\ δ_{8} = 0, \\ δ_{9} = 186974821.2229 . \end{matrix} \end{matrix}

(35)

Afterward, the relevant Hankel matrix is derived as

\hat{H} = [\begin{matrix} 14.1878 & 0 & 788.2180 & 0 & 48744.1469 \\ 0 & 788.2180 & 0 & 48744.1469 & 0 \\ 788.2180 & 0 & 48744.1469 & 0 & 3018897.7843 \\ 0 & 48744.1469 & 0 & 3018897.7843 & 0 \\ 48744.1469 & 0 & 3018897.7843 & 0 & 186974821.2229 \end{matrix}]

(36)

Following that, we try the Cholesky decomposition algorithm on matrix $\hat{H}$ to yield the Cholesky upper triangular factor $R$ whose diagonal entries are

\begin{matrix} r_{1, 1} = 3.7666, \\ r_{2, 2} = 28.0752, \\ r_{3, 3} = 70.3837, \\ r_{4, 4} = 67.1848, \\ r_{5, 5} = 9.3068 . \end{matrix}

(37)

The fact that all the diagonal components are real, shows that the Cholesky algorithm has run to completion and thus, matrix $\hat{H}$ has a Cholesky decomposition. This, in view of Corollary 3.2, classifies the dynamic model of the distillation column as stable.

Example 2. Binary distillation column with eight trays

Drgoňa et al. (2017) developed a data-driven model for a laboratory distillation column as depicted in Figure 3 for separating methanol and water. The column consisted of eight sieve trays, a reboiler, and a total condenser. An LTI model in the discrete domain was constructed by exerting a sufficient number of inputs to the experimental system and collecting the output responses. The model is expressed as

{\begin{matrix} x_{k + 1} = A x_{k} + B u_{k}, \\ y_{k} = C x_{k} . \end{matrix}

(38)

Figure 3.

Schematic view of the distillation column of Example 2.

The system matrix for the dynamics of the distillation column at a particular operation point is

A = [\begin{matrix} 0.9970 & - 0.0008 & - 0.0001 & - 0.0009 & 0.0020 & - 0.0006 & 0.0005 & - 0.0019 & - 0.0010 & - 0.0025 \\ - 0.0006 & 0.9999 & 0.0002 & 0.0006 & - 0.0002 & 0.0003 & 0.0021 & - 0.0018 & - 0.0018 & - 0.0018 \\ 0.0013 & 0.0003 & 0.9977 & 0.0018 & - 0.0002 & - 0.0007 & - 0.0023 & 0.0020 & 0.0038 & 0.0041 \\ 0.0032 & - 0.0004 & - 0.0018 & 0.9940 & - 0.0036 & 0.0318 & - 0.0008 & 0.0042 & 0.0019 & 0.0121 \\ 0.0007 & 0.0003 & 0.0002 & - 0.0163 & 0.8041 & 0.5906 & - 0.0002 & 0.0625 & 0.0191 & 0.0723 \\ 0.0007 & - 0.0006 & 0.0004 & - 0.0261 & - 0.5842 & 0.8089 & - 0.0008 & - 0.0292 & - 0.0229 & - 0.0185 \\ 0.0011 & - 0.0016 & 0.0023 & - 0.0023 & - 0.0025 & - 0.0001 & 0.9911 & 0.0141 & - 0.0086 & 0.0149 \\ 0.0019 & 0.0002 & 0.0016 & - 0.0063 & - 0.0106 & - 0.0130 & - 0.0157 & 0.2994 & - 0.9139 & 0.2570 \\ - 0.0003 & 0.0007 & - 0.0018 & 0.0062 & - 0.0254 & 0.0013 & - 0.0079 & 0.8916 & 0.1581 & - 0.3229 \\ - 0.0001 & - 0.0003 & 0.0034 & - 0.0046 & 0.0300 & 0.0145 & 0.0080 & 0.0741 & 0.0116 & - 0.4568 \end{matrix}]

(39)

The characteristic polynomial of the matrix A can be determined via the Leverrier-Takeno method (as stated in Theorem 2.2)

\begin{matrix} P (x) = x^{10} - 6.5934 x^{9} + 19.5889 x^{8} - 34.6816 x^{7} + 40.3842 x^{6} - 31.3894 x^{5} + 14.5998 x^{4} \\ - 1.4665 x^{3} - 2.8818 x^{2} + 1.8062 x - 0.3664 . \end{matrix}

(40)

To streamline the presentation, here we will skip some of the less crucial steps that were demonstrated in the previous example and list the related Markov parameters as

\begin{matrix} \begin{matrix} s_{- 1} = 1, \\ s_{0} = 13.1868, \\ s_{1} = 86.9458, \\ s_{2} = 384.3166, \\ s_{3} = 1288.1160, \\ s_{4} = 3510.3749, \end{matrix} \begin{matrix} s_{5} = 8143.9002, \\ s_{6} = 16624.7247, \\ s_{7} = 30622.5636, \\ s_{8} = 51910.4279, \\ ⋮ \end{matrix} \end{matrix}

(41)

Thus, the infinite Hankel matrix is formally determined accordingly

H = [\begin{matrix} 13.1868 & 86.9458 & 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 \\ 86.9458 & 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & \cdot \\ 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & \cdot & \cdot \\ 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & \cdot & \cdot & \cdot \\ 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & \cdot & \cdot & \cdot & \cdot \\ 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & \cdot & \cdot & \cdot & \cdot & \cdot \\ 16624.7247 & 30622.5636 & 51910.4279 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 30622.5636 & 51910.4279 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{matrix}]

(42)

Then, it is straightforward to generate the submatrix $H_{(10, 10)}$ from H as

H_{(10, 10)} = [\begin{matrix} 13.1868 & 86.9458 & 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 \\ 86.9458 & 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 \\ 384.3166 & 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 \\ 1288.1160 & 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 \\ 3510.3749 & 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 \\ 8143.9002 & 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 & 432152.5154 \\ 16624.7247 & 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 & 432152.5154 & 556173.8388 \\ 30622.5636 & 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 & 432152.5154 & 556173.8388 & 704309.5816 \\ 51910.4279 & 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 & 432152.5154 & 556173.8388 & 704309.5816 & 879685.6268 \\ 82263.2862 & 123414.9337 & 177083.3286 & 245058.0263 & 329321.0510 & 432152.5154 & 556173.8388 & 704309.5816 & 879685.6268 & 1085504.6374 \end{matrix}]

(43)

The next leading step is to compute the L matrix

L = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 13.1868 & - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 86.9458 & - 13.1868 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 384.3166 & - 86.9458 & 13.1868 & - 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1288.1160 & - 384.3166 & 86.9458 & - 13.1868 & 1 & 0 & 0 & 0 & 0 & 0 \\ 3510.3749 & - 1288.1160 & 384.3166 & - 86.9458 & 13.1868 & - 1 & 0 & 0 & 0 & 0 \\ 8143.9002 & - 3510.3749 & 1288.1160 & - 384.3166 & 86.9458 & - 13.1868 & 1 & 0 & 0 & 0 \\ 16624.7247 & - 8143.9002 & 3510.3749 & - 1288.1160 & 384.3166 & - 86.9458 & 13.1868 & - 1 & 0 & 0 \\ 30622.5636 & - 16624.7247 & 8143.9002 & - 3510.3749 & 1288.1160 & - 384.3166 & 86.9458 & - 13.1868 & 1 & 0 \\ 51910.4279 & - 30622.5636 & 16624.7247 & - 8143.9002 & 3510.3749 & - 1288.1160 & 384.3166 & - 86.9458 & 13.1868 & - 1 \end{matrix}]

(44)

Afterward, we construct ${\bar{H}}_{(10, 10)}$ through a multiplication of L by $H_{(10, 10)}$

{\bar{H}}_{(10, 10)} = [\begin{matrix} - 13.1868 & - 86.9458 & - 384.3166 & - 1288.1160 & - 3510.3749 & - 8143.9002 & - 16624.7247 & - 30622.5636 & - 51910.4279 & - 82263.2862 \\ - 86.9458 & - 762.2208 & - 3779.7901 & - 13475.7533 & - 38146.7125 & - 90767.2595 & - 188604.3561 & - 351903.1947 & - 602269.1448 & - 961374.5696 \\ - 384.3166 & - 3779.7901 & - 19938.9793 & - 73849.6257 & - 214445.2671 & - 519473.9530 & - 1093547.5783 & - 2060235.5948 & - 3552021.5615 & - 5702086.3923 \\ - 1288.1160 & - 13475.7533 & - 73849.6257 & - 280599.1050 & - 829621.4367 & - 2036288.5159 & - 4328922.1347 & - 8216738.0905 & - 14247952.9414 & - 22974835.1973 \\ - 3510.3749 & - 38146.7125 & - 214445.2671 & - 829621.4367 & - 2485481.7270 & - 6161366.2645 & - 13197836.0888 & - 25197461.8121 & - 43891818.4801 & - 71028698.4398 \\ - 8143.9002 & - 90767.2595 & - 519473.9530 & - 2036288.5159 & - 6161366.2645 & - 15390307.7464 & - 33161555.9299 & - 63604862.9929 & - 111196369.2000 & - 180462053.3207 \\ - 16624.7247 & - 188604.3561 & - 1093547.5783 & - 4328922.1347 & - 13197836.0888 & - 33161555.9299 & - 71785237.3239 & - 138190735.8521 & - 242291295.5626 & - 394125871.0347 \\ - 30622.5636 & - 351903.1947 & - 2060235.5948 & - 8216738.0905 & - 25197461.8121 & - 63604862.9929 & - 138190735.8521 & - 266800395.4622 & - 468870703.0042 & - 764110816.8748 \\ - 51910.4279 & - 602269.1448 & - 3552021.5615 & - 14247952.9414 & - 43891818.4801 & - 111196369.2000 & - 242291295.5626 & - 468870703.0042 & - 825519567.9594 & - 1347346264.2244 \\ - 82263.2862 & - 961374.5696 & - 5702086.3923 & - 22974835.1973 & - 71028698.4398 & - 180462053.3207 & - 394125871.0347 & - 764110816.8748 & - 1347346264.2244 & - 2201676769.3487 \end{matrix}]

(45)

Finally, we apply the Cholesky decomposition on $- {\bar{H}}_{(10, 10)}$ and obtain the following diagonal elements of matrix R (see Definition 2.4)

\begin{matrix} \begin{matrix} r_{1, 1} = 3.6314 i, \\ r_{2, 2} = 13.7460 i, \\ r_{3, 3} = 22.8941 i, \\ r_{4, 4} = 24.3541 i, \\ r_{5, 5} = 18.8121 i . \end{matrix} \begin{matrix} r_{6, 6} = 11.5439 i, \\ r_{7, 7} = 6.3266 i, \\ r_{8, 8} = 3.1927 i, \\ r_{9, 9} = 0.7944 i, \\ r_{10, 10} = 0.3858 . \end{matrix} \end{matrix}

(46)

Since there are non-real entries on the diagonal of the matrix R, we deduce that the considered system is unstable. We can confirm the truth of this statement by observing that at least one eigenvalue of the system matrix A has a positive real part. The eigenvalues are {–0.4466, 0.9975, 0.9996, 0.2219 ± 0.8887i, 0.8078 ± 0.5871i, 0.9913, 0.9960 ± 0.0030i}.

Let us now apply the second stability test outlined in Corollary 3.2 to system (38). As the characteristic polynomial equation (40) is even-ordered, we should expand the rational function $N (x) / M (x)$ knowing that

\begin{matrix} M (x) = x^{10} + 19.5889 x^{8} + 40.3842 x^{6} + 14.5998 x^{4} \\ - 2.8818 x^{2} - 0.3664 \end{matrix}

(47)

and

\begin{matrix} N (x) = - 6.5934 x^{9} - 34.6816 x^{7} - 31.3894 x^{5} - 1.4665 x^{3} \\ + 1.8062 x \end{matrix}

(48)

After repetitive polynomial long divisions, the 19 Markov parameters are read according to equation (11)

\begin{matrix} \begin{matrix} δ_{1} = - 6.5934, \\ δ_{2} = 0, \\ δ_{3} = - 94.4758, \\ δ_{4} = 0, \\ δ_{5} = - 1615.7982, \\ δ_{6} = 0, \\ δ_{7} = - 27931.1745, \\ δ_{8} = 0, \\ δ_{9} = - 483284.7885, \\ δ_{10} = 0, \end{matrix} \begin{matrix} δ_{11} = - 8362899.4302, \\ δ_{12} = 0, \\ δ_{13} = - 144715342.4459, \\ δ_{14} = 0, \\ δ_{15} = - 2504221129.7552, \\ δ_{16} = 0, \\ δ_{17} = - 43334203111.4795, \\ δ_{18} = 0, \\ δ_{19} = - 749875142566.8753 . \end{matrix} \end{matrix}

(49)

Next, the Hankel matrix $\hat{H}$ is formed as

\hat{H} = [\begin{matrix} - 6.5934 & 0 & - 94.4758 & 0 & - 1615.7982 & 0 & - 27931.1745 & 0 & - 483284.7885 & 0 \\ 0 & - 94.4758 & 0 & - 1615.7982 & 0 & - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 \\ - 94.4758 & 0 & - 1615.7982 & 0 & - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 \\ 0 & - 1615.7982 & 0 & - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 \\ - 1615.7982 & 0 & - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 \\ 0 & - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 & - 2504221129.7552 \\ - 27931.1745 & 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 & - 2504221129.7552 & 0 \\ 0 & - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 & - 2504221129.7552 & 0 & - 43334203111.4795 \\ - 483284.7885 & 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 & - 2504221129.7552 & 0 & - 43334203111.4795 & 0 \\ 0 & - 8362899.4302 & 0 & - 144715342.4459 & 0 & - 2504221129.7552 & 0 & - 43334203111.4795 & 0 & - 749875142566.8753 \end{matrix}]

(50)

On Cholesky decomposition, we notice that the matrix $\hat{H}$ is not positive definite because the diagonal of the upper triangular matrix R includes non-real entries

\begin{matrix} \begin{matrix} r_{1, 1} = 2.5677 i, \\ r_{2, 2} = 9.7198 i, \\ r_{3, 3} = 16.1885 i, \\ r_{4, 4} = 17.2208 i, \\ r_{5, 5} = 13.3022 i . \end{matrix} \begin{matrix} r_{6, 6} = 8.1629 i, \\ r_{7, 7} = 4.4735 i, \\ r_{8, 8} = 2.2576 i, \\ r_{9, 9} = 0.5615 i, \\ r_{10, 10} = 0.2731 . \end{matrix} \end{matrix}

(51)

As a result, system (38) at the given conditions cannot be stable.

Example 3. A hypothetical 13th-order system

As a demanding benchmark, we will test our method for a stable, hypothetical system with 13 states. The order of the system is high enough to be a fair representative of many chemical engineering unit operations. The characteristic polynomial of this system in the product form is

P (x) = Π_{j = 1}^{13} (x + j)

(52)

We can invoke Vieta’s formulas to determine the coefficients of our characteristic polynomial as follows

\begin{matrix} P (x) = x^{13} - e_{1} (- 1, - 2, \dots, - 13) x^{13} + e_{2} \\ (- 1, - 2, \dots, - 13) x^{12} + \dots + {(- 1)}^{13} \\ e_{13} (- 1, - 2, \dots, - 13) \end{matrix}

(53)

where e denotes elementary symmetric polynomials. Obviously, the system is stable as the roots of the characteristic equation, namely {–1, –2, –3, –4, –5, –6, –7, –8, –9, –10, –11, –12, –13}, all have negative real parts.

Skipping the less important steps in our procedure, we perform the Maclaurin expansion on the rational function $Q (x)$ to ultimately yield the normalized Markov parameters set forth as

\begin{matrix} \begin{matrix} \begin{matrix} s_{- 1} = - 4.8507 e - 36, \\ s_{0} = 8.828274 e - 34, \\ s_{1} = - 8.03373 e - 32, \\ s_{2} = 4.900575 e - 30, \\ s_{3} = - 2.266315 e - 28, \\ s_{4} = 8.517441 e - 27, \end{matrix} \begin{matrix} s_{5} = - 2.721844 e - 25, \\ s_{6} = 7.636194 e - 24, \\ s_{7} = - 1.926143 e - 22, \\ s_{8} = 4.448884 e - 21, \\ ⋮ \end{matrix} \end{matrix} \end{matrix}

(54)

In the interest of brevity, we will omit several intermediate steps in the procedure again. Therefore, the diagonal entries of the upper triangular Cholesky factor of $- {\bar{H}}_{(13, 13)}$ are computed as

\begin{matrix} \begin{matrix} \begin{matrix} r_{1, 1} = 6.5440 e - 35, \\ r_{2, 2} = 3.4192 e - 33, \\ r_{3, 3} = 7.8567 e - 32, \\ r_{4, 4} = 1.1486 e - 30, \\ r_{5, 5} = 1.2007 e - 29, \\ r_{6, 6} = 9.4682 e - 29, \\ r_{7, 7} = 5.7751 e - 28, \end{matrix} \begin{matrix} r_{8, 8} = 2.7464 e - 27, \\ r_{9, 9} = 1.0124 e - 26, \\ r_{10, 10} = 2.8342 e - 26, \\ r_{11, 11} = 7.3965 e - 26, \\ r_{12, 12} = 1.3160 e - 24, \\ r_{13, 13} = 3.6993 e - 24 . \end{matrix} \end{matrix} \end{matrix}

(55)

As all the $r_{j, j}$ are real, we deem our 13th-order system as stable.

Example 4. A nonlinear boiler-turbine system

The dynamics of a boiler-turbine unit has been modeled as the following system of ODEs (Tan et al., 2005)

{\begin{matrix} {\overset{\cdot}{x}}_{1} = - 0.0018 u_{2} x_{1}^{9 / 8} + 0.9 u_{1} - 0.15 u_{3}, \\ {\overset{\cdot}{x}}_{2} = (0.073 u_{2} - 0.016) x_{1}^{9 / 8} - 0.1 x_{2}, \\ {\overset{\cdot}{x}}_{3} = (141 x_{3} - (1.1 u_{2} - 0.19) x_{1}) / 85, \\ y_{1} = x_{1}, \\ y_{2} = x_{2}, \\ y_{3} = 0.05 (0.13073 x_{3} + 100 a_{cs} + q / 9 - 67.975), \end{matrix}

(56)

where the three state variables namely x₁, x₂, and x₃ denote drum pressure (kg/cm²), electric power output (MW), and fluid density (kg/m³), respectively. The three inputs u₁, u₂, and u₃ are the valve positions for the fuel, the steam, and the feed stream. The system’s output y₃ is the drum water level (m). The parameters a_cs and q are the steam quality and the evaporation rate (kg/s), which are given by

a_{cs} = \frac{(1 - 0.001538 x_{3}) (0.8 x_{1} - 25.6)}{x_{3} (1.0394 - 0.0012304 x_{1})}

(57)

q = (0.854 u_{2} - 0.147) x_{1} + 45.59 u_{1} - 2.514 u_{3} - 2.096

(58)

It is clear that system (56) is MIMO and still our criterion is valid to analyze its stability.

We want to investigate the stability of system (56) at a particular steady-state condition (equilibrium point), namely $x_{1} = 108, x_{2} = 66.65, x_{3} = 428, u_{1} = 0.34, u_{2} = 0.69,$ $u_{3} = 0.433$ .

As the first step, we derive the Jacobian matrix associated with system (56) as

J (x) = [\begin{matrix} - 0.002025 u_{2} x_{1}^{1 / 8} & 0 & 0 \\ 9 / 8 (0.073 u_{2} - 0.016) x_{1}^{1 / 8} & - 0.1 & 0 \\ - (1.1 u_{2} - 0.19) / 85 & 0 & 141 / 85 \end{matrix}]

(59)

which is evaluated at the aforementioned equilibrium point as

J (x_{e}) = [\begin{matrix} - 0.0025 & 0 & 0 \\ 0.0694 & - 0.1 & 0 \\ - 0.0067 & 0 & 1.6588 \end{matrix}]

(60)

Next, it is straightforward to obtain the characteristic polynomial of matrix $J (x_{e})$ by the Leverrier-Takeno method as follows

P (x) = x^{3} - 1.5563 x^{2} - 0.169777 x - 0.0004147

(61)

Subsequently, by a Maclaurin expansion of $Q (1 / x)$ as described before, we determined the Markov parameters as

\begin{matrix} \begin{matrix} s_{- 1} = - 1, \\ s_{0} = - 3.1126, \\ s_{1} = - 4.8441, \\ s_{2} = - 8.0682, \\ s_{3} = - 13.3802, \\ s_{4} = - 22.1955, \end{matrix} \begin{matrix} s_{5} = - 36.8178, \\ s_{6} = - 61.0735, \\ s_{7} = - 101.3087, \\ s_{8} = - 168.0510, \\ ⋮ \end{matrix} \end{matrix}

(62)

Then, the associated infinite Hankel matrix formally reads

H = [\begin{matrix} - 3.1126 & - 4.8441 & - 8.0682 & - 13.3802 & - 22.1955 & - 36.8178 & \cdot \\ - 4.8441 & - 8.0682 & - 13.3802 & - 22.1955 & - 36.8178 & \cdot & \cdot \\ - 8.0682 & - 13.3802 & - 22.1955 & - 36.8178 & \cdot & \cdot & \cdot \\ - 13.3802 & - 22.1955 & - 36.8178 & \cdot & \cdot & \cdot & \cdot \\ - 22.1955 & - 36.8178 & \cdot & \cdot & \cdot & \cdot & \cdot \\ - 36.8178 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{matrix}]

(63)

and its submatrix $H_{(3, 3)}$ is extracted as

H_{(3, 3)} = [\begin{matrix} - 3.1126 & - 4.8441 & - 8.0682 \\ - 4.8441 & - 8.0682 & - 13.3802 \\ - 8.0682 & - 13.3802 & - 22.1955 \end{matrix}]

(64)

As the next step, the lower triangular matrix L is calculated

L = [\begin{matrix} - 1 & 0 & 0 \\ - 3.1126 & 1 & 0 \\ - 4.8441 & 3.1126 & - 1 \end{matrix}]

(65)

and immediately it follows that

{\bar{H}}_{(3, 3)} = L H_{(3, 3)} = [\begin{matrix} 3.1126 & 4.8441 & 8.0682 \\ 4.8441 & 7.0097 & 11.7328 \\ 8.0682 & 11.7328 & 19.6316 \end{matrix}]

(66)

Finally, we perform the Cholesky decomposition on the matrix $- {\bar{H}}_{(3, 3)}$ and obtain the diagonal components of the upper triangular Cholesky factor R as

\begin{matrix} r_{1, 1} = 1.7643 i, \\ r_{2, 2} = 0.7275, \\ r_{3, 3} = 0.0119 . \end{matrix}

(67)

Since at least one diagonal component, that is $r_{1, 1}$ , is non-real, we deem system (56) as unstable at the mentioned equilibrium point.

Alternatively, we will test our second criterion, namely Corollary 3.2 for system (56) in the sequel. From equation (61) it is easy to identify

M (x) = - 1.5563 x^{2} - 0.0004147

(68)

N (x) = x^{3} - 0.169777 x

(69)

Since the order of the characteristic polynomial, see equation (61), is odd, we understand that we must expand $M (x) / N (x)$ through a series of polynomial long divisions. In view of equation (12), we compute the first five Markov parameters for $M (x) / N (x)$ as follows

\begin{matrix} δ_{1} = - 1.5563, \\ δ_{2} = 0, \\ δ_{3} = 0.2646, \\ δ_{4} = 0, \\ δ_{5} = - 0.0449, \end{matrix}

(70)

and from them, we build the Hankel matrix $\hat{H}$

\hat{H} = [\begin{matrix} - 1.5563 & 0 & 0.2646 \\ 0 & 0.2646 & 0 \\ 0.2646 & 0 & - 0.0449 \end{matrix}]

(71)

Afterward, we run a Cholesky decomposition to check the positive definiteness of matrix $\hat{H}$ . The upper triangular Cholesky factor matrix R has diagonal entries

\begin{matrix} r_{1, 1} = 1.2475 i, \\ r_{2, 2} = 0.5144, \\ r_{3, 3} = 0.0083 . \end{matrix}

(72)

The existence of a non-real entry in R implies that matrix $\hat{H}$ cannot be Cholesky-decomposed. Hence, $\hat{H}$ is not positive definite. Consequently, the stated equilibrium point of system (56) is not stable.

Computational performance analysis

In this section, we will compare the computational performance of our first stability criterion, namely Corollary 3.1, against that of the Routh–Hurwitz method through a CPU-time benchmarking experiment on a set of randomly generated LTI dynamic systems with orders ranging from 2 to 180. The results of this comparative investigation are depicted in Figure 4 (observe the log-scale). Based on the average CPU-time values, we can conclude that our stability analysis algorithm is approximately 12.8 times faster than the Routh–Hurwitz method. This superiority can be attributed to the fact that the most time-consuming step in our approach is the Maclaurin expansion. The rest are simple substitutions and finally, a single matrix multiplication exists. In contrast, the Routh–Hurwitz algorithm entails a large number of arithmetic operations including multiplication, division, and subtraction for the completion of the Routh table prior to reaching a conclusion on the stability of the system under study. As a side note, we may even increase the computational speed of our method for some cases by exploiting the fact that on the appearance of a non-real diagonal entry in the Cholesky pseudo code, see Definition 2.4, we can prematurely terminate the algorithm in favor of computational efficiency and conclude unstability.

Figure 4.

Computational efficiency analysis of our stability criterion and the Routh–Hurwitz method.

Moreover, in contrast to the Routh–Hurwitz method, which has certain special cases that require special handling and can significantly impact computational performance, our proposed stability criterion is free of such exceptions. These special cases in the Routh–Hurwitz method, such as zero rows or columns in the Routh table, can lead to difficulties in determining the stability of a system and may require additional calculations or adjustments to account for their presence, resulting in increased computational complexity and time.

Our criterion, however, does not suffer from any such issues, as it is applicable to a wide range of linear control systems without the need for any special treatment or exceptions, and therefore offers a more reliable and efficient approach for stability analysis, particularly in terms of computational performance, when compared to the Routh–Hurwitz method and its special cases.

Conclusion

We developed a stability criterion for LTI systems based on Markov parameters of a rational function that is defined from the characteristic polynomial of the system under study. For state-space represented dynamic systems, we used the Leverrier-Takeno algorithm to derive the coefficients of the characteristic polynomial as the initial step. Afterward, the Markov parameters were systematically calculated by a particular Maclaurin series expansion. The relevant Hankel matrices were then generated. The positive definiteness of the appropriate Hankel matrix was verified by the Cholesky decomposition algorithm as the final step of our criterion.

Furthermore, a second criterion was also developed that extracted a different set of Markov parameters of a distinct rational function from the characteristic polynomial by iterated polynomial long divisions. The latter was concluded to be computationally inferior as it entails repeated time-consuming polynomial long divisions. Hence, it was suggested suitable only for low-ordered systems.

Our methods were extended for the local stability analysis of nonlinear continuous-time systems by the concept of the Jacobian matrix and the indirect Lyapunov method of stability. The application of the proposed stability criteria was demonstrated through four illustrative examples including two distillation columns and a nonlinear boiler-turbine combination. In the end, based on the conducted CPU-time analysis for dynamic systems with orders ranging from 2 to 180, it was revealed that our first stability criterion is about 12 times faster than the classical Routh–Hurwitz method. Hence, it can possibly be implemented in fast or real-time stability analysis schemes.

As highlighted within the “Introduction,” our stability criterion is constrained by its inability to accommodate systems with time delays. Nonetheless, it is worth noting that this limitation is a common characteristic shared by the majority of previously established methods, as demonstrated in Table 1.

We intend to improve the computational efficiency of our current method by employing a more advanced positive definiteness checking algorithm as developed by Rump (2006), in a future research work.

Footnotes

Acknowledgements

The authors thank the editor and reviewers of the Transactions of the Institute of Measurement and Control for their invaluable comments and suggestions, which have significantly enhanced the quality of this paper. In addition, the first author (H.F.) thanks Raybod Fatoorehchi for diligently proofreading the paper drafts.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Hooman Fatoorehchi

References

Abolpour

Khayatian

Dehghani

, et al. (2023) An equivalent condition for stability analysis of LTI systems with bounded time-invariant delay. Applied Mathematics and Computation 438(1): 127585.

Adamy

(2022) Nonlinear Systems and Controls. Berlin: Springer.

Cheney

Kincaid

(2009) Linear Algebra: Theory and Applications. London: Jones & Bartlett Learning.

Cogan

De Paor

Quinn

(2009) PI control of first-order lag plus time-delay plants: Root locus design for optimal stability. Transactions of the Institute of Measurement and Control 31(5): 365–379.

Cohen

(1994) Is the polynomial so perfidious? Numerische Mathematik 68: 225–238.

Corriou

J-P

(2013) Process Control: Theory and Applications. London: Springer.

Datta

(1979) Application of Hankel matrices of Markov parameters to the solutions of the Routh-Hurwitz and the Schur-Cohn problems. Journal of Mathematical Analysis and Applications 68(1): 276–290.

Drgoňa

Klaučo

Janeček

, et al. (2017) Optimal control of a laboratory binary distillation column via regionless explicit MPC. Computers & Chemical Engineering 96: 139–148.

Fatoorehchi

Djilali

(2023) Stability analysis of linear time-invariant dynamic systems using the matrix sign function and the Adomian decomposition method. International Journal of Dynamics and Control 11(2): 593–604.

10.

Fatoorehchi

Ehrhardt

(2023) A combined method for stability analysis of linear time invariant control systems based on Hermite-Fujiwara matrix and Cholesky decomposition. Canadian Journal of Chemical Engineering 101: 7043–7052.

11.

Ghiasy

Boodhoo

KVK

Tham

(2012) Control of intensified equipment: A simulation study for pH control in a spinning disc reactor. Chemical Engineering and Processing: Process Intensification 55: 1–7.

12.

Katbab

Jury

(1990) On the stability of a general polynomial in the space of Markov parameters. Systems & Control Letters 14(5): 419–430.

13.

Kavitha

Ramakalyan

(2014) Simple and straight proofs of stability criteria for finite-dimensional linear time invariant systems. Transactions of the Institute of Measurement and Control 36(4): 523–528.

14.

Rasmussen

Jørgensen

(1999) Parametric uncertainty modeling for robust control: A link to identification. Computers & Chemical Engineering 23(8): 987–1003.

15.

Roy

Banerjee

(2014) Linear Algebra and Matrix Analysis for Statistics. London: Taylor & Francis.

16.

Rump

(2006) Verification of positive definiteness. BIT Numerical Mathematics 46: 433–452.

17.

Shinskey

(1996) Process Control Systems: Application, Design, and Tuning. New York: McGraw Hill.

18.

Sun

Kang

(2016) Power system observability and dynamic state estimation for stability monitoring using synchrophasor measurements. Control Engineering Practice 53: 160–172.

19.

Tan

Marquez

Chen

, et al. (2005) Analysis and control of a nonlinear boiler-turbine unit. Journal of Process Control 15(8): 883–891.

20.

Wayland

(1945) Expansion of determinantal equations into polynomial form. Quarterly of Applied Mathematics 2(4): 277–306.

21.

Yakhnin

Menzinger

(2004) High-temperature transients in catalytic fixed-bed reactors. Reviews in Chemical Engineering 20(3–4): 175–224.

22.

Yang

Hou

(2022) A generalization of Routh–Hurwitz stability criterion for fractional-order systems with order α∈ (1, 2). Fractal and Fractional 6(10): 557.