Effective control of two-dimensional Rayleigh–Bénard convection: Invariant multi-agent reinforcement learning is all you need

Deep reinforcement learning (DRL) has recently emerged as a promising methodology for performing active flow control¹ (AFC). DRL is a subclass of the much more general field of RL, where deep artificial neural networks are used to parameterize the value function in Q-learning methods² or serve to parameterize the policy or the mapping from state to actions in policy gradient methods.^3,4 The DRL “controller,” which includes the control law to be applied on the system, and algorithms around it to allow tuning of the control law and query of the control law, is called the DRL agent. The general principle of the (D)RL paradigm⁵ is based on considering the system to control as a black box (the “plant” in classical control theory or the “environment” in DRL terms) and to interact with it through just three channels of communication. At a given time step t, (i) the DRL agent can query the environment for a (possibly noisy, stochastic) state estimate s_t; (ii) the DRL agent can provide to the environment a control instruction called the action in DRL terms, a_t; and (iii) the environment can provide back a reward r_t that indicates the quality of the current state of the system. Based on these three channels of interaction, the DRL paradigm aims at defining effective algorithms to learn control strategies that maximize the combined expected rewards at all future times from the given state at time t, also known as the expected return $G t = ∑ k = 0 ∞ r t + k + 1 γ k$⁠, where γ is the discount factor, typically $γ ≈ 0.95 − 0.99$⁠. In order to make the computation of R_t tractable, DRL interaction with the environment takes place on an episode basis during training. An episode is a simulation with a given time extent (typically covering a few timescales of the slowest dynamics of the system) and corresponds to one individual trajectory of the system in its phase space. The optimization of the control applied by the DRL agent is performed through trial-and-error learning and direct online interaction with the environment, a process that involves stochastic exploration of the environment properties (referred to as the “training mode” of the DRL agent). This allows to explore and map the behavior of the environment and to determine the best trajectory that the environment should be actuated to follow in the phase space to optimize the reward. A number of algorithms have been developed to perform this optimization, the most famous belonging to either (i) the Q-learning paradigm^6,7 based on the Bellman equation; (ii) the policy-gradient method,⁸ which relies on a direct estimate of the sensitivity of the reward to the policy parameters; and (iii) a combination of both paradigms, such as the actor-critic proximal-policy optimization (PPO) method,⁹ which combines their respective strengths. Once the trial-and-error learning process has converged and the DRL agent is fully trained, the trial-and-error exploration component can be turned off, so that a deterministic controller is obtained (corresponding to the “deterministic mode” of the DRL agent).

This DRL framework is very attractive for the control of complex, high-dimensional, non-linear systems such as those obtained in AFC. Indeed, DRL does not make any assumptions on the system to control; actual DRL algorithms have been increasingly refined in the course of the last 10 years and are well adapted to solve non-linear, high-dimensional control problems from games,^7,10,11 industrial control problems,¹² or classical physical systems such as those present in fundamental fluid mechanics.^13,14 Compared with traditional approaches, such as the adjoint method,¹⁵ DRL does not need direct insight into the full state of the system and its governing equations. Moreover, while DRL training is expensive in terms of the number of interactions with the environment needed for the trial-and-error learning to take place, once the DRL agent is trained, predicting the next action to take given the current state is very cheap computationally. These aspects make DRL an attractive method for real-world control applications, including AFC.

Over the last 5 years, DRL for AFC has been applied to increasingly complex cases. Starting from simple 2D cylinder cases¹⁶ at low Reynolds number (⁠ $R e = U L / ν$⁠, which measures the relative importance of inertia and viscosity, where U is the characteristic velocity scale of the system, L its characteristic length scale, and ν the kinematic viscosity), DRL has further been used for investigating a number of different configurations involving 2D cylinders,^17–20 looking into increasingly complex stochastic wake control problems as Re is increased,^21,22 and recently, considering the control of 3D channel flows.^13,14 DRL has also been applied to a number of other AFC cases, from 2D Rayleigh–Bénard convection in a small-size channel,²³ to instabilities developing in thin fluid films,²⁴ engineering cases related to wings,^25,26 and other practical problems.^27–29 DRL has also been combined with flow-stability analysis to guide the placement of observation probes²⁰ for effective control. In this case, a stability analysis was used to reveal the structural sensitivity that determines the location of the origin of the instability, i.e., the “wavemaker” region. The structural sensitivity can also be computed using data-driven techniques.³⁰ Positioning the observation probes in the region of the wavemaker has shown to achieve the best possible learning curve with the reward maximization then taking place in the least number of episodes, compared with other configurations of observation probes.²⁰ While these studies have mostly considered numerical test cases, demonstrations of the applicability of DRL for real-world AFC have been presented in test experiments.³¹ 

Moreover, a number of peripheral challenges associated with discovering efficient flow control strategies have been solved over the last few years. Rabault and Kuhnle³² demonstrated how the multi-environment paradigm can be leveraged to obtain speedups in training speed proportional to the number of environments used. In particular, they showed that DRL training speed is directly proportional to the number of (state, action, and reward) triplets generated per wall clock time unit, and that this is in turn directly proportional to the number of CFD simulations used in parallel to train the DRL algorithm and feed it with data streams. Belus et al.²⁴ have demonstrated that the curse of dimensionality on the control dimension can be effectively solved by leveraging a multi-agent reinforcement-learning (MARL) approach (although the authors did not refer to MARL under this name at the time), in cases where some general properties of the system are invariant in space. In particular, MARL was demonstrated to be an effective way to control systems containing several generally similar structures that are far enough to be mostly independent from each other, by sharing the knowledge acquired among them. This, in turn, allows to reduce the combinatorial cost and alleviate the curse of dimensionality associated with the exploration of large action spaces.

Following these advances, DRL can now be used to control different physics as the effect of non-linearity increases,²² to perform effective control over a range of flow conditions in simple AFC test cases,³³ or to control complex three-dimensional (3D) turbulent channel flows.¹⁴ DRL has also been demonstrated in a number of test situations involving chaotic flows, such as the one arising from the one-dimensional (1D) Kuramoto–Sivashinsky system of equations.^34,35 The placement and optimization of the sensor inputs used to drive DRL algorithms has also been investigated, see Refs. 36 and 37. As the literature is growing fast, we refer the reader curious of a more detailed overview of DRL applications to fluid mechanics in general and AFC in particular to one of the several reviews that have recently been provided on the topic.^38–40 Similarly, DRL is becoming an increasingly common topic in AFC, and the DRL methodology per se has now been discussed in a number of AFC studies, so we refer the reader curious of specific DRL algorithmic details and refinements to the corresponding body of literature.⁴¹ As a consequence of this increasing popularity, DRL frameworks are starting to emerge, both in the DRL^42–44 and in the fluid-mechanics communities,^45–47 a fact that now provides frameworks for effectively coupling DRL with well-established computational-fluid-dynamics (CFD) tools. This makes it increasingly easier to apply DRL to AFC problems.

At this point, we want to emphasize that, while DRL is a promising method given that it has provided state-of-the-art results in a variety of fields as previously discussed, it is not the only method that can perform AFC by considering the flow to control as a black box to be optimized. Machine-learning control (MLC^48–50) has, indeed, been a topic of discussion for quite a long time, and the DRL method is, in this aspect, a newcomer in this field. Traditionally, methods such as genetic programing (GP^51,52) or Bayesian optimization (BO⁵³) have also been applied with great success to a variety of practical problems^54–57 and can be competitive with DRL in a number of applications.⁵⁸ However, we focus our efforts on DRL in the present study, as DRL is showing great flexibility in constraining the admissible structure of the control law (see, e.g., the discussions around invariants and MARL in Ref. 24). Moreover, the DRL methodology is promising thanks to its ability to use local (in space, but also in time) information to perform optimization. Compared to, e.g., a direct application of GP or BO, which only optimize based on the mean efficiency of a control law over one simulation containing several typical physical time scales of the underlying dynamics, and which do not exploit information about the high-frequency quality of the control law, DRL gets a reward at each interaction step. Thus, DRL can learn to both reproduce positive sequences of actions and avoid negative sequences of actions happening within a single episode. As a consequence, DRL is able to generate information about the quality of the control law at a higher time granularity than the episode itself and generates a larger volume of usable reward information than would be provided by monitoring for example the averaged reward over an episode.

In the present work, we use DRL techniques to discover effective strategies to control the instabilities and the heat flux in a two-dimensional (2D) Rayleigh–Bénard convection (RBC) problem. The RBC phenomenon is intrinsic to many industrial applications and can be found in a wide range of natural phenomena.⁵⁹ The canonical initial condition in a RBC problem is a fluid at rest that is being heated from a lower wall and/or cooled from an upper wall (see Fig. 1). Because of the induced temperature gradients and due to the fluid-buoyancy effects, natural convection occurs. The dynamical state of the RBC system can be characterized by a non-dimensional parameter, the Rayleigh number Ra, which is a ratio of the timescale of heat transport by diffusion or conduction to the timescale of heat transport by buoyancy-induced convection. Beyond a critical value of Ra = Ra_c when convective effects dominate, RB instability occurs through a supercritical bifurcation from an initial quiescent state, to time-dependent cellular motion.⁶⁰ 

The convective cells cause a non-spatially uniform heat transfer along the domain, as well as large-scale fluid motions and structures containing significant kinetic energy, which can be detrimental in many industrial applications. Thus, control of RBC is crucial to maximize the efficiency and quality in industrial processes or systems in various sectors. The control might involve, for example, modulating the spatial distribution of the bottom-plate heating in the canonical RBC configuration. Optimizing such a control is challenging topic for classical control-theory methods. Indeed, many classical linear-control techniques such as the adjoint method have only a limited domain of applicability, owing to their local, linear nature, which limits their ability to control non-linear phenomena such as RBC.

Passive flow control has been extensively used with minor improvement in the Ra_c for which instabilities start.^61–65 On the other hand, AFC has explored more sophisticated control strategies that manage to keep a stable convective cell at considerably higher Rayleigh numbers than Ra_c. The most common way of carrying out this AFC has been with the use of jet controllers, performing small perturbations of the velocity and/or temperature, or with a linear control of the temperature in the heated boundary layer.^66–80 In configurations where the RB instability is too strong to be completely eliminated, a reduction in the strength of the RBC cells and associated flow motion, as well as kinetic energy, are also key goals.²³ Recently, the use of DRL to control the RBC has been employed for the first time in Ref. 23. Comparing the results with the linear control employed in Ref. 80, the Ra_c beyond which instabilities start to emerge in the specific configuration considered was increased from 10⁴ to $3 × 10 4$⁠, and even in the case where RBC cells were present, their intensity was strongly reduced.

In the present work, we extend the work in Ref. 23 (which we will refer to in the rest of this paper as GB for “Gerben Beintema,” from the name of the lead author), to look further into DRL control of RBC. In particular, we improve on the following aspects:

• The domain used by GB is bounded by adiabatic walls on the right and left ends, and the domain itself is relatively narrow. This in turn constrains the RBC cells obtained. Moreover, the walls presented in GB's setup can be used by the DRL agent to help move around and break the RBC cells, as they provide strong flow confinement. By contrast, we use a wider domain of height H = 2 and width $L = 2 π$ (compared with the H = 1 and L = 1 used in GB) with periodic conditions on the right and left boundaries. In other words, our domain has an aspect ratio of π compared with the aspect ratio of 1 used in GB. As a consequence, our DRL agent has to control the RBC instabilities without the possibility to help itself with the walls, and the RBC cells are closer to unconstrained cells obtained in wide domains relevant for industrial use.

• One limitation in GB's work is how DRL control is applied to the RBC simulation. In order to control the RBC cells, GB modulates the temperature at the lower wall, keeping its mean temperature constant. This, naturally, means that multiple outputs, i.e., temperature perturbations on bottom-wall segments, need to be applied to the N different control segments at the bottom wall. GB formulates this control problem by simply asking a single DRL agent to generate N outputs. As discussed in Ref. 24, this approach results in the curse of dimensionality on the control space, which makes learning challenging. By contrast, we apply the invariant MARL approach suggested by Ref. 24, which allows to control an arbitrary number of heating segments without resulting in the curse of dimensionality. This has a considerable impact on our ability to rapidly learn an effective control of the RBC system.

• GB uses a lattice Boltzmann method (LBM), parallelized to run on a specific graphics-processing unit (GPU), to help accelerate the simulations. In such a way, GB partially alleviated the curse of dimensionality introduced by controlling a multiple-output DRL system without using a MARL approach. Indeed, the use of a GPU-parallelized LBM allows GB to run a very large number of DRL trial-and-error steps and, therefore, to partially compensate for the very slow and challenging learning arising due to the curse of dimensionality. Note, however, that this is only applicable for narrow channels with a limited number of control values, where the curse of dimensionality remains under control. Moreover, LBM is effective, especially as they are well adapted to the parallel nature of GPUs, but understanding their accuracy can be challenging. Finally, GB's implementation is not publicly released and may be hardware dependent. By contrast, we use a spectral Galerkin solver running on standard central-processing Units (CPUs). The solver is both open source, massively parallelizable across many CPU cores, has high order and well-understood accuracy (see Sec. II for more information). Moreover, we release all our code as open-source materials (see  Appendix A). Our code can, therefore, be used to further evaluate a wide range of domain sizes and non-dimensional parameters in future studies.

The structure of our study is as follows. In Sec. II, we describe the formulation of the RBC problem including the geometry, equations, and boundary conditions adopted (Sec. II. A.). We then describe the numerical method used to simulate the CFD (Sec. II B.), followed by the details of the DRL-based control (Sec. II C.). In Sec. III, we proceed to present the results. We compare the RBC obtained without control (baseline) to the case with control following a non-MARL method similar to GB and also to the case with control following a MARL approach. We discuss learning speed and robustness. We then show that the invariant MARL controller is able to effectively coalesce RBC cells in the case studied. Finally, we discuss this result and its implication for future works and industrial applications in the conclusions in Sec. IV.

During the entire study, the Prandtl number employed is Pr = 0.71, which corresponds to the Prandtl number of air. As described in the introduction (Sec. I), Ra is the ratio of the timescales associated with the thermal transport due to diffusion and convection. This implies that flows with a high Ra are more prone to instabilities due to buoyancy-driven convection. The Ra used in this work is $R a = 10 4$⁠, which is above the critical Ra_c when no control is applied in the flow.

The three scalar equations (5), (1a), and (1c) are implemented using a highly accurate spectral Galerkin⁸³ discretization in space and a third-order implicit/explicit (IMEX) Runge–Kutta method⁸⁴ for the temporal integration. The Galerkin method makes use of tensor product basis functions constructed from Chebyshev polynomials for the wall-normal direction and Fourier exponentials for the periodic direction. The boundary conditions in both directions are built into the basis functions and as such enforced exactly. For the wall-normal direction, this requires the use of composite Chebyshev polynomials,⁸³ whereas the Fourier exponentials for the x-direction are already periodic. Since the continuity equation cannot be used to find u for Fourier wavenumber 0, we solve for this wavenumber the momentum equation in the x direction. All other unknowns are closed through Eqs. (5), (1a), and (1c). The convection terms H and $u · ∇ T$ are computed in physical space after expanding the number of collocation points by a factor of 3/2 in order to avoid aliasing. For H, we use the rotational form $H = − u × ( ∇ × u )$⁠, with the remaining $1 / 2 ∇ u · u$ absorbed by the pressure, and for temperature, we use the divergence form $u · ∇ T = ∇ · u T$⁠.

The code is implemented using the open-source spectral Galerkin framework “shenfun,”⁸⁵ where equations can be automatically discretized through a high-level scripting language closely resembling the Mathematics. The Navier–Stokes solver has been verified by reproducing the growth of the most unstable eigenmode of the Orr–Sommerfeld equations over long time integrations. The Navier–Stokes and Rayleigh–Bénard solvers are distributed as part of the shenfun software, and there is a demonstration guide published in the documentation,⁸⁶ which also provides a much more detailed description of the numerical method. Links are provided in  Appendix A, and we refer the reader curious of all the technicalities of the CFD solver implementation to the resources detailed therein.

The observation probes are distributed over the domain as a uniform probe-mesh as can be seen in Fig. 1. There are 32 × 8 probe-mesh points in the x and y directions, respectively. The number of quadrature points in the solution approximation used by the spectral Galerkin solver is 64 × 96, with Fourier and Chebyshev basis functions in the x and y directions, respectively. The resulting solution is thus spatially continuous, and the observations are then evaluated from these global, continuous spectral Galerkin functions at the uniform-mesh probe locations. These observations are used in the DRL-based control methodology which we describe next.

Using the numerical methods described in the previous section, the flow is solved starting from an initial condition of a constant adverse temperature gradient in the vertical direction. The simulation is performed from this state, leading to the RB instability to occur, which develops into a two-cell configuration that becomes ultimately stationary. No control is applied up to this point, and this constitutes a “baseline” starting state upon which control can be applied. The baseline is saved and then loaded at the start of each episode in the learning phase. The DRL control is performed using the “Tensorforce”⁴² framework, from which we use the PPO algorithm.

As mentioned in the Sec. I, the interaction between the agent and the RBC environment occurs through three fundamental channels of communication: the state or observation s_t, reward r_t and action a_t, where the subscript t is the time step of the CFD. The state s_t consists of sets of time-averaged velocity and temperature $( u ̃ t , v ̃ t , T ̃ t )$ measured by computational “probes” distributed across the domain in a uniform mesh as seen in Fig. 1. The averaging is performed over the previous 4 time steps (where the tilde over the values above represents this average). The observations are flattened into one single-dimensional array and then fed to the agent's neural network input layer. The way in which the observations are passed to the agent differs in the single-agent and multi-agent frameworks, and will be described in the following.

Finally, to avoid nonphysical sharp discontinuities in the spatial variation of actuations across the segments, the left and right $10 %$ of each segment are smoothed relatively to their immediate neighbors following a cubic transition function. In such a way, the cubic function transitions from $T t ( j )$ to $T t ( j + 1 )$ in a continuous and differentiable way. Figure 2 shows a schematic representation of how the temperature control is applied as a whole. An illustration of the cubic transition is shown in the zoomed Fig. 2(b).

While the general DRL setup we just described is perfectly valid from a DRL point of view, and it is actually implemented and referred to as the “single-agent-reinforcement-learning” (SARL) setup in the following, as illustrated in Fig. 3, it does not exploit the intrinsic structure of the RBC problem. Indeed, the physical mechanism driving RBC does not depend per se on the position along the wall direction (i.e., along the x axis). This, combined with the periodic boundary condition we implement in the CFD solver on the left and right ends of the domain, implies that the DRL-and-RBC setup is globally invariant along the x direction. This invariance can, therefore, be exploited through implementing a MARL approach, as detailed in Ref. 24, and previously discussed.

In the following, we proceed to define the “MARL-DRL” setup based on modified versions of the states, rewards, and actions described thus far. We first regard the N segments introduced previously as separate “pseudo-environments,” each with their own independent state, action, and reward channels of communication with the PPO agent. Since our RBC setting is invariant in the x direction, each trajectory in the phase space obtained for each pseudo-environment provides knowledge about the dynamics to control to the DRL agent. As discussed in Ref. 24, this can drastically speed up and improve the quality of the training and learning.

We make a note here regarding the number of pseudo-environments N, which is used is a metaparameter of the MARL setup. Typically, N needs to be chosen so that the following needs are balanced:

1. The segments should be small enough to ensure sufficient spatial granularity in the control signal to be able to actually control the flow. As a consequence, there should be at least a few control segments per spatial feature of the flow. In our baseline case, given that there are 2 convection cells of length scale $≈ π$⁠, we should have at least a few control segments per π distance along the x axis.

2. Despite the fact that MARL is efficient at learning even in the case of many control segments,²⁴ the number of segments should not be too large to overcome certain practical limitations:

3. A large number of segments will make the control output discontinuous as a large number of points, requiring the use of more aggressive smoothing functions.

4. Each pseudo-environment corresponding to each segment will require its own stream of data to be passed to and received from the agent. From a computational viewpoint, this increases the cost of the DRL deployment.

5. Neural networks tend to learn best when learning from uncorrelated data. This means that, ideally, feeding similar data many times to the DRL agent should be avoided. Therefore, providing segments large enough so that the view on each of them is significantly different from the views on the others is beneficial for the quality of the learning.

6. Segments that are too small in size can produce an overall temperature profile (which is the control output), which is highly non-uniform. As a consequence, a Fourier decomposition of the corresponding temperature profile in space will consist of elements with a high wavenumber energy content. These high wavenumber energy content elements can make it more challenging to obtain good stability of the CFD simulation.

Considering the points above, we choose a value of N = 10 (Table I) as a reasonable trade-off. We note that this is the analogous (in space) to the considerations that are presented in Ref. 32 when choosing the action duration and episode length of the DRL algorithm (i.e., in time).

In the following, we resume the description of how the states, actions, and rewards communication channels are modified in the MARL-DRL setup. The state s_t of each pseudo-environment is defined as the global set of observations whose order has been rearranged in such a way that the observations from the probes directly above its respective control segment are in the center of the observation list. This “recentering” is further elaborated in the description of the schematic of Fig. 4.

A schematic representation of the MARL-based control methodology is presented in Fig. 4. In order to make the figure easier to read, only three representative pseudo-environments are shown there instead of N = 10. Each pseudo-environment corresponds to a localized region of the physical domain. Note that for ease of explanation, in this plot, we deviate from the notation convention used in the rest of the paper: the subscripts for s, a, and r are now the environment index j = 1, 2, and 3. We define s_j as the values of the observation from the probes directly above each control segment, r_j as the total reward from each pseudo-environment [Eq. (10)] and a_j as the control actions applied by each MARL pseudo-environment. The top diagram represents the control methodology. At the start of the cycle, the a_j are merged and communicated to the first pseudo-environment (j = 1), which is the sole pseudo-environment that runs the whole CFD simulation, while the others pseudo-environments wait until CFD simulation is finished. The CFD simulation is run for a time equal to the action duration (Table I). Once the action is finished, the entire set of observations $s = ( s 1 , s 2 , s 3 )$ is retrieved from the observation probes and communicated to all the pseudo-environments. Then, in the “Merge results” block, each pseudo-environment executes the “recenter()” function. In recenter(), a rearranging of the order of s_j is performed in each pseudo-environment, so that the s_j corresponding to pseudo-environment j is moved to the middle of the observation vector. We denote the recentered observation vectors as $s ′ j$⁠. The r_j is obtained from the $N u loc$ and $N u inst$ computed using the local s_j and the global s, respectively. This is done using Eq. (10), and it is represented for brevity as a function of g in Fig. 4. The $s ′ j$ and r_j are then passed to the respective MARL agents, which output the actions a_j and the cycle is repeated again. The MARL agents are actually implemented as data streams to a repeated DRL agent, which means that all MARL agents share the exact same policy and neural networks. This is how the control invariance across the domain is implemented in the MARL case, i.e., by re-using the exact same agent multiple times across the domain, similar to what is described in Ref. 24.

Having described the setup of SARL and MARL, we now proceed to the training process. The main parameters used for the DRL setup are summarized in Table I, and the hyperparameters used for the definition and training of the agent are in Table II. A total of 350 CFD episodes are performed for the training run in SARL, with 200 actuation updates per episode. Each actuation consists of a set of N actions or temperature control values for the N segments. In SARL, the policy update is performed every “batch size” number of episodes, which in our case is set to be 20. In MARL, we make a distinction between “CFD episodes” and “MARL episodes.” The “CFD episode” is when the CFD simulation is run for 200 actuations, as in the SARL case. Each actuation comprises N actions applied locally by each MARL pseudo-environment. Thus, each MARL pseudo-environment, in the course of one CFD episode, applies 200 control actions, which we define as a single MARL episode. Given that there are N MARL pseudo-environments, there are therefore 350  $× N$ MARL episodes in one training run for MARL. The batch size for MARL is 20 MARL episodes, or 20  $/ N$ CFD episodes. In the discussion of the results section, in the context of MARL, the word “episode” refers to “CFD episode.”

Training in both SARL and MARL is performed with a discount factor of $γ = 0.99$⁠. Wall-clock training time for 350 episodes in SARL and 350 CFD episodes in MARL is around 34 h using a single core on a modern workstation (the wall clock duration is approximately the same in both cases, as CFD is the dominating computational cost). The duration of each episode in CFD time units is action duration × number of actions per episode $= 300$ time units (Table I). This time is enough to capture the dynamics of the flow in the transition from the baseline state to the optimized controlled flow. Within the duration of one action (i.e., within a single control step), the control temperature applied to each control segment is established by the DRL agent through Eqs. (8) and (9). This control step is chosen to be short enough so that the agent can actuate under fast changes of the state of the flow, but also long enough so that the flow can develop according to the actuation performed, similar to what is recommended in Rabault et al.¹⁶ 

In Sec. III, we analyze and discuss the results obtained. In the first place, at the start of the baseline convergence simulation, three convective cells appear. However, this is an unstable state that lasts for just a few time units and the final baseline, which includes a two-cell configuration, is achieved as visible in Fig. 5 (Multimedia view) (see also the corresponding multimedia data from  Appendix C). A full video of the baseline simulation can be seen in the link provided in  Appendix C.

Figure 6 shows the evolution of the Nusselt number and the kinetic energy along the duration of the baseline convergence simulation. The flow reaches a stable state after approximately 400 time units, where the Nusselt number converges to a value of Nu = 2.68.

In the baseline flow, two convective cells are present, with a total of two counter-rotating vortex pairs. This is different from the baseline layout obtained in GB, where a single rotating vortex was obtained. The differences in flow topology with respect to the results by GB come from the adiabatic lateral walls and the aspect ratio of 1 used in GB, compared with the periodicity in x and an aspect ratio of π used in the present work. Therefore, we can expect different control strategies compared to those obtained in GB.

We now present the results obtained from DRL control of the Rayleigh–Bénard system. To compare the learning in the SARL and MARL frameworks, we plot both their learning curves in Fig. 7 (Multimedia view) (see also the corresponding multimedia data from  Appendix C). The learning was performed for 350 CFD episodes of 200 actuations each as previously highlighted, starting from the final state of the baseline (Fig. 5, see also the corresponding multimedia data from  Appendix C). This means that the same CFD duration is used to train both the SARL and MARL agents. Since CFD simulation is by far the most time- and resource-intensive part of the training process (similar to what was reported in Ref. 16), this means that the wall-clock time and resources used for performing the SARL and MARL trainings are about equal. The time-averaged Nu at the end of the baseline is plotted for reference, corresponding to a value of 2.68 in Fig. 7 (see also the corresponding multimedia data from  Appendix C). We also show the time-averaged Nu at the final episode of the MARL, which is equal to 2.07, corresponding to the dashed horizontal black line.

Figure 7 (see also the corresponding multimedia data from  Appendix C) clearly demonstrates that the MARL approach learns much faster and effectively than the SARL approach. This figure shows that, while the SARL approach plateaus after around 150 episodes and displays only very limited Nu reduction, the MARL method learns initially much faster than SARL and keeps learning for a much longer time until it reaches a far better Nu reduction. To be more specific, at the end of 350 episodes, the MARL approach is achievable a reduction in Nu of 22.7%, whereas the SARL approach achieves approximately 5.1% reduction. Moreover, a closer look at the last state of the last episodes in the SARL and MARL cases, respectively, shows that while the MARL agent has managed to change the topology of the flow and to reach a coalesced RBC convection cell state (more details below), the SARL approach has only a limited effect on the flow topology (visible from the videos in  Appendix C). This implies that, while the MARL agent has learned an effective policy and control law, the agent in the SARL framework has failed to do so. This is very similar to the findings reported in Ref. 24, as previously discussed. The reader curious of more details about the flow dynamics in each case is referred to the videos in  Appendix C to see the difference in the controlled-flow configurations. Furthermore, insight into the MARL control mechanism is provided in the next paragraphs.

The price to pay for the DRL ability to discover and optimize non-linear control laws is the lack of direct interpretability of the DRL agent policy. Given a trained DRL agent, and in particular, its policy network weights, it is difficult to derive a detailed, quantitative understanding of the control law using a bottom-up approach, and the DRL policy per se is effectively a black box. This contrasts with, e.g., linear flow-control methods, where properties such as the eigenmodes and eigenvectors of the associated operators can be used to extract some form of explainability.^87,88 However, a higher-level, phenomenological understanding of the DRL agent policy is still obtainable, by applying expert knowledge to the analysis of the control laws exhibited by the DRL agent. Such a phenomenological understanding is what we discuss in the following paragraphs.

More specifically, different phases are clearly visible when the DRL agent controls RBC, as highlighted in Fig. 8 (Multimedia view) (see also the corresponding multimedia data from  Appendix C). Starting from the baseline flow, the DRL agent starts by weakening the RBC cells through modulation of the bottom wall temperature. We will refer to this destabilization of the baseline regime as the Phase I. The destabilization is performed by the agent through applying (i) increased heating on the regions of the bottom plate where a falling plume of low temperature fluid is observed, and (ii) decreased heating (or relative cooling) on the regions of the bottom plate where a rising plume of high temperature fluid is observed. This is visible in Figs. 8(a) and 8(b) (see also the corresponding multimedia data from  Appendix C). This, in practice, weakens the driving mechanism of the baseline RBC cells, by producing excess buoyancy in the downwelling part of the RBC plume. As a result, the control on Phase I reduces the excess of buoyancy in the rising part of the RBC plume, hence, opposing the buoyancy driven flow motion.

This opposition to and weakening of the baseline flow results in an unstable situation, where the double RBC cells are broken. We recognize this transition as the phase II and it can be observed in Figs. 8(c)–8(f) (see also the corresponding multimedia data from  Appendix C) and in Fig. 9. This Phase II ends with a recombination of the two initial convective cells into a single coalesced RBC cell, thus leading to the Phase III. Phase III is distinguished by a single coalesced bubble regime, which is visible in Figs. 8(g) and 8(h) (see also the corresponding multimedia data from  Appendix C).

Therefore, it appears that the DRL agent has learnt to “navigate” in the phase space of the problem. While the phase space configuration corresponding to the baseline configuration is stable, it leads to a suboptimal reward. In particular, the DRL agent has found a strategy (applied during the transition phases I and II) that can bring the system away from the baseline configuration and into the neighborhood of a new topological flow configuration, i.e., the coalesced bubble observed in phase III, which provides higher reward.

Interestingly, the single coalesced bubble obtained in phase III is naturally stable once it is established and is in fact providing the optimal reward, even without control. This is what we define as phase IV in Fig. 9, i.e., the flow obtained upon the re-imposition of the lower-wall constant temperature boundary condition (as in the baseline) after the end of episode. A snapshot of this final state is observed in Fig. 10 (Multimedia view) (see also the corresponding multimedia data from  Appendix C).

Figure 11 shows the Nusselt number obtained through three different episodes at different stages of the learning. While at episode 45, the DRL agent was not able to coalesce the two convective cells into one, at episodes 275 and 350, it has already learned the strategy to perform the coalescence. This is indicated by the large jumps in the Nusselt number at time units 115 and 215 for episodes 350 and 275, respectively. We can also observe that the agent trained for more episodes (episode 350) is able to reach the Phase III (single-cell configuration) earlier than the agent trained for fewer episodes (episode 275). However, once the flow configuration is a single cell, no matter when it is reached, the value of the Nusselt number converges to a minimum of approximately 2.1. This is related to the fact that the single-cell configuration is a stable state with a higher reward, i.e., lower Nu, than the two-cell configuration. If once on phase III the control is turned off, i.e., in phase IV, the Nusselt number achieved gradually converges to 2.05, which is the best Nu we ever observe. Therefore, any small control applied during Phase III is noise (either exploration noise during training, or agent “uncertainty” during evaluation). The black dash-dotted horizontal line in Fig. 11 shows the value of the Nusselt number for phase IV.

In the present work, we demonstrate effective Rayleigh–Bénard Convection (RBC) control by using deep reinforcement learning (DRL). Compared with previous works, e.g., Ref. 23, we consider a domain with a higher aspect ratio, i.e., a wider domain. This leads to the existence of several rolls or RBC cells, and the need for using more separate actuators at the bottom walls. We show that a naive single-agent DRL method fails to learn an effective control policy. This is due to the difficulty presented by having to control several actuators simultaneously, which results in the curse of dimensionality in the control space. By contrast, we demonstrate that leveraging invariant multi-agent reinforcement learning (MARL), which takes advantage of the invariant structure of the underlying RBC control problem, allows to discover an effective control strategy in the case studied. We demonstrate that, in practice, “multi-agent reinforcement learning is all you need” to control our multiple-input, multiple-output fluid-dynamics problem. This finding is similar to previous results,²⁴ and we expect to keep observing it repeatedly when applying DRL to control physical systems with strong underlying structures and symmetries. Using a MARL approach results in a much faster learning than what was obtained in Ref. 23, and is key to our ability to control wider channels with more actuators.

Going into more details, our MARL controller finds a complex, non-trivial control strategy in the RBC configuration under consideration. Through trial and error, the MARL approach discovers that adequate control of the lower wall temperature profile can destabilize the double RBC cell pattern observed in our baseline flow, and force these RBC cells to coalesce into a single convection cell. This convection cell is, in turn, intrinsically stable and yields better heat exchange performance, i.e., a lower value of the Nusselt number Nu, following the definition of our optimization problem.

The present study is a first application of MARL to RBC control, and we provide a strong framework for further studies through the release of all the codes and scripts as open material. Many interesting questions remain open to discussion in future works, and we expect that more studies of the RBC controllability using MARL will take place. In particular, this framework is well adapted to studying problems such as the effect of the Rayleigh number Ra, the channel aspect ratio, the strength of the control authority, etc., on controllability and optimal states of the RBC system. We believe that this work is also a first step toward demonstrating RBC MARL control in both 3D configurations, and in geometries relevant for industrial applications.

Part of the deep-learning-model training was enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS). RV acknowledges financial support from ERC Grant No. “2021-CoG-101043998, DEEPCONTROL.” Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

The authors have no conflicts to disclose.

Colin Vignon: Data curation (equal); Investigation (lead); Validation (equal); Visualization (equal); Writing – review & editing (equal). Jean Rabault: Conceptualization (equal); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Joel Varun Vasanth: Data curation (equal); Writing – original draft (equal); Writing – review & editing (equal). Francisco Alcántara Ávila: Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Mikael Mortensen: Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Ricardo Vinuesa: Conceptualization (equal); Funding acquisition (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

All the codes, scripts, and post-processing tools used in this work are made available on Github, together with readmes and user instructions, at the following address: https://github.com/KTH-FlowAI/DeepReinforcementLearning_RayleighBenard2D_Control.

The shenfun CFD case setup is described in great details, including all the numerical implementation considerations, elements chosen, and code walk-through, on the shenfun RBC documentation page: https://shenfun.readthedocs.io/en/latest/rayleighbenard.html. Note that the conventions used in the present paper and in this documentation page are slightly different: while the present paper uses x as the horizontal (“wall parallel”) description and y as the vertical (“wall normal”) direction, the documentation page uses the opposite conventions. This has no influence on the CFD results per se. The code uses the same conventions as the shenfun documentation, and only the present paper uses separate conventions, for simplicity of the writing and conformity with the literature.

Extra materials are provided as a series of videos on YouTube, which illustrate

• The development and settling of the baseline flow (Fig. 5): see https://www.youtube.com/watch?v=aSKZ1qNgMWk

• The RBC undergoing control by the SARL controller, starting from the baseline flow, which illustrates that SARL fails to learn a topology-changing control strategy (Fig. 7, SARL episode 350): see https://www.youtube.com/watch?v=N9wTi4K2uJY

• The RBC undergoing successful control by the MARL controller, starting from the baseline flow, which illustrates the bubble coalescence phase and the large bubble phase (Fig. 8): see https://www.youtube.com/watch?v=_HYmWee3P0A

• The modified flow configuration corresponding to a single stable RBC cell, obtained by re-applying a uniform bottom wall temperature at the end of the MARL control sequence (Fig. 10): see https://www.youtube.com/watch?v=MEGl2IiybSg