A Reduced-Order RNN Model for Solving Lyapunov Equation Based on Efficient Vectorization Method

With the trend of electronization of the power system, a traditional serial numerical algorithm is more and more difficult to adapt to the demand of real-time analysis of the power system. As one of the important calculating tasks in power systems, the online solution of Lyapunov equations has attracted much attention. A recursive neural network (RNN) is more promising to become the online solver of the Lyapunov equation due to its hardware implementation capability and parallel distribution characteristics. In order to improve the performance of the traditional RNN, in this study, we have designed an efficient vectorization method and proposed a reduced-order RNN model to replace the original one. First, a new vectorization method is proposed based on the special structure of vectorized matrix, which is more efficient than the traditional Kronecker product method. Second, aiming at the expanding effect of vectorization on the problem scale, a reduced-order RNN model based on symmetry to reduce the solution scale of RNN is proposed. With regard to the accuracy and robustness, it is proved theoretically that the proposed model can maintain the same solution as that of the original model and also proved that the proposed model is suitable for the Zhang neural network (ZNN) model and the gradient neural network (GNN) model under linear or non-linear activation functions. Finally, the effectiveness and superiority of the proposed method are verified by simulation examples, three of which are standard examples of power systems.

Introduction

With the trend of the electronic power system, the scale of system computing is increasing day by day, while the demand of real-time analysis and calculation in the process of system operation remains unchanged. Traditional serial algorithms cannot solve this contradiction well, so various parallel algorithms and distributed methods appear successively. In power system state estimation, Chen et al. (2017) have used the SuperLU_MT solver to estimate the state of the actual power grid, making full use of the parallel characteristics of multicore and multi-thread solver. Liu Z. et al. (2020) have fully explored the parallelism in the calculation of continuous power flow and applied the continuous Newton method power flow model to realize the parallel solution algorithm of continuous power flow based on GPU in large scale and multiple working conditions. Moreover, a novel distributed dynamic event-triggered Newton–Raphson algorithm is proposed to solve the double-mode energy management problem in a fully distributed fashion (Li et al., 2020). Similarly, Li Y. et al. (2019) proposed an event-triggered distributed algorithm with some desirable features, namely, distributed execution, asynchronous communication, and independent calculation, which can solve the issues of day-ahead and real-time cooperative energy management for multienergy systems. Given that software algorithms are essentially run by hardware, implementing functions directly from hardware is also an option for real-time computing. For example, Hafiz et al. (2020) proposed a real-time stochastic optimization of energy storage management using deep learning–based forecasts for residential PV applications, where the key of the real-time computation is the hardware controller. It is worth pointing out that compared with the aforementioned methods, the neural dynamics method has greater potential in the field of real-time calculation of power systems (Le et al., 2019), and its time constant can reach tens of milliseconds (Chicca et al., 2014) because of its parallel distribution characteristics and the convenience of hardware implementation.

The Lyapunov equation is widely used in some scientific and engineering fields to analyze the stability of dynamic systems (He et al., 2017; He and Zhang, 2017; Liu J. et al., 2020). In addition, the Lyapunov equation plays an important role in the controller design and robustness analysis of non-linear systems (Zhou et al., 2009; Raković and Lazar, 2014). In the field of power systems, the balanced truncation method, controller design, and stability analysis are also inseparable from the solution of the Lyapunov equation (Zhao et al., 2014; Zhu et al., 2016; Shanmugam and Joo, 2021). Therefore, many solving algorithms have been proposed to solve the Lyapunov equation. For example, Bartels and Stewart proposed the Bartels–Stewart method (Bartels and Stewart, 1972), which is a numerically stable solution. Lin and Simoncini (Lin and Simoncini, 2013) proposed the minimum residual method for solving the Lyapunov equation. Stykel (2008) used the low-rank iterative method to solve the Lyapunov equation and verified the effectiveness of the method through numerical examples. However, the efficiency of these serial processing algorithms is not high in large-scale applications and related real-time processing (Xiao and Liao, 2016).

Recently, due to its parallelism and convenience of hardware implementation, recurrent neural networks have been proposed and designed to solve the Lyapunov equation (Zhang et al., 2008; Yi et al., 2011; Yi et al., 2013; Xiao et al., 2019). The RNN mainly includes the Zhang neural network (ZNN) and gradient neural network (GNN) (Zhang et al., 2008). Most of the research studies on RNN focus on the improvement of model convergence. For example, Yi et al. (2013) point out that when solving a stationary or a non-stationary Lyapunov equation, the convergence of the ZNN is better than that of GNN. Yi et al. (2011) used a power-sigmoid activation function (PSAF) to build an improved GNN model to accelerate the iterative convergence of Lyapunov equation. In (Xiao and Liao, 2016), the sign-bi-power activation function (SBPAF) is used to accelerate the convergence of the ZNN model for solving the Lyapunov equation and the proposed ZNN model has finite-time convergence, which is obviously better than the previous ZNN and GNN models. In recent years, some studies have considered the noise-tolerant ZNN model. In Xiao et al. (2019), two robust non-linear ZNN (RNZNN) are established to find the solution of the Lyapunov equation under various noise conditions. Different from previous ZNN models activated by the typical activation functions (such as the linear activation function, the bipolar sigmoid activation function, and the power activation function), these two RNZNN models have predefined time convergence in the presence of various noises.

However, both GNN and ZNN need to transform the solution matrix from the matrix form to the vector form through the Kronecker product, which is called vectorization of the RNN model (Yi et al., 2011). The use of the Kronecker product will make the scale of the problem to be solved larger. As the size of the problem increases, the scaling effect of the Kronecker product becomes more obvious. The enlargement effect of the Kronecker product on the model size will not only lead to insufficient memory when the RNN is simulated on software but also make the hardware implementation of the RNN model need more devices and wiring, which increases the volume of hardware, the complexity of hardware production, and the failure rate of hardware. However, no study has discussed the order reduction of the RNN model.

It should be pointed out that the vectorized RNN model needs to be solved using a hardware circuit. However, as the relevant research of the RNN for solving the Lyapunov equation is still in the stage of theoretical exploration and improvement, there are no reports about hardware products of the RNN solver of the Lyapunov equation. Relevant studies (Zhang et al., 2008; Yi et al., 2011; Yi et al., 2013; Xiao and Liao, 2016; Xiao et al., 2019) simulate the execution process of the RNN hardware circuit through the form of software simulation, and this study also adopts this form. It is undeniable that the results of software simulation are consistent with those of hardware implementation. Therefore, the theoretical derivation and simulation results of the RNN in this article and in the literature (Zhang et al., 2008; Yi et al., 2011; Yi et al., 2013; Xiao and Liao, 2016; Xiao et al., 2019) can be extended to the scenarios of hardware implementation.

The RNN is used to solve the Lyapunov equation, and the ultimate goal is to develop an effective online calculation model to solve the Lyapunov equation, so it is of great significance to improve the calculation speed of the RNN. Current studies focus on improving the computational speed of the RNN by improving the convergence of RNN. However, how to efficiently realize vectorization of the RNN model is also a breakthrough to improve the computational efficiency of the RNN method. At present, the Kronecker product is generally used to transform the solution matrix into the vector form (Horn and Johnson, 1991). The Kronecker product actually performs multiple matrix multiplication operations, and the time complexity of multiplying two n×n matrices is O (n^3), so the time complexity of the Kronecker product increases rapidly as the scale increases. This means that the traditional matrix vectorization method based on the Kronecker product still has room for optimization.

In summary, this article proposes an efficient method for vectorizing the RNN model based on the special structure of the vectorized matrix, which is more efficient than the traditional expansion method by the Kronecker product. Aiming at the expanding effect of vectorization on the problem scale, a reduced-order RNN model based on symmetry was proposed for solving the time-invariant Lyapunov equation, and the validity and applicability of the reduced-order RNN model were proved theoretically. The main contributions of this article are as follows.

1) An efficient method for vectorization of RNN model is proposed. Compared with the traditional vectorization method, this method has higher efficiency and less time consumption.

2) The reduced-order RNN model for solving the Lyapunov equation based on symmetry is proposed, which greatly reduces the solution scale. It is proved theoretically that the proposed model can maintain the same solution as that of the original model. Meanwhile, it is proved theoretically that the proposed model is suitable for the ZNN model and GNN model under linear or non-linear activation functions.

3) Several simulation examples are given to verify the effectiveness and superiority of the proposed efficient method for vectorization of the RNN and the reduced-order RNN model. It is also verified that the neural dynamics method is suitable for solving the Lyapunov equation of power systems through three standard examples of power systems.

In order to show the contributions of this study more clearly, the logical graph using the RNN model for solving the Lyapunov equation is shown in Figure 1, and the main novelties and differences of this article from Refs Yi et al. (2011); Yi et al. (2013); Xiao and Liao (2016); Xiao et al. (2019) are shown in Table 1.

FIGURE 1

FIGURE 1. Logical graph of using a RNN model for solving the Lyapunov equation.

TABLE 1

TABLE 1. Main novelties and differences of this article from the relevant references.

In Table1, items and numbers correspond to the three steps of Figure 1. The relevant references include Yi et al. (2011); Yi et al. (2013); Xiao and Liao (2016); Xiao et al. (2019).

In conclusion, Refs (Yi et al., 2011; Yi et al., 2013; Xiao and Liao, 2016; Xiao et al., 2019) focus on constructing a stronger RNN model to improve the convergence and noise-tolerant ability, including using different activation functions and neural networks. However, this study focuses on the vectorization method and the reduced-order RNN model.

Problem Formulation and Related Work

Problem Formulation

Consider the following well-known Lyapunov equation (Yunong Zhang and Danchi Jiang, 1995)

$\begin{array}{c}{A}^{\text{T}}X\left(t\right)+X\left(t\right)A=-C\end{array},(1)$

where A ∈ ℝ^n×n is a constant stable real matrix and C ∈ ℝ^n×n is a constant symmetric positive-definite matrix. The objective is to find the unknown matrix X(t) ∈ ℝ^n×n to make the Lyapunov matrix Eq. 1 hold true. Let $X^{\ast }\in {\mathbb{R}}^{\mathit{n\times n}}$ denote the theoretical solution of Eq. 1. In addition, two of the most relevant works (i.e., GNN and ZNN models) are presented to solve the Lyapunov Eq. 1 in the following.

GNN

According to the principle of GNN (Yi et al., 2011) and combined with the characteristics of Lyapunov equation, a corresponding GNN model can be designed to solve the Lyapunov equation. The design steps are as follows:

First, construct an energy function based on norm as follows:

$\begin{array}{c}\mathrm{\Delta }=\frac{\parallel A^{\text{T}}X\left(t\right)+X\left(t\right)A+C{\parallel }_{\text{F}}^2}{2}\end{array}(2)$

where $\parallel .{\parallel }_{\text{F}}$ means F-norm. The minimum value of the energy function is the solution of the Lyapunov equation.

Second, based on the principle of the negative gradient descent of the GNN, the following formula can be constructed:

$\begin{array}{c}\begin{array}{c}-\frac{\partial \text{Δ}}{\partial X}=-A\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)-\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)A^{\text{T}}\end{array}\end{array}(3)$

By introducing the adjustable positive parameter $\text{γ}$, the following GNN model can be obtained:

$\begin{array}{c}\begin{array}{c}\dot{X}\left(\text{t}\right)=-\text{γ}A\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)\\ -\text{γ}\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)A^{\text{T}}\end{array}\end{array}(4)$

where $\text{γ}>0$, $X\left(t\right)\in {\mathbb{R}}^{\mathit{n\times n}}$, and $X\left(0\right)\in {\mathbb{R}}^{\mathit{n\times n}}$ is the initial value of X(t).

Finally, the conventional linear GNN (Eq. 4) can be improved into the following non-linear expression by employing a non-linear activation function array $\mathrm{ℱ}(\cdot )$:

$\begin{array}{l}\dot{X}\left(\text{t}\right)=-\text{γ}(A\mathrm{ℱ}\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)\\ +\mathrm{ℱ}\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)A^{\text{T}})\end{array}(5)$

where $\mathrm{ℱ}(\cdot ):{\mathbb{R}}^{\mathit{n\times n}}\to {\mathbb{R}}^{\mathit{n\times n}}$ denotes a matrix-valued activation function array of the GNN models. In this study, the bipolar sigmoid activation function (BPAF) is selected as the representative of the non-linear activation function of the GNN model for simulation because of its strong convergence (Yi et al., 2011). The expression of BPAF is as follows:

$\begin{array}{c}ℱ\left(x\right)=\frac{1-\exp \left(-\delta x\right)}{1+\exp \left(-\delta x\right)}\end{array}(6)$

where $\delta$ is a constant and $\delta >1$.

ZNN

First, following Zhang et al.’s design method (Zhang et al., 2002), we can define the following matrix-valued error function to monitor the solution process of Lyapunov Eq. 1:

$\begin{array}{c}E\left(t\right)=A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\end{array}(7)$

Then in view of the definition of $E\left(t\right)$ and the design formula $\text{d}E\left(t\right)/dt=-\text{γφ}\left(E\left(t\right)\right)$, the dynamic equation of the ZNN model for solving the online Lyapunov Eq. 1 is derived as follows:

$\begin{array}{c}{A}^{\text{T}}\dot{X}\left(\text{t}\right)+\dot{X}\left(\text{t}\right)A=-\gamma \mathrm{\phi }\left(A^{\text{T}}X\left(t\right)+X\left(t\right)A+C\right)\end{array}(8)$

where $\mathrm{\phi }(\cdot ):{\mathbb{R}}^{\mathit{n\times n}}\to {\mathbb{R}}^{\mathit{n\times n}}$ denotes a matrix-valued activation function array of the ZNN models. The definition of $\text{γ}$ in the ZNN model is the same as that in the GNN model.

In this study, the RNZNN-1 model is selected as the representative of the non-linear activation function of the ZNN model for simulation because of its strong convergence (Xiao et al., 2019). The expression of the non-linear activation function in the RNZNN-1 model is as follows:

$\begin{array}{c}\mathrm{\phi }\left(x\right)=(a_1{\left|x\right|}^{\eta }+a_2{\left|x\right|}^{\omega })\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x\right)+a_{3}x+a_4\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x\right)\end{array}(9)$

where design parameters $0<\eta <1$, $\omega >1$, $a_1>0$, $a_2>0$, $a_3\ge 0$, $a_4\ge 0$, and $\text{sign}\left(x\right)$ denotes the signum function.

An Efficient Method for Vectorization of RNN Model

General Method of Vectorizing RNN Model

The RNN model needs to be transformed to the vector form so that it can be used for software simulation (Li X. et al., 2019) and hardware implementation.

Vectorization of the GNN Model

Yi et al. (2011) pointed out that the vectorization of GNN model is as follows:

$\begin{array}{l}\text{vec}\dot{X\left(t\right)}=-\gamma \left(\left(A\otimes I\right)\mathrm{ℱ}\left(\left(A^{\text{T}}\otimes I\right)\text{vec}X\left(t\right)+\left(I\otimes A^{\text{T}}\right)\text{ vec}X\left(t\right)+\text{vec}C\right)+\left(I\otimes A\right)\mathrm{ℱ}\left(\left(A^{\text{T}}\otimes I\right)\text{vec}X\left(t\right)+\left(I\otimes A^{\text{T}}\right)\text{vec}X\left(t\right)+\text{vec}C\right)\right)\\ =-\gamma \left(\left(A\oplus A\right)\mathrm{ℱ}\left(\left(A^{\text{T}}\oplus A^{\text{T}}\right)\text{vec}X\left(t\right)+\text{vec}C\right)\right)\end{array}(10)$

where

$\begin{array}{c}A\oplus A=A\otimes I+I\otimes A\end{array}(11)$

$\begin{array}{c}{A}^{\text{T}}\oplus A^{\text{T}}=A^{\text{T}}\otimes I+I\otimes A^{\text{T}}\end{array}(12)$

where $\otimes$ means the Kronecker product. Given $X=\left[x_{\mathit{ij}}\right]\in {\mathbb{R}}^{\mathit{n\times n}}$, we can vectorize X as a column vector, vec(X) ∈ ℝ^n2×1, which is defined as $\text{vec}\left(X\right)={\left[x_{11},\dots ,x_{1n},x_{21},\dots ,x_{n1},\dots ,x_{nn}\right]}^{\text{T}}$.

Since the order of matrix addition and matrix transpose is interchangeable (Cheng and Chen, 2017),

$\begin{array}{c}{\left(Y+Z\right)}^{\text{T}}=Y^{\text{T}}+Z^{\text{T}}\end{array}(13)$

Applying this property to Eq. 11, we can get

$\begin{array}{c}{\left(A\otimes I+I\otimes A\right)}^{\text{T}}={\left(A\otimes I\right)}^{\text{T}}+{\left(I\otimes A\right)}^{\text{T}}\end{array}(14)$

According to Chen and Zhou (2012), the relationship between the matrix transpose and Kronecker product is as follows:

$\begin{array}{c}{\left(Y\otimes Z\right)}^{\text{T}}=Y^{\text{T}}\otimes Z^{\text{T}}\end{array}(15)$

Applying this property to Eq. 14, we can get

$\begin{array}{c}{\left(A\otimes I\right)}^{\text{T}}+{\left(I\otimes A\right)}^{\text{T}}=A^{\text{T}}\otimes I^{\text{T}}+I^{\text{T}}\otimes A^{\text{T}}\end{array}(16)$

Considering $I=I^{\text{T}}$ and combining Eqs 11, 12, 14 and 16, we can get

$\begin{array}{c}{\left(A\oplus A\right)}^{\text{T}}=A^{\text{T}}\oplus A^{\text{T}}\end{array}(17)$

Vectorization of the ZNN Model

The vectorization process of the ZNN model is similar to that of the GNN. Carry out Kronecker product on Eq. 8, and we can get:

$\begin{array}{c}\begin{array}{c}\left(A^{\text{T}}\oplus A^{\text{T}}\right)\text{vec}\dot{X\left(t\right)}=-\text{γφ}\left(\left(A^{\text{T}}\otimes I\right)\text{vec}X\left(t\right)+\left(I\otimes A^{\text{T}}\right)\text{vec}X\left(t\right)+\text{vec}C\right)\\ =-\text{γφ}\left(\left(A^{\text{T}}\oplus A^{\text{T}}\right)\text{vec}X\left(t\right)+\text{vec}C\right)\end{array}\end{array}(18)$

Vectorization of the RNN Model

By comparing Eqs 10, 17 and 18, it can be seen that the key of vectorization of the RNN model is to solve $A^{\text{T}}\oplus A^{\text{T}}$.

According to Eq. 12, the calculation of $A^{\text{T}}\oplus A^{\text{T}}$ can be divided into three steps:

1) Calculate $A^{\text{T}}\otimes I$

$\begin{array}{c}{A}^{\text{T}}\otimes I=\left[\begin{array}{ccccccccc}{a}_{11}& \dots & 0& \dots & \dots & \dots & a_{n1}& \dots & 0\\ ⋮& \ddots & ⋮& \dots & \dots & \cdots & ⋮& \ddots & ⋮\\ 0& \dots & a_{11}& \dots & \dots & \cdots & 0& \dots & a_{n1}\\ ⋮& ⋮& ⋮& \ddots & & & ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& & \ddots & & ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& & & \ddots & ⋮& ⋮& ⋮\\ a_{1n}& \dots & 0& \dots & \dots & \dots & a_{nn}& \dots & 0\\ ⋮& \ddots & ⋮& \dots & \dots & \dots & ⋮& \ddots & ⋮\\ 0& \dots & a_{1n}& \dots & \dots & \dots & 0& \dots & a_{nn}\end{array}\right]\end{array}(19)$

where $\left[\begin{array}{ccc}{a}_{ij}& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & a_{ij}\end{array}\right]$ is a diagonal matrix with $n$ rows and $n$ columns. $A^{\text{T}}\otimes I$ is a matrix with $n^2$ rows and $n^2$ columns.

2) Calculate $I\otimes A^{\text{T}}$

$\begin{array}{c}\begin{array}{c}I\otimes A^{\text{T}}=\left[\begin{array}{ccccccccc}{a}_{11}& \dots & a_{n1}& \dots & \dots & \dots & 0& \dots & 0\\ ⋮& \ddots & ⋮& \dots & \dots & \cdots & ⋮& \ddots & ⋮\\ a_{1n}& \dots & a_{nn}& \dots & \dots & \cdots & 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & & & ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& & \ddots & & ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮& & & \ddots & ⋮& ⋮& ⋮\\ 0& \dots & 0& \dots & \dots & \dots & a_{11}& \dots & a_{n1}\\ ⋮& \ddots & ⋮& \dots & \dots & \dots & ⋮& \ddots & ⋮\\ 0& \dots & 0& \dots & \dots & \dots & a_{1n}& \dots & a_{nn}\end{array}\right]\\ =\left[\begin{array}{cccc}{A}^{\text{T}}& 0& \dots & 0\\ 0& A^{\text{T}}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & A^{\text{T}}\end{array}\right]\end{array}\end{array}(20)$

where $I\otimes A^{\text{T}}$ is a matrix with $n^2$ rows and $n^2$ columns.

3) Add $A^{\text{T}}\otimes I$ to $I\otimes A^{\text{T}}$

$\begin{array}{c}\begin{array}{c}{A}^{\text{T}}\otimes I+I\otimes A^{\text{T}}=\\ \left[\begin{array}{ccccccc}2a_{11}& \dots & a_{n1}& \dots & a_{n1}& \dots & 0\\ ⋮& \ddots & ⋮& \dots & ⋮& \ddots & ⋮\\ a_{1n}& \dots & a_{nn}+a_{11}& \dots & 0& \dots & a_{n1}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ a_{1n}& \dots & 0& \dots & a_{11}+a_{nn}& \dots & a_{n1}\\ ⋮& \ddots & ⋮& \dots & ⋮& \ddots & ⋮\\ 0& \dots & a_{1n}& \dots & a_{1n}& \dots & 2a_{nn}\end{array}\right]\end{array}\end{array}(21)$

An Efficient Method for Vectorization of RNN Model

According to the previous analysis, no matter how the matrix $A$ is changed, the matrix structure of $A^{\text{T}}\oplus A^{\text{T}}$ is fixed. Based on the special structure of $A^{\text{T}}\oplus A^{\text{T}}$, an efficient method for vectorization of the RNN model is proposed in this article. The steps are as follows.

1) Create a matrix with $n^2$ rows and $n^2$ columns named K and fill $K$ with the elements of $A$ according to Eq. 19.

2) Fill $K$ with the elements of $A$ according to Eq. 20.

3) Add the corresponding element of $A$ to the diagonal element of $K$ according to Eq. 21.

The vectorization method of the RNN model proposed in this article is still based on the Kronecker product, but the time complexity is greatly reduced. Because the vectorization method proposed in this article replaces matrix multiplication with assignment and addition.

The Reduced-Order RNN Model for Solving Lyapunov Equations Based on Symmetry

Since the solution of Eq. 1, $X^{\text{*}}$, is always symmetric, as long as the upper trigonometric elements of $X^{\text{*}}$ are solved, the lower trigonometric elements of $X^{\text{*}}$ can be obtained correspondently, which can greatly reduce the computational amount of solving the Lyapunov equations. Based on this idea, a reduced-order RNN model for solving the Lyapunov equation based on symmetry is proposed in this article.

The Reduced-Order ZNN Model With Linear Activation Function

Vectorization

Let’s consider a ZNN model with linear activation function after vectorization. The formula is as follows:

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{cccccc}{k}_{11}& k_{12}& \dots & k_{1,n+1}& \dots & k_{1,n^2}\\ k_{21}& k_{22}& \dots & k_{2,n+1}& \dots & k_{2,n^2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+1,n+1}& \dots & k_{n+1,n^2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,n+1}& \dots & k_{n^2,n^2}\end{array}\right]\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{n+1}}\\ ⋮\\ \dot{x_{n^2}}\end{array}\right]=\\ -\text{γ}\left(\left[\begin{array}{cccccc}{k}_{11}& k_{12}& \dots & k_{1,n+1}& \dots & k_{1,n^2}\\ k_{21}& k_{22}& \dots & k_{2,n+1}& \dots & k_{2,n^2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+1,n+1}& \dots & k_{n+1,n^2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,n+1}& \dots & k_{n^2,n^2}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n+1}\\ ⋮\\ x_{n^2}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\right)\\ \end{array}\end{array}(22)$

where $K=A^{\text{T}}\oplus A^{\text{T}}$.

For the convenience of later discussion, S∈$\mathbb{R}$^n×n is constructed. Assign the following values to $S$ as follows:

$\begin{array}{c}S=\left[\begin{array}{cccc}1& 2& \cdots & n\\ n+1& n+2& \cdots & 2n\\ ⋮& ⋮& \ddots & ⋮\\ n\left(n-1\right)+1& n\left(n-1\right)+2& \cdots & n^2\end{array}\right]\end{array}(23)$

Each element of $S$ is the index number of the element of $A$ at the same position.

We can expand $X^{\text{∗}}$ to vecX^∗ ∈ $\mathbb{R}$^n2×1. Assuming that $x_2^{\ast }$ and $x_{n+1}^{\ast }$ are, respectively, the elements of the 1st row and the n + 1th row of $\text{vec}{X}^{\text{∗}}$. Due to the symmetry, $x_2^{\ast }$ will be equal to $x_{n+1}^{\ast }$.

Reduce the Column Number of $\mathit{K}$

If the Kronecker product is directly carried out on Eq. 1, then

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{cccccc}{k}_{11}& k_{12}& \dots & k_{1,n+1}& \dots & k_{1,n^2}\\ k_{21}& k_{22}& \dots & k_{2,n+1}& \dots & k_{2,n^2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+1,n+1}& \dots & k_{n+1,n^2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,n+1}& \dots & k_{n^2,n^2}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n+1}\\ ⋮\\ x_{n^2}\end{array}\right]\\ =-\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\end{array}\end{array}(24)$

We can use $K\text{vec}X=-\text{vec}C$ to express Eq. 24. Lan (2017) points out that if $A$ is stable and $C$ is symmetrically positive definite, then the Lyapunov Eq. 1 has a unique symmetric positive definite solution. Therefore, the $K$ matrix of Eq. 24 must be invertible.

Multiply both sides of Eq. 22 by the inverse matrix of $K$, then we get

$\begin{array}{c}\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{n+1}}\\ ⋮\\ \dot{x_{n^2}}\end{array}\right]=-\mathrm{\gamma }\left(\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n+1}\\ ⋮\\ x_{n^2}\end{array}\right]+K^{-1}\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\right)\end{array}(25)$

From Eq. 25, we can see that $K^{-1}\text{vec}C$ is the solution of the Lyapunov Eq. 1, which means $K^{-1}\text{vec}C=\text{vec}{X}^{\text{*}}$. Since $X^{\text{*}}$ is a symmetric matrix, the differential equations of $x_2\left(t\right)$ and $x_{n+1}\left(t\right)$ are the same. If $x_2\left(0\right)$ and $x_{n+1}\left(0\right)$ are equal, then the time domain trajectories of $x_2\left(t\right)$ and $x_{n+1}\left(t\right)$ are the same, namely, $\dot{x_2}\left(t\right)=\dot{x_{n+1}}\left(t\right)$ and $x_2\left(t\right)=x_{n+1}\left(t\right)$. Therefore, for Eq. 22, the column $n+1$ of $K$ can be added to the second column. Similarly, the same operation of column addition can be performed on the other columns in the symmetric positions. So the column number of $K$ reduces to $0.5\left(n+1\right)n$. Eq. 22 becomes

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}\left(n+1\right)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}\left(n+1\right)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+\mathrm{1,0.5}\left(n+1\right)n}\\ ⋮& ⋮& ⋮& ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,0.5\left(n+1\right)n}\end{array}\right]\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{0.5\left(n+1\right)n}}\end{array}\right]=\\ -\text{γ}\left(\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}\left(n+1\right)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}\left(n+1\right)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+\mathrm{1,0.5}\left(n+1\right)n}\\ ⋮& ⋮& ⋮& ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,0.5\left(n+1\right)n}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{0.5\left(n+1\right)n}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\right),\\ \end{array}\end{array}(26)$

Reduce the Row Number of $\mathit{K}$

When the steady state is considered, the differential term of Eq. 26 is 0, and we can get:

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}\left(n+1\right)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}\left(n+1\right)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+\mathrm{1,0.5}\left(n+1\right)n}\\ ⋮& ⋮& ⋮& ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,0.5\left(n+1\right)n}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{0.5\left(n+1\right)n}\end{array}\right]\\ =-\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\end{array}\end{array}(27)$

As mentioned before, if $A$ is stable and $C$ is symmetrically positive definite, then Lyapunov Eq. 1 has a unique symmetric positive definite solution. Therefore, the number of equations should be the same as the number of variables (Cheng and Chen, 2017), namely, the rank of the coefficient matrix of Eq. 27 is equal to $0.5\left(n+1\right)n$, which means the row vector set of the coefficient matrix of Eq. 27 is linearly correlated.

We can construct the augmented matrix of Eq. 27 and name it as $G$. If we define the first row of $G$ as the vector ${\alpha }_1$, the second row as the vector ${\alpha }_2$, and so on, the row $n^2$ is defined as the vector ${\alpha }_{n^2}$.

According to the knowledge of linear algebra, the vector set ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_{n^2}$ is linearly dependent only if at least one of the vectors in the set can be represented linearly by the other vectors.

Let us define the vectors which can be represented linearly by the other vectors as the redundant vectors. Suppose that

$\begin{array}{c}{\alpha }_{n^2}=h_1{\alpha }_1+h_2{\alpha }_2+\dots +h_{n^2-1}{\alpha }_{n^2-1}\end{array}(28)$

where $h_1,h_2,\dots ,h_{n^2-1}$ are real numbers and at least one of them is not equal to 0. Then ${\alpha }_{n^2}$ is a redundant vector. As long as the redundant vectors are found out and the equations of their corresponding rows are deleted, a new augmented matrix with full row rank can be obtained.

According to the aformentioned analysis, as long as $A$ is stable and $C$ is symmetrically positive definite, then the redundant vectors must exist, which means there are always some vectors that satisfy Eq. 28. However, $A$ and $C$ are independent of each other. Considering there are always some vectors satisfying Eq. 28 in the case of any stable $A$ and any symmetrically positive definite $C$, there is only one possibility that for some redundant vector, there is another vector that is equal to it, and the row indexes of both of them are symmetric in the matrix $S$. Only in this way, based on the symmetric characteristics of matrix $C$, can the redundant vectors always satisfy Eq. 28 when $A$ and $C$ are independent of each other.

For a redundant vector, its own row index and the row index of another vector equal to it can form a pair of indexes. We can use these index pairs to find the redundant vectors and delete the corresponding rows. In general, in an index pair, the equation corresponding to the index whose value is larger is selected for deletion.

After the row deletion, the row number of $K$ also reduces to $0.5\left(n+1\right)n$. Eq. 26 becomes

$\begin{array}{c}\begin{array}{c}\\ \left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}\left(n+1\right)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}\left(n+1\right)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{0.5\left(n+1\right)n,1}& k_{0.5\left(n+1\right)n,2}& \dots & k_{0.5\left(n+1\right)n,0.5\left(n+1\right)n}\end{array}\right]\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{0.5\left(n+1\right)n}}\end{array}\right]=-\text{γ}\left(\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}\left(n+1\right)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}\left(n+1\right)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{0.5\left(n+1\right)n,1}& k_{0.5\left(n+1\right)n,2}& \dots & k_{0.5\left(n+1\right)n,0.5\left(n+1\right)n}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{0.5\left(n+1\right)n}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{0.5\left(n+1\right)n}\end{array}\right]\right)\\ \end{array}\end{array}(29)$

We can use $K_r\text{vec}\dot{X_r}=-\text{γ}\left(K_r\text{vec}{X}_r+\text{vec}{C}_r\right)$ to express Eq. 29.

Reduced-Order GNN Model With Linear Activation Function

Consider a GNN model with linear activation function after vectorization. The formula is as follows:

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{n+1}}\\ ⋮\\ \dot{x_{n^2}}\end{array}\right]=-\text{γ}{\left[\begin{array}{cccccc}{k}_{11}& k_{12}& \dots & k_{1,n+1}& \dots & k_{1,n^2}\\ k_{21}& k_{22}& \dots & k_{2,n+1}& \dots & k_{2,n^2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+1,n+1}& \dots & k_{n+1,n^2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,n+1}& \dots & k_{n^2,n^2}\end{array}\right]}^{\text{T}}\\ \left(\left[\begin{array}{cccccc}{k}_{11}& k_{12}& \dots & k_{1,n+1}& \dots & k_{1,n^2}\\ k_{21}& k_{22}& \dots & k_{2,n+1}& \dots & k_{2,n^2}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ k_{n+\mathrm{1,1}}& k_{n+\mathrm{1,2}}& \dots & k_{n+1,n+1}& \dots & k_{n+1,n^2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ k_{n^2,1}& k_{n^2,2}& \dots & k_{n^2,n+1}& \dots & k_{n^2,n^2}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{n+1}\\ ⋮\\ x_{n^2}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{n+1}\\ ⋮\\ c_{n^2}\end{array}\right]\right)\\ \\ \end{array}\end{array}(30)$

When the steady state is considered, the differential term of Eq. 30 is 0, and Eq. 30 is changed into Eq. 24. From the aforementioned derivation, it can be known that both $K\text{vec}X=-\text{vec}C$ and $K_r\text{vec}{X}_r=-\text{vec}{C}_r$ can obtain the solution of Lyapunov Eq. 1. On this basis, an attempt is made to construct a reduced-order GNN model with linear activation function after vectorization, as follows:

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}\dot{x_1}\\ \dot{x_2}\\ ⋮\\ \dot{x_{0.5(n+1)n}}\end{array}\right]=-\text{γ}{\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}(n+1)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}(n+1)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{0.5(n+1)n,1}& k_{0.5(n+1)n,2}& \dots & k_{0.5(n+1)n,0.5(n+1)n}\end{array}\right]}^T\\ \left(\left[\begin{array}{cccc}{k}_{11}& k_{12}& \dots & k_{\mathrm{1,0.5}(n+1)n}\\ k_{21}& k_{22}& \dots & k_{\mathrm{2,0.5}(n+1)n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{0.5(n+1)n,1}& k_{0.5(n+1)n,2}& \dots & k_{0.5(n+1)n,0.5(n+1)n}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{0.5(n+1)n}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{0.5(n+1)n}\end{array}\right]\right)\\ \end{array}\end{array}.(31)$

When the steady state is considered, Eq. 31 is changed into $K_r\text{vec}{X}_r=-\text{vec}{C}_r$, which means the solution of Lyapunov Eq. 1 can be finally obtained by solving Eq. 31.

Reduced-Order RNN Model With Non-Linear Activation Functions

Before every non-linear activation function is introduced into the linear RNN, it will be theoretically proved that their introduction can guarantee the correct convergence of RNN. However, the reduced-order RNN in this article does not change the structure of RNN, but only changes the size of RNN. Because both the problem scale and the specific values of the matrix are generally expressed in symbolic form in the theoretical derivation (Xiao and Liao, 2016; Xiao et al., 2019), and the reduced-order RNN model in this article is still applicable to the relevant theoretical proof of introducing the non-linear activation functions into the linear RNN. In other words, the reduced-order method in this article can be applied to the RNN model with non-linear activation functions.

The Generation of the Reduced-Order RNN Model

The steps for generating the reduced-order RNN model are as follows:

1) $S$ is constructed, and the construction logic is as described above.

2) Considering that the indexes in the symmetric positions of $S$ can form $\frac{n\left(n-1\right)}{2}$ index pairs, we construct a matrix $L$ and fill $L$ with $\frac{n\left(n-1\right)}{2}$ index pairs. It is important to note that all of the elements in the first column of $L$ must be the upper trigonometric elements (excluding diagonal elements) of $S$. For the convenience of presentation, suppose $M$ is the first column of L, and $N$ is the second column of matrix L.

3) For $K$ and $\text{vec}C$ of RNN model, the rows corresponding to the element values of $N$ are deleted.

4) For $K$ and $\text{vec}C$ obtained in step 3, add the symmetric columns according to the element values of $M$ and $N$, and the columns obtained by the addition will replace the columns corresponding to the element values of $M$, while the columns corresponding to the element values of $N$ will be deleted. For $\text{vec}\dot{X\left(t\right)}$ and $\text{vec}X\left(t\right)$, the rows corresponding to the element values of $N$ will be deleted.

It should be pointed out that in mathematical proof, if the order of row reduction and column reduction is exchanged, the correctness of the reduced-order RNN model cannot be proved or another proof method is needed to complete the proof. However, in the case that the mathematical proof has been completed, the order of row reduction and column reduction does not affect the final result, since we know in advance which rows and columns are to be deleted. In the aformentioned steps of generating the reduced-order RNN model, the reason why we carry out step 3 first is that it can reduce the computational amount of the column addition to achieving higher computational efficiency.

The Significance of the Reduced-Order RNN Model

In order to better explain the value and significance of the reduced-order RNN model proposed in this article, the differences before and after the order reduction are shown from the perspectives of software simulation and hardware implementation, respectively. For the convenience of discussion, the GNN is taken as an example to illustrate.

Simulation on the Software

We use the ode45 function of MATLAB to solve the GNN model after vectorization. By comparing Eq. 30 and Eq. 31, it can be seen that the memory requirement of the reduced-order GNN model is much smaller than that of the original GNN model. Therefore, the reduced-order GNN model greatly alleviates the problem of insufficient memory that may occur in the software simulation of the GNN.

Hardware Implementation

When we use the traditional GNN model, the structure of the circuit diagram is shown as Figure 2 (Yi et al., 2011).

FIGURE 2

FIGURE 2. Circuit schematics which realizes GNN.

Where$M=\left[m_{ij}\right]\in R^{n^2\times n^2}=A\oplus A$, $P={\left[p_1,p_2,\dots ,p_{n^2}\right]}^{\text{T}}=\text{vec}C\in R^{n^2\times 1},$$Y={\left[y_1,y_2,\dots ,y_{n^2}\right]}^{\text{T}}=\text{vec}X\in R^{n^2\times 1}.$

When we use the reduced-order GNN model, the structure of the circuit diagram is shown as Figure 2, except that the $n^2$ in the diagram becomes $0.5\left(n+1\right)n$.

So the reduced-order GNN model greatly reduces the number of devices and wiring required for the hardware realization of GNN model, which is conducive to reducing the volume of hardware, the complexity of hardware production, and the failure rate of hardware.

Illustrative Verification

The simulation examples in this article are all completed on the MATLAB 2013b platform. In this article, the ode45 function of MATLAB is used to simulate the iterative process of RNN (Zhang et al., 2008). The corresponding computing performance is tested on a personal computer with Intel Core i7-4790 CPU @3.2GHz and 8 GB RAM.

Since there are great differences between software and hardware in the principle of realizing the integral function, there will be a big gap between the time cost in simulating the RNN process using software and the time cost in implementing the RNN model using hardware. Considering the research on the RNN model used to solve the Lyapunov equation is still in the stage of theoretical exploration, and has not reached the stage of hardware production for the time being, this article does not discuss the influence of the proposed reduced-order RNN model on the time consuming of RNN.

The Reduced-Order RNN Model for Solving Lyapunov Equation Based on Symmetry

Example 1

Let us consider the Lyapunov Eq. 1 with the following coefficient matrices:

$A=\left[\begin{array}{ccc}-11& 2& 3\\ 4& -7& 6\\ 1& 8& -12\end{array}\right]\text{and}C=I^{3\times 3}$

where $A$ is similar to Example II in Xiao et al. (2019). However, $A$ and $C$ of Example II in Xiao et al. (2019) do not fit the definition of Eq. 1, so we change $A$ a little bit and set $C$ to be the identity matrix.

In this example, we set $\text{γ}=10$, $\eta =0.25$, $\omega =4$, $a_1=a_2=a_3=a_4=1$, and $\delta =4$ (Yi et al., 2011; Xiao et al., 2019).

In order to demonstrate the advantages of the reduced-order RNN model, this article compares the performance of the reduced-order RNN and the original-order RNN, as shown in Table 2. In Table 2, linearity and non-linearity, respectively, mean the linear activation function and the non-linear activation function. Scale means the row number of $\text{vec}X$ or $\text{vec}{X}_r$. Proportion means the row number of $\text{vec}{X}_r$ divided by the row number of $\text{vec}X$. The ZNN F-norm and GNN F-norm, respectively, mean $\parallel A^{\text{T}}X+XA+C{\parallel }_F$ at the end of the simulation of ZNN and GNN.

TABLE 2

TABLE 2. Comparison of the original-order RNN and the reduced-order RNN under Example 1.

In order to study the effect of order reduction method proposed in this article on the convergence of the RNN model, the F-norm curves of the original-order RNN model and the reduced-order RNN model are drawn, as shown in Figure 3. In Figure 3A, LAF means the ZNN model with linear activation functions. NAF means the ZNN model with the non-linear activation function of Eq. 9. In Figure 3B, LAF means the GNN model with linear activation functions. BPAF means the GNN model with the non-linear activation function of Eq. 6. In both Figure 3A and Figure 3B, F-norm refers to $\parallel A^{\text{T}}X\left(t\right)+X\left(t\right)A+C{\parallel }_F$.

FIGURE 3

FIGURE 3. Comparison of the convergence between the original-order and the reduced-order RNN models under Example 1. (A) is ZNN, (B) is GNN.

Example 2

To enlarge the scale of the example, we consider Lyapunov Eq. 1 with the following coefficient matrices:

$A=\left[\begin{array}{ccccc}-17& 3& 4& 5& 6\\ 3& -17& 3& 4& 5\\ 4& 3& -17& 3& 4\\ 5& 4& 3& -17& 3\\ 6& 5& 4& 3& -17\end{array}\right]\text{and}C=I^{5\times 5}$

where $A$ is similar to Example III in Xiao et al. (2019). However, $A$ and $C$ of Example III in Xiao et al. (2019) do not fit the definition of Eq. 1, so we change $A$ a little bit and set $C$ to be the identity matrix.

In this example, RNN’s model parameters are the same as Example 1. Similar to Example 1, we can get Table 3 and Figure 4. The definitions of all nouns in Table 3 are the same as those in Table 2, and the definitions of all nouns in Figure 4 are the same as those in Figure 3.

TABLE 3

TABLE 3. Comparison of the original-order RNN and the reduced-order RNN under Example 2.

FIGURE 4

FIGURE 4. Comparison of the convergence between the original-order and the reduced-order RNN models under Example 2. (A) is ZNN, (B) is GNN.

Example 3

$A=\left[\begin{array}{cccccccccc}-50& 5& 2& 5& 9& 5& 8& 3& 6& 10\\ 8& -51& 10& 10& 7& 1& 4& 4& 3& 8\\ 4& 3& -51& 2& 4& 10& 3& 9& 8& 5\\ 6& 10& 8& -49& 6& 10& 5& 1& 2& 5\\ 2& 2& 9& 2& -47& 5& 1& 1& 7& 5\\ 7& 9& 9& 2& 1& -47& 2& 2& 2& 4\\ 3& 6& 1& 9& 3& 4& -42& 7& 4& 6\\ 7& 10& 4& 6& 2& 10& 10& -44& 7& 6\\ 7& 1& 3& 6& 2& 4& 6& 7& -44& 9\\ 8& 5& 9& 2& 3& 2& 1& 5& 1& -44\end{array}\right]\text{and}C=I^{10\times 10}$

A 10∗10 matrix is randomly generated and then α-shift is applied to the matrix to make it stable (Yang et al., 1993), which is the generation process of $A$ of Example 3 in this article. $C$ is set to be the identity matrix.

In this example, RNN’s model parameters are the same as Example 1. Similar to Example 1, we can get Table 4 and Figure 5. The definitions of all nouns in Table 4 are the same as those in Table 2, and the definitions of all nouns in Figure 5 are the same as those in Figure 3.

TABLE 4

TABLE 4. Comparison of the original-order RNN and the reduced-order RNN under Example 3.

FIGURE 5

FIGURE 5. Comparison of the convergence between the original-order and the reduced-order RNN models under Example 3. (A) is ZNN, (B) is GNN.

Based on the information in the aforementioned three tables, we can draw the following conclusions:

a) The reduced-order RNN model has a very obvious effect, with the scale reduced by about 33–45%. Moreover, the effect of the reduced-order RNN model becomes more obvious with the increase in the size of the example. According to $\underset{n\to \infty }{\lim }\frac{0.5n(n+1)}{n^2}=0.5$, it can be seen that when the size of the example is larger, the percentage of the scale decrease is closer to 50%. The reduced-order RNN model not only greatly alleviates the problem of insufficient memory in the software simulation of RNN but also greatly reduces the number of devices and wiring required for the hardware realization of RNN model, which is conducive to reducing the volume of hardware, the complexity of hardware production, and the failure rate of hardware.

b) Under different case scales, whether it is ZNN or GNN, whether it is linear activation function or non-linear activation function, the steady-state errors of the reduced-order RNN model are very close to 0, which means the reduced-order RNN model can always converge to the correct solution of the Lyapunov equation. This indicates that the reduced-order RNN model is applicable to ZNN and GNN, as well as the scenarios of linear activation function and non-linear activation function, which is consistent with the theoretical derivation results above.

c) Under different case scales, the difference in the steady-state accuracy between the reduced-order RNN and the original-order RNN is very small, indicating that the reduced-order RNN basically does not affect the steady-state accuracy of RNN.

Based on the information in the aforementioned three figures, we can draw the following conclusions:

a) Under different case scales, the reduced-order RNN models with linear or non-linear activation functions either have a little effect on the iterative convergence characteristics or enhance the convergence at the beginning of the iteration process and have a little effect on the convergence at the end of it.

b) Under the non-linear activation functions, the convergence of the ZNN model is always stronger than that of the GNN model when other conditions are fixed.

c) Under the linear activation function, the convergence of the ZNN model is weaker than that of the GNN model when the size of the examples is small (e.g., Example 1 and Example 2). The convergence of the linear ZNN model is stronger than that of the linear GNN model when the size of the examples is large (e.g., Example 3).

d) For both ZNN and GNN, the convergence of the RNN model with non-linear activation function is always stronger than that of the linear RNN model.

e) With the increase in the size of the examples, the convergence of ZNN is basically unchanged, while the convergence of GNN will become significantly worse.

Example 4

In order to verify the applicability of neural dynamics method to the power system, the corresponding Lyapunov equation describing system controllability is generated for the IEEE three-machine nine-node system according to the principle of the balanced truncation method in (Zhao et al., 2014). The input signal is the rotor speed deviation and the output signal is the auxiliary stabilizing signal (Zhu et al., 2016). It should be noted that the IEEE standard systems used in this article come from the examples of the PST toolkit (Lan, 2017), and the linearization process of the system is realized by the svm_mgen.m of PST toolkit. We set $\text{γ}=10$. A ∈$\mathbb{R}$^15×15 and C ∈$\mathbb{R}$^15×15 of the Lyapunov equation are detailed in Supplementary Material.

In this example, the linear ZNN was selected for testing. Similar to Example 1, we can get Table 5 and Figure 6. The definitions of all nouns in Table 5 are the same as those in Table 2 and the definitions of all nouns in Figure 6 are the same as those in Figure 3.

TABLE 5

TABLE 5. Comparison of the original-order ZNN and the reduced-order ZNN under Examples 4-6.

FIGURE 6

FIGURE 6. Comparison of the convergence between the original-order and the reduced-order ZNN models under Example 4.

Example 5

Similar to Example 4, we generate the corresponding Lyapunov equation describing system controllability of the IEEE 16-machine system. A ∈ $\mathbb{R}$^35×35 and C ∈ $\mathbb{R}$^35×35 of the Lyapunov equation are detailed in Supplementary Material. We set $\text{γ}=100$. The simulation results are shown in Table 5 and Figure 7. The definitions of all nouns in Figure 7 are the same as those in Figure 3.

FIGURE 7

FIGURE 7. Comparison of the convergence between the original-order and the reduced-order ZNN models under Example 5.

Example 6

Similar to Example 4, we generate the corresponding Lyapunov equation describing system controllability of the IEEE 48-machine system. A ∈ $\mathbb{R}$^97×97 and C ∈ $\mathbb{R}$^97×97 of the Lyapunov equation are detailed in Supplementary Material. We set $\text{γ}=100$. The simulation results are shown in Table 5 and Figure 8. The definitions of all nouns in Figure 8 are the same as those in Figure 3.

FIGURE 8

FIGURE 8. Comparison of the convergence between the original-order and the reduced-order ZNN models under Example 6.

It can be seen from Table 5 and Figures 6–8 that the neural dynamics method used to solve Lyapunov equations is also suitable for solving Lyapunov equations in power systems, and the reduced-order RNN mosdels proposed in this article is effective in the example of power systems. Moreover, with the increase in the power system scale, the convergence and steady-state accuracy of ZNN model are almost unchanged, indicating the applicability of the RNN model to power systems of different scales.

It is worth mentioning that the integration between the electric power and natural gas systems has been steadily enhanced in recent decades. The incorporation of natural gas systems brings, in addition to a cleaner energy source, greater reliability and flexibility to the power system (Liu et al., 2021). Since the dynamic model of the electricity–gas coupled system can be expressed by differential-algebraic equations (Zhang, 2005; Yang, 2020), which means the dynamic model of the electricity–gas coupled system is the same as that of the power system, the aforementioned applicability analysis of the methods proposed in this article for large power systems are also applicable to large electricity–gas coupled systems.

An Efficient Method for Vectorization of RNN Model

Table 6 compares the time cost of the RNN model vectorization method proposed in this article and the traditional RNN model vectorization method. For the sake of convenience, the former is called method A and the latter is called method B. In order to better demonstrate the effect of the vectorization method of RNN model proposed in this article, four examples are added, as shown in Table 6. Four newly added examples are generated in the same way as Example 3 and are detailed in Supplementary Material, where ms means millisecond; scale means the order of $A$; and proportion refers to the time taken by method A divided by the time taken by method B.

TABLE 6

TABLE 6. Comparison of the time cost of two vectorization methods for RNN models.

It can be seen from Table 6 that method A is significantly better than method B in terms of time cost, with the decrease in time cost between 48 and 98%. With the increase in the size of the examples, the proportion of time cost improvement generally increases. It should be pointed out that when the system sizes are 15, 35, and 97, the corresponding examples are the IEEE standard systems mentioned before, which indicates that the vectorization method proposed in this article is also effective in the example of power systems.

Conclusion

1) We propose an efficient method for vectorizing RNN models, which can achieve higher computational efficiency than the traditional method of vectorizing RNN based on the Kronecker product.

2) In order to reduce the solving scale of the RNN model, a reduced-order RNN model for solving the Lyapunov equation was proposed based on symmetry. At the same time, it is proved theoretically that the proposed model can maintain the same solution as that of the original model, and it is also proved that the proposed model is suitable for both the ZNN model and GNN model under linear or non-linear activation functions.

3) Several simulation examples are given to verify the effectiveness and superiority of the proposed method, while three standard examples of power systems are given to verify that the neural dynamics method is suitable for solving the Lyapunov equation of power systems.

Because the neural dynamics method has parallel distribution characteristics and hardware implementation convenience, its convergence and computation time are not sensitive to the system scale. Considering the current development level and trend of the very large-scale integration (VLSI) chip and the ultra large-scale integration (ULSI) chip, the wide application of the neural dynamics method in large-scale systems is expected.

In addition, the research on the RNN model used to solve the Lyapunov equation is mainly in the stage of theoretical improvement and exploration, and there are few reports about hardware products. The hardware product design will be the main content of the next stage.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This work was supported in part by the Key-Area Research and Development Program of Guangdong Province (2019B111109001), the National Natural Science Foundation of China (51577071), and the Southern Power Grid Corporation’s Science and Technology Project (Project No. 037700KK52190015 (GDKJXM20198313)).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2022.796325/full#supplementary-material

References

Bartels, R. H., and Stewart, G. W. (1972). Solution of the Matrix Equation AX + XB = C [F4]. Commun. ACM 15 (9), 820–826. doi:10.1145/361573.361582

CrossRef Full Text | Google Scholar

Chen, Z., and Zhou, J. (2012). Introduction to Matrix Theory. Beijing: Beihang University Press.

Google Scholar

Chen, Q., Gong, C., Zhao, J., Wang, Y., and Zou, D. (2017). Application of Parallel Sparse System Direct Solver Library Super LU_MT in State Estimation. Automation Electric Power Syst. 41 (3), 83–88. doi:10.7500/AEPS20160607008

CrossRef Full Text | Google Scholar

Cheng, K., and Chen, X. (2017). Linear Algebra. Chongqing: Chongqing University Press.

Google Scholar

Chicca, E., Stefanini, F., Bartolozzi, C., and Indiveri, G. (2014). Neuromorphic Electronic Circuits for Building Autonomous Cognitive Systems. Proc. IEEE 102 (9), 1367–1388. doi:10.1109/JPROC.2014.2313954

CrossRef Full Text | Google Scholar

Hafiz, F., Awal, M. A., Queiroz, A. R. d., and Husain, I. (2020). Real-Time Stochastic Optimization of Energy Storage Management Using Deep Learning-Based Forecasts for Residential PV Applications. IEEE Trans. Ind. Applicat. 56 (3), 2216–2226. doi:10.1109/TIA.2020.2968534

CrossRef Full Text | Google Scholar

He, W., and Zhang, S. (2017). Control Design for Nonlinear Flexible Wings of a Robotic Aircraft. IEEE Trans. Contr. Syst. Technol. 25 (1), 351–357. doi:10.1109/TCST.2016.2536708

CrossRef Full Text | Google Scholar

He, W., Ouyang, Y., and Hong, J. (2017). Vibration Control of a Flexible Robotic Manipulator in the Presence of Input Deadzone. IEEE Trans. Ind. Inf. 13 (1), 48–59. doi:10.1109/TII.2016.2608739

CrossRef Full Text | Google Scholar

Horn, R. A., and Johnson, C. R. (1991). Topics in Matrix Analysis. Cambridge: Cambridge University Press.

Google Scholar

Lan, X. (2017). Research on Model Order Reduction Method and Predictive Control Algorithm of Grid Voltage Control System (Beijing: North China Electric Power University). [dissertation/master’s thesis].

Google Scholar

Le, X., Chen, S., Li, F., Yan, Z., and Xi, J. (2019). Distributed Neurodynamic Optimization for Energy Internet Management. IEEE Trans. Syst. Man. Cybern, Syst. 49 (8), 1624–1633. doi:10.1109/TSMC.2019.2898551

CrossRef Full Text | Google Scholar

Li, X., Yu, J., Li, S., Shao, Z., and Ni, L. (2019a). A Non-linear and Noise-Tolerant ZNN Model and its Application to Static and Time-Varying Matrix Square Root Finding. Neural Process. Lett. 50 (2), 1687–1703. doi:10.1007/s11063-018-9953-y

CrossRef Full Text | Google Scholar

Li, Y., Zhang, H., Liang, X., and Huang, B. (2019b). Event-Triggered-Based Distributed Cooperative Energy Management for Multienergy Systems. IEEE Trans. Ind. Inf. 15 (4), 2008–2022. doi:10.1109/TII.2018.2862436

CrossRef Full Text | Google Scholar

Li, Y., Gao, D. W., Gao, W., Zhang, H., and Zhou, J. (2020). Double-Mode Energy Management for Multi-Energy System via Distributed Dynamic Event-Triggered Newton-Raphson Algorithm. IEEE Trans. Smart Grid 11 (6), 5339–5356. doi:10.1109/TSG.2020.3005179

CrossRef Full Text | Google Scholar

Lin, Y., and Simoncini, V. (2013). Minimal Residual Methods for Large Scale Lyapunov Equations. Appl. Numer. Maths. 72, 52–71. doi:10.1016/j.apnum.2013.04.004

CrossRef Full Text | Google Scholar

Liu, J., Zhang, J., and Li, Q. (2020a). Upper and Lower Eigenvalue Summation Bounds of the Lyapunov Matrix Differential Equation and the Application in a Class Time-Varying Nonlinear System. Int. J. Control. 93 (5), 1115–1126. doi:10.1080/00207179.2018.1494389

CrossRef Full Text | Google Scholar

Liu, Z., Chen, Y., Song, Y., Wang, M., and Gao, S. (2020b). Batched Computation of Continuation Power Flow for Large Scale Grids Based on GPU Parallel Processing. Power Syst. Techn. 44 (3), 1041–1046. doi:10.13335/j.1000-3673.pst.2019.2050

CrossRef Full Text | Google Scholar

Liu, H., Shen, X., Guo, Q., and Sun, H. (2021). A Data-Driven Approach towards Fast Economic Dispatch in Electricity-Gas Coupled Systems Based on Artificial Neural Network. Appl. Energ. 286, 116480. doi:10.1016/j.apenergy.2021.116480

CrossRef Full Text | Google Scholar

Raković, S. V., and Lazar, M. (2014). The Minkowski-Lyapunov Equation for Linear Dynamics: Theoretical Foundations. Automatica 50 (8), 2015–2024. doi:10.1016/j.automatica.2014.05.023

CrossRef Full Text | Google Scholar

Shanmugam, L., and Joo, Y. H. (2021). Stability and Stabilization for T-S Fuzzy Large-Scale Interconnected Power System with Wind Farm via Sampled-Data Control. IEEE Trans. Syst. Man. Cybern, Syst. 51 (4), 2134–2144. doi:10.1109/TSMC.2020.2965577

CrossRef Full Text | Google Scholar

Stykel, T. (2008). Low-rank Iterative Methods for Projected Generalized Lyapunov Equations. Electron. Trans. Numer. Anal. 30 (1), 187–202. doi:10.1080/14689360802423530

CrossRef Full Text | Google Scholar

Xiao, L., and Liao, B. (2016). A Convergence-Accelerated Zhang Neural Network and its Solution Application to Lyapunov Equation. Neurocomputing 193, 213–218. doi:10.1016/j.neucom.2016.02.021

CrossRef Full Text | Google Scholar

Xiao, L., Zhang, Y., Hu, Z., and Dai, J. (2019). Performance Benefits of Robust Nonlinear Zeroing Neural Network for Finding Accurate Solution of Lyapunov Equation in Presence of Various Noises. IEEE Trans. Ind. Inf. 15 (9), 5161–5171. doi:10.1109/TII.2019.2900659

CrossRef Full Text | Google Scholar

Yang, J., Chen, C. S., Abreu-garcia, J. A. D., and Xu, Y. (1993). Model Reduction of Unstable Systems. Int. J. Syst. Sci. 24 (12), 2407–2414. doi:10.1080/00207729308949638

CrossRef Full Text | Google Scholar

Yang, H. (2020). Dynamic Modeling and Stability Studies of Integrated Energy System of Electric, Gas and Thermal on Multiple Time Scales (Hunan: Changsha University of Science & Technology). [dissertation/master’s thesis].

Yi, C., Chen, Y., and Lu, Z. (2011). Improved Gradient-Based Neural Networks for Online Solution of Lyapunov Matrix Equation. Inf. Process. Lett. 111 (16), 780–786. doi:10.1016/j.ipl.2011.05.010

CrossRef Full Text | Google Scholar

Yi, C., Chen, Y., and Lan, X. (2013). Comparison on Neural Solvers for the Lyapunov Matrix Equation with Stationary & Nonstationary Coefficients. Appl. Math. Model. 37 (4), 2495–2502. doi:10.1016/j.apm.2012.06.022

CrossRef Full Text | Google Scholar

Yunong Zhang, Z., and Danchi Jiang, M. T. J. (1995). A Recurrent Neural Network for Solving Sylvester Equation with Time-Varying Coefficients. New York: Academic Press.

Google Scholar

Zhang, Y., Jiang, D., and Wang, J. (2002). A Recurrent Neural Network for Solving Sylvester Equation with Time-Varying Coeffificients. IEEE Trans. Neural Netw. 13 (5), 1053–1063. doi:10.1109/TNN.2002.1031938

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Chen, K., Li, X., Yi, C., and Zhu, H. (2008). “Simulink Modeling and Comparison of Zhang Neural Networks and Gradient Neural Networks for Time-Varying Lyapunov Equation Solving,” in Proceedings of IEEE International Conference on Natural Computation. Jinan. IEEE, 521–525. doi:10.1109/ICNC.2008.47

CrossRef Full Text | Google Scholar

Zhang, Y. (2005). Study on the Methods for Analyzing Combined Gas and Electricity Networks (Beijing: China Electric Power Research Institute). [dissertation/master’s thesis].

Zhao, H., Lan, X., Xue, N., and Wang, B. (2014). Excitation Prediction Control of Multi‐machine Power Systems Using Balanced Reduced Model. IET Generation, Transm. Distribution 8 (6), 1075–1081. doi:10.1049/iet-gtd.2013.0609

CrossRef Full Text | Google Scholar

Zhou, B., Duan, G.-R., and Li, Z.-Y. (2009). Gradient Based Iterative Algorithm for Solving Coupled Matrix Equations. Syst. Control. Lett. 58 (5), 327–333. doi:10.1016/j.sysconle.2008.12.004

CrossRef Full Text | Google Scholar

Zhu, Z., Geng, G., and Jiang, Q. (2016). Power System Dynamic Model Reduction Based on Extended Krylov Subspace Method. IEEE Trans. Power Syst. 31 (6), 4483–4494. doi:10.1109/TPWRS.2015.2509481

CrossRef Full Text | Google Scholar