Analysis and Prediction of Electric Power System’s Stability Based on Virtual State Estimators

Abstract

The stability of bilinear systems is investigated using spectral techniques such as selective modal analysis. Predictive models of bilinear systems based on inductive knowledge extracted by big data mining techniques are applied with associative search of statistical patterns. A method and an algorithm for the elementwise solution of the generalized matrix Lyapunov equation are developed for discrete bilinear systems. The method is based on calculating the sequence of values of a fixed element of the solution matrix, which depends on the product of the eigenvalues of the dynamics matrix of the linear part and the elements of the nonlinearity matrixes. A sufficient condition for the convergence of all sequences is obtained, which is also a BIBO (bounded input bounded output) systems stability condition for the bilinear system.

1. Introduction

Stability estimators (soft sensors) are increasingly used in mode planning and supervisory control in present day electric power systems (EPS) [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. Due to the growing popularity of distributed generation and renewable energy sources (RES) in EPS, its operation modes may approach the stability limits. It is, therefore, necessary to improve both predictive modeling techniques and the methods of analysis and preventive control not only for small deviation modes, but also for essentially nonlinear ones [14,20,21,22,23]. Bilinear systems are the closest ones to the essentially nonlinear class among all nonlinear systems. Therefore, research methods for bilinear systems have been actively developed over recent decades [2,3,4,5,20,21,24,25,26,27,28]. Spectral methods of stability analysis, in particular, selective modal analysis, are widely used in the design and operation of EPS [16,18,22,23].

This article presents the results of the development of these methods and their extension over the class of bilinear EPS models. This will expand their application area [7,8,19,21,28]. To create bilinear models of discrete stationary dynamical systems, digital identification methods and associative search algorithms are used, based on the intelligent analysis of big data obtained from system operation monitoring.

The work [13] develops the Poincaré normal form method for analyzing the stability of energy systems based on continuous dynamical systems with smooth nonlinearities. This approach can be considered to be an alternative to stability estimator development. The article [21] uses a virtual model of energy system’s inertia to control the frequency in a system with a high level of microgrid penetration that shows the possibility of using stability estimators not only for stability monitoring tasks, but also for controlling the frequency of low-frequency oscillations.

In [26], Volterra equations are proposed for analyzing the stability of power systems with renewable energy sources. As against [26], Volterra equations are used by the authors for developing digital twins of bilinear models for EPS. The work [27] shows an effective method of Lyapunov stability indices for studying small-signal stability of EPS, which can be used to solve problems discussed in our article.

Separable spectral expansions of discrete Lyapunov equations are obtained for MISO LTI (multiple input single output linear time invariant) discrete dynamical systems governed by state equations in controllability and observability forms. A method and an algorithm for the element-wise solution of the generalized matrix Lyapunov equation are developed for discrete bilinear systems. The method is based on calculating the sequence of values of a fixed element of the solution matrix, which depends on the product of the eigenvalues of the dynamics matrix of the linear part and the elements of the nonlinearity matrixes. The new method is a spectral version of the iterative method used for solving this equation.

The new method changes the very computing paradigm: sequences of elements are now calculated instead of decision matrixes. In addition to this paradigm change, the method changes the approach to studying the stability of the initial nonlinear model of the power system. The convergence of all sequences of elements means BIBO stability of the original system in a wider area than the area of “small-signal stability of power systems”, as well as identify new indicators of stability in this area.

For developing bilinear models of discrete stationary dynamical systems of this class, the authors propose digital identification methods and associative search algorithms, based on the intelligent analysis of big data collected during system operation.

Computing speed is the key advantage of such identification techniques against the known methods of bilinear model development. The off-line training of the identification system is carried out in advance; further, in the course of real-time operation, the current values of model parameters are obtained using the knowledge accumulated at the training phase [29]. This is done by choosing analogs from the appropriate cluster [18]. It should be noted that the associative search method creates a new linear model at each time step. Section 2 shows how such a model can be further used for digital bilinear model development. In Section 6, an example of obtaining the values of the parameters of linear associative models is cited. In future studies, it is planned to develop a version of this method for the non-stationary case.

2. Knowledge-Based Bilinear Models of Discrete Stationary Dynamic Systems

The essence of the machine learning procedure is as follows [9]. For the current time instant k, a set of impacts $u_k \left(u_k\in {ℝ}^m\right)$ on the stationary system during the time interval $T=\left\{k-T, k-T+1,\dots ,k\right\}$, is divided into clusters (together with the corresponding values of the outputs $y_{k-i}, i=0,\dots ,T$). The clustering procedure is carried out with reference to the distance between the vectors. For the current vector $u_k$, a set of vectors $u_{k-i}$ and the corresponding outputs $y_{k-i}$ are collected within the corresponding cluster. Next, a system of linear equations is formed for the unknown coefficients and the output $y_k.$ Unlike traditional regression models, this model does not contain all the prehistory, but rather especially selected vectors (the closest to the current input vector subject to a certain criterion) named “associations”.

The least squares method provides a solution to this system of equations, which is optimal if the conditions of the Gauss-Markov theorem are met [30]. Statistical independence of the model variables is a condition of this theorem, which is not met for closed-loop systems. The transition to a system of simultaneous linear equations can be done in particular per the Moore-Penrose procedure [31,32]. As a result, a pseudo-solution of the original system of equations can be obtained such that the resulting linear model will have accuracy admissible for a wide range of applications.

It should be noted that the described identification algorithm generates point models, the best ones for the nonlinear system under investigation at a time instant. Therefore, unlike traditional identification algorithms, we do not improve a single model ad infinitum, rather we deal with a sequence of digital ad hoc models; each one is the best fit at the specific time instant subject to the chosen criterion.

Another feature of the models obtained by machine learning is the fact that if the corresponding model accuracy requirements are met, then the solution does not need to be found every time, it can be rather “found” in the cluster, which contains the current vector $u_k$. This can be the “nearest neighbor”, or a vector selected in some other way, in particular, the cluster’s centroid.

If it is nevertheless necessary to solve a system of linear equations, it will be possible in the near future using quantum algorithms.

Machine learning procedures are carried out off-line, at the training stage. Therefore, this identification algorithm demonstrates high speed.

There is a class of essentially nonlinear systems, for which the accuracy of point linear models may be insufficient. This class can include nonlinear systems described by the following equations:

$$\begin{array}{c}\dot{x}\left(t\right)=f\left(x,u,t\right)=f\left(x,t\right)+b\left(x\right)u\left(t\right), x\left(0\right)=0\\ y\left(k\right)=Cx\left(k\right),\end{array}$$

(1)

where$x\left(t\right)$is the state vector of the system, $x\left(t\right)\in {ℝ}^n$, $f \mathrm{and} b$ are nonlinear functions. Models of such systems in the form of:

$$\dot{x}\left(t\right)=ax\left(t\right)+Nx\left(t\right)u\left(t\right)+bu\left(t\right)$$

(2)

where $ax\left(t\right)$ + $bu\left(t\right)$ is the system’s linear part, are called bilinear models.

If the matrix $B\in {ℝ}^{n\times m},$ the system can be represented as

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

(3)

or, in the discrete case:

$$\begin{array}{c}x\left(k+1\right)=Ax\left(k\right)+{\displaystyle \sum _{\gamma =1}^m}{N}_{\gamma }x\left(k\right)u_{\gamma }\left(k\right)+Bu\left(k\right), x\left(0\right)=0,\\ y\left(k\right)=Cx\left(k\right),\end{array}$$

(4)

where $x\in {ℝ}^n,y\in {ℝ}^1, u\in {ℝ}^m, A,B,C, N_{\gamma }$ are matrices of appropriate dimensions.

Equation (4) can be rewritten as:

$$x\left(k+1\right)=[A ⋮ N_1 ⋮ \cdots ⋮ N_{m }]\tilde{x}\left(k\right)+Bu\left(k\right),$$

(5)

where

$$\tilde{x}\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ x\left(k\right)·u_1\left(k\right)\\ \begin{array}{c}\dots \\ x\left(k\right)·u_m\left(k\right)\end{array}\end{array}\right]\in {ℝ}^{1\times \left(n\times \left(m+1\right)\right)}.$$

(6)

Thus, we get the representation:

$$\begin{array}{c}x\left(k+1\right)=\tilde{A}\tilde{x}\left(k\right)+Bu\left(k\right), \\ \tilde{A}=[A ⋮ N_1 ⋮ \cdots ⋮ N_{m }],\end{array}$$

(7)

where:

$$\begin{array}{c}\tilde{x}\left(k\right)=D^n·x\left(k\right),\\ D^n=\left[\begin{array}{c}{D}^{0n}\\ \dots \\ D^{nn}\end{array}\right], D^{0n}=\left[\begin{array}{ccc}1& \cdots & 0\\ ⋮& 1& ⋮\\ 0& \cdots & 1\end{array}\right], D^{in}=\left[\begin{array}{ccc}{u}_i\left(k\right)& \cdots & 0\\ ⋮& u_i\left(k\right)& ⋮\\ 0& \cdots & u_i\left(k\right)\end{array}\right], i=1,\dots ,n. \end{array}$$

For simplicity, we assume that $m=n.$ Otherwise, the matrices $N_{i }$are padded with zero rows so that the number of rows is equal to $n.$

Then Equation (4) can be represented as

$$x\left(k+1\right)=\tilde{A}{D}^{n}x\left(k\right)+Bu\left(k\right).$$

(8)

Furthermore, to identify the parameters of the system, we will carry out the transition from a state-space model to a linear input-output model. This transition is a standard procedure described, e.g., in [33]. The identification is carried out using the associative search algorithm together with the Moore-Penrose procedure [31,32], which delivers a solution to a system of linear equations with the statistical dependence of the components of the vector $\tilde{x}$.

Returning to the canonical form of the model results in an estimate of all parameters, i.e., the updated bilinear model. The ability to determine the system’s state and output for various control impacts enables the usage of identification models for predicting the approach to stability boundaries in advance.

3. Controllability and Observability Gramians of Discrete Stationary Bilinear Systems

Let the model (4) of a discrete stationary dynamical system be obtained as a result of identification using the algorithms described in Section 2. We will assume that it belongs to the class of MISO LTI systems.

Consider a MISO LTI discrete stationary dynamical system in the form:

$$\begin{array}{c}x\left(k+1\right)=Ax\left(k\right)+Bu\left(k\right), x\left(0\right)=0,\\ y\left(k\right)=Cx\left(k\right),\end{array}$$

(9)

where $x\in {ℝ}^n,y\in {ℝ}^1, u\in {ℝ}^m$

We will consider real matrices of the corresponding sizes $A,B,C$. Let us assume that the system (9) is stable, fully controllable and fully observable, all eigenvalues of matrix $A$ are different. Consider discrete algebraic Lyapunov equations associated with Equation (9) in the form:

$$\begin{array}{c}A{P}^c{A}^{*}+B{B}^{*}=P^c,\\ A^{*}{P}^{o}A+C^{*}C=P^o.\end{array}$$

Consider a bilinear discrete stationary dynamical system in the form:

$$\begin{array}{c}x\left(k+1\right)=Ax\left(k\right)+{\displaystyle \sum _{\gamma =1}^m}{N}_{\gamma }x\left(k\right)u_{\gamma }\left(k\right)+Bu\left(k\right), x\left(0\right)=0,\\ y\left(k\right)=Cx\left(k\right),\end{array}$$

(10)

where $x\in {ℝ}^n,y\in {ℝ}^1, u\in {ℝ}^m, A,B,C, N_{\gamma }$ are matrices of appropriate dimensions. One of the most important properties of control systems is the controllability. In [4,5], controllability and observability Gramians of discrete bilinear dynamical systems were introduced and iterative algorithms for their computation were proposed. Let us denote:

$$\begin{array}{c}{P}_1\left(k_1\right)=A^{k_1}B,\\ P_i\left(k_1,\dots ,k_i\right)=A^{k_i}\left[N_1{P}_{i-1}{N}_2{P}_{i-1}\dots N_m{P}_{i-1}\right], i\ge 2,\end{array}$$

The controllability Gramian of a bilinear system is defined as follows:

$$P={\displaystyle \sum }_{i=1}^{\infty }{\displaystyle \sum }_{k_i=0}^{\infty }\dots \dots {\displaystyle \sum }_{k_i=0}^{\infty }{P}_i{P}_i^T.$$

(11)

It was shown in [5] that if the spectrum of A belongs to the interior of the unit circle, then, under certain additional conditions, the solutions to the following two Lyapunov generalized matrix equations:

$$AP{A}^{*}+B{B}^{*}+{\displaystyle \sum }_{j=1}^m{N}_{j}P{N}_j^{*}=P,$$

(12)

$$A^{*}QA+C^{*}C+{\displaystyle \sum }_{j=1}^m{N}_{j}Q{N}_j^{*}=Q,$$

(13)

are the Gramians of controllability and observability. The Gramian of controllability of a bilinear system is the limiting solution

$$P=\underset{i\to \infty }{\lim }{P}_i$$

(14)

obtained as a result of the implementation of the following iterative procedure

$$A{P}_1{A}^{*}-P_1=-B{B}^{*},$$

$$A{P}_i{A}^{*}-P_i+{\displaystyle \sum }_{j=1}^m{N}_j{P}_{i-1}{N}_j^{*}=0, i=2,\dots \infty .$$

(15)

Similarly, the observability Gramian of a bilinear system is the limiting solution obtained by implementing a similar iterative procedure. The disadvantage of such procedures is that the resulting limiting solution is not always the corresponding Gramian of the bilinear system.

Our goal is to create improved iterative algorithms for calculating the Gramians of bilinear systems and to develop a method and algorithms for calculating the stability indices of bilinear systems based on them. To achieve this goal, it is proposed to change the computation paradigm by transferring computations from the matrixes to their elements in the course of iterations.

The very idea of the element-wise computation of Gramians is not new: for example, the method for vectorizing the solution of generalized matrix Lyapunov equations is based on it [3,5]. However, the calculation of sequences of numeric elements of Gramian matrixes will reveal new patterns of sequences behavior, for example, the formation of geometric progressions of elements. This will allow investigating the behavior of sequences for small, medium and large matrixes and develop new approaches to approximate calculations. Another argument in favor of the element-wise approach is that this approach to calculating the spectral decompositions of Gramians of linear continuous systems was previously proposed in [7] and has shown its effectiveness.

4. Iterative Methods for Calculating the Solutions of the Generalized Lyapunov Equations for Canonical State form Equations

Consider further the spectral methods of Gramians calculating for discrete linear systems. These methods were studied in early works [1,2,3,7,11] Consider a MIMO LTI discrete system reduced using a non-degenerate coordinate transformation to the diagonal form of the dynamics matrix

$$\begin{array}{c}x=T{x}_d x_d\left(k+1\right)=\Lambda x_d\left(k\right)+B_{d}u\left(k\right), y_d\left(k\right)=C_d{x}_d\left(k\right), \\ \Lambda =T^{-1}\mathrm{A}T, B_d=T^{-1}B, C_d=CT,\end{array}$$

(16)

or

$$A=\left[\begin{array}{ccc}{u}_1& u_2\cdots & u_n\end{array}\right]\left[\begin{array}{cccc}{z}_1& 0& 0& 0\\ 0& z_2& & 0\\ 0& & \ddots & \\ 0& 0& & z_n\end{array}\right]\left[\begin{array}{c}{v}_1^{*}\\ v_2^{*}\\ \begin{array}{c}⋮\\ v_n^{*}\end{array}\end{array}\right]=T\Lambda T^{-1}, TV=VT=\mathrm{I},$$

(17)

where the matrix $T$consists of the right eigenvectors $u_i$, and matrix $T^{-1}$ consists of the left eigenvectors $v_i^{*}$ corresponding to the eigenvalues $z_i$. The last equality is a condition for eigenvectors normalization.

In particular, in [8], the spectral decompositions of Gramian controllability and observability matrixes for LTI MIMO discrete stable systems with a simple spectrum are as follows:

$$P^c={\displaystyle \sum }_{k=1}^n{\displaystyle \sum }_{\rho =1}^n{P}_{k,\rho }^c, P_{k,\rho }^c={\displaystyle \sum }_{\eta =0}^{n-1}{\displaystyle \sum }_{j=0}^{n-1}\frac{z_k^j{z}_{\rho }^{\eta }}{\dot{N}\left(z_k\right)\dot{N}\left(z_{\rho }\right)}·\frac{1}{1-z_{\rho }{z}_k}{A}_{j}B{B}^T{A}_{\eta }^T,$$

(18)

$$P^o={\displaystyle \sum }_{k=1}^n{\displaystyle \sum }_{\rho =1}^n{P}_{k,\rho }^o, P_{k,\rho }^o={\displaystyle \sum }_{\eta =0}^{n-1}{\displaystyle \sum }_{j=0}^{n-1}\frac{z_k^j{z}_{\rho }^{\eta }}{\dot{N}\left(z_k\right)\dot{N}\left(z_{\rho }\right)}·\frac{1}{1-z_{\rho }{z}_k}{A}_{\eta }^T{C}^{T}C{A}_j,$$

(19)

where $z_k, z_{\rho }$ are the roots of the characteristic equation, $A_j,A_{\eta }^T$are the Faddeev matrices. For the Lyapunov equations of the same diagonalized systems of the form:

$$\begin{array}{c}\Lambda P_d^c{\Lambda }^{*}+B_d{B}_d^{*}=P_d^c,V_d=B_d{B}_d^{*},\\ {\Lambda }^{*}{P}_d^o\Lambda +C_d^{*}{C}_d=P_d^o,W_d=C_d^{*}{C}_d, \end{array}$$

we have the following formulas for spectral decomposition:

$$P_{d,\rho k}^c=\frac{1}{1-z_{\rho }{z}_k}{R}_k{B}_d{B}_d^{*}{R}_{\rho }^{*}, \forall z:\left|z\right|<1,$$

(20)

$$P_{d,\rho k}^o=\frac{1}{1-z_{\rho }{z}_k}{R}_k^{*}{C}_d^{*}{C}_d{R}_{\rho }, \forall z:\left|z\right|<1,$$

(21)

where $R_k$are the residues of the resolvent of the matrix $\Lambda$ in the eigenvalues of the matrix$z_k$.

The elements “ρk” of the sub-Gramian matrixes (20)–(21) satisfy the formulas:

$$p_{d\rho k}^c=\frac{1}{1-z_{\rho }{z}_k}{v}_{d\rho k, }{p}_{d\rho k}^o=\frac{1}{1-z_{\rho }{z}_k}{w}_{d\rho k. }$$

(22)

When transforming Equation (4) by decomposing the matrix A in its eigenvalues (16)–(17), we obtain the equations

$$\begin{array}{c}{x}_d\left(k+1\right)=\Lambda x_d\left(k\right)+{\displaystyle \sum _{\gamma =1}^m}{N}_{d\gamma }{x}_d\left(k\right)u_{\gamma }\left(k\right)+B_{d}u\left(k\right), x_d\left(0\right)=0,\\ y\left(k\right)=C_{d}x\left(k\right).k=0,1,2\dots \end{array}$$

(23)

$$N_{d\gamma }=T{N}_{\gamma }.$$

Definition 1.

Consider the following matrix and vector identities:

$$A\equiv {\displaystyle \sum }_{i,j}{a}_{ij}{1}_{ij},\left\{a\right\}\equiv {\displaystyle \sum }_i{a}_i{1}_i, 1_{ij}=e_i{e}_j^T, 1_i=e_i,$$

where the unit vector is as follows:

$$e_i={\left[\begin{array}{ccc}0\dots & 0& \begin{array}{ccc}1& 0& \dots 0\end{array}\end{array}\right]}^T.$$

We call the above decompositions of matrixes and vectors separable decompositions. The separability property means the change in the very paradigm of solution computing: the transition from the matrix-vector consideration to the element-wise one.

Let us derive a formula for spectral decomposition of Gramian for the dynamics matrix of the system, transformed into the canonical controllability form using a linear nondegenerate transformation of coordinates with the following matrix:

$$R_c^F, x=R_c^F{x}_c.$$

We assume full-controllability and full-observability conditions are fulfilled. Furthermore, we consider the channel “γ” MISO LTI of the linear system in the canonical controllability form:

$$\begin{array}{c}{x}_c\left(k+1\right)=A_c^F{x}_c\left(k\right)+b_{\gamma }^{F}u\left(k\right), x_c\left(0\right)=0,\\ y\left(k\right)=c_{\gamma }^{F}x\left(k\right),k=0,1,2\dots \end{array}$$

(24)

$$\begin{array}{c}{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], b_{\gamma }^F={\left[\begin{array}{ccc}0& 0& \begin{array}{ccc}\dots & 0& 1\end{array}\end{array}\right]}^T, N_{c\gamma }^F={\left(R_c^F\right)}^{-1}{N}_{\gamma }.\\ c_{\gamma }^F=\left[\begin{array}{ccc}{\xi }_0& {\xi }_1& \begin{array}{ccc}\cdots & {\xi }_{n-2}& {\xi }_{n-1}\end{array}\end{array}\right].\end{array}$$

The following relations are valid:

$$\begin{array}{c} R_c^F=\left[\begin{array}{ccc}B& AB& \begin{array}{ccc}\dots & \dots & A^{n-1}B\end{array}\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(R_c^F\right)}^{-1}A{R}_c^F=A^F,{\left(R_c^F\right)}^{-1}B=B^F,C{R}_c^F=C^F.\end{array}$$

Lemma 1.

Consider a linear discrete MISO system in the form (2) represented by equations in the canonical form of controllability of the form (24). Let us further consider the decomposition of the dynamics matrix $A^F$resolvent into a segment of the Faddeev series of the form

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

where:$N\left(z\right)$is a characteristic polynomial,$A_j^F$is Faddeev matrix, $j=1,2,\dots n$.

The elements of the last column of the matrix $A_j^F$ satisfy the statements:

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

(25)

Comment. 

Note, first that the decomposition of the resolvent in the Faddeev series form does not require calculating the eigenvalues of the dynamics matrix $A^F$. Second, the transfer function of the “$\gamma$” channel of the linear part is determined by the formula:

$$V_{\gamma }^{Fln}=\left[\begin{array}{ccc}{\xi }_0& {\xi }_1& \begin{array}{ccc}\dots & {\xi }_{n-2}& {\xi }_{n-1}\end{array}\end{array}\right]{\left(Iz-A^F\right)}^{-1} b_{\gamma }^F, b_{\gamma }^F={\left[\begin{array}{ccc}0& \begin{array}{cc}\dots & 0\end{array}& 1\end{array}\right]}^T,$$

hence it follows that it is determined only by the elements of the last column of the matrix $A^F.$

Proof. 

Consider the expansion of the resolvent of the matrix $A^F$ in the form of a segment of the Faddeev series

$${\left(Iz-A^F\right)}^{-1}=\frac{{\displaystyle \sum }_{j=0}^{n-1}{A}_j^F{z}^j}{N\left(z\right)}.$$

We will accept

$$N\left(z\right)=z^n+a_{n-1}{z}^{n-1}+\cdots a_{1}z+a_0, R_j=A_{j-1}^F, j=1,2,\dots n.$$

We apply the method of mathematical induction. An iterative algorithm for calculating the Faddeev matrixes and the coefficients of the characteristic equation has the form at the first step [17]:

$$a_n=1, R_n=A_{cn-1}^F=I,$$

at the step “k”:

$$a_{n-k}=-\frac{1}{k}tr\left(A^F{R}_{n-k+1}\right), R_{n-k}=a_{n-k}I+A^F{R}_{n-k+1}, k=1,2,\dots n.$$

Consider the formation of the last column of matrixes $A_{n-k}^F.$

The first step:

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

The second step:

$$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_1& -a_2& \dots & -a_{n-2}& 0\end{array}\right], \left\{a_{cn-2,n}^F\right\}={\left[\begin{array}{ccc}0& \dots & \begin{array}{ccc}0& 1& 0\end{array}\end{array}\right]}^T.$$

Suppose that for step “k − 1”, the last column of the matrix $A_{cn-\left(k-1\right)}^F$ has the form:

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

We introduce the notation:

$$A_c^F{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 is:

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

In accordance with the Faddeev—Le Verrier algorithm, we have:

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

□

Corollary 1.

Without loss of generality, we assume m = 1. The general formulas (9) of spectral decompositions for the controllability Gramians of the linear system transformed in the canonical controllability form, taking into account the Lemma 1, acquire a simpler form:

$$P^{cF}={\displaystyle \sum }_{k=1}^n{\displaystyle \sum }_{\rho =1}^n{P}_{k,\rho }^{cF}, P_{k,\rho }^{cF}={\displaystyle \sum }_{\eta =0}^{n-1}{\displaystyle \sum }_{j=0}^{n-1}\frac{z_k^j{z}_{\rho }^{\eta }}{\dot{N}\left(z_k\right)\dot{N}\left(z_{\rho }\right)}\frac{1}{1-z_{\rho }{z}_k}{e}_{j+1}{e}_{\eta +1}^T.$$

(26)

A similar approach can be used to derive a formula for spectral decompositions for observability Gramians of the MISO system transformed in the canonical form of observability. In this case, the following formulas are valid [18]:

$$\begin{array}{c}x=R_0^F{x}_o, \\ 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{array}$$

Let us use (16) and consider the formation of the expression $A_{oj}^{FT}$ $c^{\mathrm{F}T}{c}^{\mathrm{F}}{A}_{oj}^F.$ In accordance with duality principle, due to the complete controllability and observability, the following formula is valid:

$$c^F{A}_{oj}^F={\left(A_{cj}^F\right)}^T{b}^F.$$

(27)

Substituting (27) into (19), we obtain the expression:

$$P^{oF}={\displaystyle \sum }_{k=1}^n{\displaystyle \sum }_{\rho =1}^n{P}_{k,\rho }^{oF}={\displaystyle \sum }_{\eta =0}^{n-1}{\displaystyle \sum }_{j=0}^{n-1}\frac{z_k^j{z}_{\rho }^{\eta }}{\dot{N}\left(z_k\right)\dot{N}\left(z_{\rho }\right)}\frac{1}{1-z_{\rho }{z}_k}{e}_{j+1}{e}_{\eta +1}^T, P^{oF}=P^{cF}.$$

(28)

The Gramians of the original system are related to the Gramians of the systems transformed into canonical forms as follows

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

Please note that in this case the expressions of the Gramians and sub-Gramians of controllability and observability depend only on the eigenvalues of the dynamics matrix. In addition, the proposed approach using canonical forms made it possible to simplify the general formulas significantly.

5. Separable Spectral Method and Algorithm for Solving the Generalized Lyapunov Equation

Theorem 1.

Consider a MISO (multiple input single output) discrete bilinear stationary system in the form (4) [3,4,5,24].

(1)

Let the matrix A be a Hurwitz one with the simple spectrum,

(2)

Let the vector $u\left(k\right)$ be bounded:

$$\Vert u\left(k\right)\Vert =\sqrt[2]{{\displaystyle \sum }_{i=1}^m{\left|u_i\left(k\right)\right|}^2}<M.$$

(3)

There exist real numbers $\alpha ,\rho$ such that the following inequalities hold:

$$\begin{array}{c}\Vert A\Vert \le \alpha {\rho }^i, i=0,1,2,\dots ,\\ \rho <i,\alpha >0, \\ M<\sqrt[]{1-{\rho }^2{\alpha }^{-1}}.\end{array}$$

(4)

Suppose, in addition, the following conditions are satisfied:

$$\left|\frac{p_{ij}^{C bln\left(k+1\right)}}{p_{ij}^{C bln\left(k\right)}}\right|\le \overline{N}L<1, \forall k,\nu ,\mu ,j,\eta ,\gamma ,$$

(29)

where

$$\overline{N}=\sup n^2\underset{\nu ,\mu ,j,\eta ,\gamma }{\underbrace{|n_{d\nu i}^{\gamma }||n_{d\nu j}^{\gamma }|}}, L=\left|{\left(1-z_{\nu max}{z}_{\nu max}{}^{*}\right)}^{-1}\right|.$$

(30)

In (30),$z_{\nu max}$ denotes the maximum eigenvalue of the dynamics matrix of the linear part of the system.

Then, there also exists a uniquely following separable iterative spectral solution to the generalized Lyapunov equation (12) for the diagonalized system (16):

$$p_d^{\left(k\right)ij}=\left({\displaystyle \sum }_{\nu ,\mu }{p}_{d\nu \mu }^{\left(k-1\right)ij}\right)\left[{\left(1-z_{\nu }{z}_{\mu }\right)}^{-1}{n}_{d\nu i}^{\gamma }{n}_{d\mu j}^{\gamma }\right], k=2,3,\dots \infty .$$

(31)

$$\forall \nu ,\mu ,i,j=1,2,\dots n;\gamma =1,2,\dots m.$$

Sequences of partial sums (31) converge uniformly and absolutely to the corresponding elements of the solution matrix of the generalized Lyapunov equation (12) if the conditions of the theorem are satisfied. The controllability Gramian of the original bilinear system $P^{cbln}$is related to the controllability Gramian of the diagonalized bilinear system $P_d^{cbln}$ as follows

$$T{P}_d^{cbln}{T}^T=P^{cbln}.$$

(32)

Proof. 

Consider an iterative process, which develops the solution of (31).

Step 1. Let us consider the forming of the right-hand side of the generalized Lyapunov equation for the case m = 1. We do not need the matrix of the Lyapunov equation solution of the linear part; rather we need a separable spectral decomposition of this solution in the pair spectrum of the matrix [18]:

$$P_d^{bln\left(1\right)}={\displaystyle \sum }_{\gamma =1}^m{\displaystyle \sum }_{\nu =1}^n{\displaystyle \sum }_{\mu =1}^n\frac{v_{d,\gamma \nu \mu }}{1-z_{\nu }{z}_{\mu }}{1}_{\nu \mu }.$$

(33)

Step 2. Consider the formation of the right-hand side of the generalized Lyapunov equation with the example of the matrix $N_d^{\gamma }{1}_{ij}{\left(N_d^{\gamma }\right)}^T$:

$$N_d^{\gamma }{1}_{ij}{\left(N_d^{\gamma }\right)}^T={\displaystyle \sum }_{\nu =1}^n{\displaystyle \sum }_{\mu =1}^n{n}_{d,\nu i}^{\gamma }{n}_{\mu j}^{\gamma }{1}_{\nu \mu },$$

(34)

The solution of the Lyapunov equation takes at Step 2 the form:

$$P_d^{bln\left(2\right)}={\displaystyle \sum }_{\gamma =1}^m{\displaystyle \sum }_{\nu =1}^n{\displaystyle \sum }_{\mu =1}^n\frac{1}{1-z_{\nu }{z}_{\mu }}\left(n_{d,\nu i}^{\gamma }{n}_{d\mu j}^{\gamma }\right)p_{\nu \mu }^{bln\left(1\right)ij\gamma }{1}_{\nu \mu }.$$

Proceeding in a similar way and taking into account the summation of sub-Gramians over the index “γ”, we obtain a formula for calculating the matrix of the Gramian kernel of the order “k” at step “k”.

$$P_d^{bln\left(k\right)ij\gamma }={\displaystyle \sum }_{\nu ,\mu }{r}^{\left(k\right)ij\gamma }{p}_{\nu \mu }^{\left(k-1\right)ij\gamma }{1}_{\nu \mu }, r^{\left(k\right)ij\gamma }=\left[{\left(1-z_{\nu }{z}_{\mu }\right)}^{-1}{n}_{d\nu i}^{\gamma }{n}_{d\mu j}^{\gamma }\right],$$

(35)

$$p_{d\nu \mu }^{\left(k\right)ij\gamma }=\left({\displaystyle \sum }_{\nu ,\mu }{p}_{d\nu \mu }^{\left(k-1\right)ij\gamma }\right)\left[{\left(1-z_{\nu }{z}_{\mu }\right)}^{-1}{n}_{d\nu i}^{\gamma }{n}_{d\mu j}^{\gamma }\right], k=2,3,\dots \infty ,$$

$$\forall \nu ,\mu ,i,j=1,2,\dots n;\gamma =1,2,\dots m.$$

(36)

□

This proves the theorem’s statement about the iterative spectral decomposition of the solution in the case of solution convergence. Let us show that under the conditions of Theorem 2, the convergence of the sequences is absolute and uniform. To this end, we construct a majorizing sequence for the elements of the sub-Gramian matrixes. Suppose that conditions (29)–(30) are satisfied. For all converging sequences’ elements “ij”, the following conditions must be satisfied:

$$\left|p_d^{\left(k\right)ij\gamma }\right|\le \left|p_d^{\left(k-1\right)ij\gamma }\right|\dots \le \left|p_d^{bln\left(1\right)ij\gamma }\right|\le \underset{ij}{\underbrace{max}}\left|p_d^{bln\left(1\right)ij\gamma }\right|=M_{max}^{ij\gamma }.$$

Let us introduce the notation $M_{max}=\underset{ij\gamma }{\underbrace{max}}{M}_{max}^{ij\gamma }.$ For the matrix $N_d^{\gamma },$ the exact upper bound of the products exists:

$$n^2\left|n_{d\nu j}^{\gamma }\right|\left|n_{d\mu \eta }^{\gamma }\right|\le \overline{N}, \overline{N}>0, \forall \gamma ,\nu ,i,\mu ,j: \gamma ,\nu ,i,\mu ,j=1,2,\dots n.$$

In addition, in addition, due to the stability of the linear part, the exact upper bound exists for the functions:

$$L=\underset{\nu \mu }{\underbrace{max}}\left|{\left(1-z_{\nu max}{z}_{\nu max}{}^{*}\right)}^{-1}\right|, L>0, \forall \nu ,\mu ,: \nu ,\mu =1,2,\dots n,$$

where $z_{\nu max}$ is the maximum eigenvalue of the dynamics matrix of the system’s linear part.

Therefore, the following inequality holds:

$$\left|\frac{p_{dij}^{cbln\left(k+1\right)}}{p_{dij}^{cbln\left(k\right)}}\right|\le \overline{N}L \forall \gamma ,\nu ,i,\mu ,j: \gamma ,\nu ,i,\mu ,j=1,2,\dots n.$$

We choose a single majorant for all numerical sequences in the form:

$$\begin{array}{c}{S}_{0 }, S_{1 },\dots S_{k }.\\ S_{k }=M_{max}^{ij\gamma }{n}^2{\left[\underset{\nu \mu }{\underbrace{max}} \left|{\left(1-z_{\nu max}{z}_{\nu max}{}^{*}\right)}^{-1}\right| \underset{\nu \mu ij\gamma }{\underbrace{max}} \left|n_{d\nu i}^{\gamma }{n}_{d\mu j}^{\gamma }\right| \right]}^{k-1}. \\ k=2,3,\dots \infty . \end{array}$$

Obviously, with such a choice, according to (36), the following inequality holds:

$$\left|p_{d\nu \mu }^{\left(k\right)ij\gamma }\right|<S_{k }, \forall k,i,j,\nu ,\mu ,\gamma .$$

(37)

It follows thereof that under conditions (29)–(30), the inequality

$$\left|\frac{p_{dij}^{bln\left(k+1\right)}}{p_{dij}^{bln\left(k\right)}}\right|\le \frac{S_{k+1 }}{S_{k }}<1,$$

(38)

is valid.

The majorizing sequence for all sub-Gramians of the bilinear system forms a geometric progression with positive terms. In accordance with the convergence criterion for geometric progressions, it converges if the following condition is satisfied:

$$\overline{N}L<1.$$

In accordance with the Weerstrass test (37)–(38), the sequences of partial sums $p_{dij}^{cbln\left(k\right)}$ converge uniformly and absolutely. The uniqueness of the iterative solution under conditions (1)–(3) was proved in [5].

The Gramians method can be used simultaneously for state monitoring and control of large-scale power systems, in particular, for static stability analysis, for developing stability estimators, detecting dangerous free and forced oscillations, and assessing the resonant interaction of dangerous oscillations [1,7,8,9,10].

Algorithm of the spectral iterative solution of the generalized Lyapunov equation of the form (12) is as follows:

• Step 1. Calculate the spectrum of the dynamics matrix of the linear part, check the stability of the linear part, the absence of multiple roots of the characteristic equation. Find a non-degenerate coordinate transformation that transforms the dynamics matrix of the linear part into a diagonal matrix. Let us transform the equations of the bilinear system (9) to the diagonal form.
• Step 2. Check the fulfillment of conditions (1)–(4) of Theorem 1.
• Step 3. By analyzing conditions (4), we identify the numerical sequences of elements of the matrixes of the kernels of the spectral expansion of the matrixes of the solution of the generalized Lyapunov equation, which are critical from the point of view of convergence.
• Step 4. Using algorithm (34), we compute the sequences of elements “ij” of the matrixes of the kernels of the Gramian expansion of the bilinear system at each step. We aggregate the elements of the sequences into the matrixes of the kernels of the decomposition of the bilinear system Gramian. We estimate the accuracy of the solution.
• Step 5. Using Formula (32), we calculate the Gramian matrix of the original bilinear system.

Comment. In [3,5,20,24], various versions of the generalized Lyapunov equation solutions are proposed using conditions (1)–(3) given in Theorem 1., but the similarity transformation of the dynamics matrix of the linear part to the diagonal form is not used, and the separable spectral decomposition is not used solutions of the Lyapunov equation of the linear part and the generalized Lyapunov equation for a bilinear system. Such a technique allows one to switch from calculating decision matrixes at separate iterations to calculating sequences of their elements.

As is known [2], the necessary and sufficient condition for energy stability of the system in terms of the square of the H₂ norm of the linear system transfer function $\mathrm{G}\left(\mathrm{z}\right)$has the form:

$${\Vert \mathrm{G}\left(\mathrm{z}\right)\Vert }_2^2=\mathrm{trC}{P}^c{C}^{\mathrm{T}}=\mathrm{tr}{B}^{\mathrm{T}}{P}^{o}B<+\infty .$$

Therefore, we define the stability loss risk functional of a bilinear system as:

$$J\left(z_1,z_2,\dots ,z_n\right)=\mathrm{trC}{P}^{cbln}{C}^{\mathrm{T}}=\mathrm{tr}{B}^{\mathrm{T}}{P}^{obln}B.$$

(39)

As the system approaches the stability threshold caused by the approaching of the characteristic equation roots to the imaginary axis, the risk functional approaches the infinity. Let us define the acceptable risk of stability loss of the bilinear system as:

$$J^{\left(\gamma \right)}\left(z_1,z_2,\dots ,z_n,\gamma \right)=M_{\gamma perm}, \gamma =1,2,\dots m.$$

We will consider any system as conditionally unstable if all its roots lie in a unit circle, but the functional of the stability loss risk (39) exceeds the established acceptable risk value. Accordingly, we will consider the system conditionally stable if:

$$J^{\left(\gamma \right)}\left(z_1,z_2,\dots ,z_n,\gamma \right)<M_{\gamma perm}. \gamma =1,2,\dots m.$$

(40)

The inequalities (40) define a set of energy functionals, the boundedness whereof guarantees the BIBO stability of the bilinear system. Conditions (1)–(4) of the theorem are sufficient conditions for the BIBO stability of the bilinear system and, at the same time, sufficient conditions for the boundedness of the energy functionals $J^{\left(\gamma \right)}$.

It is easy to see that inequalities (40) determine the stability conditions for a bilinear energy system in a wider range as against the traditional selective modal analysis. The analysis of expressions (37), (38) shows that the elements of the numerical sequences for the Gramian of the bilinear system converge at different rates, the guaranteed estimate whereof is specified by expressions (29)–(30).

This estimate depends on the choice of the “γ” channel, on the values of the elements of the nonlinearity matrixes for a specific channel, and on the proximity of the product of two eigenvalues of the linear part to unity. Sufficient conditions for BIBO stability of a bilinear system were obtained earlier in [3,5,24]. Theorem 1 establishes additional sufficient conditions (29)–(30) that guarantee the existence of not only a matrix for the solution of the generalized Lyapunov equation, but also of complete controllability and observability properties for the bilinear system.

6. Case Studies

The increased requirements for speed, accuracy and control capabilities under conditions of uncertainty in the presence of various kinds of disturbances in the control systems of production processes in industry and the electric power industry have demonstrated the inadequacy of the capabilities of traditional approaches to the synthesis of automatic control systems. Methods of identification synthesis, in which models are developed on the basis of data mining and machine learning, are gaining more and more popularity [18].

The authors have developed an intelligent system designed to dynamically assess the state of facilities in the power system [34]. The system is underpinned by intelligent algorithms of grid dynamics identification with automatic on-line self-tuning based on the data from monitoring systems.

State estimation models for power facilities with on-line model tuning are based on data monitoring and application of a predictive method for state estimation—the associative search method.

The acquisition, storage, processing, displaying, analysis and documenting of the information are executed in real time based on the data from automated power generation, distribution and consumption systems and supervisory control, monitoring and accounting. Figure 1 demonstrates power dynamic estimation for a certain facility in the power system.

In Figure 1, we have:

• the black line represents the dynamics of the real process;
• the red line shows the result in a linear model;
• and the blue line shows the result obtained by associative search model.

Figure 1 shows how a more accurate estimate of a real process dynamics can be obtained using the associative search model, compared with the classical linear models.

7. Conclusions

Predictive bilinear models of discrete dynamical systems are obtained using the associative search algorithm. The method is based on the use of machine learning procedures and inductive knowledge (associative patterns) extraction from historical data. The method features high algorithmic speed, since the main computational load falls on the training stage.

According to the proposed scheme, we, at first, obtain a bilinear model of a nonlinear dynamic object, and then analyze the stability. The advantages of the scheme are the accuracy and speed of the identification algorithm. Section 6 demonstrates the operation of the associative search algorithm. It shows that for nonlinear systems the models obtained through this algorithm are more accurate as against the ones obtained using traditional linear techniques.

Furthermore, according to our scheme, separable spectral expansions of discrete Lyapunov equations are obtained for MISO LTI discrete dynamical systems governed by state equations in controllability and observability forms. A method and an algorithm for the element-wise solution of the generalized matrix Lyapunov equation are developed for discrete bilinear systems. The new method is a spectral version of the well-known iterative method used for solving this equation.

A sequence of values is calculated for a fixed element of the solution matrix. The element depends on the eigenvalues product of the dynamics matrix of the linear part and the elements of the nonlinearity matrixes. A sufficient condition for the convergence of all sequences is obtained, which is also a BIBO stability condition for a bilinear system.

The article discusses MIMO, MISO, and SISO classes of bilinear systems of the form (10) but does not consider bilinear systems with distributed parameters. In the future, the authors intend to extend the new method over this class of systems. Time-variant systems will be also investigated.