Suboptimality analysis of receding horizon quadratic control with unknown linear systems and its applications in learning-based control

We consider the discrete-time linear system

$\displaystyle x_{t+1}=A_{\star}x_{t}+B_{\star}u_{t}+w_{t},$

(2)

where $w_{t}$ is a white noise signal with $\mathbb{E}(w_{t}w_{t}^{\top})=\sigma_{w}^{2}I$, $\sigma_{w}>0$, $A_{\star}\in\mathbb{R}^{n\times n}$, and $B_{\star}\in\mathbb{R}^{n\times m}$. The initial state $x_{0}$ is a zero-mean random vector with $\mathbb{E}(x_{0}x_{0}^{\top})=\sigma_{x}^{2}I$ and $\sigma_{x}\geqslant 0$, and it is uncorrelated with the noise.

In this work, the true system $(A_{\star},B_{\star})$ is unknown, and we have access to an approximate model $(\hat{A},\hat{B})$ that satisfies $\|\hat{A}-A_{\star}\|\leqslant\varepsilon_{\mathrm{m}}$ and $\|\hat{B}-B_{\star}\|\leqslant\varepsilon_{\mathrm{m}}$ for some $\varepsilon_{\mathrm{m}}\geqslant 0$. This approximate model and its error bound can be obtained, e.g., from linear system identification [22, 23].

In the above setting, we consider the nominal receding-horizon LQ controller [10]. At time step $t$, $x_{t}$ is measured, and the RHC controller solves the following optimization problem:

	$\displaystyle\min_{\{u_{k\mid t}\}_{k=0}^{N-1}}\mathbb{E}_{\{w_{k+t}\}_{k=0}^{N-1}}\Big{[}\sum_{k=0}^{N-1}l(x_{k\mid t},u_{k\mid t})+x_{N\mid t}^{\top}Px_{N\mid t}\Big{]},$		(3a)
	$\displaystyle\text{ s.t. }x_{k+1\mid t}=\hat{A}x_{k\mid t}+\hat{B}u_{k\mid t}+w_{k+t},\text{ }x_{0|t}=x_{t},$		(3b)

where $l(x_{k\mid t},u_{k\mid t})=x_{k\mid t}^{\top}Qx_{k\mid t}+u_{k\mid t}^{\top}Ru_{k\mid t}$, $R\in\mathbb{S}^{m}_{++}$, and $Q$, $P\in\mathbb{S}^{n}_{+}$. Then $u_{t}=u_{0\mid t}$ is the control input at time step $t$.

We will first consider controllable systems. The extension to stabilizable systems is in Section 5. To this end, we define $\mathcal{C}_{i}(A_{\star},B_{\star})\triangleq\begin{bmatrix}B_{\star}&A_{\star}B_{\star}&\dots&A_{\star}^{i-1}B_{\star}\end{bmatrix}$ for $i\in\{1,2,\dots\}$ and introduce the following assumptions:

For some intermediate technical results, we will relax the controllability requirement whenever possible for generality.

This RHC controller can be further characterized by the Riccati equation as follows. Given any $(A,B)\in\mathbb{R}^{n\times n}\times\mathbb{R}^{n\times m}$, we define

$$S_{B}\triangleq BR^{-1}B^{\top},$$

(4)

and the Riccati mapping $\mathcal{R}_{A,B}:\mathbb{S}_{+}^{n}\to\mathbb{S}_{+}^{n}$ as

	$\displaystyle\mathcal{R}_{A,B}(P)$	$\displaystyle\triangleq A^{\top}P(I+S_{B}P)^{-1}A+Q$		(5)
		$\displaystyle=A^{\top}PA-A^{\top}PB(R+B^{\top}PB)^{-1}B^{\top}PA+Q,$

where the last equality holds by the Woodbury matrix identity. To apply the Riccati mapping recursively, we define the Riccati iteration $\mathcal{R}^{(i)}_{A,B}$ via recursion: For any integer $i\geqslant 0$,

$$\mathcal{R}^{(i+1)}_{A,B}\triangleq\mathcal{R}_{A,B}\circ\mathcal{R}^{(i)}_{A,B},\text{ and }\mathcal{R}^{(0)}_{A,B}=\mathrm{id}.$$

Then the RHC controller in (3) is equivalent to a linear static controller $u_{t}=-K_{\mathrm{RHC}}x_{t}$ [10], where

$$K_{\mathrm{RHC}}\!=\!\Big{[}R+\hat{B}^{\top}\!\big{[}\mathcal{R}^{(N-1)}_{\hat{A},\hat{B}}(P)\big{]}\hat{B}\Big{]}^{-1}\!\hat{B}^{\top}\!\big{[}\mathcal{R}^{(N-1)}_{\hat{A},\hat{B}}(P)\big{]}\hat{A}.$$

(6)

We consider the expected average infinite-horizon performance of $K_{\mathrm{RHC}}$, applied to the true system (2). Given any controller gain $K$ with $\rho(A_{\star}-B_{\star}K)<1$, we denote its expected average infinite-horizon performance as

	$\displaystyle J_{K}\triangleq\lim_{T\to\infty}$	$\displaystyle\frac{1}{T}\mathbb{E}_{x_{0},\{w_{t}\}}\Big{[}\sum_{t=0}^{T-1}\big{(}x_{t}^{\top}Qx_{t}+u_{t}^{\top}Ru_{t}\big{)}\Big{]},$
		$\displaystyle\text{s.t. }\eqref{eq:TrueModel},\text{ }u_{t}=-Kx_{t},\text{ }t\in\mathbb{Z}^{+}.$

Under a stabilizing $K$, the closed-loop system satisfies

$$\lim_{t\to\infty}\mathbb{E}(x_{t}x_{t}^{\top})=\Sigma_{K},$$

(7)

for some $\Sigma_{K}\in\mathbb{S}^{n}_{++}$ that satisfies the Lyapunov equation:

$$\Sigma_{K}=(A_{\star}-B_{\star}K)\Sigma_{K}(A_{\star}-B_{\star}K)^{\top}+\sigma_{w}^{2}I.$$

Under Assumption 2, the optimal controller gain $K_{\star}$, which achieves the minimum $J_{K}$, is the LQR controller gain:

$$K_{\star}=(R+B_{\star}^{\top}P_{\star}B_{\star})^{-1}B_{\star}^{\top}P_{\star}A_{\star},$$

(8)

where $P_{\star}$ is the unique positive-definite solution of the discrete-time Riccati equation (DARE) [24]:

$$P_{\star}=\mathcal{R}_{A_{\star},B_{\star}}(P_{\star}).$$

(9)

We further define the closed-loop matrix $L_{\star}$ under $K_{\star}$ as

$$L_{\star}\triangleq A_{\star}-B_{\star}K_{\star}=(I+S_{B_{\star}}P_{\star})^{-1}A_{\star}.$$

(10)

It is well-known that $L_{\star}$ is Schur stable, and its decay rate can be characterized by the standard Lyapunov analysis:

The proof of this result is presented in Appendix for completeness. If Assumption 1 holds additionally, then $(A_{\star},B_{\star})$ is controllable, which implies that $(L_{\star},B_{\star})$ is controllable, and thus there exists a real number $\nu_{\star}$ such that

$$\underline{\sigma}(\mathcal{C}_{n_{\mathrm{cr}}}(L_{\star},B_{\star}))\geqslant\nu_{\star}>0.$$

(12)

In this work, we analyze the performance gap $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}$ between the RHC controller and the ideal LQR controller.

In Problem 1, a non-zero $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}$ reflects how the errors $\varepsilon_{\mathrm{m}}$, $\varepsilon_{\mathrm{p}}$, and the finite prediction horizon $N$ of $K_{\mathrm{RHC}}$ lead to a performance deviation from $K_{\star}$, which is an ideal controller with an infinite prediction horizon and the correct model. Choosing the infinite-horizon optimal controller as a baseline for performance comparison is standard in studies on learning-based control [19, 20] and MPC[4].

In the special case of a known system and $P=P_{\star}$, we have $K_{\mathrm{RHC}}=K_{\star}$ and $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}=0$ for any $N\geqslant 1$. However, in the case of $P\not=P_{\star}$, which happens when the system is unknown or a numerical error is present for the computed $P_{\star}$ even if the system is known, we can increase the prediction horizon $N$ to achieve a smaller $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}$. In the more complex situation where there is a modeling error, a larger $N$ leads to the propagation of the modeling error and thus may enlarge the performance gap. Our target is to characterize the above complex tradeoff among $N$, $\varepsilon_{\mathrm{m}}$, and $\varepsilon_{\mathrm{p}}$.

In this work, the obtained bounds hold regardless of whether $\varepsilon_{\mathrm{m}}$ and $\varepsilon_{\mathrm{p}}$ are coupled. The presence of coupling, e.g., $\varepsilon_{\mathrm{p}}=h(\varepsilon_{\mathrm{m}})$ for some function $h$, can be easily incorporated by plugging function $h$ into the bound $g$. In addition, the error bounds obtained will also depend on the system matrices $A_{\star}$ and $B_{\star}$. To simplify the algebraic expressions of the bounds, we upper bound the system matrices as

$$\max\{\|A_{\star}\|,\|B_{\star}\|,\|P_{\star}\|\}\leqslant\Upsilon_{\star},$$

where $\Upsilon_{\star}\geqslant\max\{1,\varepsilon_{\mathrm{m}},\varepsilon_{\mathrm{p}}\}$ is a real constant. Note that there always exists a sufficiently large $\Upsilon_{\star}$ such that the above holds.

We first characterize the performance gap between the LQR controller and a general linear controller that has a similar structure to the RHC controller (6). Given any $(A,B)\in\mathbb{R}^{n\times n}\times\mathbb{R}^{n\times m}$, we define function $\mathcal{K}_{A,B}:\mathbb{S}_{+}^{n}\to\mathbb{R}^{m\times n}$ as

$$\mathcal{K}_{A,B}(F)\triangleq(R+B^{\top}FB)^{-1}B^{\top}FA.$$

(13)

Then given the controller gain $\mathcal{K}_{A_{\star},B_{\star}}(F)$ for any $F\in\mathbb{S}_{+}^{n}$, the following result characterizes its performance gap based on the difference $\|F-P_{\star}\|$.

Lemma 2 directly follows from Lemma 4, which is presented later, as a special case of a known model. These results extend [19, Thm. 1] and [20, Prop. 7] by incorporating the terminal value function error. Lemma 2 shows that if $F$ is sufficiently close to $P_{\star}$, the controller gain $\mathcal{K}_{A_{\star},B_{\star}}(F)$ can stabilize the system with a near-optimal performance.

In Case 1, since the RHC controller gain (6) is a special case of (13) with $F=\mathcal{R}^{(N-1)}_{A_{\star},B_{\star}}(P)$, we will exploit Lemma 2 to characterize the performance gap of $K_{\mathrm{RHC}}$. Particularly, if we can quantify the difference

$$\|\mathcal{R}^{(N-1)}_{A_{\star},B_{\star}}(P)-P_{\star}\|,$$

Lemma 2 can be used to bound $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}$. Note that this difference captures the convergence of $\mathcal{R}^{(N-1)}_{A_{\star},B_{\star}}(P)$ to $P_{\star}$ as $N$ increases.

It is a classical result that $\mathcal{R}^{(i)}_{A_{\star},B_{\star}}(P)$ converges to $P_{\star}$ exponentially as $i$ increases [25]. Different convergence analyses have been investigated in the literature, leading to different convergence rates. The rates in [26] and [27, Thm. 14.4.1] are complex functions of the system parameters, characterized by the Thompson metric in [26] and a norm weighted by a special matrix in [27, Thm. 14.4.1]. These convergence rates are hard to interpret because of the specific metric and norm. On the other hand, the rates in [25, Ch. 4.4] and [28] utilize the standard matrix 2-norm and thus are easy to interpret; however, [25, Ch. 4.4] does not exploit controllability, leading to a slower convergence rate than [28]. To the best of our knowledge, [28] provides the most recent state-of-the-art analysis of convergence rates for controllable systems. Therefore, we develop our results based on [28]. While the result in [28] considers only a single convergence rate, we generalize it to incorporate two distinct convergence rates depending on the range of $i$. We also compute the rates as a function of $\beta_{\star}$ and thus of $P_{\star}$. The explicit dependency of the convergence rates on $P_{\star}$ will facilitate the subsequent perturbation analysis in Section 4, where the effect of the modeling error on $P_{\star}$ will be exploited.

To this end, we first define several new notations and functions. Recall $S_{B}$ defined in (4), and we define the function $\mathcal{L}_{A,B}:\mathbb{S}_{+}^{n}\to\mathbb{R}^{n\times n}$ as

$$\mathcal{L}_{A,B}(P)\triangleq(I+S_{B}P)^{-1}A=A-B\mathcal{K}_{A,B}(P),$$

(15)

which denotes the closed-loop matrix under the controller gain $\mathcal{K}_{A,B}(P)$, and thus $\mathcal{L}_{A_{\star},B_{\star}}(P_{\star})=L_{\star}$. In addition, for any integers $0\leqslant j\leqslant i$, define

$$\Phi^{(j:i)}_{A,B}(P)\triangleq\begin{cases}\mathcal{L}_{A,B}\big{(}P^{(j)}\big{)}\dots\mathcal{L}_{A,B}\big{(}P^{(i-1)}\big{)},&i>j\geqslant 0,\\ I,&i=j,\end{cases}$$

(16)

where $P^{(i)}$ is a shorthand notation for $\mathcal{R}^{(i)}_{A,B}(P)$. The definition shows $\Phi^{(0:i)}_{A_{\star},B_{\star}}(P_{\star})=L_{\star}^{i}$. Note that $\Phi^{(0:i)}_{A,B}(P)$ is a state-transition matrix of a time-varying linear system, whose asymptotic stability is reflected by the convergence of $\Phi^{(0:i)}_{A,B}(P)$ to $0$ as $i\to\infty$ [29].

The state-transition matrix is helpful for characterizing $\|\mathcal{R}^{(N-1)}_{A_{\star},B_{\star}}(P)-P_{\star}\|$ for any $P\in\mathbb{S}_{+}^{n}$ [28]:

$$\mathcal{R}^{(i)}_{A_{\star},B_{\star}}(P)-P_{\star}=\big{[}\Phi^{(0:i)}_{A_{\star},B_{\star}}(P)\big{]}^{\top}(P-P_{\star})L_{\star}^{i},$$

(17)

which can be verified by applying Lemma 7(a) recursively. Recall that $L_{\star}$ is Schur stable, satisfying $\lim_{i\to\infty}L^{i}_{\star}=0$. Therefore, (17) shows that as $i\to\infty$, the convergence of $\mathcal{R}^{(i)}_{A_{\star},B_{\star}}(P)$ to $P_{\star}$ is governed by the convergence of $L_{\star}^{i}$. Furthermore, exploiting the exponential decay of $\Phi^{(0:i)}_{A_{\star},B_{\star}}(P)$ can further contribute to the convergence of $\mathcal{R}^{(i)}_{A_{\star},B_{\star}}(P)$.

Overall, we can characterize the rate of convergence of $\mathcal{R}^{(i)}_{A_{\star},B_{\star}}(P)$ to $P_{\star}$ as follows:

Since $\overline{\beta}\leqslant\beta_{\star}$, we have $1-\beta_{\star}\leqslant\sqrt{(1-\beta_{\star})(1-\overline{\beta})}$. Therefore, as $i$ increases, Lemma 3 shows a slower decay rate²²2Given $c_{1}a^{i}$ and $c_{2}b^{i}$ with $0\leqslant a\leqslant b<1$ and $c_{1}$, $c_{2}\geqslant 0$, $i\in\mathbb{Z}^{+}$, we say $a$ is a faster decay rate than $b$ to indicate $a^{i}$ decays faster than $b^{i}$ as $i$ increases. $\sqrt{(1-\beta_{\star})(1-\overline{\beta})}$ for $i<n_{\mathrm{cr}}$ and a faster rate $1-\beta_{\star}$ for $i\geqslant n_{\mathrm{cr}}$. Recall $\Phi^{(0:i)}_{A_{\star},B_{\star}}(P)$ in (16) and (17). Intuitively, the case $i<n_{\mathrm{cr}}$ is when the state-transition matrix $\Phi^{(0:i)}_{A_{\star},B_{\star}}(P)$ is affected by the initial matrix $P$ and the consequent transient behaviors, which is evident from the dependency of $\overline{\beta}$ on $\varepsilon_{\mathrm{p}}$. Then, if $i\geqslant n_{\mathrm{cr}}$, the time-varying closed-loop systems in (16) get closer to $L_{\star}$, leading to a faster convergence rate.

The technical challenge in developing Lemma 3 is the derivation of the two distinct decay rates. Given controllability, the faster rate for $i\geqslant n_{\mathrm{cr}}$ is established from the full rank property of $\mathcal{C}_{i}(A_{\star},B_{\star})$ by exploiting Lemma 1 and [28, Lem. 4.1]. Since $\mathcal{C}_{i}(A_{\star},B_{\star})$ is not of full rank for $i<n_{\mathrm{cr}}$, the slower rate is established from stabilizability. The analysis is achieved by utilizing the stability analysis of linear time-varying systems [29], introduced in Lemma 9, and further deriving the decay rate as an explicit function of $\beta_{\star}$. Note that if $\varepsilon_{\mathrm{p}}=0$, we obtain the rate $(1-\beta_{\star})$ in both cases of (18), recovering the square of the rate $\sqrt{1-\beta_{\star}}$ in (11).

The following main result is a direct consequence of (6), Lemmas 2 and 3:

In Theorem 1, the upper bound on $\sqrt{v(N,\varepsilon_{\mathrm{p}})/c_{\star}}$ means that $N$ should be sufficiently large or $\varepsilon_{\mathrm{p}}$ should be sufficiently small, such that the RHC controller can stabilize the system with a suboptimal performance. Moreover, the performance gap satisfies $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}=$

$$\begin{cases}O\left(\left[(1-\beta_{\star})(1-\overline{\beta})\right]^{N-1}\varepsilon_{\mathrm{p}}^{2}\right),&\text{ for }N\leqslant n_{\mathrm{cr}},\\ O\Big{(}\left(1-\beta_{\star}\right)^{2(N-1)}\varepsilon_{\mathrm{p}}^{2}\Big{)},&\text{ for }N\geqslant n_{\mathrm{cr}}+1.\end{cases}$$

(20)

The above captures the tradeoff between the error $\varepsilon_{\mathrm{p}}$ and the prediction horizon $N$: The performance gap decays quadratically in $\varepsilon_{\mathrm{p}}$ and exponentially in $N$. The exponential decay obtained here is analogous to the exponential decay rate in [30] for linear RHC with estimated additive disturbances over the prediction horizon and in [6] for nonlinear RHC, while we have a novel characterization with two distinct decay rates.

The case $N\geqslant n_{\mathrm{cr}}+1$ in (20) has a faster decay rate than the rate for $N\leqslant n_{\mathrm{cr}}$, which indicates the advantage of choosing a horizon³³3In this work, the control horizon and the prediction horizon are equal. larger than the controllability index, e.g., larger than the state dimension for general controllable systems. Selecting a horizon larger than the state dimension is also suggested in [31]. This also suggests that we may choose a smaller horizon for controlling fully-actuated systems, with $\mathrm{rank}(B)=n$ and $n_{\mathrm{cr}}=1$, compared to under-actuated systems.

Given any $F\in\mathbb{S}_{+}^{n}$ and a general controller gain $\mathcal{K}_{\hat{A},\hat{B}}(F)$ in (13), formulated on the approximate model, the following result characterizes the performance gap between $\mathcal{K}_{\hat{A},\hat{B}}(F)$ and the optimal controller gain $K_{\star}$:

The above result is analogous to [19, Thm. 1] and [20, Prop. 7], with extensions to incorporate the additional error term $\|F-P_{\star}\|$. Lemma 4 shows that if the error $\varepsilon_{\mathrm{m}}+\varepsilon$ is sufficiently small, the resulting controller gain $K=\mathcal{K}_{\hat{A},\hat{B}}(F)$ can stabilize the real system $(A_{\star},B_{\star})$ with a suboptimal performance gap.

As the RHC controller (6) is a special case of $\mathcal{K}_{\hat{A},\hat{B}}(F)$ with $F=\mathcal{R}^{(N-1)}_{\hat{A},\hat{B}}(P)$, Lemma 4 can be exploited to analyze the suboptimality of the RHC controller. Particularly, Lemma 4 shows that $J_{K_{\mathrm{RHC}}}-J_{K_{\star}}$ can be upper bounded by a function of the modeling error bound $\varepsilon_{\mathrm{m}}$ and $\varepsilon$ satisfying

$$\big{\|}\mathcal{R}^{(N-1)}_{\hat{A},\hat{B}}(P)-P_{\star}\big{\|}\leqslant\varepsilon.$$

(21)

Therefore, the main challenge is to characterize $\varepsilon$ by investigating the Riccati iterations using the approximate model.

The above problem has been studied in the classical perturbation analysis of the Riccati difference equation [32]; however, the final bound in [32] does not explicitly show the dependence on $N$, $P_{\star}$, and $\varepsilon_{\mathrm{p}}$. In this work, we provide a novel upper bound as in (21) that captures this dependence.

To this end, due to the approximate model, we first define a state-transition matrix similar to (16):

$$\bar{\Phi}^{(j:i)}(P)\triangleq\begin{cases}\bar{\mathcal{L}}\big{(}\mathcal{R}^{(j)}_{A_{\star},B_{\star}}(P)\big{)}\cdots\bar{\mathcal{L}}\big{(}\mathcal{R}^{(i-1)}_{A_{\star},B_{\star}}(P)\big{)},&i>j\geqslant 0,\\ I,&i=j,\end{cases}$$

(22)

where

$$\bar{\mathcal{L}}(P)\triangleq\mathcal{W}(P)[\mathcal{H}(P)+\mathcal{L}_{A_{\star},B_{\star}}(P)],$$

with the functions $\mathcal{W}$ and $\mathcal{H}$ defined on $\mathbb{S}_{+}^{n}$:

$$\mathcal{W}(P)\triangleq I+(S_{B_{\star}}-S_{\hat{B}})P,\text{ }\mathcal{H}(P)\triangleq(I+S_{B_{\star}}P)^{-1}(\hat{A}-A_{\star}).$$

Recall the closed-loop matrix $\mathcal{L}_{A_{\star},B_{\star}}(P)$ in (15), and thus, we can interpret $\bar{\mathcal{L}}(P)$ as a perturbed $\mathcal{L}_{A_{\star},B_{\star}}(P)$: If there is no modeling error, we have $\bar{\mathcal{L}}(P)=\mathcal{L}_{A_{\star},B_{\star}}(P)$ due to $\mathcal{H}(P)=0$ and $\mathcal{W}(P)=I$. Therefore, $\bar{\Phi}^{(j:i)}(P)$ can also be interpreted as a perturbed version of the state-transition matrix $\Phi_{A_{\star},B_{\star}}^{(j:i)}(P)$ due to the approximate model.

Then, as an extension of (17) to the approximate model, the following lemma can be obtained:

Equation (23) is closely related to (17): While (17) concerns the state-transition matrix $\Phi_{A_{\star},B_{\star}}$ of the true system, (23a) contains the perturbed state-transition matrix $\bar{\Phi}$ and the state-transition matrix $\Phi_{\hat{A},\hat{B}}$ of the approximate model; moreover, (23b) contains an extra term caused by the modeling error.

To upper bound $\|\mathcal{R}^{(N-1)}_{\hat{A},\hat{B}}(P)-P_{\star}\|$, we exploit (23) with $P_{1}=P$ and $P_{2}=P_{\star}$. This motivates us to analyze the spectral norm of the perturbed state-transition matrices in (23). The analysis of these state-transition matrices is analogous to the stability analysis of the corresponding time-varying systems.

Then we establish the following results for the two perturbed state-transition matrices.