In June of 2023, I defended my thesis project to obtain both a master’s degree in Marine Engineering and Computer Science, awarded an 8.5 in both disciplines. It is adapted as a paper and published open-access in the International Shipbuilding Progress journal.

This post only discusses those things that are interesting but could not be shown by the paper. The complete thesis titled “Multi-Fidelity Kriging Extrapolation” can be found on the TU Delft repository.

Before reading this article, it might be good to read up on the following subjects (from the post series “Master’s thesis”):

• Surrogate modelling: what is surrogate modelling and what is its typical application in engineering?
• Kriging, also known outside engineering as Gaussian Process Regression, is the surrugoate modelling technique at the core of the work of this thesis due to some beneficial characteristics with respect to its use together with expensive-to-evaluate objective functions.
• Multi-fidelity Kriging: extending the Kriging surrogate to beneficially combine multiple infromation sources of different precision and cost.
• (Multi-fidelity) Efficient Global Optimisation: an algorithm that uses the characteristics of (multi-fidelity) Kriging to converge to the global optimum solution of a design problem.

The master’s thesis introduces a novel multi-fidelity methodology that could increase (design) optimisation efficiency and effectivity in deterministic computer experiments. The methodology is among others tested in a design optimisation case that mimics a lifeboat impacting with a still body of water, where the input data is acquired using Computational Fluid Dynamics simulations.

Pedagogical example

To show the use and effectiveness of the proposed method, I would like to re-iterate the example of Citation: Forrester, Sóbester & al. (2007) Forrester, A., Sóbester, A. & Keane, A. (2007). Multi-fidelity optimization via surrogate modelling. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463(2088). 3251–3269. https://doi.org/10.1098/rspa.2007.1900 , and then slightly change their multi-fidelity transformation function. The promise of Multi-Fidelity Kriging (MFK) per the contrived example of Citation: Forrester, Sóbester & al. (2007) Forrester, A., Sóbester, A. & Keane, A. (2007). Multi-fidelity optimization via surrogate modelling. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463(2088). 3251–3269. https://doi.org/10.1098/rspa.2007.1900 is shown by Fig. 1.

This result looks amazing! However, when we only slightly change the lower fidelity function (and thereby the relation between the levels), this result completely breaks down as is shown by Fig. 2

The accompanying change in functions is given by the following Python snippet:

1
2
3
4
5
6
7
8
9

def HF_function(x): # high-fidelity function
    return ((6*x - 2)**2)*np.sin((12*x - 4))  

def LF_function(x): # low-fidelity function
    return 0.5 * HF_function(x) + (x-0.5) * 10. - 5
    # or:  0.5 * ((6*x - 2)**2)*np.sin((12*x - 4.0)) + (x - 0.5) * 10. - 5

def LF_function_adapted(x):
    return 0.5 * ((6*x - 2)**2)*np.sin((12*x - 4.9)) + (x - 0.5) * 10. - 5

So, with respect to the high-fidelity, the adapted low-fidelity function only experiences a small shift in the sinusoïdal term and a resulting shift in its minimum. This should not be a problem for a robust method, right? However, if we now use the same MFK methodology, we get quite bad results as shown in Fig. 3!

Now, normally one might say the lower fidelity does not represent the high fidelity truth enough, in which case we could opt for a level with a higher fidelity than the current low fidelity. So, let’s pick a medium-fidelity, in this case a simple average between the medium- and high-fidelity. The corresponding formula is given in the Python snippet below.

1
2

def MF_function(x):
    return LF_function_adapted(x) + (HF_function(x) - LF_function_adapted(x))/2

The result is shown in Fig. 4. After adding the medium-fidelity, the MFK (left) surrogate performs as bad as seen before. What is worse: it believes it is quite accurate too! Clearly, the underlying relationships between the fidelity levels cannot be captured by the MFK method. This is where the proposed method comes in (right side of the image)! It is perfectly able to capture and leverage the relations between the fidelities to exactly re-create the high-fidelity truth. What is more: it can do so with less high-fidelity data points (cheaper!), given proper validity of the assumption the method makes on these fidelity relations.

Sure enough, we trade one required data property for the other, but this will extend the applicable range of multi-fidelity modelling. Moreover, wether or not the proposed method’s assumption on the fidelity relations is valid can be cheaply and reliably tested, as we will see later.

Experimental setup

In the experimental setup to show the efficacy of the proposed method, I use two types of test cases:

• One using synthetic optimisation (minimisation) functions, altered to include noise and transformed to multi-fidelity variants using various relationships between the fidelity levels. In total, around a 1000 variations were used of the analytical functions shown in Fig. 5.

  

• One where the objective function is retrieved using actual CFD simulations, adapted to my needs and limitations. The solver is in development and limited to 2D scenarios. Shape variations are produced using a one-to-one parameterisation using Non-Uniform Rational B-Splines (NURBS). Fidelity levels are realised by using various levels of grid resolution for simulations.

For the complete experimental setup of the experiments using the synthetic objective functions I refer to my thesis document.

In this post I will focus on the workflow of the approach using CFD simulations. Below flowchart shows the aspects involved with creating the initial surrogate model (based on the initial DoE, without entering the optimisation loop). The further sections corresponding to the experimental setup will follow this flowchart.

Lifeboat case: Problem formulation

The optimisation case mimics the shape optimisation of a free-fall lifeboat: a rescue vessel that drops off an offshore drilling platform into the sea in case of emergency. We only ‘mimic’ a complete optimisation because the used CFD solver, although specialised for these types of scenario, is under development and restricted to 2D problems. The optimisation goal is to minimise the maximum (de-)acceleration the lifeboat is subjected to when impacting with a flat body of water, to guarantee the safety of its passengers. This goal is achieved by changing the bottom shape of the lifeboat.

The main dimensions, dropping height (and thereby effective vessel speed at the moment of impact), and weight of this lifeboat are taken from the recordholding FF1200 lifeboat of Palfinger Marine (the lifeboat seen at the start of this page).

With the goal of testing the proposed method, two variations of the same scenario are used:

1. The lifeboat falls perfectly vertically onto the water.
2. The lifeboat ludicrously falls horizontally onto the water.

Setting up more scenarios was not feasible given the time constraints of the thesis.

The optimisation objective is to minimize the maximum acceleration of the wedge, taking into account the safety of its hypothetical passengers: \begin{equation}\label{eqn:objective_function} y = \min\limits_{x\in X}\left(\max(a_x)\right) \end{equation} Here, $a_x$ denotes a time trace of accelerations, at location $x$ within the full domain $X$.

Shape parameterisation

The shape parameterisation of the lifeboat is done using Non-Uniform Rational B-Splines (NURBS, see Citation: Piegl & Tiller (1997) Piegl, L. & Tiller, W. (1997). The NURBS Book 2ed . Springer . ) using 2 parameters. The limitation of 2 parameters is due to the 2D nature of the CFD solver and the requirement of one-to-one parameter-shape mappings. Fig. 6 shows the used range of possible shape variations of the 2D lifeboat bottom. To create the B-Splines I used the geomdl python package.

In the creation of these shapes, the numerical stability and feasibility of the CFD solver imposes some constraints:

• The keel cannot be too flat, to avoid flat-plate impacts that impose extreme high pressures. Although this is already something occuring in the real-world, the CFD solver has the additional factor of a discretised grid with limited (vertical) resolutions able to capture the water displacements.
• The keel cannot be too sharp, due to problems otherwise encountered during the interface reconstruction calculations.
• The bilge cannot transition to be horizontal, again to avoid flat-plate impacts. The relevant outermost angles have therefore been reduced by 10 degrees.

In-grid body reconstruction

The parameterised lifeboat shape must be mapped onto the numerical grid of the CFD solver, to enable actually simulating the impact of the falling lifeboat. The grid definition itself should be symmetric over the lifeboat and in a domain large enough to not inflict boundary effects, next to solver-specific requirements.

So, the continuous NURBS shape must be transformed into its counterpart in a discretised numerical grid. This is done through the following steps:

From a NURBS shape parameterisation a Signed Distance Function (SDF) is created. This function provides the shortest distances to the body boundary with negative values defined to be inside the body and positive values outside the body. My implementation uses the compute_sdf function of the glumpy library.

This SDF function is then used to project a body-fraction field upon the simulation grid corresponding to the NURBS shape. Lastly, based on this fraction field the body is reconstructed in the grid using the Piecewise Linear Interface Calculation (PLIC) technique.

Fidelity level definition

A necessary step to transform the process of the flowchart into its multi-fidelity counterpart is to determine how the fidelity levels are determined. By definition of the method proposed during my thesis, which assumes that the grid-convergence provides information we can use, the construction of these fidelity levels is done by varying the resolution of the numerical grid. The proposed method uses three such fidelity levels:

• The lowest grid resolution (low-fidelity) cannot be rougher than required for achieving numerical stability during the CFD simulation
• The highest grid resolution (high-fidelity) should be such that the result is converged and can be considered the ’truth’. In this choice, computational budget can also be a significant limiting factor.
• The in-between grid resolution (medium-fidelity) should be chosen such that there is some convergence of the objective value present over the fidelity levels.

Like for any multi-fidelity experimental setup, this step is crucial and turned out to be challenging.

Experimental method

Simulation starting conditions

The surrogate model can be initialised in different ways depending on the surrogate modeling method used.

• The MFK method, due to limitations in solving the underlying GLS problem, needs to start with a minimum of 3 high-fidelity samples in the initial DoE.
• The proposed surrogate modelling method on the other hand can create a surrogate model with as little as 1 high-fidelity sample. To more fairly compare the proposed method with MFK but also consider this potential competitive capability, the proposed method is started using both 1 and 3 high-fidelity samples in the initial DoE.

Surrogate performance measurement

The surrogate’s accuracy is measured by comparing the predicted value with the real sampled results by using the Root Mean Squared Error (RMSE) as done in Citation: Toal (2015) Toal, D. (2015). Some considerations regarding the use of multi-fidelity Kriging in the construction of surrogate models. Structural and Multidisciplinary Optimization, 51(6). 1223–1245. https://doi.org/10.1007/s00158-014-1209-5 :

\begin{equation} \mathrm{RMSE}=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left(y_{e_{i}}-y_{c_{i}}\right)^{2}} \end{equation}

However, I prefer the normalised RMSE% for better comparison between different optimisation cases:

\begin{equation}\label{eqn:NRMSE} \mathrm{NRMSE}=\frac{\mathrm{RMSE}}{\max(y_e)-\min(y_e)} \end{equation}

RMSEs are calculated using the high-fidelity `truth’ sampled on the full initial DoE $y_e$, where subscript $i$ indexes such that high-fidelity samples used in the surrogate model are excluded to prevent a bias for surrogate models that use a lot of high-fidelity samples.

Results

Synthetic cases & structured experiments

In short, in the synthetic toy-case experiments the proposed method consistently outperforms the reference SBGO methodology of Citation: Meliani, Bartoli & al. (2019) Meliani, M., Bartoli, N., Lefebvre, T., Bouhlel, M., Martins, J. & Morlier, J. (2019). Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design. AIAA Aviation 2019 Forum. https://doi.org/10.2514/6.2019-3236 (multi-fidelity optimisation based on MFK) in terms of optimisation cost, NRMSE of the surrogate and closeness of the found design to the optimum, even when confronted with severe noise or when deteriorating the validity of the assumption upon which the proposed method relies. This is because the proposed method, as a surrogate, turns out to be more accurate for the given circumstances than the traditional multi-fidelity MFK surrogate of Citation: Gratiet & Garnier (2014) Gratiet, L. & Garnier, J. (2014). Recursive co-kriging model for design of computer experiments with multiple levels of fidelity. International Journal for Uncertainty Quantification, 4(5). 365–386. https://doi.org/10.1615/int.j.uncertaintyquantification.2014006914 , the surrogate at the basis of the SBGO method of Citation: Meliani, Bartoli & al. (2019) Meliani, M., Bartoli, N., Lefebvre, T., Bouhlel, M., Martins, J. & Morlier, J. (2019). Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design. AIAA Aviation 2019 Forum. https://doi.org/10.2514/6.2019-3236 .

These findings confirm the potential of the proposed method under a more varied and complex set of scenarios than was seen before.

Lifeboat cases

In this section, we will be looking at two things:

• The competitiveness of optimisation of the proposed surrogate modelling method compared to optimisation using the MFK surrogate.
• The quality of the optimisation result and its found designs

In both cases, the MFK surrogate knows a large NMRSE (also visually), caused by a wrongly inferred difference trend model. Around the optimum, there is no to little convergence of the objective value (meaning a difference of $\approx 0$), while elsewhere there is convergence present (meaning differences $>0$) with differences that seem to stabilise when moving further from the optimum design. With three starting samples in the initial DoE, the MFK methodology therefore infers a parabolic difference model, and faulty maps this to the rest of the domain.

The proposed surrogate on the other hand does not suffer from this problem. This is partly because the method recognises undesirable conditions and defaults to the medium-fidelity, which in this case is - especially around the optimum design where little to no convergence is present - a stable and accurate alternative.

Best designs per optimisation case and their hydrodynamic explanation

The following shows renders of the simulations that correspond to the best designs found during the vertically and horizontally falling lifeboat cases, respectively.

Note that for the ‘vertically falling’ case, the different orientation is encoded by a higher specific mass of the body in the simulation. We therefore do not actually see a more vertically oriented structure as compared to the horizontally falling case.

This subsection will further provide a qualitative hydrodynamic explanation of the optimisation results that were provided in this section. While doing this, we consider the following points made by Citation: Eijk & Wellens (2023) Eijk, M. & Wellens, P. (2023). Experimental, numerical and analytical comparison of a 2D buoyant falling and emerging wedge. Journal of Fluids Mechanics. :

1. The buoyancy force on the wedge is insignificant during slamming.
2. The added mass increases quadratically with the wetted length between contact points of the free surface with the wedge. Added mass is the apparent increase in the mass of a body moving through a fluid due to the inertia of the fluid that must also be displaced as the body moves.
3. The jets formed in the stage before detachment from the chine of the wedge lead to a strong drop in acting force and a decrease in added mass.
4. After the separation of the free surface from the chine, the force drops even more rapidly.

The last two points are related to the high-pressure front caused by the slamming escaping from between the body of water and the rigid body, dispersing energy into the jet and air.

Roughly, the inertial effects from accelerating the surrounding water to make place for the wedge constitute the slamming force. For example, this means that when the wedge would have an almost flat bottom, water needs to be accelerated almost perpendicular to the direction of movement, instigating high corresponding pressures and thus a high slamming force. These forces lead to a deceleration of the wedge. However, this means the water we were already moving and are effectively dragging along with us, the added mass, needs to be decelerated too. Therefore the force caused by the added mass opposes the slamming forces. A higher added mass is therefore beneficial to lowering the body’s acceleration at moments when high slamming forces occur.

With this in mind, let us look at some results.

Falling wedge: best horizontally falling design

Fig. 10 shows us that the found design stabilizes the maximum force during the largest part of the impact. Force peaks, and thereby peaks in the maximum acceleration, are prevented. There cannot be gained (much) more by reducing these peaks, which is a direct indication of a good design.

The design itself is explained by the low specific mass. The buildup of added mass earlier during the impact effectively reduces the force of the more flat-plate-like impact that occurs at a later stage, where we note that the added mass is of a relatively high influence to the the low specific mass of the boat.

Falling wedge: best vertically falling design

Take a look at Fig. 9. Interestingly enough, one would expect the optimal force curve to be as flat as possible, like in Fig. 10. Instead, there exist two local maxima that are about equal in magnitude.

The first maximum is caused by the initial impact at the keel. Then, upon closer inspection of the in-grid simulation output (you could closely look at the above video here!), we observe an initial jet caused by the impact of the keel. This jet separates and thereby lowers the impact forces, recalling points 3 and 4 made by Citation: Eijk & Wellens (2023) Eijk, M. & Wellens, P. (2023). Experimental, numerical and analytical comparison of a 2D buoyant falling and emerging wedge. Journal of Fluids Mechanics. .

Then, the second force maximum emerges when the body joins again with the new wavefront.

This design is in terms of the involved hydrodynamics quite non-trivial, and shows the power of surrogate modeling!

Bibliography