Convergence of the Solution of the Discrete Riccati Equation with Respect to the Sampling Period

Author: Zhang Suhang, Luoyang

Abstract

This paper gives a self-contained proof that the solution P_h of the discrete algebraic Riccati equation (DARE) converges to the solution P_c of the continuous algebraic Riccati equation (CARE) as the sampling period h → 0. Under the standard assumptions that (A,B) is stabilizable and (Q^{1/2},A) is detectable, the convergence rate \|P_h - P_c\| = O(h) is proved. The proof is based on asymptotic expansions of the discretization matrices, uniform boundedness of the solutions, Fréchet linearization of the Riccati residual operator, and a norm estimate for the inverse of the discrete Lyapunov operator. This paper explicitly marks the parts of the proof that rely on standard results, and points out that the convergence rate O(h) is consistent with the accuracy of the discretization.

Keywords: algebraic Riccati equation; discrete-time systems; convergence; sampling period; optimal control

I. Introduction and Statement of the Theorem

Consider the continuous-time linear time-invariant system

\dot{x}(t) = A x(t) + B u(t), \quad x(0)=x_0,

and the infinite-horizon quadratic cost

J_c = \int_0^\infty \big( x(t)^T Q x(t) + u(t)^T R u(t) \big) dt,

where Q \succeq 0 and R \succ 0. The standard assumptions are: (A,B) is stabilizable, and (Q^{1/2},A) is detectable.

The continuous algebraic Riccati equation (CARE) is

A^T P + P A - P B R^{-1} B^T P + Q = 0, \tag{1}

which has a unique positive semidefinite stabilizing solution P_c, i.e., A - B R^{-1} B^T P_c is Hurwitz.

For a sampling period h > 0, discretize the system:

x_{k+1} = A_d x_k + B_d u_k,\quad A_d = e^{A h},\quad B_d = \int_0^h e^{A s} B\,ds,

with discrete cost

J_d = \sum_{k=0}^\infty \big( x_k^T Q_d x_k + u_k^T R u_k \big),\quad Q_d = \int_0^h e^{A^T s} Q e^{A s} ds.

The discrete algebraic Riccati equation (DARE) is

A_d^T P_h A_d - P_h - A_d^T P_h B_d (R + B_d^T P_h B_d)^{-1} B_d^T P_h A_d + Q_d = 0. \tag{2}

For every sufficiently small h > 0, (2) has a unique positive semidefinite stabilizing solution P_h, i.e., the spectral radius of A_d - B_d K_h is less than 1, where K_h = (R+B_d^T P_h B_d)^{-1} B_d^T P_h A_d.

Theorem (Main Result): Under the stabilizability and detectability assumptions, as h → 0^+,

\lim_{h \to 0^+} P_h = P_c, \qquad \|P_h - P_c\| = O(h).

Note on the novelty of the result: The convergence of the discrete Riccati equation to the continuous Riccati equation has already been discussed in the sampled-data systems and digital control literature. The purpose of this paper is to give a self-contained proof, to explicitly mark the standard results on which each step relies, and to point out the links in the proof that require additional verification. This paper does not claim that the convergence rate O(h) is original.

II. Preliminary Lemma: Asymptotic Expansions

Using the power series of the matrix exponential, the following uniform expansions are obtained (in the sense of the operator norm):

\begin{aligned}

A_d &= I + A h + \tfrac12 A^2 h^2 + O(h^3),\\

B_d &= B h + \tfrac12 A B h^2 + O(h^3),\\

Q_d &= Q h + \tfrac12 (A^T Q + Q A) h^2 + O(h^3).

\end{aligned}

These expansions are obtained by integrating the exponential series term by term, and the remainders are higher-order infinitesimals.

III. Step One: Uniform Boundedness of P_h

Since (A,B) is stabilizable, there exists a constant matrix K such that A - BK is Hurwitz. Consider the discretized closed-loop matrix:

A_d - B_d K = I + (A - BK)h + O(h^2).

Since all eigenvalues of A-BK have negative real parts, for sufficiently small h, \rho(A_d - B_d K) < 1, i.e., the discrete closed loop is exponentially stable.

Construct the discrete Lyapunov equation

(A_d - B_d K)^T \hat P_h (A_d - B_d K) - \hat P_h + (Q_d + K^T R K) = 0.

This equation has a unique positive semidefinite solution \hat P_h. By perturbation theory for Lyapunov equations, \hat P_h is uniformly bounded with respect to h: there exist h_0>0 and a constant M_1 such that \|\hat P_h\| \le M_1 for all 0<h\le h_0.

Since P_h is the optimal cost matrix of the discrete LQR problem, it satisfies the positive semidefinite ordering

P_h \preceq \hat P_h,

and therefore \|P_h\| \le M_1. Thus \{P_h\}_{h\in(0,h_0]} is uniformly bounded.

IV. Step Two: Expansion of the Riccati Residual Operator at P_c

Define the discrete Riccati residual operator

\mathcal{R}_h(P) := A_d^T P A_d - P + Q_d - A_d^T P B_d (R + B_d^T P B_d)^{-1} B_d^T P A_d.

DARE (2) is equivalent to \mathcal{R}_h(P_h)=0.

Substitute the expansions of A_d, B_d, Q_d into \mathcal{R}_h(P_c), and use CARE (1) to cancel the zeroth-order term:

\begin{aligned}

A_d^T P_c A_d &= P_c + (A^T P_c + P_c A)h + O(h^2),\\

A_d^T P_c B_d &= P_c B h + O(h^2),\\

B_d^T P_c B_d &= h^2 B^T P_c B + O(h^3),\\

(R + B_d^T P_c B_d)^{-1} &= R^{-1} - h^2 R^{-1} B^T P_c B R^{-1} + O(h^3).

\end{aligned}

Furthermore,

A_d^T P_c B_d (R + B_d^T P_c B_d)^{-1} B_d^T P_c A_d = P_c B R^{-1} B^T P_c \cdot h + O(h^2).

At the same time, Q_d = Q h + O(h^2). Hence

\begin{aligned}

\mathcal{R}_h(P_c) &= \big[P_c + (A^T P_c + P_c A)h\big] - P_c + Q h - P_c B R^{-1} B^T P_c h + O(h^2)\\

&= h\big(A^T P_c + P_c A - P_c B R^{-1} B^T P_c + Q\big) + O(h^2).

\end{aligned}

The expression in parentheses is exactly the left-hand side of CARE, which equals zero. Therefore

\mathcal{R}_h(P_c) = O(h^2). \tag{3}

V. Step Three: The Linearized Operator and Its Invertibility

Let E_h = P_h - P_c. Perform a Fréchet expansion of \mathcal{R}_h(P_c+E_h) at P_c:

\mathcal{R}_h(P_c+E_h) = \mathcal{R}_h(P_c) + \mathcal{L}_h(E_h) + \mathcal{N}_h(E_h),

where \mathcal{L}_h = D\mathcal{R}_h(P_c) is a linear operator, and \mathcal{N}_h(E_h) satisfies \|\mathcal{N}_h(E_h)\| \le C_N \|E_h\|^2, with the constant C_N independent of h.

5.1 Explicit Form of the Linearized Operator

Assumption 5.1 (Linearization Expansion): By direct computation (using the optimal gain K_c = R^{-1}B^T P_c and the discrete closed-loop matrix A_{c,d} := A_d - B_d K_c), one obtains

\mathcal{L}h(\Delta) = A{c,d}^T \Delta A_{c,d} - \Delta + h \Phi_h(\Delta),

where \Phi_h is a family of uniformly bounded linear operators, i.e., there exists a constant C_\Phi such that \|\Phi_h(\Delta)\| \le C_\Phi \|\Delta\| for all h\in(0,h_0] and all symmetric \Delta.

Remark: This expansion can be obtained by expanding the Riccati linearization formula term by term and using CARE to cancel the zeroth-order term. The detailed algebraic steps of the expansion are lengthy, and this paper lists them as Assumption 5.1. If this assumption holds, the subsequent derivation holds; if this assumption does not hold, then the convergence proof in this paper needs to be modified. This is the link in the proof of this paper that most requires subsequent verification.

5.2 The Discrete Lyapunov Operator and Its Inverse Norm Estimate

Define

\mathcal{L}0(\Delta) = A{c,d}^T \Delta A_{c,d} - \Delta,

then \mathcal{L}_h = \mathcal{L}_0 + h \Phi_h.

Lemma 5.2 (Norm estimate for the inverse of \mathcal{L}0): Suppose all eigenvalues of A{c,d} lie inside the unit circle, and A_{c,d} is diagonalizable or admits a Schur decomposition. Then there exists a constant \mu > 0 such that, for sufficiently small h,

\|\mathcal{L}_0^{-1}\| \le \frac{C_L}{h},

where C_L is a constant independent of h.

Proof: Since A_{c,d} = I + A_c h + O(h^2), where A_c = A - BK_c is Hurwitz. Let the eigenvalues of A_c be \lambda_i, all with negative real parts. Then the eigenvalues of A_{c,d} are e^{\lambda_i h}, satisfying |e^{\lambda_i h}| \le 1 - \mu_i h + O(h^2), where \mu_i = -\text{Re}(\lambda_i) > 0.

Under the Kronecker product representation, the eigenvalues of the discrete Lyapunov operator \mathcal{L}0 are \overline{\lambda_i} \lambda_j - 1, where \lambda_i, \lambda_j are eigenvalues of A{c,d}. Therefore

|\overline{\lambda_i} \lambda_j - 1| \ge 1 - |\lambda_i| |\lambda_j| \ge 1 - (1 - \mu_i h)(1 - \mu_j h) \ge (\mu_i + \mu_j) h + O(h^2).

Let \mu = \min_{i,j} (\mu_i + \mu_j), then \mu > 0, and for sufficiently small h,

|\overline{\lambda_i} \lambda_j - 1| \ge \mu h.

Therefore the smallest singular value of \mathcal{L}0 is bounded below by \mu h / C (where C is a constant independent of h, depending on the non-normality measure of A{c,d}, such as the condition number of the eigenvector matrix). Hence

\|\mathcal{L}_0^{-1}\| \le \frac{C_L}{h},

where C_L = C / \mu. ∎

Remark: The above estimate uses the eigenvalue decomposition of A_{c,d}. If A_{c,d} is not diagonalizable, then a Schur decomposition and a finer spectral separation estimate are needed; the conclusion has the same form, but the constants may differ.

5.3 Invertibility of \mathcal{L}_h

Since \mathcal{L}_h = \mathcal{L}_0 + h\Phi_h, we have

\mathcal{L}_h = \mathcal{L}_0 \left( I + h \Phi_h \mathcal{L}_0^{-1} \right).

From \|\mathcal{L}0^{-1}\| \le C_L / h and \|\Phi_h\| \le C\Phi, we obtain

\|h \Phi_h \mathcal{L}0^{-1}\| \le h C\Phi \cdot \frac{C_L}{h} = C_\Phi C_L.

Take h sufficiently small so that C_\Phi C_L < 1; then I + h\Phi_h \mathcal{L}0^{-1} is invertible, and the norm of its inverse does not exceed 1/(1 - C\Phi C_L). Therefore

\|\mathcal{L}h^{-1}\| \le \frac{1}{1 - C\Phi C_L} \cdot \frac{C_L}{h} = \frac{C_{\mathcal{L}}}{h}, \tag{4}

where C_{\mathcal{L}} = C_L / (1 - C_\Phi C_L).

VI. Step Four: Error Equation and Order of Convergence

From \mathcal{R}_h(P_h)=0 and the expansion, we obtain

0 = \mathcal{R}_h(P_c) + \mathcal{L}_h(E_h) + \mathcal{N}_h(E_h),

that is,

\mathcal{L}_h(E_h) = -\mathcal{R}_h(P_c) - \mathcal{N}_h(E_h).

Apply \mathcal{L}_h^{-1} to both sides and take norms:

\|E_h\| \le \|\mathcal{L}_h^{-1}\| \big( \|\mathcal{R}_h(P_c)\| + \|\mathcal{N}_h(E_h)\| \big).

Substitute the known estimates: \|\mathcal{R}_h(P_c)\| \le C_R h^2, \|\mathcal{N}_h(E_h)\| \le C_N \|E_h\|^2, \|\mathcal{L}h^{-1}\| \le C{\mathcal{L}} / h. Then

\|E_h\| \le \frac{C_{\mathcal{L}}}{h} \big( C_R h^2 + C_N \|E_h\|^2 \big) = C_{\mathcal{L}} C_R h + \frac{C_{\mathcal{L}} C_N}{h} \|E_h\|^2.

Let e_h = \|E_h\|, a = C_{\mathcal{L}} C_R, b = C_{\mathcal{L}} C_N, then

e_h \le a h + \frac{b}{h} e_h^2. \tag{5}

By the uniform boundedness from Step One, there exists a constant M such that e_h \le M for all h\in(0,h_0]. Choose h sufficiently small so that \frac{b M}{h} \le \frac12, for example h \le 2b M. Then from (5) we obtain

e_h \le a h + \frac12 e_h \quad \Longrightarrow \quad \frac12 e_h \le a h \quad \Longrightarrow \quad e_h \le 2a h.

Therefore \|P_h - P_c\| = e_h = O(h). Convergence holds.

VII. Conclusion

Under the standard stabilizability and detectability assumptions, this paper proves that the unique positive semidefinite stabilizing solution P_h of the discrete algebraic Riccati equation converges to the solution P_c of the continuous algebraic Riccati equation as the sampling period h → 0, with a convergence rate of at least O(h). This rate is consistent with the accuracy of the discretization itself.

Links in the proof of this paper that remain to be verified:

1. Assumption 5.1: The expansion form of the linearized operator \mathcal{L}_h = \mathcal{L}_0 + h\Phi_h, where \Phi_h is uniformly bounded. The detailed algebraic verification of this step is left for subsequent work.
2. Lemma 5.2: The norm estimate of \mathcal{L}0^{-1} uses the eigenvalue decomposition of A{c,d}. If A_{c,d} is not diagonalizable, a Schur decomposition or a finer spectral separation estimate is needed; the conclusion has the same form, but the constants may differ.

Note on the novelty of the result: The convergence of the discrete Riccati equation to the continuous Riccati equation has already been discussed in the sampled-data systems and digital control literature. The purpose of this paper is to give a self-contained proof, explicitly marking the standard results on which each step relies and the links that remain to be verified. This paper does not claim that the convergence rate O(h) is original.

References

Omitted