Artificial Neural Networks and Guided Gene Expression Programming to
  Predict Wall Pressure Spectra Beneath Turbulent Boundary Layers
This study evaluates the efficacy of two machine learning (ML) techniques,
namely artificial neural networks (ANN) and gene expression programming (GEP)
that use data-driven modeling to predict wall pressure spectra (WPS) underneath
turbulent boundary layers. Different datasets of WPS from experiments and
high-fidelity numerical simulations covering a wide range of pressure gradients
and Reynolds numbers are considered. For both ML methods, an optimal
hyperparameter environment is identified that yields accurate predictions. ANN
is observed to be faster and more accurate than GEP with an order of magnitude
lower training time and logarithmic mean squared error ($lMSE$), despite a
higher memory consumption. Novel training schemes are devised to address the
shortcomings of GEP. These include (a) ANN-assisted GEP to reduce the noise in
the training data, (b) exploiting the low and high-frequency trends to guide
the GEP search, and (c) a stepped training strategy where the chromosomes are
first trained on the canonical datasets followed by the datasets with complex
features. When compared to the baseline scheme, these training strategies
accelerated convergence and resulted in models with superior accuracy ($\approx
30\%$ reduction in the median $lMSE$) and higher reliability ($\approx 75\%$
reduction in the spread of $lMSE$ in the interquartile range). The final GEP
models captured the complex trends of WPS across varying flow conditions and
pressure gradients, surpassing the accuracy of Goody's model.
Nachiketa Narayan Kurhade, Nagabhushana Rao Vadlamani, Akash Haridas
Artificial Neural Networks and Guided Gene Expression Programming to Predict Wall Pressure Spectra Beneath Turbulent Boundary Layers

###### Abstract

This study evaluates the efficacy of two machine learning (ML) techniques, namely artificial neural networks (ANN) and gene expression programming (GEP) that use data-driven modeling to predict wall pressure spectra (WPS) underneath turbulent boundary layers. Different datasets of WPS from experiments and high-fidelity numerical simulations covering a wide range of pressure gradients and Reynolds numbers are considered. For both ML methods, an optimal hyperparameter environment is identified that yields accurate predictions. ANN is observed to be faster and more accurate than GEP with an order of magnitude lower training time and logarithmic mean squared error (\(lMSE\)), despite a higher memory consumption. Novel training schemes are devised to address the shortcomings of GEP. These include (a) ANN-assisted GEP to reduce the noise in the training data, (b) exploiting the low and high-frequency trends to guide the GEP search, and (c) a stepped training strategy where the chromosomes are first trained on the canonical datasets followed by the datasets with complex features. When compared to the baseline scheme, these training strategies accelerated convergence and resulted in models with superior accuracy (\(\approx 30\%\) reduction in the median \(lMSE\)) and higher reliability (\(\approx 75\%\) reduction in the spread of \(lMSE\) in the interquartile range). The final GEP models captured the complex trends of WPS across varying flow conditions and pressure gradients, surpassing the accuracy of Goody's model.

+
Footnote †: preprint: Accelerated GEP to Predict Wall Pressure Spectra beneath Turbulent Boundary Layers

## I Introduction

Turbulent boundary layers (TBLs) developing over the surfaces induce wall-pressure fluctuations. These dynamic loads are highly undesirable as they increase structural vibrations and noise. Fatigue failure due to aeroacoustic loads is a critical problem in wind turbine and gas turbine blades [1; 2; 3; 4; 5]. Aerodynamic noise due to attached or separated TBLs contributes significantly towards medium and high-frequency sound pressure levels (SPL) of the cabin [6; 7]. Acoustic loads due to wall-pressure fluctuations are much more severe over the skin panels of high-speed vehicles. The amplitude and frequency of these loads can be significantly higher in certain regions which can potentially damage the vehicle.

Estimating the wall-pressure fluctuations beneath TBLs is hence crucial to facilitate the structural and aerodynamic design process. These fluctuations can either be estimated from high-fidelity eddy-resolving simulations or can be directly measured from experiments. Both these approaches are prohibitively expensive. Although eddy-resolving methods like Direct Numerical Simulations (DNS) [8; 9; 10] and Large Eddy Simulations (LES) [11; 12; 13] can accurately quantify the spatiotemporal variation of acoustic loads, they are computationally expensive to cater to the highly iterative design process. On the other hand, most of the experiments in the literature are confined to measuring wall pressure spectra (WPS) of TBLs developing over simple geometries like a flat plate subjected to zero, favorable, or adverse pressure gradients [14; 15; 16; 17] (ZPG, FPG, and APG).

One of the widely used approaches to overcome these limitations is to couple the empirical models of WPS with the mean boundary layer characteristics obtained from the Reynolds Averaged Navier-Stokes (RANS) simulations to predict acoustic loads [19]. These models rely on scaling laws based on the inner and outer regions of the boundary layer to estimate the sound pressure levels. Several empirical models for WPS were developed based on the theoretical foundations of Lighthill [20; 21] and Kraichnan [22]. Blake [23] integrates the Poisson equation for pressure fluctuation beneath a TBL. As noted by Grasso _et al._[24], the theoretical model involves repeated integration in multiple dimensions, and the computational cost increases exponentially with the number of dimensions. An alternate approach involves correlating the non-dimensionalized SPL distribution with boundary layer parameters like thickness, wall-shear stress, freestream velocity, etc. This approach, as reflected in the works of Chase [25] and Howe [26] has an evolutionary advantage as the subsequent models continued to modify the same baseline empirical structure to account for newer observations. Goody [27] further modified the Chase-Howe model to include Reynolds number effects and capture the frequency dependence. Goody's model is designed to predict the WPS of TBLs under ZPG with a \(\omega^{2}\) growth at low frequencies (in line with Kraichnan-Phillips' theorem), \(\omega^{-1}\) decay in the inertial range (in line with Bradshaw prediction [28]) and \(\omega^{-5}\) decay at the high-frequencies (in line with Blake's observations [23]). Subsequent works of Kamruzzaman _et al._[29], Rozenberg _et al._[19], Hu [30] and Lee [31] extended Goody's model to account for pressure gradients. Ritos _et al._[32] recognised that the existing models fall short in predicting the WPS beneath supersonic and hypersonic flat plate TBLs under ZPG conditions and modified Goody's model to account for the compressibility effects. Wind tunnel experiments of Thomson and Rocha [33] measured wall pressure fluctuations of TBL on a flat plate under FPG. They further tuned the parameters of the universal spectrum model (which is inspired by Goody's model) to incorporate FPG effects and improve predictions at high-frequencies. Thomson and Rocha [34] compare the performance of Goody's model on both the wind
tunnel data and flight test data. They highlight that the model accurately predicts the WPS of the former while underpredicts the latter and accordingly updated Goody's model to accurately fit the flight test data. Accounting for the different flow conditions and pressure gradients, Fritsch _et al._[35] developed a semi-empirical WPS model using numerical optimization algorithms. The predictions are shown to be sensitive to the values of constants involved and the optimized set is bound to evolve with an ever-evolving dataset. All these empirical models, however, were developed based on the respective datasets and hence their accuracy suffers when employed for 'extrapolated' flow conditions.

With increasing computational power and the databases of wall pressure spectra data published in public domains, robust wall pressure spectra models can be developed by exploiting the entire database. Additionally one can identify the range of TBL parameters covered by these databases [36]. These insights can then be used by modern machine learning (ML) algorithms to solve this regression problem. Haridas and Vadlamani [37] and Dominique _et al._[36] used artificial neural networks (ANN) for this purpose and obtained better predictions than any of the empirical models discussed above. The study conducted by Haridas and Vadlamani [37] used subsonic and supersonic datasets to train their ANN model. This approach also enabled Dominique _et al._ to quantify the confidence in their predictions and identify the input space where more data will be helpful in improving the predictions. ANN being the universal function approximator [38] is indeed best suitable for the job as it can be retrained with little effort on the new datasets. The downside with ANN is twofold: (a) It results in a highly non-linear function that is difficult to translate to a conventional functional form and (b) It provides little insight into the underlying physics of the problem [37]. On the other hand, Gene expression programming (GEP), a symbolic machine learning algorithm, operates on the input variables through generations to evolve into an analytical expression that tries to approximate provided data. Dominique _et al._[39] showed the capabilities of this technique to predict the WPS models that are close to the established empirical correlations in addition to discovering the new dependencies that were not considered earlier. GEP however suffers from converge issues. Also, for a given set of input features, the algorithm produces different analytical expressions, although their structures and predictions can be very similar.

The objective of the present work is multifold. We first compare the strength and weaknesses of ANN and GEP (with an optimized hyperparameter environment) against the same input data assessing the time taken to converge, and the computational resources utilized. Subsequently, we propose strategies (like ANN-assisted training incorporating physical insights) to accelerate the training of the GEP algorithm and mitigate its convergence issues. Finally, we demonstrate the efficacy of the training strategies and the robustness of the models predicted using the guided GEP approach.

The manuscript is structured as follows: We introduce the regression problem and the relevant WPS models in Section II. Section III presents a brief description and analysis of the input dataset and a short overview of ANN and GEP. We present the optimum hyperparameter environment and compare the computational efficacy of ANN and GEP in Section IV. Strategies to accelerate GEP training to mitigate the convergence issues are addressed in Section V. Lastly, Section VI concludes the key findings of the study.

## II Characterizing Turbulent Boundary Layer and Wall Pressure Spectrum

As mentioned in the introduction, the empirical models of WPS rely on the mean boundary layer characteristics to predict acoustic loads. All the TBLs have a negligible zero pressure gradient in the wall-normal direction (\(\partial_{x}p=0\)) but can be subjected to a considerable gradient (\(\partial_{x}p\)) in the streamwise direction. Figure 1 sketches a typical velocity profile \(U(y)\) at a distance \(x\) downstream of the leading edge for a general case where the equilibrium velocity at the boundary layer edge changes with \(x\) as \(U_{\epsilon}(x)\). The local boundary layer thickness (\(\delta(x)\)) is usually defined as the wall-normal distance at which \(U=0.99U_{\epsilon}(x)\). In this study, following Dominique _et al._, we use an alternate definition for boundary layer thickness based on the pseudo-velocity [9]\(U^{*}(x,y)=-\int_{0}^{y}\Omega_{z}(x,\xi)\,d\xi\) to account for the effect of pressure gradients on the outer flow velocity. Here, \(\Omega_{z}\) is the spanwise component of the mean vorticity \(\Omega=\nabla\times\widetilde{U}\), \(\overline{U}\) being the mean velocity vector and \(\xi\) is the length in wall-normal direction.

At the same point \(x\) downstream of the leading edge, a probe can be placed on the wall to record the pressure fluctuations due to TBL. The unsteady pressure signals \(p(\tau)\) can be decomposed into its mean and fluctuating components as follows:

\[p(x,\tau)=\overline{p}(x,\tau)+p^{\prime}(x,\tau) \tag{1}\]

The point-wise power spectral density (PSD) \(\Phi_{pp}\) of these pressure fluctuations is obtained from the following equation:

\[\Phi_{pp}(x,\omega)=\int_{-\infty}^{+\infty}R(x,\tau)e^{-i\omega t}\,d\tau \tag{2}\]

Here, \(R\) is the auto-correlation function defined as \(R(x,\tau)=\int_{0}^{\infty}p^{\prime}(x,t)p^{\prime}(x,t+\tau)\,dt\), \(\Phi_{pp}(x,\omega)\) is the PSD obtained from

Figure 1: Typical flat plate boundary layer
the fast Fourier transform of \(R\), \(\omega\) is the angular frequency and \(i=\sqrt{-1}\). Several empirical models in the literature express the empirical correlation for \(\Phi_{pp}\) in terms of the mean boundary layer parameters:

\[\Phi_{pp}(\omega)=f(\omega,\delta,\delta^{*},\theta,U_{e},\rho,\nu,\tau_{w}, \partial_{x}p,c,\Pi) \tag{3}\]

Where \(\delta^{*}=\int_{0}^{\delta}(1-U(x,y)/U_{e})\,dy\) is the displacement thickness, \(\theta=\int_{0}^{\delta}U(x,y)/U_{e}(1-U(x,y)/U_{e})\,dy\) represents the momentum thickness of the TBL, and \(\tau_{w}\) is the wall shear stress at location \(x\). \(\rho\), \(\nu\), and \(c\) represent the density, dynamic viscosity, and speed of sound in the freestream respectively. A TBL typically consists of four distinct regimes in the wall-normal direction: viscous sublayer closer to the wall, intermediate buffer layer, logarithmic region influenced by large turbulent scales, and the wake region. The deviation of the velocity profile from the log law in the wake region is characterized by the wake strength parameter [40], \(\Pi\). The number of independent variables in Eq. 3 can be reduced using the Buckingham-Pi theorem. Following the work of Dominique _et al._[36], Eq. 3 is rewritten in terms of the non-dimensionalize variables as:

\[\widetilde{\Phi}_{pp}=f(\widetilde{\omega},\Delta,H,M,\Pi,C_{f},R_{T},\beta) \tag{4}\]

Where,

\[\text{Dimensionless PSD:}\;\widetilde{\Phi}_{pp}=(\Phi_{pp}U_{e})/(\tau_{w}^{2}\delta)\]

\[\text{Dimensionless angular frequency:}\;\widetilde{\omega}=(\omega\delta)/U_{e}\]

\[\text{Zagarola-Smits's parameter}\]

\[\text{Shape factor:}\;H=\delta^{*}/\theta\]

\[\text{Mach number:}\;M=U_{e}/c\]

\[\text{Wake strength parameter}\]

\[\text{Friction coefficient:}\;C_{f}=\tau_{w}/(\rho U_{e}^{2})\]

\[\text{Outer-to-inner-layer timescale ratio:}\;R_{T}=(\delta/U_{e})/(\nu/U_{\tau}^{2})\]

\[\text{Clauser parameter:}\;\beta=(\theta/(\rho U_{\tau}^{2}))\partial_{x}p\]

### Goody's Model

One of the earliest empirical models, the Chase-Howe model [25; 26], proposed a correlation of PSD (\(\Phi_{pp}\)) as a sole function of angular frequency \(\omega\). It however fails to predict the high frequency \(\widetilde{\omega}^{-5}\) drop observed in the wall pressure spectra of TBLs under ZPG. Goody [27] proposed the following changes to the Chase-Howe model: (a) Since the largest coherent structures in TBL are of the order of boundary layer thickness, \(\delta^{*}\) is replaced with \(\delta\) based scaling (b) Ratio of outer-layer to inner-layer timescales, \(R_{T}\) (a different form of Reynolds number) is incorporated into the formulation to control the extent of inertial subrange. Goody's model, given in Eq. 5, accurately predicts ZPG turbulent boundary layers that are homogeneous in the spanwise direction. Hence it is considered as the baseline semi-empirical model although it does not account for pressure gradients.

\[\frac{\Phi_{pp}U_{e}}{\tau_{w}^{2}\delta}=\frac{3(\omega\delta/U_{e})^{2}}{[( \omega\delta/U_{e})^{0.75}+0.5]^{3.7}+[1.1R_{T}^{-0.57}(\omega\delta/U_{e})]^{ 7}} \tag{5}\]

### Dominique's GEP Model

Dominique _et al._[39] trained the gene expression programming algorithm on the datasets described in Section III.1 to come up with the following WPS model:

\[\widetilde{\Phi}_{pp}=\frac{\left(5.41+C_{f}\left(\beta+1\right)^{5,41}\right) \widetilde{\omega}}{\widetilde{\omega}^{2}+\widetilde{\omega}+(\beta+1)M+( \widetilde{\omega}+3.6)\frac{\widetilde{\omega}^{4.76}}{C_{f}R_{T}^{5.83}}} \tag{6}\]

Although derived from scratch, Dominique's GEP model reasonably captures the frequency dependencies. Their model shows a \(\widetilde{\omega}^{1}\) dependence at low frequencies, in contrast to the \(\widetilde{\omega}^{2}\) dependence proposed by Kraichnan [22]. This is considered a strong trait for modeling realistic scenarios using GEP where the data does not always agree with the theoretical models.

## III Methodology

### Data Collection

The dataset used in the present work stems from the study by Dominique _et al._[36] (also referred to as the von Karman Institute (VKI) team henceforth). It consists of the variation of \(\widetilde{\Phi}_{pp}\) with \(\widetilde{\omega}\) at subsonic Mach numbers, in addition to the mean TBL characteristics (reported in Eq. 4), collected from the following sources:

* Experiments of Salze _et al._[14] on a flat plate under ZPG, APG, and FPG.
* High-fidelity computations of Deuse and Sandberg [9], Hao _et al._[10] and, Christophe _et al._[12] on the configuration of the flow over a controlled diffusion (CD) airfoil, where the TBL spectra under APG and mild FPG are recorded at different stations.

Figure 2 plots all the 117 datasets, each of which has 500 logarithmically spaced points resampled from the experiments listed above. A machine learning algorithm is otherwise oblivious to the different experiments from which the dataset has been extracted. Rather, it solves the regression problem modeled using Eq. 4 treating the input as a set of independent data points with respective features and labels.

### Data analysis

A closer look at Fig. 2 shows the dominance of the APG dataset over the FPG and ZPG. It is essential to design ML
algorithms to handle this skewness, else it treats the underrepresented cases as outliers and the model predictions suffer. A common practice to balance out skewed databases is to identify similarity clusters and assign the respective weights. One of the widely used methods involves clustering based on the Euclidean distance, \(d=\sqrt{\sum_{i=1}^{n}(P_{i}-Q_{i})^{2}}\) between two vectors (say \(\overline{P}\) and \(\overline{Q}\)) which represent the independent flow parameters listed in Eq. 3.

In the current study, we demonstrate the similarity by comparing the trends of \(\widetilde{\Phi}_{pp}\) vs \(\widetilde{\omega}\) across datasets. A typical WPS, shown in Fig. 3, comprises three distinct regimes: a low-frequency region where \(\widetilde{\Phi}_{pp}\) rises monotonically with \(\widetilde{\omega}\) to reach its peak value, an inertial subrange at mid-frequencies (typically noticeable at high Reynolds numbers) and a \(\widetilde{\omega}^{-5}\) roll-off at higher frequencies due to dissipation. However, not all the datasets in Fig. 2 show this exact trend. The range of frequencies encompassing each of the regimes also differs across datasets. Hence, the following slope vector \(\overline{m}\) is extracted for each dataset of length \(n\):

\[m_{i}=\frac{\widetilde{\Phi}_{pp,i+1}(dB)-\widetilde{\Phi}_{pp,i}(dB)}{ \widetilde{\omega}_{i+1}-\widetilde{\omega}_{i}},i=[1,n-1] \tag{7}\]

It should be noted that the raw WPS dataset is noisy and hence should be smoothed out to extract meaningful trends. The present work used the predictions of an ANN model for the same. Since all datasets have an equal number of points, \(\overline{m}\) also captures the extent of respective regimes. We can now compute the cosine similarity, which results in a value between [-1,1], to quantify the closeness between two slope vectors as follows:

\[CS=cos(\overline{m}_{i},\overline{m}_{j})=\frac{\overline{m}_{i}\cdot\overline {m}_{j}}{||\overline{m}_{i}||\times||\overline{m}_{j}||},(i,j)=[1,117] \tag{8}\]

In a 3D space, \(CS=1\) represents similar vectors, \(CS=0\) indicates orthogonal vectors and \(CS=-1\) represents vectors pointing in opposite directions. In a multidimensional space such as ours, cosine similarity is merely representative of how similar trends are to each other; from identical to very dissimilar respectively. Figure 4 plots the heatmap of the resulting cosine similarity matrix. \(CS\approx 1\) within the boxes aligned along the diagonal indicate that the trends within a given experiment are similar. Bands of \(CS\approx 1\) (highlighted using solid arrows) represent similar WPS trends observed across the experiments with different flow conditions. In contrast, bands of \(CS\approx-1\) (highlighted using dashed arrows) represent dissimilar trends. For example, the inset plots in Fig. 4 illustrate the similarity of WPS extracted at points A, B, and C. It is apparent that the WPS at B and A are similar with a \(CS\approx 0.92\) while the WPS at C is dissimilar to that of A with a \(CS\approx-0.54\). This similarity matrix helps in identifying clusters that can be used to assign weights to the datasets and potentially improve model predictions. This aspect is further discussed in Section IV.2.

### Artificial neural networks

Artificial neural network (ANNs) is a powerful function approximator that is commonly used to fit non-linear data using supervised learning. A brief description of ANN is given here and readers can refer to the works of Goodfellow _et al._[43] for further details. Figure 5 illustrates a typical ANN architecture comprising an input layer (which can be normalized for improved scaling of data), several hidden layers (to accurately

Figure 3: Low and high-frequency trends in TBL wall pressure spectrum, black dashed lines as eye guides

Figure 2: Wall pressure spectra collected by the von Karman Institute(VKI) team rescaled to Goody’s scales
capture the complexity of the non-linear function that fits the data), and an output layer. The hidden layers consist of neurons. Each neuron generates a weighted sum \(\Sigma\) from its inputs (\(\overline{w}\cdot\overline{x}\)) with respective biases \(b\). A predefined activation function (\(\sigma\), usually a non-linear function) operates on this weighted sum \(\Sigma\) to produce an output \(y\). The bias is a constant value that is added to the weighted sum of the inputs to shift the activation function. It is mathematically expressed as \(y=\sigma\Sigma=\sigma(\overline{w}\cdot\overline{x}+b)\). The error between the expected output and the prediction is minimized by computing the gradient of the objective function with respect to the weights and biases of each neuron through backpropagation [44]. This helps in adjusting the weights and biases of the individual neurons, resulting in a potentially improved model.

The current study considers the data used by the VKI team and hence the optimal hyperparameters from their study have been used to build the ANN model. The neural network has three hidden layers of ten neurons each in a fully-connected feed-forward fashion. These hidden layers are preceded by a normalization layer. In this layer, the input vector \(\overline{x}=[\widehat{\mathbf{\omega}},\Delta,H,Ma,\Pi,C_{f},R_{T},\beta]\) from Eq. 4 is normalized using the mean and standard deviation of each feature to accelerate the learning process [45]. The output layer of the model is made up of a single neuron that predicts the value, \(y=10log_{10}(\Phi_{pp}U_{e}/\tau_{w}^{2}\delta)\). The training is carried out using

Figure 4: Cosine similarity matrix of WPS trends observed in collected datasets. Christophe:1-78, Salze:79-88, Deuse: 89-101: Hao: 102-117

Figure 5: Example of an Artificial neural network build
Nadam optimizer [46] and Selu activation function [47] at a learning rate of 0.0001 with a batch size of 32 randomly chosen samples. The dataset is randomly split into training data and validation data with an 80:20 split respectively. Logarithmic mean squared error (_IMSE_), defined in Eq. 9, is used as the objective function. Training of the ANN is driven by the training loss while the early stopping is identified by monitoring the validation loss.

\[\textit{IMSE}=\frac{1}{N}\sum_{i=1}^{N}W_{i}(10log_{10}(y_{i}^{ expected})-10log_{10}(y_{i}^{predicted}))^{2} \tag{9}\]

### Gene expression programming (GEP)

Like a genetic algorithm GEP discovers mathematical expressions to best represent a dataset. We present a brief summary of the method here and the readers can refer to the article by Candida F. [48] for a detailed discussion on GEP. Figure 6 and Table 1 illustrate the mainstream GEP algorithm, the entities, and the operations involved. The algorithm starts with constructing several fixed-length strings called genes using combinations of the functions and terminals. Part of such a gene (ORF) is interpreted as an expression tree (ET) which is a graphical representation of a mathematical expression. Chromosomes (strings with multiple genes) with linking functions are used to improve upon the complexity captured by the mathematical expression derived from a single gene. The fitness of these chromosomes is then evaluated based on their ability to match the predictions of a given dataset. Chromosomes with desired fitness are retained, reproduced, and modified through genetic operators to result in different expressions altogether. Well-performing combinations are thus retained and mutated over multiple generations until the one with desired fitness is discovered.
\begin{table}
\begin{tabular}{l|c}
**Term** & **Description** \\ \hline Functional set (F) & Set of operators that a GEP inferred mathematical expression can be built out of. \\ \hline Terminal set (T) & Set of variables that a GEP inferred mathematical expression can be built out of. \\ \hline Gene & Fixed length string that is a combination of members from the functional and the terminal set. \\ \hline Head & Part of the gene that is a combination of entities from the functional set and the terminal set. \\ \hline Tail & Part of the gene made exclusively out of the terminal set to ensure that enough terminals are available to result in a valid mathematical expression. \\ \hline Open reading frame (ORF) & Part of the gene that can be read into an expression tree. \\ \hline Expression tree (ET) & Plane tree structure with each node as a member of the functional set and each leaf as a member of the terminal set that can be read into a mathematical expression. \\ \hline Chromosome/Individual & Fixed length string that consists multiple genes. \\ \hline Sub Expression Tree (Sub - ET) & Expression tree that is derived from a gene present in a chromosome. \\ \hline Prescribed Linking function & Function prescribed at the start of GEP training that combines Sub - ETs to produce the final mathematical expression from a chromosome. \\ \hline Fitness & Performance of GEP individuals quantified through objective functions. \\ \hline Selection & Process of selecting parents for the next generation from a pool of individuals while discarding the rest. \\ \hline Reproduction & Process of filling up the population back to its original volume by deriving daughters from parent individuals. \\ \hline Genetic operators & Operators that genetically modify chromosomes. \\ \hline Evolution & Process of performing selection and reproduction on a population while retaining the desired genetic characteristics identified through fitness. \\ \hline Generation & Pool of individuals that undergoes fitness evaluation at any point of time through GEP evolution. \\ \hline \end{tabular}
\end{table}
Table 1: GEP Terminology

Figure 6: (a) GEP algorithm flowchart with least-square based local optimization. (b) Entities involved in GEP
\[Fit=100\times\sqrt{\alpha\frac{lMSE}{A_{log}^{2}}+(1-\alpha)\frac{MSE}{A_{lin}^{2}}} \tag{10}\]

where,

\[MSE=\frac{1}{N}\sum_{i=1}^{N}W_{i}(y_{i}^{expected}-y_{i}^{predicted})^{2} \tag{11}\]

Here, \(A_{log}\) and \(A_{lin}\) represent the maximum amplitude of training points in the respective scales. \(\alpha\) is a weight assigned to the objectives, which is set to 0.5 thus giving equal weight to both contributions. The idea to have such a multi-objective function is derived from the works of Dominique _et al._[39]. They postulated that while the _MSE_ penalizes the error in the amplitude of the spectrum, the _lMSE_ captures error in the trends at low and high frequencies effectively.

A set of experiments were conducted for 100 individuals of different head lengths [3, 4, 5, 6], over 1001 generations, and with optimization probability 0.005 (every generation is being optimized, but very few individuals are optimized in each generation). Figure 8 plots the errors in exponent predictions (which is defined as \(|Exponent_{expected}-Exponent_{observed}|\)) with \(Fit\). It is evident that the \(R_{T}\) and high-frequency errors tend to drop with a decrease in \(Fit\), implying a closer agreement with Goody's model. In almost all the cases, \(R_{T}\) error

Figure 8: Exponent error observed in GEP solution while predicting the Goody’s model

Figure 7: Training dataset for hyperparameter study of GEP (a) With added noise (b) Without any noise along with a typical GEP solution observed at h = 4, 10001 generations and, \(P_{opt}=0.1\)
has dropped below 0.5 for _Fit_\(<\) 5 and below 1.0 for _Fit_\(<\) 6. Similar behavior is observed for high-frequency error where the error drops below 1 for _Fit_\(<\) 6. However, it is evident that the low frequency \(\tilde{\omega}^{2}\) errors are quite unpredictable, although the errors are \(\leq\) 1.0 for all the cases. This unpredictable behavior is attributed to the high noise at low frequencies, a feature of input data that is observable in Fig. 7 (a). As pointed out by the Dominique _et al_[39], _Fit_ by definition is a cumulative error criterion. Hence, a lower value of _Fit_ for a solution does not guarantee the local closeness to the parent function (Goody's model) in small individual subdomains.

The value of _Fit_ with the noisy data reported in Fig. 8 is \(\approx\)\(\geq\) 4. Subsequent studies were hence carried out with clean data (see Fig. 7 (b)) to determine if the GEP predictions improved. Figure 9 (a), shows a scatter plot with varying radius and color to demonstrate the effect of hyperparameter setup (head length, local optimization probability, generations) on the model performance (\(log(Fit)\)). Here, the color and radius of the circle represent the \(log(Fit)\) and the number of generations evolved in each trial respectively. Following are the key inferences that can be drawn from the figure:

* As expected, improved values of _Fit_ (below 1.0) can be achieved with clean data when compared to noisy data.
* _Fit_ drops with increasing head length, although with larger head lengths there is a risk of making the model too complex. Note that the _Fit_ values were still higher with increasing head lengths (up to 10) with 1001 generations (not shown in the plot). However, with further optimization, _Fit_ is observed to drop below 1 for h = 10 and 3001 generations (refer to marker A in Figure 9 (a)). This is because a higher number of generations allow for more mutations, thereby improving the _Fit_. Decreasing the head length also demands a larger number of generations which translates to more time spent in model evolution. Local optimization can circumvent this issue by minimizing the divergence within the population generated by the added power feature using the local least squares optimizer discussed in Section III.4.
* Increasing the local optimization probability improves the _Fit_ even at lower head lengths with less number of generations. Nevertheless, this considerably increases the time per generation as illustrated in Fig. 9 (b). Figure 9 (a) shows that increasing the optimization probability from 0.005 to 0.1 drops _Fit_ below 1 for a head length of 4 trained over 10001 generations (refer to marker B in Figure 9 (a)). The fact that a solution at the complexity of h = 4 produced comparable predictions to the ones observed for h = 10 shows the potential of the local optimizer. Further increasing the optimization probability to 1, dropped the required number of generations to 1001 resulting in a _Fit_ below 1, (refer to marker C in Figure 9 (a)).
* Higher number of generations results in a lower _Fit_ even at lower head lengths and optimization probabilities (refer to marker B in Figure 9 (a)).

All these insights lead to the conclusion that a user should prioritize mutations by evolving GEP over a greater number of generations. The subsequent diverging behavior due to mutations, if any, can be controlled with increased use of the local optimization loop. If both these solutions are ineffective, the user can consider increasing the head length as the GEP expressions with a smaller head length may not be complex enough to capture the features of the data. Figure 9 (c), shows the encouraging trend that the proposed setup that reduces the _Fit_ error also reduces the low frequency \(\tilde{\omega}^{2}\) error which was found to be unpredictable earlier. Figure 7 (b) illustrates the typical curve fit obtained with GEP (with \(h=4\), \(P_{opt}=0.05\) and 10000 generations) where the _Fit_ drops below 1.

Due to its unsupervised evolutionary nature, GEP can produce ill-performing solutions and there is always a possibility

Figure 9: (a) _Fit_ (log(_Fit_) as colorbar) as a function of head length (x-axis), local optimization probability (y-axis), and number of generations (area of circle). (b) Effect of optimization probability on the computational time required to run the algorithm through 200 generations for h = 4 (c) Predictable drop observed in low-frequency error at low _Fit_ values
of the algorithm not converging to the desired accuracy. The aforementioned hyperparameter tuning strategies aim to minimize such issues with convergence. As will be demonstrated in Section V, unlike ANN, data of reasonable quality (with minimal noise and outliers) should also be provided to ensure convergence of GEP.

### Weights and Objective functions

It is established in Section III that the database used in this study derives from different experiments with their own set of flow conditions. Table 2 consolidates different experiments and the corresponding number of datasets, \(N\). It is worth noting that the datasets collected within the same experiment can be mutually similar and differ greatly from those of other experiments which have different flow conditions. If every dataset is given equal weight/importance in the overall objective function, the model will be biased towards the over-represented data (for example, the data from Christophe). The rest of the minority datasets have little effect on the direction of the solution and suffer from poor model predictions even if the datasets are just as important [49]. One way to mitigate this problem is through weighting datasets to balance the objective function. The selection of weights is a user decision and an optimization problem in itself which is beyond the scope of this work. To illustrate the effect of a weighted objective function on the ANN and GEP predictions, a rather simple weighting strategy is considered where the weight is estimated as the ratio of the number of datasets in the least represented experiments (the data from Salze) to the experiment considered i.e. \(W=N_{Salze}/N_{Experiment}\). The weights assigned this way are listed in Table 2.

Both GEP and ANN were trained on the same input data with weighted and unweighted objective functions. Technically training can be done indefinitely and the loss calculated by the objective function may continue to drop till it reaches a minimum. However, after a steep drop of the \(lMSE\), Fig. 10 shows a marginal improvement in an ANN model performance with a further increase in epochs. For practical purposes, an early stopping criterion is therefore used to halt the training process. This can be done in a number of ways, say by setting an upper bound on epochs or by setting a lower bound on the loss calculated by the objective function. Here for this test, 20% of the randomly sampled training data is set aside as validation data. The calculated loss against validation data was then monitored and training was stopped if there is no improvement for at least 20 epochs, followed by restoring ANN weights (not to be confused with database weights) to the ones observed for the epoch with the least validation loss. Figure 10 also shows that the rate of accuracy gain (indicated by the rate of drop of error reported as weighted \(lMSE\)) is higher for the ANN models trained with \(lMSE\) as the objective function than the ones trained with \(Fit\).

The effect of different objective functions (\(lMSE\) and \(Fit\)) on the predictions of wall pressure spectra (\(\Phi_{pp}\)) using ANN is illustrated in Fig. 11 (a-d). It is observed that the predictions of ANN models trained with weighted/unweighted \(lMSE\) as the objective function are more accurate than the models trained with weighted/unweighted \(Fit\). Hence, we use \(lMSE\) as the objective function for the rest of the study unless otherwise mentioned. Since GEP predicts a different mathematical expression every run, 10 trials with 100 individuals each, having a head length = 4, and optimization probability 0.05 were evolved over 5000 generations with weighted/unweighted \(lMSE\) as the objective function. Figures 11 (e) and (f) plot the Predicted vs True values of PSD (\(\Phi_{pp}\)) obtained from the corresponding best-performing expressions which are listed below:

Weighted:

\[\widetilde{\Phi}_{pp}=\frac{0.46\widetilde{\omega}^{0.99}(\beta+1)^{1.5}}{ \frac{R_{T}^{0.55}\widetilde{\omega}^{1.02}}{\Pi^{2.76}+\Delta}+M\Delta^{6.75} \left(\frac{R_{T}+\widetilde{\omega}^{5.07}}{R_{T}^{5.84}}\right)} \tag{12}\]

and Unweighted:

\[\widetilde{\Phi}_{pp}=\frac{\widetilde{\omega}\left(\Pi^{6.71}+(\beta+1)^{1.4 }\right)}{\widetilde{\omega}^{0.87}\Delta^{0.98}M^{0.44}(1+M^{0.97})+\frac{ \widetilde{\omega}^{5.31}}{R_{T}^{7.07}}(\Pi^{5.31}+\widetilde{\omega} \Delta^{6.09})} \tag{13}\]

Interestingly, Fig. 11 shows that ANN and GEP models trained with unweighted objective functions resulted in

\begin{table}
\begin{tabular}{l l l} Experiment & N(Datasets) & Weight \\ \hline Salze & 10 & 1.00 \\ \hline Deuse & 13 & 0.77 \\ \hline Hao & 16 & 0.62 \\ \hline Christophe & 78 & 0.13 \\ \hline Total & 117 \\ \end{tabular}
\end{table}
Table 2: Weights assigned to different experiments

Figure 10: Example of the evolution of model performance in ANN training with (a) weighted \(lMSE\) and (b) weighted \(Fit\) as objective functions respectively.
Figure 11: Predictions of ANN model trained with weighted (a) and unweighted (b) _IMSE_ as the objective function; weighted (c) and unweighted (d) _Fit_ as the objective function. Predictions of GEP models, Eq. 12 and Eq. 13, trained with weighted, (e), and unweighted _IMSE_, (f), respectively as the objective functions. (Every tenth point is reported for plotting purposes)
Figure 12: Salze dataset predictions suffering (a) and (c); Hao dataset predictions improving (b) and (d), from unweighted ANN training. Salze dataset predictions suffering (e), and Christophe dataset predictions improving (f) from unweighted GEP training with \(IMSE\) as the objective function
a lower _lMSE_ error than those with the weighted ones. Although this is an unexpected trend, Table 3 delves into the performance of the models with weighted objective functions on the minority datasets. As intended, the prediction of datasets in the minority experiment (Salze) has indeed improved with weighted objective functions. In contrast, the accuracy of the datasets from a majority experiment (Christophe) has suffered due to the choice of weights. Figure 12 plots the wall pressure spectra of the majority and minority datasets to further demonstrate the effect of weighted/unweighted objective functions across different training schemes. The results show that the weighted objective functions significantly improve predictions of minority datasets (Salze). This trend is particularly observed for poorly converging training methods such as ANN with _Fit_ or the GEP. While the prediction accuracy of the dominant datasets (Christophe and Hao) decreases with the weighted training schemes, the spectra are not far off. On the other hand, unweighted schemes are observed to overfit the dominant datasets retaining their features and deteriorating the accuracy of the minority dataset predictions (refer to the predictions of inertial subrange in Fig 12(c) and (d)). These effects, however, become less pronounced in the case of an aggressively converging training scheme such as ANN with _lMSE_ where both the weighted/unweighted datasets result in competitive predictions.

The aforementioned discussion highlights the importance of optimizing the weights assigned to skewed datasets which will (a) minimize any bias in the predictions of ML models and (b) improve the performance of poorly converging models as it allows the training algorithm to focus more on minority datasets.

### ANN vs GEP

The present section compares the computational efficacy of GEP and ANN in predicting the wall pressure spectra. Both the algorithms are trained against the same input data and with the same objective function (weighted _lMSE_, Eq. 9). The trials are conducted on a machine equipped with a dual node Intel(R) Xeon(R) Gold 5120 CPU @ 2.20GHz, with 14 cores, and an Nvidia Tesla V100 32GB GPU.

Figure 13 (a) compares the convergence of ANN and GEP models with time. It is evident that the ANN achieves superior results in a shorter duration. This is in contrast to GEP which struggles to converge to the _lMSE_ achieved by ANN even after a significantly longer training. With GEP, the best-fitting model from a pool may or may not improve over the generations thereby resulting in a'stepwise' drop in _lMSE_. The mechanism for evolving a best-fitting model in GEP is mutation and the other genetic operators, but they do not necessarily result in a better model and are therefore considered diverging in nature. As a result, the next best-fitting model can appear anytime through evolution. While there are ways to control the divergence induced by mutation, such as balancing it out with aggressive selection or using a local optimization function, in general, the'stepped' convergence still follows.

Figure 13 (b) compares the virtual memory consumption with the present implementation of both algorithms. Space consumption of GEP is dependent on the size of the population, the length of each chromosome (which are basically fixed-length strings) and, additional storage required for performing actions such as evaluation, selection and, mutation. All these operations are inexpensive when compared to the ANN training. Although the memory occupied by the trainable parameters in ANN is minimal, the training process itself involves supervised learning that requires storage of gradients, and other intermediate computations that are memory intensive. Parallel computing using GPUs does help out in reducing this requirement by allowing for larger batch sizes that reduce the required number of forward and backward passes.

Figures 14 (a) and (b) compare the training time required with a reduced data resolution. As expected, both algorithms require less training time when trained against the same input data that has a lower resolution. With ANN, five trials were conducted to account for variation in the reported training time and one can observe a monotonic drop in training time with coarser input data. This is much easier to establish in the case of ANN models since they all converge aggressively resulting in competitive fits (see Fig. 14 (c)). Performing a similar analysis with GEP involves some practical complications. Different trials of the GEP algorithm require different runtimes to arrive at models with similar performance. It highlights the need to stop GEP training after a predefined number of generations although the models obtained this way can have different performances. Figure 14 (b) reports the _lMSE_ values from the 10 trials for each resolution of the sampled data. The plot shows that cases with coarser data input, in general, require less time for training and can still result in competitive predictions. To summarize, it seems beneficial to train models (specifically GEP) with coarsely sampled data, although ensuring that the local features of the data are not

\begin{table}
\begin{tabular}{l l|l l l|l l l} \hline \hline \multirow{2}{*}{Experiment} & N & \multirow{2}{*}{Weight} & \multicolumn{2}{l|}{ANN (_lMSE_ error)} & \multicolumn{2}{l}{GEP (_lMSE_ error)} & \multicolumn{1}{l}{\multirow{2}{*}{\%}} \\ \cline{4-4} \cline{6-7}  & & _lMSE_ & _lMSE_ & & _lMSE_ & _lMSE_ & _lMSE_ & _lMSE_ \\  & & & (weighted) & (unweighted) & & w.r.t & (weighted) & (unweighted) & \\  & & & training & & weighted & training & training & w.r.t \\ \hline Salze & 10 & 1.00 & 0.24 & 0.33 & +37.50 & 6.98 & 9.43 & +35.10 \\ \hline Christophe & 78 & 0.13 & 0.39 & 0.23 & -41.03 & 7.41 & 3.81 & -48.58 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Comparison of _lMSE_s in predictions across experiments
lost in the re-sampling process.

## V Guided search in gene expression programming

GEP and ANN algorithms update the model parameters through every epoch to accurately fit the data. The superior performance of ANN (in terms of speed and reliability) is attributed to training algorithms like stochastic gradient descent that guide the direction in which the ANN model parameters should be updated. On the other hand, GEP updates the models through genetic operators which rely on the evolutionary approach of trial and error. Although this approach allows for a diverse search space, it is computationally expensive. Oftentimes, with some of the known dependencies, the search space can be reduced resulting in faster convergence. This guided search approach also facilitates the discovery of complex trends in the rest of the search space. The following subsections demonstrate several guided search strategies that have been explored to accelerate GEP training and reliability.

### GEP Training Schemes

In contrast to the baseline scheme that uses the raw data, the current section introduces three additional schemes: Noise reduction with ANN filter, Omega2, and Gene4, to accelerate the GEP training.

#### Raw data scheme

This is a baseline GEP training scheme that uses the default noisy raw data and the linking function, Eq. 15, to train the GEP population.

Figure 14: Lower training times achieved with low-resolution training data - (a) ANN, (b) GEP. (c) Effect of low-resolution training data on ANN predictions for datasets with complex trends

Figure 13: (a) Model performance evolution through time, and (b) the RAM and VRAM used by the ANN and GEP models while training
#### Noise reduction with ANN filter

Noisy input datasets can misguide GEP towards sub-optimal solutions where the resulting model fails to recognize the underlying trend. The absence of feedback in mutations increases the training time required for the GEP model to converge to a global optimum. Smoothing techniques, such as kernel smoothing or moving averages, can be used to reduce the noise although these might suffer from issues like uneven noise filtering and might filter important uncertainties in parts of the data. On the other hand, Artificial neural networks are well known for their ability to handle the trade-offs between retaining information and filtering noise. Sections IV.2 and IV.3 have also demonstrated their ability to produce well-performing models of WPS with reasonable training times. Hence, we explore the strategy of ANN-filtered input data to assist the GEP algorithm.

#### Omega2 Scheme

The canonical shape of the WPS remains more or less the same across TBLs and hence a good guess of the frequency dependence of PSD in the respective regimes can be made. For example, one of the popular WPS models by Goody has three parts A, B, and C:

\[\widetilde{\Phi}_{pp}=\frac{C_{2}\widetilde{\omega}^{2}}{\left(\widetilde{ \omega}^{0.75}+C_{1}\right)^{3.7}+\left(C_{3}\widetilde{\omega}\right)^{7}}= \frac{A}{B+C} \tag{14}\]

Where,

\[C_{1}=0.5;\,C_{2}=3.0;\,C_{3}=1.1R_{T}^{-0.57}\]

The numerator \(A\) dominates at low frequencies (\(\widetilde{\omega}\to 0\)), resulting in the \(\widetilde{\omega}^{2}\) trend, which is particularly valid under zero pressure gradients (ZPG) [27]. The exponent of \(\widetilde{\omega}\), however, changes with pressure gradients as observed in Figures 2 and 3. For example, WPS of a TBL subjected to APG exhibits a \(\approx\widetilde{\omega}^{0.5}\) rise at low frequencies (see Fig. 3 from experiments of Salze _et al._[14]). While both \(B\) and \(C\) in Eq. 14 dictate the mid-frequencies, the term \(C\) governs trends at higher frequencies. Goody's model exhibits a \(\widetilde{\omega}^{-5}\) trend at higher frequencies which is consistent with those observed in Fig. 3. Unlike the other contemporary models in literature [35], Goody's model is stiff and reasonably retains its shape over a range of flow conditions.

As discussed above, the shortcomings of Goody's model which only depends on \(R_{T}\) (in addition to \(\widetilde{\omega}\)) are apparent under pressure gradients and non-canonical flows. GEP algorithm is effective in seeking additional contributing variables and accelerates the development of semi-empirical models like Goody. To assist GEP search with the inherent nonlinearity of the WPS problem, Dominique _et al._[39] prescribed the following trigenic linking function, Eq. 15 which is analogous to Goody's model:

\[Y^{GEP}=\frac{Sub-ET_{1}}{Sub-ET_{2}+Sub-ET_{3}} \tag{15}\]

Where \(Sub-ET_{i}\) represents a sub-expression interpreted from a gene in the chromosome and \(Y^{GEP}\) represents the output model.

GEP still has a tough task of working out the correct combination, especially the strong dependence on frequency (\(\widetilde{\omega}\)) across different regimes. Since the WPS at low frequencies follows a \(\widetilde{\omega}^{2}\) trend for ZPG flows, we explore the benefit of premultiplying Eq. 15 with \(\widetilde{\omega}^{2}\) to build that dependency right into the solution. The scheme further uses assistance from ANN-filtered data (modified inputs) for improved predictions. Hereafter, this strategy will be referred to as the _Omega2 scheme_ which is implemented by modifying the objective function as follows:

\[\mathit{IMSE}=\frac{1}{N}\sum_{i=1}^{N}W_{i}\left(10log_{10}\left(\widetilde{ \omega}_{i}^{2}Y_{i}^{GEP}\right)-10log_{10}Y_{i}^{true}\right)^{2} \tag{16}\]

The proposed modification facilitates the algorithm to focus on finding a better combination \(Sub-ET_{1}\) using any of the variables (\(\widetilde{\omega},\Delta,H,M,\Pi,C_{f},R_{T},\beta\) ). Table 4 lists some of the sample models obtained using the Omega2 scheme. Interestingly, for the datasets where the low-frequency trend deviates from the \(\widetilde{\omega}^{2}\) dependence, GEP is observed to predict a \(Sub-ET_{1}\) with \(\widetilde{\omega}^{i}\) (i being some index) to accurately fit the low frequencies.

#### Gene4 Scheme

Extending the omega2 approach, the _Gene4 scheme_ is aimed at further incorporating the trends at high-frequencies into the GEP models. This is achieved by modifying the original linking function in Eq. 15 with (a) an \(\widetilde{\omega}^{2}\) term in the numerator (to predict lower frequencies), (b) a \(\widetilde{\omega}^{7}\) in the denominator (to predict higher frequencies) and (c) a timescale ratio dependency \(R_{T}^{4}\) from Goody's model, a dependency also observed in other contemporary empirical models, to further assist the search. The gene4 linking function is realized using four genes as formulated below:

\[Y^{GEP}_{gene4}=\frac{\widetilde{\omega}^{2}Sub-ET_{1}+Sub-ET_{2}}{\frac{ \widetilde{\omega}^{7}}{R_{T}^{4}}Sub-ET_{3}+Sub-ET_{4}} \tag{17}\]

We have used two \(Sub-ET\)s in the numerator with only \(Sub-ET_{1}\) being pre-multiplied with \(\widetilde{\omega}^{2}\). \(Sub-ET_{2}\) allows for searching the terms that are decoupled from the \(\widetilde{\omega}^{2}\) dependence and facilitates a faster convergence than the omega2 scheme (as the algorithm does not have to look for a complex division to eliminate the effect of pre-multiplication). Table 4 lists some of the competitive models predicted using the gene4 scheme. It is interesting to note that the \(Sub-ET_{2}\) term in models 1 and 2 is independent of \(\widetilde{\omega}^{2}\) while that in model 3 is closer to \(\widetilde{\omega}^{2}\).
### Comparisons

In this section, we compare the performance of the aforementioned schemes in terms of (a) convergence rate during training and (b) model accuracy at the end of GEP evolution. For each of the four schemes, 10 trials were conducted on a pool of 100 individuals, with a head length (\(h\)) of 4 trained over 20000 generations. The input data is resampled to 25% of the original resolution and for every iteration, the population passed through the least squares optimization loop with an optimization probability of \(P_{opt}=0.05\). We probe into the temporal evolution of the training statistics with weighted _IMSE_ as the objective function.

For different training schemes, Fig. 15 (a-e) compares the temporal evolution of the resampled training statistics (see Appendix A) from the 10 trials using _IMSE_ as the objective

\begin{table}
\begin{tabular}{l l l l l} Model & \(\dfrac{\widehat{\omega}^{2}Sub-ET_{I}+Sub-ET_{2}}{\dfrac{\widehat{\omega}^{7} }{R_{T}^{T}}Sub-ET_{3}+Sub-ET_{4}}\) & _IMSE_ & \(\widehat{\omega}\to 0\) & \(\widehat{\omega}\rightarrow\infty\) \\ model1 & \(\dfrac{\widehat{\omega}^{2}\left(\dfrac{\beta^{2.45}\Pi^{5.13}}{R_{T}^{1.22}} +\dfrac{0.015}{M}\right)+\left(\dfrac{C_{f}}{M^{2.08}}+\beta^{2.84}\right) \left(\dfrac{R_{T}}{M^{0.14}}\right)}{\dfrac{\widehat{\omega}^{7}\left( \dfrac{\Delta^{2.89}H^{7.62}}{R_{T}\omega^{1.12}}\right)+\left(R_{T}^{1.29}+ \Delta^{2.03}\right)\left(M^{0.34}+\beta\right)}}\) & 10.41 & 2 & -3.88 \\ model2 & \(\dfrac{\widehat{\omega}^{2}\left(\dfrac{\Pi^{0.74}+\widehat{\omega}^{0.65}}{ \Delta^{5.02}}\right)\left(0.1R_{T}^{1.27}\right)+\left(M^{-1}\Delta^{-1}+ \beta^{1.21}\right)}{\dfrac{\widehat{\omega}^{7}}{R_{T}^{1}}14.2M^{2.96}\left( M^{1.24}+1\right)+\widehat{\omega}+5.2\widehat{\omega}^{-0.35}}\) & 10.95 & 3 & -4.35 \\ model3 & \(\dfrac{\widehat{\omega}^{2}\left(\dfrac{\beta^{2.43}\left(\Delta+\Pi\right)}{ H\Delta^{1.15}}\right)+\dfrac{\omega^{2.64}}{M\left(R_{T}+\Pi^{0.06}\right)}}{ \dfrac{\widehat{\omega}^{7}\left(\dfrac{M\Delta^{6.79}}{R_{T}^{2.39}}\right)+ \widehat{\omega}^{1.55}\left(H+\widehat{\omega}^{0.94}\right)}}\) & 12.23 & 2.64 & -4.36 \\ \end{tabular}
\end{table}
Table 5: Sample WPS models generated using gene4 scheme

\begin{table}
\begin{tabular}{l l l l l} Model & \(\widetilde{\omega}^{2}\cdot Y^{GEP}\) & _IMSE_ & \(\widetilde{\omega}\to 0\) & \(\widetilde{\omega}\rightarrow\infty\) \\ model1 & \(\widetilde{\omega}^{2}\cdot\dfrac{\widetilde{\omega}^{-1.86}\beta^{1.34}M^{-0. 55}\Delta^{-0.78}\left(R_{T}+1\right)}{\widetilde{\omega}^{1.77}+R_{T}H^{1.86}+ \widetilde{\omega}^{6.60}\left(\dfrac{\Delta}{R_{T}}\right)^{6.18}}\) & 10.61 & 0.14 & -6.46 \\ model2 & \(\widetilde{\omega}^{2}\cdot\dfrac{\widetilde{\omega}^{-2}\beta^{1.42}\left( \Pi^{3.57}\widetilde{\omega}^{0.58}+R_{T}^{1.02}\right)}{\Delta\left( \widetilde{\omega}^{5.53}\dfrac{\Delta^{5.59}}{R_{T}^{6}}+R_{T}M^{0.37}\right)}\) & 11.17 & 0.58 & -4.95 \\ model3 & \(\widetilde{\omega}^{2}\cdot\dfrac{\widetilde{\omega}^{-1.7}\beta^{1.72}M^{-0.72}}{ \widetilde{\omega}^{1.99}+\Delta\left(\dfrac{\widetilde{\omega}^{4.98}H^{8.81} }{R_{T}^{4.28}+\beta^{5.63}}+1\right)}\) & 13.67 & 0.3 & -4.68 \\ \end{tabular}
\end{table}
Table 4: Sample WPS models generated using omega2 scheme
Figure 15: Comparison of fitness distribution of the best individuals across 10 trials with \(IMSE\) as the objective function for different schemes. Evolution of following statistics: (a) \(<y>_{E}\), (b) \(median(y)\), (c) \(y_{25-75}\) (d) \(y_{max}\) (e) \(y_{min}\), of GEP trials are plotted against time. (f) Box plot comparing unweighted \(IMSE\) for the last generation, showing the scatter in \(IMSE\) from minimum to maximum across the trials. The blue box represents the trials with \(IMSE\) between the \(25-75\) percentile and the intermediate line shows the median value.)
function. These metrics include \(<y>_{E}\), \(median(y)\), \(y_{25-75}\), \(y_{max}\), and \(y_{min}\) which are briefly described in Table 6. For every trial, the best-performing model at the end of the training is recorded. Figure 15 (f) compares the _IMSE_ distributions of these best models across different schemes. Table 7 provides a quantitative overview of the relative performance of these schemes with respect to the baseline scheme, probed every 5 hours into the training period. The following inferences can be drawn from the statistics presented in Fig. 15 and Table 7:

* The gene4 scheme predicts more accurate models than other schemes throughout the evolution. When compared to the baseline scheme, Table 7 shows the superior performance of gene4 with a consistent drop of \(\approx\) 25% in its ensemble averaged _IMSE_ and \(\approx\) 20% drop in its median _IMSE_. A \(\geq\approx\) 60% drop of _IMSE_ within the interquartile range, \(y_{25-75}\) with the gene4 scheme suggests that the scheme predicts more competitive models at any instance throughout the evolution. Compared to the baseline scheme, the \(y_{max}\)_IMSE_ (worst-performing model) and \(y_{min}\)_IMSE_ (best-performing model) from the gene4 scheme trials exhibit a \(\geq\) 30% and \(\geq\) 15% lower value respectively. Trends in figure 15(a-e) also illustrate that the aforementioned improvements from the gene4 scheme are consistently observed throughout the evolution which implies that the gene4 scheme converges faster than the rest. Since the evolution of different training schemes can stop at different times, Figure 15(f) compares the _IMSE_ statistics of the last generation. Across trials, gene4 scheme has a best-performing model with \(\approx\) 30% drop in the \(y_{min}\)_IMSE_, a \(\approx\) 17% drop in the median _IMSE_ and a \(\approx\) 40% reduction in _IMSE_ within the interquartile range, \(y_{25-75}\) implying that gene4 scheme predicts more accurate models with improved reliability.
* Statistical improvements in _IMSE_ obtained using the ann-filter scheme are similar to those of the gene4 scheme albeit to a lesser extent.
* In contrast, the omega2 scheme yields either marginal or inconsistent gains over the raw data scheme throughout the training resulting in arguably worse-performing last-generation models (see Figure 15(f)). This indicates that multiplying the entire numerator (\(Sub-ET_{1}\) of Eq. 15) with \(\widetilde{\omega}^{2}\) has resulted in poorly optimized expressions. As seen in the case of model 2 in Table. 4, the algorithm has probably worked more towards decoupling the frequency terms.

Considering all the inferences, it can be concluded that the strategies involved in the development of the gene4 scheme (which include using ann-filtered input data, built-in low and high-frequency trends, and a modified linking function) are justified. Hence, the Gene4 scheme serves as a faster, and more accurate alternative to the baseline scheme while inferring models out of WPS data using the GEP algorithm.

### Validation

In this section, we evaluate the performance of the gene4 scheme on unseen data to ensure that the scheme predicts the mathematical expressions of appropriate complexity. Four datasets are discarded from each experiment from the complete database. Gene4 scheme was trained on the rest of the data with the same hyperparameter setup used in the previous section. One of the models within the interquartile range, \(y_{25-75}\), is presented below:

\[\widetilde{\Phi}pp=\frac{\widetilde{\omega}^{2}\left(\frac{\Pi^{8.05}}{\Delta }+0.22\beta^{1.54}\right)+0.69\beta^{2.66}+\Pi^{8.22}}{\frac{ \widetilde{\omega}^{7}}{R_{T}^{4}}(M^{4.07}(100(10+\Delta)))+(M^{0.44}(H^{5.81 }+\widetilde{\omega}^{2.37}))} \tag{18}\]

Equation 18 shows that the model retains the built-in frequency trends and incorporates additional complexities of the pressure gradient (\(\beta\)) and Mach number dependence. Figure 16 plots the corresponding predictions of the unseen datasets. It is evident from Figure 16(a-c) that the model presented surpasses the predictions of both Goody's model and Dominique's GEP model (which has been trained over the entire database). This ensures that the model does not overfit the data, which is in fact an encouraging step towards the development of a generalized WPS model that can reasonably predict unseen data. However, Fig. 16 (d) illustrates that, although the model performance is superior to that of Goody, there is scope to improve it further. For example, the appearance of coefficients such as 100 and 10 in this sample model shows that at times the GEP algorithm may struggle to find correct coefficients from the RNC array (Section III.4) resulting in a sub-optimal mathematical solution. However, the

\begin{table}
\begin{tabular}{l|r r r|r r r|r r r|r r r|r r r} \multirow{2}{*}{Scheme} & \multicolumn{10}{c}{Percentage _IMSE_ change w.r.t raw scheme} \\ \cline{2-13}  & \multicolumn{4}{c|}{\(<y>_{E}\)} & \multicolumn{4}{c|}{\(median(y)\)} & \multicolumn{4}{c|}{\(y_{25-75}\)} & \multicolumn{4}{c|}{\(y_{max}\)} & \multicolumn{4}{c}{\(y_{min}\)} \\ \cline{2-13}  & 5 hrs & 10 hrs & 15 hrs & 5 hrs & 10 hrs & 15 hrs & 5 hrs & 10 hrs & 15 hrs & 5 hrs & 10 hrs & 15 hrs & 5 hrs & 10 hrs & 15 hrs \\ \hline ann-filter & -10.4 & -14.4 & -15.8 & -18.5 & -19.0 & -13.8 & -57.6 & -49.0 & -58.6 & -2.7 & -18.8 & -24.4 & 5.4 & 4.9 & 4.6 \\ omega2 & 2.7 & -5.3 & -4.1 & -5.9 & -12.3 & -5.9 & -57.7 & -37.7 & -54.0 & 18.7 & -1.8 & -1.8 & 29.0 & -5.9 & -5.9 \\ gene4 & -27.6 & -27.0 & -26.1 & -29.5 & -25.3 & -19.2 & -76.0 & -52.2 & -59.4 & -33.9 & -36.1 & -39.4 & -14.9 & -19.7 & -24.4 \\ \end{tabular}
\end{table}
Table 7: Quantitative overview of the evolution of the relative performance of the proposed GEP training schemes with respect to the baseline raw data scheme. This relative performance is evaluated every five hours into the training phase.
same gene/\(Sub-ET_{3}\) does not hold a frequency term implying that the built-in trend by gene4 is valid, and with more mutations, GEP can find a superior alternative.

### Stepped schemes

We have further explored the strategy of stepped training schemes which is analogous to the idea of learning a language. One should learn the basic features of the language first, to understand the intricate constructions in the literature. Likewise learning the basic trends of the datasets that show a close resemblance to the canonical WPS is a priority. Once the GEP population learns to fit these priority datasets satisfactorily it is allowed to further improve the formulations by exposing it to a dataset that has similar canonical features with added complexities. In the present work, we train the GEP algorithm in two steps. In step 1, we prioritize exposing the population to Salze datasets which have the highest similarity to the canonical form of WPS. After the first step, the entire population, which also includes the best individual, was saved. In step 2,

Figure 16: Performance a sample GEP model trained with gene4 scheme on unseen datasets.

Figure 17: Box plot comparing unweighted \(IMSE\) for the last generation (\(40000^{th}\)), showing the scatter in \(IMSE\) from minimum to maximum across the trials. The blue box represents the trials with \(IMSE\) between the 25 - 75 percentile and the intermediate line shows the median value.
Figure 18: Dataset-wise predictions of a GEP model, Eq. 19, trained with the stepped training scheme and comparison with ANN predictions along with the predictions of Dominique’s GEP model and Goody’s model
the entire population saved in step 1 is exposed to the complete dataset for further training. For each step, the GEP algorithm was trained with the same hyperparameter environment detailed in Section V.2. Figure 17 shows the stepped gene4 scheme exhibiting superior accuracy and reliability over its non-stepped counterpart with \(\approx 10\%\) reduction in the median (from 12.2 to 10.9) and \(\approx 60\%\) drop in the \(lMSE\) spread of the interquartile range (from 3.0 to 1.1). Hence, it can be deduced that at least half of the models trained using the stepped gene4 scheme are competitive with each other. A sample model, with \(lMSE=10.61\) lying within the interquartile range, reads as follows:

\[\widetilde{\Phi}_{pp}=\frac{\widetilde{\omega}^{3.12}+\beta^{2.14}\widetilde{ \omega}^{2.41}+\left(\frac{C_{f}}{M^{4.42}\Delta^{2.34}}\right)\widetilde{ \omega}^{2}+\beta^{2.14}R_{T}^{0.29}}{\left(M^{3.2}\Delta^{6.35}+R_{T}\right) \frac{\widetilde{\omega}^{7.71}}{R_{T}^{5}}+\Delta^{0.9}(\widetilde{\omega}^{ 2.97}+H^{2})} \tag{19}\]

where the low and the high frequency exponent predictions are,

\[\widetilde{\omega}\to 0=\widetilde{\omega}^{3.12};\ \widetilde{\omega} \rightarrow\infty=\widetilde{\omega}^{-4.59}.\]

Figure 18 further compares the WPS predictions of the sample GEP model (Eq. 19) obtained using stepped gene4 scheme across different datasets. Its superior accuracy as compared to Goody's model and Dominique's GEP model is apparent.

Figure 19: Predictions of a GEP model, Eq. 19, trained with the stepped training scheme and comparison with ANN predictions along with the predictions of Dominique’s GEP model and Goody’s model. (Every tenth point is reported for plotting purposes)
In particular, Figures 18 (a,b,e) illustrate the high complexity captured by this model. Interestingly, its frequency trends and magnitude are in close agreement with those predicted by the ANN model. Figure 18 (c) assures that its predictions of the canonical WPS are unaffected and competitive with other models. Despite its superior performance, Fig. 18 (b), (d), and (f), show that the model fails to capture some local features and under or over-predicts the trends in certain datasets.

The scatter plot in Figure 19 compares the predicted vs true values for the entire dataset comprising different flow conditions (APG, FPG, and ZPG). The aforementioned GEP model trained with the stepped gene4 scheme significantly improved the predictions across the entire dataset resulting in a _IMSE_ of 10.6 in contrast to _IMSE_ of \(\approx 60\) with Dominique's GEP model and \(\approx 90\) with Goody's model. In particular, the scatter in the predictions of datasets with APG and FPG has considerably reduced when compared to Dominique's GEP model and Goody's model. Despite these improvements, the _IMSE_ of the proposed GEP model is still an order of magnitude higher than that of the ANN. It implies that there is further scope to develop a generalized WPS model using GEP through improvements in training strategies.

## VI Conclusions

The study presents a machine learning-based framework that uses data-driven modeling to predict wall pressure spectra (WPS) underneath turbulent boundary layers. Different datasets of WPS from experiments and high-fidelity numerical simulations covering a wide range of pressure gradients and Reynolds numbers are considered for training. This dataset however appears to be skewed and a cosine-similarity matrix has been used to quantify the visual resemblance of WPS trends across the experiments. The efficacy of two machine learning techniques, namely artificial neural networks (ANN) and gene expression programming (GEP) is evaluated. Firstly, an optimal hyperparameter environment is identified that yields the most accurate predictions for the respective ML methods. This includes assessing the effect of objective functions (\(LMSE/Fit\)) on the convergence rate of ANN models. Of these, \(Fit\) (which is a multi-objective function comprising both _IMSE_ and mean squared error) is shown to converge at a slower rate. Interestingly, it has been observed that the prediction accuracy of such weakly converging training methods can be improved by increasing the weight of minority datasets.

For a given input database, the computational resources (training time and memory consumption) of some of the best-performing ANN and GEP models are compared. In terms of accuracy, the results show the clear superiority of ANN over GEP with the logarithmic mean squared error (_IMSE_) of ANN being less than 1 while that of GEP being around \(O(10)\). The corresponding training time of ANN models is also 8 times lower (\(\approx 3\)_hours_) than that of the GEP (\(\approx 24\)_hours_), despite a higher memory consumption. Nevertheless, the advantage of GEP lies in predicting a realizable closed-form mathematical expression. In contrast to the unconventional structure of ANN models, these expressions from GEP can provide direct physical insight.

Novel training schemes are devised to address the shortcomings of GEP. These include (a) ANN-assisted GEP to reduce the noise in the training data (b) exploiting the physical trends in the spectra observed at low and high frequencies to guide the GEP search (c) a stepped training strategy where the chromosomes are first trained on the canonical datasets followed by the datasets with complex features. The first two training strategies are shown to accelerate the convergence and yield more consistent models (\(\approx 40\%\) reduction in the spread of _IMSE_ in the interquartile range) with superior accuracy (\(\approx 17\%\) reduction in the median _IMSE_). Stepped training strategy further improves upon the reliability and accuracy with an additional reduction of \(\approx 60\%\) and \(\approx 10\%\) in the aforementioned statistics respectively. The predictions of the resulting GEP models appear to have captured the complex trends of WPS while surpassing the accuracy of both Dominique's GEP model and Goody's model.

With the inclusion of the ever-evolving data repository of WPS, the authors believe that the methods and insights from the present study will facilitate the discovery of more consistent and generalized models. The goal is to discover GEP models that result in competitive predictions to the ones produced by the ANN models. The analytical nature of GEP models will allow for the generalization of WPS beyond the trained regimes.

## Acknowledgements

The authors wish to acknowledge NVIDIA for generously awarding Quadro P6000 and A100 GPU cards under the Academic Hardware Grant Program. The authors also wish to acknowledge the VKI team: Joachim Dominique, Jan Van den Berghe, Dr. Christophe Schram, Dr. Miguel Alfonso Mendez, Dr. Julien Christophe, and Dr. Richard D. Sandberg (from the University of Melbourne) for making their source code and data publicly available. Dr. Nagabhushana Rao Vadlamani also acknowledges the financial support from Science and Engineering Research Board (SERB), India under Mathematical Research Impact Centric Support (MATRICS) scheme (MTR/2022/000807).

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability statement

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

Here, we provide details of the resampling procedure which maps generations in GEP evolution with wall-clock time. At the end of every generation, we probed the Weighted objective function values (of the best individual in the generation) and the wall-clock time spent for the generation. Although each of the 10 trials is trained for a fixed set of generations, the training can take different wall-clock times. Hence, the statistics are re-sampled over a common time interval using a Piece-wise Cubic Hermite Interpolating Polynomial (PCHIP) interpolator [50] as illustrated in Table 8. Although this shortens the observation window to the trial with the shortest runtime to reach 20000 generations, it facilitates the comparison of statistics in wall-clock time across different trials and training schemes. This approach of choosing the shortest window is justified since the interest lies in developing methods that converge faster.

## Nomenclature

\(<y>_{E}\): Ensemble average of the batch of outputs

\(\beta\): Clauser parameter

\(\Delta\): Zagarola-Smits's parameter

\(\delta\): Boundary layer thickness

\(\delta^{*}\): Displacement thickness

\(v\): Dynamic viscosity

\(\omega\): angular frequency \(\omega=2\pi f\)

\(\Omega_{z}\): Vorticity component in \(z\) (spanwise) direction

\(\overline{p}\): Mean pressure

\(\partial_{n}\): Partial derivative with respect to \(n^{th}\) coordinate

\(\Phi_{pp}\): Power spectral density

\(\Pi\): Wake strength parameter

\(\rho\): Density at wall

\(\tau\): Time

\(\tau_{w}\): Wall shear stress

\(\theta\): Momentum thickness
