Bootstrap Policy Iteration for Stochastic LQ Tracking with Multiplicative Noise

TL;DR

This paper introduces a bootstrap policy iteration algorithm for optimal tracking control of stochastic linear systems with multiplicative noise, enabling model-free reinforcement learning and data-driven computation of control gains.

Contribution

It develops a two-phase bootstrap policy iteration method and a data-driven off-policy RL approach for stochastic LQ tracking, including a novel method for systems with state-dependent noise.

Findings

Converges to the optimal feedback gain under interval excitation.

Ensures stability and optimality in stochastic LQ tracking.

Validated through numerical examples demonstrating effectiveness.

Abstract

This paper studies the optimal tracking control problem for continuous-time stochastic linear systems with multiplicative noise. The solution framework involves solving a stochastic algebraic Riccati equation for the feedback gain and a Sylvester equation for the feedforward gain. To enable model-free optimal tracking, we first develop a two-phase bootstrap policy iteration (B-PI) algorithm, which bootstraps a stabilizing control gain from the trivially initialized zero-value start and proceeds with standard policy iteration. Building on this algorithm, we propose a data-driven, off-policy reinforcement learning approach that ensures convergence to the optimal feedback gain under the interval excitation condition. We further introduce a data-driven method to compute the feedforward using the obtained feedback gain. Additionally, for systems with state-dependent noise, we propose a shadow system-based optimal tracking method to eliminate the need for probing noise. The effectiveness of the proposed methods is demonstrated through numerical examples.

Tables

Table 1. Table 1 : Output trajectories of the desired system with different H d H_{d} .

Case	1	2	3	4	5	6	7	8
$H_d$	$[1,0,0]$	$[2,0,0]$	$[3,0,0]$	$[3,1,-1]$	$[1,1,0]$	$[1,0,1]$	$[0,1,1]$	$[0,1,0]$
$y_d\left(t\right)$

Table 2. Table 2 : Estimated FF gains corresponding to different references.

Case	1	2	3	4	5	6	7	8
${\hat{F}}^{\ast }$	${\left[\begin{array}{c}-29.898\\ 0.189\\ 0.07\end{array}\right]}^{\top }$	${\left[\begin{array}{c}-59.796\\ 0.379\\ 0.14\end{array}\right]}^{\top }$	${\left[\begin{array}{c}-89.695\\ 0.568\\ 0.21\end{array}\right]}^{\top }$	${\left[\begin{array}{c}-89.676\\ -60.932\\ 16.745\end{array}\right]}^{\top }$	${\left[\begin{array}{c}-29.887\\ -23.908\\ -7.42\end{array}\right]}^{\top }$	${\left[\begin{array}{c}-29.905\\ 37.592\\ -23.954\end{array}\right]}^{\top }$	${\left[\begin{array}{c}0.004\\ 13.305\\ -31.514\end{array}\right]}^{\top }$	${\left[\begin{array}{c}0.011\\ -24.098\\ -7.49\end{array}\right]}^{\top }$

Equations

$$\Sigma :\{\begin{array}{rl}& \mathrm{d}x=\left(Ax+Bu\right)\mathrm{d}t+(Cx+Du)\mathrm{d}w,\\ & y=Hx,x(0)=x_0\in {\mathbb{R}}^n,\end{array}$$

$${\Sigma }_d:\{\begin{array}{rl}& {\dot{x}}_d=A_d{x}_d,x_d(0)=x_{d0}\in {\mathbb{R}}^{n_d},\\ & y_d=H_d{x}_d.\end{array}$$

$$P_d(u)=\underset{T\to \infty }{\lim }\frac{1}{T}{J}_T(u),$$

$$J_T=\frac{1}{2}\mathbb{E}\left[{\int }_0^T(\Vert y-y_d{\Vert }_Q^2+\Vert u{\Vert }_R^2)\mathrm{d}t\right],$$

$${\mathcal{L}}_{[K;\mathcal{S}]}(X)=$$

$$+(C-DK{)}^{\top }X(C-DK)$$

$$\sigma ({\mathcal{L}}_{[K;\mathcal{S}]}):=\{\lambda \in \mathbb{C}:{\mathcal{L}}_{[K;\mathcal{S}]}(X)=\lambda X,X\in {\mathbb{S}}^n,X\ne 0\}.$$

$$P_d(u)=P_d(u^{\ast })+\underset{T\to \infty }{\lim }\frac{1}{T}\left(ϵ\delta J_T+\frac{{ϵ}^2}{2}{\delta }^2{J}_T\right),$$

$$\delta J_T=\mathbb{E}\left[{\int }_0^T(f_1^{\top }\delta x+f_2^{\top }\delta u)\mathrm{d}\tau \right]-\mathbb{E}\left[\mathbf{p}(T{)}^{\top }\delta x(T)\right],$$

$$\mathbf{p}(t)=P^{\ast }{x}^{\ast }(t)+{\Pi }^{\ast }{x}_d(t),$$

$$\mathbf{f}(t)=P^{\ast }A{x}^{\ast }(t)+P^{\ast }B{u}^{\ast }(t)+{\Pi }^{\ast }{A}_d{x}_d(t),$$

$$\mathbf{q}(t)=P^{\ast }(C{x}^{\ast }(t)+D{u}^{\ast }(t)).$$

$$\left|\mathbb{E}\left[\mathbf{p}(T{)}^{\top }\delta x(T)\right]\right|\le ({\lambda }_{\max }(P^{\ast })+{ϵ}^2)\mathbb{E}[\Vert x^{\ast }(T){\Vert }^2]$$

$$+{\lambda }_{max}({{\Pi }^{\ast }}^{\top }{\Pi }^{\ast })\mathbb{E}[\Vert x_d(T){\Vert }^2]+{ϵ}^2\mathbb{E}[\Vert x(T){\Vert }^2].$$

$$u^{\ast }(t)=-K^{\ast }{x}^{\ast }(t)-F^{\ast }{x}_d(t),$$

$$K^{\ast }=(R+D^{\top }{P}^{\ast }D{)}^{-1}(B^{\top }{P}^{\ast }+D^{\top }{P}^{\ast }C),$$

$$F^{\ast }=(R+D^{\top }{P}^{\ast }D{)}^{-1}{B}^{\top }{\Pi }^{\ast },$$

$$A^{\top }{P}^{\ast }+P^{\ast }A+Q-C^{\top }{P}^{\ast }C-(P^{\ast }B+C^{\top }{P}^{\ast }D)$$

$$\times (R+D^{\top }{P}^{\ast }D{)}^{-1}(B^{\top }{P}^{\ast }+D^{\top }{P}^{\ast }C)=0,$$

$${\Pi }^{\ast }{A}_d+(A-B{K}^{\ast }{)}^{\top }{\Pi }^{\ast }-H^{\top }Q{H}_d=0,$$

$$\mathcal{Z}(\alpha ):=\{K\in {\mathbb{R}}^{m\times n}|\sigma ({\mathcal{L}}_{[K;S(\alpha )]})\subset {\mathbb{C}}^{-}\}.$$

$${\mathcal{L}}_{[K;\mathcal{S}(\alpha )]}(P)+\Vert K{\Vert }_R^2+\mathcal{Q}=0,$$

$$\mathcal{K}(P):=(R+D^{\top }PD{)}^{-1}(B^{\top }P+D^{\top }PC),$$

$$\mathcal{R}(\alpha ,P,K):=\alpha +\eta \frac{{\lambda }_{\min }\left(\Vert K{\Vert }_R^2+\mathcal{Q}\right)}{{\lambda }_{\max }\left(P\right)},\eta \in (0,1).$$

$$\chi (t)=e^{-\tilde{\alpha }t}x(t),\zeta (t)=e^{-\tilde{\alpha }t}y(t).$$

$$\{\begin{array}{rl}& \mathrm{d}\chi =\left(A({\alpha }_0)\chi +Bu\right)\mathrm{d}t+\left(C\chi +Du\right)\mathrm{d}w,\\ & \zeta =H\chi ,\chi (0)=x_0.\end{array}$$

$$\varphi_{[\alpha,K]}(t)=\left[\begin{array}[]{c}\bar{\delta}_{\chi}(t)+(\alpha-\alpha_{0})\bar{I}_{\chi}(t)\\ -2\big{(}(I_{n}\otimes K)\bar{I}_{\chi}(t)+\bar{I}_{\chi u}(t)\big{)}\\ -\bar{I}_{u}(t)+\bar{I}_{K\chi}(t)\end{array}\right].$$

$$\Phi (\alpha ,K)=[{\phi }_{[\alpha ,K]}(t_1),\cdots ,{\phi }_{[\alpha ,K]}(t_l){]}^{\top },$$

$${\mathcal{I}}_{\chi }=[{\bar{I}}_{\chi }(t_1),\cdots ,{\bar{I}}_{\chi }(t_l)],{\mathcal{I}}_{\chi u}=[{\bar{I}}_{\chi u}(t_1),\cdots ,{\bar{I}}_{\chi u}(t_l)],$$

$${\mathcal{I}}_u=[{\bar{I}}_u(t_1),\cdots ,{\bar{I}}_u(t_l)],{\Delta }_{\chi }=[{\bar{\delta }}_{\chi }(t_1),\cdots ,{\bar{\delta }}_{\chi }(t_l)],$$

$$\mathcal{I}=[{\mathcal{I}}_{\chi }^{\top },{\mathcal{I}}_{\chi u}^{\top },{\mathcal{I}}_u^{\top }{]}^{\top }.$$

$$\mathrm{rank}\left(\mathcal{I}\right)=\frac{1}{2}n(n+1)+nm+\frac{1}{2}m(m+1),\forall l\ge l_1.$$

$${\mathcal{L}}_{[K_{i-1};\mathcal{S}({\alpha }_{i-1})]}(P_i)+\Vert K_{i-1}{\Vert }_R^2+\Theta =0.$$

$$K_i=\mathcal{K}(P_i).$$

$${\alpha }_i=\mathcal{R}({\alpha }_{i-1},P_i,K_i).$$

Full text

{IEEEkeywords}

Optimal tracking control, stochastic system, policy iteration, reinforcement learning

1 Introduction

Reinforcement learning (RL), as a prominent subfield of machine learning, has become increasingly influential for tackling complex optimization and decision-making problems within uncertain environments. RL enables agents to learn optimal strategies through iterative interactions with their environment, guided by feedback in the form of rewards or penalties. This capability has driven extensive research in optimal control field, resulting in notable practical applications and theoretical progress [1, 2]. Two fundamental iterative methods within RL are value iteration (VI) [3, 4, 5] and policy iteration (PI) [6, 7, 8]. VI iteratively updates the value function based on the Bellman optimality principle, ensuring convergence but often slowly. PI alternates between policy evaluation and improvement phases, typically achieving quadratic convergence, but requires an initial stabilizing policy. Enhanced algorithms, such as the $\lambda$-PI algorithm [9, 10], the homotopy-based PI algorithm [11, 12], and hybrid or composite methods [13, 14], have been developed to mitigate these limitations and improve convergence rates. Moreover, efficient exploration is essential in the development of RL-based optimal control to ensure the uniqueness of solutions to data-based iterative equations. This requirement is related to the concept of persistence of excitation (PE) condition [8, 15]. Without the PE condition, RL algorithms fail to converge to the optimal control policy. Recent research has introduced a more relaxed alternative, the interval excitation (IE) condition, implemented through filter-based PI and VI algorithms [16, 17]. Unlike PE, IE only requires sufficient excitation over certain finite intervals. In practice, both PE and IE are usually achieved by injecting probing noise into the input channel during the learning process. Designing this excitation requires carefully balancing the need for exploration with the goal of maintaining system stability and performance. This complex trade-off lies at the core of the exploration-exploitation dilemma (see [18]) and remains an open research topic in RL-based control design.

Among the various optimal control problems, the linear quadratic tracking problem (LQTP) for continuous-time (CT) systems has recently attracted renewed interest as a benchmark for evaluating RL-based control techniques. The objective of LQTP is to guide a linear system along a predetermined trajectory while minimizing a quadratic cost function. Traditional solutions in the RL-based optimal tracking control reformulate LQTPs into equivalent regulation problems by augmenting the state space and redefining the cost function, often in discounted forms. This key transformation, initially proposed in foundational works [19, 20], enables the straightforward application of existing RL methodologies but significantly increases system dimensionality and computational complexity. Subsequent research has developed diverse RL-based optimal tracking methods, including model-free approaches without discount factors [21], output-feedback controllers for incomplete state information [22, 23], off-policy algorithms for external disturbances [24, 25], and exponential tracking controllers using output regulation theory [26, 27].

While most of the existing literature on RL-based optimal control focuses on deterministic systems, stochastic systems, particularly involving multiplicative noise, offer a richer and more realistic modeling framework. This framework is relevant in diverse application areas, including financial systems [28] and networked control systems [29]. However, RL methods tailored for this type of system are scarce, largely due to the analytical complexities introduced by multiplicative noise. These complexities lead to the lack of closed-form solutions and require non-trivial extensions of classical analysis tools like the Popov–Belevitch–Hautus test. Recent foundational advances have addressed these challenges [30, 31, 32, 33, 34], enabling initial developments of on-policy RL algorithms [35, 36] and off-policy RL algorithms [37, 38] for stochastic linear quadratic problems. Further advancements integrating state augmentation techniques with RL methods have also been introduced, providing model-free solutions to the stochastic linear quadratic tracking problem (SLQTP) [39].

Despite the notable theoretical progress brought about by recent advancements in RL for stochastic control, several significant challenges remain. This paper addresses these issues by investigating the CT-SLQTP with multiplicative noise. Firstly, the drawbacks of PI-based methods (i.e., the need for initial stabilizing policies) and VI-based approaches (i.e., slow convergence rates) still exist. But the enhanced methods developed for deterministic systems (e.g. [11, 12, 14]) are difficult to extend directly to stochastic scenarios due to the complexities introduced by multiplicative noise. Secondly, in data-driven methodologies, traditional state augmentation approaches (see [19, 39]) to solving SLQTP significantly increases system dimensionality, making iterative computations more challenging and complicating the fulfillment of necessary IE conditions. Thirdly, a probing noise is typically required in data-driven methods to ensure PE/IE conditions. However, the introduction of probing noise may cause unnecessary oscillations even risk destabilization. Thus, we explore whether full-rank conditions can still be established without probing noise.

To address these challenges, we formulate the SLQTP under an average cost criterion. Utilizing the calculus of variations, we derive an optimal tracking solution that includes a feedback (FB) gain determined by solving a stochastic algebraic Ricatti equation (SARE) matching the original system’s dimension and a feedforward (FF) gain computed through a Sylvester equation. Next, we introduce a designed parameterized system to remove the need for initial stabilizing control gain. Based on this system, we propose a B-PI method that allows convergence from a zero initialization. By leveraging the inherent mean-square stability of the parameterized system, we subsequently develop a novel off-policy RL method incorporating IE conditions, which uses trajectories of the parameterized system. Furthermore, for the scenario where the system volatility depends only on the state, we introduce shadow systems that enable the data-driven matrices to satisfy the IE conditions without the need for probing noise.

Within this proposed SLQTP framework, our main contributions are summarized as follows.

(1) A novel bootstrap PI algorithm is proposed to eliminate the requirement for the initial stabilizing condition. We demonstrate analytically that the stabilizing control gain can be obtained within finite steps, subsequently guaranteeing convergence to the optimal solution.

(2) We further propose an innovative data-driven optimal tracking control algorithm that completely eliminates the need for explicit system dynamics. This method integrates an off-policy RL strategy for FB learning and a data-driven solution for FF computation.

(3) A new shadow system-based optimal tracking method specifically for systems with state-dependent noise is introduced. This methodology ensures necessary full-rank conditions even in the absence of control inputs, thereby avoiding the need for probing noise.

This paper is organized as follows. Section 2 describes the problem formulation and presents basic results for the SLQTP. In section 3, a B-PI framework for solving the SLQTP is introduced, along with convergence analysis. A data-driven optimal tracking design method is then developed, which eliminates the need for prior knowledge of the system dynamics. Finally, to avoid using probing noise while still satisfying the required IE condition, a shadow system-based learning approach is proposed. Section 4 presents simulation examples to demonstrate the effectiveness of the proposed approach. Finally, Section 5 concludes the paper.

The notations used in this paper are introduced as follows. Let $A=[a_{ij}]\in\mathbb{R}^{m\times n}$ and $R=[r_{ij}]\in\mathbb{S}^{m}$ be a symmetric matrix. We define the following operators: $\|A\|_{R}^{2}:=A^{\top}RA$; $\mathrm{vec}(A):=[a_{11},\cdots,a_{m1},a_{12},\cdots,a_{m2},$ $\cdots,a_{1(n-1)},\cdots,$ $a_{m(n-1)},a_{1n},\cdots,a_{mn}]\in\mathbb{R}^{nm}$; $\mathrm{vech}(R):=[r_{11},\cdots,r_{1m},r_{22},\cdots,r_{2m},\cdots,r_{(m-1)(m-1)},$ $r_{(m-1)m},r_{mm}]\in\mathbb{R}^{\frac{1}{2}m(m+1)}$. Let $\{x(t)\in\mathbb{R}^{n},t\geq 0\}$ and $\{y(t)\in\mathbb{R}^{m},t\geq 0\}$ be families of random variables. For a given time interval $T>0$, we define the following terms: $\delta_{x}(t):=\mathrm{vech}(2x(t+T)x(t+T)^{\top}-\mathrm{diag}\left(x(t+T)\right)^{2})-\mathrm{vech}(2x(t)x(t)^{\top}-\mathrm{diag}\left(x(t)\right)^{2})$; $\delta_{xy}(t):=x(t+T)\otimes y(t+T)$ $-x(t)\otimes y(t)$; $I_{x}(t):=\int_{t}^{t+T}\mathrm{vech}(2x(\tau)x(\tau)^{\top}-\mathrm{diag}\left(x(\tau)\right)^{2})\mathrm{d}\tau$; $I_{xy}(t):=\int_{t}^{t+T}x(\tau)\otimes y(\tau)\mathrm{d}\tau$. Let $A^{\dagger}:=(A^{\top}A)^{-1}A^{\top}$ denote the Moore-Penrose inverse of a full column-rank matrix $A$. For any random vector ${x}$, we denote its expectation by $\bar{x}=\mathbb{E}[x]$.

2 Problem formulation and preliminaries

2.1 Problem formulation

Consider a CT stochastic linear system governed by the Itô-type stochastic differential equation (SDE)

[TABLE]

where $x\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$, and $y\in\mathbb{R}^{q}$ represent the system state, control input, and system output, respectively. $w(t)$ denotes a scalar standard Brownian motion defined on the filtered probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ with $\mathcal{F}_{t}=\sigma(w(s):0\leq s\leq t)$. The constant matrices $A,B,C,D,H$ have compatible dimensions and are assumed to be unknown.

We introduce the standard assumption for system (1).

Assumption 1

System (1) or $\mathcal{S}:=[A,B|C,D]$ is mean-square stabilizable and the pair $[A,H|C]$ is exactly detectable.

Remark 1

Mean-square stabilizability guarantees the existence of a FB gain such that the closed-loop system is asymptotically mean square stable (see, e.g., [36, 40]). Exact detectability, a weaker concept than stochastic detectability, was introduced in [34, Definition 3.1], along with a necessary and sufficient condition known as the stochastic Popov–Belevitch–Hautus criterion [34, Theorem 3.1]. These properties are fundamental to designing and analyzing of optimal control strategies for stochastic linear systems [30, 41]. They also play a critical role in the subsequent theoretical developments.

The control objective is to design a control policy such that the system output $y(t)$ tracks a desired reference trajectory $y_{d}(t)$, which is generated by an autonomous linear system

[TABLE]

Here, $A_{d}\in\mathbb{R}^{n_{d}\times n_{d}}$ and $H_{d}\in\mathbb{R}^{q\times n_{d}}$ are also unknown.

Remark 2

The reference system can be designed to be either asymptotically stable or marginally stable. In this work, we consider the more general case of marginal stability, which generates bounded trajectories. This formulation accommodates a broad class of reference signals, including unit steps, sinusoids, and their combinations. Furthermore, varying $H_{d}$ allows for the synthesis of a diverse family of trajectories from the same underlying ODE, as demonstrated in the simulation section.

To formalize this structural property, we introduce the following assumption for the reference generator.

Assumption 2

All eigenvalues of the matrix $A_{d}$ have zero real parts i.e., $\mathrm{Re}(\lambda_{i}(A_{d}))=0$, $i=1,\cdots,n_{d}$, and their algebraic and geometric multiplicities are equal.

To quantify the tracking performance and penalize control effort, we introduce the infinite-horizon average cost functional

[TABLE]

where

[TABLE]

$Q=Q^{\top}>0$, and $R=R^{\top}>0$.

Definition 1** (Admissible control policy)**

A control policy is said to be admissible, denoted by $u\in\mathcal{A}$, if it satisfies the following properties: i). $u(t)=-Kx(t)-Fx_{d}(t):=\mu(x(t),x_{d}(t))$; ii). system (1) with $\mu(x,{0})$ is asymptotically mean-square stable; iii). $\lim_{T\rightarrow\infty}\frac{1}{T}\mathbb{E}\left[\int_{0}^{T}\|x(\tau)\|^{2}\mathrm{d}\tau\right]<\infty$.

For the controlled system (1) and the reference system (2), it is aimed to determine a control policy $u^{*}\in\mathcal{A}$ that minimizes the infinite-horizon average cost functional defined in (4).

2.2 Optimality analysis

We derive the optimal solution of the above SLQTP by applying the calculus of variations directly. First, given matrices from $\mathcal{S}$ and matrix $K\in\mathbb{R}^{m\times n}$, we introduce the generalized Lyapunov operator $\mathcal{L}_{[K;\mathcal{S}]}:\mathbb{S}^{n}\rightarrow\mathbb{S}^{n}$ defined as

[TABLE]

with its spectrum given by

[TABLE]

Let $u^{*}(t)\in\mathcal{A}$ denote the optimal control with associated state trajectory $x^{*}(t)$. We examine perturbed trajectories around this optimal pair $(u^{*},x^{*})$, taking the form ${u}(t)=u^{*}(t)+\epsilon\delta u(t)$ and $x(t)=x^{*}(t)+\epsilon\delta x(t)$, with small $\epsilon\in\mathbb{R}$. Here, $u(t)\in\mathcal{A}$ and $\delta x(t)$ evolves according to $\mathrm{d}\delta x=(A\delta x+B\delta u)\mathrm{d}t+(C\delta x+D\delta u)\mathrm{d}w$ with $\delta x(0)=0$. The cost functional (3) expands as

[TABLE]

where $\delta J_{T}=\mathbb{E}\left[\int_{0}^{T}\!\left((Hx^{*}\!-\!H_{d}x_{d})^{\top}QH\delta x+{u^{*}}^{\top}R\delta u\right)\mathrm{d}\tau\right]$ and $\delta^{2}J_{T}=\mathbb{E}\left[\int_{0}^{T}\!\left(\delta x^{\top}H^{\top}QH\delta x+\delta u^{\top}R\delta u\right)\mathrm{d}\tau\right]$. Given $H^{\top}QH\geq 0$ and $R>0$, it follows immediately that $\lim_{T\rightarrow\infty}\frac{1}{T}\frac{\epsilon^{2}}{2}\delta^{2}J_{T}\geq 0$.

To derive the first-order necessary conditions, introduce a costate process $\boldsymbol{p}:[0,\infty)\rightarrow\mathbb{R}^{n}$ satisfying $d\boldsymbol{p}(t)=\boldsymbol{f}(t)\mathrm{d}t+\boldsymbol{q}(t)\mathrm{d}w(t)$, where $\boldsymbol{f}(t)$ and $\boldsymbol{q}(t)$ will be specified later. Consider the term $\mathbb{E}[\int_{0}^{T}\mathrm{d}(\boldsymbol{p}(t)^{\top}\delta x(t))]$, which equals to $\mathbb{E}[\boldsymbol{p}(T)^{\top}\delta x(T)]$ due to $\delta x(0)=0$. Applying Itô’ s formula, we rewrite this term as $\mathbb{E}\left[\int_{0}^{T}\mathrm{d}(\boldsymbol{p}(t)^{\top}\delta x(t))\right]=\mathbb{E}[\int_{0}^{T}((\boldsymbol{f}+A^{\top}\boldsymbol{p}+C^{\top}\boldsymbol{q})^{\top}\delta x$ $+(B^{\top}\boldsymbol{p}+D^{\top}\boldsymbol{q})^{\top}\delta u)\mathrm{d}t]$. Thus, we have the following equality $\mathbb{E}[\int_{0}^{T}((\boldsymbol{f}+A^{\top}\boldsymbol{p}+C^{\top}\boldsymbol{q})^{\top}\delta x+(B^{\top}\boldsymbol{p}+D^{\top}\boldsymbol{q})^{\top}\delta u)\mathrm{d}t]=\mathbb{E}\left[\boldsymbol{p}(T)^{\top}\delta x(T)\right]$. Substituting this into $\delta J_{T}$ yields

[TABLE]

where f_{1}=(H^{\top}Q(Hx^{*}-H_{d}x_{d})+\boldsymbol{f}+A^{\top}\boldsymbol{p}+C^{\top}\boldsymbol{q}\big{)} and $f_{2}=\left(Ru^{*}+B^{\top}\boldsymbol{p}+D^{\top}\boldsymbol{q}\right)$. Following the sweep method [30, 42], the costate is expressed as

[TABLE]

with $P^{*}={P^{*}}^{\top}\geq 0$ and $\Pi^{*}\in\mathbb{R}^{n\times n_{d}}$, leading to

[TABLE]

Additionally, the terminal term of (6) is bounded by

[TABLE]

Since $u(t),u^{*}(t)\in\mathcal{A}$ and $x_{d}(t)$ is bounded, this gives rise to $\lim_{T\rightarrow\infty}\frac{1}{T}\mathbb{E}\left[\boldsymbol{p}(T)^{\top}\delta x(T)\right]=0$.

Setting the first variation to zero for all admissible perturbations implies the optimality conditions $f_{1}={0}$ and $f_{2}={0}$. By substituting Eqs. (7) and (9) into $f_{2}=0$, we obtain $(R+D^{\top}P^{*}D)u^{*}(t)+(B^{\top}P^{*}+D^{\top}P^{*}C)x^{*}(t)+\Pi^{*}x_{d}(t)=0$. Since the matrix $R+D^{\top}P^{*}D$ is positive definite, we can uniquely solve for the optimal control policy as

[TABLE]

where

[TABLE]

Further substituting Eqs. (7)-(9) and the obtained optimal control policy into the condition $f_{1}=0$ leads to the equality $(\mathcal{L}_{[K^{*};S]}(P^{*})+\|K^{*}\|_{R}^{2}+\|H\|_{Q}^{2})x^{*}+(\Pi^{*}A_{d}+(A-BK^{*})^{\top}\Pi^{*}-H^{\top}QH_{d})x_{d}=0$, which holds for all $[{x^{*}}^{\top},x_{d}^{\top}]^{\top}\in\mathbb{R}^{n+n_{d}}$. Consequently, we deduce the following euquations

[TABLE]

[TABLE]

where Eq. (14) is a SARE and Eq. (15) is a Sylvester equation. Next, we establish sufficient conditions for the optimal solution based on these two equations.

Theorem 1

Suppose Assumptions 1 and 2 hold. Let $P^{*}\geq 0$ and $\Pi^{*}$ be the solutions of (14) and (15), respectively. Then, we have $P_{d}(u)\geq P_{d}(u^{*}),\forall u\in\mathcal{A}$.

Proof: Under Assumption 1, by [34, Theorem 4.1] there exists a unique solution $P^{*}\geq 0$ to Eq (14) such that $\sigma(\mathcal{L}_{[K^{*};\mathcal{S}]})\subset\mathbb{C}^{-}$. It follows that the matrix $A_{c}:=A-BK^{*}$ is Hurwitz. Moreover, Assumption 2 ensures that all the eigenvalues of $A_{d}$ have zero real parts. Therefore, for all $i,j$, $\mathrm{Re}(\lambda_{j}(A_{c}))+\mathrm{Re}(\lambda_{i}(A_{d}))<0$, which implies that Eq. (15) admits a unique solution. As a result, we can conclude that $u^{*}\in\mathcal{A}$. Combing this with Eqs. (6) and (10) and noting the boundedness of $x_{d}$, we obtain $\delta J_{T}(u^{*})=0$. Finally, substituting this result into Eq. (5) yields $P_{d}(u)\geq P_{d}(u^{*})$. $\square$

Remark 3

This approach differs from the augmented formulations found in [19, 20, 39], which reformulate the tracking problem as a regulation problem on an augmented state of dimension $(n+n_{d})$. While those methods lead to an optimal control policy structurally similar to (11), they require solving a Riccati equation of size $(n+n_{d})\times(n+n_{d})$ that includes a redundant block for the construction of $u^{*}$, increasing the computational burden. In contrast, our framework separates the computation: the feedback gain $K^{*}$ is derived from a nonlinear SARE (14), while the feedforward gain $F^{*}$ comes from a linear Eq. (15). The off-policy RL algorithm we propose later is tailored to solve this SARE associated with the FB component.

3 Main results

We propose a B-PI framework to compute the optimal pair $(P^{*},K^{*})$, incorporating both model-based and data-driven algorithms. Unlike conventional PI methods [36, 43], which require an initial stabilizing gain, the B-PI method introduces a novel two-phase iterative mechanism that removes this dependency.

In the first phase, B-PI bootstraps a stabilizing control gain from the zero-value start by iteratively solving a parameterized Lyapunov equation. Once stabilization is achieved, the second phase continues iterating to optimize performance with respect to the specified cost criterion. After obtaining $(P^{*},K^{*})$, these are then used to design a data-based equation for computing $(\Pi^{*},F^{*})$, in place of solving model-based Eq. (15).

3.1 Useful lemmas for parameterized system

To develop these iterative schemes, both model-based and data-driven, we first establish several key results.

Let $\gamma>0$ be a given constant, and consider the parameterized system $\mathcal{S}(\alpha):=(A(\alpha),B|C,D)$ with $A(\alpha):=A-\frac{1}{2}(\gamma-\alpha)I_{n}$. For this system, we define the stabilizing gain set as

[TABLE]

Moreover, given $K\in\mathcal{Z}(\alpha)$, we introduce the parameterized Lyapunov equation

[TABLE]

where $\alpha\in\mathbb{R}$ and $\mathcal{Q}\in\mathbb{S}^{n}$ are design parameters that provide flexibility in shaping the iterative process. Based on the solution $P$, we define the following operators for feedback gain update and the $\alpha$ update as

[TABLE]

We next present an important result, whose proof is provided in Appendix .1. In the remainder of the paper, $\Theta$ denotes a positive definite matrix.

Lemma 1

Let $K\in\mathcal{Z}(\alpha)$ for an arbitrary $\alpha$, and suppose that Assumption 1 holds. Let ${P}$ be the unique solution of Eq. (16), and let $K^{\prime}=\mathcal{K}(P)$ and $\alpha^{\prime}=\mathcal{R}(\alpha,P,K^{\prime})$. Then the following statements hold:

If $\mathcal{Q}=\Theta$ and $\alpha\in\mathcal{(}0,\gamma)$, then $0<{P}<\infty$, $K^{\prime}\in\mathcal{Z}(\alpha^{\prime})$, and $\alpha^{\prime}>\alpha$. 2. 2.

If $\mathcal{Q}=\|H\|_{Q}^{2}$ and $\alpha=\gamma$, then ${P}\geq P^{*}$ and $K^{\prime}\in\mathcal{Z}(\gamma)$, where $P^{*}$ is the solution of Eq. (14).

Remark 4

Lemma 1 forms the analytical cornerstone for the convergence result of the B-PI algorithm along with the closed-loop asymptotic mean-square stability of the parameterized system. Specifically, when $\mathcal{Q}=\Theta$, Statement 1 of Lemma 1 is crucial for ensuring that the search for a stabilizing control gain can be completed within a finite number of steps. Furthermore, note that when $\alpha=\gamma$, the parameterized system $\mathcal{S}(\gamma)$ coincides with the original system $\mathcal{S}$. In this setting, choosing $\mathcal{Q}=\|H\|_{Q}^{2}$, corresponding to the cost functional (3), allows the iterative procedure governed by (16)-(17) to approximate the optimal solution. Statement 2 of Lemma 1 is thus key for establishing convergence.

In the subsequent two-phase scheme, the parameter $\mathcal{Q}$ will be assigned different symmetric matrices depending on the respective phase. Meanwhile, the parameter $\alpha$ will be iteratively updated during Phase I according to (18), starting from the following initial condition.

Lemma 2

If there exist $\gamma>\alpha_{0}>0$ such that $\gamma>\sigma(\mathcal{L}_{[0;\mathcal{S}]})+\alpha_{0}$, then $0\in\mathcal{Z}(\alpha_{0})$.

Proof: Given $\gamma>\sigma(\mathcal{L}_{[0;\mathcal{S}]})+\alpha_{0}$, we obtain $\sigma(\mathcal{L}_{[0;\mathcal{S}(\alpha_{0})]})=\sigma(\mathcal{L}_{[0;\mathcal{S}]})-(\gamma-\alpha_{0})\subset\mathbb{C}^{-}$. Thus, it follows directly that $0\in\mathcal{Z}(\alpha_{0})$. $\square$

Remark 5

Note that the scalar $\gamma$ can be chosen sufficiently large and $\alpha_{0}$ sufficiently small to meet the required initial condition. Consequently, the parameterized system $\mathcal{S}(\alpha_{0})$ is inherently asymptotically mean-square stable. Compared to the restrictive assumption of the existence of a stabilizing gain matrix, this condition is easy to achieve. It plays an important role not only in the development of the model-based B-PI approach but also in the establishment of the data-driven method. In particular, within the data-driven framework, this inherent stability allows for designing the input $u(t)$ solely based on probing noise, which is typically time-varying and sufficiently rich [8, 18]. Such a design ensures the collected data remains informative for learning, without requiring an explicit feedback control component.

To facilitate this, we define $\tilde{\alpha}:=\frac{1}{2}(\gamma-\alpha_{0})$, and introduce the discounted state and output variables

[TABLE]

Then the system dynamics given in (1) is transformed into the parameterized system $\mathcal{S}(\alpha_{0})$:

[TABLE]

Thus, the trajectories of $\mathcal{S}(\alpha_{0})$ can be interpreted as a discounted transformation (DT) of those of the original system $\mathcal{S}$, with the input $u(t)$ unchanged.

Given parameters $\alpha$ and $K\in\mathbb{R}^{m\times n}$, we can get the expected integrals $\bar{I}_{\chi}(t),\bar{I}_{\chi u}(t),\bar{I}_{u}(t),\bar{I}_{K\chi}(t)$, and $\bar{\delta}_{\chi}(t)$, along the trajectories (i.e., $\chi(t),u(t)$) over the interval $[t,t+T]$, where $\bar{(\cdot)}:=\mathbb{E}[(\cdot)]$, and construct

[TABLE]

Let $T_{s}>0$ be a sampling interval, and define the discrete time instants $t_{i}:=t_{1}+(i-1)T_{s}$, where $t_{1}$ is the initial time for data collection. Denote the matrices

[TABLE]

We have the following result, whose proof is provided in Appendix .2.

Assumption 3

There exists an integer $l_{1}>0$ such that

[TABLE]

Lemma 3

Suppose Assumption 3 holds. If $K\in\mathcal{Z}(\alpha)$, then the matrix $\varPhi(\alpha,K)$ has full column rank.

Remark 6

Lemma 3 provides theoretical support for the convergence of the data-driven method. The rank condition in Assumption 3 corresponds to the IE condition. Specifically, the full rank of $\mathcal{I}$ implies positive definiteness of $T_{s}\mathcal{I}\mathcal{I}^{\top}$, which is a Riemann sum approximation to the integral $\int_{t_{1}}^{t_{l}+T}\bar{I}(\tau)\bar{I}(\tau)^{\top}\mathrm{d}\tau>0$, where $\bar{I}=[\bar{I}_{\chi}^{\top},\bar{I}_{\chi u}^{\top},\bar{I}_{u}^{\top}]^{\top}$. Unlike the strict PE requirement, which demands sufficient energy in every sliding window, IE only asks that the signal be rich enough over the finite interval. Moreover, it is worth emphasizing that the integration interval $T$ and the sampling period $T_{s}$ are fundamentally different: $T$ specifies the time span over which data is integrated, while $T_{s}$ represents the sampling frequency and can be chosen arbitrarily small depending on the capabilities of the hardware.

3.2 Bootstrap policy iteration algorithm

With the theoretical foundation established by Lemma 1, we are now ready to formalize the B-PI algorithm, summarized in Algorithm 1. The algorithm operates in two phases: Phase I is outlined in lines 2-7, while Phase II is given in lines 8-12. Define the index $\mathbb{I}:=\min\{i\in\mathbb{N}|\alpha_{i}\geq\gamma\}$. Theoretical results of this algorithm are established in Theorem 2, with detailed proof given in Appendix .3.

Throughout the subsequent analyses, we consistently initialize $K_{0}$ (or $\hat{K}_{0}$) as the zero matrix and denote by $P^{*}\geq 0$ and $K^{*}$ the unique solutions to Eqs. (14) and (12), respectively.

Theorem 2

Suppose Assumption 1 holds, and there exist constants $\gamma>\alpha_{0}>0$ such that $\gamma>\sigma(\mathcal{L}_{[0;\mathcal{S}]})+\alpha_{0}$. Let the sequence $\{(P_{i},K_{i},\alpha_{i})\}_{i=1}^{\mathbb{I}}$ be generated by iteratively solving (23)-(25), and let the sequence $\{(P_{i},K_{i})\}_{i=\mathbb{I}+1}^{\infty}$ be generated by iteratively solving (24) and (26). Then the following properties are satisfied:

$\mathbb{I}<\infty$* and $K_{\mathbb{I}}\in\mathcal{Z}(\gamma)$. *

• 2)

$\lim_{i\rightarrow\infty}P_{i}=P^{*}$* and $\lim_{i\rightarrow\infty}K_{i}=K^{*}$.*

Remark 7

Classical PI methods (e.g. [35, 36, 44]) for CT stochastic linear systems require a stabilizing initial feedback gain. In contrast, our B-PI framework removes this restrictive requirement by introducing a parameterized system $\mathcal{S}(\alpha)$ and a parameterized Lyapunov equation (16). Phase I automatically drives $\alpha_{i}$ upward and returns a stabilizing gain in finite number of iterations even from $K_{0}=0$; Phase II continues on the original system to reach the optimal solution. Unlike [14], which enforces $Q>0$ to obtain a monotonic sequence $\{\alpha_{i}\}$, our design employs a tunable weighting matrix in (23) that decouples the “stabilize‑first” stage from the task cost, thereby relaxing the requirement on $Q$. As a result, the method accommodates $Q\geq 0$ and extends to indefinite $Q$ when combined with the result in [45].

Note that in practical implementations, the calculated control gains may be biased due to model uncertainties. To address this issue, we can directly extend the B-PI framework to explicitly account for these disturbances, and analyze the robustness of the algorithm under perturbed conditions by using [45, Theorem 6]. It demonstrates that if the perturbations remain small, the algorithm ensures convergence of the sequence to a neighborhood around the optimal solution.

3.3 Data-driven optimal tracking control

We now develop a data-driven method to sequentially determine the gain matrices $K^{*}$ and $F^{*}$, without requiring knowledge of the system parameters, based on Algorithm 1 and Eq. (15).

To derive a data-based iterative equation, we first recall system (20) and apply Itô’s formula to $\chi^{\top}P_{i}\chi$ to obtain

[TABLE]

Let $\mu_{i}:=K_{i}\chi$, $M_{i}:=\Upsilon_{i}K_{i}$, and $\Lambda_{i}:=D^{\top}P_{i}D$. We then substitute (24) into the above equation, and use $\mathscr{L}_{[K_{i-1};\mathcal{S}(\alpha)]}(P_{i})$ to rearrange the expression, which gives

[TABLE]

where

[TABLE]

By substituting $\alpha_{i-1}$ for $\alpha$ in (27) and applying the iterative Eq. (23) from Phase I of Algorithm 1 to the right-hand side of Eq. (27), we obtain

[TABLE]

Similarly, replacing $\alpha$ with $\gamma$ and applying Eq. (26) from Phase II of Algorithm 1, the expression becomes

[TABLE]

Next, we integrate both sides of Eqs. (28) and (29) along (20) over $[t,t+T]$, take expectations, and use Kronecker product representation, to obtain

[TABLE]

where $\theta_{i}=[\mathrm{vech}(P_{i})^{\top},\mathrm{vec}(M_{i})^{\top},\mathrm{vech}(\Lambda_{i})^{\top}]^{\top}$, $\psi_{i}(t)=\varphi_{[\alpha_{i},K_{i}]}(t)$, $\phi_{i}(t)=\varphi_{[\gamma,K_{i}]}(t)$, and $\bar{I}_{\zeta}(t)=\mathbb{E}[I_{\zeta}(t)]$.

Based on the sampling points, we construct the matrices $\Psi_{i}:=[\psi_{i}(t_{1}),\cdots,\psi_{i}(t_{l})]^{\top}$, $\Phi_{i-1}:=[\phi_{i}(t_{1}),\cdots,\phi_{i}(t_{l})]^{\top}$, and $\mathcal{I}_{\zeta}:=[\bar{I}_{\zeta}(t_{1}),\cdots,\bar{I}_{\zeta}(t_{l})]$. This leads to the matrix-form representations of Eqs. (30) and (31) as

[TABLE]

Since $\Psi_{i}=\varPhi(\alpha_{i},K_{i})$ and $\Phi_{i}=\varPhi(\gamma,K_{i})$, and based on the data-driven matrices and theoretical foundations established in Lemma 3), we are now in the position to develop an off-policy RL algorithm as an extension of Algorithm 1. In particular, lines 4-5 of the model-based algorithm are replaced by lines 5-6 of Algorithm 2.

In addition, lines 10-11 of Algorithm 1 are modified as lines 11-12 of Algorithm 2.

We have the following results, whose proof is provided in Appendix .4.

Theorem 3

Suppose the conditions of Theorem 2 and Assumption 3 are satisfied. Let the sequence $\{(P_{i},K_{i},\alpha_{i})\}_{i=1}^{\mathbb{I}}$ be generated by iteratively solving Eqs. (39),(40), and (25), and let the sequence $\{(P_{i},K_{i})\}_{i=\mathbb{I}+1}^{\infty}$ be generated by iteratively solving Eqs. (42) and (40). Then the following properties hold:

$\mathbb{I}<\infty$* and $K_{\mathbb{I}}\in\mathcal{Z}(\gamma)$.*

• 2)

$\lim_{i\rightarrow\infty}P_{i}=P^{*}$* and $\lim_{i\rightarrow\infty}K_{i}=K^{*}$.*

Remark 8

From Eq. (30) the unknowns $\theta_{i}$ also satisfies the integral closed form

[TABLE]

The update (39) is precisely the discrete normal equation obtained by stacking samples of $\psi_{i}$ and ${\bar{I}}_{\chi}$. It can therefore be viewed as a Riemann‑sum or least‑squares approximation of (34). Importantly, Eq. (39) is derived directly from the expected identity (30), rather than by numerically integrating trajectories to approximate (34). This avoids bias, which can be non‑negligible when the integral window $T$ is short or sampling is coarse.

After obtaining the optimal FB gain $K^{*}$ and matrix $\Lambda^{*}:=D^{\top}P^{*}D$, we proceed to develop a data-based equation for computing the optimal FF gain $F^{*}$. To this end, define the expected values $\bar{\chi}=\mathbb{E}(\chi),\bar{u}=\mathbb{E}(u),\bar{\zeta}=\mathbb{E}(\zeta)$. The dynamics of $\bar{\chi}$ satisfy $\dot{\bar{\chi}}=A(\alpha_{0})\bar{\chi}+B\bar{u}$.

Using Eqs. (13) and (15), the time derivative of the term $\bar{\chi}(t)^{\top}\Pi^{*}x_{d}(t)$ can be expressed as

[TABLE]

By integrating Eq. (35) over $[t,t+T]$ along the dynamics of $x_{d}$ and $\bar{\chi}$, and rearranging terms, we obtain

[TABLE]

where $\vartheta^{*}=[\mathrm{vec}(\Pi^{*})^{\top},\mathrm{vec}(F^{*})^{\top}]^{\top}$, and $\xi(t)$ is given by

[TABLE]

To facilitate a matrix-based formulation, we define the data-based matrices, ${\Delta}_{x_{d}\bar{\chi}}\!:=\![{\delta}_{x_{d}\bar{\chi}}(t_{1}),\!\cdots\!,\!{\delta}_{x_{d}\bar{\chi}}(t_{l})]$, $\mathcal{I}_{x_{d}\bar{\chi}}\!:=\![I_{x_{d}\bar{\chi}}(t_{1}),\!\cdots\!,\!I_{x_{d}\bar{\chi}}(t_{l})]$, $\mathcal{I}_{x_{d}\bar{u}}\!:=\![I_{x_{d}\bar{u}}(t_{1}),\!\cdots\!,\!I_{x_{d}\bar{u}}(t_{l})]$, and $\mathcal{I}_{y_{d}\bar{\zeta}}\!:=\![I_{y_{d}\bar{\zeta}}(t_{1}),\!\cdots\!,\!I_{y_{d}\bar{\zeta}}(t_{l})]$. We then construct the matrix $\Xi\!:=\![\xi(t_{1}),\cdots\!,\xi(t_{l})]^{\top}$, which has the following structure

[TABLE]

Therefore, the linear matrix form of Eq. (36) becomes

[TABLE]

Lemma 4

If there exists $l_{1}>0$ such that for all $l\geq l_{1}$ the following rank condition holds

[TABLE]

then $\Xi$ has full column rank.

The proof is similar to that of Lemma 3 and is omitted for brevity. As a result, under this condition, the unknown parameters can be uniquely determined by solving 43.

Remark 9

As summarized in Algorithm 2, the method consists of two main components: the FB off-policy RL (lines 3–13) and the FF computation (line 15). These components are executed sequentially while utilizing the same dataset. In contrast to most existing methods that reformulate tracking problems into regulation problems of dimension $n+n_{d}$ [19, 20, 39], our approach strategically decomposes the FB and FF control design, significantly reducing computational burden. Moreover, the convergence analysis of the iterative procedure relies on the asymptotic mean-square stability property of the parameterized system’ s closed-loop dynamics, rather than those of the original system.

3.4 Shadow system-based RL (case $D=0$)

In this subsection, we consider the special case where the matrix $D$ vanishes ($D=0$). Under this condition, the term $\Lambda_{i}$ consistently vanishes, simplifying the iterative procedure by removing its calculation. Simultaneously, the term $M_{i}$ reduces to $RK_{i}$.

As previously noted in Remark 9, although the transformed dynamics are asymptotically mean-square stable, the original system may still be inherently unstable. Therefore, injecting probing noise could amplify unbounded trajectories, posing a practical risk in real-world control applications. Conversely, as emphasized in Remark 5, the input $u(t)$ is specifically designed to include only probing noise due to this mean-square stability. In this regard, an important question arises: What would happen if no probing noise were injected into the input channel, i.e., if $u(t)=\ell(t)=0$?

To analyze this scenario, we rewrite Eqs. (39), (42), and (43) under conditions $D=0$ and $u=\bar{u}=0$, yielding

[TABLE]

Here, $\theta_{i}$ becomes $[\mathrm{vech}({P}_{i})^{\top},\mathrm{vec}({K}_{i})^{\top}]^{\top}$, and the columns of $\Psi_{i-1}$, $\Phi_{i-1}$, and $\Xi$ are given by

[TABLE]

Since $\bar{I}_{\chi u}=0$ and ${I}_{x_{d}\bar{u}}=0$, the previously stated full-rank conditions (Lemmas 3 and 4) no longer hold. Consequently, Algorithm 2 fails in the absence of probing noise.

However, certain partially model-free, on-policy methods ([19, 46]) circumvent this limitation by resetting initial states. This observation suggests a possible solution: would knowing matrix $B$ solve this issue?

Inspired by this insight, we propose a new data-driven algorithm employing shadow systems to maintain an off-policy framework. Specifically, given a known $B$, we define a shadow system, $\Sigma_{a}$, for FB learning as

[TABLE]

and a shadow system, $\Sigma_{da}$, for FF computing as

[TABLE]

From these systems, we form matrices $\Omega_{K}:=[{\omega}_{K}(t_{1}),\cdots$ $,{\omega}_{K}(t_{l})]^{\top}$ and $\Omega_{F}:=[{\omega}_{F}(t_{1}),\cdots,{\omega}_{F}(t_{l})]^{\top}$, with components given by

[TABLE]

$I_{x_{a}^{\prime}}:=\int_{t}^{t+T}\mathrm{vech}(x_{a}^{\prime})\mathrm{d}\tau$, $x_{a}^{\prime}:=2(A_{a}x_{a}x_{a}^{\top}+x_{a}x_{a}^{\top}A_{a}^{\top})-\mathrm{diag}(A_{a}x_{a}x_{a}^{\top}+x_{a}x_{a}^{\top}A_{a}^{\top})$. These matrices are embedded into the original data-based matrices:

[TABLE]

The next result confirms the effectiveness of the shadow system approach. The proof follows the same structure as that of Lemma 3 and is omitted here for brevity. Our focus is on discussing the advantages of introducing shadow systems.

Corollary 1

Let $D=0$ and $u(t)=0$. If there exists $l_{1}>0$ such that for all $l\geq l_{1}$,

[TABLE]

then the matrices $\Psi_{i}^{a}$, $\Phi_{i}^{a}$, and $\Xi^{a}$ have full column ranks.

As a result, unknown parameters can be uniquely determined by

[TABLE]

Combining the uniqueness result with the derivations of Eqs. (52), (53), and (54), we can conclude that the shadow system-based solutions are equivalent to their model-based counterparts. Therefore, using expectation approximations, the estimated version of shadow system-based equations can be substituted into Algorithm 2, which results in a shadow system-based optimal tracking algorithm. The structure of this algorithm is illustrated in Fig. 1.

Remark 10

The key advantages of the proposed approach are summarized below:

• 1).

The matrices $\mathcal{I}_{x_{a}u_{a}}$ and $\mathcal{I}_{y_{a}u_{a}}$ depend on the shadow systems, offering flexibility in generating excitation signals. The matrix $A_{a}$ can be chosen to make the pair $(A_{a},B)$ controllable. This allows $u_{a}$ to be designed to excite the full system modes and achieve the required rank conditions. Additionally, the rank conditions of $\mathcal{I}_{\chi}$ and $\mathcal{I}_{x_{d}\bar{\chi}}$ can be readily achieved by resetting the system’s initial states.

• 2).

In contrast to partially model-free methods (e.g. [19, 46]), which require applying an updated control policy to the physical system to collect data at each iteration, the proposed method is off-policy. It only requires a single data collection phase prior to learning. Furthermore, the method enables the simultaneous update of both the value function and the control gain.

• 3).

During learning, neither control input/probing noise nor shadow trajectories are injected into the physical system. Shadow trajectories enter only algebraically in the iterative equations to satisfy IE, avoiding any amplification of potentially unstable on‑plant behavior.

4 Numerical Simulations

This section presents two numerical examples to validate the effectiveness, accuracy, and robustness of the proposed data-driven algorithms for SLQTP. In both cases, the system dynamics incorporate stochastic diffusion, modeled as multiplicative noise, while the drift components are represented by two well-known benchmark systems: a spring-mass-damper system (see [20]) and a two-mass-spring system (see [27]). The simulation results constitute three parts:

• 1).

Demonstration of the algorithm’s capability to learn a stabilizing FB gain directly from data, without requiring explicit system dynamics. Furthermore, the convergence and robustness properties are verified.

• 2).

Evaluation of the flexibility and accuracy of the FF computation by applying the same learned FB gain to track multiple desired trajectories. Moreover, comparative analysis against the LQTP model-based algorithm ($C=0,D=0$) highlights the importance of considering multiplicative noise to ensure accurate tracking.

• 3).

Examination of the shadow system-based learning method. Results confirm that the rank conditions necessary for FB learning and FF computation are satisfied by generating appropriate trajectories using shadow systems, without injecting probing noise into the original input channel.

4.1 Numerical example $I$

The spring-mass-damper system with multiplicative noise are characterized by the following parameters

[TABLE]

The quadratic cost function employs weight matrices $Q=10$ and $R=0.01$. A collection of desired output trajectories is generated by a time-invariant linear system

[TABLE]

Eight different output maps $H_{d}$ for each reference case are detailed in Table 1.

Algorithm 2 parameters include $\gamma=1$, $\alpha_{0}=0.1$, $\eta=0.95$, $\Theta=\mathrm{diag}([10,10])$, and $\epsilon=10^{-5}$. The control input includes probing noise given by [8, 37] as $u(t)=10\sum_{j=1}^{50}\sin(\omega_{j}t)$, where $\omega_{j}\in[-100,100]$ are randomly selected. Data is collected over the interval $[t_{1},t_{l}]=[0,5+T]$, with sampling period $T_{s}=0.001$ and integration window $T=0.1$. The rank conditions (22) and (38) are satisfied at the end of this phase. In order to distinguish the results of data-driven algorithms from those of model-based algorithms, the results of data-driven algorithms are represented by $\hat{(\cdot)}$ in the following text.

The iterative procedure starts with $K_{0}=0$. During Phase I, the iterations continue until the parameter $\hat{\alpha}_{i}$ reaches the threshold $\gamma=1$. As shown in Fig. 2. (c), this condition is first met at iteration $i=2$, resulting in a stabilizing FB gain matrix $\hat{K}_{2}=\left[24.8937,42.5251\right]$ such that $\sigma(\mathcal{L}_{[\hat{K}_{2};\mathcal{S}]})\subset\mathbb{C}^{-}$. The iterative process converges at iteration $i=9$, with the final gain $\hat{K}_{9}=\left[25.8381,7.5651\right]$. The Evolution of $\{\hat{P}_{i}\}$ and $\{\hat{K}_{i}\}$ are illustrated in Fig. 2 (a)-(b). Accuracy of the learned gains is assessed by comparing the data-driven results with Algorithm 1 (used as ground truth). The error norms at each iteration are provided in Fig. 2 (d). The final errors are found to be small, while minor discrepancies across iterations may result from numerical approximations and data noise. It can be concluded that as long as these disturbances remain small, the parameters converge to a small neighborhood around the optimal solution.

Using the converged FB gain $\hat{K}_{9}$, FF gains are computed for all eight reference trajectories via equation (43). The results are summarized in Table 2. All FF gains are computed based on the same dataset from $\mathcal{S}(\alpha_{0})$. To assess tracking performance, two reference tracking scenarios are constructed:

• •

Scenario 1: $y_{d}^{1}\to y_{d}^{2}\to y_{d}^{3}\to y_{d}^{4}$, with the last segment lasting $10$ seconds and the others $5$ seconds each.

• •

Scenario 2: $y_{d}^{1}\to y_{d}^{5}\to y_{d}^{6}\to y_{d}^{7}\to y_{d}^{8}$, each lasting 5 seconds.

Here, the superscript represents case index. Fig. 3 shows that in both scenarios, the system achieves fast and accurate tracking over a 25-second horizon.

Additionally, to benchmark the SLQTP tracker’s effectiveness, we compare it to the LQTP tracker (assuming no stochastic diffusion i.e., $C=0$ and $D=0$). In the tracking task for case 8, Fig. 4 (a) shows the SLQTP tracking, while Fig. 4 (b) shows the LQTP result. The results indicate that while both trackers stabilize the system, the SLQTP tracker significantly outperforms the LQTP tracker, maintaining output trajectories close to the desired reference. Thus, it is important to consider multiplicative noise to ensure accurate tracking.

4.2 Numerical example $II$

In this example, we evaluate the performance of the shadow system-based learning method in no probing noise setting. The controlled system under consideration is given by

[TABLE]

The weight matrices of the quadratic cost function are $Q=100$ and $R=1$, and the desired trajectories are the same as those in Example $I$. Set $u(t)=0$.

To satisfy the rank conditions (50) and (51) in the absence of external excitation, auxiliary systems are introduced. Initial states are reset once and state trajectories are collected over the interval $[0,10]$, with each trajectory spanning $[0,5]$. The auxiliary systems are designed as

[TABLE]

The matrix $A_{a}$ is Hurwitz and controllable with $B$, and $F_{a}$ is marginally stable. All other algorithm parameters remain the same as in Example $I$.

Starting from $\alpha_{0}=0.1$ and $K_{0}=0$, the parameter $\hat{\alpha}_{i}$ meet the threshold $\gamma=1$ at iteration $i=4$, resulting in a stabilizing gain $\hat{K}_{4}=[5.2039,4.4951,-1.8831,3.3963]$. The algorithm proceeds and converges by iteration $12$, as shown in Fig. 5 (a)-(b). To evaluate the accuracy of the learning process, we compare the parameters$\{\hat{P}_{i}\}$ and $\{\hat{K}_{i}\}$ against model-based results. The error norms, provided in Fig. 5 (c), remain low throughout, demonstrating the algorithm’s robustness. Furthermore, tracking performance validation is conducted using Scenario 2 from Example I. The tracking results, shown in Fig. 5 (d), indicate that the learned control policy achieves accurate and rapid tracking performance.

5 Conclusion

This paper introduced a B-PI-based optimal tracking framework for the SLQTP in the presence of multiplicative noise. The proposed data-driven method removes the traditional requirement for an initial stabilizing control gain and explicit knowledge of the dynamical model. Furthermore, a shadow system-based optimal tracking algorithm was developed to eliminate the requirement for probing noise by assuming that volatility depends only on the system state and that the input-to-state matrix is given. The off-policy RL algorithm in this methodology can be directly applied to the stochastic linear quadratic regulation problem.

.1 Proof of Lemma 1

1). Given $K\in\mathcal{Z}(\alpha)$, we first consider Eq. (16) with $\mathcal{Q}=\Theta$. According to [47, Theorem 3.2.3], this equation admits a unique solution $P$ given by

[TABLE]

where $\{X(t);t\geq 0\}$ satisfies the matrix SDE

[TABLE]

Since $K\in\mathcal{Z}(\alpha)$ also implies exponential mean-square stability of the above system due to time invariance (see [48]), and since $R,\Theta>0$, it follows from Eq. (59) that

[TABLE]

To establish that $\alpha^{\prime}>\alpha$ and $K^{\prime}\in\mathcal{Z}(\alpha^{\prime})$, we substitute the definition of $K^{\prime}$ from (17) into Eq. (16) and complete the squares to obtain

[TABLE]

where $\Theta^{\prime}=\|K^{\prime}\|_{R}^{2}+\Theta>0$ and $\Upsilon=R+D^{\top}PD>0$. Rewriting $\alpha^{\prime}=\mathcal{R}(\alpha,P,K^{\prime})$ with $\mathcal{Q}=\Theta$ yields

[TABLE]

Since $\eta\in(0,1)$, $\Theta^{\prime}>0$, and $P>0$, it clearly implies that $\alpha^{\prime}-\alpha>0$. Moreover, we get

[TABLE]

Combing (60) and (61) gives

[TABLE]

By [31, Theorem 1], we conclude that $\sigma(\mathcal{L}_{[K^{\prime};\mathcal{S}(\alpha^{\prime})]})\subset\mathbb{C}^{-}$, and thus $K^{\prime}\in\mathcal{Z}(\alpha^{\prime})$.

2). Since $\mathcal{S}(\gamma)=\mathcal{S}$, we subtract Eq. (14) from Eq. (16) with $\mathcal{Q}=\|H\|_{Q}^{2}$ and $\alpha=\gamma$ to obtain

[TABLE]

where $\Upsilon^{*}=R+D^{\top}P^{*}D>0$. Together with $K\in\mathcal{Z}(\gamma)$, this equation implies that ${P}\geq{P}^{*}$.

Next, we aim to show that $K^{\prime}\in\mathcal{Z}(\gamma)$. Rewriting Eq. (62) using $K^{\prime}=\mathcal{K}(P)$ gives

[TABLE]

where $\Delta K:=\|K^{\prime}-{K}^{*}\|_{\Upsilon^{*}}^{2}+\|K^{\prime}-{K}\|_{\Upsilon}$. According to [34, Theorem 3.2], since $P-P^{*}\geq 0$ and $\Delta K\geq 0$ hold, it is sufficient to prove the exact detectability of the pair $(A-BK^{\prime},\Delta K|C-DK^{\prime})$ to conclude that $K^{\prime}\in\mathcal{Z}(\gamma)$.

Suppose, for contradiction, that the pair is not exactly detectable. Then, by the stochastic Popov–Belevitch-Hautus criterion ([34, Theorem 3.1]), there exists a nonzero symmetric matrix $X\neq 0$ such that

[TABLE]

From $\Delta KX=0$, it follows $K^{\prime}X={K}X$, which leads to $\mathcal{L}_{[K^{\prime};\mathcal{S}]}(X)=\mathcal{L}_{[{K};\mathcal{S}]}(X)=\lambda X$. But since $K\in\mathcal{Z}(\gamma)$, we have $\sigma(\mathcal{L}_{[K;\mathcal{S}]})\subset\mathbb{C}^{-}$, which contradicts $\mathrm{Re}(\lambda)\geq 0$. Thus, the exact detectability of $(A-BK^{\prime},\Delta K|C-DK^{\prime})$ holds. This completes the proof. $\square$

.2 Proof of Lemma 3

We prove the result by contradiction. Suppose there exists a nonzero vector $\mathcal{Y}=[\mathcal{Y}_{1}^{\top},\mathcal{Y}_{2}^{\top},\mathcal{Y}_{3}^{\top}]^{\top}\neq 0$, such that $\varPhi(\alpha,K)^{\top}\mathcal{Y}={0}$. Let $\mathcal{Y}_{1}=\mathrm{vech}(Y_{1})\in\mathbb{R}^{\frac{n}{2}(n+1)}$, $\mathcal{Y}_{2}=\mathrm{vec}(Y_{2})\in\mathbb{R}^{mn}$ with $Y_{2}\in\mathbb{R}^{m\times n}$, and $\mathcal{Y}_{3}=\mathrm{vech}(Y_{3})\in\mathbb{R}^{\frac{m}{2}(m+1)}$.

Using $\varPhi(\alpha,K)^{\top}\mathcal{Y}=0$, the definition of $\varPhi(\alpha,K)$, and Eq. (21), we obtain

[TABLE]

Next, applying Itô’s formula to $\chi^{\top}Y_{1}\chi$ along with (20), substituting $\mathcal{L}_{[K;\mathcal{S}(\alpha)]}$, and rearranging the terms yields

[TABLE]

where $F_{1}(Y_{1})=\mathcal{L}_{[K;\mathcal{S}(\alpha)]}(Y_{1})-(\alpha-\alpha_{0})Y_{1}+K^{\top}(B^{\top}Y_{1}+D^{\top}Y_{1}C)+(Y_{1}B+C^{\top}Y_{1}D)K_{\alpha}-K_{\alpha}^{\!\top}D^{\!\top}Y_{1}DK_{\alpha}$, $F_{2}(Y_{1})=B^{\top}Y_{1}+D^{\top}Y_{1}C$, and $F_{3}(Y_{1})=D^{\top}Y_{1}D.$

We then integrate over multiple time intervals $[t_{1},t_{1}+T],\cdots,[t_{l},t_{l}+T]$, take expectations, use the Kronecker product representation, and construct the matrix form of the obtained linear equations. This gives

[TABLE]

Adding Eqs. (64) and (65) results in

[TABLE]

where $M_{1}=\mathcal{L}_{[K;\mathcal{S}(\alpha)]}(Y_{1})+K^{\top}M_{2}+M_{2}^{\top}K-K^{\top}M_{3}K$, $M_{2}=B^{\top}Y_{1}+D^{\top}Y_{1}C-Y_{2}$, $M_{3}=D^{\top}Y_{1}D-Y_{3}$.

By the rank condition (22) of Assumption 3, the matrix $\mathcal{I}\mathcal{I}^{\top}>0$. It follows from Eq. (66) that $M_{1}={0}$, $M_{2}={0}$, and $M_{3}={0}$. As a result, $\mathcal{L}_{[K;\mathcal{S}(\alpha)]}(Y_{1})={0}$, together with $K\in\mathcal{Z}(\alpha)$, which implies that the only solution is $Y_{1}={0}$. Substituting into $M_{2}=0$ and $M_{3}=0$ yields $Y_{2}={0}$ and $Y_{3}={0}$. Hence, $\mathcal{Y}={0}$, contradicting the assumption $\mathcal{Y}\neq 0$. Hence, the matrix $\varPhi(\alpha,K)$ has full column rank. $\square$

.3 Proof of Theorem 2

1). We first show that $0<P_{i}<\infty$, $K_{i}\in\mathcal{Z}(\alpha_{i})$, and $\alpha_{i}>\alpha_{i-1}$ hold for all $1\leq i\leq\mathbb{I}$ by mathematical induction.

1. i)

Suppose $i=1$. Given $\gamma>\sigma(\mathcal{L}_{[0;\mathcal{S}]})+\alpha_{0}$, we have $K_{0}=0\in\mathcal{Z}(\alpha_{0})$ (Lemma 2). By Lemma 1 (Statement 1), there exist a unique positive definite solution $0<P_{1}<\infty$ to Eq. (23), and the solutions $K_{1}$ and $\alpha_{1}$ of Eqs. (24) and (25), respectively, satisfy that $K_{1}\in\mathcal{Z}(\alpha_{1})$ and $\alpha_{1}>\alpha_{0}$. 2. ii)

Assume for some $1<i<\mathbb{I}$, $K_{i}\in\mathcal{Z}(\alpha_{i})$ holds. Apply Lemma 1(Statement 1) again at step $i+1$, Eq. (23) admits a unique solution $0<P_{i+1}<\infty$, the updated feedback gain $K_{i+1}$ belongs to $\mathcal{Z}(\alpha_{i+1})$ with $\alpha_{i+1}>\alpha_{i}$.

Therefore, the sequence $\{\alpha_{i}\}$ is strictly increasing and bounded above by $\gamma$. By construction, the stopping index $\mathbb{I}$ satisfies $\gamma\leq\alpha_{\mathbb{I}}$ and $\alpha_{i}<\gamma,\forall 1<i<\mathbb{I}$. Also note that since each $P_{i}$ is bounded, there exist positive constants $c_{0},c_{1}>0$ such that $\{\lambda_{\max}(P_{i})\}\leq c_{0}$ and $\lambda_{\min}(\|K_{i}\|_{R}^{2}+\Theta)\geq c_{1}$, $\forall 1<i<\mathbb{I}$. Applying these bounds to Eq. (25), we obtain $\alpha_{i}-\alpha_{i-1}\geq\eta c_{1}/c_{0}$. A summation from $i=1$ to $i=\mathbb{I}-1$ gives

[TABLE]

which implies that

[TABLE]

This proves the finiteness of $\mathbb{I}$. Then, at iteration $i=\mathbb{I}$, Eq. (23) becomes

[TABLE]

Express it in terms of the original system $\mathcal{S}$ and ${K}_{\mathbb{I}}$ as

[TABLE]

where $\Theta_{\mathbb{I}}:=\|K_{\mathbb{I}}\|_{R}^{2}+\Theta+(\alpha_{\mathbb{I}-1}-\alpha_{\mathbb{I}})P_{\mathbb{I}}$. The positive definiteness of $\Theta_{\mathbb{I}}$ directly follows from Eq. (25) and the facts that $P_{\mathbb{I}}>0$ and $\alpha_{\mathbb{I}}>\alpha_{\mathbb{I}-1}$. Finally, since $\alpha_{\mathbb{I}}\geq\gamma$, by [31, Theorem 1], it can be immediately concluded that $K_{\mathbb{I}}\in\mathcal{Z}(\gamma)$.

2. Next, we verify that for all $i\geq\mathbb{I}+1$, $P_{i}\geq P^{*}$ and $K_{i}\in\mathcal{Z}(\gamma)$ by mathematical induction.

1. i)

Suppose $i=\mathbb{I}+1$. Since $K_{\mathbb{I}}\in\mathcal{Z}(\gamma)$, based on Lemma 1(Statement 2), we know that Eq. (26) admits a unique solution $P_{\mathbb{I}+1}\geq P^{*}$ and Eq. (24) gives $K_{\mathbb{I}+1}\in\mathcal{Z}(\gamma)$. 2. ii)

Assume for some $i>\mathbb{I}+1$, $K_{i}\in\mathcal{Z}(\gamma)$ holds. Apply Lemma 1(Statement 2) again at step $i+1$, the solutions to Eqs. (26) and (24) satisfy $P_{i+1}>0$ and $K_{i+1}\in\mathcal{Z}(\gamma)$.

To establish the convergence result, we derive a recursive identity from Eqs. (26) and (24) that

[TABLE]

Since $\sigma(\mathcal{L}_{[K_{i};\mathcal{S}]})\subset\mathbb{C}^{-}$ has been established, this equation implies $P_{i}\geq P_{i+1}$ for all $i\geq\mathbb{I}+1$. As a result, the sequence $\{P_{i}\}_{i=\mathbb{I}+1}^{\infty}$ is monotonically decreasing and bounded below by $P^{*}$. Therefore, the limit exists, i.e., there exists $P_{\infty}$ such that $\lim_{i\rightarrow\infty}P_{i}=P_{\infty}$. Under Assumption 1, the solution $P^{*}\geq 0$ of (14) is unique and satisfies (26) along with $K^{*}$ satisfying (24). Thus, it can be concluded that $P_{\infty}=P^{*}$ and $K_{\infty}=K^{*}$. $\square$

.4 Proof of Theorem 3

1. Assume $K_{i-1}\in\mathcal{Z}(\alpha_{i-1})$. We have that the unique solution triplet $(P_{i},K_{i},\alpha_{i})$ to Eqs. (23)-(25) also satisfies (39), (40), and (25). On the other hand, Lemma 3 ensures that Eq. (39) has a unique solution under Assumption 3. Therefore, the iterative procedure defined by Eqs. (39), (40), and (25) is equivalent to solving (23)-(25). From the proof of Theorem 2 (Statement 1), the condition $K_{i}\in\mathcal{Z}(\alpha_{i})$ is verified for all $1\leq i\leq\mathbb{I}$, starting from $K_{0}=0\in\mathcal{Z}(\alpha_{0})$. Hence, the conclusion of property 1 follows directly from Theorem 2 (Statement 1).

2. We now assume $K_{i-1}\in\mathcal{Z}(\gamma)$ holds. Under this assumption, Lemma 3 ensures that Eq. (42) admits a unique solution under Assumption 3. Moreover, the unique solution pair $(P_{i},K_{i})$ obtained from Eqs. (26) and (24) also satisfies (42) and (40), which implies that the policy iteration defined by (42)and (40) is equivalent to solving (26) and (24). By using the obtained result $K_{\mathbb{I}}\in\mathcal{Z}(\gamma)$ and Theorem 2(Statement 2), we immediately conclude that ${K}_{i}\in\mathcal{Z}(\gamma)$ for all $i\geq\mathbb{I}+1$, and thus the convergence results are established. $\square$