Data-driven Estimation of the Algebraic Riccati Equation for
the Discrete-Time Inverse Linear Quadratic Regulator Problem

The inverse optimal control (IOC) problem is a framework used to determine the objective function optimized by a given control system. By optimizing the objective function obtained by IOC, a different system can imitate the given control system’s behavior. For example, several studies have been conducted on robots imitating human or animal movements [1, 2, 3]. To apply IOC, in these studies, the researchers assumed that human and animal movements are optimal for unknown criteria. Such an optimality assumption holds in some cases [4].

The first IOC problem, which assumes a single-input linear time-invariant system and quadratic objective function, was proposed by Kalman [5]. Anderson [6] generalized the IOC problem in [5] to the multi-input case. For this type of IOC problem, called the inverse linear quadratic regulator (ILQR) problem, Molinari [7] proposed a necessary and sufficient condition for a linear feedback input to optimize some objective functions. Moylan et al. [8] proposed and solved the IOC problem for nonlinear systems. IOC in [5, 6, 7, 8] is based on control theory, whereas Ng et al. [9] proposed inverse reinforcement learning (IRL), which is IOC based on machine learning. Recently, IRL has become an important IOC framework, along with control theory methods [10].

We consider the ILQR problem with respect to an unknown discrete-time system. Such a problem often occurs in biological system analysis, which is the main application target of IOC. The control theory approach solves the ILQR problem by solving a linear matrix inequality (LMI) that contains the algebraic Riccati equation (ARE) [11]. However, to obtain the ARE, the system model is required. The IRL approach also has difficulty solving our problem. There is IRL with unknown systems [12], and IRL with continuous states and input spaces [13]; however, to the best of our knowledge, no IRL exists with both.

In the continuous-time case, we proposed a method to estimate the ARE from the observation data of the system state and input directly [14]. In the present paper, we estimate the ARE for a discrete-time system using a method that extends the result in [14]. Similarly to [14], our method transforms the ARE by multiplying the observed state and input on both sides. However, this technique alone cannot calculate ARE without the system model because the form of the ARE differs between continuous and discrete-time. We solve this problem using inputs from the system’s controller. Also, the use of such inputs enable us to estimate the ARE without knowing the system’s control gain. We prove that the equation obtained from this transformation is equivalent to the ARE if the dataset is appropriate. The advantage of our method is the economization of the observation data using prior information about the objective function. We conducted a numerical experiment to demonstrate that our method can estimate the ARE with less data than system identification if the prior information is sufficient.

The structure of the remainder of this paper is as follows: In Section II, we formulate the problem considered in this paper. In Section III, we propose our estimation method and prove that the estimated equation is equivalent to the ARE. In Section IV, we describe numerical experiments that confirm our statement in Section III. Section V concludes the paper.

We consider the following discrete-time linear system:

where $A\in{\mathbb{R}}^{n\times n}$ and $B\in{\mathbb{R}}^{n\times m}$ are constant matrices, and $x{\left(k\right)}\in{\mathbb{R}}^{n}$ and $u{\left(k\right)}\in{\mathbb{R}}^{m}$ are the state and the input at time $k\in\mathbb{Z}$, respectively.

Let $\mathbb{Z}_{+}=\left\{0,1,2,\cdots\right\}$. We write the set of real-valued symmetric $n\times n$ matrices as ${\mathbb{R}}^{n\times n}_{\rm sym}$. In the LQR problem [15], we determine the input sequence $u{\left(k\right)}$ $\left(k\in\mathbb{Z}_{+}\right)$ that minimizes the following objective function:

where $Q\in{\mathbb{R}}^{n\times n}_{\rm sym}$ and $R\in{\mathbb{R}}^{m\times m}_{\rm sym}$.

We assume that the system (1) is controllable and matrices $Q$ and $R$ are positive definite. Then, for any $x{\left(0\right)}$, the LQR problem has a unique solution written as the following linear state feedback law:

where $P\in{\mathbb{R}}^{n\times n}_{\rm sym}$ is a unique positive definite solution of the ARE:

For $u$ defined in (3), $J{\left(u\right)}=x{\left(0\right)}^{\top}Px{\left(0\right)}$ holds.

In the ILQR problem, we find positive definite matrices $Q\in{\mathbb{R}}^{n\times n}_{\rm sym}$ and $R\in{\mathbb{R}}^{m\times m}_{\rm sym}$ such that an input $u{\left(k\right)}=-Kx{\left(k\right)}$ with the given gain $K\in{\mathbb{R}}^{m\times n}$ is a solution of the LQR problem, that is, $u$ minimizes (2).

We can solve the ILQR problem by solving an LMI. Let the system (1) be controllable. Then, $u{\left(k\right)}=-Kx{\left(k\right)}$ minimizes (2) if and only if a positive definite matrix $P\in{\mathbb{R}}^{n\times n}_{\rm sym}$ exists that satisfies (4) and the following:

By transforming (4) and (5), we obtain the following pair of equations equivalent to (4) and (5):

Hence, we can solve the ILQR problem by determining $P$, $Q\in{\mathbb{R}}^{n\times n}_{\rm sym}$, and $R\in{\mathbb{R}}^{m\times m}_{\rm sym}$ that satisfy the following LMI:

In biological system analysis or reverse engineering, the system model and control gain is often unknown, and the ARE (6) is not readily available. Hence, we consider the following problem of determining a linear equation equivalent to (6) using the system state and input observation:

We discuss without assuming that the data is sequential, that is, always satisfies $x_{i+1}{\left(0\right)}=x_{i}{\left(1\right)}$. However, in an experiment, we show that our result can also be applied to one sequential data. The ARE has multiple solutions and we call the set of these solutions, which is a linear subspace, a solution space.

The simplest solution to Problem 1 is to identify the control system, that is, matrices $A$, $B$ and $K$. We define matrices $X\in\mathbb{R}^{n\times N_{d}}$, $X^{\prime}\in\mathbb{R}^{n\times N^{\prime}_{d}}$, $U\in\mathbb{R}^{m\times N_{d}}$, $U^{\prime}\in\mathbb{R}^{m\times N^{\prime}_{d}}$, and $D\in\mathbb{R}^{n+m\times N_{d}}$ as follows:

Let the matrices $D$ and $X^{\prime}{\left(0\right)}$ have row full rank. Then, matrices $A$, $B$, and $K$ are identified using the least square method as follows:

In ILQR problems, prior information about matrices $Q$ and $R$ may be obtained. However, such information cannot be used for system identification. We propose a novel method that uses the prior information and estimates the ARE using less observation data than system identification.

The following theorem provides an estimation of the ARE:

Equation (13) is the estimated linear equation that we propose in this paper and it can be obtained directly from the observation data without system identification. Our method obtains one scalar linear equation of $P$, $Q$, and $R$ from a pair of observation data. However, at least one of the paired data must be a transition by a linear feedback input with the given gain $K$. Equation (13) consists of $N_{d}N^{\prime}_{d}$ scalar equations obtained in such a manner. Therefore, it is expected that if $N_{d}$ and $N^{\prime}_{d}$ are sufficiently large, (13) is equivalent to (6), which is also linear for the matrices $P$, $Q$, and $R$.

Equation (12) is similar to the Bellman equation [16]. Suppose that the ARE (6) holds, that is, $u{\left(k\right)}=-Kx{\left(k\right)}$ minimizes $J{\left(u\right)}$. Then, the following Bellman equation holds:

The equation $f_{i,j}=0$ is equivalent to the Bellman equation if $x_{i}{\left(0\right)}=x_{j}{\left(0\right)}$ and $u_{i}=u_{j}=-Kx_{i}{\left(0\right)}$. The novelty of Theorem 1 is that $f_{i,j}=0$ holds under the weaker condition, that is, only $u_{i}=-Kx_{i}{\left(0\right)}$. Theorem 1 is also proved from the Bellman equation (17). See the supplemental material for the details.

The following theorem demonstrates that if the data satisfies the condition required for system identification, our estimation is equivalent to the ARE:

The advantage of our method is the data economization that results from having prior information about $Q$ and $R$. Suppose that some elements of $Q$ and $R$ are known to be zero in advance, and independent ones are fewer than the scalar equations in the ARE (6). Then, there is a possibility that our method can estimate the ARE with fewer data than $m+n$ required for system identification.

For fixed $N_{d}$, our method can provide the largest number of scalar equations when $N^{\prime}_{d}=\min\left\{n,N_{d}\right\}$. We use the following lemma to explain the reason:

Lemma 1 means that if data $d_{k}$ is the linear combination of other data, the equations obtained from $d_{k}$ are also linear combinations of other equations not using $d_{k}$. Hence, without noise, such a data as $d_{k}$ is meaningless. If $N^{\prime}_{d}>n$, at least one of $d_{1},\ldots,d_{N^{\prime}_{d}}$ is the linear combination of other data because the inputs in those data are obtained from the same gain $K$. Therefore, $N^{\prime}_{d}$ larger than $n$ reduces data efficiency. From the above and $f_{i,j}=f_{j,i}$, our method can provide $\frac{n\left(n+1\right)}{2}+n\left(N_{d}-n\right)$ equations if $N_{d}\geq n$.

For an example of prior information, we consider diagonal $Q$ and $R$. In this case, the number of the independent elements of $P$, $Q$, and $R$ is $\frac{n\left(n+1\right)}{2}+n+m$. Hence, if $N_{d}$ is $N_{d\rm min}{\left(n,m\right)}=n+1+\lceil\frac{m}{n}\rceil$ or larger, we have more equations than the independent elements. Fig. 1 compares the numbers $N_{d\rm min}$ and $n+m$ of data required for our method and system identification, respectively, if $Q$ and $R$ are diagonal. From Fig. 1, our method may estimate the ARE with less data than system identification. Additionally, this tendency becomes stronger as the number $m$ of inputs becomes larger than the number $n$ of states.

In the case of noisy data, we discuss the unbiasedness of the coefficients in (13). The coefficient of $k$th row $l$th column element $q_{k,l}$ of $Q$ in (12) is expressed as follows:

where $x_{i,k}{\left(0\right)}$ is the $k$th element of $x_{i}{\left(0\right)}$. Suppose that the data is a sum of the true value and zero-mean noise distributed independently for each data. Then, we can readily confirm that if $i\neq j$, the coefficient (23) of $q_{k,l}$ in (12) is unbiased, and otherwise, it is biased. This result is also the same for the coefficients of the elements of $P$ and $R$.

In this paper, we proposed a method to estimate the ARE with respect to an unknown discrete-time system from the input and state observation data. Our method transforms the ARE into a form calculated without the system model by multiplying the observation data on both sides. We proved that our estimated equation is equivalent to the ARE if the data are sufficient for system identification. The main feature of our method is the direct estimation of the ARE without identifying the system. This feature enables us to economize the observation data using prior information about the objective function. We conducted a numerical experiment that demonstrated that that our method requires less data than system identification if the prior information is sufficient.

We give another proof of Theorem 1 using the following Bellman equation:

To this end, we prove the following theorem equivalent to Theorem 1.

Theorem 1’: Consider Problem 1. Let $i$, $j\in\left\{1,\ldots,N_{d}\right\}$ be arbitrary natural numbers. Suppose that the following ARE holds:

Then, we have $f_{i,j}=0$ if $u_{i}=-Kx_{i}{\left(0\right)}$ holds, where $f_{i,j}$ is defined as