12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

LEARNING TURBULENCE CONTROL STRATEGIES WITH DATA-DRIVEN
REDUCED-ORDER MODELS AND DEEP REINFORCEMENT LEARNING

Kevin Zeng
Department of Chemical and Biological Engineering
University of Wisconsin - Madison
1415 Engineering Dr, Madison, WI 53706
kzeng3@wisc.edu

Alec J. Linot
Department of Chemical and Biological Engineering
University of Wisconsin - Madison
1415 Engineering Dr, Madison, WI 53706
linot@wisc.edu

Michael D. Graham
Department of Chemical and Biological Engineering
University of Wisconsin - Madison
1415 Engineering Dr, Madison, WI 53706
mdgraham@wisc.edu

ABSTRACT

Data-Driven Reduced Order Modeling

Design of active control strategies for turbulent drag reduction is a challenging task due to the complex dynamics
and the difficulty in devising good control targets. Deep
reinforcement learning (RL), an emergent machine learning
method capable of learning complex control strategies forhigh-dimensional systems from data, is able to address broad
macroscopic control objectives, such as drag minimization,
making it promising for flow control application. However, the iterative RL process requires vast amounts of data
to be generated from interactions with the target system.
For high-dimensional and computationally demanding simulations, such as direct numerical simulations (DNS), this data
generation step can be prohibitively expensive.
We mitigate this challenge in a completely data-driven
fashion by combining data-driven reduced-order models(ROM) of the flow system with deep RL. We construct ourROMs by combining an undercomplete autoencoder with aneural ordinary differential equation (ODE)–both of which aretrained directly from data. The autoencoder compresses thehigh-dimensional state representation into a low-dimensional
representation while the neural ODE predicts the system dynamics as an ODE in this new representation. This ROM issubstituted in place of the true system during RL training toaccelerate the learning process. The ROM-based control strategy is then deployed to the flow system for control validation.
We first apply our method to the Kuramoto-Sivashinsky
equation (KSE), a 1D turbulence proxy that exhibits spatiotemporal chaos, equipped with actuators. We demonstratethat we are capable of learning a ROM of the actuated dynamics and with a control goal analogous to drag reduction, weshow that the ROM-based RL strategies perform well in the
KSE. We highlight that the RL agent discovers and stabilizes
a forced equilibrium solution. Next, we apply our method to a
DNS of Couette flow with control in the form of jets that modify the wall-normal velocity at the wall. We find a ROM that
captures key short-time behavior of the underlying system using drastically fewer dimensions, and we show RL forces thesystem to relaminarize far more frequently than the underlying
system.

The first step in our procedure is to learn a ROM surrogate of the system we aim to control. We assume we haveaccess to time series data [s(t0 ), s(t1 ), ..., s(tM )] of the systemand that the dynamics of the system lie (or approximately lie)on some low-dimensional manifold M embedded in Rd (i.e.
s ∈ M ⊂ Rd ); an assumption that has been affirmed for many
dissipative systems (Temam, 1989). We can find a parameterization h ∈ Rdh that is of a much smaller dimension than
the ambient space dh ≪ d (when the dimension is minimal wewrite dM in place of dh ). We then formulate the time-evolutionof h with an ODE. Once constructed, the ROM forecasts trajectories of h (which we can map to the ambient space s) anduses these trajectories to train the RL agent. Figure 1 (a)-(c)outlines the learning procedure for the ROM.
To find the ROM we first use an undercomplete autoencoder to learn the parameterization of the manifold h. An
autoencoder consists of an encoding function h = χ(s; θE )and a decoding function s = χ̌(h; θD ). Both of these functions are constructed by neural networks (NN) with weights
θ = [θE , θD ] that are trained simultaneously to minimize the2reconstruction loss L = 1/M ∑M
i ||s(ti ) − χ̌(χ(s(ti )))|| . We
compute the gradient of this loss with respect to the NN parameters (dL/dθ ) and use stochastic gradient descent to update these parameters.
A challenge with this approach is determining the necessary degrees of freedom to fully parameterize the manifold
on which the the data lies. Whitney’s Embedding Theorem
proves as long as the encoding maps to R2dM then h is homeomorphic to s (Whitney, 1936). However, dM is unknown apriori, so we must find a way to determine it. We estimate
this dimension by tracking the error of the autoencoder as afunction of dimension, where we have empirically observed
a dramatic improvement in performance once dh agrees withdM (Linot & Graham, 2020, 2021).
After finding a reduced-dimensional representation h, wenext find a dynamical system for h. We chose to learn an ODE
for the dynamics of h such that ḣ = g(h, a; θODE ), where g isa NN and a is the action taken by the controller. The action
must be input here because the RL agent iteratively interacts

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 1: Procedure for learning a NODE-ROM from data, combining it with RL to approximate a control policy, anddeploying the approximate policy back to the true system

Example 1: Kuramoto-Sivashinsky Equation

with the system to determine the best control policy. Using
this ODE we predict new states by integrating forward in time

We consider the periodic KSE, given by∂v∂ v ∂ 2v ∂ 4v
= −v − 2 − 4 + f (x,t)
∂t
∂x ∂x
∂x

Z ti +τh̃(ti + τ) = h(ti ) +tig(h(t), a; θODE )dt.

(1)

(2)

where f is the control term defined as in Bucci et al. (2019)and corresponds to 4 evenly spaced Gaussian jets,
We train g by minimizing the difference between the prediction (h̃(ti + τ)) and the known state (h(ti + τ)), in the manifold
2
coordinates, giving the loss J = ∑M
i ||h(ti + τ) − h̃(ti + τ)||(Chen et al., 2019).

4
ai (t)(x − Xi )2
f (x,t) = ∑ √
exp −
.
2σs2
i=1 2πσs(3)

We evolve trajectories for a domain size of L = 22, which
exhibits spatiotemporal chaos with a Lyapunov time of ∼22time units. Spatial discretization is performed with Fourier
collocation on a mesh of 64 evenly spaced points and in ourformulation the state vector u consists of the solution values
at the collocation points. The system is time evolved with athird-order semi-implicit Runge-Kutta scheme (Bucci et al.,2019). From this simulation we train autoencoders with size
500 sigmoid activated hidden layers (for encoder, decoder) andneural ODEs on 40,000 snapshots of data (st = u(t)) separated
0.25 time units apart experiencing random jet actuations. Figure 2 shows the performance of autoencoders as we vary thedimension dh . For this system, the error is sufficiently low at
dh = 12, so we use this autoencoder for the mappings χ andχ̌. Furthermore, we emphasize that we achieve orders of magnitude improvement in reconstruction performance comparedto PCA while using the same number of dimensions. In thisrepresentation of the manifold coordinates, we train a neural
ODE using the same dataset. In Figure 3 we show the shorttime predictions of this model with random actuations matchthe true system well for around 1-1.5 Lyapunov times. We
also comment that in the absence of actuation, i.e. at = 0, theROM performs similarly well, indicating that the ROM is able

Learning a Control Strategy From Data
We use deep RL to learn control strategies from data. The
general deep RL framework is a cyclic interactive learning process. During each cycle, the RL agent, the embodiment of thecontrol policy represented by neural-network at = π(st ; θRL ),outputs a control action, at , given a state observation of thesystem, st . The impact of this action on the system is thenquantified by a scalar reward, rt , which is defined by the control objective. During training, the agent attempts to learn the
mapping between st and at that maximizes the cumulative longtime reward and updates accordingly each cycle. In our ROMbased RL framework, we train the agent in the reduced spaceof the ROM, such that we instead learn at = π(ht ; θRL ). During
agent deployment to the original system for control validation,
we simply precede the agent with the encoder ht = χ(st ; θE ) tomap state observations, st , to ht as this is coordinate system theagent was trained with. Fig. 1 (d)-(e) outlines the RL training
and deployment procedure. We comment that our framework
can applied with any general deep RL method.

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
st+τ , where τ = 0.25. In this work the DDPG actor and critic
networks utilized ReLU activated hidden layers of size 128
and 64, respectively, followed by tanh and linear activations to
the outputs of size 4 and 1, respectively.
To assess the performance of our ROM-based policy, thelearned control agent was applied to the NODE-ROM, with an
example controlled trajectory shown in Fig. 4a. We note thatafter a brief control transient, the control agent navigates the
NODE-ROM to an equilibrium (steady) state and stabilizes it.
The quantities targeted for minimization, D and Pf , estimatedfrom the predicted trajectory u(t), are shown in Figure 4c, revealing that this equilibrium exhibits dissipation much lower
than the natural unactuated dynamics. To assess how wellthis ROM-based control policy transfers to the original KSE(i.e. the true system), the same policy is applied to the true
KSE with the same initial condition, as shown in Fig. 4b. We
note that the controlled trajectory in the KSE yields not only
quantitatively similar transient behavior but also the same lowdissipation equilibrium state as was targeted in the NODEROM. The transient behaviors between the two are structurally
very similar, although the NODE-ROM displays slightly lessstrongly damped oscillations as it drives the trajectory to thesteady state. The values of D and Pf computed from the true
system, shown in Fig. 4d, are nearly identical to that of the
NODE-ROM in Fig. 4c.
To demonstrate the robustness of the ROM-based policy,
shown in Fig. 4e are the dissipation trajectories of the true KSE
beginning from 15 randomly sampled test initial conditions
that the ROM-based control agent has not seen before. We
highlight that the control agent is able to consistently navigate
the system to the same low-dissipation state within ∼ 150 timeunits, with one initial condition requiring ∼ 200 time units toconverge. Finally we note that although the RL training horizons were only 100 time units long, the control agent is able
to generalize to achieve and maintain control well beyond thethe horizon it was trained in.
Here we emphasize that the ROM-based policy drives the
dynamics to an equilibrium state in both the NODE-ROM andthe true KSE, indicating that not only does the NODE-ROM
capture this state, but it captures the dynamics leading to it
accurately enough such that the RL agent could discover it
during training and exploit it in a manner that still translates
to the original system. We further emphasize that both the
NODE-ROM and agent were never explicitly informed of thislow-dissipation state’s existence. Finally, we highlight that thediscovered low-dissipation state is an unstable state that is stabilized by the control agent. If control is removed, the system
returns to the natural chaotic dynamics. The ROM-based RL
agent is comparable to that of an RL agent trained directly withoriginal system via direct interactions in both performance andtargeting strategy (Zeng & Graham, 2021).
These observations indicate that the RL policy trained on
the model transfers very well to the true system. We attribute
this performance to the fact that both the NODE-ROM and RL
agent operate in Markovian fashion, i.e. even if the model has
slight inaccuracies, so long as the modeled dynamics are reasonably accurate this does not matter once the agent makes
its new state observation. Returning to the dynamical significance of the low-dissipation equilibrium state discoveredand stabilized by the RL agent, a continuation in mean forcing magnitude was performed. To do so, we Newton-solved
for equilibrium solutions to the KSE starting with the discovered equilibrium state while gradually decreasing the magnitude of the mean actuation profile to zero, as was done in Zeng
& Graham (2021). Solutions identified by this continuation

Figure 2: Mean squared error of autoencoders and PCA

on test data of the KSE for various dimensions.

(a) random actuation sequences ai (t) (b)ground truth KSE trajectory starting from a random initial conditions following actuation sequences in (a), (c)the decoded NODE-ROM trajectory following actuation
sequence (a) and the same initial condition in (b).
Figure 3:

to capture the underlying natural dynamics absent of control
inputs despite having been trained with no control-free data.
We now use this ROM as a surrogate model for training the RL agent. In this demonstration we utilize the Deep
Deterministic Policy Gradient (DDPG) RL method (Lillicrap
et al., 2016). As an analogue to energy-saving in the flow
control problem, we seek to minimize the dissipation and total power input of the KSE. This optimization is realized by
2
maximizing the reward r = −D + Pf , where D = ⟨( ∂∂ xu2 )2 ⟩ and
Pf = ⟨( ∂∂ ux )2 ⟩ + ⟨u f ⟩, respectively. Here ⟨·⟩ is the spatial average.
The control agent was trained with 1000 episodes of 100
time units long (i.e. 400 transitions per episode), with each
episode beginning from a random on-attractor initial conditionof the natural, i.e. unforced, KSE. Jet actuations implemented
by the control agent, at , were maintained constantly from st to3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 4: ROM-based RL agent applied to the same initial condition in the a) true KSE and b) data-driven reduced-ordermodel (decoded, dh = 12). The corresponding invariant quantities of dissipation and total input power for the c) true KSE
and d) learned reduced-order model. The dashed black line represents the system average of the natural KSE dynamics.
e) Controlled dissipation trajectories of the true KSE beginning from 15 randomly sampled test initial conditions.
Example 2: Turbulent Couette Flow
The second system we consider is Plane Couette flow in achannel with two periodic boundary conditions defined by thenondimensionalized incompressible Navier-Stokes Equations.
∂u+ u · ∇u = −∇P + Re−1 ∇2 u∂t
∇ · u = 0,(4)

where u = [u, v, w] is the vector of streamwise, wall-normal,
and spanwise velocities, and P is the pressure. The boundary conditions at the wall are u(±1) = ±1, ∂ v/∂ y(±1) = 0,w(±1) = 0, v(1) = 0, and v(−1) = ai (t) f (x, z). This final
boundary condition is the actuation, located only at the bottom
wall, which is in the form of two spatially localized slot-stylejets that inject or suck fluid in the wall-normal direction withzero net-flux. These slot-jets are oriented along the streamwise
direction with Gaussian profiles f (x, z) in the spanwise direction with Vmax = max(| f (x, z)|) = 1/20 and ai (t) ∈ [−1, 1]. In
Fig. 6 we show the wall-normal velocity in the case of a0 = 1
and a1 = −1.
Our code was written in Python to speed up communication with the RL code and follows the same algorithms used by Gibson et al. (2008). In particular, weconvert the field to Fourier space in x and z and Chebyshev space in y with a resolution of 32x35x32 FourierChebyshev-Fourier modes. Then we use the multistep AdamsBashforth/Backward-Differentiation 3 method described in
Peyret (2011) for time evolution using a timestep of ∆t = 0.02.
This results in a set of implicit Helmholtz equations that wesolve at each timestep, and through an influence matrix wecompute the pressure required to satisfy ∂ v/∂ y(±1) = 0. This
procedure does not result in an incompressible velocity field,
so we apply a Tau correction to exactly enforce this constraint.
The influence matrix and Tau correction methods are describedin Kleiser & Schumann (1980). Unlike the Gibson et al. code,

Figure 5: Forcing continuation from the forced equi-

librium state (under forcing f = α22 ) discovered by
NODE-ROM based RL policy (orange, dashed) to theunactuated KSE system (purple). The known equilibrium E1 of the KSE system is also provided (dots).

in forcing magnitude are shown in Fig. 5, which reveals thatequilibrium state captured by the NODE-ROM and discoveredby the RL agent is connected to a known existing solution of
the KSE known as E1 (Cvitanović et al., 2010); we obtained
a similar result with an RL agent trained on interactions withthe full system (Zeng & Graham, 2021). A similar observation
was made for RL control of 2D bluff body flow (Li & Zhang,2022). We speculate that in systems with complex nonlinear
dynamics, the discovery and stabilization of desirable underlying equilibrium solutions (or other recurrent saddle-point solutions such as unstable periodic orbits) of the system may be afairly general feature of RL flow control approaches. The nonlinear and exploratory nature of RL algorithms facilitates thediscovery of such solutions, and since the dynamics are slow
near these solutions, little control action should be required tokeep trajectories near them.

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 6: Visualization of the localized wall-normal ve-

locity of the slot-style jets located at y = −1.

Figure 9: Trajectories beginning from 25 unseen initialconditions evolving under no control (blues) and random slot actuations (reds). The mean drag achieved
across the 25 trajectories is shown in black for no control (solid) and random actuation (dashed)

unit. We performed proper orthogonal decomposition (POD)
for dimension reduction. In future work an autoencoder will be
used. After reducing the dimension with POD, we train neural ODEs with this lower-order representation. Figure 8 shows
the trajectory of the 1st POD coefficient for the true systemand the ROM where h is constructed from the first 10 POD
modes. This preliminary result indicates we can achieve qualitatively similar trajectories with very few degrees of freedom
using neural ODEs.
With the success of training RL agents with ROMs of
the KSE, and having promising preliminary results with lowdimensional ROMs for Couette flow, our next objective isto train a RL agents on Couette flow ROMs utilizing theSoft Actor-Critic (SAC) RL algorithm (Haarnoja et al., 2018).
However, before this, we first learn a control policy directlyfrom the DNS (i.e. direct interactions) to determine a benchmark control performance and timescale required for training.
We train our RL agent for several hundred episodes, with each
episode lasting for 300 time units and beginning from a turbulent initial condition that does not naturally laminarize in
300 time units. Actions input by the control agent persist for 5
time units, i.e. the agent makes 60 actions per training episode.
The reward the agent attempts to maximize is the negative ofthe drag (normalized such that laminar is 0), which is given
R RL
by r = − 0Lx 0 z ∂ u/∂ y(−1) + ∂ u/∂ y(1) − 2dxdz. Shown in
Fig. 9 is the reward vs. time of trajectories generated by our
DNS evolving from 25 new unseen initial conditions under nocontrol and with random slot actuations, i.e. the slot actuations were randomly sampled from a uniform distribution of
the allowable control range every 5 time units. We observethat in the presence of random slot actuations, the mean dragacross the 25 trajectories over a horizon of 900 time units isworse than in the presence of no control. Furthermore, someinstances where the system would naturally laminarize are disrupted by the random slot-jets. Shown in Fig. 10 is again
the total drag vs. time of Couette flow generated by our DNS
evolving from the same 25 new unseen initial conditions under no control and in the presence of our trained RL agent. We
observe that with the learned RL control policy, the mean dragacross the 25 trajectories over a horizon of 900 time units isdrastically lowered compared to no control. We highlight that
the RL agent learns to direct the chaotic turbulent dynamics tothe laminar state within this window of 900 time units, which

Figure 7: The L2 norm of the velocity field, u, of a trajec-

tory produced by Channel Flow 2.0 by Gibson and our
python-based Couette DNS code starting from the same
initial condition.

Figure 8: Trajectory of the projection of the Couette flow

field onto the first POD mode, and the prediction by theneural ODE.

our code allows us to modify the boundary conditions at thewall by solving for intermediate values in the influence matrix and tau correction methods at each timestep, as opposed
to once at the beginning. As illustrated in figure 7, we validated our code by showing that solutions of both codes tracknearly exactly for 1000 time units.
We consider Couette flow at a Re = 400 and domain sizeof [Lx , Ly , Lz ] = [7π/4, 2, 6π/5], which are the parameters considered by Gibson et al. (2008). For this system we generated 80,000 snapshots of data (st = u(t)) separated by 1 time5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
by deep reinforcement learning. Proc. R. Soc. A. 475.
Chen, Ricky T. Q., Rubanova, Yulia, Bettencourt, Jesse & Duvenaud, David 2019 Neural ordinary differential equations.
arXiv preprint arXiv:1806.07366 .
Cvitanović, P., Davidchack, Ruslan L. & Siminos, Evangelos 2010 On the State Space Geometry of the Kuramoto–Sivashinsky Flow in a Periodic Domain. SIAM
Journal on Applied Dynamical Systems 9 (1), 1–33.
Gibson, J. F., Halcrow, J. & Cvitanović, P. 2008 Visualizing
the geometry of state space in plane Couette flow. Journalof Fluid Mechanics 611 (1987), 107–130.
Haarnoja, Tuomas, Zhou, Aurick, Abbeel, Pieter & Levine,
Sergey 2018 Soft actor-critic: Off-policy maximum entropy
deep reinforcement learning with a stochastic actor.
Kleiser, L. & Schumann, U. 1980 Treatment of Incompressibility and Boundary Conditions in 3-D Numerical Spectral
Simulations of Plane Channel Flows, pp. 165–173. Wiesbaden: Vieweg+Teubner Verlag.
Li, Jichao & Zhang, Mengqi 2022 Reinforcement-learningbased control of confined cylinder wakes with stability analyses. Journal of Fluid Mechanics 932, A44.
Lillicrap, Timothy P., Hunt, Jonathan J., Pritzel, Alexander,
Heess, Nicolas, Erez, Tom, Tassa, Yuval, Silver, David &Wierstra, Daan 2016 Continuous control with deep reinforcement learning.
Linot, Alec J. & Graham, Michael D. 2020 Deep learning
to discover and predict dynamics on an inertial manifold.
Phys. Rev. E 101, 062209.
Linot, Alec J. & Graham, Michael D. 2021 Data-driven
reduced-order modeling of spatiotemporal chaos withneural ordinary differential equations. arXiv preprint
arXiv:2109.00060 .
Peyret, Roger 2011 Spectral methods for incompressible viscous flow. Springer.
Temam, R. 1989 Do inertial manifolds apply to turbulence?
Physica D: Nonlinear Phenomena 37 (1-3), 146–152.
Whitney, Hassler 1936 Differentiable Manifolds. Annals of
Mathematics 37 (3), 645–680.
Zeng, Kevin & Graham, Michael D. 2021 Symmetry reduction
for deep reinforcement learning active control of chaotic
spatiotemporal dynamics. Phys. Rev. E 104, 014210.

Figure 10: Trajectories beginning from 25 unseen initialconditions evolving under no control (blues) and the RL
control agent (reds). The mean drag achieved across the
25 trajectories is shown in black for no control (solid)
and the RL agent (dashed).

can be seen by the 19 of 25 trajectories laminarizing due to thecontrol agent. This is a stark improvement compared to the 5of 25 naturally laminarizing trajectories in the presence of nocontrol and 4 of 25 in the presence of random slot actuations.
The current accomplishment is the first, to the extent ofour knowledge, application of deep RL to learn an active control strategy from a direct numerical simulation of turbulent
Couette flow. To train such an agent required 3+ weeks ofsimulation time on 4 processors. We highlight that our preliminary NODE-ROMs is capable of simulating 500 time units of
Couette flow in less than the time it takes to simulate 1 time
unit by DNS. In forthcoming work we aim to demonstrate theapplication of training the RL agent via NODE-ROM in placeof the computationally intensive DNS for efficient estimationof control strategies.

REFERENCES
Bucci, M. A., Semeraro, O., Allauzen, A., Wisniewski, G.,
Cordier, L. & Mathelin, L. 2019 Control of chaotic systems

6