Policy iteration (PI) is a widely used Newton-type method for continuous-time linear quadratic regulator (LQR) problems [1, 2, 3]. PI is initialized with a stabilizing feedback gain and proceeds by repeatedly evaluating the current policy and updating the control law. The iteration preserves closed-loop stability, generates a monotonically nonincreasing sequence of value matrices, and converges locally at the Newton rate [4, 5]. Integral and off-policy reinforcement learning (RL) formulations implement both steps using measured state–input data without explicitly identifying the system matrices [6, 7, 8, 9, 10, 11]. Together, these properties make PI attractive for data-driven optimal control. Yet standard infinite-horizon PI is not self-starting [12, 13, 14]. The policy-evaluation step in infinite-horizon PI requires the current feedback gain to stabilize the closed-loop system. If the gain is not stabilizing, closed-loop trajectories need not decay, the infinite-horizon cost is generally not finite, and the associated Lyapunov equation no longer defines a valid value matrix. Consequently, the set of admissible policies defines the domain on which infinite-horizon policy evaluation admits a valid value-function interpretation. When the system matrices are unknown, constructing and certifying an admissible gain may require precisely the model information that data-driven PI seeks to avoid. This circular dependence poses a structural obstacle to implementing infinite-horizon PI.

Value iteration (VI) avoids the need for an initial stabilizing feedback gain by initializing a value matrix instead [15, 16, 17]. In this setting, existing continuous-time VI schemes typically rely on diminishing or small step sizes and may use bounded-set projections or reset rules [18, 19]. Their first-order updates may also converge more slowly than the local Newton iteration underlying PI. Hybrid iteration (HI) uses VI to reach the admissible-policy region and then switches to PI to exploit its fast local convergence [20, 21, 22]. The HI framework has also been applied to Markov jump systems [23, 24], networked control systems [25], and autonomous vehicles [26]. Although HI reduces reliance on an initial stabilizing feedback gain, its VI stage still requires step-size selection, boundedness verification, and an explicit VI-to-PI switching criterion.

A second class of methods modifies the policy-evaluation operator to enlarge the initialization region. Bias PI introduces a bias parameter into policy evaluation and transfers the stability requirement from the original closed-loop matrix to a shifted auxiliary matrix [27]. A smaller bias yields PI-like behavior, whereas a larger bias can expand the initialization region at the cost of making the iteration more VI-like. Related ideas have been explored for discrete-time systems [31] and zero-sum games [32, 33]. Bias PI thus relaxes the initial admissibility requirement while retaining linear matrix equations at each iteration. Implementing the method requires selecting a sufficiently large bias to ensure that auxiliary policy evaluation is well posed. Global variants additionally use boundedness tests and bias-adjustment rules.

Homotopy-based methods take a different approach. They start from an artificially stabilized system that admits a trivial or readily constructed stabilizing policy. They then gradually deform the auxiliary system into the original plant while maintaining a stabilizing policy for each intermediate problem. Homotopic PI methods move the artificial dynamics toward those of the original plant using a cumulative continuation factor and provide an explicit bound on the closed-loop pole location [28, 30]. A more recent predictor-corrector homotopy method constructs a path of stabilizing Riccati solutions and uses first-order differential information to predict the controller for the next auxiliary problem, thereby enabling faster continuation along the path [29]. Related homotopy-based methods have subsequently been extended to nonlinear systems [34], inverse RL [35], multiagent systems [36], and Markov jump systems [37]. These methods substantially relax the requirement for an initially stabilizing policy for the original plant. Nevertheless, each continuation stage still involves an infinite-horizon subproblem, and the corresponding intermediate policy must be stabilizing to ensure that policy evaluation is well posed. Thus, homotopy methods address the initialization difficulty by constructing a stability-preserving path through a sequence of auxiliary infinite-horizon problems. Table 1 summarizes the initialization requirements of the representative continuous-time data-driven RL methods discussed above and highlights the distinctions between these methods and the proposed finite-horizon bootstrap method.

The difficulty of evaluating a nonstabilizing policy stems from the unbounded evaluation interval. By contrast, over a compact interval, the terminal-value policy-evaluation equation is well posed for any bounded time-varying feedback policy, regardless of closed-loop stability. This distinction suggests using finite-horizon PI as a bootstrap stage rather than as a replacement for the original infinite-horizon problem. The finite-horizon stage generates the stabilizing policy needed to initialize infinite-horizon PI. However, implementing this bootstrap from data raises two main challenges. First, the time-varying finite-horizon value matrix and policy must be reconstructed using trajectories of the original plant alone, without collecting data from the auxiliary shifted system. Second, approximation and regression errors inevitably affect this reconstruction. The resulting constant gain must therefore be independently certified to stabilize the original plant before being used to initialize infinite-horizon PI. The main contributions are as follows:

• We introduce a finite-horizon bootstrap mechanism that eliminates the need for an initially stabilizing policy. Unlike initialization methods based on bias PI [27], homotopy [28, 30], or VI [20, 21, 22], our method allows arbitrary bounded intermediate policies because it performs policy evaluation over a compact time interval.
• We establish a finite-to-infinite-horizon connection showing that the initial-time gain generated by the finite-horizon problem approaches the optimal gain of a shifted infinite-horizon LQR problem as the horizon increases. This connection provides a principled route to a prescribed feasible stability margin.
• We develop a data-driven off-policy implementation based solely on trajectories of the original plant and characterize the effects of basis-approximation and regression errors through a perturbed finite-horizon policy-improvement recursion. We further introduce an independent data-driven Lyapunov certificate to verify the admissibility of the candidate gain before using the gain to initialize infinite-horizon PI.

Notation: Throughout this article, we use \( \mathbb{R}^{n}\) for the \( n\) -dimensional real Euclidean space and \( \mathbb{R}^{m\times n}\) for the space of real \( m\times n\) matrices. The identity matrix of order \( n\) is written as \( I_n\) , whereas \( \mathbf{0}\) represents a zero scalar, vector, or matrix whose dimensions are inferred from context. For a matrix \( \textbf{A}\) , we use \( \sigma_{\min}(\textbf{A})\) to represent its smallest singular value, while \( \textbf{A}^{\dagger}\) refers to its Moore–Penrose pseudoinverse. When \( \textbf{A}\) is square, \( \sigma(\textbf{A})\) represents its spectrum, and its spectral abscissa is given by \( \alpha(\textbf{A}) = \max_{\lambda\in\sigma(\textbf{A})} \operatorname{Re}(\lambda) \) , where \( \operatorname{Re}(\lambda)\) represents the real part of \( \lambda\) . Moreover, \( \mathbb{C}_{-} = \{z\in\mathbb{C}:\operatorname{Re}(z)<0\} \) denotes the open left-half complex plane. The minimum and maximum eigenvalues of a symmetric matrix \( \textbf{X}\) are denoted by \( \lambda_{\min}(\textbf{X})\) and \( \lambda_{\max}(\textbf{X})\) , respectively. For matrix-valued function \( \textbf{H}\) defined on \( [0,T]\) , let \( \|\textbf{H}\|_{\infty,T} = \sup_{t\in[0,T]}\|\textbf{H}(t)\|_F . \) Given a matrix \( \textbf{A}=[a_{ij}]\in\mathbb{R}^{m\times n}\) , \( \operatorname{vec}(\textbf{A})\) is formed by stacking the columns of \( \textbf{A}\) into a single vector. For a symmetric matrix \( \textbf{P}=[p_{ij}]\in\mathbb{R}^{n\times n}\) , define \( \operatorname{vecs}(\textbf{P})=[p_{11},2p_{12},…,2p_{1n},p_{22},2p_{23},…,p_{nn}]^T.\) For a vector \( x=[x_1,…,x_n]^T\in\mathbb{R}^{n}\) , define \( \operatorname{vecv}(x)=[x_1^2,x_1x_2,…,x_1x_n,x_2^2,x_2x_3,…,x_n^2]^T.\) For matrices \( X_1,…,X_q\) having the same number of columns, \( \operatorname{col}(X_1,…,X_q) = [X_1^{T},…,X_q^{T}]^{T}\) denotes their vertical concatenation.

Consider the continuous-time linear system

where \( x(t)\in \mathbb{R}^{n}\) is the system state and \( u(t)\in \mathbb{R}^{m}\) is the control signal. The system matrices \( \mathbf{A}\in\mathbb{R}^{n\times n}\) and \( \mathbf{B}\in\mathbb{R}^{n\times m}\) are constant but unavailable. The infinite-horizon quadratic cost is

where \( \mathbf{Q}=\mathbf{Q}^{T}\geq \mathbf{0}\) and \( \mathbf{R}=\mathbf{R}^{T}>\mathbf{0}\) are the state and control penalties, respectively. The pair \( (\mathbf{A},\mathbf{B})\) is assumed to be stabilizable, and \( (\mathbf{A},\mathbf{Q}^{1/2})\) is assumed to be detectable.

When \( \mathbf{A}\) and \( \mathbf{B}\) are available, the optimal control law is

where \( \mathbf{P}^{*}=(\mathbf{P}^{*})^{T}\geq\mathbf{0}\) is the unique stabilizing solution of the algebraic Riccati equation

The standard continuous-time PI for solving (5) is Kleinman’s iteration. The following standard result of Kleinman’s iteration is recalled.

Lemma 1 shows that Kleinman’s iteration is not self-starting, because its well-posedness and stability-preserving property both require an admissible initial gain. The problem addressed in this paper is therefore to construct, from measured data, a constant gain \( \mathbf{K}^{[0]} \) satisfying (3), which can serve as the initializer for infinite-horizon data-driven PI.

2.1 Why a finite-horizon bootstrap method

The central difficulty in removing the initial admissible controller is that the infinite-horizon policy evaluation operator is intrinsically defined only on stabilizing policies. For a fixed gain \( \mathbf{K}\) , the evaluation step requires the Lyapunov equation

\[ \begin{aligned} (\mathbf{A}-\mathbf{B}\mathbf{K})^{T}\mathbf{P}_{K} +\mathbf{P}_{K}(\mathbf{A}-\mathbf{B}\mathbf{K}) +\mathbf{Q} +\mathbf{K}^{T}\mathbf{R}\mathbf{K} =\mathbf{0}. \end{aligned} \]

When \( \mathbf{A}-\mathbf{B}\mathbf{K}\) is Hurwitz, its solution admits the value representation

\[ \mathbf{P}_{K} = \int_{0}^{\infty} e^{(\mathbf{A}-\mathbf{B}\mathbf{K})^{T}t} \left(\mathbf{Q}+\mathbf{K}^{T}\mathbf{R}\mathbf{K}\right) e^{(\mathbf{A}-\mathbf{B}\mathbf{K})t} \,\mathrm{d}t. \]

If \( \mathbf{K}\) is not stabilizing, the state-transition matrix \( e^{(\mathbf{A}-\mathbf{B}\mathbf{K})t}\) does not decay exponentially. Since the cost is accumulated over the unbounded interval \( [0,\infty)\) , the integral is generally not finite, and the evaluation no longer defines a valid cost matrix.

As shown in Fig. 1, a finite-horizon evaluation avoids this obstruction because it is posed on a compact interval. For any bounded time-varying gain \( \mathbf{K}(t)\) on \( [0,T] \) , the terminal-value Lyapunov equation

\[ \begin{equation} \begin{aligned} -\dot{\mathbf P}(t) =&\,(\mathbf A-\mathbf B\mathbf K(t))^T\mathbf P(t) +\mathbf P(t)(\mathbf A-\mathbf B\mathbf K(t)) \\ &+\mathbf Q+\mathbf K^T(t)\mathbf R\mathbf K(t), ~~ \mathbf P(T)=\mathbf M \end{aligned} \end{equation} \]

(8)

has a unique solution on \( [0,T]\) with a terminal weight \( \textbf{M}=\textbf{M}^{T}\geq\textbf{0}\) . This well-posedness does not require the time-varying closed-loop system to be exponentially stable. Thus, unlike the infinite-horizon evaluation, the finite-horizon evaluation is well defined for arbitrary bounded policies and can therefore be used before an admissible gain is available.

The finite-horizon stage must still produce a gain that is admissible for the original infinite-horizon problem. To this end, the evaluation is performed for the shifted dynamics \( \dot x(t)=\mathbf{A}_{\rho }x(t)+\mathbf{B}u(t)\) with \( \mathbf{A}_{\rho }=\mathbf{A}+\rho \mathbf{I}\) and \( \rho>\varepsilon\) . Let \( \mathbf{K}_{\rho,T}^\ast(t)\) denote the optimal finite-horizon gain of this shifted problem. As \( T\) increases, the initial-time gain \( \mathbf{K}_{\rho,T}^\ast(0)\) approaches the stabilizing infinite-horizon gain \( \mathbf{K}_\rho^\ast=\mathbf{R}^{-1}\mathbf{B}^T\mathbf{P}_\rho^\ast , \) where \( \mathbf{P}_\rho^\ast\) solves the shifted algebraic Riccati equation. Since \( \mathbf{A}_\rho-\mathbf{B}\mathbf{K}_\rho^\ast\) is Hurwitz, applying the same gain to the original system gives

\[ \alpha(\mathbf{A}-\mathbf{B}\mathbf{K}_\rho^\ast) = \alpha(\mathbf{A}_\rho-\mathbf{B}\mathbf{K}_\rho^\ast)-\rho <-\rho<-\varepsilon . \]

Therefore, by continuity of the spectral abscissa, \( \mathbf{K}_{\rho,T}^\ast(0)\) is also \( \varepsilon\) -admissible for all sufficiently large \( T\) .

Pan/Zoom

This section presents the model-based version of the proposed finite-horizon bootstrap method. The model-based construction clarifies the role of the finite horizon and serves as the reference algorithm for the data-driven development.

Let \( \varepsilon \geq 0\) be the desired stability margin in (3). Choose a design parameter \( \rho >\varepsilon \) and define the shifted system matrix \( \mathbf{A}_{\rho }\) . The following assumption is imposed on the shifted pair.

Assumption 1. The pair \( (\mathbf{A}_{\rho},\mathbf{B})\) is stabilizable, and \( (\mathbf{A}_{\rho},\mathbf{Q}^{1/2})\) is detectable.

Assumption 1 guarantees the existence of the stabilizing solution \( \mathbf{P}_{\rho }^{\ast }\) to the shifted algebraic Riccati equation

The control gain is \( \mathbf{K}_{\rho }^{\ast }=\mathbf{R}^{-1}\mathbf{B}^{T}\mathbf{P}_{\rho }^{\ast }\) . Since \( \mathbf{A}_{\rho }-\mathbf{B}\mathbf{K}_{\rho }^{\ast }\) is Hurwitz, the matrix \( \mathbf{A}-\mathbf{B}\mathbf{K}_{\rho }^{\ast }=\mathbf{A}_{\rho }-\mathbf{B}\mathbf{K}_{\rho }^{\ast }-\rho \mathbf{I}\) is Hurwitz with a decay margin greater than \( \rho \) . The purpose of the finite-horizon stage is to approximate this gain without requiring an admissible initial policy.

For a horizon \( T>0\) and a terminal matrix \( \mathbf{M}=\mathbf{M}^{T}\geq \mathbf{0}\) , consider the shifted finite-horizon Riccati differential equation

with terminal condition \( \mathbf{P}_{\rho ,T}^{\ast }(T)=\mathbf{M}\) . The associated finite-horizon optimal gain is

The finite-horizon gain in (11) can be computed by a finite-horizon PI without requiring a stabilizing initial policy. Let \( \mathbf{K}_{\rho ,T}^{[0]}(t)\) be any bounded time-varying matrix on \( [0,T]\) . For \( j=0,1,2,… \) , solve

with \( \mathbf{P}_{\rho ,T}^{[j]}(T)=\mathbf{M}\) , and update

Because (12) is solved over a compact interval, the initial policy \( \mathbf{K}_{\rho ,T}^{[0]}(t)\) need not stabilize the shifted system. As will be shown below, the finite-horizon PI sequence converges to \( \mathbf{K}_{\rho ,T}^{\ast }(t)\) . After convergence, with \( \bar{j}\) denoting the terminal iteration index, the initializer is taken as \( \mathbf{K}^{[0]}=\mathbf{K}_{\rho ,T}^{[\bar{j}]}(0)\) . In the model-based setting, admissibility can be checked directly by verifying

If (14) holds, \( \mathbf{K}^{[0]}\) is used to initialize Kleinman’s infinite-horizon PI (6)–(7) for the system (1). Otherwise, after the finite-horizon PI has converged for the current horizon, the horizon is updated according to

Algorithm 1 summarizes the model-based finite-horizon bootstrap procedure. The following theorem justifies the algorithm by separating the fixed-horizon finite-horizon PI property from the admissibility guarantee obtained by enlarging the shifted horizon.

The multiplicative update (15) gives \( T_{s}=\gamma ^{s}T_{0}\) . Since \( \gamma >1\) , \( T_{s}\rightarrow \infty \) as \( s\rightarrow \infty \) . Theorem 2 implies that the horizon is sufficiently large if

then \( T_{s}\geq T_{\varepsilon }\) . Therefore, for every \( s\geq s_{\varepsilon }\) , there exists \( j_{\varepsilon }(T_{s})\) such that

Thus, the finite-horizon bootstrap mechanism ensures that an \( \varepsilon\) -admissible initializer can be obtained after a sufficiently large horizon and sufficiently many finite-horizon PI iterations. In the model-based algorithm, this admissibility is verified directly by the spectral test (14).

This section develops a data-driven realization of the finite-horizon PI stage in Algorithm 1. The data are collected from the original system \( \dot{x}(t)=\mathbf{A}x(t)+\mathbf{B}u_{b}(t)\) , where \( u_{b}(t)\) is a behavior input.

For the current policy \( u(t)=-\mathbf{K}_{\rho ,T}^{[j]}(t)x(t)\) , define \( \mathbf{Q}_{\rho ,T}^{[j]}(t)=\mathbf{Q}+(\mathbf{K}_{\rho ,T}^{[j]}(t))^{T}\mathbf{R}\mathbf{K}_{\rho ,T}^{[j]}(t)\) and \( \mathbf{z}^{[j]}(t)=u_{b}(t)+\mathbf{K}_{\rho ,T}^{[j]}(t)x(t).\) The signal \( \mathbf{z}^{[j]}(t)\) is the off-policy correction term. Let \( V^{[j]}(x,t)=x^{T}\mathbf{P}_{\rho ,T}^{[j]}(t)x\) . Along the behavior trajectory, one obtains

Integrating (25) over each stored interval \( I_k=[t_k,t_{k+1}]\subset[0,T]\) , \( k=1,…,N_d\) , gives

Choose differentiable basis functions \( \phi _{\ell }(t)\) , \( \ell =1,… ,N_{p}\) with \( \phi _{\ell }(T)=0\) and basis functions \( \psi _{\ell }(t)\) , \( \ell =1,… ,N_{k}\) . Approximate \( \mathbf{P}_{\rho,T}^{[j]}(t)\) and \( \mathbf{K}_{\rho,T}^{[j+1]}(t)\) as

where \( \mathbf{P}_{\ell }^{[j]}=(\mathbf{P}_{\ell }^{[j]})^{T}\) , and

where \( \Delta _{P}^{[j]}(t)\) and \( \Delta _{K}^{[j+1]}(t)\) denote the approximation residuals induced by the chosen basis functions.

To facilitate the solution of (26) using the representations in (27)-(28), we define the following interval-wise quantities for each stored interval \( I_{k}\)

Then (26) is equivalent to

where

Stacking all stored intervals gives

where

where \( \boldsymbol{\Omega}_{\rho,T}^{[j]} \in\mathbb{R}^{N_{d}\times N_{\theta}}\) , \( \theta^{[j]}\in\mathbb{R}^{N_{\theta}}\) , and \( \mathbf{b}_{\rho,T}^{[j]}\in\mathbb{R}^{N_{d}}\) , with \( N_{\theta} = N_{p}\frac{n(n+1)}{2}+N_{k}mn\) . The residual \( \mathbf{r}_{\rho,T}^{[j]}\) is induced by the approximation errors associated with the selected basis functions. The following lemma provides an explicit upper bound on this term.

Lemma 4 allows the data-driven update to be viewed as an inexact finite-horizon PI step. For a fixed horizon \( T>0\) , define the finite-horizon policy-improvement mapping as follows. For any bounded matrix-valued function \(\mathbf K:[0,T]\to\mathbb R^{m\times n}\), let \( \mathbf{P}_{\rho ,T}(t;\mathbf{K})\) be the unique solution on \([0,T]\) of

with terminal condition \( \mathbf{P}_{\rho ,T}(T;\mathbf{K})=\mathbf{M}\) . Define

Then the exact finite-horizon PI recursion can be written as \( \mathbf{K}_{\rho ,T}^{[j+1]}=\mathcal{F}_{\rho ,T}(\mathbf{K}_{\rho ,T}^{[j]})\) . Moreover, by Lemma 4, the data-driven finite-horizon recursion has the inexact form

Corollary 6 gives a theoretical admissibility guarantee, but the neighborhood radius involved there is model dependent and cannot be checked directly from data. Therefore, we introduce a data-driven Lyapunov certificate for the candidate constant gain generated by the finite-horizon stage. Let \( \mathbf{K}_{c}=\widehat{\mathbf{K}}_{\rho ,T}^{[j+1]}(0)\) be a candidate constant gain generated by the data-driven finite-horizon recursion. Choose a matrix \( \mathbf{W}_{c}=\mathbf{W}_{c}^{T}>\mathbf{0}\) . The purpose of the following test is to determine, using the stored behavior data only, whether there exists a matrix \( \mathbf{P}_{c}=\mathbf{P}_{c}^{T}>\mathbf{0}\) such that \( (\mathbf{A}-\mathbf{B}\mathbf{K}_{c}+\varepsilon \mathbf{I})^{T}\mathbf{P}_{c}+\mathbf{P}_{c}(\mathbf{A}-\mathbf{B}\mathbf{K}_{c}+\varepsilon \mathbf{I})=-\mathbf{W}_{c}\) . If such a matrix \( \mathbf{P}_{c}\) exists, then \( \mathbf{K}_{c}\) is \( \varepsilon \) -admissible for the original system.

Define \( \mathbf{L}_{c}=\mathbf{B}^{T}\mathbf{P}_{c} \) and \( z_{c}(t)=u_{b}(t)+\mathbf{K}_{c}x(t)\) . The desired Lyapunov equation is equivalent to

Along the trajectory \( \dot{x}(t)=\mathbf{A}x(t)+\mathbf{B}u_{b}(t),\) one has

Integrating (57) over \( I_{k}=[t_{k},t_{k+1}]\) gives

For each interval \( I_{k}\) , define

Then (58) can be written as

Stacking all stored intervals yields

where

with

The following result gives a data-driven admissibility certificate for the candidate gain.

The data-driven bootstrap scheme is summarized in Algorithm 2 and illustrated in Fig. 2. It consists of two nested loops: an inner finite-horizon PI loop that generates a candidate gain at a fixed horizon, and an outer horizon-expansion loop that is activated whenever the candidate fails the data-driven admissibility test. Once the certificate succeeds, the candidate gain is returned as an initializer for the subsequent infinite-horizon PI.

Pan/Zoom

5.1 Case study of the batch reactor system

In this section, a batch reactor system [38] is used to illustrate the proposed finite-horizon bootstrap method. The system matrices are given by

\[ \mathbf{A} = \begin{bmatrix} 1.3800 & -0.2077 & 6.7150 & -5.6760 \\ -0.5814 & -4.2900 & 0 & 0.6750 \\ 1.0670 & 4.2730 & -6.6540 & 5.8930 \\ 0.0480 & 4.2730 & 1.3430 & -2.1040 \end{bmatrix}, \]

\[ \mathbf{B} = \begin{bmatrix} 0 & 0 \\ 5.6790 & 0 \\ 1.1360 & -3.1460 \\ 1.1360 & 0 \end{bmatrix}. \]

The open-loop poles are \( 1.9910\) , \( 0.0635\) , \( -8.6659\) , and \( -5.0566\) , which shows that the uncontrolled system is unstable. The weighting matrices are selected as \( \mathbf{Q}= \rm {diag} (1,1,0.1,0.01)\) , \( \mathbf{R}=I_{2}\) , and the terminal weight is chosen as \( \mathbf{M}=\mathbf{0}\) . We first validate the model-based finite-horizon bootstrap procedure. The design parameters are set as \( \varepsilon =0.2\) , \( \rho =0.6\) , \( T_{0}=0.1\) , \( \epsilon_{\mathrm{PI}}=10^{-7}\) and \( \gamma =1.25\) . The finite-horizon policy is initialized by \( \mathbf{K}_{\rho ,T_{s}}^{[0]}(t)=\mathbf{0}\) for \( t\in \lbrack 0,T_{s}]\) , and the horizon is updated according to \( T_{s+1}=\gamma T_{s}\) .

The numerical results are summarized in Table 2 and Fig. 3. It can be seen that, as the finite horizon \( T_s\) increases, the spectral abscissa of the closed-loop matrix \( \mathbf{A}-\mathbf{B}\mathbf{K}_{c}\) decreases monotonically in this example. When \( T_{s}=0.3052\) , the finite-horizon bootstrap produces a gain satisfying \( \alpha (\mathbf{A}-\mathbf{B}\mathbf{K}_{c})<-0.2. \) Therefore, the algorithm terminates and the obtained gain is accepted as an \( \varepsilon \) -admissible initializer.

Table 2. Model-based finite-horizon bootstrap results for the batch reactor system

\( s\)	\( T_s\)	Iterations	\( \alpha(\mathbf{A}-\mathbf{B}\mathbf{K}_c)\)
0	0.1000	4	1.8108
1	0.1250	4	1.6755
2	0.1563	4	1.4501
3	0.1953	5	1.0694
4	0.2441	5	0.3866
5	0.3052	5	-0.5822

The resulting stabilizing initializer is

\[ \mathbf{K}_{c} = \begin{bmatrix} 0.1086 & 0.5935 & 0.1185 & 0.0401 \\ -1.0155 & -0.0248 & -0.6106 & 0.4062 \end{bmatrix}. \]

The corresponding closed-loop poles are \( -0.5822\pm 0.5852i\) , \( -8.3664\) , and \( -7.6089\) , where \( i\) represents the imaginary unit, which confirms that the proposed finite-horizon bootstrap successfully generates an admissible initial policy.

Pan/Zoom

We next evaluate the data-driven bootstrap procedure in Algorithm 2. The certificate matrix is chosen as \( \mathbf{W}_{c}=I_{4}\) . The initial states are randomly generated with norms in \( [0.3,0.7]\) . The basis functions (27)–(28) are selected as

\[ \begin{align*} \phi_{\ell}(t) &= \left(1-\frac{t}{T_{s}}\right) \left(\frac{t}{T_{s}}\right)^{\ell-1},\\ \psi_{\ell}(t) &= \left(\frac{t}{T_{s}}\right)^{\ell-1}, ~~ \ell=1,\ldots,8. \end{align*} \]

Thus \( N_{p}=N_{k}=8\) . The data-driven PI loop is initialized by \( \widehat{\mathbf{K}}_{\rho ,T_{s}}^{[0]}(t)=\mathbf{0}\) . The data-driven results are shown in Fig. 4. For all tested horizons, the data-richness condition for the Lyapunov certificate is satisfied with \( \rm {rank} =18\) , and the finite-horizon regression matrix has full column rank, i.e., \( \rm {rank} (\boldsymbol{\Omega}_{\rho,T_s}^{[j]})=144\) . For numerical implementation, the linear equation (60) is solved in the least-squares sense. The relative residual

\[ \begin{equation} \delta_{c} = \frac{ \left\| \widetilde{\boldsymbol{\Psi}}_{c}\widehat{\theta}_{c} +\widetilde{\mathbf{b}}_{c} \right\|_{2} }{ \max\left\{ 1, \left\|\widetilde{\mathbf{b}}_{c}\right\|_{2} \right\} } \end{equation} \]

(65)

is reported only as a numerical accuracy indicator, where \( \widetilde{\boldsymbol{\Psi}}_c\) and \( \widetilde{\mathbf{b}}_c\) denote the row-scaled data matrix and right-hand side. As \( T_s\) increases, the spectral abscissa of the candidate closed-loop matrix moves toward the \( \varepsilon\) -admissible region. The candidates generated for the first five horizons are rejected by the data-driven Lyapunov certificate. At \( T_s=0.3052\) , the data-driven finite-horizon PI converges in five iterations and produces

\[ \mathbf{K}_{c}= \begin{bmatrix} 0.1085 & 0.5936 & 0.1184 & 0.0402 \\ -1.0155 & -0.0247 & -0.6106 & 0.4062\end{bmatrix}. \]

For this gain, the Lyapunov certificate (60) has relative residual \( 7.432\times 10^{-4}\) and \( \lambda_{\min}(\mathbf{P}_c)=5.938\times 10^{-2}>0\) . Therefore, the certificate succeeds and Algorithm 2 terminates. Thus the gain certified from data is an \( \varepsilon\) -admissible initializer for the subsequent infinite-horizon PI. Using \( \mathbf{K}_c\) as the initial policy, the subsequent off-policy RL converges to the optimal solution \( \mathbf{K}^\ast\) within a few iterations. As shown in Fig. 5, both \( \|\mathbf{P}^{[j]}-\mathbf{P}^\ast\|\) and \( \|\mathbf{K}^{[j]}-\mathbf{K}^\ast\|\) decrease rapidly, confirming that the bootstrap gain provides a valid initializer for the infinite-horizon data-driven PI stage.

Pan/Zoom

Pan/Zoom

5.2 Comparative study of the two-mass-spring system

We next consider the two-mass-spring system used in the homotopic PI study [28] to compare Algorithm 2 with existing data-driven RL methods under the same plant and cost. The system matrices are

\[ \mathbf{A}=\begin{bmatrix} 0 & 1 & 0 & 0 \\ -\frac{k_{1}+k_{2}}{m_{1}} & 0 & \frac{k_{2}}{m_{1}} & 0 \\ 0 & 0 & 0 & 1 \\ \frac{k_{2}}{m_{2}} & 0 & -\frac{k_{2}}{m_{2}} & 0\end{bmatrix},~ \mathbf{B}=\begin{bmatrix} 0 \\ \frac{1}{m_{1}} \\ 0 \\ 0\end{bmatrix}. \]

The physical parameters are \( k_{1}=1.5 \mathrm{N/m}\) , \( k_{2}=1 \mathrm{N/m}\) , \( m_{1}=1.1 \mathrm{kg}\) , and \( m_{2}=0.9 \mathrm{kg}\) . The weights are selected as \( \mathbf{Q}= \rm {diag} (1,1,2,1)\) and \( \mathbf{R}=I_{1}\) .

Pan/Zoom

Pan/Zoom

Pan/Zoom

For the existing data-driven homotopic PI method [28], the reported simulation uses \( \mathbf{K}^{[0]}=\mathbf{0}\) and \( \beta -\alpha ^{[0]}=0.1\) . The data are collected over \( [0,5]\) from the initial condition \( x(0)=[0.2,0.2,0.2,0.2]^{T} \) , with a probing input formed by sine and cosine signals of different frequencies. As shown in Fig. 6, the learned stabilizing gain is

\[ \mathbf{K}_{h}^{[101]} = \begin{bmatrix} 14.7218 & 6.5855 & -2.4442 & 16.4721 \end{bmatrix}, \]

which gives the closed-loop poles \( -1.51455\pm 1.66164i\) and \( -1.47887\pm 0.752358i\) .

For comparison, the hybrid PI method [20] uses \( \mathbf{P}^{[0]}=0.1I_{4}\) , \( \epsilon _{j}=1/(j+2)\) , and \( \mathcal{B}_{q}=\{\mathbf{P}>\mathbf{0}:\Vert \mathbf{P}\Vert \leq 5(q+1)\}\) , \( q=0,1,… \) . At each VI iteration, \( \mathbf{H}^{\left[ j\right] }\) and \( \mathbf{K}^{\left[ j\right] }_{v}\) are recovered from the data equation, and \( \widetilde{\mathbf{P}}^{\left[ j+1\right] }\) is generated by the VI update. If \( \widetilde{\mathbf{P}}^{\left[ j+1\right] }\notin \mathcal{B}_{q}\) , the iteration is reset to \( \mathbf{P}^{[0]}\) and \( q\) is increased by one. Otherwise, the VI stage is terminated once the data-driven Lyapunov certificate \( \mathbf{P}^{\left[ j\right] }>\mathbf{0}\) and \( \mathbf{H}^{\left[ j\right] }-2(\mathbf{K}^{[j]}_{v})^{T}\mathbf{R}\mathbf{K}^{\left[ j\right] }_{v}<\mathbf{0}\) is satisfied, where \( \mathbf{H}^{\left[ j\right] }=\mathbf{A}^{T}\mathbf{P}^{\left[ j\right] }+\mathbf{P}^{\left[ j\right] }\mathbf{A}\) . Fig. 7 shows the VI-stage learning process. After one reset caused by the bounded-set test, the spectral abscissa \( \alpha(\mathbf{A}-\mathbf{B}\mathbf{K}^{[j]}_{v})\) remains negative and decreases gradually, as shown in Fig. 7. Meanwhile, \( \lambda_{\max}(\mathbf{H}^{[j]}-2(\mathbf{K}^{[j]}_{v})^{T}\mathbf{R}\mathbf{K}^{\left[ j\right] }_{v})\) crosses zero at the \( 193\) rd VI iteration, which activates the data-driven termination certificate. The certified initializer is

\[ \mathbf{K}_{v}^{[193]} =\begin{bmatrix} 0.7714 & 1.5364 & -0.1604 & 0.9992 \end{bmatrix}, \]

and the corresponding closed-loop poles are \( -0.3376\pm1.6540i\) and \( -0.3607\pm0.7862i\) .

For Algorithm 2, we use \( \varepsilon =0\) , \( \rho =0.1\) , \( T_{0}=0.1\) , \( \gamma =1.25\) , \( \mathbf{M}=\mathbf{0}\) , and \( \mathbf{W}_{c}=I_{4}\) . The finite-horizon policy is initialized by \( \widehat{\mathbf{K}}_{\rho ,T_{s}}^{[0]}(t)=\mathbf{0}\) , and the same basis functions and stopping tolerance \( 10^{-7}\) as in the preceding Section 5.1 are used. As shown in Fig. 8, the data-driven bootstrap procedure succeeds at the first tested horizon \( T_s=0.1 \) . The resulting candidate gain is

\[ \mathbf{K}_{c} = \begin{bmatrix} -0.0057 & 0.0911 & 0.0041 & 0.0003 \end{bmatrix}. \]

For this gain, the data-driven Lyapunov certificate gives \( \lambda_{\min}(\mathbf{P}_c)=10.84354>0\) and a relative residual \( 4.732\times10^{-4}\) , which is below the prescribed tolerance. The closed-loop spectral abscissa is \( \alpha(\mathbf{A}-\mathbf{B}\mathbf{K}_{\mathrm{c}}) =-1.0445\times10^{-2}<0\) , confirming that the certified initializer is admissible. The results in Figs. 6-8 illustrate the distinct advantage of the proposed bootstrap mechanism. The homotopic PI [28] and hybrid PI [20] both succeed in generating admissible initializers, but their certification requires \( 101\) stabilizing iterations and \( 193\) iterations, respectively. By contrast, Algorithm 2 succeeds at the first tested horizon \( T_s=0.1\) and the finite-horizon PI converges in only three iterations.

This article develops a finite-horizon bootstrap method for data-driven infinite-horizon PI without requiring an initial stabilizing policy. By exploiting the well-posedness of finite-horizon policy evaluation and enlarging the horizon for a shifted system, the method generates an admissible candidate gain. An off-policy implementation reconstructs the time-varying value matrix and policy from measured trajectories, while an independent data-driven Lyapunov certificate verifies the candidate gain before it initializes infinite-horizon PI. Numerical studies demonstrate the effectiveness of the proposed method.