Data-driven optimization of a gas turbine combustor: A Bayesian approach addressing NO_x emissions, lean extinction limits, and thermoacoustic stability

Abstract

The design of gas turbine combustors for optimal operation at different power ratings is a multifaceted engineering task, as it requires the consideration of several objectives that must be evaluated under different test conditions. We address this challenge by presenting a data-driven approach that uses multiple probabilistic surrogate models derived from Gaussian process regression to automatically select optimal combustor designs from a large parameter space, requiring only a few experimental data points. We present two strategies for surrogate model training that differ in terms of required experimental and computational efforts. Depending on the measurement time and cost for a target, one of the strategies may be preferred. We apply the methodology to train three surrogate models under operating conditions where the corresponding design objectives are critical: reduction of NO_x emissions, prevention of lean flame extinction, and mitigation of thermoacoustic oscillations. Once trained, the models can be flexibly used for different forms of a posteriori design optimization, as we demonstrate in this study.

Impact statement

The development of gas turbine combustors for safe, stable, and low-emission operation under various load conditions is a highly challenging engineering task that requires extensive testing and is usually approached through iterative trial-and-error procedures. We present a data-driven approach based on multiple probabilistic surrogate models that automatically selects optimal burner designs from a large parameter space, requiring only a few experimental data points. We show how different criteria that require separate measurements can be considered in the automated routine, and thus the design process can be rendered significantly more efficient. The proposed approach is particularly suitable for the design of gas turbine burners for operational flexibility but can also be adapted to other costly design processes.

1. Introduction

The design of gas turbine (GT) combustors is a multifaceted challenge characterized by stringent requirements for emission control, flame anchoring, and thermoacoustic stability. Employing lean premix combustion in modern GTs effectively minimizes the emissions, yet it presents challenges for flame stabilization. Typically, combustors are engineered to a full-load operating point, where aerodynamic flame stabilization is achieved within the flow field. However, one of the main strengths of gas turbines lies in their exceptional load flexibility, that is, the ability to adapt their power output to meet changing demands. The part-load operation is realized by reducing the fuel flow, which reduces the heat output of the flame. The associated changes in local stoichiometry, pressure, and the flow field result in the variation of flame speed, bearing the risk of lean blowout, or extinction (Lieuwen, 2012). Furthermore, the flame’s susceptibility to acoustic perturbations varies with the operating parameters, potentially leading to thermoacoustic instability (Krebs et al., 2002; Poinsot, 2017; Beuth et al., 2023). To enable part-load operation, most technical combustion systems use multiple injection stages to adapt to changing conditions (Janus et al., 2014; Krebs et al., 2022; Pennell et al., 2023). Although each additional fuel line enhances operational flexibility, it inevitably increases system complexity. To fully exploit the capability of operational flexibility, combustor design must take into account the complexities of varying operating conditions, which further complicates the engineering process. Addressing diverse load conditions demands extensive testing across multiple operation conditions, typically requiring individual experimental measurements or computational fluid dynamics (CFD) simulations for each operating point. The effort required to design a combustion chamber with regard to operational flexibility is therefore particularly high, while tests of different designs (experiments and CFD) are limited to a few.

In order to overcome this hurdle, a data-driven approach is pursued in this work, in which probabilistic surrogate models are trained using data from automated experiments, which can then be used for prediction and design optimization. The individual surrogate models map the parameterized burner design to target variables and, as such, are system specific. The presented methodology, however, is universally applicable and particularly suited to the multi-objective design of GT burners, as will be elucidated in the following sections.

In aerodynamics, it is common to use data-driven methods for automated CFD-based design optimization (Martins, 2022), for example, for design of blades in turbo machinery (Li and Zheng, 2017) or shape design of tall buildings (Xie, 2014). In the field of reacting flows, such optimization methods are less common. This is because the associated simulations are generally more complex, hardly automatable, and the design of a combustor involves numerous degrees of freedom and different design objectives (Rogero, 2002). Yet, in some studies, optimization methods were applied to combustor design, for example, in combination with semiempirical preliminary design tools (Despierre et al., 1997; Rogero, 2002; Pegemanyfar et al., 2006; Fuligno et al., 2009; Angersbach et al., 2013) or steady flow simulations (Jeong et al., 2006; Janiga and Thévenin, 2007; Motsamai et al., 2008; Elmi et al., 2021; Yang et al., 2023), which are numerically sufficiently inexpensive to enable iterative optimization. Optimization involving combustion dynamics poses specific challenges, as these are generally hard to predict with sufficient accuracy and at affordable cost. Only a handful of studies report application of purely data-driven optimization to unsteady numerical simulations (Wankhede et al., 2011) or to low-order thermoacoustic network models (Reumschüssel et al., 2022a). Making use of the governing physical equations, adjoint-based methods have been developed for efficient computation of thermoacoustic stability with respect to parameter changes (Magri, 2019). These methods can be embedded in classic optimization routines, for example, to maximize thermoacoustic damping rate, as demonstrated by Aguilar and Juniper (2018). Another approach to incorporate physical knowledge into optimization was recently proposed by Huhn and Magri (2022), who combined model-informed neural networks for time prediction of a chaotic thermoacoustic system with Bayesian techniques for efficient optimization based on sparse data.

While low-order models and CFD simulations generally allow flexible variation of the geometric parameters of combustion chambers, they are always subject to modeling errors. Experimental optimization, in contrast, is only realizable with automatically variable burner hardware. In addition, many variables that are essential for assessing performance can only be detected indirectly and with measurement noise. Accordingly, the literature on experimental applications of data-driven optimization is very limited. Multi-objective optimization of the fuel injection system of an industrial GT combustor has been conducted by Büche et al. (2002) and Paschereit et al. (2003). In these studies, emission of nitrogen oxides (NO _x) and flame oscillations are measured simultaneously and a Pareto front is identified by iteratively varying the injection configuration through an evolutionary algorithm. Evolutionary algorithms are a well-known machine learning technique, which are known for robust detection of global minima. This robustness comes at the price of many function calls, required until convergence (Eiben et al., 2003).

A method that is more sample-efficient is Bayesian optimization (BO), which uses probabilistic surrogate models to find optima in the design space. In BO, surrogate models are usually determined from the available data using Gaussian process regression (GPR, also known as Kriging (Krige, 1951; Rasmussen and Williams, 2006)). BO is used in various areas where function calls are time-consuming or expensive, such as optimizing policies for robot navigation (Martinez-Cantin et al., 2007) or the tuning of hyperparameters of large neural networks (Snoek et al., 2012). See Shahriari et al. (2016) for an extensive review. We adopted BO for experimental combustor optimization with respect to emissions of nitric oxides (NO _x) (Reumschüssel et al., 2022b). Our work, in which a combustor with automatically variable fuel injection was used, confirmed that the method is very well suited for experimental optimization since it converges after only a few experimental measurements. In that study, full-load GT operation was simulated in the test rig and a design that leads to minimal NO _x emissions was identified automatically. Despite the high efficiency of the method, the determined optimal design proved to be insufficient in terms of flame stabilization under different operating conditions. Under full-load operation, limiting the emission of NO _x is a critical design goal. However, in part-load operation, where the combustor operates at a lower global equivalence ratio, flame stability is usually a greater problem. To simultaneously consider NO _x emissions, flame stabilization and other potential design goals within a data-driven framework, multi-objective optimization methods could be employed. Although several viable methods using BO for multi-objective problems have been proposed, for example, by Bradford et al. (2018) and Daulton et al. (2021), their application in the context of experimental combustor design presents particular challenges. This is because evaluating different design objectives requires testing under different experimental conditions. Since it is time-consuming to change the operating mode of the test rig, evaluating each design candidate against all objectives sequentially in an automated optimization loop is impractical. Therefore, instead of seeking to identify optimal designs from iterative tests in which sampling is guided by data-based surrogate models (which is the concept of BO), we first employ an automated routine to train the surrogate models of the combustion system in different operating points, which we then combine for a posteriori multi-objective optimization. We applied this approach to consider the trade-off between NO _x emissions under full-load and part-load flame stabilization in a previous study (Reumschüssel et al., 2024b). In this study, we generalize the methodology. We present a data-driven method that allows different operating points and different design criteria to be considered within a surrogate-model based optimization, allowing a sample-efficient, multi-objective design process for safe and flexible operation.

Successively, three surrogate models are trained that map the combustor design to different objectives under different operating conditions. Namely, we train one model for NO _x emissions in conditions that mimic full-load GT operation, a second model for extinction limits under part-load operation, and a third model that predicts thermoacoustic oscillations in a third operating point. We present training strategies for the models that require a minimal number of experimental data points. The models allow cost-efficient and quantitatively precise prediction of design objectives for all burner designs within the design space. After training, they can be combined and used for multi-objective optimization. In this way, our methodology addresses the intricate task of experimental GT burner design for load flexibility in an automated, efficient, and precise manner.

The remainder of this article has the following structure. First, we introduce the experimental test rig and the variable combustion system that allows automated testing of different combustor designs. Then, the strategy for the sample-efficient training of the surrogate models is presented. The subsequent results section first describes the training of the three models before combining them to demonstrate how they can be used for multi-objective burner optimization.

2. Experimental setup

The proposed approach is applied to design the pilot injection unit of an industrial swirl combustor utilized in the MGT6000 GT. In that turbine, six of these combustors are arranged in a ring around the shaft to form the combustion system. During full-load operation, the combustors operate at approximately 14 bar with a heat input of approximately 4 MW each, which corresponds to roughly 300 kW under atmospheric test conditions. For the present study, the swirl combustor is installed in a test rig that allows to operate the combustor in its original dimensions under atmospheric pressure. The test rig is equipped with a preheated air supply and provides visual access to the flame through a quartz glass combustion chamber of inner radius $ {R}_{\mathrm{g}}=112\mathrm{mm} $ . Machine-equivalent load conditions are simulated by adapting the volumetric airflow rate, the inlet air temperature and the global equivalence ratio to the corresponding values from the gas turbine. This also implies a machine-equivalent momentum ratio between injected fuel and main flow, such that mixing properties similar to those under pressure can be expected. An overview of the setup is shown in Figure 1. The test rig is equipped with a variety of sensors for monitoring and measurements. The ones relevant for this study are mentioned hereafter. The time-averaged flame shape is captured from UV images. To detect the main reaction zone of the flame, the camera is therefore equipped with a band-pass filter, transmitting light around the wavelength range of active OH* radicals. These are known to be emitted in heat-releasing reactions within the flame and thus serve as an indication for the location of heat release (Guethe et al., 2012). Assuming rotational symmetry of the time-averaged distribution of OH* emission, the line-of sight integrated intensity captured by the camera is deconvoluted to a planar distribution across the flame $ {\mathcal{I}}_{{\mathrm{OH}}^{\ast }} $ by using an inverse Abel transform (Abel, 1826; Moeck et al., 2013). To assess heat release fluctuations $ \hat{q} $ , the integral intensity of OH* chemiluminescence $ {I}_{{\mathrm{OH}}^{\ast }} $ was recorded using a photomultiplier tube at a sampling rate of 1 kHz. Denoting fluctuations with $ \hat{\Big(}.\Big) $ and time-averaged quantities by $ \overline{(.)} $ , this reads

(1) $$ \hat{q}/\overline{q}\sim {\hat{I}}_{{\mathrm{OH}}^{\ast }}/{\overline{I}}_{{\mathrm{OH}}^{\ast }}. $$

Figure 1. Test rig for combustion tests under atmospheric pressure with variable burner system.

It should be mentioned at this point that both sides of Eq. (1) are only equal if fluctuations in equivalence ratio can be excluded (Schuermans et al., 2010). This is not the case here, which is why the fluctuations of $ {I}_{{\mathrm{OH}}^{\ast }} $ only serve as a qualitative indicator for ranking different designs against each other. To analyze the composition of combustion products, a portion of the exhaust gas is directed to a gas analysis system. Within this system, carbon monoxide levels are quantified through infrared spectroscopy, while the assessment of NO _x and NO₂ is conducted via wet analysis employing a chemiluminescence detector. Furthermore, the presence of unburned hydrocarbons (UHCs) in the exhaust gas is determined using a flame ionization detector. Finally, the gas temperature in the exhaust pipe is monitored with several thermocouples.

The swirl combustion system under investigation is shown in Figure 2. The schematic shows the combustor as it is installed in the atmospheric test rig. The airflow, indicated by black arrows in the figure is fed through radial swirl generators. In the gas turbine, the airflow is directed along the outside of the combustion chamber walls prior to entering the chamber. To mimic the flow conditions from the can-combustor assembly while providing visual access to the flame through the quartz combustion chamber the test rig features a special air distribution system. Engine-equivalent inflow conditions are generated by air preheating and circumferential distribution of the airflow around the swirlers. The main fuel line injects natural gas within the swirl vanes. As shown in a previous study, this injection stage generates a high degree of premixing which leads to a uniform stoichiometry at the flame front and contributes only moderately to NO _x formations (Reumschüssel et al., 2021). The combustion system includes a secondary fuel line, referred to as the pilot fuel line, which enhances operational flexibility. The pilot injects fuel toward the center of the burner and into the central recirculation zone near the root of the flame, resulting in local enrichment (Reumschüssel et al., 2021). For orientation, the cylindrical coordinate system shown in Figure 2 is introduced, in which $ z $ describes the axial coordinate in the main flow direction and $ r $ is the radial coordinate.

Figure 2. Schematic of the swirl combustor, including the AUTOPILOT with 61 fuel injectors. Red arrows indicate fuel flow and black arrows indicate airflow. Left: Side view of the combustion system, installed in the atmospheric test rig. Right: Sectional view from downstream with AUTOPILOT highlighted. Circular markers indicate injectors aligned perpendicular to the burner base plate and rectangular markers represent inclined injectors.

In order to carry out data-driven combustor development, the investigated burner is equipped with a special, highly variable pilot fuel injection system, which was manufactured using selective laser melting. The custom piloting unit is referred to as the Automatized Experimental Optimization Pilot (Autopilot). In the series (standard) version of the GT combustor, the pilot is fed from one plenum and controlled through a single fuel line. In contrast, the Autopilot allows variable fuel injection, featuring a total of 61 injection locations, each connected to a separate fuel line with an individual magnet valve. The fuel distribution system with individual valves and fuel lines is shown in Figure 1. By controlling the 61 valves, the injectors can be opened and closed individually enabling a large variety of possible injection schemes to be tested. The injection holes are arranged in four rings around the symmetry axis of the burner, two rings with injection holes aligned perpendicular to the burner head plate, Rings I and III, and two rings with injection holes inclined inward, Rings II and IV as shown in Figure 2. Rings I, II, III, and IV have 16, 16, 15, and 14 injection holes, respectively, evenly distributed around the circumference at different radial distances—the uneven number of injectors results from spatial constraints.

2.1. Design parameterization

To enable the automatic optimization of the burner design, it is imperative to define a set of describing variables. In the present case, the potential design configurations arise from the opening and closing of the 61 fuel lines, connected to the Autopilot. The most comprehensive representation of the input space would thus comprise 61 Boolean variables. To increase the efficiency of the optimization procedure, a preselection is made by using a four-dimensional parameterization. To this end, four independent input variables are defined, $ \mathbf{x}=\left[{x}_{\mathrm{I}},{x}_{\mathrm{I}\mathrm{I}},{x}_{\mathrm{I}\mathrm{I}\mathrm{I}},{x}_{\mathrm{I}\mathrm{V}}\right] $ , corresponding to the number of open injectors on the four rings and thus ranging from 0 and $ {N}_{\mathrm{I}}=16 $ , $ {N}_{\mathrm{II}}=16 $ , $ {N}_{\mathrm{III}}=15 $ and $ {N}_{\mathrm{IV}}=14 $ , respectively. Accordingly, the design space $ {\Omega}_{\mathbf{x}} $ is defined as follows

(2) $$ {\Omega}_{\mathbf{x}}=\left\{\mathbf{x}\in {\mathrm{\mathbb{Z}}}_{16}\times {\mathrm{\mathbb{Z}}}_{16}\times {\mathrm{\mathbb{Z}}}_{15}\times {\mathrm{\mathbb{Z}}}_{14}:\mathbf{x}=\left[{x}_{\mathrm{I}},{x}_{\mathrm{I}\mathrm{I}},{x}_{\mathrm{I}\mathrm{I}\mathrm{I}},{x}_{\mathrm{I}\mathrm{V}}\right],\hskip1em \sum \limits_{i=\mathrm{I}}^{\mathrm{I}\mathrm{V}}{x}_i\ge 5\right\}, $$

where $ {\mathrm{\mathbb{Z}}}_l $ is the set of integers between 0 and $ l $ and configurations with less than five open injectors in total are excluded from the design space to keep the pressure loss across the pilot fuel line at a reasonable level. The definition of $ {\Omega}_{\mathbf{x}} $ comprises 69,360 possible combinations in total. The numbers of open injectors on the four rings are distributed across the injector positions according to the following mapping.

(3) $$ {\mathrm{f}}_{\mathrm{p}}:\mathbf{x}\to {\mathbf{x}}_{\mathrm{I}/\mathrm{O}},\hskip2em {[{\mathbf{x}}_{\mathrm{I}/\mathrm{O}}]}_{i,j}=\{\begin{array}{ll}1,& \mathrm{i}\mathrm{f}\hskip2pt j\in \mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\hskip2pt ([1,2,\dots, {x}_i]\frac{N_i}{x_i})\\ {}0,& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\operatorname{},\hskip1em i=\mathrm{I},..,\mathrm{I}\mathrm{V} $$

where $ {\left[{\mathbf{x}}_{\mathrm{I}/\mathrm{O}}\right]}_{i,j} $ denotes whether the $ j $ th injector on Ring $ i $ is opened $ \left({\left[{\mathbf{x}}_{\mathrm{I}/\mathrm{O}}\right]}_{i,j}=1\right) $ or closed $ \left({\left[{\mathbf{x}}_{\mathrm{I}/\mathrm{O}}\right]}_{i,j}=0\right) $ with $ i $ counting from inside to outside and $ j $ counting circumferentially. This parameterization spans a design space of configurations with approximately rotationally symmetric distributions, as these are considered most relevant for technical use. Asymmetrical fuel injection schemes would lead to an inhomogeneous temperature distribution and thus potentially to an excessive thermal load on the components. This preselection significantly reduces the problem dimensionality and thus contributes fundamentally to the sample-efficient feasibility of a global optimization approach for experimental tests. Figure 3 shows the configuration of the Autopilot for four exemplary burner designs x included in $ {\Omega}_{\mathbf{x}} $ . For the indicated values of x, the fuel lines connected to the injectors indicated in red are opened.

Figure 3. Examples of burner designs that can be generated from the chosen design parameterization. For the indicated values of $ \hat{\mathrm{f}}\left(\mathbf{x}\right) $ , the injectors marked in red are open.

3. Methodology for surrogate model training

We intend to apply multi-objective optimization for GT combustor design addressing thermoacoustic oscillations $ \hat{q} $ , emissions of NO _x $ {c}_{{\mathrm{NO}}_x} $ , and the risk of flame extinction, $ {r}_{\mathrm{FE}} $ . In order to consider the various target variables simultaneously without requiring sequential tests of individual candidate designs, we first train surrogate models that allow the prediction of these quantities and upon which the design optimization can be carried out a posteriori. The goal of the surrogate model training is to find mappings $ {\hat{\mathrm{f}}}_{\mathrm{c}}\left(\mathbf{x}\right) $ which minimize the maximum deviation between actual values $ \hat{\mathrm{f}}\left(\mathbf{x}\right) $ (as obtained from measurements $ {\hat{\mathrm{f}}}_{\mathrm{c}}\left({\mathbf{x}}_{\mathrm{I}/\mathrm{O}}\right) $ of parameterized configurations $ {\mathbf{x}}_{\mathrm{I}/\mathrm{O}}={\mathrm{f}}_{\mathrm{p}}\left(\mathbf{x}\right) $ ) and predicted values $ \hat{\mathrm{f}}\left(\mathbf{x}\right) $ :

(4) $$ {\hat{\mathrm{f}}}_{\mathrm{c}}(\mathbf{x})=\underset{\hat{\mathrm{f}}}{\mathrm{argmin}}\ \underset{\mathbf{x}\in {\Omega}_{\mathbf{x}}}{\max}\mid \hat{\mathrm{f}}(\mathbf{x})-\mathrm{f}(\mathbf{x})\mid, \hskip2em \mathrm{f}(\mathrm{x})={\mathrm{f}}_{\mathrm{c}}({\mathrm{f}}_{\mathrm{p}}(\mathbf{x})). $$

The concept is summarized graphically in Figure 4. The surrogate models approximate the combination of the parameterization function and the combustion experiments. Separate surrogate models are trained to predict the three target quantities, which can subsequently be used for a multi-objective design study,

(5) $$ {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}:\mathbf{x}\to {c}_{{\mathrm{NO}}_x},\hskip2em {\hat{\mathrm{f}}}_{\mathrm{FE}}:\mathbf{x}\to {r}_{\mathrm{FE}},\hskip2em {\hat{\mathrm{f}}}_{\mathrm{TA}}:\mathbf{x}\to \hat{q}. $$

Figure 4. Interpretation of the combustion system as a function from injection configuration to measured quantities and approximation through a surrogate model.

The surrogate models are acquired by Gaussian Process Regression, which is a technique to generate probabilistic predictions based on a set of training data. The selected approach for GPR-based modeling is detailed below. To generate training data for GPR, we use two different methods to sample designs for experimental tests. In the first method, all design samples are selected upfront according to a specific sampling plan. Subsequently, GPR is performed once on the complete dataset. This approach requires minimal computational overhead, but bears the risk of excessive measurement effort as the number of training points has to be estimated in advance. For the second sampling method, we were inspired by BO. Therein two steps are repeated iteratively. First, a GPR step is performed to approximate the latent function to be optimized based on the available measurement data. This involves adopting the model parameters to the data as will be described below. The second step is the selection of the next measurement points based on the GPR prediction. For this purpose, an acquisition function is used as a sampling criterion, the choice of which is also discussed below. In contrast to the first method, the second method requires more computational effort for the individual regression steps, but allows the monitoring of a convergence metric, so that only a minimum number of training samples is tested.

3.1. Gaussian process model for combustion experiments

Gaussian Process Regression (Rasmussen and Williams, 2006) is used to generate an approximation of the functional relationship between burner design $ \mathbf{x} $ and target quantities $ \mathrm{f}:{c}_{{\mathrm{NO}}_x},{r}_{\mathrm{FE}},\hat{q} $ . As a result of the burner design, the parameter vectors describing the burner design $ \mathbf{x} $ lie on the discrete domain $ {\Omega}_{\mathbf{x}} $ . Yet, to introduce a Gaussian process model for $ \mathrm{f} $ , we first consider input vector $ \boldsymbol{x} $ on the continuous equivalent of $ {\Omega}_{\mathbf{x}} $ that is referred to as $ \Omega $ and then ensure that any function evaluations (i.e., experiments) are performed exclusively on samples from the integer subdomain $ {\Omega}_{\mathbf{x}} $ .

In GPR, Bayes’ conditional probability rule is used to predict function values based on some known (measured) input–output data. GPR begins with a prior belief about the function that describes the data. This prior is represented by a stochastic process, characterized by a mean function $ \mu \left(\boldsymbol{x}\right) $ , which determines the function shape of $ \mathrm{f} $ considered most probable and a covariance function $ k\left(\boldsymbol{x},{\boldsymbol{x}}^{\ast}\right) $ that characterizes deviations from the mean,

(6) $$ \mu \left(\boldsymbol{x}\right)=\unicode{x1D53C}\left[\mathrm{f}\left(\boldsymbol{x}\right)\right], $$

(7) $$ k\left(\boldsymbol{x},{\boldsymbol{x}}^{\ast}\right)=\unicode{x1D53C}\left[\left(\mathrm{f}\left(\boldsymbol{x}\right)-\mu \left(\boldsymbol{x}\right)\right)\left(\mathrm{f}\left({\boldsymbol{x}}^{\ast}\right)-\mu \left({\boldsymbol{x}}^{\ast}\right)\right)\right]. $$

Here, $ \unicode{x1D53C} $ is the expectation operator and $ \boldsymbol{x},{\boldsymbol{x}}^{\ast}\in \Omega $ . In general, knowledge about the function under investigation that is available a priori of tests can be incorporated into the model through the choice of an appropriate prior distribution. The prior mean function will influence the posterior function prediction, particularly in regions of the parameter space where data are scarce, as, for example, explained by Garnett (2023). For BO applications in which the mean also affects the sampling policy, a well-considered choice can therefore enhance optimization performance (De Ath et al., 2020). As described below, in the present case, the mean value is not considered in the sampling criterion, and the training data are sampled throughout the design space. Therefore, and for the sake of simplicity, a constant mean value of zero is employed,

(8) $$ \mu \left(\boldsymbol{x}\right)\equiv 0. $$

The prior covariance function determines joint probability of function realizations and thus encodes assumptions about smoothness and characteristic shape of the function. In the present case, it can be reasonably assumed that the fuel injection through the four different rings will have a similar effect on the target variables under investigation. However, due to the differing distances to the flame front and different positions in the main flow field, it can be postulated that these effects will scale differently with the individual components of $ \boldsymbol{x} $ . Furthermore, we aim to model the impact of the geometry and location of the four fuel rings, rather than that of individual injectors, a paradigm which has already been integrated into the design parameterization. Consequently, we impose a high degree of smoothness of f over the number of injectors $ {x}_{\mathrm{I}} $ to $ {x}_{\mathrm{IV}} $ and select a kernel function of squared exponential type with one lengthscale per coordinate, for which the covariance function reads

(9) $$ k\left(\boldsymbol{x},{\boldsymbol{x}}^{\ast}\right)=\exp \left(-\sum \limits_{b=1}^N{\boldsymbol{\theta}}^{(b)}{\left({\boldsymbol{x}}^{(b)}-{\boldsymbol{x}}^{\ast (b)}\right)}^2\right). $$

Here, superscript $ (b) $ indicates the bth component of a vector and $ \theta $ is a vector of hyperparameters controlling the correlation in the corresponding directions. Furthermore, GPR permits the consideration of measurement noise by defining the belief in measured function values through an observation model. Within the considered measurement range, it can be reasonably assumed that no significant systematic measurement error will be introduced by the employed measurement systems. Therefore, we simply assume function realizations $ y $ (i.e., measurements) at location $ \boldsymbol{x} $ to be corrupted by additive Gaussian noise $ \unicode{x025B} $ independent of measured value,

(10) $$ y=\mathrm{f}\left(\boldsymbol{x}\right)+\unicode{x025B}, \hskip1em \unicode{x025B} \sim N\left(0,{\sigma}_{\mathrm{n}}^2\right). $$

Here, $ z\sim N\left(m,s\right) $ denotes, that quantity $ z $ is normal distributed, with mean $ m $ and covariance $ s $ and $ {\sigma}_{\mathrm{n}}^2 $ is the noise variance. Introducing process variance $ {\sigma}_{\mathrm{f}}^2 $ , this prior induces a multivariate Gaussian distribution for any finite set of $ l $ noisy observations $ \mathbf{y}={\left[{y}_1,\dots, {y}_l\right]}^T $ at corresponding locations $ {\boldsymbol{x}}_1 $ to $ {\boldsymbol{x}}_l $ ,

(11)

Here, $ \mathbf{1} $ is the unit matrix of size $ l $ , $ \mathbf{K} $ is the noise-free covariance matrix, $ \mathbf{C} $ is the covariance matrix of the noisy targets, and superscript $ \left(i,j\right) $ denotes the matrix entry in row $ i $ and column $ j $ .

The prior can be calibrated to available measurement data by adjusting the six hyperparameters of the model, $ {\boldsymbol{\theta}}_1 $ to $ {\boldsymbol{\theta}}_4,{\sigma}_{\mathrm{f}} $ and $ {\sigma}_{\mathrm{n}} $ . This is done by maximizing the log marginal likelihood of observing the measured values $ {\mathbf{y}}_{\mathrm{m}} $ at the corresponding locations $ {\mathbf{X}}_{\mathrm{m}} $ ¹, given $ \boldsymbol{\theta}, {\sigma}_{\mathrm{x}} $ and $ {\sigma}_{\mathrm{n}} $ ,

(12) $$ \log L\left({\mathbf{y}}_{\mathrm{m}}|{\mathbf{X}}_{\mathrm{m}},\boldsymbol{\theta}, \sigma, {\sigma}_n\right)=-\frac{1}{2}{\mathbf{y}}_{\mathrm{m}}^T{\mathbf{C}}_{\mathrm{m}}^{-1}{\mathbf{y}}_{\mathrm{m}}-\frac{1}{2}\log \left(\det \left({\mathbf{C}}_{\mathrm{m}}\right)\right)-\frac{N}{2}\log \left(2\pi \right), $$

where $ {\mathbf{C}}_{\mathbf{m}} $ is the covariance matrix of the noisy measurement targets, as introduced in Eq. (11). The terms in this expression allow for an interpretation: The first term quantifies the data fit. In contrast, the second term is independent of the measured values $ {\mathbf{y}}_{\mathrm{m}} $ and penalizes the model complexity while the third term represents a constant and does not affect the parameter choice. The parameter selection by maximizing the likelihood thus establishes a balance between model complexity and data fit (Rasmussen and Williams, 2006). It should be noted here that the adjustment of six parameters based on a few measurement points is not necessarily well conditioned. However, since the predictive capability of the surrogate models is regularly checked in the proposed procedure, an adequate choice of parameters is ensured. In the implementation used in this work, this auxiliary optimization task is handled by use of the Constrained Optimization by Linear Approximation algorithm (Powell, 1998).

Upon the observation of function realizations, the prior distribution is updated to a posterior distribution, which represents a refined estimate of the function in light of the acquired data. After the parameters of the Gaussian process have been adjusted to the available measurement data, it can be used for predictions of function values $ {\hat{\mathrm{f}}}_{\mathrm{p}} $ at any unseen location in design space $ {\boldsymbol{x}}_{\mathrm{p}} $ . Equation (11) implies that the distribution of noisy targets $ \mathbf{y} $ is again of Gaussian type. The joint distribution of a set of measured training values that are subject to noise $ {\mathbf{y}}_{\mathrm{m}} $ and some noise-free test value $ {\hat{\mathrm{f}}}_{\mathrm{p}}, $ for which the covariance function is $ {\mathbf{K}}_{\mathrm{p}}={\sigma}_{\mathrm{f}}k\left({\boldsymbol{x}}_{\mathrm{p}},{\boldsymbol{x}}_{\mathrm{p}}\right)={\sigma}_{\mathrm{f}} $ can therefore be written as

(13) $$ \left[\begin{array}{c}{\mathbf{y}}_{\mathrm{m}}\\ {}{\hat{\mathrm{f}}}_{\mathrm{p}}\end{array}\right]\sim N\left(\mathbf{0},\left[\begin{array}{cc}{\mathbf{C}}_{\mathrm{m}}& {\kappa}_{\mathrm{p}\mathrm{m}}\\ {}{\kappa}_{\mathrm{p}\mathrm{m}}^T& {\mathbf{K}}_{\mathrm{p}}\end{array}\right]\right), $$

where $ {\kappa}_{\mathrm{p}\mathrm{m}}={\sigma}_{\mathrm{f}}^2{\left[k\left({\boldsymbol{x}}_{\mathrm{p}},{\mathbf{x}}_{\mathrm{m},1}\right),k\Big({\boldsymbol{x}}_{\mathrm{p}},{\mathbf{x}}_{\mathrm{m},2}\Big),\dots \right]}^T $ denotes the cross-covariance function between observations and prediction (Garnett, 2023). The above form allows to apply Bayes’ rule to derive a closed expression for the posterior, that is, the conditional probability of output $ {\hat{\mathrm{f}}}_{\mathrm{p}} $ at input $ \boldsymbol{x} $ , given the measured values $ {\mathbf{y}}_{\mathrm{m}} $ for inputs $ {\mathbf{X}}_{\mathrm{m}} $ :

(14) $$ {\hat{\mathrm{f}}}_{\mathrm{p}}\mid {\mathbf{y}}_{\mathrm{m}}\left({\boldsymbol{x}}_{\mathrm{p}}\right)\sim N\left({\mu}_{\mathrm{p}\mid \mathrm{m}},{\sigma}_{\mathrm{p}\mid \mathrm{m}}^2\right), $$

(15) $$ {\mu}_{\mathrm{p}\mid \mathrm{m}}\left({\boldsymbol{x}}_{\mathrm{p}}\right)={\kappa}_{\mathrm{p}\mathrm{m}}^T{\mathbf{C}}_{\mathrm{m}}^{-1}{\mathbf{y}}_{\mathrm{m}}, $$

(16) $$ {\sigma}_{\mathrm{p}\mid \mathrm{m}}^2\left({\boldsymbol{x}}_{\mathrm{p}}\right)={\sigma}_{\mathrm{f}}-{\kappa}_{\mathrm{p}\mathrm{m}}^T{\mathbf{C}}_{\mathrm{m}}^{-1}{\kappa}_{\mathrm{p}\mathrm{m}}. $$

Equation (15) reflects how the prior mean (0 in the present case) is updated through observations $ \mathbf{y} $ . It yields the most likely function value at $ {\boldsymbol{x}}_{\mathrm{p}} $ , $ {\mu}_{\mathrm{p}\mid \mathrm{m}} $ , which serves as a prediction from GPR and is used as surrogate model output after sufficient data are sampled. Equation (16) in turn describes a data-augmented variance function, which results from the prior covariance for $ {\boldsymbol{x}}_{\mathrm{p}} $ , $ {\sigma}_{\mathrm{f}} $ and an update term that depends on the proximity of $ {\boldsymbol{x}}_{\mathrm{p}} $ to the measurements in the training data. The variance $ {\sigma}_{\mathrm{p}\mid \mathrm{m}}^2 $ provides a degree of predictive confidence, which will serve to select designs for experiments and efficiently refine the model, as discussed below.

The posterior predictive distribution (Eqs. 14–16) can be evaluated for any input $ \boldsymbol{x} $ on the continuous domain $ \Omega $ . In the context of the study, however, we intend to use it to sample integer values $ \mathbf{x}\in {\Omega}_{\mathbf{x}} $ only, as these are the only evaluable inputs for the combustion experiments. Approximations (e.g., rounding) must therefore be applied to use the Gaussian process models for sampling. Garrido-Merchán and Hernández-Lobato (2020) showed that simple rounding can cause the GPR-based sample selection to get stuck and suggested an improved methodology in which they transform the inputs before computing the covariance function. The approach is implemented in the utilized SMT software package and has been applied here (Bouhlel et al., 2019). An example of the integer-based regression is also shown in Figure 5.

Figure 5. One-dimensional illustration of the two strategies used for sampling of training data. Uncertainty sampling utilizes the GPR prediction, obtained from measurement data $ {\mathbf{x}}_{\mathrm{m}} $ to select the subsequent training sample $ {\mathbf{x}}_n $ at the location of highest $ {\sigma}_{\mathrm{p}\mid \mathrm{m}} $ , while in random sampling all $ {N}_{\mathrm{rs}} $ samples are determined beforehand.

3.2. Sampling strategy for model training

The described method of GPR can generate a surrogate model based on a given training dataset. Two sampling strategies are used to select designs for the training data. They are summarized graphically in Figure 5.

3.2.1 Random sampling for computational-efficient surrogate model training

One approach is to randomly sample a fixed number of designs $ {N}_{\mathrm{rs}} $ using Latin hypercube sampling (LHS), measure the corresponding target values and then apply GPR to obtain a surrogate model.

(17) $$ {\mathbf{X}}_{\mathrm{n},\mathrm{rs}}={\mathtt{LHS}}_{\Omega_{\mathrm{x}}}\left({N}_{\mathrm{rs}}\right), $$

where $ {\mathtt{LHS}}_{\Omega_{\mathrm{x}}}\left({N}_{\mathrm{rs}}\right) $ describes the operation of drawing $ {N}_{\mathrm{rs}} $ samples from $ {\Omega}_{\mathrm{x}} $ following the LHS method described by McKay et al. (2000). This approach requires only one GPR step (i.e., parameter adoption through maximization of the log marginal likelihood in Eq. (12)) to be performed, which saves computational costs. A disadvantage, however, is that the number of samples required must be estimated in advance. This means that there is a risk of unnecessarily high measurement effort. Accordingly, this strategy is particularly suitable for variables for which the acquisition time is less than the time required for a regression step so that the savings in computing time are likely to outweigh the disadvantage of excessive measurements. In Figure 5, the random sampling method is illustrated using a one-dimensional function with $ {N}_{\mathrm{rs}}=6 $ measurement samples selected at once.

The accuracy of the surrogate model trained this way can be tested by K-fold cross-validation. Therefore, the set of training data is divided into $ K $ randomly chosen subsets of $ {N}_k $ samples each. A surrogate model is trained K-times: $ K-1 $ subsets serve as training data, and the remaining set is used to assess the model accuracy by calculating a mean deviation between prediction and validation data. The overall relative model error $ {e}_{\mathrm{rs}} $ of the surrogate model trained on the entire data $ {\hat{\mathrm{f}}}_K $ is estimated by averaging over the $ K $ validation tests,

(18) $$ {e}_{\mathrm{rs}}\left({\hat{\mathrm{f}}}_K\right)=\frac{1}{K}\sum \limits_{k=1}^K\frac{1}{N_k}\sum \limits_{b=1}^{N_k}\left|\frac{{\hat{\mathbf{y}}}_k^{(b)}-{\mathbf{y}}_k^{(b)}}{{\mathbf{y}}_k^{(b)}}\right|,\hskip2em {\hat{\mathbf{y}}}_k={\hat{\mathrm{f}}}_{K\backslash k}\left({\mathbf{X}}_k\right). $$

Here, $ {\hat{\mathrm{f}}}_{K\backslash k} $ is a surrogate model trained on all $ K $ subsets except subset $ k $ and $ \left\{{\mathbf{X}}_k,{\mathbf{y}}_k\right\} $ is the subset $ k $ used for validation. As before, $ {\mathbf{y}}_k^{(l)} $ denotes the lth component of vector $ {\mathbf{y}}_k $ .

3.2.2 Uncertainty sampling for measurement-efficient surrogate model training

For expensive target functions, it is useful to perform a regression step (i.e., adopting the parameters through maximization of the log marginal likelihood, Eq. (12)) after each measurement in order to extract the maximum amount of information from the available data and use it to select the next measurement sample. In this way, the number of measurements can be kept to a minimum. To select samples, surrogate model predictions from GPR are used to form a sampling criterion, which is expressed in terms of an acquisition function. The next point in search space to be evaluated $ {\mathbf{x}}_{\mathrm{n},\mathrm{us}} $ is then determined by locating the maximum of the acquisition function $ a\left(\mathbf{x}\right) $ ,

(19) $$ {\mathbf{x}}_{\mathrm{n},\mathrm{u}\mathrm{s}}=\underset{\mathbf{x}\in {\Omega}_{\mathbf{x}}}{\mathrm{argmax}}\ a(\mathbf{x}). $$

In classical BO, the aim is to find an optimum of the unknown function within a minimal number of function calls, which makes BO especially suitable for objective functions that are expensive to evaluate. The strategy behind the choice of acquisition function is then to strike a balance between exploration and exploitation. Here, exploration involves taking samples in areas of the search space where the function is uncertain or poorly understood, while exploitation focuses on areas where low values can be expected based on available knowledge. Various acquisition functions have proved highly effective to fulfill this purpose like the expected improvement (Mockus, 1998; Jones et al., 1998) or the lower confident bound (Srinivas et al., 2009), among many others (Forrester and Keane, 2009; Hennig and Schuler, 2012; Blanchard and Sapsis, 2021).

Surrogate models generated during BO therefore generally only capture global trends throughout the design space and are only quantitatively precise in promising areas where more samples are selected during exploitation. This characteristic enables sample-efficient optimization with BO (Garnett, 2023). In contrast to BO, for the purpose of this study, we do not aim to directly minimize an expensive objective function. Instead, our goal is to train surrogate models by minimizing the deviation between the surrogate model prediction and measured values of the target function, using as few measurements as possible, Eq. (4). As the optimization is only carried out a posteriori in combination with other models, we cannot select promising areas during the training process and thus require the surrogate models to be quantitatively precise throughout the design space. To train the surrogate models, the objective of the sampling strategy is therefore completely geared toward exploration and the acquisition function is formed by the prediction’s variance only,

(20) $$ a\left(\mathbf{x}\right)={\sigma}_{\mathrm{p}\mid \mathrm{m}}^2\left(\mathbf{x}\right). $$

This means that the method always tests the burner design for which the predictive uncertainty of the surrogate model is greatest and subsequently adds it to the training set. The strategy is shown on the left in Figure 5. Here, an unknown function (dashed line) that is evaluated at five points (blue markers) is approximated using integer-based GPR. The prediction and uncertainty are shown as a solid line and a gray shade, respectively. The next measurement point $ {\mathbf{x}}_{\mathrm{n},\mathrm{us}} $ is selected at the point of highest uncertainty.

For sampling, the maximization task in Eq. (19) must be solved using the surrogate model from GPR. In general, this can be done within another auxiliary optimization loop as evaluation of the GP posterior is computationally cheap. In the present case, however, the number of designs $ \mathbf{x} $ to consider is finite due to the integer nature of the design space $ {\Omega}_{\mathbf{x}} $ . Therefore, the acquisition function for all designs $ \mathbf{x} $ in $ {\Omega}_{\mathbf{x}} $ is evaluated for sampling and the maximum is determined by direct comparison, which ensures that the global maximum is selected.

The use of uncertainty sampling allows for a threshold-based convergence criterion by monitoring the model accuracy during training. For this purpose, the relative training error $ {e}_{\mathrm{us}} $ is calculated in each step by comparing the predicted value of the target function with the subsequently measured one,

(21) $$ {e}_{\mathrm{us}}\left({\hat{\mathrm{f}}}_n\right)=\left|\frac{{\hat{y}}_n-{y}_n}{y_n}\right|,\hskip1em {\hat{y}}_n={\hat{\mathrm{f}}}_{n-1}\left({\mathbf{x}}_n\right). $$

Here, $ {\hat{\mathrm{f}}}_{n-1} $ is the surrogate model trained with training data available in iteration $ n-1 $ and $ \left\{{\mathbf{x}}_n,{y}_n\right\} $ is the tested design and corresponding measured value in iteration $ n $ . In practice, this error metric is calculated in each iteration and the training is considered to be converged if $ {e}_{\mathrm{us}} $ falls below a predefined threshold value $ {\unicode{x025B}}_{\mathrm{us}} $ in 10 consecutive iterations. After fulfillment of the convergence criterion, the accuracy of the models trained through uncertainty sampling is assessed on additional cross-validation data $ {\left\{{\mathbf{x}}_{\mathrm{cv},i},{y}_{\mathrm{cv},i}\right\}}_{i=1}^{N_{\mathrm{cv}}} $ . Therefore, a relative cross-validation error $ {e}_{\mathrm{cv}} $ is estimated from comparison of the measured data with predictions from the final surrogate model $ {\hat{\mathrm{f}}}_n $ ,

(22) $$ {e}_{\mathrm{cv}}=\frac{1}{N_{\mathrm{cv}}}\sum \limits_{i=1}^{N_{\mathrm{cv}}}\left|\frac{{\hat{\mathrm{f}}}_n\left({\mathbf{x}}_{\mathrm{cv},i}\right)-{y}_{\mathrm{cv},i}}{y_{\mathrm{cv},i}}\right|. $$

It should also be mentioned that, in general, approaches are also possible that fall between the extremes of only one regression step after all data are acquired (random sampling) and a regression step after each measurement (uncertainty sampling). Methods for selecting batches of measurement samples have been proposed in the context of BO, for example, by Nguyen et al. (2016) and Lyu et al. (2018). However, in order to avoid the need to select additional hyperparameters in this article, the sampling procedure is limited to the two methods mentioned above.

4. Results

The described method is used to train surrogate models for three target variables that are relevant for GT combustor design. In the following, the training process of the three surrogate models is discussed before we demonstrate how the models can be used for multi-objective design optimization.

4.1. Surrogate model training

4.1.1 NO _x emissions under full-load-equivalent conditions — $ {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}:\mathbf{x}\to {c}_{{\mathrm{NO}}_x} $

The first design objective is to minimize NO _x emissions under full-load operation. The formation of NO _x emissions is associated with the dissociation of stable N₂, which is why these emissions are particularly critical at high flame temperatures, as they occur at full load. To address this in the surrogate model approach, the test rig is operated at full-load-equivalent conditions and the concentration of NO _x in the exhaust gas $ {c}_{{\mathrm{NO}}_x} $ is measured for different burner design schemes. To be close to realistic engine operation under full load, the pilot fuel ratio (PFR, share of total fuel flowing through the pilot line), is set to 10%. The flame shape under these operating conditions is shown in Figure 6a for an exemplary pilot design. The V-shaped swirl flame stabilizes in the shear layers and the outer recirculation zone, with the heat release zone extending to the outer edge of the combustion chamber. It should be noted that NO _x formation varies with combustion chamber pressure, making direct transfer of atmospheric test values to GT conditions challenging (Biagioli and Güthe, 2007). Nevertheless, it is a common strategy to use atmospheric tests to assess burner designs for machine conditions (Sattelmayer et al., 1992). However, the extent to which the surrogate models trained under atmospheric conditions are also suitable for ranking various pilot designs in machine operation is not within the scope of this study. Instead, we emphasize at this point that the methodology presented here is also applicable to build surrogate models from tests under pressure.

Figure 6. Operating condition, model training and validation of the surrogate model mapping burner design $ \mathbf{x} $ to concentration of NO $ {}_x $ in the exhaust gas.

The change in NO _x formation with the burner design is related to the temperature distribution in the combustion chamber. A settling time of one minute after a design change is therefore found to be necessary to allow transients to fade out. To obtain reproducible values, the measured concentration is then averaged over a two minutes period such that a total evaluation time per measurement $ {t}_{\mathrm{eval}} $ of three minutes is required. In order to minimize the number of these extensive tests for model training, the above-described training method of uncertainty sampling (Section 3.2.2) is used after two measurements with randomly chosen burner designs for initialization. Convergence of the training process is detected with a threshold value $ {\unicode{x025B}}_{\mathrm{us}} $ of $ 5\% $ which is obtained after 76 iterations, as shown in the learning curve in Figure 6b. After the training process has converged, the model allows predictions for every possible burner design in the design space. For further validation of the model and the training method, model predictions and additional measurements are shown in Figure 6c. To enable visualization, the considered designs lie on two two-dimensional slices through the four-dimensional design space $ {\Omega}_{\mathbf{x}} $ —in the plots, two components of $ \mathbf{x} $ are varied while the others are set to zero. The model accurately replicates the observed trends in the measured values, demonstrating a high quantitative agreement. The average relative validation error between the measured and predicted values $ {e}_{\mathrm{cv}} $ is $ 4.4\% $ and the surrogate model can be considered reliable for use in design optimization.

4.1.2 Risk of flame extinction and emissions of unburned hydrocarbons — $ {\hat{\mathrm{f}}}_{\mathrm{FE}}:\mathbf{x}\to {r}_{\mathrm{FE}} $

If the gas turbine is to be operated at part load, the fuel mass flow is reduced and a leaner global equivalence ratio is attained. With reduced equivalence ratio, the combustion velocity is decreased and the flame is in risk of being extinguished or blown out. Under these conditions, reliable operation can only be ensured by an increased pilot fuel ratio and an appropriately designed pilot injection unit.

To investigate the influence of the pilot injection on the combustion process under part-load operation, the operating parameters (volumetric airflow, equivalence ratio, and preheating temperature) are adjusted to the respective conditions from the GT. The operating conditions adapted accordingly are expressed in the following in terms of the relative load $ {l}_{\mathrm{r}} $ . Figure 7 shows the characteristics of the combustion process for different burner designs and pilot fuel ratios. A series of OH* images at the top of the figure illustrates the change in flame shape with load variation. The two plots below show emissions of unburned hydrocarbons (top) and corresponding exhaust gas temperature (bottom) for a load variation for two different burner designs and two different pilot fuel ratios. To be able to compare the exhaust gas temperatures despite the different operating conditions, the values have been normalized to the respective theoretic adiabatic flame temperature, $ {\mathrm{T}}_{\mathrm{ad}} $ . The left end of each graph corresponds to the minimal relative load at which the configuration could be operated, such that the length of the graphs indicates the width of the operating window. It is evident that the burner design has a significant influence on the operating window (and thus the risk of extinction at an operating point) and potentially enables operation with a lower pilot fuel ratio. Furthermore, the graphs show that flame extinction is preceded by incomplete combustion, which manifests by a strongly shortened flame and increased concentration of UHCs in the exhaust gas. The incomplete combustion and thus the proximity to the lean extinction limit is also indicated by a reduced exhaust gas temperature, as shown in the bottom plot of Figure 7, which sharply drops to below 0.5 $ {\mathrm{T}}_{\mathrm{ad}} $ when approaching flame extinction. The proximity to the extinction limit and thus the risk of flame extinction $ {r}_{\mathrm{FE}} $ can therefore be measured equally by the concentration of emissions of UHCs and the exhaust gas temperature,

Figure 7. Influence of design and pilot fuel ratio on part-load combustion characteristics. Top: Time-averaged, Abel deconvoluted flame images for $ \mathbf{x}=[0,0,0,14] $ , $ \mathrm{P}\mathrm{F}\mathrm{R}=45\% $ . Bottom left: Concentration of unburned hydrocarbons $ {c}_{\mathrm{UHC}} $ and exhaust gas temperatures $ {\mathrm{T}}_{\mathrm{exh}} $ , normalized by adiabatic flame temperature $ {\mathrm{T}}_{\mathrm{ad}} $ . Graphs end where flame is extinguished. Bottom right: Exhaust temperature as a function of UHC emission for different burner designs.

(23) $$ {r}_{\mathrm{FE}}\sim {c}_{\mathrm{UHC}}\sim {\left({\mathrm{T}}_{\mathrm{exh}}/{\mathrm{T}}_{\mathrm{ad}}\right)}^{-1}. $$

This insight is used to efficiently train a surrogate model that represents the ability of individual designs $ \mathbf{x} $ to avoid extinction. Unlike emission measurements, exhaust gas temperature responds rapidly to changes of the fuel injection scheme and does not require a long settling or measurement time. The time required to evaluate a single burner design $ {t}_{\mathrm{eval}} $ is therefore around 10 seconds. The short evaluation time makes a sophisticated iterative sampling search impractical, which is why the training data are acquired form random sampling, as described in Section 3.2.1. The model is trained at a fixed operating point which is representative of the width of the operating window. A relative load $ {l}_{\mathrm{r}} $ of $ 25\% $ and a pilot fuel ratio of 45% are selected for this purpose. This operating point has been chosen to ensure that the burner design has sufficient influence on the exhaust gas temperature while the flame is not extinguished for any configuration so that the tests can be run in a fully automated manner. Additional evidence of the indicative significance of the exhaust gas temperature under these conditions is provided at the bottom right of Figure 7, in which the emission of UHC is plotted against the measured exhaust gas temperature for various burner designs. The plot underlines that incomplete combustion (high UHC concentration) can be detected by low exhaust gas temperature. We use $ {N}_{\mathrm{rs}}=600 $ measurement points to train the Gaussian process surrogate model. The model’s accuracy is evaluated by $ K $ -fold cross validation, Eq. (18), where $ K $ is set to 10. The test revealed a validation error $ {e}_{\mathrm{rs}} $ of $ 4.9\% $ . We therefore conclude that the surrogate model is able to rank burner designs with respect to the associated operating window with sufficient accuracy. As an example, Figure 8 shows the predictions for the two slices through $ {\Omega}_{\mathbf{x}} $ already depicted above. The plots show a range of values between 0.4 and 0.6, which appears to correspond to conditions close to the extinction limit and complete combustion based on a comparison with Figure 7. Furthermore, similar to the NO _x emissions shown in Figure 6c, both slices largely show an increase in exhaust gas temperature for more open injectors. This suggests that the goals of low NO _x emissions under full load and complete combustion in part-load operation are at least partially contradictory, so that careful consideration is required in order to strike a balance.

Figure 8. Predictions of exhaust gas temperature for two cuts through the four dimensional parameter space; injection through Rings I and IV (left) and through Rings II and III (right).

4.1.3 Thermoacoustic oscillations — $ {\hat{\mathrm{f}}}_{\mathrm{TA}}:\mathbf{x}\to \hat{q} $

In addition to a wide operating window and low NO _x emissions, many other design goals play a role in burner development; one of them is thermoacoustic stability. The following section will therefore demonstrate how the proposed approach can be extended to include the thermoacoustic behavior of different burner designs.

It should be noted at this point that thermoacoustic stability does not only depend on the flame properties that potentially change with operating pressure (Sabatino et al., 2018), but also on the eigenfrequencies and thus the geometry of the combustion chamber. Thus, in order to obtain an optimization result that is transferable to the GT, the acoustic boundary conditions in the test rig must be adapted to those of the GT—either by modifying the geometry or by active control, as suggested by Bothien et al. (2008). Another option would be to measure the flame response separately from the chamber acoustics and optimize it in combination with a model of the GT chamber acoustics, as suggested by Reumschüssel et al. (2024a). However, within the scope of the present study, none of the above measures were undertaken. Instead, we limit ourselves to the optimization of the thermoacoustic oscillations in the test rig in order to investigate the applicability of the optimization method in this context. We emphasize that the adjustment of the chamber geometry will not lead to any methodological differences.

Thermoacoustic instability describes the coupling of the unsteady heat release and acoustic oscillations in the combustion chamber, which degrades the performance of the system. One way to alter the flame’s response to acoustic oscillation and thereby prevent such unsteady phenomena is to change the location of fuel injection. To incorporate thermoacoustics as an objective into the combustor design, a surrogate model is trained that maps the pilot injection schemes to heat release oscillation amplitudes in the flame. The latter are quantified by measuring fluctuations of intensity of OH* chemiluminescense, $ {I}_{{\mathrm{OH}}^{\ast }} $ by means of the photomultiplier tube as described in Section 2.

In general, no operating range can be specified in which thermoacoustic fluctuations are most likely to occur, as this is specific to the combustion system. To demonstrate how this objective can be incorporated into the suggested surrogate model framework, the operating parameters of the test rig were therefore adapted to conditions, where thermoacoustic fluctuations are particularly prominent. For this purpose, the preheating temperature of the air was increased above the corresponding value for full load ( $ {l}_{\mathrm{r}}\sim 110\% $ ) and the pilot fuel ratio is set to 20%. After experimental acquisition, the $ {I}_{{\mathrm{OH}}^{\ast }} $ signal is transformed into an amplitude spectrum $ {\hat{I}}_{{\mathrm{OH}}^{\ast }}(f) $ using a Fourier transform. The maximum of the spectrum between 50 and 500 Hz, normalized by the mean value acts as the target quantity. Figure 9 shows four time signal examples for different burner designs and the corresponding amplitude spectra. The figure also shows a histogram of the captured time signals. Although the signals associated with low amplitudes in the spectrum show a normal distribution around the mean value, which can be associated with noise in a (marginally) stable flame, the histogram for design $ \mathbf{x}=\left[\mathrm{0,0,0,14}\right] $ resembles the typical bimodal distribution of limit cycle oscillations (Noiray and Denisov, 2016). It can therefore be concluded that the burner design can switch the flame between weakly stable thermoacoustics and limit cycle oscillations at the chosen operating point.

Figure 9. Data processing from OH* chemiluminescence signal captured by photomultiplier tube. Top: Exemplary section of a normalized time signal and histogram. Bottom: Corresponding amplitude spectrum for four different burner designs.

Repeatable and accurate detection of self-sustained flame oscillations is challenging, especially in the marginally stable regime where oscillations are mainly driven by noise. A total measurement time $ {t}_{\mathrm{eval}} $ of 20 seconds per design is therefore used for model training. The uncertainty sampling method was selected for sampling with an accuracy threshold $ {\unicode{x025B}}_{\mathrm{us}} $ of $ 10\% $ . Compared to the NO _x model, the threshold has been increased to account for the stochastic nature of the oscillations. The learning curve for the thermoacoustic model is shown in Figure 10a. The difficulty of training a model for thermoacoustics can also be seen from this graph. During training, the error reaches significantly higher values than for the NO _x emissions (note the logarithmic scale in Figure 10a) and the convergence criterion is only fulfilled after 178 measurements.

Figure 10. Model training and validation of the surrogate model mapping burner design $ \boldsymbol{x} $ to thermoacoustic oscillations $ {\hat{I}}_{{\mathrm{OH}}^{\ast }} $ .

As before, the surrogate model is validated with additional measurements, which were not included during model training. Figure 10b shows the model predictions in the form of surfaces and validation measurements using spheres. The varying distribution of the spheres in the diagram is indicative for the partially random nature of $ {\hat{I}}_{{\mathrm{OH}}^{\ast }} $ . However, the surrogate model captures the main trends and the average relative deviation between measurement and prediction $ {e}_{\mathrm{cv}} $ is $ 12.4\% $ and therefore considered sufficiently accurate.

4.1.4 Summary of surrogate model training

Before proceeding with a description of the surrogate models’ application to multi-objective design optimization, a summary of the overall training process of the three models is provided. Table 1 therefore lists the operating conditions for the experiments and the measurement effort required for the training. Specifically, the physical quantities used as indicators for the design targets, the associated times required for the evaluation of a single burner design and the number of designs used for model training are specified. As explained above, the choice of sampling strategy for each of the objectives was made in consideration of the evaluation time. The highest evaluation time per sample is required for the $ {\mathrm{NO}}_x $ tests. Accordingly, uncertainty sampling was applied. This involved an additional computing time of about 10 seconds on a workstation for each evaluation. However, the method allowed to limit the number of training samples to 76 such that the model could be trained within less than 4.5 hours in total. The model for flame extinction is based on measurements of exhaust gas temperature, which can be rapidly measured ( $ {t}_{\mathrm{eval}}=10\hskip0.1em \mathrm{s} $ ). Therefore, random sampling was applied as a training strategy, in which 600 samples were selected upfront for training. Despite the high number of samples, a total training time of 100 minutes could be achieved, leading to a prediction accuracy similar to that of the NO _x model. For thermoacoustic oscillations, the evaluation time is $ 20\hskip0.1em \mathrm{s} $ . Including the computing time of around 10 seconds for each step, the model took around 1.5 hours to train. With a total amount of around 8 hours, the training of all models could thus be realized within 1 day of automated tests.

Table 1. Operating conditions and parameters of surrogate model training. Abbreviations us and rs refer to uncertainty sampling and random sampling

4.2. Surrogate model-based analysis and optimization

With the surrogate models trained as described above, the individual pilot designs can be evaluated with regard to NO _x emissions, the width of the operating window, and the influence on the thermoacoustic heat release oscillations, each under the respective most critical operating conditions. Thus, three target parameters can be assigned to each pilot design x. The fast evaluation time of the models enables efficient simultaneous consideration of the parameters. In the following, we first make use of the surrogate models to gain some general insights into the combustor characteristics and the influence of the individual design parameters. Subsequently, we discuss how a multi-objective combustor optimization can be conducted using the surrogate models.

4.2.1. Influence of individual design parameters

To evaluate the influence of the four design parameters, that is, the number of opened fuel injectors on the four rings, on the three target variables Figure 11 shows model predictions $ {\mu}_{\mathrm{p}\mid \mathrm{m}} $ with error bars indicating a confidence interval of $ \pm 2{\sigma}_{\mathrm{p}\mid \mathrm{m}} $ . Each parameter $ {x}_i $ is examined by considering designs for which all other rings are closed ( $ {x}_k=0 $ for $ k\ne i $ ) while the number of injectors on Ring $ i $ is varied. The left panel of Figure 11 shows the predictions of NO _x emissions. All four parameter variations indicate that more open injectors lead to a higher NO _x concentration in the exhaust gas. This trend is most probably associated with a reduced quality of mixing of the pilot fuel with the passing gas mixture prior to reaction, which increases the emission of NO _x. Fuel momentum and thus mixing quality are reduced by increasing the number of open injectors while maintaining a constant total pilot flow. The rating between the four injector rings can also be associated with fuel-air mixing. The injectors on Ring I introduce fuel at the furthest distance from the flame front, which enables a high level of homogenization prior to the reaction leading to the lowest NO _x values. Injection through Ring IV results in the shortest mixing path lengths before the fuel reacts in the flame front and thus also generates the highest NO _x emissions. Despite their different radial positioning, the emission values for the variation of $ {x}_{\mathrm{II}} $ and $ {x}_{\mathrm{III}} $ are very similar. The angled injectors on Ring III presumably play a role here, shortening the length of the mixing path. Overall, the NO _x model thus provides a quantitative means of analyzing the influence of the fuel momentum, the injector position, and the injection angle on the emissions.

The aforementioned trends can be observed in a similar form for the part-load model. The corresponding model predictions are shown in the center panel of Figure 11. Before analyzing the plot in detail, it should be recalled that, in contrast to NO _x emissions, high values of T $ {}_{\mathrm{exh}} $ are desirable as they correspond to a most complete combustion and thus a low risk of extinction. With this in mind, it is quite obvious that the two targets of low NO _x emissions and secure part-load operation are contradicting. For an increasing number of open injection nozzles, an increase in the exhaust gas temperature is predicted for most of the parameter range of all four variations. The comparison of the four injector rings is also approximately inverse to the rating of NO _x emissions, with Ring IV offering the best performance, while injection through Ring I appears to produce a flame which is most susceptible to extinction. The contradiction to the NO _x emissions is presumably due to the fact that a less homogeneous fuel-air mixture leads to a richer local stoichiometry and thus to higher combustion temperatures. This increases the anchoring of the flame and reduces the risk of extinguishing. However, unlike the NO _x concentration, significantly different values of exhaust gas temperature are predicted for Rings II and III at part load. This difference could be explained by changes in flame shape at part load, which presumably result in different mixing path lengths compared to full-load operation. In any case, the surrogate model allows for accurate prediction without the need for further physical insight.

The predictions for thermoacoustic oscillations are shown in the right panel of Figure 11. Although higher oscillations are mostly predicted for more open injectors, the comparison of the individual fuel rings shows a significantly different characteristic than for the other two models. Oscillations above 10% of the mean value (and therefore probably related to instability, see discussion of Figure 9) are only predicted for fuel injection through Ring IV.

4.2.2. Design optimization through combined model analysis

The predictive capabilities of the surrogate models go beyond the variation of individual parameters. They enable the quantitative forecast of the target variables for each design in the parameter space and thus allow for multi-objective optimization. Several parameters can be taken into account in different ways.

One possible approach is a weighted optimization with respect to all three objectives predicted from the surrogate models $ {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x},{\hat{\mathrm{f}}}_{\mathrm{FE}},{\hat{\mathrm{f}}}_{\mathrm{TA}} $ , for example, in the form

(24) $$ {\mathbf{x}}_{\mathrm{p}}=\underset{x\in {\Omega}_{\mathbf{x}}}{\mathrm{argmin}}\hskip1em {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}\left(\mathbf{x}\right)-{\alpha}_1{\hat{\mathrm{f}}}_{\mathrm{FE}}\left(\mathbf{x}\right)+{\alpha}_2{\hat{\mathrm{f}}}_{\mathrm{TA}}\left(\mathbf{x}\right),\hskip2em {\alpha}_1,{\alpha}_2\in {\mathrm{\mathbb{R}}}^{+}. $$

In this case, the optimum $ {\mathbf{x}}_{\mathrm{p}} $ depends on the choice of the weights $ {\alpha}_1 $ and $ {\alpha}_2 $ . The entirety of all optimal designs for every possible combination of weights is contained in the set of Pareto-optimal designs $ {\Omega}_{\mathrm{p}} $ , which is commonly defined as the set of all non-dominated designs. Here, a design $ {\mathbf{x}}_1 $ is said to dominate design $ {\mathbf{x}}_2 $ if it is superior or equal with regard to all target variables considered. We introduce the symbol $ \succ $ to denote this relation in the context of surrogate-model based combustor design:

(25) $$ {\mathbf{x}}_1\succ {\mathbf{x}}_1\hskip1em \iff \hskip1em {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}\left({\mathbf{x}}_1\right)\le {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}\left({\mathbf{x}}_2\right),\hskip1em {\hat{\mathrm{f}}}_{\mathrm{FE}}\left({\mathbf{x}}_1\right)\ge {\hat{\mathrm{f}}}_{\mathrm{FE}}\left({\mathbf{x}}_2\right),\hskip1em {\hat{\mathrm{f}}}_{\mathrm{TA}}\left({\mathbf{x}}_1\right)\le {\hat{\mathrm{f}}}_{\mathrm{TA}}\left({\mathbf{x}}_2\right). $$

The set of Pareto-optimal designs $ {\Omega}_{\mathrm{p}} $ can then be defined as all designs $ \mathbf{x} $ in $ {\Omega}_{\mathbf{x}} $ , for which no design $ {\mathbf{x}}_d $ exists $ {\Omega}_{\mathbf{x}} $ , that dominates $ \mathbf{x} $ ,

(26) $$ {\Omega}_{\mathrm{p}}=\{\mathbf{x}\in {\Omega}_{\mathbf{x}}:\{{\mathbf{x}}_{\mathrm{d}}\in {\Omega}_{\mathbf{x}}:{\mathbf{x}}_{\mathrm{d}}\ne \mathbf{x},{\mathbf{x}}_{\mathrm{d}}\succ \mathbf{x}\}=\varnothing \}. $$

Here, $ \varnothing $ denotes the empty set. The entirety of all target quantities associated with the Pareto-optimal designs is referred to as the Pareto front $ P $ :

(27) $$ P=\{Y(\mathbf{x}):\mathbf{x}\in {\Omega}_{\mathrm{p}}\},\hskip1em Y(\mathbf{x})=[{\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}(\mathbf{x}),{\hat{\mathrm{f}}}_{\mathrm{FE}}(\mathbf{x}),{\hat{\mathrm{f}}}_{\mathrm{TA}}(\mathbf{x})]. $$

For three target variables, the optimization result in Eq. (24) depends on two weighting factors and thus the Pareto front can be approximated by a two-dimensional hyperplane in three-dimensional target space. Figure 12 shows the predictions from the three surrogate models for all 69,360 designs comprised in $ {\Omega}_{\mathbf{x}} $ as black markers. The axes in the diagram are orientated such that optimal configurations are found near the rear corner. Because of the integer steps of the components of x, $ {\Omega}_{\mathbf{x}} $ contains a finite number of designs, so that the Pareto curve can be determined by direct comparison. Since the prediction process from the surrogate models is computationally cheap, it is not necessary to apply an additional a posteriori optimization routine to identify optima. The Pareto front considering all three target quantities consists of 379 designs. The corresponding predictions $ P $ are marked with an orange outline in the three-dimensional plot. To further visualize the Pareto front, a two-dimensional spline was fitted to $ P $ . It is shown in the form of a blue surface in Figure 12 and represents a limit to the achievable performance. Taking only two target variables into account results in a smaller set of Pareto-optimal designs. This is shown on the projection planes in the figure. Here, the gray dots represent the target values for all possible designs, corresponding to the black three-dimensional markers. The Pareto fronts under consideration of the respective two variables are indicated by orange outlines.

Figure 11. Surrogate model predictions ( $ {\mu}_{\mathrm{p}\mid \mathrm{m}} $ ) and 95% confidence interval ( $ \pm 2{\sigma}_{\mathrm{p}\mid \mathrm{m}} $ ) of target values for fuel injection through individual rings ( $ {x}_k=0 $ for $ k\ne i $ ).

Figure 12. Model prediction for all designs in $ {\Omega}_{\boldsymbol{x}} $ (scattered black markers) and Pareto front (orange). Blue surface indicates three-dimensional Pareto front. Orange markers on projection planes show two-dimensional Pareto fronts considering only corresponding two target quantities.

Depending on the extent to which the targets contradict each other, the Pareto front contains a large number or only a few burner designs. For example, the burner designs for which low NO _x emissions are predicted also predominantly exhibit low flame oscillations, so that only four burners lie on the corresponding Pareto front (bottom projection plane in Figure 12). Aiming for a high exhaust gas temperature under part-load operation, on the other hand, contradicts the other two targets strongly. The Pareto fronts therefore include more designs, and a weighting must be chosen for a final design selection (left/right projection planes in Figure 12). Regardless of which target variables are to be considered, the surrogate models can be used to efficiently identify Pareto-optimal designs.

Another option to use the surrogate models is for constrained optimization. If some of the target variables do not need to be optimized, but instead are required to fulfill a certain threshold limit, the corresponding surrogate model can be used to limit the parameter space accordingly. For GT burners, the emission of NO _x emissions is typically restricted legally. A fixed limit value must be adhered to during full-load operation. In the following, we will therefore consider the Pareto optimization between the part-load operation and thermoacoustic stability under the constraint of limited NO _x emissions under full load, yielding the following constrained optimization problem:

(28) $$ \mathbf{x}=\underset{x\in {\Omega}_{\mathbf{x}}}{\mathrm{argmin}}\hskip1em {\hat{\mathrm{f}}}_{\mathrm{TA}}\left(\mathbf{x}\right)-\alpha {\hat{\mathrm{f}}}_{\mathrm{FE}}\left(\mathbf{x}\right),\hskip2em {\hat{\mathrm{f}}}_{{\mathrm{NO}}_x}\left(\mathbf{x}\right)<{t}_{{\mathrm{NO}}_x},\hskip2em \alpha \in {\mathrm{\mathbb{R}}}^{+}, $$

where $ {t}_{{\mathrm{NO}}_x} $ denotes the constraining threshold. This scenario is of interest, as reduced pulsations increase the service life of the components and load flexibility increases the machine’s range of application. From a manufacturing point of view, both objectives should therefore be realized as thoroughly as possible while fulfilling the legal restrictions for NO _x emissions. Figure 13 shows the predicted values of flame oscillations plotted against the corresponding predicted values of $ {\mathrm{T}}_{\mathrm{exh}}/{\mathrm{T}}_{\mathrm{ad}} $ . The colored markers indicate Pareto-optimal configurations with respect to these two design goals, which are identified taking into account different NO _x emission constraints. The Pareto boundaries are additionally approximated by dashed lines. The front without NO _x constraints is shown in blue and corresponds to the one drawn on the right projection plane of Figure 12. If the limit value for NO _x emissions is tightened, this comes at the expense of lower part-load flexibility, which is evident from the lower maximum achievable $ {\mathrm{T}}_{\mathrm{exh}} $ values. The selection algorithm increasingly sorts out designs that have a large number of open injectors on the inner rings, as exemplified in the figure. The achievable level of flame oscillations in the Pareto front is not influenced by the threshold value. This is consistent with the finding in Figure 12 that for the investigated combustion system, similar configurations result in a minimum of NO _x emissions and a minimum of thermoacoustic oscillations. Overall, this analysis shows another strength of the surrogate model-based design approach: once the models have been trained and validated, the requirements for the machine can be flexibly adapted during the design process. One possible application would be, for example, when GTs are developed for several countries in which different emission limits must be complied with, requirement-specific designs can easily be identified without the need for new measurements.

Figure 13. Surrogate model predictions for part-load performance and thermoacoustic oscillations. All burner designs in $ {\Omega}_{\mathbf{x}} $ are shown in gray dots and Pareto-optimal configurations under different levels of constraint for maximum NO_x emissions are marked in color.

5. Conclusions

We have proposed a data-driven approach to tackle the intricate engineering challenge of designing gas turbine combustors for safe, stable and low-emission operation under various load conditions. This design task is characterized by different targets, which have to be evaluated in complex measurements under different operating conditions and thus in separate experiments. In order to account for these constraints, the methodology is based on multiple, separately trained surrogate models that express the functional relationship between the burner design and the target variable using probabilistic correlations. The applicability of the method was demonstrated on full-scale atmospheric tests of an industrial burner. However, the method could be equally advanced to automated high-pressure tests. We presented two strategies for surrogate model training: random sampling, which minimizes computational effort but carries the risk of excessive measurement effort, and uncertainty sampling, which minimizes the measurement effort for training but entails a higher computational overhead. Depending on how much the time and cost required to measure a target variable limits the number of measurements, one of the strategies may be preferred. Specialized hardware, which was created using additive manufacturing, enabled the training of the surrogate models to be carried out automatically. Three models were trained representing design objectives under operating conditions where they are critical: one that predicts NO _x emissions under conditions that represent full-load operation of the GT, one model that predicts the ability of different burner designs to prevent lean extinction in part-load operation, and a model that predicts thermoacoustic oscillations under conditions, in which the system was found to be particularly susceptible. All three targets can be mapped using the proposed methodology. In the present case, the training of surrogate models that represent a design space comprising more than 69,000 possible burner designs could be completed within 1 day of automated testing time. The models can be used flexibly in different ways for multi-objective optimization. They can be employed to filter out designs from the large design space that are optimal in terms of differently weighted combinations of objectives. We also showcased an application in an a posteriori constrained optimization, in which one of the objectives limits the design space while the others are optimized.

The utilization of such an automated approach enables the identification of optimal designs within the design space, which is defined by the test geometry. This necessitates the a) choice of meaningful measurement metrics to represent the design goals and b) creation of specialized hardware to efficiently test different design variants. Although these two prerequisites must be created by experienced engineers, the identification of optima can be carried out in a purely data-driven manner.