An Efficient Off-Policy Reinforcement Learning Algorithm for the Continuous-Time LQR Problem

Consider the integer $N\in\mathbb{N}$ and the positive scalar $T~{}\in~{}\mathbb{R}_{+}$. Let $\xi:[0,NT]\rightarrow\mathbb{R}^{\sigma}$, with $[0,NT]\subset\mathbb{R}$, denote a continuous-time signal of length $NT$. Using the trajectory $\xi$, we define the following matrix

$$\mathcal{H}_{T}(\xi(t)):=\left[\begin{array}[]{cccc}\xi(t)&\xi(t+T)&\cdots&\xi(t+(N-1)T)\end{array}\right]$$

(1)

for $0\leq t\leq T$. Notice that (1) is a time-varying matrix defined on the interval $t\in[0,T]$.

Now, consider the following CT-LTI system

$$\dot{x}(t)=Ax(t)+Bu(t),$$

(2)

where $x\in\mathbb{R}^{n}$ and $u\in\mathbb{R}^{m}$ are the state and input vectors of the system, respectively. The pair $(A,B)$ is assumed to be controllable throughout the paper.

Suppose that the input signal $u:[0,NT]\rightarrow\mathbb{R}^{m}$ is applied to (2), and the resulting state trajectory $x:[0,NT]\rightarrow\mathbb{R}^{n}$ is collected. From (2) and the definition in (1), we can write

$$\mathcal{H}_{T}(\dot{x}(t))=A\mathcal{H}_{T}(x(t))+B\mathcal{H}_{T}(u(t)).$$

Since it is unusual to have the state derivative $\dot{x}$ available as a measurement, integrate the expression above to obtain

$$\mathcal{H}_{T}(x(T))-\mathcal{H}_{T}(x(0))\\ =A\int_{0}^{T}\mathcal{H}_{T}(x(\tau))d\tau+B\int_{0}^{T}\mathcal{H}_{T}(u(\tau))d\tau.$$

For convenience of notation, define the matrices

$$\tilde{X}=\mathcal{H}_{T}(x(T))-\mathcal{H}_{T}(x(0)),$$

(3)

$$X=\int_{0}^{T}\mathcal{H}_{T}(x(\tau))d\tau,\quad U=\int_{0}^{T}\mathcal{H}_{T}(u(\tau))d\tau.$$

Notice that the matrix $X$ (and similarly $U$) only requires the computation of integrals of the form ${\int_{0}^{T}x(\tau+jT)d\tau}$, $j=0,\ldots,N-1$. This is simpler than the integrals computed in the existing RL literature [8, 6, 7].

By definition, the following expression holds

$$\tilde{X}=AX+BU.$$

(4)

Define the integer constants $L,N\in\mathbb{N}$. The Hankel matrix of depth $L$ of a discrete-time sequence $\{\mu_{k}\}_{k=0}^{N-1}=\{\mu_{0},\,\mu_{1},\,\ldots,\,\mu_{N-1}\}$, $\mu_{k}\in\mathbb{R}^{m}$, is defined as

$$H_{L}(\mu):=\left[\begin{array}[]{cccc}\mu_{0}&\mu_{1}&\cdots&\mu_{N-L}\\ \mu_{1}&\mu_{2}&\cdots&\mu_{N-L+1}\\ \vdots&\vdots&\ddots&\vdots\\ \mu_{L-1}&\mu_{L}&\cdots&\mu_{N-1}\end{array}\right].$$

In [16], the following definition of a PE input for discrete-time systems is made.

It is important to highlight the fact that Definition 1 provides a condition that enables a straightforward design of a PE input and that is easy to verify for any discrete sequence.

It is shown in [21] that a piecewise constant input designed by exploiting Definition 1 is persistently exciting for the continuous-time system (2). This class of inputs is formally described in the following definition.

When a PCPE input is applied to system (2), the obtained input-state data set satisfies an important rank condition, as shown below.

For a CT-LTI system (2), the infinite-horizon LQR problem concerns determining the control input $u$ that minimizes a cost function of the form

$$J(x(0),u):=\int_{0}^{\infty}\left(x^{\top}(t)Qx(t)+u^{\top}(t)Ru(t)\right)dt,$$

(8)

where $Q\succeq 0$ and $R\succ 0$. Throughout the paper, we assume that the pair $(A,Q^{1/2})$ is observable. This, together with the assumed controllability of $(A,B)$, implies that the optimal control input is given by $u^{*}(x)=-K^{*}x$, where

$$K^{*}=R^{-1}B^{\top}P^{*}$$

and the matrix $P^{*}\succ 0$ solves the algebraic Riccati equation

$$Q+P^{*}A+A^{\top}P^{*}-P^{*}BR^{-1}B^{\top}P^{*}=0.$$

In [22], Kleinman proposed a model-based iterative algorithm to solve the LQR problem. This algorithm starts by selecting an initial stabilizing matrix $K_{0}$, i.e., a matrix such that $A-BK_{0}$ is Hurwitz stable. At every iteration $i$, the Lyapunov equation

$$P_{i}(A-BK_{i})+(A-BK_{i})^{\top}P_{i}+Q+K_{i}^{\top}RK_{i}=0$$

(9)

is solved for $P_{i}$. Then, a new feedback matrix is defined as

$$K_{i+1}=R^{-1}B^{\top}P_{i}.$$

(10)

The algorithm iterates the equations (9) and (10) until convergence. With the main drawback of being a model-based method, Kleinman’s algorithm otherwise possesses highly attractive features. Namely, at each iteration the matrix $K_{i+1}$ is stabilizing, the algorithm converges such that

$$\lim_{i\rightarrow\infty}K_{i+1}=K^{*},$$

and convergence occurs at a quadratic rate [22].

The following section presents the main developments of this paper.

Combining (9) and (10), we readily obtain the following expressions

$$P_{i}A-K_{i+1}^{\top}RK_{i}+A^{\top}P_{i}-K_{i}^{\top}RK_{i+1}+Q+K_{i}^{\top}RK_{i}=0$$

and $B^{\top}P_{i}-RK_{i+1}=0$. Therefore, the matrix equation

$$\left[\begin{array}[]{cc}A&B\\ -RK_{i}&-R\end{array}\right]^{\top}\left[\begin{array}[]{c}P_{i}\\ K_{i+1}\end{array}\right]\left[\begin{array}[]{cc}I_{n}&0\end{array}\right]\\ +\left[\begin{array}[]{c}I_{n}\\ 0\end{array}\right]\left[\begin{array}[]{c}P_{i}\\ K_{i+1}\end{array}\right]^{\top}\left[\begin{array}[]{cc}A&B\\ -RK_{i}&-R\end{array}\right]\\ +\left[\begin{array}[]{cc}Q+K_{i}^{\top}RK_{i}&0\\ 0&0\end{array}\right]=0$$

(11)

holds, where $I_{n}$ is an $n\times n$ identity matrix and $0$ represents a matrix of zeros with appropriate dimensions.

Denoting the fixed matrices as

	$$\displaystyle\Phi_{i}:=\left[\begin{array}[]{cc}A&B\\ -RK_{i}&-R\end{array}\right],\quad E:=\left[\begin{array}[]{cc}I_{n}&0\end{array}\right],$$		(15)
	$$\displaystyle\bar{Q}_{i}:=\left[\begin{array}[]{cc}Q+K_{i}^{\top}RK_{i}&0\\ 0&0\end{array}\right]$$		(18)

and the unknown matrix as

$$\Theta_{i+1}:=\left[\begin{array}[]{c}P_{i}\\ K_{i+1}\end{array}\right],$$

(19)

we can write (11) in the compact form

$$\Phi_{i}^{\top}\Theta_{i+1}E+E^{\top}\Theta_{i+1}^{\top}\Phi_{i}+\bar{Q}_{i}=0.$$

(20)

The matrix $\Theta_{i+1}\in\mathbb{R}^{(n+m)\times n}$ consists of the unknown matrices in Kleinman’s algorithm, $P_{i}$ and $K_{i+1}$. It is of our interest to design a method in which solving a matrix equation as in (20) for $\Theta_{i+1}$ corresponds to solving both (9) and (10) simultaneously. However, it can be noted that (20), as it is formulated, in general does not have a unique solution $\Theta_{i+1}$. To address this issue, first express the unknown submatrices of $\Theta_{i+1}$ as

$$\Theta_{i+1}=\left[\begin{array}[]{c}\Theta_{i+1}^{1}\\ \Theta_{i+1}^{2}\end{array}\right],$$

(21)

with $\Theta_{i+1}^{1}\in\mathbb{R}^{n\times n}$ and $\Theta_{i+1}^{2}\in\mathbb{R}^{m\times n}$. In the following lemma, we show that there exists only one matrix $\Theta_{i+1}$ that solves (20) such that the submatrix $\Theta_{i+1}^{1}$ is symmetric.

Lemma 2 implies that constraining the solution of (20) to include a symmetric submatrix $\Theta_{i+1}^{1}$ leads to the desired solution (19). The following lemma shows that we achieve this by properly modifying $\Phi_{i}$ in (18).

Using this result, we formulate Algorithm 1 below. As in any policy iteration procedure, Algorithm 1 is initialized with a stabilizing matrix $K_{0}$. Using this matrix (as well as model knowledge), (23) is solved for $\Theta_{i+1}$. Then, partitioning $\Theta_{i+1}$ as in (21), a new feedback matrix is obtained as $K_{i+1}=\Theta_{i+1}^{2}$.

Using the results obtained so far, we conclude that Algorithm 1 is equivalent to Kleinman’s algorithm in the sense that, starting from the same initial matrix $K_{0}$, they provide the same updated policies $K_{i+1}$ at every iteration. This implies that Algorithm 1 preserves all the properties of Kleinman’s algorithm. In the following, we use this result to design a data-based algorithm.

To avoid the need for model knowledge in Algorithm 1, we collect persistently excited data from the system (2) as described in Section II-C. Using this data, we define the constant matrices $X$, $U$ and $\tilde{X}$ as in (3).

Lemma 1 showed that the collected data set satisfies the rank condition (7). In the following lemma, we extend this result to the matrices $X$ and $U$.

Define $Z=[X^{\top}\quad U^{\top}]^{\top}$. Since $Z$ has full row rank by Lemma 4, we can select $n+m$ linearly independent columns from it. Let $z_{k}$ represent the $k$th column of $Z$, and let $\eta=\{k_{1},\ldots,k_{n+m}\}$ be a set of indices such that

$$Z_{\eta}:=\left[z_{k_{1}}\quad\cdots\quad z_{k_{n+m}}\right]$$

(26)

is a nonsingular matrix. Then, $\Theta_{i+1}$ is a solution of (23) if and only if it is a solution of

$$Z_{\eta}^{\top}(\Phi_{i}^{-})^{\top}\Theta_{i+1}EZ_{\eta}+Z_{\eta}^{\top}E^{\top}\Theta_{i+1}^{\top}\Phi_{i}^{+}Z_{\eta}+Z_{\eta}^{\top}\bar{Q}_{i}Z_{\eta}=0.$$

(27)

From the definitions in (18) and (24), and using the expression (4), we have the following

$$\Phi_{i}^{+}Z_{\eta}=\left[\begin{array}[]{c}AX_{\eta}+X_{\eta}+BU_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right]=\left[\begin{array}[]{c}\tilde{X}_{\eta}+X_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right],$$

$$\Phi_{i}^{-}Z_{\eta}=\left[\begin{array}[]{c}AX_{\eta}-X_{\eta}+BU_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right]=\left[\begin{array}[]{c}\tilde{X}_{\eta}-X_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right],$$

$Z_{\eta}^{\top}\bar{Q}_{i}Z_{\eta}=X_{\eta}^{\top}(Q+K_{i}^{\top}RK_{i})X_{\eta}$ and $EZ_{\eta}=X_{\eta}$, where the subindex $\eta$ represents a matrix constructed using the columns specified by the set $\eta$ from the corresponding original matrix. Substituting in (27), we obtain

$$(Y_{i}^{-})^{\top}\Theta_{i+1}X_{\eta}+X_{\eta}^{\top}\Theta_{i+1}^{\top}Y_{i}^{+}+X_{\eta}^{\top}(Q+K_{i}^{\top}RK_{i})X_{\eta}=0.$$

(28)

where

	$\displaystyle Y_{i}^{-}$	$\displaystyle:=\left[\begin{array}[]{c}\tilde{X}_{\eta}-X_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right],$		(29)
	$\displaystyle Y_{i}^{+}$	$\displaystyle:=\left[\begin{array}[]{c}\tilde{X}_{\eta}+X_{\eta}\\ -RK_{i}X_{\eta}-RU_{\eta}\end{array}\right].$

Now, (28) is a data-based equation that does not require any knowledge about the system model. Algorithm 2 uses this expression to solve the LQR problem. For convenience, for Algorithm 2 we define

$$Q_{i}:=X_{\eta}^{\top}(Q+K_{i}^{\top}RK_{i})X_{\eta}.$$

(30)

The following theorem states the main properties of this algorithm.

Algorithm 2 is a purely data-based, off-policy method to solve the continuous-time LQR problem. Using Definition 2, we are able to guarantee the existence of a solution $\Theta_{i+1}$ of (31) at every iteration for data trajectories of fixed length. This contrasts with the methods in the literature that must keep collecting data until a matrix gets full rank, such as, e.g., [7, 8]. Moreover, we avoid the use of the Kronecker product and its resulting large matrices in Algorithm 2. As stated in Remark 4, methods to efficiently solve a Sylvester-transpose equation as in (31) are well known.

In this subsection, we analyze the theoretical computational complexity of Algorithm 2. Moreover, we compare this complexity with that of the algorithm proposed in [8]. This is because [8] is also an off-policy data-based method that shares many of the characteristics of Algorithm 2.

The most expensive steps in Algorithm 2 are obtaining the solution of (31) and selecting $n+m$ linearly independent vectors from $[X^{\top}\quad U^{\top}]^{\top}$. Methods to solve the Sylvester-transpose equation (31) with a complexity of $\mathcal{O}((n+m)^{3})$ are known [19]. The selection of linearly independent vectors can be performed using a simple procedure like Gaussian elimination to transform the matrix of interest into row echelon form. This method has a complexity of $\mathcal{O}((n+m)^{2}N)$ operations [24]. This step, however, only needs to be performed once in Algorithm 2 (in Step 4). Thus, we conclude that Algorithm 2 requires once $\mathcal{O}((n+m)^{2}N)$ and then in each iteration $\mathcal{O}((n+m)^{3})$ floating point operations.

The algorithm in [8] was also shown to be equivalent to Kleinman’s algorithm at every iteration. However, their method uses a Kronecker product formulation that yields matrices of large dimensions. Let $N_{\otimes}$ be the amount of data samples used in [8]. Then, the most expensive step at each iteration of their algorithm is the product of a matrix with dimensions $(\frac{1}{2}n(n+1)+mn)\times N_{\otimes}$ times its transpose. This product, and hence each iteration of the algorithm, requires $\mathcal{O}((\frac{1}{2}n(n+1)+mn)^{2}N_{\otimes})$ floating point operations [25]. Clearly, as the dimension of the system increases, the difference in performance of both algorithms becomes more significant. Moreover, we notice from [8] that the amount of collected data must satisfy $N_{\otimes}\geq\frac{1}{2}n(n+1)+mn$ for the algorithm to yield a unique solution at every iteration. Compare this with the bound $N\geq(n+1)m+n$ in Algorithm 2. In Section V, we test this theoretical comparison using numerical examples.

In [26, Remark 2], a procedure to design a stabilizing controller for continuous-time systems using only measured data was described. This method is based on the solution of a linear matrix inequality (LMI). The authors in [7] proposed to use a similar LMI-based procedure to determine the initial stabilizing gain for a Q-learning algorithm. Since one of the goals in this paper is computational efficiency, we would like to avoid the computationally expensive step of solving an LMI. In this subsection, we present an alternative method to determine the initial stabilizing matrix $K_{0}$ for Algorithm 2. The following development follows closely a procedure proposed in [27, Section IV] for discrete-time systems.

Let $F$ be the Moore-Penrose pseudoinverse of the matrix $X$ in (3). Since $X$ has full row rank (see Lemma 4), $F$ is a right inverse of $X$. Furthermore, let $G$ be a basis for the null space of $X$, such that $X(F-G\bar{K})=I$ for any matrix $\bar{K}$ of appropriate dimensions. Using the matrices $F$, $G$ and $U$ from (3), we propose to compute the initial stabilizing gain $K_{0}$ for Algorithm 2 as

$$K_{0}=-U(F-G\bar{K})$$

(32)

where $\bar{K}$ is a matrix to be determined.

From (4) and (32), notice that

	$\displaystyle\tilde{X}(F-G\bar{K})$	$\displaystyle=\left[A\quad B\right]\left[\begin{array}[]{c}X\\ U\end{array}\right](F-G\bar{K})$
		$\displaystyle=\left[A\quad B\right]\left[\begin{array}[]{c}I\\ -K_{0}\end{array}\right]$

Therefore, by designing the poles of the matrix $\tilde{X}(F-G\bar{K})$, we also set the poles of $A-BK_{0}$. Since $(A,B)$ is controllable and hence the poles of $A-BK_{0}$ can be assigned arbitrarily, also the poles of $\tilde{X}(F-G\bar{K})$ can be placed arbitrarily by a suitable choice of $\bar{K}$. Moreover, since $\tilde{X}$, $F$ and $G$ are matrices obtained from data, we can operate with them without any need of model knowledge. This procedure is summarized in the following theorem. The proof of this theorem is straightforward considering the procedure described in this subsection and is hence omitted.