Artificial intelligence control of a turbulent jet

An artificial intelligence (AI) control system is developed to maximize the mixing rate of a turbulent jet. This system comprises six independently operated unsteady minijet actuators, two hot-wire sensors placed in the jet, and genetic programming for the unsupervised learning of a near-optimal control law. The ansatz of this law includes multi-frequency open-loop forcing, sensor-feedback and nonlinear combinations thereof. Mixing performance is quantified by the decay rate of the centreline mean velocity of jet. Intriguingly, the learning process of AI control discovers the classical forcings, i.e. axisymmetric, helical and flapping achievable from conventional control techniques, one by one in the order of increased performance, and finally converges to a hitherto unexplored forcing. Careful examination of the control landscape unveils typical control laws, generated in the learning process, and their evolutions. The best AI forcing produces a complex turbulent flow structure that is characterized by periodically generated mushroom structures, helical motion and oscillating jet column, all enhancing the mixing rate and vastly outperforming others. Being never reported before, this flow structure is examined in various aspects, including the velocity spectra, mean and fluctuating velocity fields and their downstream evolution, and flow visualization images in three orthogonal planes, all compared with other classical flow structures. Along with the knowledge of the minijet-produced flow and its effect on the initial condition of the main jet, these aspects cast valuable insight into the physics behind the highly effective mixing of this newly found flow structure. The results point to the great potential of AI in conquering the vast opportunity space of control laws for many actuators and sensors and in optimizing turbulence.


Introduction
The turbulent jet is one of classical shear flows discussed in virtually every textbook. Its control finds important industrial applications, including dilution jets in combustors, fuel injection of combustion engines, noise mitigation of sub-and supersonic jets for civil and military aircrafts, thrust augmenting ejectors, thrust vector control, etc. The key to control entrainment and mixing processes in a turbulent jet is, as in other shear layers, to manipulate the coherent motions. When a jet issues from a round nozzle, a free shear layer is formed from the nozzle lip and develops downstream. Kelvin-Helmholtz instability (Ho & Huerre 1984) inherent in the shear layer rapidly grows, resulting in the formation of axisymmetric ring vortices. The vortices, along with their subsequent interaction (e.g. merging and breakdown), dominate the shear-layer growth and entrainment (Crow & Champagne 1971). Shortly downstream of the nozzle exit, three-dimensionality becomes an important feature of the flow structure; streamwise vorticity contributes predominantly the entrainment of fluid from the surroundings (Liepmann & Gharib 1992). These motions, formed near the nozzle exit, are featured by a wide range of scales, varying convection velocity and a rich set of three-dimensional patterns (Garnaud et al. 2013); they are sensitive to initial conditions (e.g. the turbulence level, boundary layer thickness, nozzle geometry) and external periodic disturbances (Vlasov & Ginevskii 1973), thus highly susceptible for control.
Many active control studies of turbulent jets involve the open-loop periodic forcing of a pre-specified form, e.g. axisymmetric, flapping or helical forcing. Broze & Hussain (1994) deployed four speakers upstream of the nozzle to add a longitudinal component of perturbation to the mean flow. The acoustic source produced axisymmetric forcing which was found to amplify vortex ring structures and subsequent vortex pairing. Koch et al. (1989) generated helical forcing on a turbulent round air jet using four speakers, each being 90 • out of phase with the adjacent speaker.  used two unsteady radial minijets separated by 60 • or 120 • to produce flapping jet column, which enhanced greatly jet mixing. The combination of individual forcings is also investigated. Juvet (1987) optimized experimentally the combinations of axisymmetric axial and helical forcing to augment mixing. The axial excitation was produced by a loudspeaker placed upstream of the jet and perpendicularly to the centreline. The helical excitation was generated by four external loudspeakers. Their acoustic excitations were directed via waveguides at an angle around the jet exit lip to the shear layer where the flow is most sensitive to acoustic forcing. While the axial excitation led to the early formation of large-scale vortices, the helical caused the shear layer to roll up into staggered vortex structures. This combined excitation generated a bifurcating jet with a much larger spreading angle than the single excitation when the frequency ratio of the axial to that of the helical excitation was equal to 2 (Reynolds et al. 2003). Hilgers & Boersma's (2001) three-dimensional direct numerical simulation of a turbulent jet demonstrated that the superposition of two counter-rotating helical modes of the same excitation frequency f e and one axial excitation of 2f e produced a bifurcating jet whose centreline mean velocity and scalar concentration decayed faster than those of the counter-rotating helical actuation alone. Tyliszczak & Geurts (2015) and Tyliszczak (2018) simulated highly mixed multi-armed bifurcating jets using axial and helical excitations. However, it would be very difficult or time-consuming for conventional active controls to find the globally optimal solution for the combined excitations where many control parameters are generally involved. For instance, the control optimization of a turbulent jet has so far typically involved up to two control parameters, such as the actuation amplitude and frequency. Then, the optimization of combined modes, like axisymmetric forcing and flap-ping forcing, may involve at least four independent control parameters, i.e. the amplitude and frequency of each mode (e.g. Hilgers & Boersma (2001)). The search for its optimal solution is then already a challenge. If the control parameters for each mode is increased to three or four such as the amplitude, frequency, duty cycle and diameter ratio of an unsteady jet (e.g. Perumal & Zhou (2018)) or multiple independent actuators are deployed, the search for the globally optimal solution of the combined modes would be a daunting task. Koumoutsakos et al. (2001) and Hilgers & Boersma (2001) have pioneered the jet mixing optimization with 3 and 4 actuation parameters using Rechenbergs 1973 evolutionary strategy.
Model-based control comes, if doable, with the deep understanding of actuation dynamics, regardless of open or closed loops. In simulations, the linear dynamics can be accurately resolved by discretized Navier-Stokes (N-S) equations (Kim & Bewley 2007;Sipp et al. 2010). In experiments, linear stochastic estimation (Tinney et al. 2006) has been successfully applied to resolve the flow physics from measurement signals and PIV measurements. The linearized N-S dynamics can be encapsulated in reducedorder models employing several dominant non-normal global stability eigenmodes. The downstream evolution of wave-packets can be real-time estimated in a high-Reynoldsnumber turbulent jet thanks to the development of transfer functions based on the parabolized stability equations (Sasaki et al. 2017). So can the closed-loop control of fluctuations in a low-Reynolds-number shear layer (Sasaki et al. 2018). These controloriented models have significantly contributed to the understanding of the manipulated jet dynamics.
Model-free approaches may yield performance benefits from nonlinear dynamics which is too complex for control-oriented models. A new model-free self-learning approach for general non-linear control laws has been developed by Dracopoulos (1997) for commanding satellite motion and was re-discovered in fluid mechanics as Machine Learning Control or MLC (Gautier et al. 2015). A review of dozens of MLC experiments and simulations is provided by Noack (2019). The first MLC experiment was set to enhancing shear layer mixing with 96 jet actuators driven in unison and 25 hot-wire sensors for feedback control (Parezanović et al. 2016). The optimization of shear-layer mixing resulted in destabilizing phasor control, i.e. the feedback excitation of the dominant frequency. The control enhanced and synchronized downstream large-scale vortices with a frequency of one sixth of the initial Kelvin-Helmholtz (K-H) instability. Li et al. (2017) deployed four Coanda jet actuators placed at the trailing edge of an Ahmed body and 12 pressure sensors at the back side, achieving a drag reduction by 22% (where about 10% was attributed to the Coanda effect) when the excitation frequency was much higher than the predominant and even the shear-layer frequencies in the wake. Using an unsteady single jet actuator driven by an online PIV-based sensing, Gautier et al. (2015) cut short the reattachment length of flow over a backward-facing step. They observed surprisingly the enhancement of a low-frequency flapping mode, instead of the excitement of the dominant K-H vortex shedding. MLC has matched with or outperformed existing control strategies and solved the combined task of picking the nonlinear mechanisms for performance optimization and selecting the best sensors. These model-free control studies show that the actuation mechanism can be very difficult to anticipate, thus implying a challenge to any model-based control.
It is worth pointing out that MLC has never been applied to multiple independently operated actuators, resembling a distributed actuation, in experiments so far. MLC laws, previously developed, have been of small to moderate complexity, e.g. the phasor control, threshold-level based control, periodic or two-frequency forcing , as the actuators are typically driven by a single actuation command. Indeed, the use of independent actuators may increase dramatically the level of control complexity. For example, assume that one unsteady minijet, used to maximize jet mixing, involves three parameters, i.e. the actuation frequency f a , velocity U a and duty cycle α. Then, if the number is increased to up to say six independent minijets spatially distributed around main jet, the independent control parameters will be tremendously increased. Then one naturally wonders what the globally optimal solution of the problem is and whether an AI control technique could be developed to find this solution. Furthermore, what turbulent flow structure might this global optimal solution or forcing produce?
This experimental work aims to address the issues raised above and to optimize jet entrainment/mixing with six independently unsteady minijets placed upstream of nozzle exit, extending the MLC jet control using a single unsteady minijet (Wu et al. 2018a). The manuscript is organized as follows. The experimental setup and minijet-produced flow, along with its effect on the jet initial conditions, are described in § 2 and § 3, respectively. The following three sections 4, 5 and 6 describe the AI control system developed, the outcome of the AI-based learning and the resulting turbulent flow structures. The work is concluded in § 7.

Jet facility
Experiments were conducted in a round air jet facility, as schematically shown in figure 1. The facility was placed in an air-conditioned laboratory where the room temperature remains constant within ±0.5 • C, centrally deployed in an area of approximately 2.5 m in width and 2 m in height, enclosed by fabric walls. In order to minimize the effects of the wall on the jet, the nozzle exit is 4.0 m away from the fabric partition wall and the distance is well over 70 times jet exit diameter required for neglecting the wall effects (Malmstrom et al. 1997). As the jet is highly sensitive to background noise, careful measures are taken to avoid any external interference to airflow.
The compressed air of the round jet comes from a constant 5 bar gauge pressure, mixed with seeding particles in the mixing chamber in the case of the particle image velocimetry (PIV) or flow visualization measurements, and then enters into a plenum chamber, composed of a 300 mm long diffuser of 15 • in half-angle and a 400 mm long cylindrical settling chamber with an inner diameter of 114 mm. The flow passes two screens before entering the smooth contraction nozzle (Perumal & Zhou 2018), which is extended by a 47 mm long smooth tube of the same diameter as the nozzle exit D (= 20 mm). The Reynolds number Re D = U j D/ν of the main jet is fixed at 8000, where U j is the centreline velocity measured at the nozzle exit, the overbar denotes time-averaging and ν is the kinematic viscosity of air. A Cartesian coordinate system (x, y, z) is defined in figure 1(a, c), with its origin at the centre of the jet exit and the x-axis pointing in the direction of flow. Measurements were conducted in the x − z, x − y and y − z planes of the main jet. The instantaneous and fluctuating velocities in the x, y and z directions are denoted by (U, V, W ) and (u, v, w), respectively. Six unsteady control minijets issued from orifices with a diameter of 1 mm are equidistantly placed around the extension tube at x i = −0.85D, y i = (D/2) cos θ i , zi = (D/2)θ i , where θ i = (i − 1)2π/6, i = 1, 2 . . . , 6 (figure 1b, c). Their mass flow rate is determined by a flow-limiting valve and monitored by a mass flow meter, with a measurement uncertainty of 1 %, and the frequencies and duty cycles are independently controlled by individual electromagnetic valves that are operated in an ON/OFF mode. The maximum operating frequency of the valves is 500 Hz, exceeding three times the preferred-mode frequency, f 0 = 135 Hz, of the unforced jet, the corresponding Strouhal number being S t = f 0 D/U j = 0.45, where f 0 is obtained from the power spectral density function E u of streamwise fluctuating velocity u measured in the absence of control .

Flow measurements
The fluctuating flow velocities are monitored by two tungsten wire sensors of 5 µm in diameter, operated on a constant temperature circuit (Dantec Streamline) at an overheat ratio of 0.6, one placed at (x/D, y/D, z/D) = (0, 1/4, 0) and the other at (x/D, y/D, z/D) = (5, 0, 0). The time-averaged velocity at the latter position is denoted by U 5D . This choice is based on following considerations. Firstly, Zhou et al. (2012) demonstrated that thus defined K is correlated approximately linearly with an equivalent jet half width R eq = [R H R V ] 0.5 , where R H and R V are the jet half-widths in two orthogonal planes, implying that K is directly connected to the entrainment rate of the manipulated jet. Secondly, Fan et al. (2017) found that the difference K between the K values with and without control reaches the maximum at x/D ≈ 5, that is, the centreline mean streamwise velocity measured at x/D = 5 is most sensitive to the change in the control parameters.
Finally, the variation in K is almost linear from x/D = 0 to x/D ≈ 7 under control (Fig 7, Fan et al. (2017)), that is, a single value of K may be used to describe reasonably well the jet decay rate in the near field under control. Both hot-wires are calibrated at the jet exit using a pitot-static tube connected to a micromanometer (Furness Controls FCO510). The cutoff and sampling frequencies are 3 kHz and 6 kHz for open-loop control experiments, respectively. The experimental uncertainty of the hot-wire measurement is estimated to be less than 2%.
A planar high-speed PIV system, with a high-speed camera (Dantec Speed Sence 90C10, 2056 × 2056 pixels resolution) and a pulsed laser source (Litron LDY304-PIV, Nd: YLF, 30 mJ/pulse), is deployed for velocity field measurements in the x−z, x−y and y − z planes. An oil droplet generator (TSI MCM-30) is used to generate fog from olive oil with an averaged particle size of 1 µm for flow seeding. Flow illumination is provided by a laser sheet of 1 mm in thickness generated by the pulsed laser via a cylindrical lens. For velocity measurements in the x − z and x − y planes, the captured image covers the area of x/D ∈ [0, 6] and y/D, z/D ∈ [−2, 2]. The longitudinal and lateral image magnifications are identical, 0.09 mm per pixel. The time interval between two consecutive images is presently chosen to be 25 µs, which is found to yield satisfactory results. There are 253×253 velocity vectors, the same for the two planes. A total of 200 pairs of flow images are captured at a sampling rate of 405 Hz for each set of PIV data. In post-processing, a built-in adaptive correlation function of the Flow Map Processor (PIV 2001 type) is applied with an interrogation window of 32 × 32 pixels and a 75% overlap along both directions. The same PIV system is used for flow visualization in the three orthogonal planes. So are the seeding particles, though their concentration is higher than in the PIV measurements to provide a clear picture for the flow structure. The captured images cover the area of x/D ∈ [0, 6] and y/D or z/D ∈ [−2, 2] in the x − y and x − z planes and the area of y/D = z/D ∈ [−2, 2] at x/D = 0.25 in the y − z plane.

Real-time system
A National Instrument PXIe-6356 multifunction I/O device, connected to a computer, is used in experiments to generate the real-time control command at a sampling rate of F rf = 1 kHz. A LabVIEW Real-Time module is used to execute the program. Sensor data acquisition and control command generation for the AI control experiments are operated under the same sampling frequency of 1 kHz. It has been confirmed that the ON/OFF command lasts at least 1 ms to ensure the actuators to work effectively. The available f a can be derived from f a = F rf /N sp , where N sp is the number of sampling points in one period 1/f a . The working frequency range of actuators ([0, 500 Hz]) imposes a minimum value for N sp , i.e. N sp 2. For a given frequency, α can be deduced from m/N sp , m = 1, . . . , N sp -1. The m range ensures a response time of 1 ms for the effective working of the actuators, which is adequate as the maximum sampling rate F rf is 1 kHz due to the limitation of hardware. Thus, the number of possible duty cycles N α for a given f a is N α = N sp − 1 = F rf /f a − 1, which increases with F rf and decreases with f a . This process is similar to the one used by Li et al. (2017) and Wu et al. (2018a).

Minijet-produced flow
It is important to document the flow produced by a minijet and the effect of minijets on or the initial condition of main jet. This information is crucial for understanding physically the manipulated jet. The instantaneous velocity U a of a single radial minijet is first examined in the absence of main jet. A hot-wire is placed 17 mm or x/D = -0.85 upstream of the main jet exit and 3 mm radially from the exit of minijet 1 (figure 1c). The hot-wire is oriented normal to the minijet axisrecording the signal U a , which changes with α (figure 2). For α = 0.1, U a displays sharp peaks which are periodic and clearly separated. But these peaks are less pronounced at α = 0.4. U a is almost steady at α = 1, though showing a small variation, as observed by Johari et al. (1999). Apparently, a small α produces a large instantaneous velocity, implying a large penetration depth into main jet.
Consider the simultaneous injection of minijets 1 and 4 (figure 1c) without main jet. Two hot-wires are placed perpendicularly to the x − y plane at x/D = -0.85 and 3 mm from each of the corresponding measured minijet exit. The two minijets are injected with a phase shift Φ, which may be varied by changing the phase shift between the two square wave signals of input voltages. At Φ = 0 • , the U a1 signal exhibits a very sharp peak value, with a magnitude of close to 0 at the off-state of the minijet and about 13 at the on-state (figure 3a). Note that, even after the electromagnetic valve is closed, there may be some fluid injecting into main jet (Sailor et al. 1999). A similar observation can be made for Φ = 60 • and 180 • (figure 3b, c). The characteristics of U a2 resemble those of U a1 , regardless of the Φ value. It may be inferred that each of the minijets does not depend on Φ and is rather independent of each other.

Penetration depth and minijet number
The penetration depth of control jets may have a pronounced impact upon jet mixing (Davis 1982). Thus, its influence on main jet is examined for various minijet numbers and configurations with given U j and control parameters, i.e. α, f a /f 0 and mass flow ratio C m = m mini /m j , where m mini and m j are the mass flow rates of a single minijet and main jet, respectively. The minijet penetration depth could be approximately estimated from the U signals along the y direction, measured at x/D = 0.05 using a hot-wire placed perpendicularly to the x − y plane, as shown in figure 4, where the scale of the abscissa or ordinate is made the same for all cases to facilitate comparison. The U signals are essentially constant throughout the range of y/D ∈ [-0.3, 0.3] for the unforced jet. The periodic fluctuations of U appear at y/D = -0.3 for one minijet injection (N = 1), Figure 4. Typical hot-wire signals of instantaneous streamwise velocity U/Uj along x−yplane at x/D = 0.05 for minijet number N = 1, 3, 4 and 6 (fa/f0 = 0.5, α = 0.15, Cm =1.2%). The same scale is applied for all signals. The dots within the circle represent the hot-wire measurement points. and its magnitude grows first and then retreats with increasing y/D. The fluctuations remain discernible at y/D = -0.1. Note that the minijet is issued along the y direction. Beyond y/D = -0.1, the velocity fluctuations are negligibly small and in fact comparable to that in the unforced jet. These observations indicate that the minijet has reached a penetration depth of y/D = -0.1. With three adjacent minijets on (N = 3), the velocity fluctuations are appreciably larger in magnitude than their counterparts of N = 1, and the maximum amplitude is shifted to a deeper position, i.e. from y/D = -0.1 at N = 1 to y/D = 0 at N = 3. The fluctuations are now discernible at y/D = 0.2, indicating an increased penetration depth, though the minijets clearly have not impinged on the wall opposite to the injecting minijets. With N increasing to 4, the maximum magnitude of the velocity fluctuations is appreciably larger than that of N = 3, and again occurs at the centre (y/D = 0) where all minijets contribute to an increase in the velocity fluctuations. Furthermore, the fluctuations are now even discernible at y/D = 0.3. It is worth pointing out that we did not move the hot-wire closer to the wall because of its high fragility; therefore, we could not tell whether the minijets have penetrated through the main jet in this case. At N = 6, the velocity fluctuations display symmetry about the centre, the maximum magnitude exceeding all other cases and taking place at the centre.
3.3. Power spectral density function and minijet number  of the preferred mode structure. When the minijets as well as main jet are operated, E u exhibits more pronounced peaks at f /f a = 1 and its harmonics. These observation results from the interaction between main jet and minijet, referred to as the parametric resonance by Huang & Hsiao (1999). Evidently, the unsteady injection produces the periodic structures upstream of the nozzle exit, as noted by Zhou et al. (2012). With increasing N , the peaks become more pronounced and occur at more harmonics, echoing the enhanced periodic structures (figure 4) and hence the enhanced excitation of the shear layer. The predominant frequencies do not vary with N though. Note that E u is normalized by u 2 so that its integration over the entire frequency range is always equal to unity. As a result, E u drops appreciably over the low frequency range.

Fluctuating velocity and minijet number
The number and configuration of minijets may profoundly affect the main jet issuing from the nozzle, as the case of passive delta tabs (Zaman et al. 1994). This effect plays an important role in the downstream evolution of flow. As such, the radial profiles of the hot-wire measured root mean square (rms) velocity u rms at x/D = 0.05 are examined in the manipulated jets, along the y and z axes, respectively, for N = 1, 2, 3, 6. The data of the unforced jet are also presented for the purpose of comparison. Given the symmetrically arranged minijets (N = 2, 6) about the x − z plane (figure 6b1,d1), the u rms distributions along the y-axis exhibit reasonable symmetry. The u rms displays a pronounced peak about y/D = ±0.45, where the shear layer is expected, in the injection or x − y plane for N = 2 (figure 6b1) but remains unchanged in the orthogonal or x − z plane (figure 6b2), indicating that the shear layer between the two minijets is essentially undisturbed. Being symmetrical about the z axis, the u rms distributions are given only for z/D 0 in figure 6(a2 − d2). A broad bump is evident at y/D ≈ 0.2 for N = 2 (figure 6b1). The flow structure induced by an unsteady injecting minijet is similar to a pulsed jet in cross flow, which forms a series of periodical vortex rings (M'closkey et al. 2002). It seems that these minijet-produced periodic vortices may occur most likely at y/D ≈ 0.2, accounting for the broad bump. For N = 3 and 6, this bump moves to near the centre, with a significantly increased magnitude (figure 6c1, d1). Two factors may be responsible for this increase. Firstly, as the separation angle θ decreases from 180 • to 120 • and then 60 • , two neighbouring minijets become close and their induced unsteady flows interact more and more intensely. Zaman et al. (1994) noted that, as the neighboring delta tabs approach each other, streamwise vortices interact more vigorously, resulting in the jet core fluid ejection. Secondly, as demonstrated in figure 4, every minijet generates a velocity fluctuation at the centre. For N = 6, θ is smallest and all six minijets contribute to flow perturbations, thus producing the most pronounced bump at y/D = 0.

Artificial intelligence control system
Artificial intelligence methods allow us to explore the rich universe of nonlinear actuation mechanisms opened by independent spatially distributed actuators. Hence, we see the actuation and sensing hardware and control logic as intimately interwoven. The AI control system is sketched in figure 7. Generally, a control system facilitates a control goal for a plant by control hardware and a control logic/controller. The control hardware includes sensors and actuators as discussed in § 2. This hardware monitors the plant output (velocity signals) and executes instructions from the controller. The open-loop arrangement is shown in figure 1a for calculating the cost value J = U 5D /U j .
A minimized cost J corresponds to the maximized decay rate K = 1 − J of jet centreline mean velocity, which is an indicator of the mixing efficacy of a jet (Perumal & Zhou 2018).

Control optimization using linear genetic programming
The six-dimensional vector b = [b 1 , b 2 , . . . , b 6 ] † comprises all actuation commands. The ith minijet is 'ON' if the actuation command b i is positive and is 'OFF' otherwise. In the sequel, we assume b i = 1 for 'ON', and b i = 0 for 'OFF'. Following Wu et al. (2018a), we search for a control law including sensor-feedback with hot-wire signals s, multifrequency open-loop forcing with harmonic functions contained in where the vector function K = [K 1 . . . , K 6 ] comprises the actuation laws for each minijet. The time-averaged duty cycle of the ith minijet is determined by the control law K and the arguments, i.e. sensor signals and harmonic functions. For open-loop forcing b = K(h), the duty cycle of the ith minijet becomes the sensor-independent time-average of the actuation command K i (h). Thus, helical forcing may have a particularly simple representation, e.g. b i = h i . In general, only two harmonic functions, typically sin ω a t and cos ω a t, are sufficient for harmonic functions with arbitrary phases. Following Paschereit et al. (1995) and others, we add the cosine and sine components of ω a /2 and ω a /4, yielding a ten-dimensional vector h = [h 1 , h 2 . . . , h 10 ] † . The nonlinear function K can create arbitrary higher harmonics, arbitrary phase relationships between ω a , ω a /2, ω a /4, and higher harmonics, e.g. 1 − 2h 10 = cos(10ω a t), as well as arbitrary sum and difference frequencies. The control optimization searches for a law of form 4.1 that minimizes the cost, The regression problem implies a search for a mapping from multiple inputs to a multipleoutput signal. Even in case of a linear function this implies the optimization of a large number of parameters. We employ the powerful linear genetic programming (LGP) as a regression solver and take the same parameters for the control law representation and for the genetic operations as Wu et al. (2018a). The first generation of LGP, n = 1, contains N i = 100 random control laws, also called individuals. Each individual 'i' is experimentally tested for 5 seconds to yield the measured cost J n i , where superscript 'n' denotes the generation number. Subsequent generations are produced from the previous ones with genetic operations (elitism, crossover, mutation and replication) and tested analogously. Elitism pass directly the top-ranking individuals to next generation. Replication copies a stochastically selected number of individuals into next generation, which acts to preserve some well performing individuals. Crossover involves two selected individuals and then produces two individuals, with part of their elements exchanged. This operation tends to generate better individuals by exploitation. For the mutation operation, the instructions of a selected individual are randomly changed. Both crossover and mutation serve to explore potentially new and better minima of J. After the in situ performance measurements, the individuals are re-numbered in order of performance, J n 1 J n 2 . . . J n Ni , where subscript i represents the individual index and N i and n denote the size and number of generations. We have noted in trial tests that all the winning individuals always involve every actuator. Therefore, when generating the 100 individuals of the first generation, we exclude the possibility of permanently inactive actuator to accelerate the learning process, that is, as a plant-specific rule, we discard and replace any individual for testing if one or more actuators are not active. It is worth mentioning that the present jet control is formulated as a model-free regression problem: determine the law which minimizes the given cost function. The considered search space of control laws significantly extends hitherto considered actuations. First, general multiple-input actuation is allowed without any imposed symmetry constraints. Thus, actuations with arbitrary combinations of minijets thereof can be realized. Second, the search space includes broadband multi-frequency actuation. Third, nonlinear sensorfeedback is included, which is made by nonlinear operations with the sensor signal s, e.g. b = log 10 (s) (Wu et al. 2018a). However, this feature is not found improving appreciably the control performance and is therefore removed eventually in the learning process. Fourth, the control law may include nonlinear combinations of multi-frequency forcing and sensor-feedback. The key enabler for the control optimization in this search space is genetic programming as powerful regression solver. Genetic programming may be considered as an example for the many powerful regression solvers of AI.

Parameters and control landscape
The LGP parameters for this study are displayed in Table 1. These values are adopted from a previous MLC jet mixing study in the same facility with a single minijet (Wu et al. 2018a). The parameters are identical or close to the ones employed in numerous experimental studies as summarized by Duriez et al. (2016) and Noack (2019). Elitism is set to N e = 1, i.e. the best individual of a generation is copied to the next one. The replication, crossover and mutation probabilities are 10%, 70% and 20%, respectively. The individuals on which these genetic operations are performed come from a tournament selection of size N t = 7. The instruction number varies from 10 to 50. The operations comprise +, −, ×, ÷, sin, cos, tanh, log 10 and g 2 , where g is the input argument. The operations ÷ and log 10 are protected to prevent an undefined expression with a vanishing argument; for example, log 10 (g) is modified to log 10 (|g|). In addition, LGP uses three random constants in the range [-1, 1].
The evolution of control laws is depicted with a proximity map following Duriez et al. (2016). The main idea is that the considered ensemble of K i (h) is represented as points in a two-dimensional feature plane γ i = (γ i,1 , γ i,2 ), where i = 1, 2 . . . , N i × n, so that the difference between the control laws is optimally indicated by distance between feature vectors. The key is the definition of a metric D ij between the control laws K i (h) and K j (h). For the considered open-loop actuation, this metric is the root-mean-square averaged Euclidean difference between the actuation command vectors accounting for a potential time-delay, given by In the employed metric, we incorporate also the control performance J i by a penalization term, i.e.
The parameter β is chosen so that the maximum actuation distance M ij is equal to the maximum difference in the performance terms: Given the resulting configuration matrix D = (D ij ) (i, j = 1, 2 . . . , N i × n), classical multi-dimensional scaling (Cox & Cox 2000) uniquely determines feature vectors γ i , i = 1, 2 . . . , N i × n, so that the distances are optimally preserved: The translational degree of freedom is removed by centering the feature vectors Ni×n i=1 γ i = 0. The feature vectors are sorted and rotated so that the first coordinate has the largest variance, the second coordinate the second largest, etc. The coordinates are indeterminate by a sign (mirroring), like POD modes and their amplitudes.
Finally, a control landscape J(γ) is interpolated from the three-dimensional data points (γ i,1 , γ i,2 , J i ), i = 1, 2 . . . , N i ×n. The two-dimensional feature vectors γ i are connected by an unstructured grid from a Delaunay (1934) triangulation. This triangulation guarantees that the mesh triangles are optimally equilateral. The J-values in each mesh triangle i 1 , i 2 , i 3 ∈ {1, 2 . . . , N i × n} are interpolated from the known values at the vertices J i1 , J i2 , J i3 . These control landscapes have been employed in several AI-based control schemes . They indicate the complexity of the actuation response and the learning progress of AI-based control. Often, the feature coordinates can be linked with the physical properties of actuation a posteriori, thus providing additional insights.

Representative reference actuations
A few well-known reference forcings are firstly presented to facilitate the understanding of the AI learning process and highlight the uniqueness of this method. In our earlier studies, turbulent jet mixing has been optimized for the same cost function and experimental conditions. For single unsteady minijet forcing, the optimal f a is found to be 67 Hz (Wu et al. 2018a), 0.5f 0 , and the optimal C m is 1.2% based on a dual-input-and-one-output closed-loop control technique (Wu et al. 2018b). As such, we choose the same f a or ω a = 2πf a and C m = 1.2% for every minijet. With C m fixed for each minijet, the overall mass flow of injected fluid in one actuation period T a is the same for all actuations, that is, the input/actuation energy is the same, irrespective of control modes or laws. Consider three reference forcings (e.g. Hilgers & Boersma (2001), , Yang (2017)), viz.
The constants α a , α h and α f correspond to the duty cycles and have been optimized with respect to the cost. As mentioned before, actuation is performed only when b i > 0. The cost functions are found to be J a = 0.665, J h = 0.568 and J f = 0.423 for the optimized axisymmetric, helical and flapping forcings (Table 2), respectively, based on the conventional open-loop control, which provide the benchmarks for the AI control performance to be discussed below.

Learning process of AI control
In the initial stage of the learning process, we included a feedback signal s(t) = [u 3D , u 5D ] as one input (eq. 4.1 for the AI system, where u 3D and u 5D are the fluctuating velocity signals measured at x/D = 3 and 5, respectively. It is found from dozens of experiments that the search for the optimal solution benefits neither from the subharmonic components of h 7 , . . . , h 10 nor from the feedback signals. Therefore, the AI control laws 4.1 may be cast in the periodic open-loop form, b i = K i (h 1 , h 2 . . . , h 6 ), i = 1, 2, . . . , 6. Hence, we restrict our following discussion to this open-loop actuation.
The learning curve of the AI control is presented in figure 8, where the square symbol marks the first and best individual of each generation with N i = 100 control laws. The remaining costs grow monotonously with their indices, and the 100 individuals of each generation form a color bar. The square symbol curve unveils the best performance from generation n = 1 to 30. The best individual of the first generation or stage 1 is characterized by an axisymmetric control law (see Eq. A 1 in Appendix A).
This law is equivalent to 5.1a except for a time shift, reflected by 4/6π. The performance J 1 1 = 0.626 (see figure 8 and table 2) is slightly better, about 5.8% lower, than the benchmark of axisymmetric forcing 5.1a, though much higher than that (J u = 0.947 or K ≈ 0.05) of the unforced jet. Note that the centreline mean velocity for calculating K or J a is measured over a duration of 60 seconds in the benchmark experiments, but only 5 seconds for estimating J i as the measured J i is used to evaluate control laws and does not need to be very accurate in the learning process of AI control. An accurate longtime evaluation of J is performed only in the last generation n = 30. This difference in evaluating the cost function could account for the deviation between J 1 1 and J a . However, the AI control or specifically genetic programming breeds several copies of the winning individual, covering all possible combinations of the control parameters, and then takes the best performing one. This process differs in essence from the searching process of the conventional open loop control which optimizes one control parameter first and then moves to next with the first parameter fixed. The advantage of the former over the latter is obvious. Therefore, different search strategies cannot be excluded from the mechanisms behind the deviation, which will be further substantiated by the fact that all the best cost functions of the different stages in the learning curve are less, albeit slightly, than their corresponding benchmarks produced from the conventional open-loop control.
Stage 2 starts with the second generation when the AI control discovers a better performing helical forcing (Eq. A 2 in Appendix A). This forcing differs in form from 5.1b, but clearly shows a uniformly traveling wave in the azimuthal direction (to be demonstrated in § 5.3), its cost J 2 1 being again slightly lower than J h (Table 2). Helical forcing reduces J further as found from the numerical simulation study of a similar jet mixing optimization (Hilgers & Boersma 2001). Local spatial stability analysis indicates that, unlike axisymmetric forcing, helical perturbations are spatially amplified downstream of the potential core (Garnaud et al. 2013).
Flapping forcing takes place in stage 3, starting from the fifth generation. The law (Eq. A 3 in Appendix A) is similar to 5.1c but incorporates an asymmetry. An optimized asymmetry yields a reproducibly better mixing, again J 5 1 J f (table 2). The eleventh generation marks the emergence of stage 4. AI control discovers a very sophisticated control law. See Eq. A 4 in Appendix A. This forcing significantly outperforms the flapping forcing found in generation 5, the corresponding J 11 1 plunging to 0.305, a drop of 27% compared with the smallest J 5 1 in generation 5 and less than 1/3 of the unforced jet. The actuation mechanism does not change any more in following generations with little variation in costs, pointing to the convergence of the AI learning process. It is worth highlighting that this actuation mechanism is reproducible, that is, approximately the same converged cost has been observed in all experiments, notwithstanding a change in the initial parameters of the first generation. However, not all AI learning curves go through the stages of axisymmetric, helical and flapping forcings; some AI experiments may find only two of the three stages in the learning process.

Representative control laws and flow responses
The control mechanisms may be elucidated from the analysis of the spatio-temporal actuations, shown in figure 9, extracted from the control laws of n = 1, 2, 5 and 11. Each circular pie corresponds to one sixth of the excitation period, while its six sectors represent the six minijets. The arrow indicates that the minijet is on and the radial depth of the blue area is proportional to the duration when the minijet is injecting. The spatio-temporal actuation is found to be internally consistent with the cross-sectional and streamwise flow structure shown in figures 10 and 11. Axisymmetric forcing (figure 9a1 − a6) is characterized by simultaneous blowing of all minijets and a small α of 13.3%. As a result, the cross-sectional flow structure (figure 10a1 − a6) is axisymmetric, and the ring vortex is evident.  discussed in detail the distortion, formation of longitudinal structures and generation of mushroom-like structures in the braid region between ring vortices, and presented a scenario on how the interactions between the longitudinal structures, mushroom-like structures and ring vortices enhance entrainment and mixing. The six synchronized minijet excitations greatly strengthen the ring vortices, as shown in figure 11a (cf. figure 11u). For helical forcing (figure 9h1 − h6), two or three minijets are blowing simultaneously at any instant, with α reaching 40%. These blowing actions rotate clockwise from (h1) to (h6). The greatly increased α, probably required for the generation of helical motion, may act to inhibit the occurrence of mushroomlike structures (Perumal & Zhou 2018), which are absent in figure 10(h1 − h6). The jet (figure 11h) exhibits more spread than axisymmetric forcing (figure 11a). In case of flapping forcing (figure 9f 1 − f 6), three adjacent minijets are blowing simultaneously at one instant with α = 13.3% and are switched to the other three after a phase shift of π, thus creating the asymmetric flapping jet column (figures 10f 1 − f 6, 11f 1 − f 2). The actuation configuration of the rightmost column (n = 11) in figure 9 is complex. Firstly, the number of the simultaneously injecting minijets can be 1, 2, 3 or 4, mostly adjacent to each other. Secondly, the injecting minijets tend to rotate clockwise. The resulting effect is to produce both helical and flapping motions (figures 10c1 − c6, 11c). For convenience, hereinafter we refer to this forcing as the combined mode. Thirdly, α varies from one blowing minijet to another, from about 6.7% to 53%. The effect could be twofold. On one hand, the varying duty cycles of the blowing minijets may yield a resultant blowing force not going through the jet centreline ( figure 9c1−c6). This may produce a precession effect, causing additional jet column oscillation (Wong et al. 2004). On the other hand, whilst a small α facilitates the generation of mushroom-like structures (figure 10c1−c6), a large value enhances the strength of the flapping motion (figure 11c), as noted by Perumal & Zhou (2018). All the features, confirmed by more detailed analysis in section 6, act to promote mixing, consistent with the observed minimum J (figure 8). The generation of such a sophisticated control mode, along with the generation of a complex turbulent flow structure, would have been extremely challenging for conventional control techniques, be it open-or closed-loop and model-based or model-free approaches. . Actuation associated with the best individuals of generation n = 1, 2, 5 and 11 of figure 8. From top to bottom: six instances associated with φ = Ta/6, 2Ta/6, 3Ta/6 . . . , Ta in one excitation period. Each of the six sectors of the circles corresponds to one minijet.

Control landscape: cartographing all actuation laws
Proximity maps ( § 4.3) provide a very revealing illustration of the evolution process of the control laws The underlying metric between two control laws b and b is given by D ij (Eq. 4.4). Figure 12 presents the proximity map of the control laws in a two-dimensional plane such that this metric is optimally preserved. This plane is spanned by the feature Figure 10. Sequential photographs of the cross-sectional flow structure captured at x = 0.25. From top to bottom: six instances at ti = iTa/6 (i = 1, 2 . . . , 6) in one actuation period Ta (= 1/fa).
coordinates γ 1 and γ 2 , which are derived from a mathematical optimization process. The details of deriving γ 1 and γ 2 from individual generations are given in Cox & Cox (2000). Physically, the distance between two points, which are given in terms of (γ 1 , γ 2 ), in the plane is directly linked to the extent how closely similar to each other two control laws are.
The subfigures display the feature coordinates of the four discussed generations (n = 1, Figure 11. Flow visualization of (u) unforced jet (a) axisymmetric, (h) helical, (f ) flapping, and (c) combined forcings from the figure 8, respectively. White ellipses and crooked arrows indicate vortex rings and their rotation, respectively, and yellow arrows highlight the jet column oscillation.
2, 5 and 11). Interesting observations and inferences could be made from the subfigures. For n = 1 (figure 12a), the points appear rather randomly distributed. The optimal solution of axisymmetric forcing, indicated by an open circle, occurs at (γ 1 , γ 2 ) = (-0.04, 0.54), where γ 2 is the largest of all points. At n = 2 (figure 12b), these points appear forming three separated regions numbered 1 through 3, as enclosed by the elliptic contours, where most of the control laws fall in. The optimal solution that appears in Region 2 and corresponds to helical forcing now occurs at (γ 1 , γ 2 ) = (0.39, -0.01), where γ 1 is the largest of all points. At n = 5, the individual points tend to populate along discrete curves (figure 12c), which is a commonly observed phenomenon of the AI control, e.g. Li et al. (2017) and Wu et al. (2018a). The optimal solution now takes place at (γ 1 , γ 2 ) = (0.49, -0.01), where γ 1 is again the largest of all, within Region 3 and corresponds to flapping forcing. It has been confirmed that most of the data points within Regions 1, 2 and 3 correspond to axisymmetric, helical, and flapping forcings, respectively. By n = 11 (figure 12d), the optimal solution occurs at (γ 1 , γ 2 ) = (0.31, 0.17) within Region 2, corresponding to the combined mode. It is noteworthy that those points in this region now correspond to either helical or combined forcing; this very fact may suggest that the so-called combined forcing may have been developed from and more closely related to helical forcing. Interestingly, Region 1 where axisymmetric forcing takes place contracts to a single point, while Region 2 or 3 is reduced to a curve. The cost J associated with each point is color-coded from red (J = 0) to pink (J = 1) as indicated by the color bar. The cost values within each colored area are interpolated from the 100 individuals in each generation. The control landscape at n = 1 displays quite a few minima, as indicated by green colored areas, in the feature plane (figure 12a), suggesting the complexity of the learning task in the early stage, internally consistent with the rather random distribution of the points. The landscape becomes simpler towards the rightmost boundary of the generation. With increasing number of generations, the individuals tend to line up on the ridge-curves marking the cost valleys.
The feature coordinates have technically no a priori meaning. However, an approximate physical meaning of γ 1 and γ 2 may be inferred from the careful analysis of the control laws. The coordinate γ 1 is correlated with the degree of asymmetric forcing, while γ 2 is linked to the number of simultaneously injecting minijets. Consider six arbitrarily chosen control laws A − F in the first generation (figure 12a). Figure 13 shows the on-and offstates of the minijets corresponding to the six control laws. Evidently, the maximum number of simultaneously injecting minijets increases from one (law A) to six (law F ). Axisymmetric forcing at n = 1 (figure 11a) leads to an early disruption of the potential core and a significantly improved entrainment, as suggested by the substantially reduced 'neck' due to ambient fluid (dark color) brought into the jet by the greatly increased strength of vortices near the nozzle exit. The flow structure exhibits an appreciably increased lateral spread than the unforced jet ( figure 11u). The E u of the centreline u (figure 14a) measured at x/D = 1.0 displays the most pronounced peaks of all at f /f a = 1.0 and its harmonics, though the peaks at higher harmonics vanish rapidly from x/D = 1.0 to 4.0. These peaks become very sharp because of the minijet excitation. The ring vortices remain axisymmetric (figure 10a1 − a6), which is corroborated by the spectral phase Φ 12 , about zero over the entire range of f /f a = 0−2.0, between two simultaneously measured hotwire signals u 1 and u 2 at x/D = 1 and y/D = ±0.3 (figure 15a). The Φ 12 is calculated by tan −1 (Q 12 /Co 12 ), where Co 12 and Q 12 represent the cospectrum and quadrature spectrum of u 1 and u 2 , respectively (Zhou et al. 2002).
An even larger lateral spread is achieved for helical forcing at n = 2 (figure 11h) due to the rotating flow structure ( figure 10h1 − h6). The E u measured on the centreline shows less pronounced peaks at f /f a = 1.0 and its harmonics than its counterpart for axisymmetric forcing. This is because there are only two or three injecting minijets at any instant for helical forcing (figure 9h1 − h6), which produce considerably less velocity fluctuations at the jet centre than six minijets (figures 5 and 6). Further, the peaks at f /f a 2.0 disappear at x/D 1. The Φ 12 is about π over a range of frequencies about f /f a = 1.0 (figure 15h), as is expected based on figure 9(h1 − h6).
At n = 5, the jet column wobbles right and left, as indicated by yellow arrows in figure 11f 1, in the flapping plane (the x − y plane) but not in the orthogonal x − z plane ( figure 11f 2), which is symmetrical about y = 0. The lateral spread appears exceeding appreciably that for helical forcing, echoing the considerably improved mixing shown in figure 8. Interestingly, the peaks at f /f a = 1.0 and its harmonics in E u are less pronounced than their counterparts of helical forcing ( figure 14). The peaks of the former decay more rapidly, completely vanished by x/D = 4.0, than the latter. The observation is consistent with the perception that it is the flapping motion, not the large-scale vortices, that plays a predominant role in enhancing mixing in this case. The flapping motion of the jet is characterized by a negative correlation between the two fluctuating streamwise velocities obtained on the opposite side of the jet (Goldschmidt & Bradshaw 1973). Indeed, Φ 12 is about π over a very narrow frequency band about f /f a = 1.0 in the flapping plane but zero over a rather broad range of frequencies in the non-flapping plane ( figure 15f 1, f 2), as observed by . Note that the peak at f /f a = 2.0 is larger than at f /f a = 1.0 for flapping forcing (figure 14f ). This behaviour is ascribed to the flapping motion (figure 11f 1) caused by two separate excitations with a phase shift of π within each excitation cycle (figure 9), which are captured by the hotwire. The combined mode is distinct from all other forcings. Firstly, its spread shown in figure 11(c) is clearly the largest of all, due to the presence of both flapping and helical motions (figures 9c1−c6 and 10c1−c6), internally consistent with the smallest J in figure  8. Secondly, its E u (figure 14c) displays a number of differences from other forcings. The peak at f /f a = 1.0 grows in amplitude from x/D = 1 to 3, while its counterpart for other three forcings all decay quickly. Furthermore, the peaks at the higher harmonics of f /f a = 1.0 decay little for the same range of x/D, in distinct contrast to their counterparts of other forcings where these peaks retreat rapidly. Naturally, compared with other three forcing modes, there are many peaks at the higher harmonics of f /f a = 1.0, especially at x/D = 3. When manipulating main jet using a single unsteady minijet, Perumal & Zhou (2018) made a similar observation, which was ascribed to the use of a small duty cycle. The small duty cycle occurs in every phase of the combined mode (figure 9c1 − c6). Another note is that the peak at f /f a = 2.0 is larger than at f /f a = 1.0 for x/D = 1 and remains very pronounced downstream, similarly to the flapping forcing case (figure 14f ). Thirdly, its Φ 12 (figure 15c) approaches anti-phase at f /f a = 1.0 for both orthogonal planes examined. A rather broad plateau occurs about f /f a = 1.0, where Φ 12 ≈ -0.86 π. This phase shift differs appreciably in value from flapping or helical forcing where Φ 12 ≈ π, and is probably connected to the presence of the oscillating component in this mode. Note that the combined mode produces the nearly anti-phased behavior for all planes through the x axis. However, in flapping forcing, this anti-phased behavior takes place only in the flapping plane.

Momentums impinging upon main jet and jet centre trajectory
Additional insight may be gained into the flow physics of the combined mode by examining the sum of the momentums due to individual injecting minijets impinging upon the main jet and the direction of their resultant momentum. As the averaged mass flux is the same for all minijets, the maximum actuation velocity scales roughly with the inverse of the duty cycle α, as demonstrated in figure 2. As such, the maximum actuation velocity U a,i of the ith minijet is proportional to the product of b i (t) and 1/α i , viz.
where b i (t) is a signal generated by K i (s(t), h(t)), given in eq. 4.1, its values 0 and 1 corresponding to the 'on' and 'off' states, respectively. The time-averaged mass flow rate is the same, 1.2%, for every minijet. The sum of the momentums due to individual injecting minijets impinging upon the main jet is parameterized by By definition, the time-averaged amplitude is the number of active minijets, i.e. 0 for unforced flow and A(t) = N = 6 for actuation, implying a total time-averaged mass-flow rate of 7.2%. On the other hand, the resultant momentum vector of the momentums associated with individual injecting minijets is given by where e i = -(cos θ i , sin θ i ) is a unit vector in the direction of the ith minijet located at angle θ i in the y − z plane. The minus sign refers to an inward velocity towards the jet centre. N is the number of injecting minijets. Figure 16 shows the variation in angle θ(t) of Q(t), with respect to the y axis, and A(t) with time t over two actuation periods T a for the axisymmetric, helical, flapping and combined forcing modes (n = 1, 2, 5 and 11). For axisymmetric forcing, θ(t) is undetermined as Q(t) vanishes identically, as b 1 = b 2 =. . .= b 6 and e i = 0 ( figure 9a1). Therefore, A(t), albeit large (figure 16a2), would not make the jet column oscillate (figure 11a). In the case of helical forcing, θ varies essentially linearly with t, as indicated by the red dashed line (figure 16h1). The stepwise behaviour is caused by the discontinuous on-off actuation b i (t). A large A(t) or one half of the strength of axisymmetric forcing (figure 16a2, f 2) occurs at t/T a = 0.2 and 0.7 for flapping forcing, which correspond to a phase shift of π (figure 16f 1). The behaviours of both A(t) and θ(t) are fully consistent with our understanding of axisymmetric, helical and flapping forcings, thus providing a validation for applying A(t) and θ(t) to describe the forcing on main jet.
For the combined mode, the variations in both A(t) and θ(t) with t/T a are more complicated. Nevertheless, a number of features can be identified. Firstly, after reaching the first maximum as highlighted by a circle, θ(t) decreases, albeit not monotonically, over a duration of 1.3π (figure 16c1), similarly to helical forcing, as indicated by the red arrow. This feature implies a swirling forcing on main jet. Secondly, the maxima of A exceed those of helical forcing, suggesting a stronger swirl actuation. Thirdly, the phase shift between the second local maximum θ(t) and a local minimum θ(t) at t/T a = 0.67 is π, as highlighted by the vertical arrow in figure 16c1, pointing to the signature of flapping forcing. The local maximum A at t/T a = 0.67 (figure 16c2) is associated with a phase change of π (figure 16c1), highlighted by a circle. All the observations suggest that the combined mode is rather unique, featured by both helical and flapping motions, and further by the stronger strengths of flapping and vortical motions than other cases, thus accounting for the largest entrainment and mixing of all.
Evidently, the jet centre r c = (y c , z c ) responds to Q(t i ) (eq. 6.3), and it is plausible to assume the jet centre trajectory to be correlated with Q(t i ). This centre can be Here, r = (y, z) represents the coordinate in the cross-stream plane. For simplicity, the x-dependency of (y c , z c ) will be dropped out hereafter. Over one excitation period with a time step t = T /6, the jet centre takes six positions, i.e. r c (t i ) = (y c (t i ), z c (t i )), t i = i t, i = 0, 1 . . . , 5. We make the most simple assumption that the jet centre displacement ( y c , z c )(t i ) = (y c (t i+1 )−y c (t i ), z c (t i+1 )−z c (t i )) over t is proportional to the actuation momentum, viz.
for i = 0, . . . , 5, where the proportionality constant is set to unity again for simplicity. The discrete time dynamics 6.5 can be considered to a rough discretization of the dynamics equation d(r c )/dt = Q, which describes the jet centre motion under the external momentum. Eq. 6.5 represents 10 equations for 12 unknown jet centre coordinates. The remaining equations are obtained from the observation that the control law is periodic in time and the time-averaged actuation momentum vanishes based on experimental constraints, i.e. the same averaged mass flow through each minijet. Hence, the average jet centre position can be expected to vanish: 5 i=0 r c (t i ) = 0. (6.6) Equations 6.5 and 6.6 constitute 12 linear equations for 12 unknowns, describing the motion of the jet centre over one excitation period. The jet centre dynamics is most easily solved by starting the integration at the origin with x c (t 0 ) = y c (t 0 ) = 0, iteratively computing the positions at t i , i = 1 . . . , 5 with (eq. 6.5) and adding a translation consistent with vanishing averaged jet centre (eq. 6.6). Figure 17 presents the trajectories of the jet centre within each excitation period for the four forcings, which are calculated based on the control laws shown in figure  9 or Eqs. (A 1-A 4) and (6.5-6.6). Apparently, the jet centre vanishes identically for axisymmetric forcing where Q(t i ) ≡ 0 and oscillates along the y−direction between two extremes in the xy plane for flapping forcing where Q(t i ) changes from the positive y−direction at one phase to the negative at next phase or vice versa. Helical forcing, i.e. a uniformly rotating Q vector, leads to a uniformly processing jet and the jet centre moves along a circle around the axis of symmetry. The result conforms to previous reports. Koenig et al. (2016) experimentally investigated the turbulent jet under the helical mode excitation and observed a precessing jet column when the helical structures were spatially amplified in the shear layer. Zhang & Turner (2016) found in a similar experiment that the jet centreline under helical excitation was offset slightly and precessed around the initial axis of the core flow. For the combined forcing, the motion of the jet centre is more complicated. Its trajectory is apparently ellipse-like, suggesting the occurrence of a precession jet. In contrast to helical forcing, the distance of the jet centre from the centre of symmetry varies, along with the separation between the centres of two consecutive phases, suggesting the speed of swirling changes with time. Furthermore, this ellipse-like path indicates an oscillating jet column, a feature of the flapping motion.

Velocity field
The velocity field may provide us with the crucial information on the flow structure. Figure 18 presents the radial profiles of the hot-wire measured U /U j at x/D = 0.05. A number of observations could be made. Firstly, for all controlled cases except helical forcing, the mean velocity profile displays an overshoot at the nozzle exit, as noted by Wu et al. (2018a) who used a single unsteady minijet to manipulate main jet. Andreopoulos & Rodi (1984) made the same experimental observation in case of a circular jet in a cross-stream. They explained that the cross-stream fluid acted like a partial cover over the jet exit, causing the jet flow to bend around and to accelerate so that the velocity of the bent-over jet was somewhat higher than the cross-stream velocity. Secondly, the widths of the mean velocity profiles under control are larger than that of the unforced jet, indicating that the shear layer grows laterally, and the main jet becomes wider right at the nozzle exit. Thirdly, the mean velocity profiles of the axisymmetric and helical forcing at jet exit are almost symmetrical about y/D or z/D = 0, in general with its maximum at the centre ( figure 18a, h). The velocity profile of flapping forcing is also symmetric about y/D and z/D = 0 (figure 18f ), though displaying two peaks off the centre in the x − y or flapping plane and indicating the occurrence of bifurcation. As shown in figure 11(f 1), two consecutive rings are locally connected, forming a zigzag flow structure, as observed in Carlos & Olivier (2002) numerical investigation of a bifurcating jet (their figure 4a). Furthermore, the cross-flow distributions U /U j at x/D = 3 − 5 (figure 18f ) display two peaks in the bifurcation (x − y) plane, while those in the bisection (x − z) plane show only one peak on the centreline. Also, the jet grows slowly in width in the bisection plane. All the features are similar to Lee & Reynolds (1985) data where loudspeakers were used to produce a bifurcating water jet at Re D = 4300. The present data deviate quantitatively from Lee & Reynolds (1985) as a result of distinct actuation techniques, experimental setups and Re D between the two investigations. Finally, the velocity profile under combined forcing mode (figure 18c) is distinct from others and a little tilted at the nozzle exit. Further downstream, U /U j remains asymmetrical about the jet centre but, unlike the flapping mode (figure 18f ), does not show the twin-peak distribution. Wong et al. (2003) produced a precession jet by issuing a jet into a cylindrical chamber with a small axisymmetric inlet at one end and an exit lip at the other. The inlet flow separates at the abrupt inlet expansion and reattaches asymmetrically to the wall of the chamber. Asymmetry of the flow within the chamber causes the reattaching flow to precess around the inner wall of the chamber, resulting in a precessing exit flow. An asymmetric and rotating pressure field is thus established so that the entire flow-field, including the emerging jet, precesses (Nathan et al. 1998). The precessing jet proves to be highly effective in increasing the near-field spreading. Interestingly, the present distributions of U /U j exhibit a similarity to their counterparts of the precession jet (please refer to figure 10 in (Wong et al. 2003)) at the nozzle exit and downstream development.
The radial distributions of u rms /U j ( figure 19) are reasonably symmetrical about y/D or z/D = 0 under control, except under combined forcing. The u rms rises greatly at x/D = 0.05 throughout the main jet, as compared to the uncontrolled jet (figure 19u) where u rms /U j = 0.3% at y/D = 0 and 2.5% at y/D ≈ ±0.43 due to the shear layer. The peak in the shear layer becomes very pronounced. The result suggests a turbulent jet at the nozzle exit, which is internally consistent with flow visualization data (figure 11), and the shear layer instabilities are significantly amplified. The u rms /U j for axisymmetric forcing in the centre region is largest of all, due to the simultaneous injection of six minijets into the main jet, which causes a strong disturbance in the central region ( figure  19a). The u rms /U j distribution displays twin peaks for helical forcing (figure 19h). This is reasonable as helical forcing may produce a hurricane-like helical motion with a centre that is more stable than the surrounding motion. The u rms /U j at x/D = 0.05 under flapping forcing in the x−y plane is larger than in the x−z plane in the shear layer (figure 19f ), as observed by Hussain & Husain (1989) and Zaman (1996). The twin-peak behavior is evident along the z axis due to the flapping motion, but not so along the y axis. The u rms /U j distributions under combined forcing are asymmetrical about the centre at the jet exit (figure 19c), and again resemble the procession jet (Mi & Nathan 2005). Unlike the case under flapping forcing, the u rms /U j profiles under combined forcing display marked twin peaks along the y and z axes, where the right peak is more pronounced than the left, which is probably linked to the helical motion. Furthermore, the u rms /U j peaks at x/D = 3 are substantially higher along the y axis than along the z axis, which remains discernible at x/D = 5. All the features have been observed in the precession jet. Mi & Nathan (2005) investigated the streamwise development of u rms in a procession jet. As shown in their figure 18, the u rms of the precessing jet is asymmetric compared with a non-precession jet and exhibits two peaks, located tangentially 'in front of' and 'behind' the jet centre. The 'front' peak is more pronounced than the 'rear' peak. The similar behaviours between the jet under combined forcing and a precessing jet may suggest that the so-called combined forcing may have produced a precession jet. This suggestion is further corroborated by the downstream development of the centreline mean and fluctuating velocities U cl /U j and u cl,rms /U j presented below.
The variations in U cl /U j , u cl,rms /U j and U j /U cl are presented in figure 20 for various forcing modes as well as the unforced jet. It is worth pointing out that our contraction nozzle is extended by a 47 mm long smooth tube of the same diameter as the nozzle exit D where the minijet assembly is mounted. In spite of this difference, the unforced jet displays the well-known features. Firstly, f 0 D/U j is 0.45, falling in the expected range 0.24 ∼ 0.64 (e.g. Gutmark & Ho (1983), Zhou et al. (2012)). Secondly, U cl /U j and U j /U cl (figure 20a, c) agree both qualitatively and quantitatively with Mi & Nathan (2005) (Re D = 8050) and Seidel et al. (2005) (Re D = 8800) measurements. Following Todde et al. (2009), we may fit the data to U cl /U j = B[(x − x 0 )/D] −1 , where x 0 and B denote the virtual origin and decay constant, respectively. Then, x 0 /D = -0.49 and U cl /U j decays at a rate of x −1 beyond x/D = 5, as reported by previous investigations (e.g. Mi et al. (2001)). Thirdly, the streamwise distribution of u cl,rms exhibits one peak at x/D = 3.5 and another x/D = 7.5 (figure 20b). The former is connected to the breakdown of the primary ring vortices, and the latter is due to early transition to turbulence (Mi et al. 2013). Similar observations were made by Todde et al. (2009) (their figure 6) and by Mi et al. (2013) whose data are included in figure 20 (b). Departures between the present and other's data are not unexpected in view of differences in, inter alia, experimental setup and Re D among the investigations. The potential core length of the unforced jet is approximately 5D, beyond which U cl appears dropping approximately linearly. Note that U cl /U j for axisymmetric forcing exceeds 1.0 given x/D 2, higher than those of the other modes. This is due to an increase in the mass flow rate by 7.2%. A similar observation is made by Seidel et al. (2005) who observed, with 16 minijets blowing, an increase in the centreline mean velocity near the jet exit. Under all control modes, U cl /U j decays rapidly right from the beginning. The minijet actuation reduces U cl /U j significantly at x/D 8 (figure 20a), demonstrating the efficacy of minijet injections. This efficacy can be attributed to the earlier increase in the turbulence levels when the radial injections are present; a large turbulence level (figure 20b) is correlated with a rapid drop in U cl /U j (figure 20a). Interestingly, the variations in u cl,rms /U j for helical and combined forcings are quite similar to each other, growing almost linearly first and then experiencing a small drop before fluctuating slightly around 0.15 (figure 20b). This similarity is ascribed to the common feature of the two control modes, i.e. the swirling motion. However, it is combined forcing that maintains the rapid and linear growth further downstream, up to x/D = 3 and retreats little by x/D = 3.5. It is noted earlier in figure 14c that those peaks at f /f a = 1, 2, 3, 4 and 5 in E u also remain pronounced up to x/D = 3. This coincidence suggests that the excited coherent structures account for the linear growth in u cl,rms /U j (figure 20b) and the rapid decay in U cl /U j (figure 20a). In contrast, the rapid growth of u cl,rms /U j is only up to x/D = 2 for helical forcing and even only up to 1 for axisymmetric and flapping forcings (figure 20b). As such, U cl /U j keeps decaying rapidly and almost linearly until x/D = 4.5 for the combined forcing mode and remains well below other cases further downstream (figure 20a). It is worth pointing out that the precessing jet is also characterized by a substantially faster decay than the non-precessing jet (Mi & Nathan 2005). Jet spreading rate and the overall entrainment rate may be well quantified by the downstream variation of the jet width (e.g. Zhou et al. (2012)). Following Hussain & Husain (1989), we define an equivalent jet width by R eq = [R H R V ] 0.5 , where R H and R V denotes the mean-velocity half-widths in the x − z and x − y planes, respectively. The half-width is defined as the distance between the jet centreline and the location at which U = 0.5U cl . Figure 21 shows the downstream evolution of R V , R H and R eq . While changing little for unforced jet, R V , R H and R eq grow appreciably in controlled jet. Evidently, R eq is the largest for the combined mode, followed by flapping, helical and axisymmetric, though the latter two do not differ much. The results provide additional support for our choice of J as a measure for the mixing efficacy.
To understand further the predominant flow structures under the four forcings, we examine in figure 22 typical instantaneous V -or W -contours, measured using PIV, in the x−y and x−z planes. In the unforced jet, the positive and negative velocity concentrations occur in pair and are mirrored by another pair, though with swapped signs, on the other side of the centreline (figure 22u). Apparently, the two pairs of velocity concentrations are associated with the two vortical structures, as indicated by symbols '+' and '×', of one ring vortex. Axisymmetric forcing leads to the topologically unchanged velocity contours but a significantly increased size in the velocity concentrations which now start to occur in the near proximity of the nozzle exit (figure 22a). Under helical forcing, the vortical structures above and below the centreline are stagger-arranged and their samesigned velocity contours are connected, forming inclined strips alternately in sign (figure 22h). Under flapping forcing (figure 22f 1), the V -contours in the x − y plane exhibit alternate upward and downward motions which are absent in the W -contours in the x − z plane. Some interesting observations can be made for combined forcing. Firstly, the inclined strips of velocity contours, alternately in sign, are seen in both planes ( figure  22c1, c2). Secondly, the velocity concentrations reach in general the maximum contour level of ±0.5, greater than those associated with other forcings (±0.4), indicating a stronger spread/entrainment or mixing. Thirdly, the pattern of velocity contours and their maximum strength persists considerably further downstream, reaching x/D ≈ 5, whilst the upward and downward motions under flapping forcing disappear at x/D ≈ 3, suggesting a prolonged entrainment and more thorough mixing. The observation is again internally consistent with the rapid drop until x/D ≈ 4.5 in the centreline mean velocity decay (figure 20a).
More insight may be gained into the predominant flow structures under control by examining typical instantaneous velocity vectors superimposed with the corresponding iso-contours of streamwise vorticity ω x = ω x D/U j in the y − z plane at x/D = 0.25 (figure 23). Under axisymmetric forcing (figure 23a1, a2), the vectors show the inward or outward motions associated with the ring-like structures, which are axisymmetric and highly repeatable. There are six pairs of counter-rotating ω x concentrations arranged about the centre (figure 23a), apparently generated by the six axisymmetrically placed minijets. The rotational motion under helical forcing is evident and the fluid moves inward along the circumference ( figure 23h1, h2). The core region appears rather stagnant. The phase of injecting minijets is clockwise incremented by 60 • (figure 9h), producing a corkscrew type of structure (figure 10h1 − h6), as shown by Koch et al. (1989). Figure  23(f 1, f 2) shows the cross-flow motion that switches from one direction at one moment to the opposite at another under flapping forcing, accompanied by one pair of counterrotating ω x concentrations, as shown by . The velocity vectors in Figure 23(c1, c2) exhibit the clockwise rotational motion under combined forcing. The area of rotational motion, as indicated by the arrows in figure 23c1 is in general larger than that under helical forcing (figure 23h1 − h2). Note that ambient fluid may be entrained into the jet core area from various circumferential locations, e.g. the upper left and lower right corners as indicated by the elliptic contours in figure 23c1, while under helical forcing ambient fluid comes into the core area largely from only one location as highlighted by the elliptic contours ( figure 23h1, h2). Furthermore, there are many vorticity concentrations of both signs in figure 23(c1 − c2). The core area is dominated by the vorticity concentrations of negative sign, while the region surrounding the core is populated with those of both signs. This is very different from flapping forcing where there is only one pair of opposite-signed vorticity concentrations. This is also markedly different from helical forcing ( figure 23h1, h2) where the cross-sectional plane is characterized by the vorticity concentrations of a single sign. The observations reconfirm that the combined forcing mode is associated with a much better mixing and furthermore probably also small-scale mixing.

Insight into the 3D flow structure
To gain insight into the three-dimensional (3D) flow structure, we deploy a linear stochastic estimation method to reconstruct the predominant flow structure from the PIV data, captured in the y − z plane at x/D = 0.25 with a total of 200 images. This technique is introduced in detail by e.g. Adrian & Moin (1988) and is briefly described below. Thus, the periodic flow response may be given by q(y, z, t) = q(y, z) + q 1 (y, z)cosφ(t) + q 2 (y, z)sinφ(t), φ = 2πf a t. (6.9) The residual of this phase-averaged flow q consists of higher harmonics and an uncorrelated stochastic contribution. Note the actuation commands b i (t) are the functions of the pointer (cosφ, sinφ). Hence, the temporal Fourier component optimally represents the actuation response and no flow-intrinsic phase needs to be constructed. The iso-surfaces of reconstructed V and W are presented in figure 24 for four control modes. The flow structures of axisymmetric, helical and flapping forcing modes show excellent agreement with the control laws (figure 9). Under axisymmetric forcing (figure 24a1 − a2), V and W indicate clearly that jet fluids on the two sides of the centreline move either inward or outward simultaneously (figure 24a). Under helical forcing (figure 9h1 − h6), the helical motion is evident (figure 24h). For flapping forcing, as indicated by the iso-surfaces of W , one sector of fluid moves in one direction for one instant and the adjacent sectors move in the opposite direction ( figure 24f 2). On the other hand, the iso-surfaces of V is anti-symmetrical about the centreline (figure 24f 1). The observations are fully consistent with the characteristics of the flapping motion documented in e.g. . The results provide a validation for the presently reconstructed V and W . For the combined mode (figure 24c), the flow structure appears much more complicated. However, the swirling motion is still discernible. Furthermore, the iso-surfaces of V and W show unequivocally the occurrence of many more small parcels of fluids, suggesting a much better jet mixing than all other flow modes, which is fully consistent with the finding from figure 8 as well as figure 23(c).

Conclusions and outlook
An artificial intelligence (AI) control system has been developed for the control of turbulence. The system consists of a control plant, a sensing unit, an actuation unit, and a control or 'thinking' unit. A round jet (control plant) is manipulated to illustrate the potential of this system. Two hot-wires are deployed for sensing online the information on the turbulent jet. The control unit deploys a linear genetic programing, and six independent unsteady radial minijets placed around the nozzle prior to the issue of main jet for executing control laws generated from the control unit. The search space for control laws is extremely large, including the minijet number N , geometric configuration, frequency f a , duty cycle α and phase shift ϕ ij (i, j = 1, 2, . . . , 6) between minijets. The target is to maximize the decay rate K of the jet centreline mean velocity, which is correlated with the mixing rate.
It has been demonstrated that the AI control system can learn automatically how to optimize the spatially distributed actuators and thus the turbulent jet for the targeted cost. Like virtually all other control strategies of nonlinear dynamics, AI control solutions do not come with a proof of global optimality. Yet, the results for jet mixing optimization demonstrate several highly desirable features. First, AI control has identified a few typical and well-known control laws, i.e. axisymmetric, helical and flapping forcings, in its learning process in the order of increased control performances, and eventually converged to an unexpected spatio-temporal forcing, referred to as the combined mode, which has never been reported previously. The learning time of 1100 individuals or 2 hours windtunnel testing is remarkably short for such a complex solution. It is worth pointing out that the presently developed AI control method may not suit for numerical simulation where faster learners, e.g. the reduced-order-model-based models, would be required (once the winning control mechanism is known) to reduce significantly the testing time for N-S-based simulations. Inspirations may be gained from control studies of a turbulent boundary layer (Sasaki et al. 2020) or turbulent separation (Nair et al. 2019). Second, the AI-learned combined mode is reproducible with other initial generations. The control laws may analytically differ but produce almost identical actuation commands. Third, the parameters of the underlying linear genetic programming are taken verbatim from Duriez et al. (2016) and were already proven useful in many other experiments. No sensitive dependence on the parameters has been observed so far and AI control can be expected to yield near-optimal results in its first application to a new plant. Finally, the unique advantage of the AI control over conventional techniques, be open or closed-loop and linear or nonlinear, is its capability to find the apparently global optimum solution when the search space for control laws is extremely large. This advantage will be lost in case of a single actuator involving few control parameters such as frequency and duty cycle (e.g. Fan et al. (2017)).
The control landscape is studied by examining the feature vectors of generations n = 1, 2, 5 and 11. Several interesting findings are made. Firstly, jet mixing benefits from the increasing asymmetry of forcing in the learning process. Secondly, the feature vectors appear randomly distributed at n = 1 but evolve gradually to three distinct regions with increasing n (figure 12), characterized by axisymmetric, helical and flapping forcings, respectively. While the region of axisymmetric forcing contracts to one single point, the other two shrink essentially to two curves in the evolution process. Intriguingly, the converged globally optimized control law at n = 11 takes place in the region of helical forcing which is now populated with either helical or combined mode. Thirdly, the feature coordinate γ 1 is found to be correlated with the degree of asymmetric forcing, while γ 2 is linked to the number of simultaneously injecting minijets. The best control law is characterized by the largest γ 2 for axisymmetric forcing and largest γ 1 for the other two forcings. However, the globally optimized control law (combined mode) is associated with neither the largest γ 1 nor largest γ 2 .
This combined forcing has produced a novel turbulent flow structure characterized by strong oscillation and swirling motions, along with the generation of mushroom-like structures, all acting to enhance jet mixing. As a result, the combined mode vastly outperforms the well-known optimal benchmark forcings, increasing the entrainment or mixing rate by 54%, 47% and 28% compared with the axisymmetric, helical and flapping forcings, respectively. Extensive measurements are conducted in three orthogonal planes to understand this novel flow structure, which is compared with the flow structures under the benchmark axisymmetric, helical and flapping forcings. It has been found that this flow structure is characterized by a number of features, including (1) the helical motion, (2) 3D oscillating jet column, (3) strong coherent structures, (4) asymmetrical cross-flow distributions of mean and fluctuating velocities, the latter displaying twin-peak behavior in every plane through the x-axis, (5) spiral behavior of the swirling jet centre, (6) changing speed of swirling, (7) spectral phase shift by π over a considerable frequency band about f /f a = 1 between two simultaneously captured hotwire signals placed at (x/D, y/D) = (1, ± 0.3), and (8) many vorticity concentrations of opposite signs over the entire cross-sectional jet plane, which is distinct from the flow structures under other forcings and suggests better and smaller scale mixing. It is further found that features (1) through (5) resemble those of a precession jet generated by a passive device.
We expect that AI control will be commonly applied to discover the unknown winning nonlinear actuation mechanism of multi-input multi-output flow control experiments in the very foreseeable future. Conventional model-based or model-free control design may be then deployed to refine and optimize the AI-based actuation mechanism, provided the control law is sufficiently simple. AI control may be improved in numerous aspects. Examples include a human-interpretable control law, an increased learning speed in experiment, robustness against varying operating conditions and the inclusion of prior knowledge and expectations of control laws. One can safely assume that AI will be an essential tool in future turbulence control applications, just as AI is indispensable in robotics now.
YZ wishes to acknowledge support given to him from NSFC through grants 11632006, 91752109 and 91952204. This work is supported by the French National Research Agency (ANR) via the grants ANR-11-IDEX-0003-02 (iCODE), 'ACTIV ROAD' and 'FlowCon', and by the OpenLab Fluidics consortium (Fluidics@poitiers) of PSA Peugeot-Citroë n and Institute Pprime.