Optimization of the Shape of a Hydrokinetic Turbine’s Draft Tube and Hub Assembly Using Design-by-Morphing with Bayesian Optimization

Renewable energy is fundamental to meeting our energy demands in a sustainable fashion. Although there exist several renewable sources for clean energy (wind, wave, solar, etc.), hydroelectric power is perhaps the most pertinent as a viable replacement for petroleum, natural gas and other fossil fuels. In fact, hydropower was one of the largest renewable energy sources in the United States in 2021 [1]. As such, it is important to consider ways in which the hydroelectric power plant (HEPP) can be more efficient, and thereby more cost-effective.

The relative size of HEPP’s is usually classified by the amount of power they generate. Although there exists no unique definition of “small” hydropower, it is often accepted as a generating capacity at or below 10 megawatt electrical (MWe) [2]. The demand for “small” hydropower is steadily increasing [3] despite concerns about its potential adverse environmental impacts (but impact studies and how the impacts scale with the size of the hydropower plant remain controversial [4, 5]). Currently, there are numerous undeveloped sites across the globe that have large potentials for efficient and sustained power generation via small hydropower [6]. The sites with the most potential and that are easiest to exploit are those with low-impact stream-reaches, existing non-powered dams, and sites with existing conduits [2]. Motivated by this potential for inexpensive and sustainable energy, we propose here a new methodology for the design of some of the parts of a small hydrokinetic turbine. As a specific demonstration of our methodology, we provide a new, more efficient design of the draft tube and hub assembly of a small, low-head Kaplan turbine.

The turbine component of the HEPP works by converting the potential and kinetic energy of the water entering the turbine into mechanical work, and then producing electricity via a generator [7]. As one of the oldest and largest sources of renewable energy, there exists a wide variety of hydrokinetic turbines, but a concept common to all of them is the dynamic pressure or hydraulic head $P+\rho{\bf v}^{2}/2$ (where $P$ is the static pressure, $\rho$ the density and ${\bf v}$ the velocity) of the water entering the turbine. Using Bernoulli’s principle (see $\S$ 2) we can relate this head to the amount of power an ideal turbine can produce and the physical properties of the hydropower system such as the relative heights of the dam, penstock, and tail water discharge [8]. For high head ranges, impulse turbines (e.g., Pelton) exploit only the velocity of the fluid across the runner (see Fig. 1) to create the mechanical of the turbine blade; whereas in medium and low head ranges, reaction turbines (e.g., Francis and, with later developments for low head applications, Kaplan) exploit both the fluid’s velocity and the fluid pressure or enthalpy across the runner. It is the conversion of enthalpy that allows the draft tube and hub assembly to enhance the performance of a reaction turbine [9, 10]. The draft tube and hub assembly make up only part of the overall hydropower system’s performance, but they are most important for low-head turbines, which we consider here to be those under 20 m.

The draft tube is a diffuser, or several diffusers joined together, that sits beneath the runner of the turbine and directs the flow downstream of the turbine blades to the tailwater pool. It therefore has a large role in determining the dimensions of the lower section of the power plant [11]. The draft tube increase the efficiency of the turbine by adjusting the dynamic pressure $(P+\rho{\bf v}^{2}/2)$ such that the static pressure $P$ just downstream of the turbine blade is decreased. The adjustment is done by decelerating the velocity of the fluid passing through the draft tube, and it is this property (hereafter referred to as the “pressure recovery”) that affects the power-generating capacity of Kaplan and other reaction turbines [12, 13].

The hub is a conically shaped part that extends just past the inlet of the draft tube, centered about the axis turbine’s of rotation and connecting its blades. Water flows into the draft tube through an annular-shaped region between the hub and the draft tube wall. The hub rotates with the same angular velocity as the turbine blade, while the outer boundary of the draft tube is non-rotating. (See Fig. 1.) It is important to optimize the hub because it modifies the inlet flow to the draft tube, and therefore has a large effect on the pressure recovery.

Poorly designed parts of the draft tube or hub can significantly decrease a turbine’s efficiency by exacerbating turbulence and increasing friction losses. For example, the draft tube’s elbow (Fig.1), which is necessary for redirecting the tailwater flow, can promote flow separation due to excessive centrifugal force at its inner radius. Similarly, a poorly designed hub can allow the swirl flow at the inlet due to the turbine blades create instabilities in the flow that lead to noise, vibration (prompting failure due to fatigue), and even the reversal of flow through the center of the draft tube (causing sudden changes in power output of the turbine) [14]. A well-designed hub can help prevent these problems, and, in addition, allow a larger opening angle of the diffuser (and therefore higher pressure recovery). As such, we not only optimize the draft tube, but also the hub, which is usually neglected in prior optimization studies [15, 16, 17].

Much of the optimization effort of Kaplan draft tubes in recent years was focused on improving the sharp-heeled elbow draft tube [17]. This shape was first installed in 1949 at the Hölleforsen Hydro Power Station in Vattenfall, northern Sweden, utilising a 25m head and with a power generating capacity of 150 MWe [18]. Gubin [17] and Dahlbäck [19] independently argued that this draft tube shape needed improvement, and subsequently there have been many proposed design changes. Marjavaara and Lundström [20] and Marjavaara [15] use a Response Surface Method (RSM) surrogate modeling strategy, as well as a commercial CFD code (ANSYS CFX4.4) to create new designs with different parametrizations (circular and elliptical, respectively) of the elbow section. More recent improvements include those of Daniels et al. [21, 22] who use multi-objective Bayesian optimization to maximize pressure recovery using a series of subdividing curves, optimizing over the inflow cone, outer-heel, and secondary straight diffuser. Other improved designs focus on the optimization of the draft tube for low-head applications while retaining the sharp-heel shape [23, 24].

Design search space creation or parametrization is an integral part of any shape optimization process. Parameterization determines both the design space and the complexity of the optimization problem. A desirable parameterization technique must cover a wide design space within a limited number of design parameters[25, 26, 27], necessary during the early design stage when minimum geometric constraints are placed and radical changes during the optimization process are required.

For a set number of parameters, parametrization techniques are judged on the different fidelity and ranges of control they offer [28, 27, 25]. Discrete methods [29], where design parameters are exactly the discrete surface points that define the airfoil shapes. The displacement of each point can be adjusted [30], and very fine local control can be achieved. However, to describe a shape accurately and with high fidelity, a large number of surface points are needed, which increases the number of parameters of the optimization problem. To accommodate the large number of design variables, gradient-based method are usually used to guide the optimization which can easily get stuck at a local optimum.

High-fidelity features can be captured by B-splines [31, 32], and nonuniform rational B-spline (NURBS) [33] which form curves by connecting low-order Bézier segments defined by the control points. The fidelity of shape representation with these techniques depend on the number of parameters used for curve definition, increasing the computational complexity. To reduce the number of the design parameters, the control points can be grouped together, free-form deformation (FFD) method [34, 35]. Other techniques include spectral representation techniques which use some basis functions or modes, such as proper orthogonal decomposition (POD) [36, 37], Hicks-Henne’s approach [38], and class/shape function transformation (CST) method [39, 40].

To address the design challenges for improved hydrodynamic and aerodynamic surfaces, we introduced a new framework, Design-by-Morphing (DbM), for creating a design space that is versatile enough to include old and new designs and that is sufficiently free of human bias to produce a radical and counter-intuitive design search space. The strategy cannot be called fully free from designer bias due to the choice of baseline shapes by designer, but the negative weights and radical baseline shapes, described later, minimize this effect. DbM uses the shape space as the basis for the parametrization scheme, thus allowing infinite fidelity without increasing the design space parameters for any shape, since it is the weight on the shape itself that is optimized. This is not possible with any conventional scheme discussed above and DbM has been shown to reproduce the entire UIUC database to very fidelity using just 25 parameters in [41]. DbM was first used by Oh et al. [42] and later used in other optimization problems [41, 43]. An $N$-dimensional design space for a shape is created by choosing $N$ baseline shapes. The shape within our design space is specified by the choice of the $N$ weights of these baselines from which the design is morphed. The bounds of the weights are sufficiently large that the morphed shapes can be not only interpolations among the shapes, but also extrapolations. Furthermore, any of the weights can be negative so that features of poorly performing baseline shapes can either be suppressed or entirely avoided. Negative weights and large positive weights allow for unintuitive shapes and for extrapolations. Many design techniques only allow for small departures from existing designs [44] or allow only local changes at one or a few specific locations, rather than global changes to the overall shape [45, 46]. This is especially problematic with most methods that use control points. CAD, and/or NURBS. Generally, methods with parametric control [28, 45, 47, 48, 49] and the adjoint method [44, 46, 50] limit the amount of change that can be made to a design so that radical new designs are not possible. DbM also allows a more extensive design space where both spatially local and global changes can be made to the design, and those changes can be subtle and/or radical [41]. This allows us to find an optimum which may be a completely new, unconventional shape [42].

Once a design space is chosen, many engineering optimization problems require the repeated numerical (or laboratory) evaluation of an expensive multi-modal black-box function to determine the performance or cost of a particular design. In optimizations of a surface that interacts with a fluid, generally the fluid flow must be computed with an expensive Computational Fluid Dynamics (CFD) program to determine the quantitative performance of each candidate design. Those quantitative results are then fed into an optimizing algorithm to determine the best performing design. The expense of computing the performance function with CFD makes many optimization problems intractable. Furthermore, the highly nonlinear behaviour of fluid flows often leads to a design space where the performance function has many local maxima, and it is difficult for the search algorithm to find the global maximum.

Optimizing draft tubes is an example of a search requiring an expensive multi-modal, black-box, performance function because the efficiency and pressure recovery of each candidate draft tube must be computed with CFD. Gradient based optimization techniques are unsuitable for multi-modal optimization since they get stuck in local optima and might require a number of function calls to either multi-start or get out of the local optimum. Global optimization schemes, such as evolutionary algorithm on the other hand require a large number of function calls to reach a global optimum, making them unsuitable for expensive black-box functions. Compared to Bayesian optimization, a genetic algorithm requires at least two orders of magnitude greater number of function calls to find the global optimum [51]. Bayesian Optimization (BO) is an efficient search method for this type of optimization because it requires far fewer evaluations of the performance function to find an optimum than most other optimization algorithms [52, 53]. BO techniques have been successful in the design of architected meta-materials [54, 55, 56, 57, 58, 59], hyper-parameter tuning for machine learning algorithms [60, 61, 62], drug design [63, 64], and controller sensor placement [65]. In this study, we use an improvement that we made to BO, that we call Mixed-variable Multi-Objective Bayesian Optimization (MixMOBO) [51]. Previously, we used MixMOBO to design an architected meta-material that has maximum strain-energy density [59].

The design space is created from five baseline draft tube shapes and two baseline hub shapes. The shapes are morphed with the weights chosen by MixMOBO to create a new draft tube and hub shape. These baselines are shown in Figure 2. The five baseline draft tube shapes are homeomorphic to each other as are the two baseline hub shapes. Some of the baselines that we chose were used previously in literature so that we were able to exploit any inherent advantageous features they might have, while other baselines were chosen to have non-intuitive features that would lead to radical designs, rather than incremental improvements. In the past, the fundamental design features of the Kaplan draft tube were formed through experimental observation and quasi-empirical formulae derived from geometries already installed and in use in HEPPs [21]. The works of Gubin [17], Cervantes [11], Mulu [67] and Nilsson et. al. [14] [68]) provide insights into draft tube geometries.

Our first baseline shape for the draft tube is the sharp heel draft tube, which has been the focus of extensive optimization attempts. It is also the subject of extensive experimental and numerical studies, the majority of which were completed through the European Research Community On Flow, Turbulence And Combustion (ERCOFTAC) Turbine-99 Workshop series [69, 70, 66]. For this reason, this shape also provides excellent validation of our CFD setup.

The second and third baseline draft tube shapes are based on designs cited by Gubin [17] and are particularly well-suited to low-head applications (i.e., for Kaplan turbines). The fourth and fifth baseline tube shapes were both devised in order to create features which would expand the design space. These shapes were generated by taking the average shape of the first three baseline shapes and then adding radically different features such as a rounded outlet and a circular diffuser shape. The addition of these baseline shape can be thought of in a similar manner to adding mutation within a generation in genetic algorithms. All the baseline draft tubes have the same shape at their inlets, and in all cases the planes containing inlet and outlet are perpendicular to each other.

Both baseline hubs have the same radius at the inlet, the same radius at the end of the hub, and the same length. The first baseline hub shape is based on one currently used in a Kaplan turbine [67, 71]. The second baseline hub is a cone shape, and which was historically used in a wide variety of low to medium head reaction turbines [17]. Morphing these two hubs allows for a variation in the inlet geometry, which in turn modifies the flow profile near the inlet and therefore the pressure recovery. The baseline hub shapes are shown in Figure 2.

Design-by-Morphing (DbM) works by creating a one-to-one correspondence between the set of baseline shapes. The shapes are then “morphed” i.e. linearly combined together using “weights”. These weights are the independent parameters that form the search space for optimizing the shape. DbM requires that the set of shapes be homeomorphic or topologically equivalent so that a one-to-one correspondence can be created between the surface elements. This one-to-one correspondence is created by defining a collocation strategy for the baseline shapes that ensures all constraints for our design are fulfilled (for example inlet and outlet orientation) even if individual shapes are radically different from each other.

For the draft tube baseline shapes (shown in Figs. 3 and 4), the inlet is constrained to a plane parallel to the horizontal ($x$-$y$) plane and the outlet is constrained to a plane parallel to the vertical ($y$-$z$) plane for all the baseline shapes. Note that this ensures that all the morphed shape inlets and outlets are similarly constrained i.e. inlets are in the same horizontal plane and the outlets are in the same vertical plane. Fig 3 explains how we construct a single coordinate system $(r,\varphi,s)$ for all of the baseline and morphed draft tubes that is mapped from the cylindrical coordinates $(r,\varphi,z)$. The $z$ axis is mapped to the origin curve in the figure with with arc length $s$ and local unit vector $\hat{\bf s}$. The boundary of the draft tube is defined by its intersections with each of the $N_{s}$ perpendicular planes shown in Fig. 3. Within each perpendicular plane, the intersection of the plane with the draft tube boundary is defined by the closed curve $r_{j}(\varphi)$, the radial distance of the draft tube boundary from the origin curve in the $j^{th}$ perpendicular plane. We note here that, because the same origin line is used for all of the baseline shapes, the start and end points of the origin line are always constrained to the same location in the inlet and outlet planes. The draft tube shape, however, changes around the origin line, meaning that the average path length for a fluid particle to travel from inlet to outlet, can vary for each draft tube.

In each of the $N_{s}$ perpendiclar planes, we discretize $\varphi$ into $N_{\varphi}$ equally-spaced collocation angles $\varphi_{k}$ with $k=0,1,2,\cdots,N_{\varphi}-1$, with $\varphi_{0}=0$ and $\varphi_{N_{\varphi}}=2\pi(N_{\varphi}-1)/N_{\varphi}$. The radius of the $p^{th}$ draft tube as a function of the arc length $s$ and polar angle $\varphi$ is completely defined by the radial location matrix $R^{p}_{k,j}\equiv r^{p}(\varphi_{k},s_{j})$, with $k=0,1,2,\cdots,N_{\varphi}-1$, and $j=0,1,2,\cdots,N_{s}-1$. Note that the $s_{j}$ perpendicular planes and the $\varphi_{k}$ collocation angles are the same for all of the the baseline and the morphed draft tubes.

For hub shapes, the same collocation strategy is used as we used for the draft tube shapes, with the origin curve of the hubs passing through the center of both of the baseline hubs because the baseline hub are axi-symmetric.

After morphing, each morphed radial matrix can be projected back into a 3-D shape as in Figure 4.

Because all of the five baseline draft tube shapes are homeomorphic, they can easily be combined into a new morphed shape, given by $r^{morph}(\varphi_{k},s_{j})\equiv R^{morph}_{k,j}$, once the weights $w_{p}$, $p=1,2,\cdots,5$ of each baseline are chosen (see Fig. 4):

Note that negatives weights and weights greater than unity allow for extrapolations; negative weights also allow us to “avoid” some baselines. This means that the only existing constraints are those imposed by choosing the baselines shapes themselves. That said, the lack of CAD parameterization means that DbM may produce non-physical, self-intersecting shapes when negative weights are applied. If the morphed shape is not physical because it has intersecting radial boundaries, we set the pressure recovery function of that morphed draft tube function to be zero so that the optimization method avoids that region of design space. For draft tube optimization, we limit the allowable range of the weights of each of the draft tubes such that: $w_{p}\in[-0.5,1.0]$.

The morphed shape of the hub is given by an equation similar to eq. (5), but the sum is over only 2 baseline radii. We denote the weights for hubs with $\alpha$ to differentiate it from the draft tube weights. Furthermore, the sum weight for the second hub baseline is constrained by $\sum_{l=1}^{2}\alpha_{l}=0.5.$, and for the morphed hub shape, $\alpha_{1}\in[-0.5,1]$. Due to the fact that sums of the weights of the draft tube baseline shapes and of the hub baseline shapes are fixed, there are only 6 degrees of freedom in choosing the values of the weights of the baseline shapes, so our design space has 6 dimensions.

The coefficient of pressure recovery $C_{prm}$ of each morphed draft tube shape that is physically allowable is determined with CFD. We validated the CFD code by comparing our simulations of the sharp heel draft tube with previously published results. In order to compute the flow in a morphed hub and draft tube shape, a mesh of that shape is created from a surface point cloud of the shape using Gmsh [72]. The average statistics of a typical mesh is given in Table 1, where D is the inlet diameter. The software used is OpenFOAM, and the solver used is pimpleFoam.

For turbulence modeling, we use the $k$-$\omega$ $SST$ turbulence model. This choice was made based on its success in previous draft tube studies [13, 73]. Based on the kinematic viscosity of water, the distance between the outer edge of the hub and the inlet wall, and the average streamwise velocity at the inlet, the Reynolds number is 5.56e+05, and based on the azimuthal velocity of the hub, it is 9.48e+05. No-slip conditions are applied at the rotating, inner hub wall ($62.3~{}rad~{}s^{-1}$) and at the non-rotating draft tube boundary.²²2The optimization studies in the past for draft tubes [23, 24], did not consider a rotating hub inlet condition, which makes the CFD boundary conditions less close to the actual operating or experimental conditions The inlet velocity is axisymmetric with non-zero values of the azimuthal, radial, and streamwise components. Our CFD simulations use the experimentally measured inlet velocities found by Engström et al. [70]. The outlet pressure is the pressure at location ${\bf(2)}$ in Fig. 1, and is fixed gauge pressure, but the pressure at the inlet, needed to determine the $C_{prm}$, is computed by the CFD code as a function of $r$, $\varphi$, and time. We run the CFD for 128,571 time steps for each morphed hub/draft tube shape. Our choice of the number of time steps is based on how long it takes the solution to reach a statistical steady state – see Fig. 5. Note that $C_{prm}$, given by eq. (2), is not determined at a single time step, but rather it is averaged over the final 28,500 time steps of the computation (and note that as shown in Fig. 5, the solution has converged to a statistically-steady equilibrium during those final time steps). The time step was chosen to be $0.00134~{}D/u_{avg}$, where $D$ is the outer diameter of the hub and $u_{avg}$ is the average streamwise inlet velocity. This time step was chosen based on a time resolution study (see below).

The CFD is validated in two ways. The first compares our CFD-computed pressures with the experimentally-measured values along the top and bottom center lines of the baseline 1, Sharp-Heel, draft tube. The second uses baseline 1 to test the convergence of the CFD code by refining the numerical time step and grid size. We note here that the experimental studies did not provide the $C_{prm}$ values, thus no $C_{prm}$ comparison could be made.

The local pressure recovery along the top and bottom lines of the first baseline shape (i.e., the sharp heel draft tube) is shown (and defined) in Figure 6. The figure compares our numerically-computed local pressure recovery values to the experimentally-measured values [70]. Figures similar to Fig. 6 appear throughout the ERCOFTAC Turbine-99 Workshop series, and the deviation between our CFD and experimental values are consistent with the previous numerical studies.

The second way we validated our use of our OpenFoam code and gridding scheme was by time step and mesh refinement. We chose our time step size by decreasing its value in the CFD code until the late time-averaged $C_{prm}$ values of the sharp heel draft tube with the first hub baseline shape changed by less than $0.01\%$. In particular, once we set the value of the time step as $0.00134~{}D/u_{avg}$, we repeated the calculation of the local pressure recovery $C_{p}(j)$ for both the top and bottom as a function of $j$ as defined in Fig. 6 – once with a time step $20\%$ greater than that used in the figure and once with it $20\%$ lower. Plots of these new curves computed with these time steps are indistinguishable from the curves (error less than 0.5%) shown in Fig. 6 (at the resolution at which we created the figure), validating that our time step is sufficiently small. We determined the size of the spatial resolution of the grid for the CFD in a similar manner: we decreased the grid size (or equivalently, increased the square of the number of grid points) until the late time-averaged values $C_{prm}$ values of the sharp heel draft tube changed by less than $0.4\%$. We also recomputed the two curves of the $C_{p}(j)$ in Fig. 6 with a spatial resolution of the spatial gridding mesh with a size that was $20\%$ greater than that used in the figure and once with it $20\%$ lower. Again, we found that the plots of these new curves computed with these grids are indistinguishable (error less than 0.5%) from the curves shown in Fig. 6, validating that our grid is sufficiently fine.

It is important to note that Fig. 5 shows that it takes more than 50,000 time steps for our numerically-computed $C_{prm}$ to reach a statistical equilibrium. The 50,000 time steps corresponds to a time that is approximately equal to 10 “advective times”, where the latter is the streamwise length of the draft tube divided by the average streamwise velocity at the inlet $u_{avg}$. Previous numerical studies of draft tubes were often validated by comparing the numerically-computed local pressure recovery $C_{p}(j)$ at the bottom or top of the draft tube with the experimentally measured values. However, in many of those studies the codes were run for only one advective time (5000 of our time steps) before the pressure recovery factor $C_{prm}$ was evaluated [70]. Figure 5 clearly shows that although that time be sufficient to numerically capture the final statistically steady state of $C_{p}(j)$, it is insufficient for the full flow and $C_{prm}$ to have come to a statistical steady state.

In our previous work, we developed a Mixed-variable, Multi-Objective Bayesian Optimization (MixMOBO) algorithm [51], a framework for optimizing mixed-variable, multi-objective problems with noisy black-box function using parallel batch updates of the data set, that is, letting the acquisition function pick several points to evaluate at each iteration to allow their CFD evaluations to be carried out in parallel. MixMOBO was proven to be more efficient in terms of number of black-box function evaluations compared to other algorithms in small data settings. In our previous studies, we used MixMOBO to design a microlattice structure with the objective of maximizing its strain-energy density [59]. We are also applying MixMOBO for the optimization of vertical axis wind turbines [43].

Our current study adapts the generalized MixMOBO algorithm for the specific case of 6 continuous variables with parallel batch updates using the HedgeMO algorithm [51]. Our HedgeMO strategy uses the Upper Confidence Bound, Expected Improvement, Probability of Improvement, and Stochastic Monte-Carlo [52, 51] acquisition functions. We apply our Bayesian optimization algorithm in tandem with DbM here for the shape optimization of a draft tube/hub. This section details our adaptation of the MixMOBO algorithm for our current study.

The search in our 6-dimensional design space (consisting of the 5 weights of the baseline draft tubes and the 1 independent weight of the two baseline hubs) for the draft tube/hub design with the maximum $C_{rpm}$ is an example of an optimization problem that can generically be posed as:

where $f(\vec{w})$ is the objective to be maximized (in this case $C_{prm}$), and $\vec{w}$ is a $D$-dimensional (in this case 6-dimensional) variable vector, defined over a bounded set $\mathcal{W}\subset\mathbb{R}^{D}$ (in this case, the weights of the baselines). For many practical engineering problems, $f(\vec{w})$ is expensive to evaluate (as is the case here, which requires repeated CFD evaluations of a complex draft tube designs). In such cases, a Bayesian optimization algorithm is often the best choice [52]. Bayesian optimization is a sequential optimization technique, specifically designed to find the optimal solution of a noisy black-box function $f$ with the fewest possible evaluations or function calls to $f$. At every iteration of the algorithm, a surrogate model $g$, usually a Gaussian process (GP), is fit over the data set $\mathcal{D}=\{[w_{1},f(w_{1})],\ldots,[w_{i},f(w_{i})]\}$.

Here $i$ is the total number of evaluated points (in this case, morphed draft tube/hub designs) and $\vec{w}$ is the vector of continuous variable vector.

The MixMOBO algorithm is agnostic to the choice of kernel function. For our current study, we use a simple modified squared exponential kernel due to its good performance for continuous variable problems:

where $\vec{\theta}=(\{{\underline{\underline{M}}}\},\epsilon_{f})$ is a vector containing all the hyper-parameters, $\{{\underline{\underline{M}}}\}=diag(\vec{h})^{-2}$ is the covariance hyper-parameter matrix and $\vec{h}$ is the vector of covariance lengths and distance between the variables defined to be their Euclidean distances.

Once the GP surface has been fit, the next query point in the design space that should be chosen to evaluation must explore the space that most probably contains the optimum (e.g., the design with the largest $C_{prm}$), and must also try to improve the fit for the next iteration, i.e. choose a region of space where information is sparse. These two determinations are not the same, and the former is called “exploitation” and the latter “exploration”. Bayesian optimization strives to create a balance of exploitation and exploration that continues until there is either evidence that a global optimum has been found or a maximum prescribed number of iterations of the algorithm is reached. An acquisition function, which we represnt here by $\mathcal{A}^{l}$, determines the next query point (or points) $\vec{w}_{i+1}$ in the design space to be evaluated with $f$ (in this case, with the costly CFD), by balancing the competing needs of exploitation and exploration. Once the next query point (or points) has been determined , $f$ is evaluated for that point (or points) and is then appended to the data set, $\mathcal{D}=\mathcal{D}\cup(w_{i+1},f(w_{i+1}))$. To optimize the exploitation/exploration balance, we use a GA to optimize the acquisition functions, which, although expensive to use on an actual black-box function, is an ideal candidate for optimizing the acquisition function working on the surrogate surface.

MixMOBO uses HedgeMO, a hedging strategy for acquisition functions, as part of its framework. Hedge strategies use a portfolio of acquisition functions, rather than a single acquisition function [74]. HedgeMO also allows parallel ‘$Q$-batches’ of query points. The detailed algorithmic flowchart of our adapted MixMOBO algorithm is shown in Algorithm 1. We note that the convergence criterion presented by [74] holds for our portfolio as well.

The process is repeated until the evaluation budget is reached. Our choice for the evaluation budget was based on several optimization experiments on test functions and design spaces, described below, that we believe to be representative of the draft tube/hub optimization. To determine the number of iterations of the Bayesian optimization algorithm, or epochs (where each epoch determines the $C_{prm}$ of 5 morphed draft tube/hub designs in parallel), required to likely find the optimal design within the design space, we carried out optimizations of a suite of test functions with different properties and whose maximum values could be determined analytically. We optimized the Spherical, Rastringin, Syblinski-Tang and Amalgamated functions, all of which are standard functions used to test optimization schemes [75] with the exception of the Amalgmated function, which is novel and created by us to mimic some of the properties of the draft tube/hub design space. The test functions that we used to determine the number of epochs that are needed for MixMOBO to likely find the maximum $C_{prm}$ are defined in Table 3. Similar to the design space of the draft tube/hub, each test function was tested with 6 dimensions.

Figure 7 shows how MixMOBO approaches the global maximum of each of the four test functions as a function of epoch. Two of the four test functions find the global maximum within 75 epochs.

The optimization was terminated after 75 epochs because the test functions reached reasonable convergence, other than Rastringin function, which is know to have a a multitude of local optima, as shown in Fig. 7. We could only afford 75 epochs worth of CFD evaluations of candidate draft tube designs in our computational budget so the evaluations were not carried past 75 epochs. As depicted in our Results, the draft tube optimization converged to the final design in 30 epochs. MixMOBO, like most Bayesian optimization schemes, needs to be initialized with evaluations of random designs. Based on the results shown on Test Functions, we use 50 random evaluations, and 75 epochs with 5 parallel batch evaluations per epoch to optimize our 6-dimensional design space.