2.1. High-fidelity model based on the unsteady Navier–Stokes equations

The simulation of the fluid–solid interaction between the airfoil and the surrounding flow is performed with STAR-CCM+. The numerical approach is similar to the one used in Jardin and Doué (Reference Jardin and Doué2019). The evolution of the flow follows the unsteady incompressible Navier–Stokes equations. These are solved with a cell-centered finite volume method. The numerical schemes are second-order both in time and space. The flat plate is placed in the center of a circular domain, whose radius is 20 times the chord length of the airfoil. The latter is modeled as a no slip surface and the far-field is set through a null velocity condition.

The flat plate rotation is imposed by the controller through the pitch rate $ \dot{\beta} $ . However, the airfoil is free to move by translation in both directions of the 2D space, based on its interaction with the fluid. This fluid–solid interaction is simulated with the 6-DOF solver of STAR-CCM+. The airfoil is submitted to the lift $ \overrightarrow{L} $ and drag $ \overrightarrow{D} $ aerodynamic forces, as well as gravity $ \overrightarrow{g} $ , as displayed in Figure 2. All angles are measured as displayed in Figure 2 and signed counter-clockwise.Footnote ³

Figure 2. Lift $ \overrightarrow{L} $ , drag $ \overrightarrow{D} $ _, and gravity $ \overrightarrow{g} $ forces applied to the center of mass of the flat plate. The flat plate velocity is denoted $ \overrightarrow{V} $ , the pitch angle $ \beta $ and the angle of attack $ \alpha $ . All angles are measured as displayed in this figure and signed counter-clockwise: $ \alpha $ is the angle going from the flat plate to the velocity vector so that $ \alpha >0 $ in the figure and $ \beta $ is the angle going from the horizontal ( $ -\overrightarrow{x} $ ) to the flat plate so that $ \beta <0 $ in the figure.

Reducing the flat plate to its center of mass, the fundamental principle of dynamics is expressed as

(1) $$ m\cdot \overrightarrow{a}=m\cdot \overrightarrow{g}+\overrightarrow{L}+\overrightarrow{D} $$

where $ m $ is the flat plate mass, $ \overrightarrow{a} $ is the flat plate acceleration, $ \overrightarrow{g} $ is the acceleration of gravity, $ \overrightarrow{L}=\left({L}_x,{L}_y\right) $ is the lift force, and $ \overrightarrow{D}=\left({D}_x,{D}_y\right) $ is the drag force. Projecting Eq. (1) onto the $ x $ and $ y $ axes of a reference frame attached to the ground leads to the following equations of motion:

(2) $$ {\displaystyle \begin{array}{c}m\frac{d^2x}{dt^2}={L}_x+{D}_x,\\ {}m\frac{d^2y}{dt^2}=- mg+{L}_y+{D}_y.\end{array}} $$

In this high-fidelity physics model, aerodynamic forces $ \overrightarrow{L} $ and $ \overrightarrow{D} $ are computed by integration of the pressure and viscous forces over the surface of the airfoil. These forces, along with the gravitational force are then used to compute the translational motion of the flat plate by integration of the equations of motion (Eq. 2) over small time steps of duration $ dt $ .

The simulation setup, mesh, and time step used to solve the flow are chosen in accordance with the work of Jardin and Doué (Reference Jardin and Doué2019), where the numerical approach has been validated for highly unsteady cases. The flow is solved on a structured mesh. Mesh convergence is assessed by computing the flow for a constant angle of attack equal to $ \alpha ={40}^{\circ } $ and a constant translation velocity equal to $ \overrightarrow{U_0}=-10\overrightarrow{x}\mathrm{m}\;{\mathrm{s}}^{-1} $ , with three meshes of increasing resolution. The cell size at the wall, as well as the number of cells for each mesh, are given in Table 1.

Table 1. Mesh characteristics

The airfoil lift and drag coefficients obtained with the two most refined meshes showed variations of values below 1.5% compared to the coarse mesh. Therefore, the coarse mesh is selected for this study. Similarly, simulations performed with a time step equal to $ dt=5\times {10}^{-5}\mathrm{s} $ , $ dt=2.5\times {10}^{-5}\mathrm{s} $ or $ dt=1.25\times {10}^{-5}\mathrm{s} $ , exhibited differences on the aerodynamic coefficients below 0.1%. The larger time step of $ dt=5\times {10}^{-5}\mathrm{s} $ is thus retained to perform uncontrolled simulations at fixed angles of attack. However, to ensure a correct temporal resolution of the pitch motion, and for stability reasons, a smaller time step $ dt=2\times {10}^{-5}\mathrm{s} $ is set for simulations controlled by the DRL algorithm (i.e., when the pitch evolves during the simulation). Controlled computations are initialized with the converged solution obtained at a fixed angle of attack $ \alpha ={0}^{\circ } $ and a Reynolds number $ \mathit{\operatorname{Re}}=300 $ .

2.2. Low-order physics model: Quasi-steady empirical approach

In order to decrease the computational cost associated with the simulation of the flat plate motion, the high-fidelity model described above can be replaced by a simpler model to compute the aerodynamic forces involved in Eq. (2). This quasi-steady surrogate model replaces the expensive computation of the flow by empirical correlations between aerodynamic forces and the airfoil angle of attack. Such correlations are obtained by solving the Navier–Stokes equations at various fixed angles of attack with STAR-CCM+ and extracting the corresponding lift and drag coefficients. The evolution of the time-averaged lift $ {c}_L $ and drag $ {c}_D $ coefficients with the angle of attack is shown in Figure 3. These quasi-steady data reasonably agree with numerical data from Taira and Colonius (Reference Taira and Colonius2009) obtained for the same Reynolds number and a 2D flat plate whose thickness is 3.66% of the chord length.Footnote ⁴

Figure 3. Evolution of the lift $ {c}_L $ and drag $ {c}_D $ coefficients with the angle of attack $ \alpha $ . Results from Taira and Colonius (Reference Taira and Colonius2009) (■), present study using STAR-CCM+ (○), and fitting of the present study data using Eq. (3) (- - -).

As mentioned in Wang et al. (Reference Wang, Birch and Dickinson2004), data from Figure 3 can be approximated with trigonometric basis functions. To do so, the simplex method described by Lagarias et al. (Reference Lagarias, Reeds, Wright and Wright1998) is used. The fitting curves (- - - in Figure 3) are given by

(3) $$ {\displaystyle \begin{array}{r}{c}_L=1.3433\sin \left(2\alpha \right),\\ {}{c}_D=1.5055-1.3509\cos \left(2\alpha \right).\end{array}} $$

In order to solve the equations of motion (Eq. 2), these coefficients can then be used to compute the lift and drag forces applied to the flat plate at each time step $ dt $ :

(4) $$ {\displaystyle \begin{array}{l}L=\frac{1}{2}{d}_f{U}^2\mathcal{S}\cdot {c}_L,\\ {}D=\frac{1}{2}{d}_f{U}^2\mathcal{S}\cdot {c}_D.\end{array}} $$

where $ U $ is the norm of the instantaneous velocity, and $ \mathcal{S} $ is equal to the chord length times the span of the airfoil. This low-order model permits the simulation of the system’s dynamics at a reduced computational cost compared to the resolution of the Navier–Stokes equations. Such cost reductions are critical to perform the detailed study of both the control problem formulation and the DRL algorithms parametrization presented in this article.

3.1. Control problem definition

In Section 2, the flat plate motion is described through the values of its position variables $ x $ and $ y $ in a reference frame attached to the ground. However, the DRL algorithm takes a self-centered perspective: it should not build a control policy based on absolute positions in space, but instead on the relative position of the airfoil to the target point B. For that purpose, a change of variables is performed to define the systems’ state $ {s}_t\in S $ at any time step. The flat plate position and velocity are expressed in a polar coordinate system $ \left(\rho, \theta \right) $ , relative to point B, as illustrated in Figure 4.

Figure 4. Airfoil polar coordinates $ \left(\rho, \theta \right) $ defined relative to point B.

The normalized state at time step $ t $ then reads:

(5) $$ {s}_t=\left(\frac{\rho_t}{k_1},\sin \left({\theta}_t\right),\cos \left({\theta}_t\right),\frac{{\dot{\rho}}_t}{k_2},\frac{{\dot{\theta}}_t}{k_3},\sin \left({\beta}_t\right)\right) $$

where $ {k}_1 $ , $ {k}_2 $ , and $ {k}_3 $ are normalizing factors chosen in order to maintain the same order of magnitude within the input variables of control policies and value functions, to avoid scaling effects which are detrimental to stochastic gradient descent learning procedures. These state variables are chosen because they define an intuitive first-order dynamics on variables that are relative to the target position, where $ {s}_{t+1} $ can be predicted from $ {s}_t $ and the current pitch rate $ {\dot{\beta}}_t $ . At each discrete control time step $ t $ , the RL algorithm observes the state $ {s}_t $ and controls the flat plate trajectory by setting the pitch rate $ {\dot{\beta}}_t $ . This pitch rate defines the controller’s action variable $ {a}_t\in A $ at each time step, and may take any value in the interval $ \left[-{\dot{\beta}}_0;{\dot{\beta}}_0\right] $ . The time span between two control steps is $ \Delta t $ . The low-order model of Section 2.2 induces a first-order dynamic on the state $ {s}_t $ such that the state at time step $ t+1 $ is only function of the state and action at time step $ t $ : $ {s}_{t+1}=f\left({s}_t,{a}_t\right) $ , for a constant action equal to $ {a}_t $ over the $ \Delta t $ time interval between time steps.Footnote ⁵

Under a certain control strategy $ \pi $ prescribing an action $ {a}_t=\pi \left({s}_t\right) $ at each time step, the system follows a trajectory $ \left({s}_0,{a}_0,{s}_1,{a}_1,\dots \right) $ . To devise optimal control strategies, one introduces a quantitative criterion $ J\left(\pi \right) $ that evaluates the induced trajectories. This criterion is defined as a sum of instantaneous rewards $ {r}_t=r\left({s}_t,{a}_t,{s}_{t+1}\right) $ obtained at each control time step along the trajectory, averaged across possible starting states $ {s}_0 $ . For a given starting state $ {s}_0 $ , this sum of rewards is called the return under policy $ \pi $ . Introducing a discount factor $ \gamma <1 $ , quantifying the importance between immediate and future rewards, $ J\left(\pi \right)={\unicode{x1D53C}}_{s_0}\left[{\sum}_{t=0}^{\infty }{\gamma}^t{r}_t|{s}_0\right] $ is the average over-starting state of the discounted return when controlled with the strategy $ \pi $ . Therefore, the definition of the reward model shapes the criterion which is optimized by the RL algorithm.

In the airfoil trajectory problem, the goal is to reach point B. A first solution, when defining the reward model is to give a higher value to a trajectory that strongly reduces $ \rho $ , and thus choose a model under the form:

(6) $$ r\left({s}_t,{a}_t,{s}_{t+1}\right)=-\Delta \rho \left({s}_t,{s}_{t+1}\right) $$

where $ \Delta \rho $ is the variation of $ \rho $ between $ {s}_t $ and $ {s}_{t+1} $ . Such a formulation evaluates whether the airfoil moved closer or farther away from B over the last control time step. It is called dense since a non-null feedback signal is given at each step. For such a reward model, the return $ J\left(\pi \right) $ , over a trajectory, indicates the overall expected reduction in the distance to B, that is $ -\Delta \rho \left({s}_0,{s}_{\infty}\right) $ , where $ {s}_{\infty } $ is the state reached by the flat plate as $ t\to \infty $ .

An alternative reward model consists in giving a positive reward only when the flat plate reaches the vicinity of B, which takes the form

(7) $$ r\left({s}_t,{a}_t,{s}_{t+1}\right)=\left\{\begin{array}{ll}1& \hskip0.1em \mathrm{if}\hskip0.22em \rho \left({s}_{t+1}\right)<\varepsilon, \\ {}0& \hskip0.1em \mathrm{otherwise}.\end{array}\right. $$

where $ \varepsilon $ is a tolerance constant. Note that this reward definition will produce a null signal most of the time, since non-zero values will be obtained only when the airfoil reaches B: this reward is therefore sparse and can be challenging for most DRL algorithms.

In this work, the reward model includes both a dense and a sparse term:

(8) $$ r\left({s}_t,{a}_t,{s}_{t+1}\right)={c}_1\cdot {r}_d\left({s}_t,{a}_t,{s}_{t+1}\right)+{c}_2\cdot {r}_s\left({s}_t,{a}_t,{s}_{t+1}\right), $$

where

(9) $$ {r}_d\left({s}_t,{a}_t,{s}_{t+1}\right)=-\frac{\Delta \rho \left({s}_t,{s}_{t+1}\right)}{\rho_0}, $$

(10) $$ {r}_s\left({s}_t,{a}_t,{s}_{t+1}\right)=\left\{\begin{array}{ll}1& \hskip0.1em \mathrm{if}\hskip0.22em \rho \left({s}_{t+1}\right)<\varepsilon, \\ {}-1& \hskip0.1em \mathrm{if}\hskip0.22em \mid \theta \mid \ge \pi /2,\\ {}0& \hskip0.1em \mathrm{otherwise},\end{array}\right. $$

$ {c}_1 $ and $ {c}_2 $ are constant parameters and the tolerance $ \varepsilon $ is equal to $ \varepsilon ={10}^{-2}/{\rho}_0 $ . Note that compared to Eq. (7), a penalty has been added in the sparse term (Eq. 10) if $ \mid \theta \mid \ge \pi /2 $ , that is, when the airfoil goes left beyond point B. In Eq. (8), the respective weights of the sparse and dense terms depend on the values assigned to $ {c}_1 $ and $ {c}_2 $ . In particular, setting $ {c}_1=0 $ (respectively, $ {c}_2=0 $ ), leads to a purely sparse (resp. dense) reward model. Various values of these parameters are discussed in this paper.

Whatever the reward parametrization, an episode (or trajectory) may terminate in three different ways:

• • $ \rho \left({s}_{t+1}\right)<\varepsilon $ means that the airfoil reached B. This situation is called won (the “game” is won).

• • $ \mid \theta \mid \ge \pi /2 $ means that the airfoil surpassed B. This situation is called lost.

• • the maximum number of steps, set to 150, has been reached. The episode is also considered as lost.

If a won or lost condition is reached, the episode ends, and the airfoil state and parameters are reset to their initial values.

To sum up, at the beginning of an episode, the flat plate position is set to point A and the velocity to the initial velocity. The RL algorithm observes the state $ {s}_t $ and chooses an action $ {a}_t $ . Then, either the low-order model or the Navier–Stokes equations are solved to obtain the new position and velocity of the flat plate $ {s}_{t+1} $ and the instantaneous reward, associated with the transition is computed. These are fed to the RL algorithm, which again chooses an action and so on until termination. The control time step $ \Delta t $ is equal to the physical time step $ dt $ used to integrate the equations of motion of the airfoil.

Note that, in a practical application, one might desire a behavior that reaches point B with low velocity or in a specific configuration. This would require expressing a more convoluted reward model, but would not change fundamentally the practical steps of deploying RL algorithms on the problem. In the present work, we retain the simplicity and representativity of the model introduced above.

The theory of discrete-time optimal control (Bertsekas and Shreve, Reference Bertsekas and Shreve1996) guarantees that for such a valuation of trajectories and such first-order dynamics, there exists an optimal control function $ a=\pi (s) $ that only depends on $ s $ . Such dynamical systems described via the 4-tuple $ \left(S,A,f,r\right) $ are called Markov decision processes (MDP) (Puterman, Reference Puterman2014) and extend, in the general case, to stochastic transition models that follow Markov’s property. Note that, in the case where the dynamics are obtained by solving the Navier–Stokes equations with STAR-CCM+, the system still obeys a first-order differential equation, but the latter involves the full pressure and velocity fields surrounding the flat plate, which are not observed here for the control policy. In such a case, the control policy is still defined as a function of $ s $ , but without formal guarantees of optimality.

3.2. Reinforcement learning algorithms

RL algorithms solve the optimization problem $ {\max}_{\pi }J\left(\pi \right) $ , by exploiting samples $ \left({s}_t,{a}_t,{r}_t,{s}_{t+1}\right) $ obtained by interacting with the physical system. Doing so, the algorithm has no explicit knowledge of the underlying dynamics. To solve the optimization problem, one defines the value function $ {V}^{\pi } $ which represents the return from a given state of a trajectory controlled by $ \pi $ : $ {V}^{\pi }(s)={\sum}_{t=0}^{\infty }{\gamma}^t{r}_t $ provided that $ {s}_0=s $ and actions are picked according to $ \pi $ . Similarly, one can define the state-action value function $ {Q}^{\pi}\left(s,a\right) $ , which is the return when applying $ a $ in $ s $ , and thereafter applying the control policy $ \pi $ . In DRL, neural networks are used to learn approximations of the value functions $ {V}^{\pi }:S\to \mathrm{\mathbb{R}} $ and/or $ {Q}^{\pi }:S\times A\to \mathrm{\mathbb{R}} $ , and the control functions $ \pi :S\to A $ . In the present work, all neural networks are multilayer perceptrons (MLPs) as illustrated in Figures 5 and 6. Activation functions are ReLU functions for all hidden layers. Output neurons of value function networks use identity activation functions, while output neurons for policy networks use $ \tanh $ activation functions as per the common practice of each algorithm.

Figure 5. Generic architecture of the policy $ \pi $ neural network.

Figure 6. Generic architecture of the state-action value function $ {Q}^{\pi}\left(s,a\right) $ neural network.

Methods that directly solve the $ {\max}_{\pi }J\left(\pi \right) $ problem from samples fall under the umbrella of direct policy search. When they are gradient-based, such methods exploit the policy gradient theorem (Sutton et al., Reference Sutton, McAllester, Singh and Mansour1999) which necessitates repeated collection of samples under the current control policy’s stationary distribution (they are then called on-policy). The seminal proximal policy optimization (PPO) (Schulman et al., Reference Schulman, Wolski, Dhariwal, Radford and Klimov2017) algorithm, used for instance in the work of Rabault et al. (Reference Rabault, Kuchta, Jensen, Reglade and Cerardi2019), is one such policy gradient method. This makes such methods quite sample inefficient and quickly becomes a handicap for simulation-costly problems such as CFD control. On the other hand, the family of approximate dynamic programming methods (Bertsekas, Reference Bertsekas2012) solves for the optimal value function, which is the unique solution to the Bellman equation $ Q\left(s,a\right)=r\left(s,a,f\left(s,a\right)\right)+\gamma {\max}_{a^{\prime }}Q\left(f\left(s,a\right),{a}^{\prime}\right),\forall s,a $ , and then derives an optimal policy as one that is greedy with respect to this optimal value function. The Bellman equation is solved without knowledge of the $ r $ and $ f $ functions, and relies on a set of collected samples $ \left({s}_t,{a}_t,{r}_t,{s}_{t+1}\right) $ , called a replay memory or replay buffer (Lin, Reference Lin1992). This buffer is filled by following an exploration behavior that selects actions at each time step. This behavior policy should often choose actions that are deemed promising by the current value function $ Q $ , but also balance them with exploring new actions (this is coined the exploration versus exploitation trade-off). The replay memory is sampled to solve the Bellman equation through a sequence of supervised learning problems. Such methods are off-policy by construction, decoupling the behavior policy from the one being optimized, which enables sample reuse along the optimization process and makes them promising candidates for simulation-costly problems such as CFD control. In this work, three state-of-the-art approximate dynamic programming algorithms are compared, namely DDPG (Lillicrap et al., Reference Lillicrap, Hunt, Pritzel, Heess, Erez, Tassa, Silver and Wierstra2016), TD3 (Fujimoto et al., Reference Fujimoto, Hoof and Meger2018), and SAC (Haarnoja et al., Reference Haarnoja, Zhou, Abbeel and Levine2018).

All three methods share the same architecture, where an actor neural network $ \pi :S\to A $ approximates the solution of $ \arg {\max}_aQ\left(s,a\right) $ in all s. Then, a critic network $ Q:S\times A\to \mathrm{\mathbb{R}} $ approximates the mapping $ {s}_t,{a}_t\mapsto {r}_t+\gamma {Q}^{-}\left({s}_{t+1},\pi \left({s}_{t+1}\right)\right) $ for all experience samples $ \left({s}_t,{a}_t,{r}_t,{s}_{t+1}\right) $ , where $ {Q}^{-} $ is called a target network which is periodically replaced by $ Q $ . Each of these networks approximate the true greedy policy and value function by taking gradient steps towards them. This process is repeated until convergence (convergence is guaranteed in the absence of approximation errors). In DDPGFootnote ⁶, the actor network is a deterministic policy, and the behavior policy is the actor’s policy when noise is injected artificially. TD3 (twin-delayed DDPG) is an attempt to fix some weaknesses of DDPG, by delaying actor updates (TD3 performs two critic updates for a single actor update), and introducing a second critic function to account for possible over-estimation of the value function in stochastic environments. Both DDPG and TD3, despite their overall good behavior on average, still exhibit specific flaws in their optimization process. For instance, Matheron et al. (Reference Matheron, Perrin, Sigaud, Farkaš, Masulli and Wermter2020) illustrated how sparse rewards in deterministic environments could prevent altogether the convergence of these methods by inducing value function plateaus and zero-gradient updates. A state-of-the-art alternative is SAC, which makes the actor policy stochastic (the behavior policy is then the actor policy), forces this policy to imitate the soft-max of $ Q $ (instead of the hard-max in DDPG and TD3), and introduces an entropy regularization term in the resolution of the Bellman equation to make the optimization landscape smoother (Geist et al., Reference Geist, Scherrer and Pietquin2019). DDPG, TD3, and SAC are fairly close from a formal foundations perspective and often perform on-par, outperforming each other depending on the practical problem at hand. Surprisingly, in a number of applications, approximation errors or specific characteristics (such as deterministic dynamics) still make DDPG very competitive.

Common practice in RL has demonstrated that methods designed to perform better in principle and on average, fail on specific benchmarks. DRL are known to be very sensitive to the fine tuning of their hyperparameters and fair evaluation of methods is an open problem (Henderson et al., Reference Henderson, Islam, Bachman, Pineau, Precup and Meger2018). This sensitivity is an undesirable aspect for the fluid mechanics practitioner, which we try to quantify on the benchmark case described in Section 2.

4.1. Baseline RL controller

First, the case illustrated in Figure 7 is considered, where A and B are fixed during training and testing. The starting and target point coordinates are, respectively, equal to $ A=\left(0,0\right) $ and $ B=\left(-2c,0.2c\right) $ , and the flat plate initial velocity $ \overrightarrow{U_0} $ is purely horizontal and directed to the left.

Figure 7. Illustration of the first task: A and B are fixed.

Since B is higher than A, the controller has to learn how to create lift in order to gain altitude. The baseline reward model follows Eq. (8) with $ {c}_1=100 $ and $ {c}_2=10 $ . Therefore, for each transition, the agent receives a reward:

(11) $$ r\left({s}_t,{a}_t,{s}_{t+1}\right)=-100\frac{\Delta \rho \left({s}_t,{s}_{t+1}\right)}{\rho_0} $$

and an additional term at the end of the episode:

(12) $$ r\left({s}_t,{a}_t,{s}_{t+1}\right)=\left\{\begin{array}{ll}10& \hskip0.1em \mathrm{if}\hskip0.22em \rho \left({s}_{t+1}\right)<\varepsilon, \\ {}-10& \hskip0.1em \mathrm{if}\hskip0.22em \mid \theta \mid \ge \pi /2,\\ {}0& \hskip0.1em \mathrm{otherwise}.\end{array}\right. $$

For a path reaching the vicinity of B, the reward associated with each transition is slightly below 1 and the number of steps around 115. As a result, return values around 90 mean that the airfoil missed B by a small margin, and values close to 110 correspond to the flat plate managing to reach point B.

For the first experiment, the physics is solved with the low-order model and the replay buffer capacity is set to 500,000, which allows to keep track of all transitions experienced during training. The policy is trained for 1,000 episodes and then tested without changing the position of A or B. The learning process for one run performed with SAC is illustrated in Figure 8. During training, the learned policy is monitored by performing a test (also called an evaluation) at the end of each episode. The evolution of the return obtained from these evaluations is shown in Figure 8a, while all trajectories explored during learning are plotted in Figure 8b. These trajectories allow the visualization of all different positions contained in the transitions stored in the replay buffer. Recall that the state is six-dimensional (Eq. 5) and that plotting the traces in the $ \left(x,y\right) $ plane does not fully reflect the trajectories in the state space. In Figure 8, curves from both plots are colored by episode number to link the evaluation return and the explored trajectories.

Figure 8. (a) Evolution of the evaluation returns during the training phase. This return is obtained through the evaluation of the policy at the end of each episode during training. To get a sense of the learning “history” and convergence, the return is plotted at each episode and no smoothing or average is performed. (b) Visualization of the 1,000 different trajectories explored during training. To link the explored trajectories and the policy learning process, curves from both plots are colored by the episode number.

The learning process may be decomposed into four phases, denoted from (i) to (iv) in Figure 8a:

1. (i) Warm-up phase: During the first eight episodes, actions are chosen randomly at each step and the policy is not trained. The corresponding trajectories stay close to the horizontal or slightly fall from point A towards the left under the combined effects of the initial velocity and gravity. This phase allows to prefill the replay buffer with various transitions.

2. (ii) Discovery phase: The agent explores different trajectories, learns how to fly, and slowly manages to control the flat plate to fly closer to the target point B.

3. (iii) Stabilizing phase: The spread between trajectories decreases and the agent produces trajectories that reach the vicinity of B. The return jumps up and down between approximately 90 (B was missed by a small margin) and 110 (B was reached).

4. (iv) Stable phase: The trajectory dispersion becomes very low. Returns stabilize just below 110: the agent has robustly learned to perform the task.

Among all the paths shown in Figure 8b, the one from the first and “best” episodes of the training are presented in Figure 9a, as well as the trajectory obtained by testing the policy at the end of the 1,000 episodes. The so-called “best” path is the one that led to the highest return and was reached at the 393th episode of the training. The actions $ \dot{\beta} $ corresponding to these three specific episodes are shown in Figure 9b and the resulting pitch angle $ \beta $ in Figure 9c.

Figure 9. (a) Trajectory, (b) normalized action $ \dot{\beta}/{\dot{\beta}}_0 $ , and (c) pitch angle $ \beta $ for three selected episodes of the training and testing: the first and best episode (i.e., with the highest return) of training and the test episode. Corresponding returns are respectively equal to 80.99, 109.78, and 109.69.

As mentioned above, actions are chosen randomly for the first episode, which can be observed in Figure 9b where the evolution of actions is stochastic. The resulting pitch angle stays close to zero, and the airfoil trajectory slightly falls while going left. The action signal associated with the best episode still appears noisy since SAC learns a Gaussian probability density function. At each time step, the action is drawn from this distribution without correlation with the action taken at the previous time step, hence the noise and the high-frequency changes in the green curve. As the physical system acts as a low pass filter, these high-frequency actions have a very moderate (but non-zero) effect on the state trajectory. Yet, the information about the effects of actions is present in the samples and helps guide the optimization towards better trajectories. On the contrary, for the test trajectory, the action is deterministically selected as the mean of the probability density function learned by SAC. This leads to a smooth sequence of actions with localized quick variations enabled by the non-linear neural network controller. For both cases, the pitch rate is negative on average over an episode. The resulting pitch angle globally decreases over the episode,Footnote ⁷ allowing the flat plate to gain lift to counteract the effect of gravity and to move upwards towards the target. Despite the differences in the policies, the returns obtained for the best and test episodes are close (109.7 and 109.5, respectively).

The reward model defined here does not strongly discriminate between policies and several quasi-optimal policies may be found for this case. This is better evidenced in Figure 10: 20 different control policies are independently trained and then tested after 1,000 episodes. All parameters remain the same for these 20 trials except the random numbers generator’s seed.

Figure 10. (a) Trajectory, (b) normalized action $ \dot{\beta}/{\dot{\beta}}_0 $ , and (c) pitch angle $ \beta $ obtained by testing 20 identically parametrized agents. Returns from all the tests are in the range $ \left[109.12;109.58\right] $ _.

Figure 10 shows different trajectories, policies, and pitch angles, which illustrates the diversity of quasi-optimal control policies for the problem at hand. However, some trends can be observed:

• • At the beginning, most policies apply a negative pitch rate, leading to a negative pitch angle which allows the airfoil to gain lift and move upwards,

• • then, the value of the angle globally stabilizes,

• • finally, close to B, the angle values of the different policies spread as the controller adapts the pitch rate to reach B.

Owing to the model of the dynamics as well as the reward formulation, a wide range of policies are quasi-optimal. As a result, small differences in the states visited during the training as well as in the actions chosen at the beginning of each path may lead to different final strategies. This explains the high spread of the test trajectories obtained within the 20 trials despite the fact that they have the same parameters. Alternative physical constraints, reward models, or parametrization of the problem may lead to more constrained strategies. Moreover, while the learning process shown in Figure 8 is smooth and while the training curves presented in Figure 10 are all successful in performing the task, other learning trials have displayed more erratic behaviors, most often due to overall instabilities during the policy learning process (e.g., catastrophic forgetting, McCloskey and Cohen, Reference McCloskey, Cohen and Bower1989; Ratcliff, Reference Ratcliff1990, or unlikely pathological exploration). The next three sections aim to explore the influence of training parameters and control problem definition on the learning reliability and sample efficiency, as well as on the resulting policy. Specifically, we isolate the influence of the RL algorithm’s parameters, the sensitivity to the reward model’s definition, and the influence of variations in the dynamics model.

4.2. Sensitivity analysis to DRL hyperparameters

In DRL, hyperparameters greatly affect policies and learning performance (Henderson et al., Reference Henderson, Islam, Bachman, Pineau, Precup and Meger2018). Besides, their effects vary with the algorithm and the control problem. Therefore, this section investigates the influence of RL hyperparameters on the learning reliability and sample efficiency of DDPG, TD3, and SAC on the first task considered in Section 4. First, the influence of the network size is assessed by varying the number of hidden layers of the MLPs, as well as the number of neurons on each layer. For each case, these variations are applied to all neural networks (policy network, value function network, and corresponding target networks, if any) at the same time. Then, the effect of the replay buffer capacity and the number of warm-up episodes is assessed. The former is of particular importance because, in a CFD context where the whole field (e.g., pressure, velocity) would be stored as the state in the replay buffer, keeping all acquired samples might cause memory issues. Note that a preliminary study investigated other parameters such as the batch size or specific techniques to compensate for poor distribution of samples in the replay buffer (Schaul et al., Reference Schaul, Quan, Antonoglou and Silver2016; Lahire et al., Reference Lahire, Geist and Rachelson2022). None of these evidenced a significant impact on the reliability and sample efficiency of the algorithms and the results are not reported in this paper.

To establish proper metrics to evaluate the reliability of RL algorithms, the inherent stochasticity of the learning process has to be taken into account. For this reason, for each set of parameter values tested in this study, 20 learning trials, with different random numbers generator’s seeds are performed. Each of them lasts 1,000 episodes. The policy obtained at the end of the 1,000th episode is then tested. The trial is classified as won if the policy is able to reach B (i.e., $ \rho <\varepsilon $ ), and lost otherwise. These 20 trials provide a statistical estimate of each configuration’s characteristics. This permits highlighting global tendencies and allows formulating guidelines for the application of RL algorithms to the present benchmark, and to some extent to other flow control problems. Table 2 summarizes the values of the parameters and the number of won trials over the 20 computed for each set of parameters.

Table 2. Influence of RL hyperparameters: cases characteristics and number of won trials over 20 realizations

Note. Number of won trials smaller or equal to 15 are highlighted in bold red, and in black bold when equal to 20.

Among the neural networks (NN) sizes tested here, cases $ {NN}^{2\times 8} $ and $ {NN}^{3\times 8} $ induce a degraded reliability for DDPG and TD3, while SAC is remarkably robust to the network definition. A small neural network might be unable to accurately encode the optimal policy to reach point B, and increases the likelihood of finding local minima during the optimization process (hence not reaching the optimal policy), while the regularization term of SAC provides a smoother optimization landscape and consistently escapes local minima. For the benchmark considered, a neural network with two hidden layers of 32 neurons is sufficient to correctly represent the policy. There is no clear-cut effect of increasing the number of layers beyond 2, or the number of neurons on the layers beyond 32. Similarly, increasing the number of warm-up episodes from 8 to 80 does not affect the reliability significantly.

Regarding the replay buffer, given the number of episodes and their average length, a replay buffer maximum capacity equal to 500,000 means that every collected sample stays in the replay buffer during the full learning. On the contrary, maximum capacities equal to 5,000 or 500 lead to a progressive replacement of the first samples by newer ones, on a first-in-first-out basis. Such configurations induce a non-stationary distribution of samples in the replay buffer, which is detrimental to the learning process and may lead to unstable training. While the reliability is only weakly affected for case $ {RB}^{5000} $ , it is strongly degraded, for DDPG and TD3, when the replay buffer capacity is set to 500. More generally, Table 2 shows that SAC is more robust to parametric variations than DDPG and TD3, even where the other two algorithms perform particularly badly.

Besides the reliability, sample efficiency is of particular interest for CFD applications because samples are often computationally expensive to collect. A sample $ \left({s}_t,{a}_t,{r}_t,{s}_{t+1}\right) $ is collected and stored in the replay buffer at each control step. Therefore, to compare the sample efficiency between algorithms and cases, Figure 11 shows the evolution of the ensemble average and standard deviation of the evaluation return computed over all trials that led to a winning test trajectory as a function of the ensemble average number of control steps (or equivalently, the number of collected samples). In all cases, SAC and DDPG display similar learning curves, while TD3 is systematically slower to reach the maximum return. This can be attributed to the inner workings of TD3, which updates the policy’s network only once every two updates of the value function, hence trading off accuracy of the policy gradient for number of gradient steps. This makes TD3 not suited to the present problem and possibly limits its usefulness for CFD problems where sample efficiency may be critical. Looking at DDPG and SAC, the variance in score of the latter decreases sooner across episodes in most cases (shaded area in Figure 11), meaning that the learning of SAC is more stable. Cases with a lot of dispersion with SAC correspond to especially hard situations, where DDPG and TD3 are anyhow less reliable. Despite these differences in the sample efficiency, the three algorithms reach the same asymptotic maximal average score in most cases. The few exceptions happen with TD3 where the algorithm would most probably reach the maximum return if the learning lasted more than 1,000 episodes. Note that the associated policies obtained from all trials discussed in this section lead to similar envelopes of optimal test trajectories (not shown), which suggests that all algorithms produce similar sets of optimal policies. Finally, Figure 11 shows a sharp decrease at the beginning of learning with SAC in most cases. While DDPG and TD3 only seek to maximize the objective term, SAC has an additional entropy term in the loss. Since Figure 11 does not show the full term which is maximized but only the objective term, it is possible to conjecture that the decrease at the beginning is due to the algorithm maximizing the entropy term at first.

Figure 11. Ensemble average and standard deviation of the evaluation return computed with only trials that led to a winning test trajectory (the average may therefore be performed with different number of runs for each algorithm). The curves were obtained by evaluating the policy at the end of each episode during training.

4.3. Variations in the reward model definition

Besides the raw usability of RL algorithms for fluid mechanics problems, and the hyperparameter sensitivity of such algorithms, one may wish to study how these methods cope with alternative definitions of the training objective, that is the reward model. Reward formulation can impact both the learning performance and the learned policy (Ng et al., Reference Ng, Harada and Russell1999; Henderson et al., Reference Henderson, Islam, Bachman, Pineau, Precup and Meger2018). In particular, two pitfalls must be avoided. First, the agent may find a workaround instead of actually solving the task (see, for instance, Rabault et al., Reference Rabault, Kuchta, Jensen, Reglade and Cerardi2019). Second, the agent may be very slow to find a solution, or even never find one, if the reward provides too little information (e.g., sparse rewards).

In this section, various formulations are investigated. These include a dense and/or a sparse component, as discussed in Section 3.1. The sparse part stands for the task that is actually solved (and on which the agent is tested), whereas the dense part guides the airfoil in the right direction. For all cases, the discount factor is set to $ \gamma =0.99 $ (see Section 3.1). As in Section 4.2, 20 learning trials are run for each set of parameters. The reward formulations and the results on reliability are summarized in Table 3, while the learning progress of selected cases is shown in Figure 12.

Table 3. Influence of the reward model (Eq. 8): cases characteristics and number of won trials over 20 realizations

Note. Number of won trials smaller or equal to 15 are highlighted in bold red, and in black bold when equal to 20.

Figure 12. Ensemble average and standard deviation of the evaluation return computed with only trials that led to a winning test trajectory (the average may therefore be performed with different number of runs for each algorithm). The curves were obtained by evaluating the policies at the end of each episode during training.

First, let us focus on cases Baseline (where $ {c}_2=10 $ ), $ {c}_2^0 $ , $ {c}_2^{100} $ , $ {c}_2^{1000} $ , and $ {R}^0 $ , leaving aside case $ {R}^{-1} $ . In case $ {R}^0 $ , a non-null reward is only obtained at the end of the episode if a won or lost situation is reached: the reward is purely sparse. On the opposite side, case $ {c}_2^0 $ features a reward signal that tells the algorithm, at each time step, whether the action has drawn the airfoil closer or farther away from the target; the reward is purely dense. In this last case, the reward model does not fit explicitly the criterion used for evaluating trajectories a posteriori (won when the plate enters the ball of radius $ \varepsilon $ around point B and lost otherwise). Conversely, the reward model encourages the algorithm to learn how to fly as close as possible to B. Cases Baseline, $ {c}_2^{100} $ and $ {c}_2^{1000} $ use a combination of both a dense and a sparse term. The three cases vary in the weight assigned to each term. Among all these cases, the most reliable option when using DDPG or TD3 is the Baseline one, where the sparse term value is approximately 10 to 25% of the sum of the rewards collected, at each time step, over an episode. On the contrary, among the less reliable cases for DDPG and TD3 is the fully sparse $ {R}^0 $ case. One should note that piecewise constant reward models and deterministic dynamics fall out of the scope of application of deterministic policy gradient algorithms (Silver et al., Reference Silver, Lever, Heess, Degris, Wierstra and Riedmiller2014). DDPG and TD3 belong to this class of algorithms so that they cannot be expected to perform well in sparse reward cases such as $ {R}^0 $ . Besides, even when theoretically applicable, sparse rewards often make the problem hard to solve since they induce a difficult uninformed exploration problem for the RL algorithm. In case $ {c}_2^{1000} $ , the sparse term is approximately 10 times larger than the sum of the rewards collected over an episode. The dense part might be too low for DDPG and TD3 to perform well. Overall, while DDPG and TD3 are strongly affected by the reward formulations discussed here, SAC is way more robust to these variations, owing to both the regularization term in its objective function and the better exploration process of the stochastic policy.

Finally, in case $ {R}^{-1} $ , a reward of $ -1 $ is applied at each time step, as well as a sparse term at the end of each episode. Such a formulation defines a signal that directly minimizes the time taken to reach B, but does not inform the behavior policy in terms of progress towards the goal: it boils down to the difficulties of a piecewise constant and sparse reward problem. Indeed, as long as the algorithm does not meet point B at least once by chance, given the data within the replay buffer, the most straightforward way to maximize the return is to minimize the number of time steps per episode by terminating the episode as fast as possible. Interestingly, the policy may therefore remain stuck in a local minimum corresponding to a “suicidal” policy (which terminates each episode as fast as possible), and never explore enough to discover there is a better way to maximize the return. While SAC is more reliable than DDPG or TD3 in all other cases, Table 3 evidences the very poor performance of all three algorithms for case $ {R}^{-1} $ .

Changing the reward model implies that different control problems are actually solved, with possibly different optimal trajectories for each scenario. Yet, as in Section 4.2, similar test paths are obtained whatever the reward model (not shown here). Consequently, this benchmark can be used to assess the sample efficiency in regard to the reward model applied. Despite the different scales, cases $ {c}_2^{1000} $ and $ {R}^0 $ display very similar learning curves in Figure 12. The dense term in the reward of case $ {c}_2^{1000} $ may be too small to be accurately taken into account in the learning so that this case is almost equivalent to the $ {R}^0 $ case. Regarding cases including both a dense and a sparse part, the smaller the sparse term, the quicker the maximum return value is reached (or equivalently, the better the sample efficiency). For all cases presented here, once again, similar learning curves are obtained for DDPG and SAC, whereas TD3 is slower to learn, likely due to its over-caution and delayed updates.

4.4. Variations in the dynamics model

The control problem considered until now has simple dynamics: the environment is deterministic and the task is permissive such that various policies can be quasi-optimal. Now, we investigate the behavior of algorithms in potentially more challenging situations, which could be encountered when moving towards real-world applications. Four different variations from the baseline are considered: decreased actuation frequency, increased airfoil density, added environment noise on the action, and, finally, replacing the low-order model with the resolution of the Navier–Stokes equations to account for unsteady aerodynamic forces. Such variations modify the dynamics of the problem and may therefore alter the learning performance and the associated policy. Both are analyzed in this section.

The first parameter investigated is the control frequency. In this work, this is treated as a tunable parameter. This may not be the case in physical problems that require a specific frequency to control a given phenomenon, or if the goal is to mimic a specific actuation device. Decreasing the actuation frequency leads to a lower-horizon control problem (fewer control steps in each episode) and could therefore induce easier convergence. However, this can hit a threshold frequency under which the problem is under-actuated and no successful policy can be found. The second parameter is the airfoil density. As the density increases, the airfoil inertia increases as well, and the reachable states are increasingly limited. One could then expect that with an increased density, the control problem to solve is harder because the policies have less latitude. Next, environmental noise is added to the action prescribed by the controller: the environment is not deterministic anymore. Note that unlike in robustness studies, here the noise is added during both the learning and testing phases so that the policy learns how to deal with noise. This added noise can resemble the uncertainty in the prescription of the action that can occur in an experimental setting. However, here, the noise value is selected randomly according to a Gaussian distribution centered in zero and of standard deviation $ \sigma \cdot {\dot{\beta}}_0 $ , where $ \sigma $ is given in Table 4. Note that this is different from the expected noise produced by experimental devices. Nevertheless, this allows a first assessment of the algorithms’ behavior in the presence of environment noise.

Table 4. Influence of the transition model: cases characteristics and number of won trials over 20 realizations

Note. Number of won trials smaller or equal to 15 are highlighted in red bold, and black bold when equal to 20.

Finally, the low-order model of the flying flat plate physics used until now is replaced by the resolution of the Navier–Stokes equations with the STAR-CCM+ solver, described in Section 2.1. This introduces more complex dynamics since added mass, virtual camber, and wake effects on the aerodynamic loads are captured, and therefore this case represents a step towards more realistic applications. The state components fed to the RL algorithm are kept unchanged from the low-order model case. As a consequence, the algorithm has now access to only partial information about the complete state describing the underlying MDP: the state observed by the algorithm does not fully describe the environment state, since it has no information on the whole flow fields. This situation leads to ambiguous situations where two identical measurements fed to the policy might actually correspond to two different environment states (e.g., same flat plate position, velocity, and pitch angle, but different surrounding flow field, such as a different unsteady wake topology). Such a control problem moves away from the MDP theoretical framework on which the three analyzed algorithms are built, but also better represents most real-world problems. While the training lasted 1,000 episodes with the low-order model, here, only 500 episodes were performed before testing the controller. This is sufficient to obtain a converged solution and allows to reduce the computational cost of the trials with the high-fidelity model. Besides, a well-known unlearning phenomenon (where the learned policy forgets after being able to perform the task) was observed in a few cases when running the training for longer: a detailed analysis of this phenomenon is beyond the scope of this current work. For all these four variations from the Baseline, Table 4 summarizes the values of the parameters investigated and the corresponding learning reliability, while Figure 13 shows the learning history of some selected cases.

Figure 13. Ensemble average and standard deviation of the evaluation return computed with only trials that led to a winning test trajectory (the average may therefore be performed with different number of runs for each algorithm). The curves were obtained by evaluating the policy at the end of each episode during training.

In case $ {f}^{1/2} $ (respectively, $ {f}^{1/5} $ ), the actuation frequency is decreased compared to the Baseline case: the physical time step $ dt $ is kept constant, while the control time step $ \Delta t $ is increased so that $ \Delta t=2 dt $ (respectively, $ \Delta t=5 dt $ ). As a result, for the same number of episodes, the number of collected samples and neural network updates would be decreased by 2 (respectively, 5). To keep these constant across cases, the number of episodes is increased accordingly so that the total number of collected samples remains of the same order of magnitude. However, this increases the cost associated with sample collection, and, incidentally, with each training trial. Results from Table 4 show that decreasing the control frequency has no significant impact on reliability. However, comparing the curves obtained for the Baseline and the low-frequency actuation case in Figure 13 reveals that a faster learning is achieved for the latter. Indeed, case $ {f}^{1/5} $ is a lower-horizon control problem: the average number of control steps between the initial state and the end of an episode is five times smaller. As a result, the sparse earning at the end of a winning episode is faster to propagate towards initial states, resulting in a faster convergence.

To ensure that the global trend obtained when assessing the reliability of the RL algorithms is not specific to the chosen setup, an additional flat plate configuration is tested: case $ {d}^{150} $ corresponds to a different airfoil density compared with Baseline, whose density is $ {d}_s=30\mathrm{kg}\;{\mathrm{m}}^{-3} $ . DDPG and TD3’s reliability slightly decreases when increasing the flat plate density. This is consistent with the fact that the problem is more challenging. Besides, the sample efficiency of DDPG and SAC is not significantly affected, while TD3 requires more samples to converge toward the maximum value of the return.

In cases $ {\sigma}^{0.1} $ to $ {\sigma}^{0.75} $ , where environmental noise is added to the action prescribed by the agent, all three algorithms remain reliable whatever the noise intensity. Because the actuation frequency is high, and the noise is random, its average over time equals zero, which yields a null global effect. A correlated noise may lead to different results. Here, the sample efficiency remains similar to that of the Baseline for all three algorithms.

Finally, in the case of CFD, whatever the algorithm, reliability stays high and sample efficiency is slightly higher, compared to the Baseline. When using RL algorithms in conjunction with CFD solvers, computing the samples is expensive in terms of CPU usage. A trial with the low-order model lasts less than an hour on a single computing core. With the CFD environment, the total computational cost required to perform a trial (RL training and sample collection) is 1,200 times larger. This corresponds to the large CPU time required to collect samples when coupling CFD and RL algorithms, the cost of the RL part being negligible. This significant increase exemplifies the need for reliable and sample-efficient off-policy RL algorithms, which is at the core of the present study.

While the variations in the dynamics model considered here have a limited impact regarding the learning performance, their influence on policies is evidently high. Figure 14a shows all the winning test trajectories and Figure 14b shows the corresponding controller actions, obtained with SAC, for each case presented in Figure 13. Moreover, for each case, in an attempt to compare visually the policies obtained by SAC, TD3, and DDPG, we count how many times each action is taken along a trajectory, and average this count over 20 trials. This provides a probability density function (PDF) of actions along trajectories which is shown in Figure 14c. The trajectories and action signals of the Baseline case obtained with SAC have already been discussed on Figure 10: policies vary a lot from one test to another, but most policies consist in applying a negative pitch rate at the beginning, then stabilizing on average around a value close to zero. Finally, near the end of the episode, close to point B, the different policies cover the whole range of possible values. This leads to the PDF of the actions observed in Figure 14c for SAC, which peaks slightly below zero. The PDFs obtained for DDPG and TD3 are, however, quite different. While a peak below zero can still be observed for DDPG, TD3 evidences a much more evenly distributed pitch rate over the range of possible values. Besides, DDPG and TD3 show a small peak of the PDF at the extrema, which is not observed for SAC. This evidences some periods of action saturation with the former algorithms that are not experienced with SAC. Variations of the PDF are strong between the algorithms for this Baseline case. On the contrary, all three algorithms present similar PDFs for the low frequency actuation case $ {f}^{1/5} $ . Apart from the peaks at the extrema observed only for DDPG and TD3, the PDF curves overlap quite well. Because of the lower control frequency, this case leads to a graphically rougher action sequence signal than the Baseline case but the trajectory envelope is similar. Then, when the plate density is increased, the controller has much less room to maneuver the possible actions at the beginning of the episode. All policies able to reach B consist of applying a significant upward pitch rate from the very beginning until approximately half of the episode. Otherwise, the strong inertia of this case prevents the airfoil from reaching the target. The associated PDF evidences a strong peak close to the minimal pitch rate allowed by the problem. The test trajectories resulting from these resembling policies are very similar. Next, case $ {\sigma}^{0.75} $ depicts a very noisy action sequence due to the additional environment noise. However, the actuation frequency is so high that the trajectories are barely affected by the noise and are very similar to that of the Baseline. The three PDF shapes differ from one algorithm to another but they all appear flattened.

Figure 14. (a) Trajectory, (b) normalized action $ \dot{\beta}/{\dot{\beta}}_0 $ , and (c) normalized action probability density function (PDF) for selected cases.

Finally, the action sequences in the CFD case are more saturated than in the Baseline case and do not evidence the preference for a negative pitch rate at the beginning of the episode, as observed in all other cases. Instead, the pitch rate values spread over the whole range at the beginning of the episode. The action PDFs are even more flattened than in the noisy case, and do not evidence any pattern apart from the peaks at the extrema already observed for DDPG. Corresponding trajectories are less smooth and way more spread than the ones obtained with the low-order model. This is most likely associated with the richer dynamics involved and the more complex feedback coupling between the airfoil and the fluid when the flow is solved with the Navier–Stokes equations. In particular, because the airfoil can experience large accelerations in pitch, added mass forces (which are not included in the low-order model) are one to two orders of magnitude larger than quasi-steady forces induced by translation, and the controller relies on such forces (rather than quasi-steady forces) to drive the airfoil to point B. In fact, although not shown here for the sake of conciseness, the effective angle of attack can be mostly negative throughout the maneuver, hence mostly producing negative quasi-steady lift. Rather, gain and loss in altitude are correlated with rapid pitch-up and pitch-down (i.e., angular accelerations) of the airfoil. The control strategy is therefore fundamentally different from that obtained with the simplified model, which most likely explains the differences observed in Figure 14. We add that virtual camber and wake effects may also play a role in the control strategy, but their order of magnitude is significantly lower than that of added mass effects. An example of the vorticity field that develops around the airfoil during the maneuver is provided in Appendix A.

4.5. Summary on the sensitivity analysis and reasons to choose SAC over DDPG or TD3

The overall performance of DDPG, TD3, and SAC on the first task in terms of reliability is summarized in Figure 15. For each set of parameters tested in the previous sections, and detailed in Tables 2–4, the percentage of won trials is reported in green as a positive number, and the percentage of lost trials is reported in red as a negative number.

Figure 15. For all sets of parameters reported in Tables 2–4 the percentage of won trials is reported in green as a positive number and the percentage of lost trials is reported in red as a negative number.

Figure 15 highlights clearly that SAC is more reliable than DDPG and TD3, which both evidence strong reliability variations when parameters vary. On the other hand, SAC is less sensitive to parameter tuning and especially “hard” problems for DRL algorithms (e.g., small network, low replay buffer capacity). The only exception is the $ {R}^{-1} $ case (see Table 3), which leads to catastrophic performance for all three algorithms.

On this benchmark, and without special optimization of any of the algorithms, SAC was observed to be slightly more CPU consuming than DDPG or TD3 (which is consistent with its inner adjustment mechanisms), but the computational cost gap between algorithms goes to zero when the physics is solved with CFD. Indeed, in such a case, the cost associated to the RL algorithm is negligible compared to the cost associated to the CFD computation.

To sum up, on the first task at hand here, SAC is more reliable, has the same cost when used with CFD, and has a sample efficiency comparable to DDPG and superior than TD3. It is therefore more suitable to set up this new control problem for which we had no prior knowledge in terms of parametrization. The knowledge acquired on this simple task can then be used as a basis when moving towards more complexity.

Interpreting why SAC seems more robust and efficient is an open question in the RL literature, for which we provide some elements here. First, and maybe overall, the entropy regularization term in SAC is known to make the optimization landscape smoother and to drive the optimization more consistently towards the same minima (Geist et al., Reference Geist, Scherrer and Pietquin2019). In the present case, this is reflected through the consistent performance of SAC. Additionally, SAC provides a more principled way of performing exploration, as its behavior policy benefits from an adaptive temperature parameter, while DDPG and TD3 rely on adding Gaussian noise to a deterministic policy. This induces a possibly better coverage of the relevant states and actions, and, in turn, better gradients for the value function and policy. Finally, SAC has been recently analyzed under the lens of robust MDPs (Eysenbach and Levine, Reference Eysenbach and Levine2022) where the policy is optimized to perform well across a span of control problems (instead of a single one). The robustness of SAC policies to variations in the control problem dynamics makes it less sensitive to both environment identification errors and hyperparameter tuning. Nevertheless, SAC and TD3 are mathematically very close (both are approximate value iteration algorithms, using a stochastic behavior policy) and claiming pure dominance of an algorithm over another in RL is a delicate endeavor (Henderson et al., Reference Henderson, Islam, Bachman, Pineau, Precup and Meger2018). Consequently, the take-away message of this section should not be that SAC will always dominate over TD3 or even DDPG. It is actually very likely that there exists a set of hyperparameters for which TD3 or DDPG performs as well or even better than SAC. However, in the present case—and possibly in other fluid dynamics cases such as those studied by Bucci et al. (Reference Bucci, Semeraro, Allauzen, Wisniewski, Cordier and Mathelin2019), Novati et al. (Reference Novati, Mahadevan and Koumoutsakos2019), and Rabault et al. (Reference Rabault, Kuchta, Jensen, Reglade and Cerardi2019) for instance—SAC appears less sensitive to hyperparameter tuning and is thus a good candidate for a fluid dynamics practitioner wishing to set up a new case study.