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

In this paper, we propose a method for estimating the algebraic Riccati equation (ARE) with respect to an unknown discrete-time system from the system state and input observation. The inverse optimal control (IOC) problem asks, ``What objective function is optimized by a given control system?'' The inverse linear quadratic regulator (ILQR) problem is an IOC problem that assumes a linear system and quadratic objective function. The ILQR problem can be solved by solving a linear matrix inequality that contains the ARE. However, the system model is required to obtain the ARE, and it is often unknown in fields in which the IOC problem occurs, for example, biological system analysis. Our method directly estimates the ARE from the observation data without identifying the system. This feature enables us to economize the observation data using prior information about the objective function. We provide a data condition that is sufficient for our method to estimate the ARE. 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.


I. INTRODUCTION
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.

II. PROBLEM FORMULATION
We consider the following discrete-time linear system: where A ∈ R n×n and B ∈ R n×m are constant matrices, and x(k) ∈ R n and u(k) ∈ R m are the state and the input at time k ∈ Z, respectively. Let We write the set of real-valued symmetric n × n matrices as R n×n sym .In the LQR problem [15], we determine the input sequence u(k) (k ∈ Z + ) that minimizes the following objective function: where Q ∈ R n×n sym and R ∈ R m×m sym .We assume that the system (1) is controllable and matrices Q and R are positive definite.Then, for any x(0), the LQR problem has a unique solution written as the following linear state feedback law: where P ∈ R n×n sym is a unique positive definite solution of the ARE: For u defined in (3), J(u) = x(0) ⊤ P x(0) holds.
In the ILQR problem, we find positive definite matrices We can solve the ILQR problem by solving an LMI.Let the system (1) be controllable.Then, u(k) = −Kx(k) minimizes (2) if and only if a positive definite matrix P ∈ R n×n 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 ∈ R n×n sym , and R ∈ R m×m 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: Problem 1 (ARE estimation problem): Consider controllable system (1) and controller u(k) = −Kx(k) with unknown A, B, and K. Suppose N d observation data (x i (0), u i , x i (1)) (i ∈ {1, . . ., N d }) of the system state and input are given, where Let N ′ d ≤ N d .N ′ d inputs in the data are obtained from the unknown controller as follows: Determine a linear equation of P , Q ∈ R n×n sym , and R ∈ R m×m sym with the same solution space as (6).We discuss without assuming that the data is sequential, that is, always satisfies x i+1 (0) = x i (1).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.

III. DATA-DRIVEN ESTIMATION OF ARE
The simplest solution to Problem 1 is to identify the control system, that is, matrices A, B and K.We define matrices , and D ∈ R n+m×N d as follows: Let the matrices D and X ′ (0) 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: Theorem 1: Consider Problem 1.Let F be an N ′ d × N d matrix whose ith row jth column element f i,j is defined as ) Then, the following condition is necessary for the ARE (6) to hold.
The ARE ( 6) is equivalent to the combination of G 1 = 0 and G 2 = 0. Let i ∈ {1, . . ., N ′ d } and j ∈ {1, . . ., N d } be arbitrary natural numbers.From the ARE, we obtain Using the assumption ( 8) and ( 9), we can transform the lefthand side of (15) into the following: Therefore, we obtain f i,j = 0 from ( 15) and ( 16), which proves Theorem 1. 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 ′ d scalar equations obtained in such a manner.Therefore, it is expected that if N d and N ′ 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(k) = −Kx(k) minimizes J(u).Then, the following Bellman equation holds: The equation f i,j = 0 is equivalent to the Bellman equation if x i (0) = x j (0) and u i = u j = −Kx i (0).The novelty of Theorem 1 is that f i,j = 0 holds under the weaker condition, that is, only u i = −Kx i (0).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: Theorem 2: Consider Problem 1. Suppose that matrices D and X ′ (0) defined in (10) have row full rank.Then, the estimation ( 13) is equivalent to the ARE (6).
Proof: From Theorem 1, the ARE (6) implies the estimation (13), that is, F = 0.As shown in the proof of Theorem 1, matrix F is expressed as follows using the matrices defined in (10) and ( 14): Because the matrix D has row full rank, F = 0 implies Because X ′ (0) has row full rank, we obtain G 1 = 0 and G 2 = 0 from (19), which proves Theorem 2. 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 ′ d = min {n, N d }.We use the following lemma to explain the reason: Lemma 1: Suppose that there exist k ∈ {1, . . ., N ′ d } and where for any i ∈ {1, . . ., N d }.Then, the following holds for any j ∈ {1, . . ., N d }: Proof: We can readily prove Lemma 1 using following derived from (18) 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 ′ d > n, at least one of d 1 , . . ., d N ′ d is the linear combination of other data because the inputs in those data are obtained from the same gain K. Therefore, N ′ d larger than n reduces data efficiency.From the above and f i,j = f j,i , our method can provide 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 n(n+1) n ⌉ or larger, we have more equations than the independent elements.Fig. 1 compares the numbers N dmin 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 kth row lth column element q k,l of Q in ( 12) is expressed as follows: where x i,k (0) is the kth element of x i (0).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 ̸ = 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.

IV. NUMERICAL EXPERIMENT A. Distance Between Solution Spaces
In our experiments, we evaluated the ARE estimation using a distance between the solution spaces of the ARE and estimation.Let s ∈ R Nv be a vector generated by independent elements in the matrices P , Q, and R in any fixed order, where N v is defined according to the prior information.Because the ARE ( 6) is linear for s, we can transform (6) into the following equivalent form: where Θ ∈ R NARE×Nv is the coefficient matrix and + mn.Similarly, the estimation ( 13) is transformed into the following equivalent form: where Θ ∈ R Nest×Nv is the coefficient matrix and We define the solution spaces S and Ŝ of the ARE (24) and estimation (25) as follows: (26) We define the distance between the solution spaces using the approach provided [17].Assume that Θ and Θ have the same rank.Let Π, Π ∈ R Nv×Nv be the orthogonal projections to S and Ŝ, respectively.The distance between S and Ŝ is defined as follows: where ∥•∥ 2 is the matrix norm induced by the 2-norm of vectors.Distance d S, Ŝ is the maximum Euclidian distance between ŝ ∈ Ŝ and Πŝ ∈ S when ∥ŝ∥ 2 = 1, as explained in Fig. 2 and by the following equation: Hence, the distance satisfies 0 ≤ d S, Ŝ ≤ 1.
The orthogonal projections Π and Π are obtained from the orthogonal bases of S and Ŝ.To obtain these orthogonal bases, we use singular value decompositions of Θ and Θ.
In our experiments, we performed the computation using MATLAB 2023a.

B. Experiment 1
In this experiment, we confirmed Theorem 2 by solving Problem 1.We set the controllable pair (A, B) as follows: Using this setting, each of the ARE (6) and our estimation (13) contained 12 scalar equations for N v = 15 variables, that is, independent elements in the symmetric matrices P , Q, and R. Fig. 2 shows the singular values of Θ and Θ, in descending order, that we computed using the default significant digits.As shown in Fig. 2, there is no zero singular value of Θ, Θ ∈ R 12×15 .Hence, the ranks of Θ and Θ are 12, and the solution spaces S, Ŝ ⊂ R 15 are three-dimensional subspaces.To investigate the influence of the computational error, we performed the same experiment with various numbers of significant digits.Fig. 3 shows the relationship between the distance d S, Ŝ and the number of significant digits.As shown in Fig. 3, the logarithm of the distance decreased in proportion to the number of significant digits.Therefore, we concluded that non-zero d S, Ŝ was caused by the computational error and that our method correctly estimated the ARE.

C. Experiment 2
In this experiment, we confirmed that, if Q and R are diagonal, Theorem 1 can estimate the ARE with fewer observation data than system identification.We randomly set the matrices A ∈ R 100×100 and B ∈ R 100×50 by generating each element from the uniform distribution in [−1, 1].After the generation, we confirmed that (A, B) is controllable.The given K ∈ R 50×100 is the solution of the LQR problem for the diagonal matrices Q and R whose elements were generated from the uniform distribution in [0.01, 1].We used N d = N dmin = 102 observation data that satisfy (9) for N ′ d = n = 100.We generated the elements of x 1 (0), . . ., x 102 (0) ∈ R 100 and u 101 , u 102 ∈ R 50 from the uniform distribution in [−1 , 1].The number of used data, 102, is less than n + m = 150 required for system identification (11).
Using this setting, the ARE (6) and our estimation (13) contained + n + m = 5200 independent elements in the symmetric matrix P and diagonal matrices Q and R.
Because the arbitrary-precision arithmetic was too slow to perform for the scale of this experiment, we performed the computation using the default significant digits.We indexed the singular values of Θ ∈ R 10050×5200 and Θ ∈ R 5250×5200 in descending order; Fig. 4 shows the 5180-5200th singular values.From Fig. 4, we assumed that, without computational error, Θ and Θ each had a zero singular value.Then, the solution spaces S, Ŝ ⊂ R 5200 were one-dimensional subspaces, and the distance d S, Ŝ was 4.3 × 10 −10 .We compared our method with system identification.For system identification, we generated additional data x 103 , . . ., x 150 and u 103 , . . ., u 150 in the same manner as x 101 and u 101 .Using the same approach as that for S, we defined the solution space ŜSI ⊂ R Nv of the ARE (6) using A, B, and K identified by (11).Then, the distance d S, ŜSI was 2.8 × 10 −10 .Therefore, our method estimated the ARE with the same order of error as system identification using fewer observation data than system identification.

D. Experiment 3
In this experiment, we compared our method and system identification under a practical conditions: a single noisy and sequential data and sparse but not diagonal Q and R. We set the matrices A ∈ R 40×40 and B ∈ R 40×20 in the same manner as Experiment 2. The control gain K ∈ R 20×40 is the solution of the LQR problem for the sparse matrices Q and R. We performed three operations to generate Q and R. First, we generated matrices M Q ∈ R 40×40 and M R ∈ R 20×20 in the same manner as A and set Second, we randomly set 800 non-diagonal elements of Q to 0 to make it sparse while maintaining symmetry.Similarly, we set 200 elements of R to 0. Finally, add the product of a constant and the identity matrix to Q and R so that the maximum eigenvalues are 10 times as large as the minimum ones, respectively.
To generate data, we used the following system: x where x * (k) ∈ R n and u * (k) ∈ R m are the state and input at time k ∈ Z + , respectively.If ∥x * (k)∥ 2 ≤ 1, the vector v(k) ∈ R n with norm 0.2 is the product of a scalar constant and a random vector whose elements follow a uniform distribution in [−1, 1].Otherwise, v(k) = 0. We ran the system (33) through N d = 200 time steps with a random initial state x * (0) satisfying ∥x * (0)∥ 2 ≤ 1 and obtained data as follows: where ε state ∈ R n and ε input ∈ R m are the observation noise whose elements follow a normal distribution with a mean 0 and variance σ 2 .We explain the value of σ 2 later.Throughout the simulation, v(k) = 0 holds 121 times.Hence, N ′ d = 121.We sorted the data (x i (0), u i , x i (1)) (i ∈ {1, . . ., N d }) so that (9) holds if noise does not exist.
We indexed the singular values of Θ ∈ R 1620×1350 and Θ ∈ R 16940×1350 in descending order.Fig. 5 shows the 1330-1350th singular values when σ 2 = 10 −8 .Although there is no zero singular value due to noise, the last singular value is noticeably small as shown in Fig. 5. Hence, we can conclude that the solution spaces S, Ŝ ⊂ R 1350 were onedimensional subspaces.
We conducted experiments with different noise variance σ 2 from 10 −6 to 10 −16 .Also, the noise seed differs for each experiment, but the other seeds are the same.Fig. 6 shows the relationship between the noise variance σ 2 and distances d S, Ŝ and d S, ŜSI .As shown in Fig. 6, our method outperformed system identification in almost all experiments.In the experiments of variances from 10 −8 to 10 −16 , the ratio of d S, Ŝ to d S, ŜSI is approximately constant, and the average ratio is 0.17.Because the maximum value of the distances is 1, we conclude that the estimations using both methods failed with a larger noise variance than 10 −8 .

V. CONCLUSIONS
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.

Fig. 1 .
Fig. 1.Ratio N dmin n+m when 5 ≤ n ≤ 200 and 5 ≤ m ≤ 200.The contour lines are not smooth, but this is not a mistake.

2 
2 −0.4 −0.6 0.4 −0.7 −0.3 −1.0 −0.8 −0.The given gain K ∈ R 2×3 is the solution of the LQR problem for the following Q and R: following N d = n + m observation data with N ′ d = n that satisfy the condition of Theorem 2:

Fig. 3 .
Fig. 3. Relationship between distance d S, Ŝ and number of significant digits in the computation in Experiment 1.