Spectral Decompositions of Gramians of Continuous Stationary Systems Given by Equations of State in Canonical Forms

Abstract

The application of transformations of the state equations of continuous linear and bilinear systems to the canonical form of controllability allows one to simplify the computation of Gramians of these systems. In this paper, we develop the method and obtain algorithms for computation of the controllability and observability Gramians of continuous linear and bilinear stationary systems with many inputs and one output, based on the method of spectral expansion of the Gramians and the iterative method for computing the bilinear systems Gramians. An important feature of the concept is the idea of separability of the Gramians expansion: separate computation of the scalar and matrix parts of the spectral Gramian expansion reduces the sub-Gramian matrices computation to calculation of numerical sequences of their elements. For the continuous linear systems with one output the method and the algorithm of the spectral decomposition of the controllability Gramian are developed in the form of Xiao matrices. Analytical expressions for the diagonal elements of the Gramian matrices are obtained, and by making use of which the rest of the elements can be calculated. For continuous linear systems with many outputs the spectral decompositions of the Gramians in the form of generalized Xiao matrices are obtained, which allows us to significantly reduce the number of calculations. The obtained results are generalized for continuous bilinear systems with one output. Iterative spectral algorithms for computation of elements of Gramians of these systems are proposed. Examples are given that illustrate theoretical results.

1. Introduction

The matrix continuous differential and algebraic Lyapunov and Sylvester equations play an important role in a modern control theory [1,2,3,4,5,6,7]. The first spectral Gramian expansions for linear continuous systems with a simple spectrum were derived from the spectral expansions of the Lyapunov integral representation of the solution of Lyapunov or Sylvester equations [1,2]. In [3], analytical solutions of discrete and continuous Lyapunov equations based on the transformation of the dynamics matrix to Jordan form were obtained. Modern electric power generation and distribution systems are changing rapidly due to the need for the reduction of greenhouse gas emissions, the proliferation of renewable energy sources, the emergence of new storage systems, and the active entry of the consumers into the energy market [8]. New mathematical modelling techniques, such as generation and load forecasting, optimal control of energy storage systems, and new methods for energy systems control will play a key role in 21st century power engineering and in realizing the Internet of Energy concept [9]. The development of new innovative technologies requires the involvement of a modern mathematical control theory, in which the study of the structural properties and computational methods of energy systems Gramians is an important task [10,11].

The Gramian theory of linear dynamical systems is closely related to the problem of calculating the norms of transfer functions and developing the models of systems with reduced dimensionality, which makes it possible to significantly reduce the dimensionality of the approximating model. Among these methods we would like to mention the balanced cutoff, singular decomposition, Krylov subspace method, methods of synthesis of simplified model optimal by criterion H₂-norm of Gramian, and also hybrid methods [1,4,12]. Important promising results have been obtained in the field of computation of Gramians for systems whose models were represented in canonical controllability and observability forms. In [13], the methods for computing Gramians of linear systems given by equations in controllability and observability forms, based on the use of periodic structure matrices, were first proposed. It is shown that using canonical forms of controllability and observability in linear systems allows one to identify and study the periodic structure of the Gramian matrices not only for their computation but also for stability analysis of the systems themselves based on algebraic Routh criteria. Note that for linear systems the transformation of the equations of state into canonical forms of controllability and observability allowed us to develop new methods for computing Gramians in the form of Xiao matrices [6] based on the use of Routh tables [13,14]. In [5,14], the new approach was developed in the direction of using the properties of the impulse transient function and the zero-plaid structure of the controllability Gramian. Bilinear systems, due to their linearity in state and control, are the closest class with respect to linear systems, so their studies paved the way for the study of complex nonlinear systems, primarily systems with smooth nonlinearities. Studies of these systems provide the key to solving many unsolved problems of linear dynamical systems, such as systems with variable parameters [15]. Significant scientific advances in model approximation, monitoring and control of bilinear systems make it possible to extend methods of linear control theory to areas where this is possible. Research in the field of bilinear control systems are closely related to the problem of lowering the order of the model by constructing an approximating model of lower dimension. Solutions of generalized Lyapunov equations using Kronecker products and vectorization method are obtained; however, this method leads to a sharp increase in dimensionality at each iteration step. “The curse of dimensionality” requires the use of matrix elements aggregation methods in the computation process [12,16,17,18,19,20,21]. The structural properties of bilinear systems provide the key to solving of their control problems [20,21,22,23]. Bilinear systems in mathematical control theory can be considered within linear spaces, called Lie algebras, on which bilinear operations are defined in addition to the usual operations [19]. General controllability conditions can be defined for such systems. In [24], it is shown that these conditions can be obtained using the metrics given by the controllability Gramians. An alternative approach for the Gramian spectral decomposition is to apply the apparatus of linear matrix inequalities and semi-definite programming [25].

Main Contribution

The formulation in the framework of a united concept of the problem in Section 2 considers several topics of computation of controllability and observability Gramians. Its important feature is the idea of the separability of Gramian decompositions: separate computation of scalar and matrix parts of a spectral decomposition which allows the reduction of the computation of matrices of sub-Gramians to the computation of numerical sequences of their elements. The use of canonical forms of controllability allow to propose a pioneering approach to computation of Gramians based on the use of Routh tables and Xiao matrices. This paper proposes to improve this approach by using spectral decompositions of the Gramians and the representation of the resolvent of the dynamics matrix by extending the scope to multivariable linear and bilinear control systems [11,26,27,28,29].

In Section 3, we propose to use the representation of the dynamics matrix resolvent of continuous linear stationary systems with many inputs and one output (MISO LTI) in the form of a Faddeev–Levereux series segment for Gramian decomposition [27,30]. Conversion of the state equations to the canonical form of controllability allowed us to exclude the right part of the Lyapunov equations and the Fadeev matrix from the spectral expansions of the Gramians. This allowed us to further simplify the scalar part of the spectral expansions and to link the localization of the Gramian elements with the residues of the scalar transfer function of the linear system. As a result, we obtained new Gramian expansion, where the matrix part of the expansion is the product of unit vectors.

Section 4 presents the main theoretical results of the paper. For continuous MISO LTI systems, the method and algorithm for spectral decomposition of the controllability Gramian in the form of Xiao matrices have been developed. It allowed us to reduce, by almost four times, the number of calculations of the elements of the Gramian matrices and to get rid of the calculations of the Faddeev matrices, which are sensitive to rounding errors. The identities for real stable polynomials of finite order are proved, which allowed us to generalize the application of separable spectral decompositions for Gramians of linear plants with many inputs and many outputs in the form of a generalized Xiao matrix. The analytical expressions for the diagonal elements of the Gramian matrices, which can be used to compute the remaining elements, were obtained. For the continuous linear stationary systems with many inputs and many outputs (MIMO LTI) systems, the spectral expansions of the Gramians in the form of the generalized Xiao matrices were obtained. They can significantly reduce the number of calculations. The obtained results are generalized for the continuous bilinear stationary systems with many inputs and one output (MISO BTI). Iterative spectral algorithms for computing the Gramian elements of these systems are proposed. Sufficient conditions for convergence of iterative algorithms are established. Illustrative examples are given at the end of Section 4.

2. Problem Statement

Consider a linear MISO LTI continuous stationary dynamical system of the form

$$\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right),\hfill \\ y\left(t\right)=Cx\left(t\right)\hfill \end{array}$$

(1)

where x ∈ $R^n$, y ∈ $R^1$, u ∈ $R^m$, x—state vector, y—output vector, u—control vector; A, B, C are real matrices of appropriate dimensions. We will consider real matrices of appropriate dimensions A, B, C. Assume that system (1) is stable, fully controllable, and observable, and that all eigen numbers of matrix A are distinct. Consider continuous algebraic Lyapunov equations, related to Equation (1), of the following form:

$$A{P}^c+P^c{A}^T=-B{B}^T,$$

Consider a bilinear MISO LTI continuous stationary dynamical system of the form

$$\begin{gathered} \dot{x}\left(t\right)=Ax\left(t\right)+{\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{N}_{\gamma }x\left(t\right)u_{\gamma }\left(t\right)+Bu\left(t\right), \\ y\left(t\right)=Cx\left(t\right), \end{gathered}$$

(2)

where $x\in {ℝ}^n, y\in {ℝ}^1, u\in {ℝ}^m, A, B, C, N_{\gamma }$ are real matrices of appropriate dimensions. The linear system (1) is the linear part of the bilinear system (2) If we use a nondegenerate transformation of the variables with matrix ${\displaystyle R_c^F}$. We can consider the MISO LTI system defined by the equation in canonical form of controllability. The matrix B for the MISO system can be represented as

$$B=\left[\begin{array}{ccc}{b}_1& \dots & \begin{array}{ccc}{b}_{\gamma }& \dots & b_m\end{array}\end{array}\right]$$

For dynamical systems of the form (1) and (2) consider ordinary and generalized Lyapunov equations (GLE) of the form

$${\displaystyle A_c^F}{\displaystyle P_c^{Fln}}+{\displaystyle P_c^{Fln}}{\left({\displaystyle A_c^F}\right)}^T=-{\displaystyle b_{\gamma }^F}{\left({\displaystyle b_{\gamma }^F}\right)}^T,$$

(3)

$${\displaystyle A_c^F}{\displaystyle P_{c\gamma }^{Fbln}}+{\displaystyle P_{c\gamma }^{Fbln}}{\left({\displaystyle A_c^F}\right)}^T+{\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{N}_{\gamma }{\displaystyle P_{c\gamma }^{Fbln}}{\displaystyle N_{\gamma }^T}=-{\displaystyle b_{\gamma }^F}{\left({\displaystyle b_{\gamma }^F}\right)}^T,$$

(4)

Let us further consider the channel “γ” MISO LTI of the linear system in the canonical form of controllability

$$x={\displaystyle R_c^F}{x}_c.$$

$$\begin{gathered} {\dot{x}}_c\left(t\right)={\displaystyle A_c^F}{x}_c\left(t\right)+{\displaystyle b_{\gamma }^{Fc}}{u}_{\gamma }\left(t\right), \\ {\displaystyle y_c^F}\left(t\right)={\displaystyle c_{\gamma }^{Fc}}{x}_c\left(t\right), k=0,1,2\dots \end{gathered}$$

(5)

$$\begin{gathered} {\displaystyle A_c^F}=\left[\begin{array}{ccccc}0& 1& 0& & 0\\ 0& 0& 1& & 0\\ 0& 0& & 1& 0\\ 0& 0& 0& 0& 1\\ -a_0& -a_1& -a_2& \dots & -a_{n-1}\end{array}\right], {\displaystyle b_{\gamma }^{Fc}}={\left[\begin{array}{ccc}0& 0& \begin{array}{ccc}..& 0& 1\end{array}\end{array}\right]}^T, \\ a=\left[\begin{array}{ccc}-a_0& -a_1& \begin{array}{ccc}..& -a_{n-2}& -a_{n-1}\end{array}\end{array}\right],{\displaystyle c_{\gamma }^{Fc}}=\left[\begin{array}{ccc}{\xi }_0& {\xi }_1& \begin{array}{ccc}..& {\xi }_{n-2}& {\xi }_{n-1}\end{array}\end{array}\right]. \\ {\dot{x}}_c\left(t\right)={\displaystyle A_c^F}{x}_c\left(t\right)+{\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{N}_{\gamma }{x}_c\left(t\right)u_{\gamma }\left(t\right)+{\displaystyle b_{\gamma }^{Fc}}u\left(t\right), x_c\left(0\right)=0, \\ {\displaystyle y_c^F}\left(t\right)={\displaystyle c_{\gamma }^{Fc}}{x}_c\left(t\right),k=0,1,2\dots \end{gathered}$$

(6)

The following relations are valid [2].

$$\begin{gathered} R_{c\gamma }^F=\left[\begin{array}{ccccc}{b}_{\gamma }& {\displaystyle A_c^F}{b}_{\gamma }& ..& \dots & {\left({\displaystyle A_c^F}\right)}^{n-1}{b}_{\gamma }\end{array}\right]\left[\begin{array}{ccccc}{a}_1& a_2& & a_{n-1}& 1\\ a_2& a_3& a_{n-1}& 1& \\ & a_{n-1}& & & 0\\ a_{n-1}& 1& & 0& 0\\ 1& 0& 0& 0& 0\end{array}\right] \\ {\left({\displaystyle R_{c\gamma }^F}\right)}^{-1}A{\displaystyle R_{c\gamma }^F}=A^F,{\left({\displaystyle R_{c\gamma }^F}\right)}^{-1}b={\displaystyle b_{\gamma }^{Fc}}, c{\displaystyle R_{c\gamma }^F}={\displaystyle c_{\gamma }^F}, {\displaystyle N_{c\gamma }^F}={\left({\displaystyle R_c^F}\right)}^{-1}{N}_{\gamma } \end{gathered}$$

For Gramians we have

$${\displaystyle P_c^{Fln}}={\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{\displaystyle P_{c\gamma }^{Fln}},$$

$${\displaystyle P_{c\gamma }^{bln}}={\displaystyle R_{c\gamma }^F}{\displaystyle P_{c\gamma }^{Fbln}}{\left({\displaystyle R_{c\gamma }^F}\right)}^T,$$

where the matrix ${\displaystyle P_{c\gamma }^F}$ is a solution of the generalized Lyapunov Equation (4).

Consider the Lyapunov equations of the MIMO LTI system

$$A{P}^c+P^c{A}^T=-B{B}^T,$$

and MISO LTI subsystems formed by splitting the right hand side matrix into products of vectors

$$b_{\gamma }{\displaystyle b_{\gamma }^T}$$

$$B{B}^T={\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{b}_{\gamma }{\displaystyle b_{\gamma }^T}, A{\displaystyle P_{\gamma }^c}+{\displaystyle P_{\gamma }^c}{A}^T=-b_{\gamma }{\displaystyle b_{\gamma }^T},$$

Hence the equality follows

$$P^c={\displaystyle {{\displaystyle \sum }}_{\gamma =1}^m}{\displaystyle P_{\gamma }^c}$$

The use of canonical forms of controllability and observability to solve the Lyapunov and Sylvester equations has been proposed in [5,13,14]. This approach is based on reducing the Lyapunov equation to a solution of a linear algebraic equation in the form of Kronecker products and then applying a vectorization method to solve this equation. However, by applying the equations of state in canonical forms of controllability, the matrix $B{B}^T$ of the right-hand side of the Lyapunov equation turns into a near-zero matrix with a single nonzero element, which allowed us to radically change the method itself for solving the linear algebraic equation. As a result, the controllability Gramian takes the form of the Xiao matrix (for matrices of even size) [13]

$$P^{cF}=\left[\begin{array}{ccccccc}{x}_1& 0& -x_2& 0& x_3& & 0\\ 0& x_2& 0& -x_3& & & 0\\ -x_2& 0& x_3& & & & x_{n-2}\\ 0& -x_3& & & & -x_{n-2}& 0\\ x_3& & & & x_{n-2}& 0& -x_{n-1}\\ & 0& & -x_{n-2}& 0& x_{n-1}& 0\\ 0& & x_{n-2}& 0& -x_{n-1}& 0& x_n\end{array}\right]$$

(7)

The Gramian structure has also been called the zero-plaid structure [16]. The main advantage of the new approach is the radical reduction of calculations in the transition to the canonical form of controllability. Instead of calculating $n^2$ elements of the Gramian matrix, it is enough to calculate only n elements. The price for this advantage is the complexity of calculating the elements themselves, which are computed through the elements of the Routh table.

$$\{\begin{array}{c}{x}_n=\frac{1}{2R_{n1}},\\ x_{n-k}=\frac{-{\displaystyle {{\displaystyle \sum }}_{i=1}^{m-1}}{\left(-1\right)}^i{R}_{n-k,i+1}{x}_{n-k+i}, }{R_{n-k,1}}\end{array} k=1,\overline{n-1},$$

(8)

where $R_{ij }$—is an element of Routh table for the system, standing at the intersection of row “i” and column “j” (details can be found at https://openaccess.city.ac.uk/id/eprint/19115/ (accessed on 1 June 2022)).

The aim of the paper is to develop an alternative approach to computation of the controllability and observability Gramians for linear and bilinear systems based on spectral properties of the Faddeev series expansion of the dynamics matrix and study properties of these expansions.

3. Spectral Decompositions of the Controllability and Observability Gramians of the Linear System Represented by the State Equations in Canonical Forms of Controllability and Observability

Consider the spectral decomposition of the resolvent of the dynamics matrix ${\displaystyle A_j^F}$ in the form

$${\left(Is-A^F\right)}^{-1} ={\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle A_j^F}{s}^j}{N\left(s\right)} ,$$

(9)

where $N\left(z\right)$-characteristic polynomial, ${\displaystyle A_j^F}$—Faddeev matrices, $j=1,2,\dots n$ [26,27,29].

Lemma 1.

Consider a linear continuous MISO system of the form (3) represented by equations in the canonical controllability form of the form (5). Consider the expansion of the resolvent of the dynamics matrix A^F in the Fadeev series.

For elements of the last column of matrix${\displaystyle A_j^F}$the following statements are true

$${\left\{{\displaystyle a_{n-k,n-kn}^F}\right\}}^T={\displaystyle e_{n-k}^T}, k=1,2\dots n-1.$$

(10)

Proof. 

Consider the resolvent decomposition of matrix A^F as a segment of Fadeev series [26,29].

Assume $N\left(s\right)=s^n+a_{n-1}{s}^{n-1}+\dots a_1{s}^{}+a_0$, $R_j={\displaystyle A_{j-1}^F}, j=1,2,\dots n;$ Consider Fadeev–Levereux algorithm for system (1)

$${\left(Is-A\right)}^{-1}=\frac{{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}{R}_j{s}^j}{N\left(s\right)}$$

$$a_n=1, R_n=I,$$

$$a_{n-k}=-\frac{1}{k}tr(A{R}_{n-k+1}), k=1,2\dots n.$$

Let us apply the method of mathematical induction for the system (5). The iterative algorithm for calculating the Fadeev matrices and the coefficients of the characteristic equation is as follows

Step one: $a_{n-1}=1, R_n={\displaystyle A_{cn-1}^F}=I,$

Step “k”: Consider forming of the last column of matrices ${\displaystyle A_{n-k}^F}$.

First step: ${\displaystyle A_{cn-\left(k-1\right)}^F}={\displaystyle A_{cn-1}^F}=I.$ $\left\{{\displaystyle a_{n-1,n}^F}\right\}={\left[\begin{array}{ccc}0& 0& \begin{array}{ccc}\dots .& 0& 1\end{array}\end{array}\right]}^T$

Second step:

$${\displaystyle A_{cn-2}^F}=\left[\begin{array}{ccccc}0& 1& 0& & 0\\ 0& 0& 1& & 0\\ 0& 0& & 1& 0\\ 0& 0& 0& 0& 1\\ -a_0& -a_1& \dots & -a_{n-2}& 0\end{array}\right], \left\{{\displaystyle a_{cn-2,2n}^F}\right\}={\left[\begin{array}{ccc}0& \dots & \begin{array}{ccc}0& 1& 0\end{array}\end{array}\right]}^T$$

Propose that for «k−1» step the last column of the matrix ${\displaystyle A_{n-k}^F}$ has the form

$$\left\{{\displaystyle a_{cn-\left(k-1\right),k-1n}^F}\right\}={\left[\begin{array}{ccccc} \underset{1}{\underbrace{0}}& \dots & \underset{n-\left(k-2\right)}{\underbrace{0}}& \underset{n-\left(k-1\right)}{\underbrace{1..}}& \underset{n}{\underbrace{0}}\end{array}\right]}^T$$

Let us introduce the notation. ${\displaystyle A_c^F}{\displaystyle A_{cn-\left(k-1\right)}^F}=S, S=\left[\begin{array}{ccc}\left\{s_1\right\}& \left\{s_2\right\}& \begin{array}{cc}\dots & \left\{s_n\right\}\end{array}\end{array}\right]$.

The last column of the matrix has the form

$$\left\{s_n\right\}={\left[\begin{array}{ccccc}\underset{1}{\underbrace{0}}& \dots & \underset{n-k}{\underbrace{1}}& \underset{n-\left(k-1\right)}{\underbrace{0.}}& \underset{n}{\underbrace{-{\displaystyle a_{cn-\left(k-1\right)}^F}}}\end{array}\right]}^T.$$

According to the Fadeev–Levereux algorithm, we have

$$\left\{{\displaystyle a_{cn-k,n-kn}^F}\right\}={\left[\begin{array}{ccccc} \underset{1}{\underbrace{0}}& \dots & \underset{n-k}{\underbrace{1}}& \underset{n-\left(k-1\right)}{\underbrace{0}}& \underset{n}{\underbrace{0}}\end{array}\right]}^T.$$

□

Note. Note, first, that the resolvent expansion in the form of a Fadeev series does not require calculating the eigenvalues of the dynamics matrix ${\displaystyle A_c^F}$. Secondly, the transfer function of channel “γ” of the linear part is defined by formula

$$\left\{{\displaystyle a_{cn-k,n-kn}^F}\right\}={\left[\begin{array}{ccccc} \underset{1}{\underbrace{0}}& \dots & \underset{n-k}{\underbrace{1}}& \underset{n-\left(k-1\right)}{\underbrace{0}}& \underset{n}{\underbrace{0}}\end{array}\right]}^T$$

where it follows that it is determined only by elements of the last column of the matrix ${\displaystyle A_{n-k}^F}.$

Corollary 1.

The general formulas of spectral decompositions for the controllability Gramians of the canonically transformed controllability of the linear system, taking into account the lemma, take a simpler form [7]

$$\begin{gathered} P^{cF}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}{\displaystyle P_{k,\rho }^{cF}}, {\displaystyle P_{k,\rho }^{cF}}={\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}{\displaystyle A_j^F}{\displaystyle b_{\gamma }^F}{\left({\displaystyle b_{\gamma }^F}\right)}^T{\left({\displaystyle A_j^F}\right)}^T \\ {\displaystyle P_{k,\rho }^{cF}}={\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}{\mathbf{1}}_{j+1\eta +1}, \end{gathered}$$

(11)

where the notation $e_{j+1}{\displaystyle e_{j+1}^T}={\mathbf{1}}_{j+1j+1}$ is used. A similar approach can be applied to derive the formula for spectral decompositions for observability Gramians of the canonically transformed observability of the MISO system. In this case, the formulas [13] are valid

$$\begin{gathered} x={\displaystyle R_0^F}{x}_o, \\ {\dot{x}}_o\left(t\right)={\displaystyle A_c^F}{x}_o\left(t\right)+{\displaystyle b_{\gamma }^{Fo}}u\left(t\right), \\ {y}_o\left(t\right)={\displaystyle c_{\gamma }^{Fo}}x\left(t\right),k=0,1,2\dots \\ {\displaystyle A_o^F}=\left[\begin{array}{ccccc}0& 0& 0& 0& -a_0\\ 1& 0& 0& 0& -a_1\\ 0& 1& 0& 0& ⋮\\ 0& 0& 1& 0& -a_{n-2}\\ 0& 0& 0& \dots 1& -a_{n-1}\end{array}\right], {\displaystyle c_{\gamma }^{Fo}}={\left[\begin{array}{ccc}0& 0& \begin{array}{ccc}..& 0& 1\end{array}\end{array}\right]}^T, {\displaystyle N_{o\gamma }^F}={\left({\displaystyle R_o^F}\right)}^{-1}{N}_{\gamma }. \\ {\displaystyle b_{\gamma }^{Fo}}=\left[\begin{array}{ccc}{\eta }_0& {\eta }_1& \begin{array}{ccc}..& {\eta }_{n-2}& {\eta }_{n-1}\end{array}\end{array}\right]. \\ {\displaystyle R_0^F}={\left\{\left[\begin{array}{ccccc}{a}_1& a_2& & a_{n-1}& 1\\ a_2& a_3& a_{n-1}& 1& \\ & a_{n-1}& & & 0\\ a_{n-1}& 1& & 0& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}c\\ \begin{array}{c}cA\\ ⋮\end{array}\\ c{A}^{n-1}\end{array}\right]\right\}}^{-1}. \end{gathered}$$

(12)

We use Formula (7) and consider the formation of expressions${\displaystyle A_{oj}^{FT}}$$·c^{FT}{c}^F{\displaystyle A_{oj}^F}$. According to the duality principle, we obtain expressions

$$P^{oF}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}{\displaystyle P_{k,\rho }^{oF}},{\displaystyle P_{k,\rho }^{cF}}={\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}{\mathbf{1}}_{j+1\eta +1},$$

(13)

As a result, we obtain the following algorithm for the separable computation of the “jη“ elements of the matrix

$${\displaystyle p_{j\eta }^{oF}}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}\frac{-1}{s_{\rho }+s_k}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}$$

By substituting the identity

$${\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}\equiv {\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}\frac{{\displaystyle s_k^j}{\left(-s_k\right)}^{\eta }}{\dot{N}\left(s_k\right)N\left(-s_k\right)}$$

into Formulae (11) and (13), we obtain new formulas for the expansion of the controllability and observability Gramians

$${\displaystyle P_{k,\rho }^{cF}}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\left(-s_k\right)}^{\eta }}{\dot{N}\left(s_k\right)N\left(-s_k\right)}{\mathbf{1}}_{j+1\eta +1},$$

(14)

$${\displaystyle P_{j\eta }^{oF}}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\left(-s_k\right)}^{\eta }}{\dot{N}\left(s_k\right)N\left(-s_k\right)}{\mathbf{1}}_{j+1\eta +1}$$

(15)

The Gramians of the original system are related to the Gramians of the systems transformed to canonical forms of controllability and observability by equations of the form [28]

$${\displaystyle R_c^F}{P}^{cF}{\displaystyle R_c^{FT}}=P^c, {\left({\displaystyle R_o^{FT}}\right)}^{-1}{P}^{oF}{\left({\displaystyle R_o^F}\right)}^{-1}=P^o$$

(16)

4. Separable Spectral Expansion of the Controllability Gramian of the Linear System in Controllability Canonical Form

Let us write down the Lyapunov equation for system (5)

$${\displaystyle A_c^F}P+P{\left({\displaystyle A_c^F}\right)}^T=\left[\begin{array}{ccc}0& 0& \begin{array}{cc}\dots & 0\end{array}\\ 0& 0& \begin{array}{cc}\dots & \end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ \dots \end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right]$$

(17)

$$\left[\begin{array}{ccc}0& 0& \begin{array}{cc}\dots & 0\end{array}\\ 0& 0& \begin{array}{cc}\dots & \end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ \dots \end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right]=\left[\begin{array}{cc}{0}_{\left(n-1\right)\times \left(n-1\right)}& 0_{\left(n-1\right)\times \left(1\right)}\\ 0_{\left(1\right)\times \left(n-1\right)}& 1_{\left(1\right)\times \left(1\right)}\end{array}\right]$$

(18)

For the matrix on the left-hand side, let us introduce the notation

$$\begin{gathered} {\displaystyle A_c^F}P+P{\left({\displaystyle A_c^F}\right)}^T\doteq P^{\&}, \\ {P}^{\&}=\left[\begin{array}{cc}{P}^{\&}{}_{\left(n-1\right)\times \left(n-1\right)}& p^{\&}{}_{\left(n-1\right)\times \left(1\right)}\\ p^{\&}{}_{\left(1\right)\times \left(n-1\right)}& p^{\&}{}_{nn}\end{array}\right]. \end{gathered}$$

(19)

Let us call an upper (lower) diagonal of a square matrix any diagonal leading from the upper-left corner to the lower-right corner and located above (below) the main diagonal. We form the following sets of matrix elements.

The “1st” upper odd diagonal. The set $U^{\left(1\right)odd}$. Elements

$$\begin{gathered} p_{12} p_{23}\dots p_{nn}, \mathrm{if} n=2k \\ {p}_{12} p_{23}\dots p_{nn-1}, \mathrm{if} n=2k-1 \end{gathered}$$

The “2nd” upper odd diagonal. The set $U^{\left(2\right)odd}$. Elements

$$\begin{gathered} p_{14} p_{25}\dots p_{n-1n}, \mathrm{if} n=2k \\ {p}_{14} p_{25}\dots p_{n-2n}, \mathrm{if} n=2k-1 \end{gathered}$$

The “k-th” upper odd diagonal. The set $U^{\left(k\right)odd}$. Elements

$$\begin{gathered} p_{2k}\hspace{1em}\hspace{1em}\mathrm{if} n=2k \\ {p}_{2k-22k-1} \mathrm{if} n=2k-1 \end{gathered}$$

«1st» lower odd diagonal. The set $V^{\left(1\right)odd}$. Elements

$$\begin{gathered} p_{21} p_{32}\dots p_{nn}, \mathrm{if} n=2k \\ {p}_{21} p_{32}\dots p_{nn-1}, \mathrm{if} n=2k-1 \end{gathered}$$

«2nd» lower odd diagonal. The set $V^{\left(2\right)odd}$. Elements

$$\begin{gathered} p_{41} p_{52}\dots p_{n-3n}, \mathrm{if} n=2k \\ {p}_{41} p_{52}\dots p_{n-1n}, \mathrm{if} n=2k-1 \end{gathered}$$

«k-th» lower odd diagonal. The set $V^{\left(k\right)odd}$. Elements

$$\begin{gathered} p_{n1}\mathrm{if} n=2k \\ {p}_{n2}\mathrm{if} n=2k-1. \end{gathered}$$

Similarly form sets of even diagonals, considering the main diagonal as even zero. Form the sets ${\mathsf{\Omega }}_0$ and ${\mathsf{\Omega }}_{\varnothing }$ of elements of the Gramian P as

$${\mathsf{\Omega }}_0=U^{\left(1\right)odd}{{\displaystyle \cup }}^{ }{U}^{\left(2\right)odd}{{\displaystyle \cup }}^{ }\dots {{\displaystyle \cup }}^{ }{U}^{\left(k\right)odd}{{\displaystyle \cup }}^{ }{V}^{\left(1\right)odd}{{\displaystyle \cup }}^{ }{V}^{\left(2\right)odd}{{\displaystyle \cup }}^{ }\dots {{\displaystyle \cup }}^{ }{V}^{\left(k\right)odd}$$

(20)

$${\mathsf{\Omega }}_{\varnothing }=U^{\left(0\right)even}{{\displaystyle \cup }}^{ }{U}^{\left(1\right)even}{{\displaystyle \cup }}^{ }{U}^{\left(2\right)even}\dots U^{\left(k\right)even}{{\displaystyle \cup }}^{ }{V}^{\left(1\right)even}{{\displaystyle \cup }}^{ }{V}^{\left(2\right)even}\dots V^{\left(k\right)even}$$

(21)

From the above, there are two indications that the elements of a Gramian P belong to the sets ${\mathsf{\Omega }}_0$ and ${\mathsf{\Omega }}_{\varnothing }$

$$p_{j\eta }\in {\mathsf{\Omega }}_0 \mathrm{if} \forall j,\eta :j+\eta =2k, k=1,2,\dots n.$$

(22)

$$p_{j\eta }\in {\mathsf{\Omega }}_{\varnothing } \mathrm{if} \forall j,\eta :j+\eta =2k-1, k=1,2,\dots n.$$

(23)

The structure of the matrix (19) and symmetry of the Gramian makes it possible to reduce the solution of Equation (17) to the solution of four simple linear algebraic equations. Let us write down the first system of the matrix equality (19) as

$$p^{\&}{}_{ij}=\{\begin{array}{c}2p_{i+1,j}=0, if i=j, i=\overline{2,n;} j=\overline{1,n-1}\\ p_{i,j+1 }+p_{i+2,j-1 ,}\forall i\ne j, i=\overline{2,n;} j=\overline{1,n-1}\end{array}$$

(24)

From the symmetry of the Gramian it follows that

$$p^{\&}{}_{ij}=p^{\&}{}_{ji}=p_{j,i+1 }+p_{j+2,i-1 ,}\forall i\ne j, i=\overline{2,n;}j=\overline{1,n-1}$$

(25)

$$p^{\&}{}_{1n}=p_{2n}+\left(a·p_1\right)=0, p^{\&}{}_{2n}=p_{3n}+\left(a·p_2\right)=0, \dots p^{\&}{}_{nn}=\left(a·p_n\right)=0,$$

(26)

$$p^{\&}{}_{n1}=p_{2n}+\left(a·p_1\right)=0, p^{\&}{}_{n2}=p_{3n}+\left(a·p_2\right)=0, \dots p^{\&}{}_{nn}=\left(a·p_n\right)=0,$$

(27)

$$p^{\&}{}_{nn }=2\left(a·p_n\right)=-1$$

(28)

Thus, the system of linear algebraic Equations (24)–(28) implements an alternative method to determine elements of the Gramian control matrices in a simpler way than solving the linear algebraic equation in the form of Kronecker products and applying vectorization method (8) [16]. In addition, Formulas (10) and (11) give analytical solution of Lyapunov equation based on Gramian method. If a linear system is stable and the eigen numbers of its dynamics matrix are different, the solution of the Lyapunov equation is singular. Hence, the two solutions considered are the same.

Let us discuss properties of the matrix and solution elements on the basis of analysis of solutions to the system of Equations (24)–(28) and Formulas (10) and (11).

Assertions.

(1)

All elements of upper and lower even diagonals of the Gramian matrix

${\mathsf{\Omega }}_0$ are zero. Indeed. From (12) and Gramian symmetry, this is true for the elements of the sets $U^{\left(1\right)odd}, V^{\left(1\right)odd}$ odd of the first upper and lower diagonals of the matrix. For the elements of other diagonals of the matrix (both zero and nonzero), the following identities follow from Formulas (10) and (11)

$$p_{j\eta }={\left(-1\right)}^{\eta }{p}_{j+1\eta -1}={\left(-1\right)}^{\eta +2}{p}_{j+2\eta -2}={\left(-1\right)}^{\eta -2}{p}_{j-2\eta +2}=\cdots .$$

(29)

(2)

If the element $p_{j\eta }=0,$ it follows from (27) that

$$p_{j+1\eta -1}=p_{j+2\eta -2}=p_{j-2\eta +2}=\cdots =0$$

(30)

If $p_{j\eta }\ne 0,$ first calculate all diagonal elements of the Gramian by Formulas (10) and (11)

$$p_{11}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{1}{\dot{N}\left(s_k\right)N\left(-s_k\right)},$$

(31)

$$p_{22}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{-{\left(s_k\right)}^2}{\dot{N}\left(s_k\right)N\left(-s_k\right)},$$

(32)

$$p_{nn}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\left(-1\right)}^{n-1}{\left(s_k\right)}^{2\left(n-1\right)}}{\dot{N}\left(s_k\right)N\left(-s_k\right)}.$$

(33)

If j$\ne \eta$

$$p_{j\eta }={\left(-1\right)}^{\frac{j-\eta }{2}}{p}_{ll}, j+\eta =2l$$

(34)

We have proved the theorem.

Theorem 1

[13]. Consider a continuous MISO LTI system (3). Suppose the system is stable and all roots of its characteristic equation are different. Then the elements of its controllability Gramian are solutions of the system of Equations (24)–(28).

For even n, the solution matrix looks like

$$P^{cF}=\left[\begin{array}{ccccccc}{p}_{11}& 0& -p_{22}& 0& p_{33}& & 0\\ 0& p_{22}& 0& -p_{33}& & & \\ -p_{22}& 0& p_{33}& & & & p_{n-2n-2}\\ 0& -p_{33}& & & & -p_{n-2n-2}& 0\\ p_{33}& & & & p_{n-2n-2}& 0& -p_{n-1n-1}\\ & 0& & -p_{n-2n-2}& 0& p_{n-1n-1}& 0\\ 0& & p_{n-2n-2}& 0& -p_{n-1n-1}& 0& p_{nn}\end{array}\right]$$

For odd n we have the formula

$$P^{cF}=\left[\begin{array}{ccccccc}{p}_{11}& 0& -p_{22}& 0& p_{33}& & p_{\frac{n+1}{2}\frac{n+1}{2}}\\ 0& p_{22}& 0& -p_{33}& & -p_{\frac{n+1}{2}\frac{n+1}{2}}& \\ -p_{22}& 0& p_{33}& & & & p_{n-2n-2}\\ 0& -p_{33}& & & & -p_{n-2n-2}& 0\\ p_{33}& -p_{\frac{n+1}{2}\frac{n+1}{2}}& & & p_{n-2n-2}& 0& -p_{n-1n-1}\\ & & & -p_{n-2n-2}& 0& p_{n-1n-1}& 0\\ p_{\frac{n+1}{2}\frac{n+1}{2}}& & p_{n-2n-2}& 0& -p_{n-1n-1}& 0& p_{nn}\end{array}\right]$$

The analytical expressions for calculating the leading elements of the Gramian matrices are defined as

$$\begin{gathered} p_{11}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{1}{\dot{N}\left(s_k\right)N\left(-s_k\right)}, \\ {p}_{22}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\left(-1\right)}^1{\left(s_k\right)}^2}{\dot{N}\left(s_k\right)N\left(-s_k\right)}, \\ {p}_{nn}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\left(-1\right)}^{n-1}{\left(s_k\right)}^{2\left(n-1\right)}}{\dot{N}\left(s_k\right)N\left(-s_k\right)}. \end{gathered}$$

(35)

The analytical expressions for calculating the slave elements of the Gramian matrices through the leading elements are defined as

$$p_{j\eta }={\left(-1\right)}^{\frac{j-\eta }{2}}{p}_{ll}, j+\eta =2l$$

(36)

The rest elements of the Gramian matrices elements are zero.

Corollary 2.

For a continuous stationary stable linear dynamical system MISO LTI with a simple spectrum the following identities are true

$${\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}\equiv 0, j=0,1,\dots n-1; \eta =0,1,\dots n-1,$$

(37)

$${\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\displaystyle s_k^j}{\left(-s_k\right)}^{\eta }}{\dot{N}\left(s_k\right)N\left(-s_k\right)}\equiv 0, j=0,1,\dots n-1; \eta =0,1,\dots n-1,$$

(38)

if the conditions are satisfied

(a) 

For even matrices of size n = 2k, the sum of indices j + η belongs to the one of the set ${\Omega }_0$

$$\forall j, \eta : j+\eta \in \{\begin{array}{c}\left[2,3,\dots \dots .2n-1\right],\\ \begin{array}{c}or\left[3,4,\dots 2n-3\right]\\ \dots \dots \dots \dots ..\end{array}\\ or\left[2k+1\right]\end{array}$$

(39)

(b) 

For matrices of odd size n = 2k−1, the sum of indices j + η belongs to one of the sets ${\Omega }_0$

$$\forall j, \eta : j+\eta \in \{\begin{array}{c}\left[2,3,\dots \dots .2n-1\right],\\ \begin{array}{c}or\left[3,4,\dots 2n-3\right]\\ \dots \dots \dots \dots ..\end{array} \\ or\left[2k-1,2k+1\right]\end{array}$$

(40)

The consequence follows from Formulae (35) and (36) of the theorem. It expresses one of the important properties of Xiao matrices: all the elements of odd diagonals are zero. The conditions (22) and (23) formally indicate that the elements of odd diagonals belong to the set ${\Omega }_0$.

At the same time, identities (37) and (38) are valid for stable real polynomials of order “n”, all roots of which are different and situated in the left half-plane of the complex plane. Clearly, the identities are invariant with respect to any nondegenerate transformation of the coordinates of the dynamical system.

Corollary 3.

Consider a stable continuous stationary linear MIMO LTI dynamical system with simple spectrum with many inputs and many outputs

$$\begin{gathered} \dot{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right), \\ y\left(t\right)=Cx\left(t\right), \end{gathered}$$

(41)

where$x\left(t\right)\in R^n, u\left(t\right)\in R^m, y\left(t\right)\in R^m$.

The controllability sub-Gramian${\displaystyle P_{j,\eta }^c}$ of a continuous stationary stable linear MIMO LTI system of the form (1) with a simple spectrum is a matrix of the form [7]

$$P^c={\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle P_{j,\eta }^c}, {\displaystyle P_{j,\eta }^c}={\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\omega \left(s_k,s_{\rho },j,\eta \right)A_{j}B{B}^T{\left(A_{\eta }\right)}^T$$

(42)

$$\omega \left(s_k,s_{\rho },j,\eta \right)=\{\begin{array}{c}0,if index j\eta belongs to odd diagonals of the subGramian\\ {\displaystyle {{\displaystyle \sum }}_{л=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}\frac{-1}{s_{\rho }+s_k}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}if index j\eta belongs to the rest of diagonals of the sub-Gramian\end{array}$$

Proof. 

As we know, the spectral decomposition of controllability Gramian in consequence of conditions has the form [7]

$$P^c={\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}\frac{-1}{s_{\rho }+s_k}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}{A}_{j}B{B}^T{\left(A_{\eta }\right)}^T$$

(43)

Note that identities (37) and (38) are invariant for any nondegenerate coordinate I inser numbettransformation. Substitute the newly introduced scalar function $\omega \left(s_k,s_{\rho },j,\eta \right)$ into this formula and obtain Formula (42). □

Let’s define the bilinear system controllability Gramian by means of the Volterra matrix series of the form [10], namely

$$\begin{array}{l}{P}_1\left(t_1\right)=e^{A{t}_1}B,\hfill \\ P_i\left(t_1,..t_i\right)=e^{A{t}_i}\left[N_1{P}_{i-1}{N}_2{P}_{i-1}\dots N_m{P}_{i-1}\right], i=2,3,4,\dots ,\hfill \\ P={\displaystyle {{\displaystyle \sum }}_{i=1}^{\infty }}{\displaystyle {{\displaystyle \int }}_0^{\infty }}\dots {\displaystyle {{\displaystyle \int }}_0^{\infty }}{P}_i\left(t_1,..t_i\right){\displaystyle P_i^T}\left(t_1,..t_i\right)d{t}_1\dots t_i.\hfill \end{array}$$

(44)

For system (1) it is possible to define GLE as

$$AP+P{A}^T+{\displaystyle {{\displaystyle \sum }}_{j=1}^m}{N}_{j}P{\displaystyle N_j^T}=-B{B}^T,$$

(45)

The Volterra series (44) is a solution to equation (45), in the case where this solution exists. The solution matrix in this case can be called the Gramian controllability of the bilinear system [11,12,13].

Theorem 2

[11]. If the dynamics matrix of the linear part A is stable and the controllability Gramian of the bilinear system is the only solution of the generalized Lyapunov equation, then the solution is the matrix determined by the following iterative procedure.

$$\begin{array}{ll}\hfill & A{P}_1+P_1{A}^T=-B{B}^T,\hfill \\ \hfill & A{P}_i+P_i{A}^T+{\displaystyle {{\displaystyle \sum }}_{j=1}^m}{N}_j{P}_{i-1}{\displaystyle N_j^T}=0, i=2,3,4\dots ..,\hfill \\ \hfill & P=P_1+{\displaystyle {{\displaystyle \sum }}_{i=2}^{\infty }}{P}_i .\hfill \end{array}$$

(46)

If the controllability Gramian of a bilinear system exists, then it is the limit solution resulting from realization of the iterative procedure (46) [16,30].

Let us introduce the notations $N{\displaystyle A_j^F}\doteq {\mathcal{N}}_j,N{\displaystyle A_{\eta }^F}\doteq {\mathcal{N}}_{\eta }, {\mathcal{N}}_{j\gamma }\doteq \left[{\displaystyle n_{j\nu \mu }^{\gamma }}\right], {\mathcal{N}}_{\eta \gamma }\doteq \left[{\displaystyle n_{\eta \nu \mu }^{\gamma }}\right]$.

Let us consider the formation of the right-hand side of the Lyapunov equation at step “k” Let us write the element « $j\eta$ » of the matrix of the right-hand side in the form

$${\displaystyle e_j^T}{\mathcal{N}}_{j\gamma }{\mathbf{1}}_{j\eta }{\displaystyle p_{}^{Fbln\left(k-1\right)j\eta }}{\displaystyle {\mathcal{N}}_{\eta \gamma }^T}{e}_{\eta }={\displaystyle {{\displaystyle \sum }}_{\nu =1}^n}{\displaystyle {{\displaystyle \sum }}_{\mu =1}^n}-{\left(s_{\nu }+s_{\mu }\right)}^{-1} {\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}{p}^{Fbln\left(k-1\right)j\eta }.$$

(47)

Theorem 3

[16]. Consider a MISO continuous bilinear stationary system (multiple-input one output) represented in the canonical controllability of the form (1). Let the linear part of the bilinear system be fully controllable. Let the matrix A^F be Hurwitzian and have a simple spectrum. Let the inequalities be true.

$$(n^2-{\displaystyle {{\displaystyle \sum }}_{k=1}^{\left(n+1\right)/2}}{2}^k)\underset{\nu \mu }{\underbrace{ max }}\left|{\left(s_{\nu }+s_{\mu }\right)}^{-1}\right|\underset{\nu \mu ij\gamma }{\underbrace{max }}\left|{\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}\right|\underset{k,\rho ,j,\eta }{\underbrace{ max }}<1 ,$$

(48)

Then the Gramian of controllability exists and is unique. The elements of the Gramian matrix can be defined using an iterative procedure of the form [30,31]

$$P^F=P^{Fln\left(1\right)}+P^{Fbln},$$

(49)

$$\begin{gathered} P^{Fln\left(1\right)}={\displaystyle {{\displaystyle \sum }}_{k=1}^n}{\displaystyle {{\displaystyle \sum }}_{\rho =1}^n}{\displaystyle {{\displaystyle \sum }}_{\eta =0}^{n-1}}{\displaystyle {{\displaystyle \sum }}_{j=0}^{n-1}}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}{\mathbf{1}}_{j+1\eta +1}, \\ P^{Fbln\left(k\right)j\eta \gamma }={{\displaystyle \sum }}_{\nu ,\mu }{r}^{\left(k\right)j\eta \gamma }{\displaystyle p_{\nu \mu }^{Fbln\left(k-1\right)ij\gamma }}{\mathbf{1}}_{\nu \mu }, r^{\left(k\right)ij\gamma }= \\ =-\left[{\left(s_{\nu }+s_{\mu }\right)}^{-1}{\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}\omega \left(s_k,s_{\rho },j,\eta \right)\right], \\ k=2,3,..\infty . \forall \nu ,\mu ,i,j=1,2,\dots n;\gamma =1,2,\dots H, \end{gathered}$$

(50)

The initial Gramian of controllability$P^{bln}$of the bilinear system is related to the matrix$P^{cFbln}$ by the Equation (16).

Proof of Theorem 2.

According to the theorem, the formation of the solution matrix is reduced to the construction of the sequence of elements «$j\eta$» and subsequent aggregation of elements into a single matrix. All partial sequences in the general case are complex-valued. To prove the convergence of the sequence of partial sums, let us apply the comparison feature and construct a major sequence from the modules of the sequence members. For each step “k” and each matrix $N_{\gamma }$ the iterative relations (46) take place. Let us construct a comparison series for the elements of the subGramian “ijγ” for step “k”. It follows from formulae (46) that each element of the sequence is a weighted sum of all leading elements in the previous step.

Step 1. Consider the formation of the right part of the generalized Lyapunov equation at the first step for the case γ = 1. A separable spectral decomposition of this solution by a pairwise spectrum of the matrix, was obtained above (10) [32].

k-th step. Consider the formation of the right-hand side of the Lyapunov equation at step “k”. Let us write the “jη” element of the right-hand side matrix as

$${\displaystyle e_j^T}{\mathcal{N}}_{j\gamma }{\mathbf{1}}_{j\eta }{\displaystyle p_{}^{Fbln\left(k-1\right)j\eta }}{\displaystyle {\mathcal{N}}_{\eta \gamma }^T}{e}_{\eta }={\displaystyle {\sum }_{\nu =1}^n}{\displaystyle {\sum }_{\mu =1}^n}-{\left(s_{\nu }+s_{\mu }\right)}^{-1} {\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}{p}^{Fbln\left(k-1\right)j\eta }.$$

(51)

Given the summation of sub-Gramians by the index “γ”, we obtain the formula for calculating the kernel matrix of the Gramian of order “k” at step “k”.

$$\begin{gathered} P^{Fbln\left(k\right)j\eta \gamma }={{\displaystyle \sum }}_{\nu ,\mu }{r}^{\left(k\right)j\eta \gamma }{\displaystyle p_{\nu \mu }^{Fbln\left(k-1\right)ij\gamma }}{\mathbf{1}}_{\nu \mu },r^{\left(k\right)ij\gamma }=-\left[{\left(s_{\nu }+ \\ {s}_{\mu }\right)}^{-1}{\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}\omega \left(s_k,s_{\rho },j,\eta \right)\right],k=2,3,..\infty . \forall \nu ,\mu ,i,j=1,2,\dots n;\gamma =1,2,\dots m, \end{gathered}$$

(52)

Formula (52) expresses an algorithm for elemental computation of the Gramian matrix of a bilinear system. At each step, the algorithm makes it possible to compute the kernel matrix of the Gramian of order “k”.

The global major of all sequences converges if the condition

$$(n^2-{\displaystyle {{\displaystyle \sum }}_{k=1}^{\left(n+1\right)/2}}{2}^k)\underset{\nu \mu }{\underbrace{ max }}\left|{\left(s_{\nu }+s_{\mu }\right)}^{-1}\right|\underset{\nu \mu ij\gamma }{\underbrace{max }}\left|{\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}\right||\underset{k,\rho ,j,\eta }{\underbrace{max \omega (s_k,s_{\rho },j,\eta )}}|$$

Since the linear part is stable, there exists an exact upper bound for the inverse of the moduli of the sum of any eigen numbers of its matrix, equal to

$$\underset{\nu \mu }{\underbrace{max}}\left|{\left(s_{\nu }+s_{\mu }\right)}^{-1}\right|$$

For any elements of matrices ${\displaystyle n_{j\nu \mu }^{\gamma }}{\displaystyle n_{\eta \nu \mu }^{\gamma }}$ there exists an exact upper bound of the products of the moduli of its elements. Hence, if the conditions of the theorem are satisfied, inequality (48) and the signs of convergence of series with positive terms are satisfied. Hence, complex-valued sequences (52) converge uniformly and absolutely. The uniform convergence of all sequences is equivalent to the uniform convergence of matrices in (46) [16]. □

Example 1.

Validation of spectral algorithms for computing controllability Gramians of linear continuous stationary systems given by equations of state in the canonical form of controllability.

$$\begin{gathered} \dot{x}=A^{F}x+b^{F}u, \\ y=c^{F}x, \\ n=3, N\left(x\right)=x^3+a_2{x}^2+a_{1}x+a_0. a_2=4,5, a_1=6,5, a_0=3, \\ {s}_1=-1, s_2=-2, s_3=-1.5. \\ {\displaystyle A_{}^F}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -3& -6.5& -4.5\end{array}\right], {\displaystyle b_{}^F}={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T, \\ {\displaystyle c_{\gamma }^F}=\left[\begin{array}{ccc}{\xi }_0& {\xi }_1& {\xi }_2\end{array}\right]. \end{gathered}$$

(53)

Formula for calculating a Gramian elements

$$p_{j\eta }={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\displaystyle s_k^{j-1}}{\left(-s_k\right)}^{\eta -1}}{\dot{N}\left(s_k\right)N\left(-s_k\right)},$$

$$\dot{N}\left(s_1\right)N\left(-s_1\right)=7.5; \dot{N}\left(s_2\right)N\left(-s_2\right)=21; \dot{N}\left(s_3\right)N\left(-s_3\right)=-6.5625.$$

The formula calculation gives the following results

$$p_{j\eta }={\displaystyle {{\displaystyle \sum }}_{k=1}^n}\frac{{\displaystyle s_k^{j-1}}{\left(-s_k\right)}^{\eta -1}}{\dot{N}\left(s_k\right)N\left(-s_k\right)},$$

$$\dot{N}\left(s_1\right)N\left(-s_1\right)=7.5; \dot{N}\left(s_2\right)N\left(-s_2\right)=21; \dot{N}\left(s_3\right)N\left(-s_3\right)=-6.5625.$$

$$p_{21}=p_{12}=\frac{\left(-1\right)1}{7.5}+\frac{\left(-2\right)1}{21}+\frac{\left(-1.5\right)1}{-6.5625}=0, p_{32}=p_{23}=\frac{\left(1\right)1}{7.5}+\frac{{\left(2\right)}^{2}2}{21}+\frac{{\left(1.5\right)}^{2}1.5}{-6.5625}=0,$$

$$p_{11}=\frac{1}{7.5}+\frac{1}{21}+\frac{1}{-6.5625}=0.02855, p_{22}=\frac{\left(-1\right)}{7.5}+\frac{\left(-4\right)1}{21}+\frac{\left(-2.25\right)1}{-6.5625}=0.01906,$$

$$p_{33}=\frac{1}{7.5}+\frac{16}{21}+\frac{5.0625}{-6.5625}=0.1238, p_{31}=\frac{1}{7.5}+\frac{4}{21}+\frac{2.25}{-6.5625}=0.01906=-p_{22},$$

$$p_{13}=\frac{1}{7.5}+\frac{4}{21}+\frac{2.25}{-6.5625}=0.01906=-p_{22}.$$

The controllability Gramian has the form

$$P^F=\left[\begin{array}{ccc}0.02855& 0& 0. 01906\\ 0& -0. 01906& 0\\ 0. 01906& 0& 0.1238\end{array}\right]$$

Substituting the resulting matrix into the Lyapunov equation, we obtain

$$\begin{gathered} \left[\begin{array}{ccc}0& 0.01906& 0\\ -0,01906& 0& 0.1238\\ 0.00003& -0.1238& -0.4999\end{array}\right]+\left[\begin{array}{ccc}0& -0.01906& 0.00003\\ 0,01906& 0& -0.1238\\ 0& 0.1238& -0.4999\end{array}\right] \\ =\left[\begin{array}{ccc}0& 0& -0.00003\\ 0& 0& 0\\ -0.00003& 0& -0.99998\end{array}\right]\cong \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right] \end{gathered}$$

Example 2.

Consider further the generalized Lyapunov equation of the second order in the canonical form of controllability of the form

$$A^F{P}^{Fbln}+P^{Fbln}{\left(A^F\right)}^T+N^F{P}^{Fbln}{(N^F)}^T=-b^F{\left(b^F\right)}^T,$$

(54)

and the Lyapunov equation for its linear part

$$A^F{P}^{Fln}+P^{Fln}{\left(A^F\right)}^T=-b^F{\left(b^F\right)}^T,$$

(55)

Let us also give the corresponding initial Lyapunov equations in eigenvalue decomposition (EVD) canonical form

$$A{P}^{bln}+P^{bln}{\left(A\right)}^T+N{P}^{bln}{N}^T=-b{b}^T,$$

$$A{P}^{ln}+P^{ln}{\left(A\right)}^T=b{b}^T,$$

In the numerical example, let us take

$$A=\left[\begin{array}{cc}-1& 0\\ 0& -2\end{array}\right], b=\left[\begin{array}{c}\sqrt{3}\\ \sqrt{3}\end{array}\right], N=\left[\begin{array}{cc}\sqrt{3}& 1.5\sqrt{3}\\ 0.5\sqrt{3}& \sqrt{3}\end{array}\right]$$

$$A^{Fbln}=\left[\begin{array}{cc}0& 1\\ -2& -3\end{array}\right], b^F=\left[\begin{array}{c}0\\ 1\end{array}\right], N^F=\left[\begin{array}{cc}0.5& 0.5\\ 0& 0.5\end{array}\right]$$

To determine the Gramian of controllability of the linear part we apply the formula

$$\begin{gathered} P^{Fln}={\displaystyle {\displaystyle \sum }_{k=1}^2}{\displaystyle {\displaystyle \sum }_{\rho =1}^2}{\displaystyle {\displaystyle \sum }_{\eta =0}^1}{\displaystyle {\displaystyle \sum }_{j=0}^1}\frac{{\displaystyle s_k^j}{\displaystyle s_{\rho }^{\eta }}}{\dot{N}\left(s_k\right)\dot{N}\left(s_{\rho }\right)}\frac{-1}{s_{\rho }+s_k}{\mathbf{1}}_{j+1\eta +1}, \\ {P}^{Fln}=\left[\begin{array}{cc}{\displaystyle p_{11}^{Fln}}& {\displaystyle p_{12}^{Fln}}\\ {\displaystyle p_{21}^{Fln}}& {\displaystyle p_{22}^{Fln}}\end{array}\right], \end{gathered}$$

(56)

$${\displaystyle p_{11}^{Fln}}=\frac{6+3-4-4}{12}=\frac{1}{12},$$

(57)

$${\displaystyle p_{12}^{Fln}}=\frac{-6+8+4-6}{12}=0,$$

(58)

$${\displaystyle p_{11}^{Fln}}=\frac{-6+8+4-6}{12}=0,$$

(59)

$${\displaystyle p_{11}^{Fln}}=\frac{6-8-8+12}{12}=\frac{2}{12},$$

(60)

The controllability Gramian of the linear part is a Xiao matrix and simultaneously a diagonal matrix with positive diagonal elements. As we will see later, these properties are not inherited by the matrices of the kernels of the spectral expansion of the bilinear Gramian. Note that if we do not look for spectral expansions of solutions of the Lyapunov equations of the form (57)–(60), but look for the solutions themselves, then because of the simple structure of the Frobenius matrices, these solutions in our case can be obtained by solving simple systems of linear algebraic equations (SLAE). Let us introduce the notation

$$P^{Fln}=\left[\begin{array}{cc}{\displaystyle x_1^{ln}}& {\displaystyle x_4^{ln}}\\ {\displaystyle x_3^{ln}}& {\displaystyle x_2^{ln}}\end{array}\right], P^{Fbln}=\left[\begin{array}{cc}{\displaystyle x_1^{bln}}& {\displaystyle x_4^{bln}}\\ {\displaystyle x_3^{bln}}& {\displaystyle x_2^{bln}}\end{array}\right].$$

Given the notations, the Lyapunov equation of the linear part takes the form of element-by-element equality.

$$\left[\begin{array}{cc}{\displaystyle x_3^{ln}}+{\displaystyle x_4^{ln}}& {\displaystyle x_2^{ln}}-2{\displaystyle x_1^{ln}}-3{\displaystyle x_4^{ln}}\\ {\displaystyle x_2^{ln}}-2{\displaystyle x_1^{ln}}-3{\displaystyle x_3^{ln}}& 2{\displaystyle x_4^{ln}}-2{\displaystyle x_3^{ln}}-6{\displaystyle x_2^{ln}}\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& -1\end{array}\right],$$

(61)

These equations can be written as a SLAE system

$$\begin{gathered} {\displaystyle x_3^{ln}}+{\displaystyle x_4^{ln}}=0, \\ {\displaystyle x_3^{ln}}={\displaystyle x_4^{ln}}, \\ {\displaystyle x_2^{ln}}-2{\displaystyle x_1^{ln}}-3{\displaystyle x_4^{ln}}=0, \\ 2{\displaystyle x_4^{ln}}-2{\displaystyle x_3^{ln}}-6{\displaystyle x_2^{ln}}=-1. \end{gathered}$$

(62)

Substituting solutions (62) into system (61), we can see that they are solutions of the Lyapunov equation of the linear part. The same method can be applied to solve the generalized Lyapunov equation, but it is far from obvious. Let us write the corresponding SLAE of the bilinear system as

$$\begin{gathered} -N^F{P}^{Fln}{(N^F)}^T=\left[\begin{array}{cc}-{\displaystyle q_1^{\left(2\right)}}& -{\displaystyle q_4^{\left(2\right)}}\\ -{\displaystyle q_3^{\left(2\right)}}& -{\displaystyle q_2^{\left(2\right)}}\end{array}\right] \\ {\displaystyle x_3^{bln}}+{\displaystyle x_4^{bln}}+0,01{\displaystyle {{\displaystyle \sum }}_{i=1}^4}{\displaystyle x_i^{bln}}=0, {\displaystyle x_3^{ln}}={\displaystyle x_4^{ln}}, \\ 2{\displaystyle x_4^{bln}}-2{\displaystyle x_3^{bln}}-6{\displaystyle x_2^{bln}}+0,01{\displaystyle x_2^{bln}}=-1, \\ {\displaystyle x_2^{bln}}-2{\displaystyle x_1^{bln}}-3{\displaystyle x_3^{bln}}+0,01{\displaystyle {{\displaystyle \sum }}_{i=2}^3}{\displaystyle x_i^{bln}}=0 \end{gathered}$$

(63)

Its solution has the form

$${\displaystyle x_1^{bln}}=-\frac{11}{12}{\displaystyle q_1^{}}-\frac{1}{12}{\displaystyle q_2^{}}-\frac{1}{2}{\displaystyle q_3^{}}, {\displaystyle x_3^{bln}}={\displaystyle x_4^{bln}}=\frac{1}{2}{\displaystyle q_1^{}},$$

(64)

$${\displaystyle x_2^{bln}}=-\frac{1}{6}\left({\displaystyle q_2^{}}+2{\displaystyle q_1^{}}\right)$$

(65)

Let us move on to calculating the bilinear part matrices.

Second step

Let us calculate a matrix for the right part of the Lyapunov equation at the second step

$$-N^F{P}^{Fln}{(N^F)}^T=\left[\begin{array}{cc}-{\displaystyle q_1^{\left(2\right)}}& -{\displaystyle q_4^{\left(2\right)}}\\ -{\displaystyle q_3^{\left(2\right)}}& -{\displaystyle q_2^{\left(2\right)}}\end{array}\right]=\left[\begin{array}{cc}-0.002499& -0.001666\\ -0.001666& -0.001666\end{array}\right]$$

Let us use formulas (64) and (65) to calculate the elements of the solution matrix of this equation

$${\displaystyle x_1^{bln\left(2\right)}}=-\frac{11}{12}{\displaystyle q_1^{\left(2\right)}}-\frac{1}{12}{\displaystyle q_2^{\left(2\right)}}-\frac{1}{2}{\displaystyle q_3^{\left(2\right)}}=0.003677,$$

$${\displaystyle x_2^{bln\left(2\right)}}=-\frac{1}{6}\left({\displaystyle q_2^{\left(2\right)}}+2{\displaystyle q_1^{\left(2\right)}}\right)=0.001110,$$

$${\displaystyle x_3^{bln\left(2\right)}}={\displaystyle x_4^{bln\left(2\right)}}=\frac{1}{2}{\displaystyle q_1^{\left(2\right)}}=0.001249.$$

Third step

Calculate the matrix for the right-hand side of the Lyapunov equation in step three

$$-N^F{P}^{Fbln\left(2\right)}{(N^F)}^T=\left[\begin{array}{cc}-{\displaystyle q_1^{\left(3\right)}}& -{\displaystyle q_4^{\left(3\right)}}\\ -{\displaystyle q_3^{\left(3\right)}}& -{\displaystyle q_2^{\left(3\right)}}\end{array}\right]=\left[\begin{array}{cc}0.0004926& 0.00002359\\ 0.0002359& 0.0001110\end{array}\right]$$

Let us use Formulas (64) and (65) to calculate the elements of the solution matrix of this equation

$${\displaystyle x_1^{bln\left(3\right)}}=-\frac{11}{12}{\displaystyle q_1^{\left(3\right)}}-\frac{1}{12}{\displaystyle q_2^{\left(3\right)}}-\frac{1}{2}{\displaystyle q_3^{\left(3\right)}}=0.0007035,$$

$${\displaystyle x_2^{bln\left(3\right)}}=-\frac{1}{6}\left({\displaystyle q_2^{\left(3\right)}}+2{\displaystyle q_1^{\left(3\right)}}\right)=0.0002463,$$

$${\displaystyle x_3^{bln\left(3\right)}}={\displaystyle x_4^{bln\left(3\right)}}=\frac{1}{2}{\displaystyle q_1^{\left(3\right)}}=0.0001825.$$

Fourth step

Calculate the matrix for the right-hand side of the Lyapunov equation in step three

$$-N^F{P}^{Fbln\left(3\right)}{(N^F)}^T=\left[\begin{array}{cc}-{\displaystyle q_1^{\left(4\right)}}& -{\displaystyle q_4^{\left(4\right)}}\\ -{\displaystyle q_3^{\left(4\right)}}& -{\displaystyle q_2^{\left(4\right)}}\end{array}\right]=\left[\begin{array}{cc}0.00001899& 0.000004288\\ 0.000004288& 0.00001825\end{array}\right]$$

Let us use Formulas (64) and (65) to calculate the elements of the solution matrix of this equation

$${\displaystyle x_1^{bln\left(4\right)}}=-\frac{11}{12}{\displaystyle q_1^{\left(4\right)}}-\frac{1}{12}{\displaystyle q_2^{\left(4\right)}}-\frac{1}{2}{\displaystyle q_3^{\left(4\right)}}=0.00001956,$$

$${\displaystyle x_2^{bln\left(4\right)}}=-\frac{1}{6}\left({\displaystyle q_2^{\left(4\right)}}+2{\displaystyle q_1^{\left(4\right)}}\right)=0.0000931,$$

$${\displaystyle x_3^{bln\left(4\right)}}={\displaystyle x_4^{bln\left(4\right)}}=\frac{1}{2}{\displaystyle q_1^{\left(4\right)}}=0.00000949.$$

Calculate the sum of the bilinear system matrices after four iterations

$$P^{Fbln}=P^{Fln}+P^{Fbln\left(2\right)}+P^{Fbln\left(3\right)}+P^{Fbln\left(4\right)}=\left[\begin{array}{cc}0.08774& 0.001504\\ 0.001504& 0.1679\end{array}\right]$$

Let us check the obtained solution by substituting it into the generalized Lyapunov equation

$$\begin{gathered} \left[\begin{array}{cc}0.001504& 0.1679\\ -0.1799& -0.5006\end{array}\right]+\left[\begin{array}{cc}0.001504& -0.1799\\ 0.1679& -0.5006\end{array}\right]+\left[\begin{array}{cc}0.002586& 0.001694\\ 0.001694& 0.001679\end{array}\right]= \\ \left[\begin{array}{cc}0.005666& 0.010306\\ 0.010306& -0.9995\end{array}\right]\cong \left[\begin{array}{cc}0& 0\\ 0& -1\end{array}\right] \end{gathered}$$

Obviously, the matrix$P^{Fbln}$ > 0.

5. Conclusions

The following results are obtained in the paper:

• A new method and new algorithms for computing the Xiao matrices.
• Identities for the spectral expansions of the stable polynomials of finite order.
• Computation formulas for the Gramians of continuous stable systems with many inputs, given by the equations of state in the canonical forms of controllability and observability.
• Methods and algorithms for computing sequences of elements of Volterra matrices for Gramians of continuous bilinear systems.

These methods allow one to significantly reduce the computational complexity of solving the above problems. A new separable iterative spectral method to compute the elements of controllability Gramians for stable linear and bilinear stationary continuous systems with many inputs and one output given by equations in the canonical form of controllability has been developed. For the linear part of the novel method includes analytical formulas to compute all the elements of the Gramian based on the eigenvalues of the dynamics matrix resolvent and its residues. To compute the Gramian of a bilinear system, a new iterative spectral algorithm is proposed which preserves the structure of the generalized Xiao matrix at each step, which allows the computations to be reduced significantly and the analysis of convergence of iterative algorithms to be simplified [13]. The role of the diagonal elements of the Gramians is very important. They not only determine the remaining nonzero elements of the Gramian matrices of the linear part but also form the scalar multiplier of the spectral decomposition of the bilinear part. The novel method also has the advantage to novel algorithm for calculating the Gramian based on the known method of calculating the Gramian matrix using the Xiao matrix and Routh tables. However, the last method is highly sensitive to rounding errors, so it is not recommended for high dimensional systems.

The emerging areas for future research are the analysis of the stability of energy-saving control systems, development a reduced-order observer in the modal control systems and the findings using geometric theory and Lie groups.