The search functionality is under construction.
The search functionality is under construction.

Open Access
A Small-Data Solution to Data-Driven Lyapunov Equations: Data Reduction from O(n2) to O(n)

Keitaro TSUJI, Shun-ichi AZUMA, Ikumi BANNO, Ryo ARIIZUMI, Toru ASAI, Jun-ichi IMURA

  • Full Text Views

    255

  • Cite this
  • Free PDF (3.5MB)

Summary :

When a mathematical model is not available for a dynamical system, it is reasonable to use a data-driven approach for analysis and control of the system. With this motivation, the authors have recently developed a data-driven solution to Lyapunov equations, which uses not the model but the data of several state trajectories of the system. However, the number of state trajectories to uniquely determine the solution is O(n2) for the dimension n of the system. This prevents us from applying the method to a case with a large n. Thus, this paper proposes a novel class of data-driven Lyapunov equations, which requires a smaller amount of data. Although the previous method constructs one scalar equation from one state trajectory, the proposed method constructs three scalar equations from any combination of two state trajectories. Based on this idea, we derive data-driven Lyapunov equations such that the number of state trajectories to uniquely determine the solution is O(n).

Publication
IEICE TRANSACTIONS on Fundamentals Vol.E107-A No.5 pp.806-812
Publication Date
2024/05/01
Publicized
2023/10/24
Online ISSN
1745-1337
DOI
10.1587/transfun.2023MAP0010
Type of Manuscript
Special Section PAPER (Special Section on Mathematical Systems Science and its Applications)
Category

1.  Introduction

When a mathematical model is not available for a dynamical system, it is reasonable to use a data-driven approach for analysis and control of the system. Thus, various data-driven solutions have been recently developed for several stability [1], [2], controllability [3], [4], observability [4], optimal control [5], [6], and so on.

Along this direction of research, the authors have recently developed a data-driven solution to Lyapunov equations in the form of \(PA+A^\top P= -Q\), where \(A\) and \(Q\) are given constant matrices and \(P\) is the unknown [7]. The Lyapunov equations are known to play an important role in control engineering, such as stability analysis and controllability analysis. The solution is based on the so-called data-driven Lyapunov equation, which is defined by not the matrix \(A\) but the state trajectories of the system \(\dot{x}(t)=Ax(t)+Bu(t)\). Using this framework, one can solve a Lyapunov equation even when a model is not available to the system.

On the other hand, the number of state trajectories to uniquely determine the solution is \(O(n^2)\) in [7], where \(n\) is the size of the matrix \(A\). In other words, the amount of data quadratically grows with \(n\). This prevents us from applying this method to a large-scale system, and thus it is preferable to reduce the amount of data.

This paper proposes a novel class of data-driven Lyapunov equations such that the number of state trajectories to uniquely determine the solution is \(O(n)\). In the previous method [7], one scalar equation, which corresponds to a part of the Lyapunov equation, is formulated by one state trajectory. On the other hand, in the proposed method, three scalar equations are formulated by any combination of two state trajectories. This idea drastically reduces the amount of data to construct a data-driven Lyapunov equation.

Finally, we note that our data-driven approach has an advantage that prior knowledge can be easily incorporated. In fact, our approach transforms a Lyapunov equation into an equation described by data, which allows us to take prior knowledge into account as additional equations. On the other hand, the model-based approach, which is based on a mathematical model, uniquely provides the solution to a Lyapunov equation from the model. Thus, it is hard to directly incorporate prior knowledge into the problem of solving a Lyapunov equation.

This paper is organized as follows. In Sect. 2, we introduce our problem and review the existing method and the amount of required data [7]. Our solutions are presented in Sects. 3 and 4, respectively. Finally, Sect. 5 concludes this paper.

Notation: (i) For a vector \(x \in \mathbb{R}^n\) and symmetric matrix \(P \in \mathbb{R}^{n \times n}\), let \(\mathrm{lvec}(x) \in \mathbb{R}^{1 \times \frac{1}{2}n(n+1)}\) and \(\mathrm{rvec}(P) \in \mathbb{R}^{\frac{1}{2}n(n+1)}\) be the pair of a row vector and column vector, such that \(x^\top P x = \mathrm{lvec}(x)\mathrm{rvec}(P)\). For example,

\[\begin{aligned} \mathrm{lvec}(x) &= \begin{bmatrix} x_1^2 & 2x_{1}x_{2} & 2x_{1}x_{3} & x_2^2 & 2x_{2}x_{3} & x_3^2 \end{bmatrix} , \\ \mathrm{rvec}(P) &= \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{22} & p_{23} & p_{33} \end{bmatrix}^\top \end{aligned}\]

for

\[x = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} , P = \begin{bmatrix} p_{11} & p_{12} & p_{13} \\ p_{12} & p_{22} & p_{23} \\ p_{13} & p_{23} & p_{33} \\ \end{bmatrix} ,\]

where \(x_i \in \mathbb{R}\) is the \(i\)-th component of \(x\), and \(p_{ij} \in \mathbb{R}\) is the \((i,j)\)-th component of \(P\).

(ii) Consider the equations \(f(x)=0\) and \(g(x) =0\), where \(x \in \mathbb{R}^n\). If their solution sets are equal, then the equations are said to be equivalent.

2.  Problem Formulation and Existing Result

Consider the linear system

\[\begin{equation*} \dot{x}(t) = Ax(t) , \tag{1} \end{equation*}\]

where \(x(t) \in \mathbb{R}^n\) is the state and \(A \in \mathbb{R}^{n \times n}\) is a Hurwitz matrix.

For the system, a Lyapunov equation is given by

\[\begin{equation*} PA + A^\top P = -Q , \tag{2} \end{equation*}\]

where \(Q \in \mathbb{R}^{n \times n}\) is a given symmetric matrix and \(P \in \mathbb{R}^{n \times n}\) is the unknown. There exists a unique solution \(P \in \mathbb{R}^{n \times n}\) to (2) if \(A\) is Hurwitz [8].

If the matrix \(A\) is known to us, the solution can be easily obtained; otherwise, it is reasonable to use a data-driven solution to the equation. The problem of finding the solution \(P\) from the behavioral data of (1) is formulated as follows.

Problem 1: Consider the system in (1) and the Lyapunov equation in (2). Assume that \(A\) is unknown but Hurwitz. Suppose that the data of state trajectories \(x(t,x_{0i}),\; t \in [0,T]\;(i=1,2,\ldots,q)\) are given, where \(T\) is a positive real number and \(x(t,x_{0i})\) is the state \(x(t)\) of the system (1) for the initial state \(x(0) = x_{0i}\). Find the solution of (2). \(\tag*{◻}\)

For Problem 1, the following solution [7] has been presented. Let us briefly review the existing result. Consider the Lyapunov equation in (2). Let us multiply the both sides of (2) by \(x(t,x_{0i})^\top\) from the left and by \(x(t,x_{0i})\) from the right and integrate them on the interval \([0,T]\). As the result, we obtain the following equation.

\[\begin{multline*} x^\top(T,x_{0i}) P x(T,x_{0i}) - x_{0i}^\top P x_{0i} = -\int_0^T \!\!\! x^\top(t,x_{0i}) Q x(t,x_{0i})dt\\ (i=1,2,\ldots,q) \tag{3} \end{multline*}\]

This is called the data-driven Lyapunov equation.

By using the functions \(\mathrm{lvec}\) and \(\mathrm{rvec}\) introduced at the end of Sect. 1, (3) is equivalently transformed into

\[\begin{multline*} [\mathrm{lvec}(x(T,x_{0i})) - \mathrm{lvec}(x_{0i})]\mathrm{rvec}(P) \\ = -\int_0^T x^\top(t,x_{0i}) Q x(t,x_{0i})dt\\ (i=1,2,\ldots,q) , \tag{4} \end{multline*}\]

which is represented as

\[\begin{equation*} D\,\mathrm{rvec}(P) = d \tag{5} \end{equation*}\]

for

\[\begin{align} D &:= \begin{bmatrix} \mathrm{lvec}(x(T,x_{01})) - \mathrm{lvec}(x_{01})\\ \mathrm{lvec}(x(T,x_{02})) - \mathrm{lvec}(x_{02})\\ \vdots\\ \mathrm{lvec}(x(T,x_{0q})) - \mathrm{lvec}(x_{0q})\\ \end{bmatrix} \in \mathbb{R}^{q \times \frac{1}{2}n(n+1)} , \tag{6} \\ d &:= \begin{bmatrix} -\int_0^T x^\top(t,x_{01}) Q x(t,x_{01})dt \\ -\int_0^T x^\top(t,x_{02}) Q x(t,x_{02})dt\\ \vdots\\ -\int_0^T x^\top(t,x_{0q}) Q x(t,x_{0q})dt \end{bmatrix} \in \mathbb{R}^{q}. \tag{7} \end{align}\]

Equation (5) is in the standard form of a linear equation with respect to the unknown vector \(\mathrm{rvec}(P) \in \mathrm{R}^{\frac{1}{2}n(n+1)}\). Thus, if

\[\begin{equation*} \mathrm{rank}(D) = \frac{1}{2}n(n+1) , \tag{8} \end{equation*}\]

then (5) has a unique solution. Based on this idea, the solution to Problem 1 in [7] is given as follows.

Lemma 1: Consider the equation in (3). If (8) holds, then (3) has a unique solution, and it is equal to the solution to the Lyapunov equation in (2). \(\tag*{◻}\)

From this result, it turns out that the solution to the Lyapunov equation in (2) is obtained by solving a linear equation in (3) defined by the data of state trajectories. Since \(D \in \mathrm{R}^{q \times \frac{1}{2}n(n+1)}\) and (8) imply \(q \geq (1/2)n(n+1)\), the data-driven method [7] requires at least \((1/2)n(n+1)\) state trajectories, i.e., it requires data of \(O(n^2)\).

Remark: In Problem 1, \(A\) is assumed to be Hurwitz. Thus, the problem is not solved to determine whether the system in (1) is stable or not. One application of Problem 1 is to analyze the observability of (1), since the solution of (2) for \(Q=C^{\top}C\) is equal to the observability Gramian of (1) [8]. Another application is to obtain a Lyapunov function of a stable system, which can be used, for example, to construct an event-triggered controller [9], [10].

3.  Small-Data Solution to Problem 1

3.1  Main Results

In the data-driven Lyapunov equation in (3), each scalar equation is constructed by one state trajectory. However, as will be shown later, three scalar equations can be formulated by any combination of two distinct state trajectories. This fact gives the following result.

Theorem 1: (i) Consider the following equation:

\[\begin{multline*} x^\top(T,x_{0i}) P x(T,x_{0j}) - x_{0i}^\top P x_{0j} = -\int_0^T \!\!\!\!\! x^\top(t,x_{0i}) Q x(t,x_{0j})dt\\ (i=1,2,\ldots,q,\;j=i,i+1,\ldots,q) \tag{9} \end{multline*}\]

which is a modified version of (3). If

\[\begin{equation*} \mathrm{rank}(X_0) = n , \tag{10} \end{equation*}\]

then (9) has a unique solution, and it is equal to the solution to (2), where \(X_0 := \begin{bmatrix} x_{01} & x_{02} & \cdots & x_{0q} \end{bmatrix}\) is the matrix defined by the initial states of the data \(x(t,x_{0i}),\;t \in [0,T]\;(i=1,2,\ldots,q)\).

(ii) The minimum number of state trajectories to determine the unique solution to (9) is \(n\), i.e., \(O(n)\). \(\tag*{◻}\)

This means that the solution to the Lyapunov equation can be obtained with much less data than the conventional method [7] (Lemma 1). This fact is summarized in Table 1.

Table 1  Proposed and conventional methods with \(q\) state trajectories.

Proof of Theorem 1: (i) It is the consequence of the following facts.

  1. Equation (2) has a unique solution if \(A\) is Hurwitz [8].
  2. Equation (2) is equivalent to
    \[\begin{equation*} e^{A^\top T} P e^{AT} - P = -\int_0^T e^{A^\top t} Q e^{At}dt . \tag{11} \end{equation*}\]
  3. If (10) holds, then (11) is equivalent to
    \[\begin{equation*} X_0^\top e^{A^\top T} P e^{AT} X_0 - X_0^\top P X_0 = -\int_0^T X_0^\top e^{A^\top t} Q e^{At} X_0dt . \tag{12} \end{equation*}\]
  4. Equation (12) is equivalent to (9).

Let us prove Facts (II)-(IV).

(II) First, we multiply the both sides of (2) by \(e^{A^\top t}\) from the left and by \(e^{A t}\) from the right and integrate them on the interval \([0,T]\). As the result, we obtain (11). Thus, if \(P\) is a solution to (2), then it is also a solution to (11). Meanwhile, the converse holds if both (2) and (11) are unique solutions. Since (11) is a discrete-time Lyapunov equation with the Schur matrix \(e^{AT}\) (because \(A\) is Hurwitz), it has a unique solution. These facts and (I) imply (II).

(III) Under (10), the matrix \(X_0\) has full row rank and \(X_0^\top\) has full column rank, which, together with Lemma 2 in the appendix, proves (III).

(IV) It is trivial from the facts that (9) is the element-wise expression of (12) and \(x(t,x_{0i})=e^{At}x_{0i}\).

(ii) Equation (10) implies \(q \geq n\). This proves (ii). \(\tag*{◻}\)

3.2  Examples
3.2.1  Case with Noise-Free data

Consider the linear system

\[\begin{equation*} \dot{x}(t)= \begin{bmatrix} 0 & 1\\ -2 & -3\\ \end{bmatrix} x(t) \tag{13} \end{equation*}\]

and the Lyapunov equation in (2) for \(Q := I\), which has the unique solution

\[\begin{equation*} P^*= \begin{bmatrix} 1.250 & 0.250\\ 0.250 & 0.250\\ \end{bmatrix}. \tag{14} \end{equation*}\]

Then we address Problem 1 for the dataset shown in Fig. 1, where the state trajectories are generated by \(x_{01} = [1 0]^\top\) and \(x_{02} = [0 1]^\top\).

Fig. 1  Dataset in Sect. 3.2.1.

Let us apply Theorem 1 to the problem. In this case, we have \(q=n=2\) and (10) holds for \(x_{01}\) and \(x_{02}\). Then the data-driven Lyapunov equation in (9) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 0.600 & -0.465 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.600 \\ -0.465 \end{bmatrix}\\ & - \begin{bmatrix} 1.000 & 0.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 0.000 \end{bmatrix} = -0.885 , \tag{15} \\ &\nonumber \begin{bmatrix} 0.233 & -0.097 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.600 \\ -0.465 \end{bmatrix}\\ & - \begin{bmatrix} 0.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 0.000 \end{bmatrix} = -0.106 , \tag{16} \\ &\nonumber \begin{bmatrix} 0.233 & -0.097 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.233 \\ -0.097 \end{bmatrix}\\ & - \begin{bmatrix} 0.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.000 \\ 1.000 \end{bmatrix} = -0.191 , \tag{17} \end{align}\]

which provides

\[\begin{equation*} P= \begin{bmatrix} 1.250 & 0.250\\ 0.250 & 0.250\\ \end{bmatrix} . \tag{18} \end{equation*}\]

This agrees with the solution to the original Lyapunov equation, i.e., \(P^*\) in (14).

3.2.2  Case with Noisy Data

We consider the case where noisy data is given as shown in Fig. 2. The data are the state trajectories in Fig. 1 corrupted by Gaussian noise with a mean of 0 and a variance of 0.05.

Fig. 2  Noisy dataset.

The data-driven Lyapunov equation in (9) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 0.551 & -0.444 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.551 \\ -0.444 \end{bmatrix}\\ & - \begin{bmatrix} 0.943 & 0.036 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.943 \\ 0.036 \end{bmatrix} = -0.878 , \tag{19} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.551 \\ -0.444 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.943 \\ 0.036 \end{bmatrix} = -0.110 , \tag{20} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.227 \\ -0.098 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} -0.049 \\ 0.960 \end{bmatrix} = -0.199 . \tag{21} \end{align}\]

which provides

\[\begin{equation*} P= \begin{bmatrix} 1.300 & 0.317\\ 0.317 & 0.306\\ \end{bmatrix} \tag{22} \end{equation*}\]

It turns out that this approximates the true value \(P^*\) in (14).

On the other hand, the existing method [7] requires more data, i.e., three state trajectories, because three scalar equations are constructed by three state trajectories, as shown in (3). For example, for the three state trajectories in Figs. 2 and 3, the data-driven Lyapunov equation in (3) is given by (19), (21), and

\[\begin{align} &\nonumber \begin{bmatrix} 1.036 & -0.683 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ &- \begin{bmatrix} 0.969 & 1.977 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -2.066, \tag{23} \end{align}\]

which provides

\[\begin{equation*} P = \begin{bmatrix} 1.384 & 0.229\\ 0.229 & 0.305 \end{bmatrix}. \tag{24} \end{equation*}\]

Although the solution is derived by more data, the accuracy is similar to the proposed method.

Finally, we show that the proposed method provides better performance with the data as shown in Figs. 2 and 3. In this case, the data-driven Lyapunov equation is given by (19), (20), (21), (23),

\[\begin{align} &\nonumber \begin{bmatrix} 0.551 & -0.444 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ & - \begin{bmatrix} 0.943 & 0.036 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -1.085 , \tag{25} \\ &\nonumber \begin{bmatrix} 0.227 & -0.098 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.036 \\ -0.683 \end{bmatrix}\\ & - \begin{bmatrix} -0.049 & 0.960 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 0.969 \\ 1.977 \end{bmatrix} = -0.494 . \tag{26} \end{align}\]

The least-squares solution of these equations is

\[\begin{equation*} P= \begin{bmatrix} 1.289 & 0.252\\ 0.252 & 0.278\\ \end{bmatrix} , \tag{27} \end{equation*}\]

which is closer to \(P^*\) than (24).

Fig. 3  Additional noisy data.

3.2.3  Application to Input Node Design of a Large-Scale System

Consider the model of a 100-order gene regulatory network system as shown in Fig. 4, which is developed in [11]. For the system, we address the problem of finding a gene to which the control input is applied so that the trace of the controllability Gramian is maximized.

Fig. 4  Network structure.

As shown in [7], the solution is the index of the maximum diagonal element of the solution \(P\) of (2) for \(Q = I\). Thus, let us solve the corresponding data-driven Lyapunov equation. Our result provides gene 34 as the solution with the data of 100 state trajectories. On the other hand, the existing method [7] requires 5050 state trajectories to obtain the solution. This suggests that our method is more practical than the existing method.

4.  Extension to the Case with Discrete-Time Data

4.1  Main Results

Suppose that the state trajectory \(x(t,x_{0})\) is sampled for the period \(h\) and the discrete data \(\bar{x}_k := x(kh,x_{0})\) is given. Then we consider the discrete-time Lyapunov equation

\[\begin{equation*} \bar{A}^\top P \bar{A} - P = -Q , \tag{28} \end{equation*}\]

where \(\bar{A}:= e^{Ah}\), \(Q \in \mathbb{R}^{n \times n}\) is a given symmetric matrix, and \(P \in \mathbb{R}^{n \times n}\) is the unknown. There exists a unique solution \(P \in \mathbb{R}^{n \times n}\) to (28) if \(A\) is Hurwitz, i.e., \(\bar{A}\) is Schur [8].

Now, we consider the following problem.

Problem 2: Consider the system in (1) and the discrete-time Lyapunov equation in (28). Assume that \(\bar{A}\) is unknown and Schur, but the sampled data \(\bar{x}_k := x(kh,x_{0})\;(k=0,1,\dots,s)\) of the state trajectories of (1) are given, where \(h \in (0,\infty)\) is the sampling period. Find the solution of (28). \(\tag*{◻}\)

To derive a solution, we represent the data \(\bar{x}_k\) \((k=0,1,\dots,s-1)\) in the following matrix form:

\[\begin{equation*} \bar{X}_0 := \begin{bmatrix} \bar{x}_0 & \bar{x}_1 & \cdots & \bar{x}_{s-1} \end{bmatrix}. \tag{29} \end{equation*}\]

Then a solution to Problem 2 is obtained as follows.

Theorem 2: (i) Consider the following equation:

\[\begin{multline*} \bar{x}_k^\top P \bar{x}_l - \bar{x}_{k-1}^\top P \bar{x}_{l-1} = -\bar{x}_{k-1}^\top Q \bar{x}_{l-1}\\ (k=1,2,\ldots,s,\;l=k,k+1,\ldots,s). \tag{30} \end{multline*}\]

If

\[\begin{equation*} \mathrm{rank}(\bar{X}_0) = n, \tag{31} \end{equation*}\]

then (30) has a unique solution, and it is equal to the solution to (28).

(ii) The minimum number of the sampled data to determine the unique solution to (30) is \(n+1\), i.e., \(O(n)\). \(\tag*{◻}\)

Proof of Theorem 2: (i) It is the consequence of the following facts.

  1. Equation (28) has a unique solution if \(A\) is Hurwitz [8].
  2. If (31) holds, then (28) is equivalent to
    \[\begin{equation*} \bar{X}_0^\top \bar{A}^\top P \bar{A} \bar{X}_0 - \bar{X}_0^\top P \bar{X}_0 = -\bar{X}_0^\top Q \bar{X}_0 . \tag{32} \end{equation*}\]
  3. Equation (32) is equivalent to (30).

Let us prove Facts (II) and (III).

(II) Under (31), the matrix \(\bar{X}_0\) has full row rank and \(\bar{X}_0^\top\) has full column rank, which, together Lemma 2 in the appendix, proves (II).

(III) It is trivial from the facts that (30) is the element-wise expression of (32) and \(\bar{x}_k =\bar{A}\bar{x}_{k-1}\).

(ii) Equation (31) implies \(s \geq n\). This proves (ii). \(\tag*{◻}\)

4.2  Example

Consider the linear system in (13) and the discrete-time Lyapunov equation in (28) for \(Q := I\) and \(h=0.1\). The solution is given by

\[\begin{equation*} P^*= \begin{bmatrix} 13.000 & 2.508 \\ 2.508 & 3.050 \end{bmatrix} . \tag{33} \end{equation*}\]

Then we consider Problem 2 for the dataset composed of

\[\begin{equation*} \bar{x}_{0} = \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix},\, \bar{x}_{1} = \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix},\, \bar{x}_{2} = \begin{bmatrix} 1.116 \\ 0.225 \end{bmatrix}. \tag{34} \end{equation*}\]

Let us apply Theorem 2 to the problem. In this case, (31) holds for \(\bar{x}_{0}\) and \(\bar{x}_{1}\). Then the data-driven discrete-time Lyapunov equation in (30) is given by

\[\begin{align} &\nonumber \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix}\\ & - \begin{bmatrix} 1.000 & 1.000 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix} = -2.000 , \tag{35} \\ &\nonumber \begin{bmatrix} 1.116 & 0.225 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix}\\ & - \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.000 \\ 1.000 \end{bmatrix} = -1.637 , \tag{36} \\ &\nonumber \begin{bmatrix} 1.116 & 0.225 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.116 \\ -0.225 \end{bmatrix}\\ & - \begin{bmatrix} 1.077 & 0.560 \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{bmatrix} \begin{bmatrix} 1.077 \\ 0.560 \end{bmatrix} = -1.474 , \tag{37} \end{align}\]

which provides

\[\begin{equation*} P= \begin{bmatrix} 13.000 & 2.508 \\ 2.508 & 3.050 \end{bmatrix} . \tag{38} \end{equation*}\]

This agrees with the solution to the original discrete-time Lyapunov equation, i.e., \(P^*\) in (33).

5.  Conclusions

In this paper, we have developed a new class of data-driven Lyapunov equations, which are defined by data of amount \(O(n)\). Since the amount of data is much less than that of the existing result [7], i.e., \(O(n^2)\), this result broadens the applicability of the framework of data-driven Lyapunov equations.

On the other hand, it is still open whether \(O(n)\) is the minimum or not to determine the unique solution to a Lyapunov equation in a data-driven manner. In the future, we plan to answer to this question.

Acknowledgments

This research was supported by JST Moonshot R&D #JPMJMS2021.

References

[1] U.S. Park and M. Ikeda, “Stability analysis and control design of lti discrete-time systems by the direct use of time series data,” Automatica, vol.45, no.5, pp.1265-1271, 2009.
CrossRef

[2] F. Zhang and D. Söffker, “A data-driven online stability monitoring method for unknown discrete-time nonlinear systems,” 50th IEEE Conference on Decision and Control and European Control Conference, pp.7356-7361, 2011.
CrossRef

[3] V.K. Mishra, I. Markovsky, and B. Grossmann, “Data-driven tests for controllability,” IEEE Control Syst. Lett., vol.5, no.2, pp.517-522, 2021.
CrossRef

[4] Z. Wang and D. Liu, “Data-based controllability and observability analysis of linear discrete-time systems,” IEEE Trans. Neural Netw., vol.22, no.12, pp.2388-2392, 2011.
CrossRef

[5] D. Vrabie, O. Pastravanu, M. Abu-Khalaf, and F. Lewis, “Adaptive optimal control for continuous-time linear systems based on policy iteration,” Automatica, vol.45, no.2, pp.477-484, 2009.
CrossRef

[6] Y. Jiang and Z.P. Jiang, “Computational adaptive optimal control for continuous-time linear systems with completely unknown dynamics,” Automatica, vol.48, no.10, pp.2699-2704, 2012.
CrossRef

[7] I. Banno, S. Azuma, R. Ariizumi, T. Asai, and J. Imura, “Data-driven estimation and maximization of controllability Gramians,” 2021 60th IEEE Conference on Decision and Control, pp.5053-5058, 2021.
CrossRef

[8] B.N. Datta, Numerical Methods for Linear Control Systems: Design and Analysis, Academic Press, 2004.
CrossRef

[9] I. Banno, S. Azuma, R. Ariizumi, and T. Asai, “Data-driven sparse event-triggered control of unknown systems,” 2021 American Control Conference, pp.3383-3388, 2021.
CrossRef

[10] W.P.M.H. Heemels, K.H. Johansson, and P. Tabuada, “An introduction to event-triggered and self-triggered control,” IEEE Conference on Decision and Control, pp.3270-3285, 2012.
CrossRef

[11] X. Shen, M. Morishita, J. Imura, M. Oku, and K. Aihara, “Low-sample-size data-driven re-stabilization of gene network systems,” 10th IFAC Symposium on Robust Control Design ROCOND 2022, vol.55, no.25, pp.241-246, 2022.
CrossRef

[12] R. Piziak and P.L. Odell, “Full rank factorization of matrices,” Mathematics Magazine, vol.72, no.3, pp.193-201, 1999.
CrossRef

Appendix: On Equivalence of Two Equations

For a matrix-valued function \(f: \mathbb{R}^{n \times n} \to \mathbb{R}^{n \times n}\) and matrices \(A \in \mathbb{R}^{m \times n}\), \(B \in \mathbb{R}^{n \times l}\), we consider the two equations

\[\begin{align} &f(X) = 0_{n \times n} , \tag{A$\cdot $1} \\ &Af(X)B = 0_{m \times l} , \tag{A$\cdot $2} \end{align}\]

where \(0_{n \times n}\) is an \(n \times n\) zero matrix and \(0_{m \times l}\) is similarly defined. Let \(S \subset \mathbb{R}^{n \times n}\) be the set of the solutions to (A\(\cdot\)1) and \(\tilde{S} \subset \mathbb{R}^{n \times n}\) be that to (A\(\cdot\)2), i.e.,

\[\begin{align} S &:= \{ X \in \mathbb{R}^{n \times n} \mid f(X) = 0_{n \times n} \} , \tag{A$\cdot $3} \\ \tilde{S} &:= \{ X \in \mathbb{R}^{n \times n} \mid Af(X)B = 0_{m \times l} \} . \tag{A$\cdot $4} \end{align}\]

Lemma 2: Consider the sets in (A\(\cdot\)3) and (A\(\cdot\)4). If \(A\) has full column rank and \(B\) has full row rank, then \(S=\tilde{S}\). \(\tag*{◻}\)

Proof: The relation \(S \subset \tilde{S}\) is trivial from \(X \in S \Rightarrow f(X) = 0_{n \times n} \Rightarrow Af(X)B = 0_{m \times l} \Rightarrow X \in \tilde{S}\).

Next, we prove \(S \supset \tilde{S}\). If \(X \in \tilde{S}\), then

\[\begin{equation*} Af(X)B = 0_{m \times l} . \tag{A$\cdot $5} \end{equation*}\]

If \(A\) has full column rank, then \(A^\top A\) is non-singular [12]. Similarly, if \(B\) has full row rank, then \(B B^\top\) is non-singular. Therefore, multiplying the both sides of (A\(\cdot\)5) by \((A^\top A)^{-1}A^\top\) from the left and by \(B^\top(B B^\top)^{-1}\) from the right, we get

\[\begin{equation*} (A^\top A)^{-1}A^\top Af(X)B B^\top(B B^\top)^{-1} = 0_{n \times n} , \tag{A$\cdot $6} \end{equation*}\]

which is equal to

\[\begin{equation*} f(X) = 0_{n \times n} , \tag{A$\cdot $7} \end{equation*}\]

i.e., \(X \in S\). This proves \(S \supset \tilde{S}\).

Since \(S \subset \tilde{S}\) and \(S \supset \tilde{S}\), we obtain \(S=\tilde{S}\). \(\tag*{◻}\)

Authors

Keitaro TSUJI
  Nagoya University

received his B.E. degrees in engineering from Nagoya University, Nagoya, Japan in 2022. He is now M.E. student in engineering at Nagoya University.

Shun-ichi AZUMA
  Kyoto University

received the Ph.D. degrees in control engineering from Tokyo Institute of Technology in 2004. He was a research fellow of the Japan Society for the Promotion of Science from 2004 to 2005, an Assistant Professor and an Associate Professor in the Graduate School of Informatics, Kyoto University from 2005 to 2011 and from 2011 to 2017, respectively, and a Professor in the Graduate School of Engineering, Nagoya University. He is currently a Professor in the Graduate School of Informatics, Kyoto University. He serves/served as Associate Editors of several journals, such as IEEE Transactions on Automatic Control, IEEE Transactions on Control of Network Systems, IEEE Control Systems Letters, and Automatica. His research interests include analysis and control of network systems and hybrid systems.

Ikumi BANNO
  Nagoya University

received his B.E. and M.E. degrees in engineering from Nagoya University, Nagoya, Japan in 2019 and 2021, respectively. He is now Ph.D. student in engineering at Nagoya University.

Ryo ARIIZUMI
  Nagoya University

received his B.E., M.E., and Ph.D. degrees from Kyoto University, Kyoto, Japan, in 2010, 2012, and 2015, respectively. He was a research fellow of the Japan Society for the Promotion of Science from 2014 to 2015. He is currently an Assistant Professor at Nagoya University, Nagoya, Japan. His research interests include the control of redundant robots and the optimization of robotic systems.

Toru ASAI
  Nagoya University

received his B.E., M.E. and Ph.D. degrees from Tokyo Institute of Technology, Japan, in 1991, 1993 and 1996, respectively. He has worked as a Research Fellow of JSPS between 1996 and 1998. In 1999, he joined the faculty of Osaka University. He is an Associate Professor of the sub-department of Mechatronics between 2015 and 2016 and the Department of Mechanical Systems Engineering since 2017, Nagoya University. His research interests include robust control, switching control, parameter estimation and industrial applications. He is a member of IEEE, SICE, ISIJ and ISCIE. He is a member of IFAC TC 9.4 Control Education.

Jun-ichi IMURA
  Tokyo Institute of Technology

received the M.E. degree in applied systems science and the Ph.D. degree in mechanical engineering from Kyoto University, Kyoto, Japan, in 1990 and 1995, re- spectively. He served as a Research Associate at the Department of Mechanical Engineering, Kyoto University, from 1992 to 1996, and as an Associate Professor in the Division of Machine Design Engineering, Faculty of Engineering, Hi- roshima University, Hiroshima, Japan, from 1996 to 2001. From May 1998 to April 1999, he was a Visiting Researcher at the Faculty of Mathematical Sciences, University of Twente, Enschede, The Netherlands. Since 2001, he has been with Tokyo Institute of Technology, Tokyo, Japan, where he is currently a Professor at the Department of Systems and Control Engineering. His research interests include modeling, analysis, and synthesis of nonlinear systems, hybrid systems, and large-scale network systems with applications to power systems, ITS, biological systems, and indus- trial process systems. He is a member of the Society of Instrument and Control Engineers (SICE), The Institute of Systems, Control and Information Engineers (ISCIE), and The Robotics Society of Japan.

Keyword