<figure><img src="image_1.png" /></figure>

HAL
open science

# A Multistep Reinforcement Learning Control of Shear Flows in Minimal Input-Output Plants Under Large Time-delays

Amine Saibi, Lionel Mathelin, Onofrio Semeraro

► To cite this version:

Amine Saibi, Lionel Mathelin, Onofrio Semeraro. A Multistep Reinforcement Learning Control of Shear Flows in Minimal Input-Output Plants Under Large Time-delays. Flow, Turbulence and Combustion, 2025, 115 (3), pp.1379-1402. 10.1007/s10494-025-00697-w. hal-05379450

HAL Id: hal-05379450

https://hal.science/hal-05379450v1

Submitted on 24 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A multistep reinforcement learning control of
shear flows in minimal input–output plants under
large time-delays

Amine Saibi¹,²*, Lionel Mathelin² and Onofrio Semeraro²

¹Inst. Jean Le Rond d'Alembert, Sorbonne Université, Paris, 75005, France.

²Lab. Interdisciplinaire des Sciences du Numérique (LISN), CNRS, Univ. Paris-Saclay, Orsay, 91400, France.

*Corresponding author(s). E-mail(s): amine.saibi@dalembert.upmc.fr;
Contributing authors: lionel.mathelin@cnrs.fr; onofrio.semeraro@cnrs.fr;

Abstract

Flow control has attracted research for its potential role in reducing drag, sup-
pressing turbulence, and enhancing mixing in fluid systems. The emergence of
data-driven modeling and machine learning techniques has sparked new inter-
est in designing control strategies that can adapt in real time to complex,
high-dimensional flow environments. However, fluid systems remain particularly
challenging testbeds for control design due to their nonlinear and convective
nature, which introduces large time delays. In active control, additional difficul-
ties arise from practical constraints, such as the use of localized sensors in limited
number. In this work, we investigate a reinforcement learning framework based
on a suitable actor-critic algorithm designed to address these challenges. Two
test cases representative of transitional shear flows are considered: a linearized
version of the Kuramoto-Sivashinsky equation and the control of instabilities in
a two-dimensional boundary-layer flow over a flat plate, using a minimal but
realistic sensor-actuator configuration. This choice reflects our focus on the limi-
tations that arise from plants of experimental interest. Time delays are identified
during a pretraining stage, while the control algorithm employs multistep returns
during value iteration. This approach improves both the convergence rate and
stability of learning. Furthermore, we show that the lookahead in the multistep
formulation provides a non-trivial beneficial effect in plants where the control
task is characterized by a severe credit-assignment issue.
# 1 Introduction

Reinforcement learning is gaining attention due to its indisputable successes across a wide range of scientific fields, from robotics to the optimization of nonlinear, complex systems such as those found in fluid mechanics, where control has long been a major research focus. Optimizing aerodynamic flows is well recognized for its potential to reduce pollutant emissions, mitigate acoustic noise, and control separation or stall. An introduction to the topic can be found in [1]; among others, reviews focusing on specific aspects include [2] for linear control, [3] for input-output design and reduced-order modeling, and [4] for the influence of flow stability on closed-loop dynamics. A broader perspective on nonlinear control applied to fluid flows is provided in [5]. Flow optimization can be achieved through passive or active manipulation; in this work, we focus on active control, where localized actuators act as transducers to manipulate the system [6], modulated by control policies based on sensor measurements. Depending on the flow regime, control mechanisms may aim to attenuate disturbances responsible for transition to turbulence at low Reynolds numbers [4], or to modify large-scale structures in turbulent flows at higher Reynolds numbers [7].

Standard applications in flow control typically rely on reduced-order modeling as a preliminary step in control design. Common approaches include approximate balanced truncation [8, 9], system identification techniques such as the auto-regressive with moving average (ARMAX) method and related algorithms [see 10, Chapters 4 and 7], and subspace iteration methods such as the subspace state-space system identification (N4SID) algorithm [see 11, Chapters 4 and 5]. Examples of applications in fluid mechanics are reported in [12, 13]. These techniques often rely on restrictive assumptions—such as linearized dynamics, low-dimensional approximations, or perfect state information—that can introduce significant modeling uncertainties or degrade performance in practice. Consequently, when applied to more realistic flow configurations, it becomes essential to explicitly account for robustness to model inaccuracies and measurement noise, to ensure the stability of the resulting closed-loop system in feedback configurations [4]. These limitations can, in principle, be addressed using data-driven tools based on measurements from local sensors, without relying on *a priori* physical models, by introducing reinforcement learning (RL) policies. This framework encompasses numerous techniques, including algorithms rooted in control theory through dynamic programming (DP), a nonlinear optimization protocol based on the Bellman equation, whose solution is the value function. In principle, once the value function is known, the optimal policy can be determined. However, the direct solution of the Bellman is impractical in most cases; a large body of literature is dedicated to its solution, and in particular the approximation based on iterative methods that can be applied with and without a prior model at hand: from a theoretical viewpoint, these are the premises on top of which RL algorithms stand. For comprehensive overviews, we refer the reader to the classic works by Bertsekas [14] and Sutton and Barto [15]. In short, a large part of the RL algorithms are iterative, data-driven and solely relying on limited measurements. The control policy is determined from the exploration of the state-action space, the set of all state-action pairs, each consisting of a possible state of the environment and a valid action the agent can take in that state. In the limit of full knowledge of this space, the resulting policy is obtained by optimizing the
<figure><img src="image_3.png" /><figcaption>**Fig. 1:** Sketch of the control configuration, where both inputs and outputs are localized along the streamwise direction. This setup is used to control the dynamics governed by the linearized Kuramoto–Sivashinsky equation, as well as disturbances evolving in the two-dimensional boundary layer flow over a flat plate. The system inputs, denoted by $d$ and $u$, represent external disturbances and actuation, respectively. The outputs, $y$ and $z$, correspond to the controller’s observation and the reward or cost function. The spatial arrangement of sensors and actuators introduces intrinsic time delays that significantly affect the behavior of the closed-loop system. Further details are provided in section 2.</figcaption></figure>

value function, also denoting the cost function when minimized or the reward when maximized [14, 16]. Interestingly, this framework also provides a basis for generalizing optimal control problems across linear and nonlinear settings. Linear-quadratic controllers [16, Chapters 3 and 11], model-predictive controllers (MPC) [17, Chapter 16], and adaptive controllers [16, see Chapter 11, following 18] can all be regarded as special cases of this broader formulation. A recent overview from a control perspective is provided by [19], which highlights the links between reinforcement learning and standard control, addresses the challenges they pose, and discusses issues often overlooked due to prevailing misconceptions.

Although far from being exhaustive, this quick sketch provides a motivation for the success of these applications in fluid mechanics; among the early works, it is worth mentioning the ones by [20–23], followed shortly after by a large amount of work that pushed the use of these techniques to progressively more complex environments [24, 25]. A recent review on the latest development is provided by [26], while a perspective on the advantages of the RL framework in fluid mechanics, when compared to standard techniques, is given by the analysis in [27]. The main aim of this article is to demonstrate how RL can be applied to flow control in the presence of large time delays, characteristic of convective flows with localized sensing and actuation. The contribution is directly inspired by [28] and, more recently, [29], where the
convective nature of shear flows was shown to critically impact the input-output con-
figuration of controllers. We focus on the same test cases: the control of a modified
Kuramoto-Sivashinsky (KS) equation and of Tollmien-Schlichting waves in a two-
dimensional boundary layer over a flat plate. A sketch is provided in fig. 1, showing the
layout of the input-output configuration; the system inputs are the forcing signals *d*
and *u*, modeling disturbances and actuation, respectively; the outputs are denoted by
*y* and *z*, corresponding to observations available to the controller and the reward/cost
function. These configurations serve as a testbed for the control of perturbations with
the goal of delaying transition to turbulence [3, 30–32]. With localized sensors and
actuators, their relative placement along the flow direction has a decisive impact on
closed-loop behavior. As observed in experimental studies [4, 28, 33–35], the standard
configuration depicted in fig. 1 leads to a dynamically open-loop system, effectively
equivalent to a disturbance-feedforward control strategy [36], well documented in both
robust control [37, Chapter 12] and fluid mechanics [4, 38]. In this regime, the con-
troller relies only on disturbance measurements and the design model, whereas true
feedback arises when the sensor captures the actuator’s influence, which directly lim-
its performance and robustness—constraints imposed by the plant itself, especially in
noisy or uncertain environments. Several approaches have been proposed to address
this issue. One option is to preserve the dynamics between estimation and actuation,
as in [35]; while this retains robustness, performance quickly deteriorates. Another
approach recovers information through a secondary path using algorithms from active
noise cancellation, such as adaptive filters [28, 30]. More equipment-intensive strategy
employs multiple sensors upstream and downstream of the actuator.

It follows that RL algorithms will be similarly affected when applied to realis-
tic input-output configurations, particularly in view of experimental implementation.
These aspects are often overlooked; most prior studies have focused on bluff bodies or
shear flows with dense sensing: in both cases, the impact of time delays is less criti-
cal. Here we deliberately retain a minimal and experimentally feasible configuration
(as in fig. 1) and introduce a formal control framework based on an actor-critic algo-
rithm. While the governing equations are simplified through linearization, the adopted
input-output configuration is more realistic than approaches relying on dense sens-
ing or impractically long training horizons. In doing so, we isolate the effects of time
delays, sensor limitations, and their interaction with multistep policy iteration and
credit assignment. Recent works have introduced multistep methods without explicit
reference to physical considerations [39, 40], while treatments of time delays appear
in [41] and partially in [29], where deep deterministic policy gradient (DDPG) [42, 43]
is combined with reward shaping to account for both time delays and the disturbance
frequency bandwidth. However, to our knowledge, no prior study has examined the
interplay between delays, credit assignment, and convergence in flow control.

The remainder of the paper is organized as follows. In section 2, we introduce
the two demonstrators—the KS equation and the two-dimensional boundary-layer
flow—along with a concise description of the numerical setup. The reinforcement
learning algorithm is detailed in section 3 and section 4, including an analytical con-
vergence analysis extending [44] within the framework proposed in [45]. An in-house
code is used. We provide only the main ideas and results here to maintain clarity for
a fluid mechanics audience. The control of instabilities in both test cases, achieved with a minimal number of sensors and actuators, is presented in section 5. The paper concludes with a discussion of the findings and perspectives in section 6.

## 2 Numerical setup and plant

In this section, we introduce the systems used for the numerical demonstrations. The linearized version of the Kuramoto–Sivashinsky (KS) equation is introduced and discussed in [28], where a fine-tuning of the flow parameters is carried out for obtaining an input–output behavior comparable to the ones observed in two-dimensional boundary-layer flows developing over a flat plate [3, 29]. We consider a realistic input–output configuration by modeling localized and limited sensors and actuators, and by explicitly accounting for time delays. Although the environment consists of linearized models, the chosen plant more accurately reflects practical experimental conditions — particularly for the specific task of delaying transition to turbulence — from an input–output viewpoint as compared to highly idealized scenarios that assume an excessive number of sensors.

### 2.1 Kuramoto–Sivashinsky equation and input–output system

Given the state $q: \mathbb{R}^+ \times [0, L] \to \mathbb{R}$, the stochastic partial differential equation associated to a modified version of the KS equation to be controlled reads as

$$dq(t, x) = (\mathcal{L}q(t, x) + \phi(x; x_a)u(t)) dt + \sigma_a \phi(x; x_a) dW_a(t) + \phi(x; x_d) dW_d(t), \quad (1)$$

with the differential operator $\mathcal{L} = -\mathcal{V}\partial_x - \frac{1}{R}(\mathcal{P}\partial_x^2 + \partial_x^4)$. The adopted parameters are $R = 0.25$, $\mathcal{P} = 0.05$, and $\mathcal{V} = 0.4$, corresponding to a convective velocity, following [28]. The boundary conditions are defined as $q(t, 0) = \partial_x q(t, 0) = \partial_x q(t, L) = 0$ for all $t \ge 0$ with a fixed domain length $L = 800$. $W_a(t)$ and $W_d(t)$ are two independent standard Wiener processes adapted to a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \in \mathbb{R}^+}, \text{Pr})$.

The inputs are applied to eq. (1) using a spatial function $\phi(x; x_i) = e^{-(x-x_i)^2}$, localizing the actuator and external-disturbance effects at the spatial points set to $x_a = 400$ and $x_d = 40$, respectively. The outputs are defined as

$$\left\{
\begin{aligned}
y(t) &= \int_{[0, L]} \phi(x; x_s) q(t, x) dx + \nu_y(t), \\
z(t) &= \int_{[0, L]} \phi(x; x_z) q(t, x) dx + \nu_z(t),
\end{aligned}
\right.
\qquad (2)$$

representing, respectively, the observation (or estimation) and cost sensors, centered on the space points $x_s = 300$ and $x_z = 450$, with $\mathcal{F}_t$-adapted additive Gaussian white noises $\nu_y(t) \sim \mathcal{N}(0, \sigma_s)$ and $\nu_z(t) \sim \mathcal{N}(0, \sigma_z)$. We set the standard deviations of the noises to $\sigma_a = \sigma_y = \sigma_z = 10^{-1}$.

The partial differential equation is approximated numerically using Fourier discretization based on an equispaced partition of $[0, L]$ with $\Delta x = L/N$ and $N = 1600$
discretization points. The domain is made periodic using a fringe region [46], covering 20% of the domain length and leading to a shorter physical domain. The time discretization is performed by introducing a Crank–Nicolson scheme with a time-step $\Delta t = 0.1$. Given these choices for the numerical approximation, by defining the Fourier transform of a function $f \in L^2([0, L])$ as

$$ \mathcal{F}(f)(k) = \hat{f}(k) = \int_{[0,L]} f(x)e^{-\frac{2\pi k}{L}kx}dx, \quad k \in \mathbb{Z}, \quad (3) $$

the differential operator $\mathcal{L}$ writes as follows

$$ l(k) = \frac{1}{\mathcal{R}} \left[ \left( \frac{2\pi k}{L} \right)^4 - \mathcal{P} \left( \frac{2\pi k}{L} \right)^2 \right] - i\mathcal{V} \frac{2\pi k}{L}. \quad (4) $$

The application of the Fourier transform in space combined with the Crank–Nicolson scheme yields

$$ \hat{q}(t + \Delta t, k) = \left(1 - \frac{\Delta t}{2} l(k)\right)^{-1} \left[ \left(1 + \frac{\Delta t}{2} l(k)\right) \hat{q}(t, k) + \hat{\phi}(k, x_a)(u(t) + \nu_a(t)) \right. \\ \left. + \hat{\phi}(k, x_d)\nu_d(t) \right], \quad (5) $$

where the noise applied to the right-hand side is discretized as

$$ \nu_a(t) = \sigma_a (W_a(t + \Delta t) - W_a(t)) \sim \mathcal{N}(0, \sigma_a \sqrt{\Delta t}), $$

and

$$ \nu_d(t) = W_d(t + \Delta t) - W_d(t) \sim \mathcal{N}(0, \sqrt{\Delta t}). $$

The resulting dynamics is shown in fig. 2, where the spatio-temporal contour plot shows the evolution of the open-loop dynamics governed by the linearized, convective version of the KS equation.

### 2.1.1 Input–output system

Taking the inverse Fourier transform of eq. (5), the discrete operator can be written as

$$ Aq = \mathcal{F}^{-1} \left[ \left(1 - \frac{\Delta t}{2} l\right)^{-1} \left(1 + \frac{\Delta t}{2} l\right) \mathcal{F}(q)(t, \cdot) \right], $$

for a one-step time marching. By denoting in the following for all $n \in \mathbb{N}$, $q_n$ as the random vector $(q(n\Delta t, i\Delta x))_{i \in \{1..N\}}$, the control variable $u_n = u(n\Delta t) + \nu_a(t)$, and the observation and cost variables $y_n$ and $z_n$, we get the input–output discrete system

$$ \left\{ \begin{array}{l} q_{n+1} = Aq_n + B_u u_n + B_d \nu_d(n), \\ y_n = C_y q_n + \nu_y(n), \\ z_n = C_z q_n + \nu_z(n), \end{array} \right. \quad (6) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2: Spatio-temporal contour plot showing the evolution of the dynamics governed by the linearized, convective version of the Kuramoto–Sivashinsky equation.</figcaption></figure>

where the discretization of the integrals in eq. (2) yields 1 × N matrices $C_y, C_z \in \mathcal{M}_{1,N}$ to model the sensors, while the inputs are given by

$$B_i = \mathcal{F}^{-1} \left[ \left( 1 - \frac{\Delta t}{2} l \right)^{-1} \hat{\phi}(\cdot; x_i) \right] \in \mathcal{M}_{N,1}(\mathbb{R}),$$

with $x_i$ the local position of the elements. This is the mathematical representation of the input-output system sketched in fig. 1.

## 2.2 Boundary-layer flow

Formally, it is possible to obtain an input-output system as eq. (6) by discretizing the linearized 2-D incompressible Navier-Stokes (NS) equations. Given $\Omega = ]0, L[ \times ]0, H[$ with $L = 1000$ and $H = 40$, we introduce the perturbation velocity $\mathbf{u}' : [0, T] \times \Omega \to \mathbb{R}^2$; the linearized NS governing the perturbation dynamics around a base flow $\mathbf{U}$ reads

$$ \left\{ \begin{array}{l} \mathrm{d}\mathbf{u}'(t, x_0, x_1) = (\mathcal{D}\mathbf{u}'(t, x_0, x_1) + \phi(x_0, x_1; x_a)u(t))\mathrm{d}t \\ \qquad + \sigma_a\phi(x_0, x_1; x_a)\mathrm{d}W_a(t) + \phi(x_0, x_1; x_d)\mathrm{d}W_d(t), \quad (7) \\ \mathrm{div}\,\mathbf{u}'(t, x_0, x_1) = 0, \end{array} \right. $$

where the perturbation velocity $\mathbf{u}'(t, x) = (u'(t, x), v'(t, x))$ has two components representing the streamwise and normalwise velocities, respectively. $\mathcal{D}$ denotes the differential operator action defined by

$$\mathcal{D}\mathbf{u}' = -\mathbf{U} \cdot \nabla\mathbf{u}' - \mathbf{u}' \cdot \nabla\mathbf{U} + \frac{1}{Re_{\delta_0^*}} \Delta\mathbf{u}' - \nabla p,$$
where the Reynolds number $Re_{\delta_0^*} = U\delta_0^*/\nu = 1000$ is introduced based on the displacement thickness $\delta_0^*$ of the boundary layer computed at $x = 0$, the inlet velocity $U$ and the kinematic viscosity $\nu$. The variable $p$ indicates the perturbation pressure.

The function $\phi(x_0, x_1; p) = (1, 1)e^{-(x_0-p_0)^2-(x_1-p_1)^2}$ in eq. (7) is extended to the two-dimensional setting. The adaptation of the observation and cost sensors, as defined in eq. (2), follows naturally. In analogy with the input-output system based on the KS equation, we set the positions as $x_d = (40, 0)$, $x_s = (300, 0)$, $x_a = (400, 0)$, and $x_z = (450, 0)$, with noise-scaling parameters $\sigma_a = \sigma_y = \sigma_z = 10^{-4}$ also for the boundary-layer case.

Numerical simulations of eq. (7) are carried out by a code using the framework provided by FEniCS [47], based on finite elements for the spatial discretization and marched in time using Crank-Nicolson scheme. The boundary conditions imposed are Dirichlet along the flat plate $u'(t, x_0, 0) = (0, 0) \forall (t, x_0) \in [0, T \times ]0$, $L$, and the top boundary $u'(t, x_0, H) = (0, 0) \forall (t, x_0) \in [0, T \times ]0$, $L$; null perturbations are imposed at the inlet and outflow conditions and null pressure at the outlet. Note that the baseflow $U$ is obtained by running the nonlinear version of the code, where a solution of the Blasius equation is prescribed at the inlet [see 48, chapter 8].

## 3 Markov Decision Processes and the Bellman operator framework

We denote by $\mathcal{P}(S)$ the set of all probability distributions over a set $S$, and assume that all functions considered are measurable. Let $\mathcal{X} \subset \mathbb{R}^d$ and $\mathcal{A} \subset \mathbb{R}^p$ be subsets representing the state space and action space, respectively, where $d$ and $p$ two positive integers. We model the sequence of state-action pairs $(q_n, u_n)_{n \in \mathbb{N}}$ of stochastic processes as a Markov decision process (MDP) taking values in $\mathcal{X} \times \mathcal{A}$. The time evolution of this process is governed by a transition kernel $P: \mathcal{X} \times \mathcal{A} \to \mathcal{P}(\mathcal{X})$ which is derived via the inverse Fourier transform of equation (5). This kernel describes the stochastic law that determines the next state $q_{n+1}$ resulting from the application of an action $u_n$ when being in a state $q_n$ at $n \in \mathbb{N}$, and is formally given by:

$$ \forall \psi: \mathcal{X} \to \mathbb{R} \text{ bounded, } \quad \mathbb{E}[\psi(q_{n+1} | q_n, u_n)] = P\psi(q_n, u_n) := \int_{\mathcal{X}} P(q_n, u_n, dx)\psi(x). $$

We call control strategy or policy, a map $\pi: \mathcal{X} \to \mathcal{P}(\mathcal{A})$ that samples the action $u_n$ to input at each $n \in \mathbb{N}$ while being in a state $q_n$, i.e.:

$$ \forall \zeta: \mathcal{A} \to \mathbb{R} \text{ bounded, } \quad \mathbb{E}[\zeta(u_n) | q_n] = \pi\zeta(q_n) := \int_{\mathcal{A}} \pi(q_n, da)\zeta(a), $$

with $\mathbb{E}[\cdot | Y]$ denoting the conditional expectation with respect to a variable $Y$. The data of an initial state $q_0 \in \mathcal{X}$ and the couple of maps $(P, \pi)$ uniquely determine the process $(q_n, u_n)_{n \in \mathbb{N}}$, we will denote by $\mathbb{P}_{q_0, P, \pi}$ its probability distribution. We will also often write $\mathbb{P}_{q_0, u_0, P, \pi}$ to designate the distribution of the MDP with an initial state $q_0 \in \mathcal{X}$ and a first action $u_0 \in \mathcal{A}$, i.e., $u_0 = a_0$ almost surely.
This definition implies that the process $q_n$ is a Markov chain with transition law $\pi P(q, dx) := \int_A \pi(q, da) P(q, a, dx)$, and $(q_n, u_n)$ is also a Markov chain with transition law $P\pi(q, u, dx, da) := P(q, u, dx)\pi(x, da)$.

To set an optimal control problem, we define the (state-action) value function of a policy $\pi$, as the expected cumulated (discounted) costs:

$$Q_{\pi}(q, u) = \mathbb{E}_{\mathbb{P}_{q,u,P,\pi}} \left[ \sum_{n \ge 0} \gamma^n z_{n+1}^2 \right], \quad (q, u) \in \mathcal{X} \times \mathcal{A},$$

with $\gamma \in ]0, 1[$ a discount factor. The cost variable $z_{n+1}^2 = C(q_n, u_n, q_{n+1})$ is seen as a function of the transition $(q_n, u_n, q_{n+1})$. The optimal control problem in reinforcement learning consists in finding a policy $\pi^*$ that minimizes the value function of an MDP, i.e., for all $q, u \in \mathcal{X} \times \mathcal{A}$:

$$Q_{\pi^*}(q, u) = \inf_{\pi} Q_{\pi}(q, u) := Q^*(q, u),$$

where the infimum is taken along all policies $\pi$.

Let $\zeta: \mathcal{X} \times \mathcal{A} \to \mathbb{R}$ be a bounded function and $\pi$ a policy. We define the Bellman operators $\mathcal{T}$ and $T_\pi$ on $L^\infty(\mathcal{X} \times \mathcal{A})$:

$$\mathcal{T}\zeta(q, u) := c(q, u) + \gamma P \left[ \min_{a \in \mathcal{A}} \zeta(\cdot, a) \right] (q, u), \quad (8)$$

$$T_{\pi}\zeta(q, u) := c(q, u) + \gamma P\pi\zeta(q, u), \quad (9)$$

where $c(q, u) = PC(q, u)$. In the following, we will also write $c_n$ to designate the cost process $c(q_n, u_n)$. The Bellman operators are $\gamma$-Lipschitz with respect to the supremum norm $\|\cdot\|_{\infty}$, i.e., taking any two bounded functions $\zeta$ and $\psi$ we have

$$\begin{aligned} \|T_{\pi}\zeta - T_{\pi}\psi\|_{\infty} &= \sup_{(q,u) \in \mathcal{X} \times \mathcal{A}} \left| \gamma \int P(q, u, dx) \pi(x, da) (\zeta(x, a) - \psi(x, a)) \right|, \\ &\le \gamma \|\zeta - \psi\|_{\infty}. \end{aligned} \quad (10)$$

Since $\gamma < 1$, the Banach fixed-point theorem then implies that there exists a unique fixed point to a Bellman operator. Note that we can write the value function

$$Q_{\pi} = \sum_{n \ge 0} \gamma^n (P\pi)^n c, \quad (11)$$

with $(P\pi)^n$ the $n$-iterate of $P\pi$, where for all bounded $\psi: \mathcal{X} \times \mathcal{A} \to \mathbb{R}$,

$$P\pi\psi(q, u) := \int_{\mathcal{X} \times \mathcal{A}} P(q, u, dx)\pi(x, da)\psi(x, a).$$
Hence, it is possible to derive the Bellman equations $T_\pi Q_\pi = Q_\pi$ and $\mathcal{T}Q^* = Q^*$ identifying precisely the fixed points of the Bellman operators. The optimal control problem can then also be described as finding $\pi^*$ such that

$$T_{\pi^*} Q_{\pi^*} = \mathcal{T} Q_{\pi^*}. \qquad (12)$$

The processes $(q_n, y_n, u_n)_{n \in \mathbb{N}}$ generally describe a partially observable Markov decision process (POMDP) when considering control processes $u_n$ that are output feedbacks, i.e., measurable w.r.t. $y_n$. The observability of eqs. (1) and (2) is determined solely by the (discretized) differential operator $\mathcal{L}$ and the function $\phi(\cdot, x_s)$ allowing the system to be decomposed into fully observable and unobservable dynamics [see 49, Chapter 14]. The observable part of the state, denoted $q^\theta$, evolves over time, with the influence of unobserved states appearing as an additional forcing term. Its transition dynamics can therefore be modeled using a kernel $P_\theta(q_n^\theta, u_n, dq_{n+1}^\theta)$. Let $\mathcal{O} \subset \mathcal{X}$ denote the observable subspace. Takens' theorem [50] guarantees the existence of a time-delay embedding of observables $\theta: \mathcal{O} \to \mathbb{R}^{\tau+1}$, which maps $q_n^\theta$ to the vector $(y_{n-\tau}, \dots, y_{n-1}, y_n)$ for a delay $\tau \in \mathbb{N}$ that can be estimated from measurements. The same work also shows that the delays used in the embedding need not be equally spaced; any set of generic lags with cardinality $\tau+1$ suffices. Subsequent refinements in [51] extend these ideas to non-autonomous (forced) and stochastic systems, establishing that the delay-coordinate reconstruction remains valid in the presence of external forcing and measurement noise. We then consider policies of the form $\pi(y, du)$, where $y = \theta(q)$, so that the policy depends on observations. The pair $(P_\theta, \pi)$ thus defines a Markov decision process (MDP). A key observation is the need to adapt the value function accordingly. This is done by accounting for the change in the probability law of the process $z_{n+1} = C(q_n, u_n, q_{n+1}, \nu_z)$, resulting from replacing $q_n$ with $q_n^\theta$. The lost information is treated as additional noise, effectively reducing the POMDP to an MDP and aligning with the introduced optimal control formulation.

# 4 An actor-critic algorithm for time-delayed systems

The control design described in this section consists of two parts. The first is a pre-training procedure that estimates the time delay between the sensor and actuator using a generalized cross-correlation method, as outlined in section 2. The second part involves an actor-critic algorithm enhanced with a multistep return strategy; in section 4.2, we present a formal analysis with theoretical derivations to highlight the role of the multistep approach from a convergence perspective. Practical considerations and algorithmic details are provided in section 4.2.1.

## 4.1 Estimating time delay

Time delays characterizing the input-output system can be estimated solely relying on available data. To this end, let $(y(t))_{t \in [0,T]}$ and $(z(t))_{t \in [0,T]}$ be two finite-time trajectories of observations and cost stochastic processes, respectively. As an illustration, we consider the KS equation presented in section 2.1. We define, in continuous time,
<figure><img src="image_5.png" /><figcaption>Fig. 3: Sample signals from the observation sensor $y(t)$ (black line) and the cost sensor $z(t)$ (red line), shown in an open-loop configuration (left) and under active control based on the observation signal (right).</figcaption></figure>

the cross-correlation of $(y, z)$ as

$$ \chi_{y,z}(\tau) = \int_{\tau}^{T} y(t)z(t - \tau)dt. \qquad (13) $$

An estimation of the time delay between $y$ and $z$ is then given by (see [52])

$$ \hat{\tau} = \arg \max_{\tau \in [0,T]} \chi_{y,z}(\tau). \qquad (14) $$

In our setting, to estimate the time delay between the upstream sensor and the actuator we generate two sets of trajectories $(y^{(i)}, u^{(i)}, z^{(i)})_{i \in \{0,1\}}$ with $u^0 \sim \delta_0$ being a null control and $u^{(1)} \sim \mathcal{N}(y^{(1)}, \sigma)$ reproducing the sensor signal $y^{(1)}$ with a standard deviation set to $\sigma = 10^{-1}$. The generated signals are shown in fig. 3, where it is also possible to observe the time shifts between the signals. The estimated time delay $\hat{\tau}$ between the sensor and actuator is then given by

$$ \hat{\tau} = \hat{\tau}_0 - \hat{\tau}_1, \qquad (15) $$

with $\hat{\tau}_i = \arg \max_{\tau \in [0,T]} \chi_{y^{(i)},z^{(i)}}(\tau)$, $i \in \{0,1\}$.

The rationale for taking the difference is that the first scenario (null control) establishes a baseline delay $\hat{\tau}_0$ between the observation and cost sensors, while the second scenario (noisy actuation on sensing) allows us to pinpoint how the upstream signal directly influences the downstream cost upon actuation, thus establishing the delay between actuator and cost sensor. In fig. 4, the results of the extrapolation of the time delays based on this procedure are shown. Statistics are performed on 10 sample trajectories for each sensor position, ranging from $x = 50$ to $x = 390$ with a space step of $\Delta x = 20$, and compared to the reference delay $\tau = (x_a - x_s)/V$, where $V$ is the velocity defined before. We consider a setup where the actuator is fixed in $x_a = 400$
<figure><img src="image_6.png" /><figcaption>Fig. 4: On the left, the relative error between the estimated time delay $\hat{\tau}$ (see section 4.1) and the actual time delay $\tau$ is shown as a function of the observation sensor position $x_s$. The evaluation is based on 10 sample trajectories with different random seeds and shows consistent time-delay estimates up to $x_s \approx 330$. The increase in relative error beyond this point is attributed to the combined effects of time resolution and the smaller magnitude of $\tau$. Nonetheless, as shown on the right, where both the estimated and actual time delays are plotted as functions of $x_s$, the estimation remains well-behaved.</figcaption></figure>

and the reward sensor is in $x_z = 450$. In fig. 4, left, the relative error between the estimation of the time delay $\hat{\tau}$ and the actual time delay $\tau$ is shown as a function of the position of the observation sensor placed in $x_s$. Results are consistent up to $x_s \approx 330$ when the relative error increases due to time resolution and the smaller value of $\tau$; nonetheless, the estimated time delays plotted as a function of $x_s$ on the right show good agreement with the real time delays. A similar procedure is applied for the case of the boundary layer in section 5, where the convective velocity used for the estimation of the time delay corresponds to the group velocity of the propagating Tollmien-Schlichting wavepackets.

## 4.2 Approximate Policy Iteration in Actor-Critic methods: monotonicity and error Bounds

To address the optimal control problem for the MDP, we place ourselves in the actor-critic framework which consists in building a sequence $(\pi_k, Q^k)_{k \in \mathbb{N}}$ of pairs with the aim of eventually converging to the optimal policy $\pi^*$ and optimal value function $Q_{\pi^*}$. Dynamic programming provides many ways of building such a sequence [see 15, chapter 4]. A classical approach alternates between two procedures, policy improvement and policy evaluation [45].

Consider the sequence $(\pi_k, Q_{\pi_k})$ generated given an initial $\pi_0$ such that the policy improvement condition

$$ \pi_{k+1}Q_{\pi_k} \le \pi_kQ_{\pi_k} \qquad (16) $$
holds. We consider the following identity

$$ (I - \gamma P\pi_{k+1})^{-1} - (I - \gamma P\pi_k)^{-1} = \gamma (I - \gamma P\pi_{k+1})^{-1} (P\pi_{k+1} - P\pi_k) (I - \gamma P\pi_k)^{-1}. \quad (17) $$

*Proof.* This identity can easily be checked, for instance by rearranging the right-hand side:

$$ \begin{aligned} \text{RHS} &= \gamma (I - \gamma P\pi_{k+1})^{-1} (P\pi_{k+1} - P\pi_k) (I - \gamma P\pi_k)^{-1}, \\ &= (I - \gamma P\pi_{k+1})^{-1} (\gamma P\pi_{k+1} - \gamma P\pi_k + I - I) (I - \gamma P\pi_k)^{-1}, \\ &= (I - \gamma P\pi_{k+1})^{-1} (- (I - \gamma P\pi_{k+1}) + (I - \gamma P\pi_k)) (I - \gamma P\pi_k)^{-1}, \\ &= -(I - \gamma P\pi_k)^{-1} + (I - \gamma P\pi_{k+1})^{-1}. \end{aligned} $$

□

Having $Q_{\pi_k} = (I - \gamma P\pi_k)^{-1} c$ from the Bellman equation, and with the identity (17) defined above applied to $c$, we get

$$ Q_{\pi_{k+1}} - Q_{\pi_k} = \gamma (I - \gamma P\pi_{k+1})^{-1} (P\pi_{k+1}Q_{\pi_k} - P\pi_kQ_{\pi_k}). \quad (18) $$

Furthermore, the operator $\gamma P\pi_{k+1}$ is non-negative, i.e., if $\psi$ is a nonnegative function,

$$ \forall (q, u) \in \mathcal{X} \times \mathcal{A}, \quad P\pi\psi(q, u) = \int_{\mathcal{X} \times \mathcal{A}} P(q, u, dx)\pi(x, da)\psi(x, a) $$

is nonnegative as well. Hence the operator $(I - \gamma P\pi_{k+1})^{-1} = \sum_{n=0}^{\infty} \gamma^n (P\pi_{k+1})^n$ also is nonnegative, so that $(Q_{\pi_k})_{k \in \mathbb{N}}$ is nonincreasing.

Thus, the construction above illustrates the principle of exact policy iteration, in which, the monotonicity property (16) is sufficient to ensure the convergence of this policy iteration scheme to optimality.

In practice, we would have to evaluate the policy at hand, i.e. compute $Q_{\pi_k}$, to improve on it. This is generally infeasible, since it requires solving the Bellman equation exactly, for instance by iterating indefinitely over $T_{\pi_k}$. To address this limitation, one resorts to Approximate Policy Iteration (API) algorithms. These methods replace the exact policy evaluation step with an approximation, obtained in our case by truncating the number of iterations of $T_{\pi_k}$, and by employing function approximation schemes. The resulting algorithms trade off exactness for tractability.

Suppose that we have initial data $(\pi_0, Q^0)$ and given $(\pi_k, Q^k)$, we improve the policy on an approximation $Q^k$ of $Q_{\pi_k}$ seeking $\pi_{k+1}$ such that:

$$ \pi_{k+1} Q^k \le \pi_k Q^k, \quad (19) $$

then update by approximately evaluating the new policy

$$ Q^{k+1} = T_{\pi_{k+1}}^{m} Q^{k}, \quad (20) $$
for a fixed number of iterations $m \in \mathbb{N}^*$, we get a sequence $(\pi_k, Q^k)_{k \in \mathbb{N}}$. In the literature, this index is known as the **optimistic order** [see 45, Section 2.5], and the resulting modified policy iteration is also referred to as **optimistic policy iteration**. In the following, however, we will simply refer to the index $m$ as multistep lookahead.

Neural networks provide a practical means of approximately carrying out both procedures. However, the monotonicity condition (19) in this case is not sufficient to guarantee that $Q_{\pi_{k+1}} \le Q_{\pi_k}$, but we can express

$$ 
\begin{align} 
Q_{\pi_{k+1}} - Q_{\pi_k} &= (I - \gamma P \pi_{k+1})^{-1} \gamma P (\pi_{k+1} Q_{\pi_k} - \pi_k Q_{\pi_k}), \tag{21} \\
&= (I - \gamma P \pi_{k+1})^{-1} \gamma P (\pi_{k+1} Q_{\pi_k} - \pi_{k+1} Q^k + \pi_{k+1} Q^k \nonumber \\
&\qquad - \pi_k Q^k + \pi_k Q^k - \pi_k Q_{\pi_k}), \nonumber \\
&= (I - \gamma P \pi_{k+1})^{-1} \gamma P \left[ (\pi_{k+1} - \pi_k)(Q_{\pi_k} - Q^k) + (\pi_{k+1} Q^k - \pi_k Q^k) \right]. \tag{22}
\end{align} 
$$

Thus, by positivity of $\gamma(I - \gamma P \pi_{k+1})^{-1}P$, the left-hand side is non-positive if and only if

$$ (\pi_{k+1} Q^k - \pi_k Q^k) \le -(\pi_{k+1} - \pi_k)(Q_{\pi_k} - Q^k). \quad (23) $$

The policy improvement step (19) in the algorithm is then effective only when the improvement compensates for the error in policy evaluation, this happens for instance when:

$$ \pi_{k+1} Q^k - \pi_k Q^k \le -2 \|Q_{\pi_k} - Q^k\|_{\infty}. \quad (24) $$

Thus to make policy improvement feasible, the policy evaluations need to be accurate enough. Now consider an arbitrary sequence of policies $\pi_k$ to assess the influence of policy variations and the role of the lookahead order on the errors of policy evaluation. By induction, we can write from identity (20)

$$ 
\begin{align} 
T_{\pi_{k+1}}^m Q &= c + \gamma P \pi_{k+1} (T_1^{m-1} Q), \nonumber \\
&= \sum_{i=0}^{m-1} \gamma^i (P \pi_{k+1})^i c + \gamma^m (P \pi_{k+1})^m Q. \tag{25}
\end{align} 
$$

It follows that $Q_{\pi_{k+1}} - Q^{k+1} = \gamma^m (P \pi_{k+1})^m (Q_{\pi_{k+1}} - Q^k)$. Then, plugging identity (22)

$$ 
\begin{align} 
Q_{\pi_{k+1}} - Q^{k+1} &= \gamma^m (P \pi_{k+1})^m [(Q_{\pi_{k+1}} - Q_{\pi_k}) + (Q_{\pi_k} - Q^k)], \tag{26} \\
&= \gamma^m (P \pi_{k+1})^m \Biggl[ (I - \gamma P \pi_{k+1})^{-1} (\gamma P \pi_{k+1} - \gamma P \pi_k)(Q_{\pi_k} - Q^k) \nonumber \\
&\qquad + (I - \gamma P \pi_{k+1})^{-1} (\gamma P \pi_{k+1} Q^k - \gamma P \pi_k Q^k) + (Q_{\pi_k} - Q^k) \Biggr], \nonumber \\
&= \gamma^m (P \pi_{k+1})^m \Biggl[ [I + (I - \gamma P \pi_{k+1})^{-1} (\gamma P \pi_{k+1} - \gamma P \pi_k)](Q_{\pi_k} - Q^k) \nonumber \\
&\qquad + (I - \gamma P \pi_{k+1})^{-1} \gamma P (\pi_{k+1} Q^k - \pi_k Q^k) \Biggr], \nonumber \\
&= \gamma^m (P \pi_{k+1})^m (I - \gamma P \pi_{k+1})^{-1} \Biggl[ (I - \gamma P \pi_k)(Q_{\pi_k} - Q^k) \nonumber \\
&\qquad + \gamma P (\pi_{k+1} Q^k - \pi_k Q^k) \Biggr]. \tag{27}
\end{align} 
$$
**Algorithm 1** Actor-critic with lookahead value update

1: Estimate the delay $\tau$ (section 4.1)

2: Set the lookahead order $m$

3: Initialize parameters $\theta_0$ of $\pi_0$ and $\psi_0$ of $Q^0$ networks; an empty memory $\mathcal{M}$

4: **for** $k \leftarrow 0$ **to** number of episodes **do**

5: Initialize environment and observe initial output $y_0$;

6: **for** $t \leftarrow 0$ **to** time horizon $\Delta t$ **do**

7: Sample action $a_t \sim \pi_k(y_t, da)$

8: Apply action to the environment and get a transition $(y_n, u_n, y_{n+1}, c_{n+1}) = \text{Env}(n, u_n)$

9: Add $(y_n, u_n, y_{n+1}, c_{n+1})$ to $\mathcal{M}$

10: **end for**

11: Sample uniformly enough batches of transitions $\{\mathcal{B}_i\}_{0 \le i \le N} \subset \mathcal{M}$ to get satisfactory coverage of $\mathcal{M}$

12: **for** $i \leftarrow 0$ **to** N **do**

13: Update $Q^k$ by gradient descent in the direction of

$$ \frac{1}{|\mathcal{B}_i|} \partial_{\psi}|_{\psi=\psi_k} \sum_{\substack{(y_{-\tau}, u_0, y_{m-\tau}, [c_0:c_{m-1}]) \in \mathcal{B}_i \\ u' \sim \pi_k(y', du)}} \left| \sum_{p=0}^{m-1} \gamma^p c_p + \gamma^m Q^k(y_{m-\tau}, u_m) - Q(\psi; y_{-\tau}, u_0) \right| \quad (30) $$

14: **end for**

15: **for** $i \leftarrow 0$ **to** N **do**

16: Update $\pi_k$ by gradient descent in the direction of

$$ \partial_{\theta}|_{\theta=\theta_k} \frac{1}{|\mathcal{B}_i|} \sum_{\substack{(y, u, y', c) \in \mathcal{B}_i \\ u' \sim \pi_k(\theta; y', du)}} Q^k(y', u') $$

17: **end for**

18: **end for**

Taking the norm yields the estimate

$$ \|Q_{\pi_{k+1}} - Q^{k+1}\|_{\infty} \le \frac{\gamma^m}{1-\gamma} \left[ (1+\gamma) \|Q_{\pi_k} - Q^k\|_{\infty} + \gamma \delta_k \right], \quad (28) $$

with $\delta_k = \|\pi_{k+1}Q^k - \pi_k Q^k\|_{\infty}$ the policy improvement. Taking the limit superior we get

$$ \left(1 - \frac{\gamma^m(1+\gamma)}{1-\gamma}\right) \limsup \|Q_{\pi_{k+1}} - Q^{k+1}\|_{\infty} \le \frac{\gamma^{m+1}}{1-\gamma} \limsup \delta_k. \quad (29) $$

An easy calculation shows that when $m > \frac{\ln(1-\gamma)-\ln(1+\gamma)}{\ln \gamma}$ the inequality
$$ \limsup_{\|Q_{\pi_{k+1}} - Q^{k+1}\|_{\infty}} \le \frac{\gamma^{m+1}}{1 - \gamma - \gamma^m (1 + \gamma)} \limsup_{\delta_k} \delta_k, \quad (31) $$

holds and shows the asymptotic behavior of the error in fixed-point approximations. This process is not guaranteed to converge to an optimal policy in general, but it converges in practice to suboptimal policies. This is due to the propagation of errors due to function approximation for example or to the application of approximate Bellman operators, errors that are not addressed in these bounds, the interested reader is referred to [53].

### 4.2.1 Practical considerations

The proposed actor-critic algorithm with lookahead and delay estimation is outlined in the pseudocode provided in algorithm 1. An initial warm-up phase is dedicated to estimating the time delays characterizing the input-output behavior of the system. The main body of the algorithm adopts an actor-critic architecture in which a cost lookahead strategy is employed during the policy evaluation stage. In the next section, using a system governed by the Kuramoto-Sivashinsky (KS) equations, we analyze the advantages introduced by the preliminary delay estimation step, highlighting its positive influence on the accuracy of credit assignment within the learning process. Moreover, the control of the boundary layer developing over a flat plate is shown as an illustration of a case governed by Navier-Stokes equations. For both test cases, we use the same hyperparameters reported in table 1. This choice is motivated by the fine tuning of the KS equation to closely mimic the time delays characterizing the chosen Navier-Stokes test cases. In fact, since the dynamics are dominated by the presence of time delays, we opted not to change the hyperparameters of the system. Shallow networks are generally preferred over deeper ones; unlike tabular methods, policy iteration with neural networks does not guarantee monotonic improvement because approximation errors introduce both bias (due to limited network capacity) and variance (due to stochastic updates), undermining stability. In our derivation, we set aside the critic's approximation errors; however, because the actor depends on the critic's value estimates to update its policy, a poorly approximated critic can produce misleading gradients and thus lead to suboptimal or unstable learning. To account for policy improvement errors, we allow arbitrary policy sequences parameterized by neural networks. In practice, this aligns well with the standard trade-off between computational cost and expressivity: we restrict ourselves to relatively shallow networks, recognizing that while deeper models can represent more complex policies (reducing bias), they demand more samples and training time (increasing variance).

Regarding the ratio $N = \frac{\text{episode length}}{\text{batch size}}$, we choose $N$ such that every transition in memory is sampled, on average, once per epoch. In lines 12–15 of algorithm 1 we replace the customary single-step cost $c_t$ with an $m$-sample window $\{c_{t+p}\}_{p=0}^{m-1}$ (see $[c_0 : c_{m-1}]$ in (30)). Because each $c_{t+p}$ is already produced by the environment, assembling the window merely keeps a rolling buffer of $m$ scalars and performs operations that are negligible compared with the forward/backward passes of the neural networks. In practice, the wall-clock cost of forming the window is negligible for both
<table><caption>**Table 1:** Parameters used in the implementation of the numerical examples.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Actor hidden layers</td><td>2</td></tr><tr><td>Actor units per layer</td><td>[32, 64]</td></tr><tr><td>Critic hidden layers</td><td>2</td></tr><tr><td>Critic units per layer</td><td>[16, 32]</td></tr><tr><td>Actor non-linearity</td><td>None</td></tr><tr><td>Critic non-linearity</td><td>ReLU</td></tr><tr><td>Actor leaning rate</td><td>10<sup>-4</sup></td></tr><tr><td>Critic leaning rate</td><td>10<sup>-4</sup></td></tr><tr><td>Batch size</td><td>256</td></tr><tr><td>Multistep / lookahead order (m)</td><td>1–2000</td></tr><tr><td>Discount rate γ</td><td>0.99</td></tr><tr><td>Memory size</td><td>20000</td></tr><tr><td>Number of time-steps per episode</td><td>20000</td></tr></tbody></table>

the KS and Navier-Stokes test cases, while the smoothed return noticeably stabilises training and speeds up convergence.

# 5 Results

In this section, we present the results of the closed-loop systems based on the two setups described in section 2. The algorithm, which was introduced from a theoretical perspective in section 3 and section 4, is summarized as pseudocode in algorithm 1.

## 5.1 Impact of time delays on the training process

Accurate knowledge of time delays in input-output systems is crucial for designing effective control strategies, particularly in flow control settings, where actuators and sensors may be spatially separated and signals require finite time to propagate. When explicit models or exact transport velocities are unavailable or unreliable, a purely data-driven approach to estimating time delays becomes invaluable. Here, we demonstrate such an approach using the KS equation introduced in the previous section, though the methodology readily extends to more general configurations. We account for time delays, as discussed in earlier sections, by leveraging Takens' theorem and direct delay computations.

This estimation procedure is necessary, as illustrated in fig. 5(left), where we compare the average training cost of a model that includes a delay-estimation pretraining phase with the same test case lacking such estimation. In the absence of a dynamic model, neglecting sensor-to-actuator time delays can critically undermine the performance of any learned control policy. If the RL agent is not provided with accurate temporal information, the training phase yields no significant improvement in the cost function, effectively halting meaningful policy updates. In contrast, when the estimated delay is incorporated into the RL framework, the agent can properly align
<figure><img src="image_7.png" /><figcaption>**Fig. 5:** Control of the dynamics governed by the KS equation using a policy trained with an lookahead order $m = 2000$ (see fig. 6). This figure compares the effect of incorporating a delay estimation procedure on convergence. On the left panel, the average cost is shown as a function of the number of episodes: by explicitly estimating time delays during the pre-training step, convergence is achieved within approximately 10 episodes. On the right, the maximal policy variations across all data batches are shown. (see section 4).</figcaption></figure>

observations with actions in time. This temporal alignment enables effective control decisions, as evidenced by a steady cost decrease across ten distinct training runs (initialized with different random seeds); this cost reduction coincides with policy updates, which are theoretically bounded as shown in prior sections. These learned policies effectively damp perturbation energy, as shown later in section 5.3. Consequently, incorporating time-delay estimation is a prerequisite for reliable training and robust performance when the underlying system dynamics remain unknown.

## 5.2 Lookahead order and the credit assignment problem

In fig. 5, the test cases adopt a look-ahead order of $m = 2000$. In this section, we analyze the impact of this choice from a twofold perspective: convergence and the credit assignment problem. The latter is one of the main challenges in reinforcement learning (RL), involving the task of identifying which actions taken at earlier times contribute to changes in a later cost signal [see 15, chapters 1.3 and 6]. In fluid systems, where sensors, actuators, and the region of interest may be spatially separated, significant transport delays can exist. Tracing the causal chain from individual control inputs to fluctuations in performance then becomes notably complex. These credit assignment issues, combined with the time delays discussed in the previous section, can hinder learning or result in suboptimal policies if the algorithm is not equipped to handle delayed costs effectively. RL architectures may therefore struggle to converge or may converge only to suboptimal solutions.

We investigate this problem by comparing two configurations shown in fig. 6. On the left, we consider a plant with a cost sensor placed at $x_z = 405$, while on the right,
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** Evolution of the average cost per episode during the model's training on the KS equation, shown as a function of the lookahead order *m*. Solid lines represent the average over 10 different random seeds, with lightly shaded areas indicating the variance. The left panel shows results for a cost sensor located at *x*<sub>*z*</sub> = 405, while the right panel presents the same analysis for a different plant configuration, where the cost sensor is positioned at *x*<sub>*z*</sub> = 450. The right panel corresponds to a scenario with more pronounced credit assignment: the positive impact of a higher *m* is demonstrated by the improvement in average performance observed for *m* ≥ 2000 after 25–30 episodes. Lower values of *m* lead to slower convergence, whereas further increases beyond this threshold do not yield significant improvements. This effect is less pronounced in the left panel, where convergence occurs already at *m* = 500.</figcaption></figure>

the sensor is positioned further downstream at $x_z = 450$. In the latter case, more pro-
nounced credit assignment challenges arise due to the longer delay introduced by the
downstream placement of the cost sensor. However, this configuration is particularly
relevant, as it reflects potential geometric constraints that may arise from the coupling
of specific actuator types — such as plasma actuators [54] — with sensor placement.
Both configurations are tested using policies trained with varying $m$ to assess how the
look-ahead influences the average cost per episode during training on the KS equation.
For the second configuration, depicted on the right, a significant improvement only
emerges once $m$ reaches or exceeds 2000, illustrating how insufficient values of $m$ result
in slower or incomplete convergence. Notably, the drop observed for $m = 2000$ aligns
approximately with the convective time delay between a wavepacket triggered at the
actuator location $x_a = 400$ and its detection at the cost sensor. This interpretation
is further supported by the left figure, which represents the configuration less con-
strained by geometric limitations: placing the cost sensor closer, at $x_z = 405$, alleviates
the credit assignment issue, allowing convergence to occur with lower values of $m$.
These results highlight how sensor placement, time delays, and reinforcement learning
hyperparameters jointly determine the agent's ability to assign credit effectively over
extended temporal horizons.
As a further point, we note — based on the analysis in section 4.2 — that incorporating the look-ahead parameter *m* also helps in bounding the error during the policy improvement step and enhances convergence, all while using the same amount of data for training the transition dynamics. In this sense, the algorithm remains data-efficient compared to standard actor-critic algorithms, as it achieves better convergence with no increase in data. However, compared to typical applications [44], the higher number of *m* iterations required here appears primarily due to the difficulty of the credit assignment problem; as the experiments indicate, stronger or more persistent credit assignment challenges require higher values of *m* to ensure that the agent maintains the necessary temporal span to connect actions to their eventual outcomes. The identified policies are capable of effectively damping the perturbations, as demonstrated in the following two sections for both the KS and boundary-layer flow cases.

## 5.3 Control of the Kuramoto-Sivashinsky equation

The dynamics governed by the linearized KS equation are discussed in the following. In fig. 6, the convergence of the policy is analyzed as a function of the episodes by monitoring the maximal policy variations over all the data batches ($\hat{\delta}$) and the average costs, which is comparable to the value function defined with discounted costs (see [55]). The curves are obtained by averaging statistics collected by running 10 different random seeds. By perusing the behavior of the two curves, it is possible to observe that convergence is attained in about 10 episodes. The policies are then applied for the control of the dynamics. The results are illustrated in fig. 7 where on the left is shown the cost signal for the uncontrolled (black solid line) and the controlled (red solid line) cases. More than one order of magnitude of reduction is observed in the cost, coinciding with effective control of the dynamics of the system as depicted in the spatio-temporal contour plot shown in fig. 7(right).

Finally, the energy behavior is monitored for a quantitative analysis. In fig. 8, the results obtained from the uncontrolled simulations are shown as a black solid line, while in red are reported the results related to the controlled case. On the left, the energy averaged in time is shown as a function of the streamwise direction $x$; more than one order of magnitude reduction is shown after the action of the controller, for $x > x_a = 400$; on the right, the average of the energy in space is shown as a function of time $t$, using the same legend.

## 5.4 Control of the Tollmien-Schlichting waves

Finally, we discuss the control of the Tollmien-Schlichting waves, excited by a disturbance $d$ and evolving in a boundary layer developing on a flat plate at $Re_{\delta^*} = 1000$. The algorithm is tuned using the same hyper-parameters already considered for the KS case, see also table 1. We observed that for the boundary layer case as well the convergence is achieved in about 10 episodes. The resulting controlled dynamics is shown in fig. 9(lower inset) and compared with the open-loop dynamics (top inset). The actuation is introduced at $x_a = 400$, where it is possible to observe the control action. From a quantitative viewpoint, the control of the wavepackets corresponds to
<figure><img src="image_9.png" /><figcaption>**Fig. 7:** Control of the dynamics governed by the KS equation. In the left panel, the cost signal $z^2(t)$ is shown for both the uncontrolled (black solid line) and controlled (red solid line) cases. The spatio-temporal contour plot in the right panel illustrates the dynamics of the controlled system, where a significant reduction in amplitude is observed at $x_a = 400$, corresponding to the actuator location.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8:** Control of the dynamics governed by the KS equation. On the left panel, the time-averaged energy is shown as a function of the streamwise direction $x$, with the uncontrolled dynamics represented by a black solid line and the controlled dynamics by a red solid line. Using the same legend, the space-averaged energy as a function of time $t$ is displayed on the right panel.</figcaption></figure>

a reduction of both the cost signal and the energy averaged in space: these two quantities are shown as a function of time $t$ in fig. 10. Together with the KS test case, these results highlight the effectiveness of incorporating time-delay estimation into the training process, leading to a pronounced suppression of flow perturbations relative to the free uncontrolled dynamics.
<figure><img src="image_11.png" /><figcaption>**Fig. 9:** Control of the boundary layer flow. Two flow fields are compared: the uncontrolled dynamics (top) is characterized by the growth of Tollmien-Schlichting (TS) wave packets traveling from left to right, with perturbations continually excited by white noise. The controlled dynamics (bottom) shows the effect of actuation introduced at $x_a = 400$, where the incoming wave packet is suppressed by the control action. Both snapshots are taken at $t = 5000$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10:** Comparison between controlled (red) and uncontrolled (black) dynamics. Left: cost signal $z(t)$. Right: total perturbation kinetic energy.</figcaption></figure>

# 6 Conclusions

In this article, we considered the control of convective flows using a reinforcement learning strategy. To this end, two configurations were introduced as demonstrators: a linearized version of the Kuramoto-Sivashinsky equation and the two-dimensional
boundary-layer flow developing over a flat plate. An appropriate choice of flow param-
eters allows for comparable input–output behavior in the two systems once equipped
with sensors and actuators, serving respectively as outputs and inputs. While we con-
sidered simplified environments, the choice of a minimal input–output configuration
stems from practical considerations. Indeed, environments equipped with an unreal-
istically large number of sensors are potentially limited from an experimental and
practical standpoint. In that sense, from the input–output configuration perspective,
localized sensors and actuators provide a more realistic setup for practical experiments
than a dense sensor array, which is common in many reinforcement learning studies.
However, the combination of local sensing, geometry, and flow features induces large
time delays that influence control design decisions. In this limit, the relative positioning
of inputs and outputs along the streamwise direction can also significantly impact the
system dynamics. Building on analogies with classical control approaches such as feed-
forward control and adaptive filtering — successfully employed in shear flow control
— we adopt a similar setup in combination with an actor–critic algorithm. Compared
to current flow control practices using reinforcement learning, the proposed algorithm
features: i) improved convergence due to the introduction of the lookahead order *m*;
ii) efficient data usage with respect to vanilla actor–critic algorithms (*m* = 1); iii) an
application of Takens’s theorem integrated with time-delay identification during pre-
training. Theoretical considerations are thoroughly discussed regarding convergence
and error bounds.

The resulting performance aligns with literature, though notably, only $O(10)$ episodes are needed for effective training. It is emphasized that the critic update employs *m*-step costs, which require no additional data processing and thus do not affect the overall computational time. Interestingly, the lookahead cost rollout positively influences the credit assignment problem posed by time delays between actuators and cost sensors; we observe that the number of required iterations is directly related to this delay. To the best of our knowledge, this is the first time such a strategy has been applied in a fluid mechanics context with physical considerations on the interplay between multi-step strategies, credit assignment, and limitations from the plant due to time delays. Moreover, we observed that applications in the RL literature based on standard benchmarks typically employ a lookahead rollout order below $m \approx 100$, which is well below the $O(1000)$ suggested by the physical constraints of the system considered here. This discrepancy becomes particularly pronounced in scenarios with stronger credit assignment effects.

Future work will focus on enhancing data efficiency during training through the
integration of exploration strategies and on extending the theoretical convergence
bounds of the proposed actor-critic algorithm to partially observable Markov decision
processes. From a fluid mechanics perspective, the current strategy is not limited to
simplified environments and can be extended to more realistic settings, in conjunction
with the minimal input-output plant considered in this study. We hope the present
study can encourage the fluid mechanics community to further engage with reinforce-
ment learning approaches, without forgetting the lessons from standard flow control
such as the critical requirements addressed in robust control.
**Acknowledgements.** The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under the grant ANR-REASON (ANR-21-CE46-0008).

**Author contribution.** Amine Saibi: investigation, methodology, implementation and data curation, writing; Lionel Mathelin: conceptualization and writing; Onofrio Semeraro: conceptualization, funding acquisition and writing.

**Funding.** The authors declare that they have no financial interests.

**Conflict of interest/Competing interests.** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

**Ethics approval and consent to participate.** Not applicable.

**Data availability/Materials availability/Code availability.** Data and training models used for this article will be available upon reasonable request.

## References

[1] Gad-el-Hak, M.: Flow Control: Passive, Active, and Reactive Flow Management. Cambridge University Press, Cambridge, UK (2000)

[2] Kim, J., Bewley, T.R.: A linear systems approach to flow control. Annu. Rev. Fluid Mech. **39**, 383–417 (2007)

[3] Bagheri, S., Henningson, D.S., Hoepffner, J., Schmid, P.J.: Input-output analysis and control design applied to a linear model of spatially developing flows. App. Mech. Rev. **62**(2) (2009)

[4] Sipp, D., Schmid, P.J.: Linear closed-loop control of fluid instabilities and noise-induced perturbations: A review of approaches and tools. Appl. Mech. Rev. **68**(2), 020801 (2016)

[5] Brunton, S.L., Noack, B.R.: Closed-loop turbulence control: Progress and challenges. Appl. Mech. Rev. **67**(5), 050801 (2015)

[6] Cattafesta III, L.N., Sheplak, M.: Actuators for active flow control. Annu. Rev. Fluid Mech. **43**, 247–272 (2011)

[7] Kühnen, J., Song, B., Scarselli, D., Budanur, N.B., Riedl, M., Willis, A.P., Avila, M., Hof, N.: Destabilizing turbulence in pipe flow. Nat. Phys. **14**(4), 386–390 (2018)

[8] Rowley, C.W.: Model reduction for fluids, using balanced proper orthogonal decomposition. Int. J. Bifurcation Chaos **15**(03), 997–1013 (2005)
[9] Ma, Z., Ahuja, S., Rowley, C.W.: Reduced-order models for control of fluids using the eigensystem realization algorithm. Theor. Comp. Fluid Dyn. **25**, 233–247 (2011)

[10] Ljung, L.: System Identification. Wiley Online Library, Hoboken, NJ-US (1999)

[11] Overschee, P., De Moor, B.: Subspace Identification for Linear Systems: Theory, Implementation, Applications. Springer, New York, NY (1996).
https://link.springer.com/book/10.1007/978-1-4613-0465-4

[12] Hervé, A., Sipp, D., Schmid, P.J., Samuelides, M.: A physics-based approach to flow control using system identification. J. Fluid Mech. **702**, 26–58 (2012)

[13] Juillet, F., McKeon, B., Schmid, P.J.: Experimental control of natural perturbations in channel flow. J. Fluid Mech. **752**, 296–309 (2014)

[14] Bertsekas, D.P.: Dynamic Programming and Optimal Control. Athena Scientific, Boston, MA, USA (1995)

[15] Sutton, R.S., Barto, A.G.: Reinforcement Learning: An Introduction. MIT Press, Boston, MA-US (2018)

[16] Lewis, F.L., Vrabie, D., Syrmos, V.L.: Optimal Control. John Wiley & Sons, Hoboken, NJ-US (2012)

[17] Glad, T., Ljung, L.: Control Theory. Taylor & Francis, Oxford, UK (2000)

[18] Åström, K.J., Wittenmark, B.: Adaptive Control. Courier Corporation, CA-US (2008)

[19] Recht, B.: A tour of reinforcement learning: The view from continuous control. Annu. Rev. Control Robot. Auton. Syst. **2**(1), 253–279 (2019)

[20] Guéniat, F., Mathelin, L., Hussaini, M.Y.: A statistical learning strategy for closed-loop control of fluid flows. Theor. Comp. Fluid Dyn. **30**, 497–510 (2016)

[21] Rabault, J., Kuchta, M., Jensen, A., Réglade, U., Cerardi, N.: Artificial neural networks trained through deep reinforcement learning discover control strategies for active flow control. J. Fluid Mech. **865**, 281–302 (2019)

[22] Bucci, M.A., Semeraro, O., Allauzen, A., Wisniewski, G., Cordier, L., Mathelin, L.: Control of chaotic systems by deep reinforcement learning. Proc. R. Soc. A **475**(2231), 20190351 (2019)

[23] Fan, D., Yang, L., Wang, Z., Triantafyllou, M.S., Karniadakis, G.E.: Reinforce-ment learning for bluff body active flow control in experiments and simulations. PNAS **117**(42), 26091–26098 (2020)
[24] Guastoni, L., Rabault, J., Schlatter, P., Azizpour, H., Vinuesa, R.: Deep reinforcement learning for turbulent drag reduction in channel flows. Eur. Phys. J. E **46**(4), 27 (2023)

[25] Sonoda, T., Liu, Z., Itoh, T., Hasegawa, Y.: Reinforcement learning of control strategies for reducing skin friction drag in a fully developed turbulent channel flow. J. Fluid Mech. **960**, 30 (2023)

[26] Vignon, C., Rabault, J., Vinuesa, R.: Recent advances in applying deep reinforcement learning for flow control: Perspectives and future directions. Phys. Fluids **35**(3) (2023)

[27] Pino, F., Schena, L., Rabault, J., Mendez, M.A.: Comparative analysis of machine learning methods for active flow control. J. Fluid Mech. **958**, 39 (2023)

[28] Fabbiane, N., Semeraro, O., Bagheri, S., Henningson, D.S.: Adaptive and model-based control theory applied to convectively unstable flows. App. Mech. Rev. **66**(6), 060801 (2014)

[29] Xu, D., Zhang, M.: Reinforcement-learning-based control of convectively unstable flows. J. Fluid Mech. **954**, 37 (2023)

[30] Sturzebecher, D., Nitsche, W.: Active cancellation of Tollmien-Schlichting instabilities on a wing using multi-channel sensor actuator systems. Intl J. Heat and Fluid Flow **24**, 572-583 (2003)

[31] Semeraro, O., Bagheri, S., Brandt, L., Henningson, D.S.: Transition delay in a boundary layer flow using active control. J. Fluid Mech. **731**, 288-311 (2013)

[32] Fabbiane, N., Simon, B., Fischer, F., Grundmann, S., Bagheri, S., Henningson, D.S.: On the role of adaptivity for robust laminar flow control. J. Fluid Mech. **767** (2015)

[33] Rathnasingham, R., Breuer, K.S.: Active control of turbulent boundary layers. J. Fluid Mech. **495**, 209-233 (2003)

[34] Lundell, F.: Reactive control of transition induced by free-stream turbulence: an experimental demonstration. J. Fluid Mech. **585**, 41-71 (2007)

[35] Belson, B.A., Semeraro, O., Rowley, C.W., Henningson, D.S.: Feedback control of instabilities in the two-dimensional Blasius boundary layer: The role of sensors and actuators. Phys. Fluids **25**(5), 054106 (2013)

[36] Doyle, J.C., Glover, K., Khargonekar, P.P., Francis, B.A.: State-space solutions to standard $H_2$ and $H_\infty$ control problems. IEEE Trans. Autom. Control **34**, 831-847 (1989)

[37] Zhou, K., Doyle, J.C., Glover, K.: Robust and Optimal Control. Prentice Hall,
Upper Saddle River, NJ (1996)

[38] Joslin, R.D., Miller, D.N.: Fundamentals and Applications of Modern Flow Control. American Institute of Aeronautics and Astronautics, ??? (2009)

[39] Lee, T., Kim, J., Lee, C.: Turbulence control for drag reduction through deep reinforcement learning. Physical Review Fluids 8(2), 024604 (2023)

[40] Zhou, Z., Zhang, M., Zhu, X.: Reinforcement-learning-based control of turbulent channel flows at high reynolds numbers. Journal of Fluid Mechanics 1006, 12 (2025)

[41] Monroc, C.B., Bušić, A., Dubuc, D., Zhu, J.: Actor critic agents for wind farm control. In: 2023 American Control Conference (ACC), pp. 177-183 (2023). IEEE

[42] Silver, D., Lever, G., Heess, N., Dęgris, T., Wierstra, D., Riedmiller, M.: Deterministic Policy Gradient Algorithms. In: International Conference on Machine Learning (2014)

[43] Lillicrap, T.P., Hunt, J.J., Pritzel, A., Heess, N., Erez, T., Tassa, Y., Silver, D., Wierstra, D.: Continuous control with deep reinforcement learning. arXiv preprint arXiv:1509.02971 (2015)

[44] Scherrer, B., Ghavamzadeh, M., Gabillon, V., Lesner, B., Geist, M.: Approximate modified policy iteration and its application to the game of tetris. Journal of Machine Learning Research 16(49), 1629-1676 (2015)

[45] Bertsekas, D.P.: Abstract Dynamic Programming, 2nd edn. Athena Scientific, Belmont, MA (2022)

[46] Nordström, J., Nordin, N., Henningson, D.: The fringe region technique and the fourier method used in the direct numerical simulation of spatially evolving viscous flows. SIAM-SISC 20(4), 1365-1393 (1999)

[47] Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N.: The fenics project version 1.5. Archive of numerical software 3(100) (2015)

[48] Kundu, P.K., Cohen, I.M., Dowling, D.R., Capecelatro, J.: Fluid Mechanics. Elsevier, San Diego, CA (2024)

[49] Trélat, E.: Contrôle Optimal: Théorie & Applications. Vuibert, Paris (2005)

[50] Sauer, T., Yorke, J.A., Casdagli, M.: Embedology. J. Stat. Phys. 65:3-4 (1991)

[51] Stark, J., Broomhead, D.S., Davies, M.E., Huke, J.: Takens embedding theorems for forced and stochastic systems. Nonlinear Anal. Theory Methods Appl. 30(8), 5303-5314 (1997)
[52] Azaria, M., Hertz, D.: Time delay estimation by generalized cross correlation methods. IEEE Trans. Acoust., Speech, Signal. Process. **32**(2), 280–285 (1984)

[53] Farahmand, A., Szepesvári, C., Munos, R.: Error propagation for approximate policy and value iteration. In: Advances in Neural Information Processing Systems (2010)

[54] Kriegseis, J., Simon, B., Grundmann, S.: Towards in-flight applications? A review on dielectric barrier discharge-based boundary-layer control. App. Mech. Rev. **68**(2), 020802 (2016)

[55] Tsitsiklis, J.N., Benjamin, V.R.: On average versus discounted reward temporal-difference learning. Mach. Learn. **49**, 179–191 (2022)