3.1. Optimization via BO and LIPO

We assume that the policy is a predefined parametric function $\boldsymbol {a}={\rm \pi} (\boldsymbol {s}_t;\boldsymbol {w}^{{\rm \pi} })\in \mathbb {R}^{n_a}$ with a small number of parameters (say $n_w \sim {O}(10)$). The dimensionality of the problem enables efficient optimizers such as BO and LIPO; other methods are illustrated by Duriez et al. (Reference Duriez, Brunton and Noack2017).

3.1.1. Bayesian optimization

The classic BO uses a GPr as surrogate model of the function that must be optimized. In the ‘on-policy’ approach implemented in this work, this is the cumulative reward function $R(\boldsymbol {w})$; from (2.3) and (2.4), this is $R(\boldsymbol {w})=V^{{\rm \pi} }(\boldsymbol {s}_0)=Q(\boldsymbol {s_0},\boldsymbol {a}^{{\rm \pi} }_0)$.

Let $\boldsymbol {W}^*:=\{\boldsymbol {w}_1,\boldsymbol {w}_2\ldots \boldsymbol {w}_{n_*}\}$ be a set of $n_*$ tested weights and $\boldsymbol {R}^*:=\{R_1,R_2\ldots {R}_{n_*}\}$ the associated cumulative rewards. The GPr offers a probabilistic model that computes the probability of a certain reward given the observations $(\boldsymbol {W}^*,\boldsymbol {R}^*)$, i.e. $p(R(\boldsymbol {w})|\boldsymbol {W}^*,\boldsymbol {R}^{*})$. In a GPr this is

(3.1)\begin{equation} p(R(\boldsymbol{w})|\boldsymbol{R}^*,\boldsymbol{W}^*)= \mathcal{N}(\boldsymbol{\mu},\boldsymbol{\varSigma}), \end{equation}

where $\mathcal {N}$ denotes a multivariate Gaussian distribution with mean $\boldsymbol {\mu }$ and covariance matrix $\boldsymbol {\varSigma }$. In a Bayesian framework, (3.1) is interpreted as a posterior distribution, conditioned to the observations ($\boldsymbol {W}^*,\boldsymbol {R}^*$). A GPr is a distribution over functions whose smoothness is defined by the covariance function, computed using a kernel function. Given a set of data $(\boldsymbol {W}^*,\boldsymbol {R}^*)$, this allows for building a continuous function to estimate both the reward of a possible candidate and the uncertainties associated with it.

We are interested in evaluating (3.1) on a set of $n_E$ new samples $\boldsymbol {W}:=\{\boldsymbol {w}_1,\boldsymbol {w}_2\ldots \boldsymbol {w}_{n_E}\}$ and we denote as $\boldsymbol {R}:=\{R_1,R_2\ldots R_{n_E}\}$ the possible outcomes (treated as random variables). Assuming that the possible candidate solutions belong to the same GPr (usually assumed to have zero mean (Rasmussen & Williams Reference Rasmussen and Williams2005)) as the observed data $(\boldsymbol {W}^*,\boldsymbol {R}^*)$, we have

(3.2)\begin{equation} \begin{pmatrix}\boldsymbol{R}^* \\ \boldsymbol{R}\end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{0}, \begin{pmatrix}\boldsymbol{K}_{**} & \boldsymbol{K}_* \\ \boldsymbol{K}_*^T & \boldsymbol{K}\end{pmatrix} \right), \end{equation}

where $\boldsymbol {K}_{**}=\kappa (\boldsymbol {W}^*,\boldsymbol {W}^*) \in \mathbb {R}^{n_*\times n_*}$, $\boldsymbol {K}_{*}=\kappa (\boldsymbol {W},\boldsymbol {W}^*)\in \mathbb {R}^{n_E\times n_*}$, $\boldsymbol {K}=\kappa (\boldsymbol {W},\boldsymbol {W})\in \mathbb {R}^{n_E\times n_E}$ and $\kappa$ a kernel function.

The prediction in (3.1) can be built using standard rules for conditioning multivariate Gaussian, and the functions $\boldsymbol {\mu }$ and $\boldsymbol {\varSigma }$ in (3.1) become a vector $\boldsymbol {\mu _*}$ and a matrix $\boldsymbol {\varSigma _*}$,

(3.3)$$\begin{gather} \boldsymbol{\mu_*} = \boldsymbol{K}_*^T \boldsymbol{K}_R^{{-}1} \boldsymbol{R}^*\quad\in \mathbb{R}^{n_E}, \end{gather}$$

(3.4)$$\begin{gather}\boldsymbol{\varSigma_*} = \boldsymbol{K} - \boldsymbol{K}_*^T \boldsymbol{K}_R^{{-}1} \boldsymbol{K}_*\quad\in \mathbb{R}^{n_E\times n_E}, \end{gather}$$

where $\boldsymbol {K}_R=\boldsymbol {K}_{**}+\sigma _{R}^2\boldsymbol {I}$, with $\sigma _{R}^2$ the expected variance in the sampled data and $\boldsymbol {I}$ the identity matrix of appropriate size. The main advantage of BO is that the function approximation is sequential, and new predictions improve the approximation of the reward function (i.e. the surrogate model) episode after episode. This makes the GPr-based BO one of the most popular black-box optimization methods for expensive cost functions.

The BO combines the GPr model with a function suggesting where to sample next. Many variants exist (Frazier Reference Frazier2018), each providing their exploration/exploitation balance. The exploration seeks to sample in regions of large uncertainty, while exploitation seeks to sample at the best location according to the current function approximation. The most classic function, used in this study, is the expected improvement, defined as (Rasmussen & Williams Reference Rasmussen and Williams2005)

(3.5)\begin{equation} \operatorname{EI}(\boldsymbol{w}) = \begin{cases} (\varDelta - \xi)\varPhi(Z) + \sigma(\boldsymbol{w})\phi(Z) & \text{if}\ \sigma(\boldsymbol{w}) > 0,\\ 0 & \text{if}\ \sigma(\boldsymbol{w}) = 0, \end{cases} \end{equation}

with $\varDelta =\mu (\boldsymbol {w}) - R(\boldsymbol {w}^+)$ and $\boldsymbol {w}^+=\operatorname {arg\,max}_{\boldsymbol {w}} \tilde {R}(\boldsymbol {w})$ the best sample so far, where $\varPhi (Z)$ is the cumulative distribution function, $\phi (Z)$ is the probability density function of a standard Gaussian and

(3.6)\begin{equation} Z = \begin{cases} \dfrac{\varDelta- \xi}{\sigma(\boldsymbol{w})} & \text{if}\ \sigma(\boldsymbol{w}) > 0 ,\\ 0 & \text{if}\ \sigma(\boldsymbol{w}) = 0. \end{cases} \end{equation}

Equation (3.5) balances the desire to sample in regions where $\mu (\boldsymbol {w})$ is larger than $R(\boldsymbol {w}^+)$ (hence, large and positive $\varDelta$) versus sampling in regions where $\sigma (\boldsymbol {w})$ is large. The parameter $\xi$ sets a threshold over the minimal expected improvement that justifies the exploration.

Finally, the method requires the definition of the kernel function and its hyperparameters, as well as an estimate of $\sigma _y$. In this work the GPr-based BO was implemented using the Python API scikit-optimize (Head et al. Reference Head, Kumar, Nahrstaedt, Louppe and Shcherbatyi2020). The selected kernel function was a Mater kernel with $\nu =5/2$ (see Chapter 4 from Rasmussen & Williams Reference Rasmussen and Williams2005), which reads

(3.7)\begin{equation} \kappa (\boldsymbol{x},\boldsymbol{x}')=\kappa (r)=1+\frac{\sqrt{5} r}{l}+\frac{5 r^2}{3 l^2}\exp-\frac{\sqrt{5} r}{l}, \end{equation}

where $r=\|\boldsymbol {x}-\boldsymbol {x}'\|_2$ and $l$ is the length scale of the process. We report a detailed description of the pseudocode we used in Appendix A.1.

3.1.2. Lipschitz global optimization

Like BO, LIPO relies on a surrogate model to select the next sampling points (Malherbe & Vayatis Reference Malherbe and Vayatis2017). However, LIPO's surrogate function is the much simpler upper bound approximation $U(\boldsymbol {w})$ of the cost function $R(\boldsymbol {w})$ (Ahmed et al. Reference Ahmed, Vaswani and Schmidt2020). In the dlib implementation by King (Reference King2009), used in this work, this is given by

(3.8)\begin{equation} U(\boldsymbol{w})= \min_{i=1\ldots t}\left(R(\boldsymbol{w}_{i}) + \sqrt{\sigma_i + (\boldsymbol{w}-\boldsymbol{w}_i)^TK(\boldsymbol{w}-\boldsymbol{w}_i)} \right), \end{equation}

where $\boldsymbol {w}_i$ are the sampled parameters, $\sigma _i$ are coefficients that account for discontinuities and stochasticity in the objective function, and $K$ is a diagonal matrix that contains the Lipschitz constants $k_i$ for the different dimensions of the input vector. We recall that a function $R(\boldsymbol {w}):\mathcal {W}\subseteq \mathbb {R}^{n_w}\rightarrow \mathbb {R}$ is a Lipschitz function if there exists a constant $C$ such that

(3.9)\begin{equation} \lVert R(\boldsymbol{w}_1) - R(\boldsymbol{w}_2)\rVert \leq C\lVert\boldsymbol{w}_1 - \boldsymbol{w}_2\rVert, \quad \forall\ \boldsymbol{w}_1,\boldsymbol{w}_2\in\mathcal{W}, \end{equation}

where $\lVert {\cdot }\rVert$ is the Euclidean norm on $\mathbb {R}^{n_w}$. The Lipshitz constant $k$ of $R(\boldsymbol {w})$ is the smallest $C$ that satisfies the above condition (Davidson & Donsig Reference Davidson and Donsig2009). In other terms, this is an estimate of the largest possible slope of the function $R(\boldsymbol {w})$. The values of $K$ and $\boldsymbol {\sigma }_i$ are found by solving the optimization problem

(3.10)\begin{equation} \left. \begin{aligned} \min_{K,\sigma} & \quad \lVert K\rVert_F^{2} + 10^6\sum_{i=1}^{t}\,\sigma_i^2, \\ \text{s.t.} & \quad U(\boldsymbol{w}_i)\geq R(\boldsymbol{w}_i),\quad \forall i\in[1\cdots t], \\ & \quad \sigma_i \geq 0,\quad \forall i\in[1\cdots t], \\ & \quad K_{i,j} \geq 0,\quad \forall i,j\in[1\cdots d], \\ & \quad \text{K} = \{k_1,k_2,\ldots,k_{n_w}\}, \end{aligned} \right\} \end{equation}

where $10^6$ is a penalty factor and $\lVert\ {\cdot }\ \rVert _F$ is the Frobenius norm.

To compensate for the poor convergence of LIPO in the area around local optima, the algorithm alternates between a global and a local search. If the iteration number is even, it selects the new weights by means of the maximum upper bounding position (MaxLIPO),

(3.11)\begin{equation} \boldsymbol{w}_{k+1} = \mathop {{\rm arg}\ {\rm max}}\limits_{\boldsymbol{w}}(U(\boldsymbol{w})), \end{equation}

otherwise, it relies on a trust region (TR) method (Powell Reference Powell2006) based on a quadratic approximation of $R(\boldsymbol {w})$ around the best weights obtained so far $\boldsymbol {w}^*$, i.e.

(3.12) \begin{equation} \left. \begin{gathered} \boldsymbol{w}_{k+1} =\arg\max_{\boldsymbol{w}} \overbrace{(\boldsymbol{w}^* + g(\boldsymbol{w}^*)^{T}\boldsymbol{w} + \tfrac{1}{2}\boldsymbol{w}^T\boldsymbol{H} (\boldsymbol{w}^*)\boldsymbol{w})}^{m(\boldsymbol{w};\boldsymbol{w}^*)}, \\ \text{s.t.} \|\boldsymbol{w}_{k+1}\|< d(\boldsymbol{w}^*), \end{gathered} \right\} \end{equation}

where $g(\boldsymbol {w}^*)$ is the approximation of the gradient at $\boldsymbol {w}^*$ ($g(\boldsymbol {w}^*)\approx \boldsymbol {\nabla } R(\boldsymbol {w}^*))$, $\boldsymbol {H}(\boldsymbol {w}^*)$ is the approximation of the Hessian matrix $(\boldsymbol {H}(\boldsymbol {w}^*))_{ij} \approx {\partial ^2 R(\boldsymbol {w}^*)}/{\partial \boldsymbol {w}_i\partial \boldsymbol {w}_j}$ and $d(\boldsymbol {w}^*)$ is the radius of the trust region. If the TR-method converges to a local optimum with an accuracy smaller than $\epsilon$,

(3.13)\begin{equation} |R(\boldsymbol{w}_{k})-R(\boldsymbol{w}^*)|<\varepsilon, \quad \forall \boldsymbol{w}_{k}, \end{equation}

the optimization goes on with the global search method until it finds a better optimum. A detailed description of the pseudocode we used can be found in Appendix A.2.

3.2. Genetic programming

In the GP approach to optimal control, the policy $\boldsymbol {a}={\rm \pi} (\boldsymbol {s};\boldsymbol {w})$ is encoded in the form of a syntax tree. The parameters are lists of numbers and functions that can include arithmetic operations, mathematical functions, Boolean operations, conditional operations or iterative operations. An example of a syntax tree representation of a function is shown in figure 2. A tree (or program in GP terminology) is composed of a root that branches out into nodes (containing functions or operations) throughout various levels. The number of levels defines the depth of the tree, and the last nodes are called terminals or leaves. These contain the input variables or constants. Any combination of branches below the root is called a subtree and can generate a tree if the node becomes a root.

Figure 2.

Syntax tree representation of the function $2 x\sin (x)+\sin (x)+3$. This tree has a root ‘+’ and a depth of two. The nodes are denoted with orange circles while the last entries are leafs.

Syntax trees allow encoding complex functions by growing into large structures. The trees can adapt during the training: the user provides a primitive set, i.e. the pool of allowed functions, the maximum depth of the tree and set the parameters of the training algorithm. Then, the GP operates on a population of possible candidate solutions (individuals) and evolves it over various steps (generations) using genetic operations in the search for the optimal tree. Classic operations include elitism, replication, cross-over and mutations, as in genetic algorithm optimization (Haupt & Ellen Haupt Reference Haupt and Ellen Haupt2004). The implementation of GP in this work was carried out in the distributed evolutionary algorithms in Python (DEAP) (Fortin et al. Reference Fortin, De Rainville, Gardner, Parizeau and Gagné2012) framework. This is an open-source Python library allowing for the implementation of various evolutionary strategies.

We used a primitive set of four elementary operations ($+,-,/,\times$) and four functions ($\exp,\log,\sin,\cos$). In the second test case, as described in § 5.2, we also include an ephemeral random constant. The initial population of individuals varied between $n_I=10$ and $n_I=80$ candidates depending on the test case and the maximum depth tree was set to $17$. In all test cases the population was initialized using the ‘half-half’ approach, whereby half the population is initialized with the full method and the rest with the growth method. In the full method trees are generated with a predefined depth and then filled randomly with nodes and leafs. In the growth method trees are randomly filled from the roots: because nodes filled with variables or constant are terminals, this approach generates trees of variable depth.

Among the optimizers available in DEAP, in this work we used the $(\mu +\lambda )$ algorithm for the first two test cases and eaSimple (Banzhaf et al. Reference Banzhaf, Nordin, Keller and Francone1997; Vanneschi & Poli Reference Vanneschi and Poli2012; Kober & Peters Reference Kober and Peters2014; Bäck, Fogel & Michalewicz Reference Bäck, Fogel and Michalewicz2018) for the third one. These differ in how the population is updated at each iteration. In the $(\mu +\lambda )$ both the offsprings and parents participate to the tournament while in eaSimple no distinction is made between parents and offsprings and the population is entirely replaced at each iteration.

Details about the algorithmic implementation of this approach can be found in Appendix A.3.

3.3. Reinforcement learning via DDPG

The DDPG by Lillicrap et al. (Reference Lillicrap, Hunt, Pritzel, Heess, Erez, Tassa, Silver and Wierstra2015) is an off-policy actor–critic algorithm using an ANN to learn the policy (direct approach, figure 1a) and an ANN to learn the $Q$ function (indirect approach, figure 1b). In what follows, we call the $\varPi$ network the first (i.e. the actor) and the $Q$ network the second (i.e. the critic).

The DDPG combines the DPG by Silver et al. (Reference Silver, Lever, Heess, Degris, Wierstra and Riedmiller2014) and the DQN by Mnih et al. (Reference Mnih, Kavukcuoglu, Silver, Graves, Antonoglou, Wierstra and Riedmiller2013,  Reference Mnih2015). The algorithm has evolved into more complex versions such as the twin delayed DDPG (Fujimoto, van Hoof & Meger Reference Fujimoto, van Hoof and Meger2018), but in this work we focus on the basic implementation.

The policy encoded in the $\varPi$ network is deterministic and acts according to the set of weights and biases $\boldsymbol {w}^{{\rm \pi} }$, i.e. $\boldsymbol {a}={\rm \pi} (\boldsymbol {s}_t,\boldsymbol {w}^{{\rm \pi} })$. The environment is assumed to be stochastic and modelled as a MDP. Therefore, (2.3) must be modified to introduce an expectation operator,

(3.14)\begin{equation} Q^{\rm \pi}(\boldsymbol{s}_t,\boldsymbol{a}_t)=\mathbb{E}_{\boldsymbol{s}_{t},\boldsymbol{s}_{t+1}\sim E}\left[r(\boldsymbol{s}_t,\boldsymbol{a}_t)+\gamma Q^{\rm \pi}(\boldsymbol{s}_{t+1},\boldsymbol{a}^{\rm \pi}_{t+1})\right], \end{equation}

where the policy is intertwined in the action state relation, i.e. $Q^{{\rm \pi} }(\boldsymbol {s}_{t+1},\boldsymbol {a}_{t+1})=Q^{{\rm \pi} }(\boldsymbol {s}_{t+1},\boldsymbol {a}^{{\rm \pi} }(\boldsymbol {s}_{t+1}))$ and having used the shorthand notation $\boldsymbol {a}^{{\rm \pi} }_{t+1}={\rm \pi} (\boldsymbol {s}_{t+1},\boldsymbol {w}^{{\rm \pi} })$. Because the expectation operator in (3.14) solely depends on the environment ($E$ in the expectation operator), it is possible to decouple the problem of learning the policy ${\rm \pi}$ from the problem of learning the function $Q^{{\rm \pi} }(\boldsymbol {s}_t,\boldsymbol {a}_t)$. Concretely, let $Q(\boldsymbol {s}_t,\boldsymbol {a}_t;\boldsymbol {w}^{Q})$ denote the prediction of the $Q$ function by the $Q$ network, defined with weights and biases $\boldsymbol {w}^{Q}$ and let $\mathcal {T}$ denote a set of $N$ transitions $(\boldsymbol {s}_t,\boldsymbol {a}_t,\boldsymbol {s}_{t+1},r_{t+1})$ collected through (any) policy. The performances of the $Q$ network can be measured as

(3.15) \begin{equation} J^Q(\boldsymbol{w}^{Q})=\mathbb{E}_{\boldsymbol{s}_t,\boldsymbol{a}_t, \boldsymbol{r}_t\sim \mathcal{T}}[(Q(\boldsymbol{s}_t,\boldsymbol{a}_t; \boldsymbol{w}^{Q})-y_t)^2], \end{equation}

where the term in the squared brackets, called temporal difference (TD), is the difference between the old $Q$ value and the new one $y_t$, known as the TD target,

(3.16)\begin{equation} y_t=r(\boldsymbol{s}_t,\boldsymbol{a}_t)+\gamma Q (\boldsymbol{s}_{t+1},\boldsymbol{a}_{t+1};\boldsymbol{w}^{Q}). \end{equation}

Equation (3.15) measures how closely the prediction of the $Q$ network satisfies the discrete Bellman equation (2.3). The training of the $Q$ network can be carried out using standard stochastic gradient descent methods using the back-propagation algorithm (Kelley Reference Kelley1960) to evaluate $\partial _{\boldsymbol {w}^{Q}} J^Q$.

The training of the $Q$ network gives the off-policy flavour to the DDPG because it can be carried out with an exploratory policy that largely differs from the final policy. Nevertheless, because the training of the $Q$ network is notoriously unstable, Mnih et al. (Reference Mnih, Kavukcuoglu, Silver, Graves, Antonoglou, Wierstra and Riedmiller2013,  Reference Mnih2015) introduced the use of a replay buffer to leverage accumulated experience (previous transitions) and a target network to under-relax the update of the weights during the training. Both the computation of the cost function in (3.15) and its gradient are performed over a random batch of transitions $\mathcal {T}$ in the replay buffer $\mathcal {R}$.

The DDPG combines the $Q$ network prediction with a policy gradient approach to train the $\varPi$ network. This is inherited from the DPG by Silver et al. (Reference Silver, Lever, Heess, Degris, Wierstra and Riedmiller2014), who have shown that, given

(3.17)\begin{equation} J^{\rm \pi}(\boldsymbol{w}^{\rm \pi})=\mathbb{E}_{\boldsymbol{s}_t\sim E,\boldsymbol{a}_t\sim {\rm \pi}}\left[(r(\boldsymbol{s}_t,\boldsymbol{a}_t))\right] \end{equation}

is the expected return from the initial condition, the gradient with respect to the weights in the $\varPi$ network is

(3.18)\begin{equation} \partial_{\boldsymbol{w}^{\rm \pi}}J^{\rm \pi}=\mathbb{E}_{\boldsymbol{s}_t\sim E,\boldsymbol{a}_t\sim {\rm \pi}}[\partial_{\boldsymbol{a}} Q(\boldsymbol{s}_t,\boldsymbol{a}_t;\boldsymbol{w}^Q)\,\partial_{\boldsymbol{w}^{\rm \pi}} \boldsymbol{a}(\boldsymbol{s}_t;\boldsymbol{w}^{\rm \pi})]. \end{equation}

Both $\partial _{\boldsymbol {a}} Q(\boldsymbol {s}_t,\boldsymbol {a}_t; \boldsymbol {w}^Q)$ and $\partial _{\boldsymbol {w}^{{\rm \pi} }} \boldsymbol {a}(\boldsymbol {s}_t;\boldsymbol {w}^{{\rm \pi} })$ can be evaluated via back propagation on the $Q$ network and the $\varPi$ network, respectively. The main extension of DDPG over DPG is the use of DQN for the estimation of the $Q$ function.

In this work we implement the DDPG using Keras API in Python with three minor modifications to the original algorithm. The first is a clear separation between the exploration and exploitation phases. In particular, we introduce a number of exploratory episodes $n_{Ex}< n_{Ep}$ and the action is computed as

(3.19)\begin{equation} \boldsymbol{a}(\boldsymbol{s}_t)=\boldsymbol{a}(\boldsymbol{s}_t;\boldsymbol{w}^{\rm \pi})+\eta(\mbox{ep})\mathcal{E}(t;\theta,\sigma^2), \end{equation}

where $\mathcal {E}(t;\theta,\sigma )$ is an exploratory random process characterized by a mean $\theta$ and variance $\sigma ^2$. This could be the time-correlated (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930) noise or white noise, depending on the test case at hand (see § 4). The transition from exploration to exploitation is governed by the parameter $\eta$, which is taken as $\eta (\mbox {ep})=1$ if $\mbox {ep}< n_{Ex}$ where $d^{{ep}-n_{Ex}}$ if $\mbox {ep}>n_{Ex}$. This decaying term for $\mbox {ep}>n_{Ep}$ progressively reduces the exploration and the coefficient $d$ controls how rapidly this is done.

The second modification is in the selection of the transitions from the replay buffer $\mathcal {R}$ that are used to compute the gradient $\partial _{\boldsymbol {w}^{Q}} J^Q$. While the original implementation selects these randomly, we implement a simple version of the prioritized experience replay from Schaul et al. (Reference Schaul, Quan, Antonoglou and Silver2015). The idea is to prioritize, while sampling from the replay buffer, those transitions that led to the largest improvement in the network performances. These can be measured in terms of the TD error

(3.20)\begin{equation} \delta = r_t + \gamma Q(\boldsymbol{s}_{t+1}, \boldsymbol{a}^{\rm \pi}_{t+1};\boldsymbol{w}^{Q})- Q(\boldsymbol{s}_t, a_t;\boldsymbol{w}^{Q}). \end{equation}

This quantity measures how much a transition was unexpected. The rewards stored in the replay buffer ($r_{t}^{RB}$) and used in the TD computation are first scaled using a dynamic vector $r_{log}=[r^{RB}_1,r^{RB}_2,\ldots,r^{RB}_t]$ as

(3.21)\begin{equation} r_{t}^{RB} = \frac{r_t - \bar{r}_{log}}{{\rm std}(r_{log})+1\times 10^{{-}10}}, \end{equation}

where $\bar {r}_{log}$ is the mean value and $\textrm {std}(r_{log})$ is the standard deviation. The normalization makes the gradient steeper far from the mean of the sampled rewards, without changing its sign, and is found to speed-up the learning (see also van Hasselt et al. Reference van Hasselt, Guez, Hessel, Mnih and Silver2016).

As discussed by Schaul et al. (Reference Schaul, Quan, Antonoglou and Silver2015), it can be shown that prioritizing unexpected transitions leads to the steepest gradients $\partial _{\boldsymbol {w}^{Q}} J^Q$ and, thus, helps overcome local minima. The sampling is performed following a triangular distribution that assigns the highest probability $p(n)$ to the transition with the largest TD error $\delta$.

The third modification, extensively discussed in previous works on RL for flow control (Rabault & Kuhnle Reference Rabault and Kuhnle2019; Rabault et al. Reference Rabault, Ren, Zhang, Tang and Xu2020; Tang et al. Reference Tang, Rabault, Kuhnle, Wang and Wang2020), is the implementation of a sort of moving average of the actions. In other words, an action is performed for $K$ consecutive interactions with the environment, which in our work occur at every simulation's time step.

We illustrate the neural network architecture employed in this work in figure 3. The scheme in the figure shows how the $\varPi$ network and the $Q$ network are interconnected: intermediate layers map the current state and the action (output by the $\varPi$ network) to the core of the $Q$ network. For plotting purposes, the number of neurons in the figure is much smaller than the one actually used and indicated in the figure. The $\varPi$ network has two hidden layers with $128$ neurons each, while the input and output depends on the test cases considered (see § 4). Similarly, the $Q$ network has two hidden layers with $128$ neurons each and intermediate layers as shown in the figure. During the exploration phase, the presence of the stochastic term in the action selection decouples the two networks.

Figure 3.

The ANN architecture of the DDPG implementation analysed in this work. The illustrated architecture is the one used for the test case in § 4.3. During the exploration phase, the two networks are essentially decoupled by the presence of the stochastic term $\mathcal {E}$ that leads to exploration of the action space.

We detail the main steps of the implemented DDPG algorithm in Appendix A.4. It is important to note that, by construction, the weights in this algorithm are updated at each interaction with the system. Hence, $k=n$ and $N=1$ in the terminology of § 2. The notion of an episode remains relevant to control the transition between various phases of the learning process and to provide comparable metrics between the various algorithms.

4.1. A 0-D frequency cross-talk problem

The first selected test case is a system of nonlinear ordinary differential equations (ODEs) reproducing one of the main features of turbulent flows: the frequency cross-talk. This control problem was proposed and extensively analysed by Duriez et al. (Reference Duriez, Brunton and Noack2017). It essentially consists in stabilizing two coupled oscillators, described by a system of four ODEs, which describe the time evolution of four leading proper orthogonal decomposition modes of the flow past a cylinder. The model is known as the generalized mean field model (Dirk et al. Reference Dirk, Günther, Noack, King and Tadmor2009) and was used to describe the stabilizing effect of low frequency forcing on the wave flow past a bluff body (Pastoor et al. Reference Pastoor, Henning, Noack, King and Tadmor2008; Aleksic et al. Reference Aleksic, Luchtenburg, King, Noack and Pfeifer2010). The set of ODEs in the states $\boldsymbol {s}(t)=[s_1(t),s_2(t),s_3(t),s_4(t)]^\textrm {T}$, where ($s_1,s_2$) and ($s_3,s_4$) are the first and second oscillator, reads

(4.1)\begin{equation} \dot{\boldsymbol{s}} = \boldsymbol{F}(\boldsymbol{s})\boldsymbol{s} + \boldsymbol{A}\boldsymbol{a}, \end{equation}

where $\boldsymbol {a}$ is the forcing vector with a single scalar component interacting with the second oscillator (i.e. $\boldsymbol {a}=[0,0,0,a]^\textrm {T}$) and the matrix $\boldsymbol {F}(\boldsymbol {s})$ and $\boldsymbol {A}$ are given by

(4.2a,b) \begin{align} \boldsymbol{F}(\boldsymbol{s}) = \begin{bmatrix} \sigma(\boldsymbol{s}) & -1 & 0 & 0\\ 1 & \sigma(\boldsymbol{s}) & 0 & 0 & \\ 0 & 0 & -0.1 & -10 \\ 0 & 0 & 10 & -0.1\!\!\!\!\!\!\!\!\end{bmatrix},\quad \boldsymbol{A} = \left[\!\!\begin{array}{llll@{}l@{}} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1\end{array}\right]. \end{align}

The term $\sigma (\boldsymbol {s})$ models the coupling between the two oscillators:

(4.3)\begin{equation} \sigma(\boldsymbol{s}) = 0.1 - E_1 - E_2, \end{equation}

where $E_1$ and $E_2$ are the energy of the first and second oscillator given by

(4.4a,b)\begin{equation} E_1 = s_1^2+s_2^2 \quad E_2 = s_3^2 + s_4^2. \end{equation}

This nonlinear link is the essence of the frequency cross-talk and challenges linear control methods based on linearization of the dynamical system. To excite the second oscillator, the actuation must introduce energy to the second oscillator, as one can reveal from the associated energy equation. This is obtained by multiplying the last two equations of the system by $s_3$ and $s_4$, respectively, and summing them up to obtain

(4.5)\begin{equation} \tfrac{1}{2}\dot{E_2} ={-}0.2E_2 + s_4 u, \end{equation}

where $u\,s_4$ is the production term associated to the actuation.

The initial conditions are set to $\boldsymbol {s}(0)=[0.01,0,0,0]^\textrm {T}$. Without actuation, the system reaches a ‘slow’ limit cycle involving the first oscillator $(s_1,s_2)$, while the second vanishes ($(s_3,s_4)\rightarrow 0$). The evolution of the oscillator $(s_1,s_2)$ with no actuation is shown in figure 4(a); figure 4(b) shows the time evolution of $\sigma$, which vanishes as the system naturally reaches the limit cycle. Regardless of the state of the first oscillator, the second oscillator is essentially a linear second-order system with eigenvalues $\lambda _{1,2}=-0.1\pm 10\mathrm {i}$; hence, a natural frequency $\omega =10\,\textrm {rad}\,\textrm {s}^{-1}$.

Figure 4.

Evolution of the oscillator $(s_1,s_2)$ (a) of the variable $\sigma$ (4.3) (b) in the 0-D test case in the absence of actuation ($a=0$). As $\sigma \approx 0$, the system naturally evolves towards a ‘slow’ limit cycle.

The governing equations (4.1) were solved using scipy's package odeint with a time step of $\Delta \,t={\rm \pi} /50$. This time step is smaller than the one by Duriez et al. (Reference Duriez, Brunton and Noack2017) ($\Delta t={\rm \pi} /10$), as we observed this had an impact on the training performances (aliasing in LIPO and BO optimization).

The actuators’ goal is to bring to rest the first oscillator while exiting the second, leveraging on the nonlinear connection between the two and using the least possible actuation. In this respect, the optimal control law, similarly to Duriez et al. (Reference Duriez, Brunton and Noack2017), is the one that minimizes the cost function

(4.6)\begin{equation} \left.\begin{gathered} J = J_a + \gamma J_b = \overline{s_1^2 + s_2^2} + \alpha \overline{a^2}, \\ \text{where} \quad \overline{f(t)} = \frac{1}{40{\rm \pi}}\int_{20{\rm \pi}}^{60{\rm \pi}} f(t')\, {\rm d}t', \end{gathered} \right\} \end{equation}

where $\alpha$, set to $\alpha =10^{-2}$, is a coefficient set to penalize large actuations. Like the original problem in Duriez et al. (Reference Duriez, Brunton and Noack2017), the actions are clipped to the range $a_k\in [-1,1]$.

The time interval of an episode is set to $t\in [20{\rm \pi},60{\rm \pi} ]$; thus, much shorter than that used by Duriez et al. (Reference Duriez, Brunton and Noack2017). This duration was considered sufficient, as it allows the system to reach the limit cycle and to observe approximately $20$ periods of the slow oscillator. To reproduce the same cost function in a RL framework, we rewrite (4.6) as a cumulative reward, replacing the integral mean with the arithmetic average and setting

(4.7)\begin{equation} J = \frac{1}{n_t}\sum_{k=0}^{n_t-1} s^2_{1k}+s^2_{2k}+\alpha a^2_k={-} \sum_ {k=0}^{n_t-1}r_t={-}R, \end{equation}

with $r_t$ the environment's reward at each time step. For the BO and LIPO optimizers, the control law is defined as a quadratic form of the four system's states,

(4.8)\begin{equation} {\rm \pi}(\boldsymbol{s};\boldsymbol{w}) := \boldsymbol{g}_w^T\boldsymbol{s} + \boldsymbol{s}^T\boldsymbol{H}_w\boldsymbol{s}, \end{equation}

with $\boldsymbol {g}_w\in \mathbb {R}^4$ and $\boldsymbol {H}_w\in \mathbb {R}^{4x4}$. The weight vector associated to this policy is thus $\boldsymbol {w}\in \mathbb {R}^{20}$ and it collects all the entries in $\boldsymbol {g}_w$ and $\boldsymbol {H}_w$. For later reference, the labelling of the weights is as follows:

(4.9a,b)\begin{equation} \boldsymbol{g}_{w}= \begin{bmatrix}w_1\\w_2\\w_3\\w_4 \end{bmatrix} \quad \text{and}\quad \boldsymbol{H}_{w}= \begin{bmatrix}w_5 & w_9 & w_{13} & w_{17}\\w_6 & w_{10} & w_{14} & w_{18}\\w_7 & w_{11} & w_{15} & w_{19}\\w_8 & w_{12} & w_{16} & w_{20} \end{bmatrix}. \end{equation}

Both LIPO and BO seek for the optimal weights in the range $[-3,3]$. The BO was set up with a Matern kernel (see (3.7)) with a smoothness parameter $\nu =1.5$, a length scale of $l=0.01$, an acquisition function based on the expected improvement and an exploitation–exploration (see (3.5)) trade-off parameter $\xi =0.1$. Regarding the learning, 100 episodes were taken for BO, LIPO and DDPG. For the GP, the upper limit is set to 1200, considering 20 generations with $\mu =30$ individuals, $\lambda =60$ offsprings and a ($\mu +\lambda$) approach.

The DDPG experiences are collected with an exploration strategy structured into three parts. The first part (until episode 30) is mostly explorative. Here the noise is clipped in the range $[-0.8,0.8]$ with $\eta =1$ (see (3.19)). The second phase (between episodes 30 and 55) is an off-policy exploration phase with a noise signal clipped in the range $[-0.25,0.25]$, with $\eta =0.25$. The third phase (from episode 55 onwards) is completely exploitative (with no noise). As an explorative signal, we used a white noise with a standard deviation of 0.5.

4.2. Control of the viscous Burgers's equation

We consider the Burger's equation because it offers a simple 1-D problem combining nonlinear advection and diffusion. The problem set is

(4.10)\begin{equation} \left. \begin{aligned} \partial_tu + u\partial_xu & = \nu\partial_{xx}u + f(x,t) + c(x,t),\\ u(x,0) & = u_0,\\ \partial_xu(0,t) & = \partial_xu(L,t)=0, \end{aligned} \right\} \end{equation}

where $(x,t)\in (0,L)\times (0,T]$ with $L=20$ and $T=15$ is the episode length, $\nu =0.9$ is the kinematic viscosity and $u_0$ is the initial condition, defined as the developed velocity field at $t=2.4$ starting from $u(x,0)=0$. The term $f(x,t)$ represents the disturbance and the term $c(x,t)$ is the control actuation, which are both Gaussian functions in space, modulated by a time-varying amplitude,

(4.11)$$\begin{gather} f(x,t)=A_f\sin{(2{\rm \pi} f_p t)}\ {\cdot}\ \mathcal{N}(x-x_f,\sigma), \end{gather}$$

(4.12)$$\begin{gather}c(x,t)=a(t)A_c\ {\cdot}\ \mathcal{N}(x- x_c,\sigma), \end{gather}$$

taking $A_f= 100$ and $f_p=0.5$ for the disturbance's amplitude and frequencies and $A_c=300$ being the amplitude of the control and $a(t)\in [-1,1]$ the action provided by the controller. The disturbance and the controller action are centred at $x_f=6.6$ and $x_c=13.2$, respectively, and have $\sigma =0.2$. The uncontrolled system produces a set of nonlinear waves propagating in both directions at approximately constant velocities. The objective of the controller is to neutralize the waves downstream of the control location, i.e. for $x>x_c$, using three observations at $x=8,9,10$. Because the system's characteristic is such that perturbations propagate in both directions, the impact of the controller propagates backwards towards the sensors and risks being retrofitted in the loop.

To analyse how the various agents deal with the retrofitting problem, we consider two scenarios: a ‘fully closed-loop’ approach and a ‘hybrid’ approach, in which agents are allowed to produce a constant action. The constant term allows for avoiding (or at least limiting) the retrofitting problem. For the BO and LIPO controllers, we consider linear laws; hence, the first approach is

(4.13)\begin{equation} a_A(t;\boldsymbol{w})=w_0 \,u(8,t)+w_1\, u(9,t)+w_2\, u(10,t), \end{equation}

while the second is

(4.14)\begin{equation} a_B(t;\boldsymbol{w})=w_0 \,u(8,t)+w_1\, u(9,t)+w_2\, u(10,t)+w_3. \end{equation}

For the GP, we add the possibility of a constant action using an ephemeral constant, which is a function with no argument that returns a random value. Similarly, we refer to ‘A’ and ‘B’ as agents that cannot produce a constant and those that do. For the DDPG, the ANN used to parametrize the policy naturally allows for a constant term; hence, the associated agent is ‘hybrid’ by default, and there is no distinction between A and B.

One can get more insights into the dynamics of the system and the role of the controller from the energy equation associated with (4.11). This equation is obtained by multiplying (4.10) by $u$,

(4.15)\begin{equation} \partial_t \mathcal{E} + u \partial_x \mathcal{E} = {\nu}[\partial_{xx}\mathcal{E}-\left(\partial_x u\right)^2] +2 u\,f(x,t) + 2 u\,c(x,u), \end{equation}

where $\mathcal {E}=u^2$ is the transported energy and $u\,f(x,t)$ and $u\,c(x,u)$ are the production/destruction terms associated to the forcing action and the control action. Because $f$ and $c$ do not act in the same location, the controller cannot act directly on the source, but must rely either on the advection (mechanism I) or the diffusion (mechanism II). The first mechanism consists of sending waves towards the disturbing source so that they are annihilated before reaching the control area. Producing this backward propagation in a fully closed-loop approach is particularly challenging. This is why we added the possibility of an open-loop term. The second mechanism generates large wavenumbers, that is, waves characterized by large slopes so that the viscous term (and precisely the squared term in the brackets on the right-hand side of (4.15)) provides more considerable attenuation. This second mechanism cannot be used by a linear controller, whose actions cannot change the frequency from the sensors’ observation.

The controller's performance is measured by the reward function

(4.16)\begin{equation} r(t) ={-}(\ell_2(u_t)_{\varOmega_r} + \alpha{\cdot} a(t)^2), \end{equation}

where $\ell _2(\ {\cdot }\ )_{\varOmega _r}$ is the Euclidean norm of the displacement $u_t$ at time step $t$ over a portion of the domain $\varOmega _r = \{x\in \mathbb {R}|15.4\leq x\leq 16.4\}$ called the reward area, $\alpha$ is a penalty coefficient and $a_t$ is the value of the control action selected by the controller. The cumulative reward is computed with a discount factor $\gamma =1$ while the penalty in the actions was set to $\alpha =100$. This penalty gives comparable importance to the two terms in (4.16) for the level of wave attenuation achieved by all agents. Figure 5 shows the evolution of the uncontrolled system in a contour plot in the space–time domain, recalling the location of perturbation, action, observation and reward area.

Figure 5.

Contour plot of the spatio-temporal evolution of normalized $\hat {u}=u/\max (u)$ in (4.10) for the uncontrolled problem, i.e. $c(x,t)=0$ in the normalized space–time domain ($\hat {x}=x/L$, $\hat {t}=t/T$). The perturbation is centred at $\hat {x}=0.33$ (red continuous line) while the control law is centred at $\hat {x}=0.66$ (red dotted line). The dashed black lines visualize the location of the observation points, while the region within the white dash-dotted line is used to evaluate the controller performance.

Equation (4.10) was solved using Crank–Nicolson's method. The Neumann boundary conditions are enforced using ghost cells, and the system is solved at each time step via the banded matrix solver solve_banded from the python library scipy. The mesh consists of $n_x=1000$ points and the time stepping is $\Delta t=0.01$, thus leading to $n_t=1500$ steps per episode.

Both LIPO and BO optimizers operate within the bounds $[-0.1, 0.1]$ for the weights to avoid saturation in the control action. The overall set-up of these agents is the same as that used in the 0-D test case. For the GP, the selected evolutionary strategy is $(\mu +\lambda )$, with the initial population of 10 individuals $\mu = 10$ and an offspring $\lambda = 20$ trained for 20 generations. The DDPG agent set-up relies on the same reward normalization and buffer prioritization presented for the previous test case. However, the trade-off between exploration and exploitation was handled differently: the random noise term in (3.19) is set to zero every $N=3$ episodes to prioritize exploitation. This noise term was taken as an Ornstein–Uhlenbeck, time-correlated noise with $\theta = 0.15$ and $\textrm {d}t = 1\times 10^{-3}$ and its contribution was clipped in the range $[-0.3, 0.3]$. Regarding the learning, the agent was trained for 30 episodes.

4.3. Control of the von Kármán street behind a 2-D cylinder

The third test case consists in controlling the 2-D viscous and incompressible flow past a cylinder in a channel. The flow past a cylinder is a classic benchmark for bluff body wakes (Zhang et al. Reference Zhang, Fey, Noack, König and Eckelmann1995; Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003), exhibiting a supercritical Hopf bifurcation leading to the well-known von Kármán vortex street. The cylinder wake configuration within a narrow channel has been extensively used for computational fluid dynamics benchmark purposes (Schäfer et al. Reference Schäfer, Turek, Durst, Krause and Rannacher1996) and as a test case for flow control techniques (Rabault et al. Reference Rabault, Kuchta, Jensen, Réglade and Cerardi2019; Tang et al. Reference Tang, Rabault, Kuhnle, Wang and Wang2020; Li & Zhang Reference Li and Zhang2021).

We consider the same control problem as in Tang et al. (Reference Tang, Rabault, Kuhnle, Wang and Wang2020), sketched in figure 6. The computational domain is a rectangle of width $L$ and height $H$, with a cylinder of diameter $D=0.1$ m located slightly off the symmetric plane of the channel (cf. figure 6). This asymmetry triggers the development of vortex shedding.

Figure 6.

Geometry and observations probes for the 2-D von Kármán street control test case. The 256 observations used by Tang et al. (Reference Tang, Rabault, Kuhnle, Wang and Wang2020) are shown with black markers. These are organized in three concentric circles (diameters $1+0.002/D$, $1+0.02D$ and $1+0.05D$) around the cylinder and three grids (horizontal spacing $c_1 = 0.025/D$, $c_2 = 0.05/D$ and $c_3 = 0.1/D$). All the grids have the same vertical distance between adjacent points ($c_4 = 0.05/D$). The five observations used in this work (red markers) have coordinates $s_1(0,-1.5)$, $s_2(0,1.5)$, $s_3(1,-1)$, $s_4(1,1)$ and $s_5(1,0)$. Each probe samples the pressure field.

The channel confinement potentially leads to different dynamics compared with the unbounded case. Depending on the blockage ratio ($b=D/H$), low-frequency modes might be damped, promoting the development of high frequencies. This leads to lower critical Reynolds and Strouhal numbers (Kumar & Mittal Reference Kumar and Mittal2006; Singha & Sinhamahapatra Reference Singha and Sinhamahapatra2010), the flattening of the recirculation region and different wake lengths (Williamson Reference Williamson1996; Rehimi et al. Reference Rehimi, Aloui, Nasrallah, Doubliez and Legrand2008). However, Griffith et al. (Reference Griffith, Leontini, Thompson and Hourigan2011) and Camarri & Giannetti (Reference Camarri and Giannetti2010) showed, through numerical simulations and Floquet stability analysis, that for $b=0.2$ ($b\approx 0.24$ in our case), the shedding properties are similar to those of the unconfined case. Moreover, it is worth stressing that the flow is expected to be fully three dimensional for the set of parameters considered here Kanaris, Grigoriadis & Kassinos (Reference Kanaris, Grigoriadis and Kassinos2011); Mathupriya et al. (Reference Mathupriya, Chan, Hasini and Ooi2018). Therefore, the 2-D test case considered in this work is a rather academic benchmark, yet characterized by rich and complex dynamics (Sahin & Owens Reference Sahin and Owens2004) reproducible at a moderate computational cost.

The reference system is located at the centre of the cylinder. At the inlet ($x=-2D$), as in Schäfer et al. (Reference Schäfer, Turek, Durst, Krause and Rannacher1996), a parabolic velocity profile is imposed,

(4.17)\begin{equation} u_{inlet} = \frac{-4U_m}{H^2}(y^2 - 0.1Dy - 4.2D^2), \end{equation}

where $U_m=1.5\,\textrm {m}\,\textrm {s}^{-1}$. This leads to a Reynolds number of $Re=\bar {U} D/\nu =400$ using the mean inlet velocity $\bar {U}=2/3 U_m$ as a reference and taking a kinematic viscosity of $\nu =2.5\times 10^{-4}\,\textrm {m}^2\,\textrm {s}^{-1}$. It is worth noting that this is much higher than $Re=100$ considered by Jin et al. (Reference Jin, Illingworth and Sandberg2020), who defines the Reynolds number based on the maximum velocity.

The computational domain is discretized with an unstructured mesh refined around the cylinder, and the incompressible Navier–Stokes equations are solved using the incremental pressure correction scheme method in the FEniCS platform (Alnæs et al. Reference Alnæs, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015). The mesh consists of 25 865 elements and the simulation time step is set to $\Delta t=1e-4[s]$ to respect the Courant–Friedrichs–Lewy condition. The reader is referred to Tang et al. (Reference Tang, Rabault, Kuhnle, Wang and Wang2020) for more details on the numerical set-up and the mesh convergence analysis.

In the control problem every episode is initialized from a snapshot that has reached a developed shedding condition. This was computed by running the simulation without control for $T = 0.91$ s $= 3T^*$, where $T^*= 0.303$ s is the vortex shedding period. We computed $T^*$ by analysing the period between consecutive pressure peaks observed by probe $s_5$ in an uncontrolled simulation. The result is the same as that found by Tang et al. (Reference Tang, Rabault, Kuhnle, Wang and Wang2020), who performed a discrete Fourier transform of the drag coefficient.

The instantaneous drag and lift on the cylinder are calculated via the surface integrals:

(4.18a,b)\begin{equation} F_D = \int(\sigma{\cdot} n)\ {\cdot}\ e_x\, {\rm d}S, \quad F_L = \int\,(\sigma{\cdot} n)\ {\cdot}\ e_y\, {\rm d}S, \end{equation}

where $S$ is the cylinder surface, $\sigma$ is the Cauchy stress tensor, $n$ is the unit vector normal to the cylinder surface, $e_x$ and $e_y$ are the unit vectors of the $x$ and $y$ axes, respectively. The drag and lift coefficient are calculated as $C_D = {2F_D}/({\rho \bar {U}^2D})$ and $C_L = {2F_L}/({\rho \bar {U}^2D})$, respectively.

The control action consists in injecting/removing fluid from four synthetic jets positioned on the cylinder boundary as shown in figure 7. The jets are symmetric with respect to the horizontal and vertical axes. These are located at $\theta =75^{\circ }, 105^{\circ }, 255^{\circ }, 285^{\circ }$ and have the same width $\Delta \theta = 15^{\circ }$. The velocity profile in each of the jets is taken as

(4.19)\begin{equation} u_{jet}(\theta) = \frac{\rm \pi}{\Delta \theta D}Q_i^*\cos{\left(\frac{\rm \pi}{\Delta \theta}(\theta - \theta_i)\right)}, \end{equation}

where $\theta _i$ is the radial position of the $i$th jet and $Q^*_i$ is the imposed flow rate. Equation (4.19) respects the non-slip boundary conditions at the walls. To ensure a zero-net mass injection at every time step, the flow rates are mean shifted as $Q_i^* = Q_i - \bar {Q}$ with $\bar {Q}=\tfrac {1}{4}\sum _{i}^{4}\,Q_i$ the mean value of the four flow rates.

Figure 7.

Location of the four control jets for the 2-D von Kármán street control test case. These are located at $\theta =75^{\circ }, 105^{\circ }, 255^{\circ }, 285^{\circ }$ and have width $\Delta \theta = 15^{\circ }$. The velocity profile is defined as in (4.19), with flow rate defined by the controller and shifted to have zero-net mass flow.

The flow rates in the four nozzles constitute the action vector, i.e. $\boldsymbol {a}=[Q_1,Q_2,Q_3,Q_4]^\textrm {T}$ in the formalism of § 2. To avoid abrupt changes in the boundary conditions, the control action is kept constant for a period of $T_c=100\Delta t= 1\times 10^{-2}\,\textrm {s}$. This is thus equivalent to having a moving average filtering of the controller actions with an impulse response of length $N=10$. The frequency modulation of such a filter is

(4.20)\begin{equation} H(\omega)=\frac{1}{10}\left|\frac{\sin(5\omega)}{\sin(\omega/2)}\right|, \end{equation}

with $\omega =2{\rm \pi} f/f_s$. The first zero of the filter is located at $\omega =2{\rm \pi} /5$, thus $f=f_s/5=2000\,\textrm {Hz}$, while the attenuation at the shedding frequency is negligible. Therefore, this filtering allows the controller to act freely within the range of frequencies of interest to the control problem, while preventing abrupt changes that might compromise the stability of the numerical solver. Each episode has a duration of $T=0.91\,\textrm {s}$, corresponding to $2.73$ shedding periods in uncontrolled conditions. This allows for having 91 interactions per episode (i.e. 33 interactions per vortex shedding period).

The actions are linked to the pressure measurements (observations of the flow) in various locations. In the original environment by Tang et al. (Reference Tang, Rabault, Kuhnle, Wang and Wang2020), 256 probes were used, similarly to Rabault et al. (Reference Rabault, Kuchta, Jensen, Réglade and Cerardi2019). The locations of these probes are shown in figure 6 using black markers. In this work we reduce the set of probes to $n_s=5$. A similar configuration was analysed by Rabault et al. (Reference Rabault, Kuchta, Jensen, Réglade and Cerardi2019) although using different locations. In particular, we kept the probes $s_1$ and $s_2$ at the same $x$ coordinate, but we moved them further away from the cylinder wall to reduce the impact of the injection on the sensing area. Moreover, we slightly moved the sensors $s_3, s_4, s_5$ downstream to regions where the vortex shedding was stronger. The chosen configuration has no guarantee of optimality and was heuristically defined by analysing the flow field in the uncontrolled configuration. Optimal sensor placement for this configuration is discussed by Paris, Beneddine & Dandois (Reference Paris, Beneddine and Dandois2021).

The locations used in this work are recalled in figure 6. The state vector, in the formalism of § 2, is thus the set of pressure at the probe locations, i.e. $\boldsymbol {s}=[p_1,p_2,p_3,p_4,p_5]^\textrm {T}$. For the optimal control strategy identified via the BO and LIPO algorithms in §§ 3.1.1 and 3.1.2, a linear control law is assumed, hence $\boldsymbol {a}=\boldsymbol {W}\boldsymbol {s}$, with the 20 weight coefficients labelled as

(4.21)\begin{equation} \begin{bmatrix}Q_1 \\Q_2\\Q_3\\Q_4 \end{bmatrix} = \begin{bmatrix}w_1 & w_2 & w_3 & w_4 & w_5 \\w_6 & w_7 & w_8 & w_9 & w_{10} \\w_{11} & w_{12} & w_{13} & w_{14} & w_{15} \\w_{16} & w_{17} & w_{18} & w_{19} & w_{20} \end{bmatrix} \begin{bmatrix}p_1 \\p_2\\p_3\\p_4\\p_5 \end{bmatrix} . \end{equation}

It is worth noting the zero-net mass condition enforced by removing the average flow rate from each action could be easily imposed by constraining all columns of $\boldsymbol {W}$ to add up to zero. For example, setting the symmetry $w_1=-w_{11}$, $w_6=-w_{16}$, etc.(leading to $Q_1=-Q_3$ and $Q_2=-Q_4$) allows for halving the dimensionality of the problem and, thus, considerably simplifying the optimization. Nevertheless, one has infinite ways of embedding the zero-net mass condition and we do not impose any, letting the control problem act in $\mathbb {R}^{20}$.

Finally, the instantaneous reward $r_t$ is defined as

(4.22)\begin{equation} r_t = \langle F_{D}^{base}\rangle_{T_c} -\langle F_D\rangle_{T_c} - \alpha|\langle F_L\rangle_{T_c}|, \end{equation}

where $\langle \bullet \rangle _{T_c}$ is the moving average over $T_c=10\Delta t$, $\alpha$ is the usual penalization parameter set to $0.2$ and $F_{D}^{base}$ is the averaged drag due to the steady and symmetric flow. This penalization term prevents the control strategies from relying on the high lift flow configurations Rabault et al. (Reference Rabault, Kuchta, Jensen, Réglade and Cerardi2019) and simply blocking the incoming flow. The cumulative reward was given with $\gamma =1$. According to Bergmann et al. (Reference Bergmann, Cordier and Brancher2005), the active flow control cannot reduce the drag due to the steady flow, but only the one due to the vortex shedding. Hence, in the best case scenario, the cumulative reward is the sum of the averaged steady state drag contributions:

(4.23)\begin{equation} R^* = \sum_{t=1}^T r_t = \sum_{t=1}^T \langle F_{D}^{base}\rangle_{T_c} = 14.5. \end{equation}

The search space for the optimal weights in LIPO and BO was bounded to $[-1, 1]$. Moreover, the action resulting from the linear combination of such weights with the states collected in the $i$th interaction was multiplied by a factor $2\times 10^{-3}$, to avoid numerical instabilities. The BO settings are the same as in the previous test cases, except for the smoothness parameter that was reduced to $\nu = 1.5$. On the GP side, the evolutionary strategy applied was the eaSimple's (Bäck et al. Reference Bäck, Fogel and Michalewicz2018) implementation in DEAP – with hard-coded elitism to preserve the best individuals. To allow the GP to provide multi-outputs, four populations of individuals were trained simultaneously (one for each control jet). Each population evolves independently (with no genetic operations allowed between them), although the driving reward function (4.23) values their collective performance. This is an example of multi-agent RL. Alternative configurations, to be investigated in future works, are the definition of a multiple-output trees or cross-population genetic operations.

Finally, the DDPG agent was trained using the same exploration policy of the Burgers’ test case, alternating 20 exploratory episodes with $\eta =1$ and 45 exploitative episodes with $\eta =0$ (cf. (4.22)). During the exploratory phase, an episode with $\eta =0$ is taken every $N=4$ episodes and the policy weights are saved. We used the Ornstein–Uhlenbeck time correlated noise with $\theta = 0.1$ and $\textrm {d}t = 1\times 10^{-2}$ in (3.19), clipped in the range $[-0.5, 0.5]$.