# Efficient multi-objective constrained Bayesian optimization of bridge girder

Heine Havneraas Røstum, Joseph Morlier, Sebastien Gros, Ketil Aas-Jakobsen

► To cite this version:

Heine Havneraas Røstum, Joseph Morlier, Sebastien Gros, Ketil Aas-Jakobsen. Efficient multi-objective constrained Bayesian optimization of bridge girder. Structural and Multidisciplinary Optimization, 2025, 68 (10), pp.214. 10.1007/s00158-025-04157-x. hal-05489211

HAL Id: hal-05489211

https://ut3-toulouseinp.hal.science/hal-05489211v1

Submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Efficient multi-objective constrained Bayesian optimization of bridge girder

Heine Havneraas Røstum<fnref n="1" />,<fnref n="2" /> · Joseph Morlier<fnref n="3" /> · Sebastien Gros<fnref n="1" /> · Ketil Aas-Jakobsen<fnref n="2" />

Received: 22 May 2025 / Revised: 11 September 2025 / Accepted: 15 September 2025 / Published online: 4 October 2025
© The Author(s) 2025

## Abstract

The buildings and construction sector is a significant source of greenhouse gas emissions, with cement production alone contributing 7 % of global emissions and the industry as a whole accounting for approximately 37 %. Reducing emissions by optimizing structural design can achieve significant global benefits. This article introduces an efficient multi-objective constrained Bayesian optimization approach to address this challenge. Rather than attempting to determine the full set of non-dominated solutions with arbitrary trade-offs, the approach searches for a solution matching a specified trade-off. Structural design is typically conducted using computationally expensive finite element simulations, whereas Bayesian optimization offers an efficient approach for optimizing problems that involve such high-cost simulations. The proposed method integrates proper orthogonal decomposition for dimensionality reduction of simulation results with Kriging partial least squares to enhance efficiency. Constrained expected improvement is used as an acquisition function for Bayesian optimization. The approach is demonstrated through a case study of a two-lane, three-span post-tensioned concrete bridge girder, incorporating fifteen design variables and nine constraints. A comparison with conventional design methods demonstrates the potential of this optimization approach to achieve substantial cost reductions, with savings of approximately 10% to 15% in financial costs and about 20% in environmental costs for the case study, while ensuring structural integrity.

**Replication of Results** Results from the finite element analyses can be made available upon request by e-mailing the corresponding author. The methods used for the optimization and creation of surrogate models were performed as described in the article.

**Keywords** Gaussian process regression · Multi-objective optimization · Bayesian optimization · Machine learning · Bridge engineering

## 1 Introduction

Traditionally, bridge engineers have focused on minimizing financial costs while ensuring compliance with all structural and geometric requirements. However, in recent years, the environmental impact of bridge construction has gained increasing attention. The buildings and construction sector

is a major source of greenhouse gas emissions, with cement production alone accounting for 7 % of global emissions, while the sector as a whole contributes approximately 37 % of global emissions (United Nations Environment Programme 2023). Reducing the emissions associated with constructions may therefore have a significant global impact. The problem is that the traditional bridge design process is heavily based on experience and manual calculations performed by highly skilled professionals. Given the numerous factors involved in bridge design, even before considering greenhouse gas emissions, finding the optimal balance between these factors has become an increasingly challenging task for human experts. As a result, many bridge designs today may be considered suboptimal.

Numerical optimization is well established (Nocedal and Wright 1999), but it is often difficult to apply it directly to optimization of engineering design problems. It requires

Responsible Editor: Makoto Ohsaki.

<figure><img src="image_2.png" /></figure>

Heine Havneraas Røstum
hrh@aaj.no

<footnote n="1">¹ Engineering Cybernetics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway</footnote>

<footnote n="2">² Dr Ing A Aas-Jakobsen AS, Oslo, Norway</footnote>

<footnote n="3">³ ICA, Université de Toulouse, ISAE-SUPAERO, Toulouse, France</footnote>
solving the system for multiple parameter values, which in engineering design consists of expensive data generation methods like numerical simulations, where gradient information is usually not directly available. This further adds complexity, making the process impractical due to excessive computational time and memory requirements. Surrogate-based optimization addresses these challenges by employing approximations of the underlying functions, known as surrogates, which are much faster to evaluate than the underlying functions during the numerical optimization process (Martins and Ning 2021). This approach has been shown to greatly reduce computational time in bridge design (Penadés-Plà et al. 2019). A similar approach is Bayesian optimization, which builds probabilistic surrogates, where the probabilistic element is leveraged to accelerate the optimization process (Frazier 2018). The probabilistic surrogate is typically constructed using Gaussian Process Regression (GPR). Bayesian optimization has been successfully applied to various engineering design problems; materials (Frazier and Wang 2016), civil structures (Mathern et al. 2021), aerospace (Tfaily et al. 2024; Saves et al. 2024), and bridge design (Røstum et al. 2025). It ensures fast convergence in the optimization process for problems that typically involve less than 20 decision variables (Frazier 2018).

A challenge in applying GPR is its poor scalability with dataset size, as the training time increases cubically with the number of data points. For high-dimensional simulations, such as those obtained from *Finite Element (FE)* analyses, this can make training computationally intractable if attempting to approximate the full high-dimensional output using a single surrogate model. To mitigate this issue, *Reduced Order Modeling (ROM)* has been widely adopted (Lucia et al. 2004; Peherstorfer and Willcox 2015; Chiplunkar et al. 2018; Drakoulas et al. 2024). ROM replaces the full-order model with a lower-dimensional approximation, significantly reducing both memory requirements and computational costs. A widely used technique for ROM is the *Proper Orthogonal Decomposition (POD)* (Liang et al. 2002; Amsallem and Farhat 2012; Paul-Dubois-Taine and Amsallem 2015; Guo and Hesthaven 2018; Zhao 2021; Rocha et al. 2023). Guo and Hesthaven (2018) used POD with GPR to predict unobserved design variable combinations of the full solution field of non-linear FE analyses. Their approach showed good performance in both accuracy and computational efficiency. A challenge with surrogate models is that the number of design variables, or input parameters, affects its accuracy and training time. To address this, Bouhlel et al. (2016) proposed a method called *Kriging Partial Least Squares (KPLS)* where *Partial Least Squares (PLS)* is applied for the construction of the GPR. Their results showed that with the same number of training points, they could significantly reduce training time while achieving similar or improved accuracy compared to

standard GPR when applied to problems with a moderate to high number of design variables.

Using numerical optimization in the design of new structures can be an effective way to reduce global emissions. However, for this approach to be practical, it must be efficient, both in terms of training time and the number of iterations required for the optimizer to converge. Reducing the number of iterations reduces the number of potentially time-consuming numerical simulations, significantly impacting the overall duration of the optimization process. To further enhance the appeal of this approach, the optimization process should consider both the environmental impact and the monetary cost of the structure. By striking a balance between these objectives, it is possible to achieve trade-offs that reduce both, offering improvements over traditional design methods also in the monetary sense. To achieve this, this article presents a method for multi-objective constrained Bayesian optimization using a predefined trade-off function.

Early extensions of Bayesian optimization to multi-objective problems were based on the Efficient Global Optimization (EGO) framework (Jones et al. 1998). One of the first was *ParEGO* (Knowles 2005), which employs randomly generated weight vectors in a Tchebycheff scalarization to enable the use of expected improvement within a single-objective EGO loop. Simultaneously, Jeong and Obayashi (2005) applied EGO directly to multi-objective aerodynamic design by formulating expected improvements for each objective and identifying non-dominated solutions. In recent years, the *Expected Hypervolume Improvement (EHVI)* criterion has become increasingly popular for guiding Bayesian multi-objective optimization (Yang et al. 2019; Mathern et al. 2021).

However, while these approaches are effective in exploring the Pareto front, they do not directly address the practical scenario in which only one final design will ultimately be implemented. To this end, our method targets a single solution aligned with a predefined trade-off, thereby requiring fewer iterations and number of numerical simulations than traditional multi-objective optimization. The approach is conceptually similar to the work of Jeong and Obayashi (2005), who selected three representative points on the Pareto front approximation based on the objectives' expected improvements (the extremities and the midpoint). We extend this idea by incorporating constrained expected improvement and selecting the point that maximizes the aggregated improvement over the current best feasible design. Furthermore, the proposed method utilizes the well-known POD to reduce the high-dimensional space from numerical simulations. This is combined with KPLS to obtain a computationally efficient approach to optimize engineering designs based on simulations.

<figure><img src="image_3.png" /></figure>
The approach is applied to a case study of a two-lane, three-span post-tensioned concrete bridge girder. This case study considers fifteen design variables and nine high-dimensional constraints derived from simulations, making it a realistic and relevant example for demonstrating the method's applicability to practical engineering problems. The objectives are to minimize the monetary and environmental cost of the structure, while satisfying all constraints. To assess the effectiveness of using a predefined trade-off function for the objectives, as opposed to the traditional arbitrary trade-off, the approach is compared to the constrained multi-objective Bayesian optimization scheme constrained *expected hypervolume improvement* (Daulton et al. 2020). Results are benchmarked against a conventional bridge girder design approach, which is based on the authors' experience in the bridge engineering industry. This comparison highlights the potential benefits of applying the constrained multi-objective optimization approach, demonstrating its ability to improve both efficiency and performance in bridge design.

## 2 Methodology

This section presents the constrained multi-objective optimization methodology which was applied to optimize the design of the bridge girder in the case study. Though the method is here applied to a bridge girder, the methodology is applicable to any problem that involves;

* *Expensive Data Generation* (e.g., simulation results) The high computational cost of each simulation limits the available data. Reducing the number of required simulations is crucial due to time and memory constraints.

* *Small to Moderate Number of Decision Variables* Bayesian optimization struggles to scale with a large number of decision variables due to the *curse of dimensionality*. A higher-dimensional space demands more data, which GPR handles inefficiently. KPLS improves accuracy in moderate-dimensional problems with limited data, enhancing the effectiveness of the proposed approach.

These conditions are common in many engineering design problems. The section starts with a description of the initial sampling using *Latin Hypercube Sampling*, followed by a description of POD, which is used to create a low-dimensional approximation of the high-dimensional, also known as functional, simulation outputs from the initial observations. To predict unobserved variable combinations, GPR is used to interpolate the low-dimensional POD bases. Additionally, PLS is employed to reduce the number of GPR

kernel parameters, with this reduction of kernel parameters referred to as KPLS. As GPR is a probabilistic regressor, the probabilistic element of it may be exploited in a Bayesian optimization scheme using a Bayesian acquisition function. If the constraints are also approximated using GPR, they can be included to formulate a constrained Bayesian optimization scheme, which is done in this article. Since multiple objective functions are involved, each one is defined using a constrained Bayesian acquisition function and incorporated into a *multi-objective optimization* (MOO) framework, which will be detailed at the end of this section. Rather than approximating the entire set of non-dominated solutions, this MOO approach focuses on efficiently identifying a single solution that best satisfies a predefined trade-off between the objectives, aiming to do so in as few iterations as possible.

The MOO process is illustrated in Fig. 1. The figure illustrates that the POD and KPLS can be used together or separately. The POD is applied when one is approximating the full solution field of a large dimensional output space, which could be the results from *Finite Element Method* (FEM) or *Computational Fluid Dynamics* (CFD) simulations. If a moderate to high number of design variables are considered, KPLS should be used as a surrogate model instead of ordinary GPR. A multi-objective optimization is subsequently conducted, using a preferred trade-off between the objectives, which will be further described at the end of this section. If the exit criterion is satisfied, the best observed solution according to the specified trade-off is returned; otherwise, the next candidate solution is incorporated as a new sample and the process is repeated.

In this article, we use bold uppercase letters to denote matrices, bold lowercase letters to denote vectors, and plain lowercase letters to denote a scalar or variable.

### 2.1 Initial sampling

The initial sampling method, also sometimes known as Design of Experiments (DoE), is the process of determining the initial combinations of design variables, $\mathbf{X} = \{\mathbf{x}_j\}_{j=1}^{N_s}$, where $\mathbf{x}_j \in \mathbb{R}^d$, from which to acquire observations. An observation could be the results of FEM analyses or from other calculations using a combination of design variables $\mathbf{x}_j$. $N_s$ is the number of combinations of design variables to sample, or the sampling budget, and $d$ is the number of design variables. From the initial combination of design variables, an observation set $\mathcal{D} = \{(\mathbf{x}_j, \mathbf{y}(\mathbf{x}_j))\}_{j=1}^{N_s}$ is acquired. A popular sampling strategy is the *Latin Hypercube Sampling* (LHS) strategy (Negrin et al. 2023). LHS was first described by McKay et al. (1979). The LHS strategy implemented in this article is to generate a LHS with a normalized distribution in the direction of each variable. The
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Illustration of the proposed efficient multi-objective Bayesian optimization process</figcaption></figure>

technique is similar to that described in the book by Martins and Ning (2021) for LHS with normalized distribution. An

illustration of a two-variable LHS with distributions projected to the axes is shown in Fig. 2.

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 2 Illustration of two-variable Latin hypercube sampling with distribution projected to axes</figcaption></figure>

A normalized distribution is chosen to ensure that not multiple design variables are in the variable bounds' extremities for the same variable combination, as this may result in analytical problems during the FEM analyses. This depends on the type of analysis and software that is used for the simulations. However, it should be noted that in high-dimensional spaces, much of the volume of the design domain lies near the boundaries. Under-sampling these regions may lead to over-exploration by the acquisition function in later iterations due to the high associated uncertainty. A uniform distribution may typically be more applicable to most situations and ensure that more information is obtained from the initial sampling.

## 2.2 Proper orthogonal decomposition (POD)

Approximating the entire output space from FEM simulations, such as, e.g., displacements for the entire structure, is typically computationally demanding using a single GPR due to the high dimensionality of the results. POD is a data analysis method with the purpose of obtaining a low-dimensional approximation of a high-dimensional process (Chatterjee 2000). If each result for the entire output space for different combinations of design variables is stored in a vector $\mathbf{y}(\mathbf{x}) \in \mathbb{R}^{N_h}$, where $N_h$ is the number of datapoints of the high-dimensional output. The results from the initial sampling can be stored in a matrix $\mathbf{Y} \in \mathbb{R}^{N_h \times N_s}$;

$$ \mathbf{Y} = \begin{bmatrix} y_{1,1} & \cdots & y_{1,N_s} \\ \vdots & \ddots & \vdots \\ y_{N_h,1} & \cdots & y_{N_h,N_s} \end{bmatrix} \quad (1) $$

y and $\mathbf{Y}$ are here used as general terms to denote observed responses, which could be constraints or objectives. In the context of POD, the $\mathbf{Y}$ matrix is often known as a snapshot matrix, where the number of simulations $N_s$ is the number

of snapshots. We assume that all simulations share the same underlying mesh. If they do not, their outputs must first be interpolated or projected onto a common reference grid to maintain consistency. POD is commonly achieved using the Singular Value Decomposition (SVD) (Guo and Hesthaven 2018). Applying SVD to the matrix $\mathbf{Y}$ yields:

$$ \mathbf{Y} = \mathbf{U} \Sigma \mathbf{V}^{\top} \quad (2) $$

Here $\mathbf{U} \in \mathbb{R}^{N_h \times N_s}$ and $\mathbf{V} \in \mathbb{R}^{N_s \times N_h}$ are orthogonal matrices, and $\Sigma = \text{diag}\{\sigma_1, \sigma_2, \dots, \sigma_{N_s}\}$ are the singular values (eigenvalues) in descending order. A low-dimensional truncated POD basis, $\mathbf{L} \in \mathbb{R}^{N_h \times N_{POD}}$, which is derived by preserving the first $N_{POD}$ modes of $\mathbf{U}$ that corresponds to the largest singular values, can be used to calculate the projection coefficients $\mathbf{A} \in \mathbb{R}^{N_{POD} \times N_s}$, also sometimes known as the amplitude matrix (Buljak and Maier 2011; Gooijer et al. 2021), using the relationship:

$$ \mathbf{A} = \mathbf{L}^{\top} \cdot \mathbf{Y} \quad (3) $$

A common method of selecting the number of POD modes to include is by evaluating the projection error $\epsilon_{POD}$ of including the first $N_{POD}$ left singular vectors of $\mathbf{Y}$ (Eckart and Young 1936; Guo and Hesthaven 2018; Brunton and Kutz 2022). Thus $N_{POD}$ can be defined as the minimum integer such that:

$$ \frac{\sum_{i=1}^{N_{POD}} \sigma_i^2}{\sum_{i=1}^{N_s} \sigma_i^2} > 1 - \epsilon_{POD}^2 \quad (4) $$

In this article, an $\epsilon_{POD}^2$ of 0.01 is used for selecting the number of POD modes, $N_{POD}$. Due to the orthogonality of $\mathbf{U}$, an approximation of the snapshot matrix $\mathbf{Y}$ can be made from the low-dimensional POD basis $\mathbf{L}$ and the projection coefficients $\mathbf{A}$ as:

$$ \mathbf{Y} \approx \mathbf{L} \cdot \mathbf{A} \quad (5) $$

Using the training set $\mathcal{D} = \{(\mathbf{X}, \mathbf{A}_i)\}_{i=1}^{N_{POD}}$, one regression model $\hat{\mathbf{A}}_i(\mathbf{x})$ is created to approximate each of the $N_{POD}$ projection coefficients in $\mathbf{A} \in \mathbb{R}^{N_{POD} \times N_s}$, resulting in the predicted projection coefficients:

$$ \hat{\mathbf{A}}(\mathbf{x}) = [\hat{\mathbf{A}}_1(\mathbf{x}), \hat{\mathbf{A}}_2(\mathbf{x}), \dots, \hat{\mathbf{A}}_{N_{POD}}(\mathbf{x})]^{\top} \quad (6) $$

Approximations of the full set solution $\hat{\mathbf{Y}}$ can be retrieved for new variable combinations $\mathbf{x}$ using the equation:

$$ \hat{\mathbf{Y}}(\mathbf{x}) = \mathbf{L} \cdot \hat{\mathbf{A}}(\mathbf{x}) \quad (7) $$

In this article, Gaussian process regression (GPR) is applied as the regression model, which will be further described in the following section.
## 2.3 Gaussian process regression (GPR)

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution (Williams and Rasmussen 2006). A real process, $h(\mathbf{x})$ can be formulated completely as a Gaussian process by its mean function, $\mu(\mathbf{x})$, and covariance function $k(\mathbf{x}, \mathbf{x}')$, with the Gaussian process formulated as:

$$h(\mathbf{x}) \sim GP(\mu(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')) \quad (8)$$

The mean is often assumed to be zero for notational simplicity, as will be done in this article, and this can be approximated by mean centering the distribution. The covariance function, $k(\mathbf{x}, \mathbf{x}')$, typically determines the GP's smoothness, and characteristic length-scale. In this article, the Matérn 5/2 kernel is used, which is formulated as (Williams and Rasmussen 2006):

$$k(\mathbf{x}, \mathbf{x}') = \sigma_f^2 \prod_{i=1}^{d} \left(1 + \sqrt{5}r_i + \frac{5}{3}r_i^2\right) \exp\left(-\sqrt{5}r_i\right) \quad (9)$$

where $d$ is the number of design variables, $\sigma_f^2$ is the signal variance, and $r_i$ is a distance metric for the $i$th design variable component. In this article, a distance metric with an anisotropic length scale $\theta_i$ is used:

$$r_i = \theta_i |x_i - x_i'| \quad (10)$$

When a set of $N_s$ observations are available, the observations $\mathcal{D} = \{(\mathbf{x}_j, y_j)\}_{j=1}^{N_s}$ can be used to in order to infer a joint Gaussian posterior distribution for unobserved data, $\mathbf{y}_*(\mathbf{X}_*)$:

$$\mathbf{y}_*|\mathbf{X}_*, \mathcal{D} \sim \mathcal{N}(\hat{\mu}(\mathbf{X}_*), \hat{\sigma}^2(\mathbf{X}_*)) \quad (11)$$

where the posterior mean, $\hat{\mu}(\mathbf{X}_*)$, and variance, $\hat{\sigma}^2(\mathbf{X}_*)$, are calculated as:

$$\hat{\mu}(\mathbf{X}_*) = K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I}]^{-1} \mathbf{y} \quad (12a)$$

$$\hat{\sigma}^2(\mathbf{X}_*) = K(\mathbf{X}_*, \mathbf{X}_*) - K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I}]^{-1} K(\mathbf{X}, \mathbf{X}_*) \quad (12b)$$

For compactness, $K(\mathbf{X}, \mathbf{X})$ will be abbreviated as $\mathbf{K}_{\text{ff}}$, and $K(\mathbf{X}_*, \mathbf{X}_*)$ as $\mathbf{K}_{**}$, and so forth. The $\sigma_n^2 \mathbf{I}$ term is the deviation between the observations and the regressor. $\sigma_n^2$ is the deviation variance, which is often referred to as the noise variance, and $\mathbf{I} \in \mathbb{R}^{N_s \times N_s}$ is an identity matrix. Though the results from the FEM simulations are noise free, $\sigma_n^2$ is included as a regularization term.

The marginal likelihood is commonly used to determine the model hyperparameters $\beta = \{\sigma_f^2, \theta_1, \dots, \theta_d, \sigma_n^2\}$, which

consist of the noise variance, signal variance, and the parameters associated with the covariance function. The logarithmic marginal likelihood is formulated as:

$$\log p(\mathbf{y} | \mathbf{X}) = -\frac{1}{2} \log \det(\mathbf{K}_{\text{ff}} + \sigma_n^2 \mathbf{I}) \\ -\frac{1}{2} \mathbf{y}^{\top} (\mathbf{K}_{\text{ff}} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{y} - \frac{n}{2} \log(2\pi) \quad (13)$$

The optimal hyperparameters $\beta_{\text{opt}}$ are then typically obtained by the maximization problem:

$$\beta_{\text{opt}} = \arg \max_{\beta} \log p(\mathbf{y} | \mathbf{X}) \quad (14)$$

Which is also used in this article. The first term in Eq. (13) penalizes high complexity, in the sense that large values of $\mathbf{K}$ is discouraged in the covariance function, as it results in non-smooth behavior. The second term ensures a good data fit, and the last term is a normalization constant. Optimizing the hyperparameters using the marginal likelihood thus ensures a good fit for the observed values, and reduces the risk of overfitting by penalizing complex models.

## 2.4 Kriging partial least squares (KPLS)

Bouhlel et al. (2016) showed that weights derived from Partial Least Squares (PLS) could be used to efficiently reduce the number of kernel hyperparameters associated with the design variables from $\{\theta_i\}_{i=1}^d$ to $\{\theta_l\}_{l=1}^h$, where $h < d$. They called this method KPLS and showed that the method resulted in a significant reduction in the CPU training time, and that the method in general provided comparable or better predictive metrics for problems involving a medium to high number of hyperparameters. The method relies on PLS, which is a method for determining an optimum linear relationship between inputs $\mathbf{X}$ and outputs $\mathbf{y}$ (Helland 1988). The number of components, $h$, to include in the parameter reduction is chosen a priori, but a low number of weights is typically selected. By applying PLS with $h$ number of components to observed inputs $\mathbf{X} \in \mathbb{R}^{N_s \times d}$ and outputs $\mathbf{y}(\mathbf{X}) \in \mathbb{R}^{N_s}$, one obtains the PLS weights $\mathbf{W} \in \mathbb{R}^{d \times h}$ and loadings $\mathbf{P} \in \mathbb{R}^{d \times h}$. A new projected coordinate system, $\mathbf{T} \in \mathbb{R}^{N_s \times h}$, obtained by rotating the original system, can be formulated as (Bouhlel et al. 2016):

$$\mathbf{T} = \mathbf{X}\mathbf{W}_* \quad (15)$$

where $\mathbf{W}_*$ is known as the projection or rotation matrix and can be calculated using the relationship:

$$\mathbf{W}_* = \mathbf{W}[\mathbf{P}^{\top}\mathbf{W}]^{-1} \quad (16)$$

The rotation matrix $\mathbf{W}_* = [\mathbf{w}_*^{(1)}, \dots, \mathbf{w}_*^{(h)}]$, is used to reduce the number of hyperparameters associated with the
length-scales $\theta$ in the kernel from $d$ to $h$. The KPLS version of the Matérn 5/2 kernel is given in Eq. (17).

$$k(\mathbf{x}, \mathbf{x}') = \sigma_f^2 \prod_{l=1}^{h} \prod_{i=1}^{d} \left[ 1 + \sqrt{5} r_i^{(l)} + \frac{5}{3} (r_i^{(l)})^2 \right] \exp(-\sqrt{5} r_i^{(l)}) \quad (17)$$

Here, the rotation coefficient $w_{*i}^{(l)}$ is included in the distance metric $r_i^{(l)}$ as given in Eq. (18).

$$r_i^{(l)} = \theta_i \left| w_{*i}^{(l)} (\mathbf{x}_i - \mathbf{x}'_i) \right| \quad (18)$$

## 2.5 Optimization process

A typical multi-objective optimization problem of determining the optimum design variables $\mathbf{x}^*$ subjected to inequality constraints can mathematically be formulated as:

$$\mathbf{x}^* = \arg \min_{\mathbf{x} \in \mathbf{X}} \left( f_1(\mathbf{x}), f_2(\mathbf{x}), \dots, f_{n_f}(\mathbf{x}) \right) \quad (19a)$$

$$\text{s.t} \quad g_i(\mathbf{x}) \le \lambda_i, \quad i = 1, 2, \dots, n_c \quad (19b)$$

where $f(\mathbf{x}) = (f_1(\mathbf{x}), f_2(\mathbf{x}), \dots, f_{n_f}(\mathbf{x}))$ are the objective functions, subjected to $n_c$ number of inequality constraints expressed by the $i^{th}$ constraint function as $g_i(\mathbf{x})$. The goal is often to identify the set of feasible solutions for which no other feasible solution dominates them. A solution $\mathbf{x}^* \in \mathbf{X}$ is non-dominated if there exists no other feasible solution $\mathbf{x}' \in \mathbf{X}$ such that $f(\mathbf{x}') < f(\mathbf{x}^*)$. The corresponding set of feasible non-dominated solutions is known as the Pareto front approximation $\mathcal{P} = \{f^{(1)}, \dots, f^{(n_p)}\} \subset \mathbb{R}^{n_p}$, where $n_p$ is the number of non-dominated solutions. If the objective and constraint functions are approximated by probabilistic surrogate models, such as GPR, the probabilistic element can be exploited in a Bayesian optimization framework. In the Bayesian optimization framework, the original constrained multi-objective problem defined in Eq. (19) is reformulated using an acquisition function that guides the selection of the next design (or sampling) point. This acquisition function typically transforms the multi-objective minimization problem into a scalar maximization problem. Whether the goal is to approximate the entire Pareto front or to target a specific trade-off solution, the optimization process is generally performed one point at a time.

We first present a popular acquisition function for single-objective constrained Bayesian optimization. This is followed by extensions of it for multi-objective constrained Bayesian optimization. In the framework of multi-objective constrained Bayesian optimization, we consider two distinct methodologies:

1. **The method proposed in this work**, which targets a specific solution or trade-off among the objective functions. The goal is to identify a single desirable design as efficiently as possible, rather than deriving the entire Pareto front.

2. **The constrained expected hypervolume improvement approach**, which represents a standard, off-the-shelf multi-objective Bayesian optimization method. It is designed to explore the objective space broadly and incrementally improve the entire Pareto front. This serves as a reference for comparison.

The purpose of the proposed approach is to rapidly converge to a single balanced trade-off point by maximizing an acquisition function based on normalized objective values. Since only one structural design will ultimately be implemented, our method avoids the computational overhead of constructing a full Pareto front and instead prioritizes sample efficiency, i.e., in converging to a high-quality solution with a few number of iterations. Although scalarized single-objective approaches (e.g., a weighted sum of constrained expected improvements) could produce similar acquisition functions, they would essentially reduce to a simplified version of our framework. In contrast, cEHVI follows a fundamentally different optimization strategy aimed at broadly covering the Pareto front. This comparison highlights the distinction between traditional exploratory multi-objective methods, as exemplified by cEHVI, and the solution-focused trade-off optimization implemented in this work.

### 2.5.1 Single objective constrained Bayesian optimization

A popular acquisition function for constrained Bayesian optimization with a single objective was proposed by Gardner et al. (2014) and is called *expected constrained improvement* ($EI_c$). For $EI_c$, it is assumed that each constraint is independent and is also approximated by Gaussian processes. The $EI_c$ is the *expected improvement* (EI) of the approximation of the objective function $\hat{f}(\mathbf{x})$ relative to the best feasible objective value $f^*$ observed so far, weighted by the probability that $\mathbf{x}$ is a feasible design (Gardner et al. 2014). A feasible design is a variable combination $\mathbf{x}$ that satisfies $g_i(\mathbf{x}) \le \lambda_i, \forall i \in \{1, \dots, n_c\}$. The probability that the predictions of each constraint $\hat{g}_i(\mathbf{x})$ is lower than the limits $\lambda_i$, which is known as the probability of feasibility (PF) for design $\mathbf{x}$.

$$EI_c(\mathbf{x}) = EI(\mathbf{x})PF(\mathbf{x}) \quad (20)$$
The equations for calculating the components of EI<sub>c</sub>(x) are given in Eq. (20).

$$
EI(\mathbf{x}) = (f^* - \hat{\mu}_f(\mathbf{x}))\Phi\left(\frac{f^* - \hat{\mu}_f(\mathbf{x})}{\hat{\sigma}_f(\mathbf{x})}\right) + \hat{\sigma}_f(\mathbf{x})\phi\left(\frac{f^* - \hat{\mu}_f(\mathbf{x})}{\hat{\sigma}_f(\mathbf{x})}\right) \quad (21a)
$$

$$
\mathrm{PF}(\mathbf{x}) = \prod_{i=1}^{n_c} \Phi \left( \frac{\lambda_i - \hat{\mu}_{gi}(\mathbf{x})}{\hat{\sigma}_{gi}(\mathbf{x})} \right) \qquad (21\mathrm{b})
$$

Here, it is assumed that each of the constraints $\hat{g}_i(x)$ are
independent GPs, with mean $\hat{\mu}_{gi}(x)$ and variance $\hat{\sigma}_{gi}^2$. Variable
$f^*$ is the best observed objective value so far that satisfies
the constraints, $\hat{\mu}_f(x)$ and $\hat{\sigma}_f^2(x)$ are, respectively, the
predictive mean and variance of the objective function
approximation $\hat{f}(x)$. Function $\Phi$ is the cumulative density
function (CDF) of a standardized normal distribution, and
function $\phi$ is the probability density function (PDF) of a
standardized normal distribution.

2.5.2 Multi-objective constrained Bayesian optimization

2.5.3 Adopted approach (PTMOO)

Within our constrained multi-objective optimization framework, we compute an acquisition function EI<sub>c</sub><sup>(k)</sup> for each of the n<sub>f</sub> objective functions. The constrained expected improvement is evaluated relative to a specified preferred trade-off, which we define as the feasible design that yields the most desirable balance among objectives based on a normalized aggregate measure. We use the term preferred to emphasize that the optimization process is directed toward a specific solution in the objective space, rather than aiming to reconstruct the full Pareto front. Importantly, this preference is not explicitly defined by the user, e.g., via weights, but arises from a default strategy that assigns equal importance to all objectives. Specifically, we seek a preferred combination of objective function values f<sup>*</sup>(x<sub>N</sub><sup>*</sup>) = [f<sub>1</sub><sup>*</sup>(x<sub>N</sub><sup>*</sup>), f<sub>2</sub><sup>*</sup>(x<sub>N</sub><sup>*</sup>), ..., f<sub>n<sub>f</sub></sub><sup>*</sup>(x<sub>N</sub><sup>*</sup>)], where x<sub>N</sub><sup>*</sup> &isin; X denotes the preferred design at iteration N. This is defined as the feasible combination of design variables that minimizes the sum of normalized objective function values:

$$
\mathbf{x}_N^* = \arg \min_{\mathbf{x} \in \mathbf{X}} \sum_{k=1}^{n_f} \frac{f_k(\mathbf{x}) - f_k^{\min}}{f_k^{\max} - f_k^{\min}} \quad (22a)
$$

$$
\text{s.t} \quad g_i(\mathbf{x}) \le \lambda_i, \quad i = 1, 2, \dots, n_c \qquad (22b)
$$

<figure><img src="image_7.png" /><figcaption>Fig. 3 Illustration of best observed design for process with 2 objec-
tive functions</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 4 Pareto front constructed using two acquisition functions and the PTMOO candidate solution on the Pareto front using Eq. (23)</figcaption></figure>

Here, $f_k^{\min}$ and $f_k^{\max}$ are the minimum and maximum observed values of objective function $f_k$, evaluated across all sampled designs, both feasible and infeasible. Normalization ensures that all objectives are treated with equal importance when computing the aggregate. In cases where differing priorities are desired, weighted normalization can be applied to reflect the relative importance of each objective. Figure 3 illustrates the selection of the best, or preferred, feasible design in a two-objective scenario with equal weighting, as used in Eq. (22).

A Pareto front approximation is then calculated using the multi-objective genetic algorithm NSGA-II implemented in the Python package *pymoo* (Blank and Deb 2020). See the article by Deb et al. (2002) for further information on the NSGA-II algorithm. The objectives of the applied NSGA-II are to maximize the constrained acquisition function EI<sub>c</sub><sup>(k)</sup> for each of the n<sub>f</sub> number of objectives. An illustration of a Pareto front for two acquisition functions is shown with blue dots in Fig. 4. When the Pareto front has been determined, the next candidate design, x<sub>N+1</sub>, is determined as the variable combination x with the maximum sum of squared expected
relative improvement of the variable combinations in the Pareto front $X_p$, which is calculated as:

$$x_{N+1} = \arg \max_{x \in X_p} \sum_{k=1}^{n_f} \left( \frac{EI_c^{(k)}(x)}{f_k^*(x_N^*)} \right)^2 \quad (23)$$

The next candidate design, $x_{N+1}$, is chosen to maximize the aggregate relative improvement over the current preferred solution. It is important to note that $x_{N+1}$ denotes the next candidate design to be evaluated, not necessarily the new preferred design. The actual preferred design at iteration $N + 1$, denoted $x_{N+1}^*$, is determined only after evaluating the true objectives and constraints of $x_{N+1}$. We will refer to this selection rule as PTMOO for the remainder of the paper. By using relative improvement, each objective is treated with equal weight while driving the largest possible advancement from the existing preferred design. In this article, a single point is selected. Thus, one could also skip calculating the Pareto front approximation and select the next candidate design by simply maximizing Eq. (23) for the entire possible design set $X \in \mathbb{R}^d$. The two-stage approach offers three key advantages over “blind” single-objective maximization:

1. Decomposes a hard global search into tractable pieces $EI_c^{(k)}(x)$ may be highly multi-modal, so direct global optimization may get stuck in poor local maxima. By first using NSGA-II to approximate the Pareto front of the problem we efficiently identify a diverse set of high-quality candidate points.

2. Multiple candidate solutions Even though we ultimately choose one $x$ as the next candidate design in this article, the Pareto front approximation can be used to determine multiple candidate points for which to perform new analyses.

3. Insight into trade-offs at each iteration The approximate Pareto front reveals how improvements in one objective affect others, offering valuable insights into trade-offs and guiding early stopping decisions.

Figure 4 illustrates the Pareto front approximation of a process with two objective functions and the location of the next candidate design determined using Eq. (23).

Algorithm 1 summarizes the proposed efficient multi-objective constrained Bayesian optimization procedure, which incorporates a user-specified preferred trade-off for the objective functions. In our implementation, the full set of hyperparameters is trained in each iteration when approximating the constraints and objectives. We do not define a hard threshold on the preference of KPLS over GPR without hyperparameter reduction. Bouhlel et al. (2016) showed that when more than ten design variables were considered, KPLS generally matched or outperformed standard GPR. Similarly, we do not prescribe a strict cutoff for applying POD based on output dimensionality, yet we recommend its use when seeking to approximate the full solution set from FEM or CFD analyses.

<figure><img src="image_9.png" /><figcaption>Fig. 5 Illustration of a 3-span post-tensioned concrete road bridge with 2 lanes</figcaption></figure>
**Algorithm 1 Efficient Multi-Objective Constrained Bayesian Optimization Process**

1: **Input:** Bounds of design variables [X<sub>min</sub>, X<sub>max</sub>], initial sampling budget N<sub>s</sub>, constraint limits λ

2: **Output:** Preferred design x<sup>*</sup>

3: Generate initial sample variables X<sub>0</sub> using Latin Hypercube Sampling (LHS) with budget size N<sub>s</sub>

4: Evaluate initial observations: D<sub>0</sub> = {(*x*<sub>j</sub>, *f*(*x*<sub>j</sub>), *g*(*x*<sub>j</sub>))}N<sub>s</sub>j=1 = (*X*<sub>0</sub>, *Y*<sub>D<sub>0</sub></sub>)

5: Determine initial preferred design x<sub>0</sub><sup>*</sup> using Eq. (22)

6: Set converged ← **False**

7: Set iteration counter N ← 0

8: **while** converged = **False** **do**

9: &nbsp;&nbsp;&nbsp;&nbsp;**for** each objective k = 1 to n<sub>f</sub> **do**

10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** Output dimension of *f*<sub>k</sub>(*X*<sub>*N*</sub>) is large **then**

11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Apply Proper Orthogonal Decomposition (POD) to *f*<sub>*k*</sub>(*X*<sub>*N*</sub>) for dimension reduction

12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**

13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** Number of design variables is high **then**

14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Approximate *f*<sub>*k*</sub>(*X*<sub>*N*</sub>) using KPLS

15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**

16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Approximate *f*<sub>*k*</sub>(*X*<sub>*N*</sub>) using GPR

17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**

18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**

19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each constraint i = 1 to n<sub>c</sub> **do**

20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** Output dimension of *g*<sub>*i*</sub>(*X*<sub>*N*</sub>) is large **then**

21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Apply Proper Orthogonal Decomposition (POD) to *g*<sub>*i*</sub>(*X*<sub>*N*</sub>) for dimension reduction

22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**

23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** Number of design variables is high **then**

24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Approximate *g*<sub>*i*</sub>(*X*<sub>*N*</sub>) using KPLS

25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**

26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Approximate *g*<sub>*i*</sub>(*X*<sub>*N*</sub>) using GPR

27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**

28: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**

29: Compute Pareto front using Eq. (20) as acquisition functions for the n<sub>f</sub> objectives

30: Select next candidate design x<sub>*N*+1</sub> from Pareto front using Eq. (23)

31: Evaluate objectives and constraints: {f<sub>1</sub>(x<sub>*N*+1</sub>), ..., f<sub>n<sub>f</sub></sub>(x<sub>*N*+1</sub>), g<sub>1</sub>(x<sub>*N*+1</sub>), ..., g<sub>n<sub>c</sub></sub>(x<sub>*N*+1</sub>)}

32: Update observations: D<sub>*N*+1</sub> ← D<sub>*N*</sub> ∪ (x<sub>*N*+1</sub>, f(x<sub>*N*+1</sub>), g(x<sub>*N*+1</sub>))

33: Determine preferred design x<sub>*N*</sub><sup>*</sup> using Eq. (22)

34: **if** Convergence criteria met **then**

35: &nbsp;&nbsp;&nbsp;&nbsp;converged ← **True**

36: &nbsp;&nbsp;&nbsp;&nbsp;**end if**

37: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;N ← N + 1

38: &nbsp;&nbsp;&nbsp;&nbsp;**end while**

39: **Return** x<sub>*N*+1</sub><sup>*</sup>

### 2.5.4 Constrained expected hypervolume improvement (cEHVI)

To compare the results obtained using Algorithm 1, the results obtained using the constrained Expected Hypervolume Improvement (cEHVI) acquisition function are applied. The cEHVI is an extension of the cEI defined in Eq. (21) (Daulton et al. 2020). Given a finite approximation set of n<sub>P</sub> points to a Pareto front P = {f<sup>(1)</sup>, ..., f<sup>(n<sub>P</sub>)</sup>} ⊂ ℝ<sup>n<sub>f</sub></sup>, the hypervolume (HV) indicator measures the size of the dominated subspace bounded above by a reference point r ∈ ℝ<sup>n<sub>f</sub></sup>:

$$ \mathrm{HV}(\mathcal{P} \Leftrightarrow \nabla) = \lambda_{n_f} (\cup_{f \in \mathcal{P}} [f, r]), \quad (24) $$

where λ<sub>n<sub>f</sub></sub> is the Lebesgue measure. For a vector f(x) ∈ ℝ<sup>n<sub>f</sub></sup>, the hypervolume improvement (HVI) is calculated as:

$$ \mathrm{HVI}(\mathbf{x}, \mathcal{P}, r) = \mathrm{HV}(\mathcal{P} \cup \{f(\mathbf{x})\}, r) - \mathrm{HV}(\mathcal{P} \Leftrightarrow \nabla) \quad (25) $$

When the objectives f(x) are approximated by a predictive distribution, e.g., GPR, then for a given set of predictive distributions p(f|x) = N(μ(x), σ<sup>2</sup>(x)), and the Pareto front approximation P, the expected hypervolume improvement (EHVI) is defined as:

$$ \mathrm{EHVI}(\mathbf{x}, \mathcal{P}, r) = \int_{\mathbb{R}^{n_f}} \mathrm{HVI}(\mathbf{x}, \mathcal{P}, r) \cdot p(f|\mathbf{x})df(\mathbf{x}) \quad (26) $$
<table><caption>**Table 1** Variables for the MOO and associated upper and lower bounds</caption><thead><tr><th>Variable</th><th></th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>Tendon strand section area</td><td></td><td></td><td></td></tr><tr><td>Main span</td><td>(A<sub>TSM</sub>)</td><td>70.0 mm<sup>2</sup></td><td>180.0 mm<sup>2</sup></td></tr><tr><td>Side spans</td><td>(A<sub>TSS</sub>)</td><td>70.0 mm<sup>2</sup></td><td>180.0 mm<sup>2</sup></td></tr><tr><td>Tendon inflection point - Side spans*</td><td>(α<sub>L,TSS</sub>)</td><td>0.50</td><td>0.80</td></tr><tr><td>Girder cross section height</td><td>(H)</td><td>0.80 m</td><td>2.00 m</td></tr><tr><td>Reinforcement area</td><td></td><td></td><td></td></tr><tr><td>Upper section - Part 1</td><td>(A<sub>T1</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Upper section - Part 2</td><td>(A<sub>T2</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Upper section - Part 3</td><td>(A<sub>T3</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Lower section - Part 1</td><td>(A<sub>B1</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Lower section - Part 2</td><td>(A<sub>B2</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Lower section - Part 3</td><td>(A<sub>B3</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Lower section - Part 4</td><td>(A<sub>B4</sub>)</td><td>3 000 mm<sup>2</sup></td><td>10 000 mm<sup>2</sup></td></tr><tr><td>Reinforcement distribution length*</td><td></td><td></td><td></td></tr><tr><td>Upper section - Part 1</td><td>(α<sub>L,T1</sub>)</td><td>0.10</td><td>0.95</td></tr><tr><td>Upper section - Part 3</td><td>(α<sub>L,T3</sub>)</td><td>0.10</td><td>0.95</td></tr><tr><td>Lower section - Part 2</td><td>(α<sub>L,B2</sub>)</td><td>0.10</td><td>0.95</td></tr><tr><td>Lower section - Part 4</td><td>(α<sub>L,B4</sub>)</td><td>0.10</td><td>0.95</td></tr></tbody></table>

The cEHVI may then be calculated by combining the PF

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Girder cross-sectional variables and geometry</figcaption></figure>

<figure><img src="image_11.png" /></figure>

from Eq. (21b) with that of EHVI in Eq. (26):

$$ \text{cEHVI}(\mathbf{x}, \mathcal{P}, r) = \text{PF}(\mathbf{x}) \cdot \text{EHVI}(\mathbf{x}, \mathcal{P}, r) \quad (27) $$

The next candidate design to be sampled, $\mathbf{x}_{N+1}$ is then calculated by maximizing the cEHVI:

$$ \mathbf{x}_{N+1} = \arg \max_{\mathbf{x} \in \mathcal{X}} \text{cEHVI}(\mathbf{x}, \mathcal{P}, r) \quad (28) $$

The cEHVI and the PTMOO frameworks share a common goal, guiding the optimization process toward designs with the most promising improvement. However, they differ fundamentally in how they define and pursue improvement.

* PTMOO focuses on a single preferred design. The acquisition function is computed as the expected improvement, conditioned on feasibility, relative to a specific target, the preferred trade-off point $f(\mathbf{x}_N^*)$. This leads to a point-based improvement strategy, more exploitative and aligned with user preferences.

* cEHVI, on the other hand, computes the expected increase in hypervolume over the entire current approximation of the Pareto front, conditioned on feasibility. This makes it inherently set-based and more diversity-oriented, as it rewards candidate solutions that extend or fill gaps in the front.

<footnote>*The lengths are relative to the respective span lengths</footnote>
<table><caption>**Table 2** Materials applied in the bridge girder</caption><thead><tr><td>Material</td><td>Quality</td></tr></thead><tbody><tr><td>Concrete</td><td>B45</td></tr><tr><td>Reinforcement bars</td><td>B500NC</td></tr><tr><td>Post-tension tendons</td><td>Y1860S7</td></tr></tbody></table>

The objective function values of cEHVI are typically normalized to [0, 1]¹ⁿⱼ and a reference point $r = 1$¹ⁿⱼ used (Mathern et al. 2021). This approach is also followed in this article. Alternatively, a reference point using the objective values determined using Eq. (22) could be applied, which would limit the search space and avoid the extreme non-dominated points. As the cEHVI in this article is intended only to compare against the approach outlined in Algorithm 1, the typical reference point of $r = 1$¹ⁿⱼ is used.

<table><caption>**Table 3** Applied loads</caption><tbody><tr><td colspan="2">Permanent loads</td></tr><tr><td>Selfweight concrete</td><td>25.00 kN/m<sup>3</sup></td></tr><tr><td>Surfacing</td><td>2.50 kN/m<sup>2</sup></td></tr><tr><td>Edge beams</td><td>8.75 kN/m</td></tr><tr><td>Railing</td><td>0.50 kN/m</td></tr><tr><td>Creep</td><td>Annex B.1 in Standard Norge (2021)</td></tr><tr><td>Shrinkage</td><td>Annex B.2 in Standard Norge (2021)</td></tr><tr><td>Post-tension</td><td>Section 5.10 in Standard Norge (2021)</td></tr><tr><td colspan="2">Variable loads</td></tr><tr><td>Uniform traffic - Lane 1</td><td>5.40 kN/m<sup>2</sup></td></tr><tr><td>Uniform traffic - Lane 2</td><td>2.50 kN/m<sup>2</sup></td></tr><tr><td>Uniform traffic - Remaining area</td><td>2.50 kN/m<sup>2</sup></td></tr><tr><td>Tandem traffic - Lane 1</td><td>2x300 kN</td></tr><tr><td>Tandem traffic - Lane 2</td><td>2x200 kN</td></tr><tr><td>Max. air temperature</td><td>36.0 °C</td></tr><tr><td>Min. air temperature</td><td>-30.0 °C</td></tr></tbody></table>

<table><caption>**Table 4** Ultimate Limit State (ULS) load factors</caption><thead><tr><th rowspan="2">Dominant load</th><th colspan="3">Permanent loads*</th><th colspan="2">Variable loads</th></tr><tr><th>Selfweight</th><th>Deformation</th><th>Tendons</th><th>Traffic</th><th>Temperature</th></tr></thead><tbody><tr><td>Selfweight</td><td>1.35/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>0.95</td><td>0.84</td></tr><tr><td>Traffic</td><td>1.20/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>1.35</td><td>0.84</td></tr><tr><td>Temperature</td><td>1.20/1.00</td><td>1.00/0.00</td><td>1.10/0.90</td><td>0.95</td><td>1.20</td></tr></tbody></table>

# 3 Case study

The bridge considered is a three-span, post-tensioned concrete road bridge with 2 traffic lanes. The side spans are assumed to be 20 m each and the central span is 30 m; thus, the entire bridge is 70 m. A perspective view and a view of the girder section are illustrated in Fig. 5. In this article, only the bridge girder is considered during the optimization.

## 3.1 Variables and geometry

For the case study, 15 design variables are considered, as indicated in Table 1, where the upper and lower bounds of the variables are also provided. These design variables are considered to be the most important for the girder design.

* Tendons provide a compressive axial force to the girder and also a moment which counteracts the girder self-weight and additional live loads.

* Concrete has a relatively low tensile strength, and the reinforcement is used to strengthen the areas of the girder which is subjected to tensile loads, in addition to reducing crack widths.

* The distribution of the reinforcement is partitioned into what is considered as reasonable distribution areas, depending on what part of the section is most subjected to tension or compression.

* The girder's height is the most important parameter for the girder's bending stiffness about the weak axis.

The position of the variables along the bridge girder is illustrated in Fig. 6, and the cross-sectional geometry and the position of the variables in the girder's cross section are illustrated in Fig. 7.

In Fig. 7a, the geometric values of the girder are shown. As indicated in the figure, only the height of the cross section is considered as a variable. The upper width must be constant to accommodate a certain number of traffic lanes, and the bottom width must be constant to ensure adequate spacing for the tendons. Remaining geometry is kept fixed as lateral force on the girder is not considered, and the applied values ensure a realistic design.

<footnote>*The load factors are given as unfavorable/favorable</footnote>
<table><caption>**Table 5** Serviceability Limit State (SLS) load factors</caption><thead><tr><th rowspan="2">Comb.</th><th rowspan="2">Dominant load</th><th colspan="3">Permanent loads</th><th colspan="2">Variable loads</th></tr><tr><th>Selfweight</th><th>Deformation*</th><th>Tendons</th><th>Traffic</th><th>Temperature</th></tr></thead><tbody><tr><td rowspan="2">Char.</td><td>Traffic</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>1.00</td><td>0.70</td></tr><tr><td>Temperature</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>0.70</td><td>1.00</td></tr><tr><td rowspan="2">Freq.</td><td>Traffic</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>0.70</td><td>0.20</td></tr><tr><td>Temperature</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>0.20</td><td>0.60</td></tr><tr><td rowspan="2">Perm.</td><td>Traffic</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>0.50</td><td>0.20</td></tr><tr><td>Temperature</td><td>1.00</td><td>1.00/0.00</td><td>1.00</td><td>0.20</td><td>0.50</td></tr></tbody></table>

<table><caption>**Table 6** Constraints and their limits considered for the MOO of the bridge girder</caption><thead><tr><th>Ultimate Limit State (ULS)</th><th>Reference</th></tr></thead><tbody><tr><td>Tendon strain</td><td>10.0 %</td></tr><tr><td>Concrete compressive strain</td><td>-3.5 %</td></tr><tr><td>Reinforcement strain</td><td>30.0 %</td></tr><tr><td>Serviceability Limit State (SLS)</td><td></td></tr><tr><td>Concrete compressive stress (SLS-Char.)</td><td>-27.0 MPa</td></tr><tr><td>Concrete compression at tendons (SLS-Perm.)</td><td>0.0 MPa</td></tr><tr><td>Crack width (SLS-Freq.)</td><td>0.13 mm</td></tr><tr><td>Maximum deflection due to traffic*</td><td>L<sub>i</sub>/350</td></tr></tbody></table>

### 3.1.1 Materials

The material quality assumed in the girder of the concrete, reinforcement bars, and post-tension system is rendered in Table 2.

### 3.1.2 Loading

The girder is subjected to permanent and variable loads. The permanent loads are comprised of self-weight, creep, shrinkage, super dead loads, and the post-tension force. The variable loads comprise of traffic and temperature loads. The total carriage width of the girder is 7.0 m, which enables two notional lanes. The edge beams, on which the railing is mounted, are assumed suspended to the girders cantilevers. Table 3 shows the applied loads. As the cross section and amount of post-tension varies, references to the applied regulations are given when relevant.

## 3.2 Load combinations

The internal loads from the FEM simulations are combined to what is known as load combinations, which reflects load combination to which the structure is, with different likelihood and duration, subjected to during its lifetime. In Table 4, the load factors for the Ultimate Limit State (ULS) are shown.

<table><caption>**Table 7** Material monetary and environmental costs</caption><thead><tr><th>Material</th><th>Environmental cost</th><th>Monetary cost</th></tr></thead><tbody><tr><td>Concrete (grade B45)</td><td>360 kgCO<sub>2</sub>-eq/m<sup>3</sup></td><td>4 200 NOK/m<sup>3</sup></td></tr><tr><td>Normal reinforcement (grade B500NC)</td><td>4 443 kgCO<sub>2</sub>-eq/m<sup>3</sup></td><td>295 945 NOK/m<sup>3</sup></td></tr><tr><td>Tendons (grade Y1860)</td><td>12 kgCO<sub>2</sub>-eq/mMN</td><td>700 NOK/mMN</td></tr></tbody></table>

For the Serviceability Limit State (SLS), which are load combinations associated with the serviceability of the structure, three main load combinations are considered:

* Char: Characteristic

* Freq: Frequent

* Perm: Quasi-permanent

The load factors for the load combinations with associated dominant load are shown in Table 5.

## 3.3 Constraints

The constraints considered for the MOO of the bridge girder are given in Table 6. The constraints are associated with both the Ultimate and Serviceability Limit State (ULS and SLS) load combinations, in addition to deflections solely due to the traffic loads. For the constraints concerning the

<footnote>* The load factors are given as unfavorable/favorable</footnote>

<footnote>*L<sub>i</sub> is the span length of the span where the deflection is considered</footnote>
reinforcement, it is separated into reinforcement at the upper and lower part of the girder section; thus, a total of nine constraints are considered.

The constraints are derived from numerical analyses and the results from along the entire bridge are used. Thus, the constraints are represented by high-dimensional, functional outputs.

## 3.4 Objectives

The goal of this MOO process is to minimize both the monetary and environmental costs associated with the materials of the bridge girder. Though this case study is limited to two objectives, expanding a MOO process to additional objectives is straightforward; see, e.g., Yepes Llorente et al. (2024). The environmental cost is for the case study determined based on CO₂ equivalent emissions derived from the life cycle assessment (LCA) of the materials used in the structure, with data sourced from the Norwegian Public Roads Administration's LCA tool, VegLCA (Norwegian 2022). Similarly, the monetary cost is calculated using material costs considered representative of Norwegian infrastructure projects. The values applied to the different materials are provided in Table 7.

# 4 Optimizing the girder

## 4.1 Numerical analyses

The bridge girder is analyzed using a Finite Element Method (FEM) model. The FEM model is a beam element model, which is modeled and analyzed using a software called NovaFrame, which is an in-house software at Dr. Ing. A.

(2020). The girder is modeled using 60 beam elements (61 nodes). The analyses are performed using linear elastic theory to derive sectional forces. The sectional forces from the different permanent and variable loads are then combined for the various limit states and automatically sorted to determine the most adverse load combinations. The software then has a module that allows the stresses and strains in the concrete, reinforcement, and tendons at the section to be calculated by determining the equilibrium between the external and internal forces of a monotone stress field using a discrete approximation to Green's theorem.

## 4.2 Approximating simulation results

LHS with normalized distributions for each design variable is carried out using 30 samples. Hence, $N_s = 30$ independent analyses are performed, each corresponding to a unique combination of design variables $x$, to determine the initial observation set $D_0$. As described in Section 3.3, nine constraints are considered. For each constraint $i$, the numerical simulations produce results at the 61 nodal points and the 60 element midpoints, yielding $N_h = 61 + 60 = 121$ spatial evaluations. The constraints along the girder are given by:

$$g_i = \begin{cases} UR_i = \frac{F_{Ed,i}}{F_{Rd,i}}, & i \neq \text{concrete compression at tendons} \\ F_{Ed,i}, & i = \text{concrete compression at tendons} \end{cases} \qquad (29)$$

where $UR_i$ is used as an abbreviation of the utilization ratio for constraint $i$. $F_{Ed,i}$ is the subjected load for constraint $i$, $F_{Rd,i}$ its capacity (see Table 6), and for the constraint

<figure><img src="image_12.png" /><figcaption>Fig. 8 Tendon strain and corresponding POD mode shapes of the 30 initial observations along the bridge. The 4 POD mode shapes correspond to the necessary number of modes with $\epsilon_{POD}^2 = 0.01$</figcaption></figure>
“concrete compression at tendons,” we simply take the raw analysis result. The constraints are limited as:

$$g_i \le \lambda_i, \quad \lambda_i = \begin{cases} 1.0, & i \ne \text{concrete compression at tendons} \\ 0.0, & i = \text{concrete compression at tendons} \end{cases} \qquad (30)$$

This signifies that the subjected load must be lower than the capacity, and for the constraint “concrete compression at tendons,” the concrete must be in compression. The constraints are then approximated using POD with GPR to interpolate for unobserved variable combinations. As detailed in Sect. 2.2, the number of POD modes is selected using Eq. 4 with $\epsilon_{POD}^2 = 0.01$. Because each constraint has different POD modes, the required number of POD modes varies across constraints. Figure 8 illustrates the tendon strain along the bridge for the initial 30 observations, together with the first four POD modes which were sufficient to keep the projection error below 0.01 for the initial observations of the tendon strain. In Fig. 8b, the modes are arranged in descending order of significance, with “Mode 1” denoting the most significant mode.

To assess the effectiveness of the KPLS method, the constraints are approximated using POD with GPR and no kernel hyperparameter reduction, and using POD with KPLS where the number of components included in the KPLS ranges from 2 to 5.

The objectives f, i.e., the monetary and environmental cost, are obtained for the complete girder as scalar values, and are thus not approximated using POD model reduction, though they are approximated using GPR with and without kernel hyperparameter reduction in the same manner as for the constraints. The objective functions have a relatively simple shape and could also be accurately modeled using a surrogate modeling technique that is simpler than GPR. The benefit of using GPR for the objective functions is, as mentioned before, that it allows the optimization process to leverage uncertainties between observations. GPR is also relatively fast to construct when a limited number of simulations are performed, which is typical in engineering design.

The calculations are implemented in Python using TensorFlow (Abadi et al. 2015). The estimation of hyperparameters is performed using the gradient-based L-BFGS algorithm from Scipy, where the gradients are calculated using algorithmic differentiation in TensorFlow (Virtanen et al. 2020).

## 4.3 Optimization routine

The MOO processes described in Sect. 2.5 are performed for 50 iterations, increasing the sample size $N_s$ with one new observation for each iteration. This is repeated 20 times, each run using the same initial set of 30 samples, to assess

its consistency. The hyperparameter setting for the genetic MOO algorithm NSGA-II were;

* Population size: 200

* Number of generations: 100

To also assess the effect of the KPLS hyperparameter reduction, the following GPR approximations are considered:

* GPR: No hyperparameter reduction

* KPLS<sub>h</sub>: KPLS with $\forall h \in \{2, 3, 4, 5\}$ components

To implement cEHVI, we employed box decomposition (Yang et al. 2019). Optimization of Eq. (28) was carried out with the genetic algorithm (GA) provided by Pymoo (Blank and Deb 2020), ensuring both MOO routines shared the same solution strategy. The hyperparameters were configured using the same hyperparameters as for the NSGA-II. It was also considered optimizing Eq. (28) directly using a gradient-based optimizer (GBO). This was performed using the L-BFGS GBO algorithm in Scipy (Virtanen et al. 2020). The GBO employed a multi-start strategy, where 1000 initial designs were generated at each iteration using uniform Latin hypercube sampling (LHS), and the best candidate was selected as the starting point for the optimization.

For the cEHVI and GBO, only GPR with no hyperparameter reduction and KPLS with 5 components were considered.

## 4.4 Benchmark

As described in the introduction, the traditional method of bridge design is based upon expert opinion and best guesses or rules of thumb which is then followed by multiple iterations to achieve the final design. The traditional design process for a concrete road bridge is illustrated in Fig. 9. The process is based on the authors’ experience from working several years in the bridge design industry. The process follows a set of consecutive steps in which one considers one variable or component at a time. When an initial design has been determined, simulations are performed and the utilization ratio of the relevant constraints is calculated. The designer must then evaluate whether certain areas need additional reinforcement or tendons to meet the constraints, or if this can be reduced in certain areas. Decreasing the amount of materials is considered beneficial as it reduces the cost of the structure.

This process was performed for the bridge girder presented in the case study of Sect. 3. A design was determined in which the minimum reinforcement required by the Eurocodes and national regulations was applied, and the tendon
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Traditional bridge girder design process</figcaption></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 10** Environmental and monetary cost of the preferred design $\mathbf{x}_N^*$ for the PTMOO. Solid line represents the median value, and the shaded area represents the 50 % confidence interval of the 20 repetitions</figcaption></figure>

<table><caption>**Table 8** Variable combination of the preferred design at the last iteration $\mathbf{x}_{50}^*$ of the PTMOO</caption><thead><tr><th>Variable</th><th>Unit</th><th>ROT</th><th>GPR*</th><th>KPLS<sub>5</sub>*</th></tr></thead><tbody><tr><td colspan="5">Tendon strand section area</td></tr><tr><td>Main span</td><td>(A<sub>TSM</sub>) [mm<sup>2</sup>]</td><td>127.8</td><td>130.5 ± 12.8</td><td>139.8 ± 1.6</td></tr><tr><td>Side spans</td><td>(A<sub>TSS</sub>) [mm<sup>2</sup>]</td><td>62.2</td><td>74.1 ± 15.0</td><td>76.0 ± 1.3</td></tr><tr><td>Tendon inflection point - Side spans</td><td>(α<sub>TSS</sub>) [-]</td><td>0.60</td><td>0.54 ± 0.04</td><td>0.64 ± 0.02</td></tr><tr><td>Girder cross section height</td><td>(H) [m]</td><td>1.50</td><td>1.10 ± 0.10</td><td>1.01 ± 0.01</td></tr><tr><td colspan="5">Reinforcement area</td></tr><tr><td>Upper section - Part 1</td><td>(A<sub>T1</sub>) [mm<sup>2</sup>]</td><td>15 134</td><td>4 366 ± 1 613</td><td>6 793 ± 74</td></tr><tr><td>Upper section - Part 2</td><td>(A<sub>T2</sub>) [mm<sup>2</sup>]</td><td>12 382</td><td>6 491 ± 2 407</td><td>8 792 ± 360</td></tr><tr><td>Upper section - Part 3</td><td>(A<sub>T3</sub>) [mm<sup>2</sup>]</td><td>16 234</td><td>4 908 ± 1 464</td><td>5 745 ± 131</td></tr><tr><td>Lower section - Part 1</td><td>(A<sub>B1</sub>) [mm<sup>2</sup>]</td><td>5 928</td><td>5 223 ± 2 175</td><td>5 454 ± 527</td></tr><tr><td>Lower section - Part 2</td><td>(A<sub>B2</sub>) [mm<sup>2</sup>]</td><td>5 928</td><td>6 339 ± 2 539</td><td>7 197 ± 52</td></tr><tr><td>Lower section - Part 3</td><td>(A<sub>B3</sub>) [mm<sup>2</sup>]</td><td>5 928</td><td>6 538 ± 1 647</td><td>8 120 ± 153</td></tr><tr><td>Lower section - Part 4</td><td>(A<sub>B4</sub>) [mm<sup>2</sup>]</td><td>5 928</td><td>8 259 ± 1 901</td><td>9 101 ± 148</td></tr><tr><td colspan="5">Reinforcement distribution length</td></tr><tr><td>Upper section - Part 1</td><td>(α<sub>T1</sub>) [-]</td><td>0.85</td><td>0.81 ± 0.09</td><td>0.76 ± 0.02</td></tr><tr><td>Upper section - Part 3</td><td>(α<sub>T3</sub>) [-]</td><td>0.70</td><td>0.53 ± 0.32</td><td>0.89 ± 0.04</td></tr><tr><td>Lower section - Part 2</td><td>(α<sub>B2</sub>) [-]</td><td>0.50</td><td>0.66 ± 0.20</td><td>0.64 ± 0.02</td></tr><tr><td>Lower section - Part 4</td><td>(α<sub>B4</sub>) [-]</td><td>0.50</td><td>0.61 ± 0.19</td><td>0.63 ± 0.02</td></tr></tbody></table>

quantity was reduced until no further reduction was possible which was governed by the constraint *Concrete compression at tendons* of Table 6. This design is abbreviated as ROT in the results and serves as a benchmark for the MOO process.

# 5 Results

## 5.1 Comparing PTMOO to benchmark design

The environmental and monetary cost of the preferred design $\mathbf{x}_N^*$ according to Eq. (22), at each iteration of the PTMOO process is shown in Fig. 10. The figure shows the median of the objective functions $f(\mathbf{x}_N^*)$ in each iteration with solid lines, and the 50 % confidence interval as the corresponding shaded area. The figure shows that the PTMOO have an environmental cost which is about 20 % lower than the ROT,

<footnote>* Values are given as mean ± one standard deviation</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 11** Variable combination of the preferred design at the last iteration $\mathbf{x}_{50}^*$ of the PTMOO relative to ROT. The solid lines show the mean value and the shaded areas are the mean ± one standard deviation</figcaption></figure>

and simultaneously a monetary cost that is between 10 % and 15 % lower. The figure also shows that the PTMOO using KPLS converge faster, i.e., reaching a plateau at fewer iterations, than the method that utilizes only POD with full GPR. The KPLS reduces the number of hyperparameters, and since there are limited observations at the start of the optimization process, it seems to be more capable of providing accurate predictions at the early phase. The full GPR appears to converge to a design with a slightly higher environmental cost, albeit slightly lower monetary cost. In the case study, the KPLS with three weight components converge the fastest, while increasing the number of weights beyond three results in a slightly lower convergence rate. All KPLS methods seem to converge to approximately the same optimal environmental and monetary costs.

In Table 8, the variable combination of the preferred design at the last iteration $\mathbf{x}_{50}^*$ for ROT, GPR, and KPLS$_5$ is shown. For GPR and KPLS$_5$, the mean ± one standard deviation for the 20 evaluations is shown. The table shows that the ROT has a higher girder section with less tendon area, and a higher reinforcement area in the upper section than the PTMOO methods. The minimum reinforcement requirement of the Eurocodes and national regulations is not considered in the PTMOO processes, making it possible that the results do not comply with these requirements. In practice, such requirements can be satisfied through minor post-optimization adjustments, which are not expected to significantly alter the trade-off behavior, cost or performance of the designs. Since the minimum reinforcement has a linear relationship with the section height, the lower section height

in the PTMOO process will result in a lower minimum reinforcement compared to the ROT. It should also be noted that in practice, the reinforcement and tendon amounts must be according to the suppliers' available dimensions.

In Table 8, we also note that the GPR has a significantly higher standard deviation than the KPLS$_5$ for the preferred feasible design variables, which is also apparent in the wider confidence interval of Fig. 10. In Fig. 11, the variable combination of the preferred design at the last iteration $\mathbf{x}_{50}^*$ relative to ROT is illustrated. The solid lines show the mean value, and the shaded areas are the mean ± one standard deviation. The likelihood function of Eq. (13), which is used to determine the hyperparameters, is often multi-modal. Because GPR must learn many more hyperparameters, its training is more susceptible to getting stuck in local optima, which can degrade prediction accuracy. In contrast, KPLS$_5$ is limited to just five hyperparameters, greatly reducing its flexibility and thus its vulnerability to suboptimal fits. This difference in model complexity likely accounts for the larger standard deviation observed in the preferred feasible design variables when using GPR.

Figure 12 presents the utilization ratio of the constraints for the preferred solution along the bridge, referred to as the *station* in the figure, at the final iteration of the PTMOO process. The results are shown for a randomly selected case from the 20 repetitions, comparing the process using KPLS$_5$, GPR, and the benchmark ROT. Constraints are satisfied as long as they are below 1.0 or 0.0 for the *SLS: Tendons* case. As the figure shows, all constraints are satisfied. What is interesting to note is that the KPLS$_5$ and the GPR have relatively similar utilizations, sometimes the KPLS$_5$ is slightly higher and for others, it is the opposite. In all cases, the ROT has the lowest utilization, except for the *SLS: Tendons*, where utilization levels are relatively similar. This constraint governed the tendon amount in the ROT design process and appears to be quite active in the PTMOO processes as well. Higher utilization is not inherently disadvantageous, as long as it remains within the allowable limit. In fact, it appears that increasing the overall level of utilization improves the efficiency of the load bearing capacity of the structure, ultimately reducing both environmental and financial costs. The figure demonstrates that considering multiple variables simultaneously, as is done in the PTMOO processes, can lead to greater structural utilization, which further contributes to cost and environmental benefits.

## 5.2 Comparing cEHVI to PTMOO

Figure 13 depicts, at each iteration and in accordance with Eq. (22), the monetary and environmental cost of the preferred design $\mathbf{x}_N^*$ found by both optimization schemes. The solid
<figure><img src="image_17.png" /><figcaption>**Fig. 12** Constraints along the bridge for the preferred solution at final iteration for ROT, GPR, and KPLS<sub>5</sub> using PTMOO. SLS and ULS are the limit states described in Sect. 3.2. “Rein.” is normal reinforcement,</figcaption></figure>

“bot.” is at the bottom of the girder, and “top” is at the top of the girder. Limit is 1.0 for all constraints except “SLS: Tendons” where it is 0.0

curves show the PTMOO results, with the shaded region indicating the 50% confidence interval. The dashed curves show the cEHVI results, with the hatched region indicating the corresponding 50% confidence interval. Figure 13 demonstrates that, over 50 iterations, PTMOO identifies a design incurring both lower monetary and environmental costs. Although KPLS<sub>5</sub> appears to plateau after roughly 20 iterations, cEHVI is still improving at iteration 50. This slower convergence is expected. cEHVI explores multiple feasible solutions under arbitrary objective trade-offs, whereas PTMOO targets a single solution with its trade-off fixed by Eq. (22).

## 5.3 Comparing GBO to PTMOO

In Fig. 14, the results of PTMOO are compared with GBO, showing the environmental and monetary costs of the preferred design $x_N^*$ as defined in Eq. (1). The figure reports the median performance, with solid lines for PTMOO and dashed lines for GBO, while the 50% confidence intervals are indicated by shaded regions for the PTMOO and hatched regions for the GBO.
<figure><img src="image_18.png" /><figcaption>**Fig. 13** cEHVI compared to PTMOO. Environmental and monetary cost of the preferred design $\boldsymbol{x}_N^*$ at each iteration according to Eq. (3). Solid line represents the median value, and the shaded area</figcaption></figure>

represents the 50 % confidence interval of the 20 repetitions for the PTMOO. The dashed line and the hatched shaded area represent the corresponding results for the cEHVI approach,

GBO without KPLS (GPR-GBO) outperforms the convergence rate of PTMOO without KPLS (GPR-TMOO), converging more rapidly to designs with lower environmental costs, albeit a higher monetary cost. When KPLS is included, the KPLS appears to have a degrading effect on the optimization process for the GBO approach, which is opposite to what is observed for the PTMOO. The GBO tends to converge to slightly higher monetary costs than the PTMOO. The faster convergence of GBO reflects its reliance on gradient information to efficiently exploit promising descent directions, whereas PTMOO maintains a broader exploration of the design space, enabling the identification of better optima. The use of KPLS appears to improve the PTMOO by simplifying the design space through dimensionality reduction, while it appears to limit the search space of the GBO.

## 6 Conclusion

A multi-objective constrained Bayesian optimization process with a preferred trade-off for the objectives has been presented and applied to a 3-span post-tensioned concrete road bridge girder. Fifteen design variables were considered in the optimization process to minimize both the environmental and the monetary cost of the structure, while also fulfilling nine constraints. Optimization was performed using EI\_c Bayesian acquisition functions, where the improvement was measured relative to the current preferred feasible solution. Both constraints and objectives were approximated using GPR. Due to the high dimensionality of the constraints, which were derived from numerical simulations, they were approximated by using POD with GPR to approximate unobserved design variable combinations. In

addition, the use of KPLS was explored to assess its potential to enhance the effectiveness of the optimization process.

As a benchmark for the case study, design variables and corresponding environmental and monetary costs were determined using a traditional bridge girder design process. The results showed that multi-objective optimization methods with a preferred trade-off for objectives were able to reduce the environmental cost with about 20 % and the monetary cost with between 10 % and 15 % while complying to all constraints, when compared to the benchmark. The results also showed that applying KPLS can result in a faster convergence, i.e., fewer number of iterations before plateauing, than GPR with a full number of kernel parameters. For the case study, a minimum of three weight components were necessary to achieve this. The KPLS and GPR with a full number of kernel parameters converged to design variables that resulted in about the same monetary and environmental costs.

The MOO routine with a preferred trade-off was compared to the commonly used cEHVI MOO routine, where the trade-off is arbitrary. The results showed that the approach with a preferred trade-off was able to obtain a solution where both the environmental and monetary cost were lower than that achieved using cEHVI with the same number of iterations. This is as expected, as the cEHVI expends effort exploring the full set of non-dominated solutions. For the cEHVI, the objectives were normalized and bounded by a reference point encompassing the full possible Pareto front. A more aggressive approach, if the extremities of the Pareto front is not of interest, would be to use the preferred trade-off solution, e.g., Eq. (3), as a reference point.

The MOO method presented in this article is applicable to any optimization problem where data are expensive or time-consuming to generate, resulting in a limited amount
<figure><img src="image_19.png" /><figcaption>**Fig. 14** GBO compared to PTMOO. Environmental and monetary cost of the preferred design $\boldsymbol{x}_N^*$ at each iteration according to Eq. (22). Solid line represents the median value, and the shaded area</figcaption></figure>

represents the 50 % confidence interval of the 20 repetitions for the PTMOO. The dot-dashed line and the hatched shaded area represent the corresponding results for the GBO approach

of available data, and a moderate number of design variables are considered. These are conditions that are common in many engineering design problems that involve numerical simulations. The article has also shown that considering multiple variables simultaneously can lead to higher structural utilization, ultimately resulting in lower environmental and monetary costs. A bottleneck for the applicability of the presented MOO is generating enough observations to make accurate predictions. Since simulations can be time-consuming, reducing their number would be beneficial. Future research will explore the potential of using deep Gaussian process regressors for the MOO method to determine whether they can improve prediction accuracy with less data.

**Acknowledgements** A special thanks to assistant professor Frederico Afonso at Instituto Superior Técnico (IST) - University of Lisbon, for his thoughts and advises on the multi-objective optimization routine.

**Author Contributions** The authors' contribution is as follows: Heine Røstum was responsible for the conceptualization, methodology, formal analyses, data curation, software writing, and original draft work. Joseph Morlier was responsible for supervision, conceptualization, and methodology. Sebastien Gros and Ketil Aas-Jakobsen were responsible for the supervision and project administration.

**Funding** Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital - Trondheim University Hospital). This work is funded by the Research Council of Norway [project number 332126] and by Dr Ing A Aas-Jakobsen AS.

**Data Availability** The results from the finite element analyses derived from the case study of Sect. 3 can be made available upon request by e-mailing the corresponding author.

## Declarations

**Conflict of interest** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Abadi M, Agarwal A, Barham P, Brevdo E, Chen Z, Citro C, Corrado GS, Davis A, Dean J, Devin M, Ghemawat S, Goodfellow I, Harp A, Irving G, Isard M, Jia Y, Jozefowicz R, Kaiser L, Kudlur M, Levenberg J, Mané D, Monga R, Moore S, Murray D, Olah C, Schuster M, Shlens J, Steiner B, Sutskever I, Talwar K, Tucker P, Vanhoucke V, Vasudevan V, Viégas F, Vinyals O, Warden P, Wattenberg M, Wicke M, Yu Y, Zheng X (2015) TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems. Software available from tensorflow.org. https://www.tensorflow.org/

Amsallem D, Farhat C (2012) Stabilization of projection-based reduced-order models. Int J Numer Meth Eng 91(4):358–377

Blank J, Deb K (2020) Pymoo: multi-objective optimization in python. IEEE Access 8:89497–89509

Bouhlel MA, Bartoli N, Otsmane A, Morlier J (2016) Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction. Struct Multidiscip Optim 53:935–952
Brunton SL, Kutz JN (2022) Data-driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, Cambridge, England

Buljak V, Maier G (2011) Proper orthogonal decomposition and radial basis functions in material characterization based on instrumented indentation. *Eng Struct* 33(2):492–501

Chatterjee A (2000) An introduction to the proper orthogonal decomposition. Current science, 808–817

Chiplunkar A, Bosco E, Morlier J (2018) Gaussian process for aerodynamic pressures prediction in fast fluid structure interaction simulations. In: Advances in Structural and Multidisciplinary Optimization: Proceedings of the 12th World Congress of Structural and Multidisciplinary Optimization (WCSMO12) 12, 221–233. Springer

Daulton S, Balandat M, Bakshy E (2020) Differentiable expected hypervolume improvement for parallel multi-objective bayesian optimization. *Adv Neural Inf Process Syst* 33:9851–9864

Deb K, Pratap A, Agarwal S, Meyarivan T (2002) A fast and elitist multiobjective genetic algorithm: nsga-ii. *IEEE Trans Evol Comput* 6(2):182–197

Dr. Ing. A. Aas-Jakobsen A/S: NovaFrame: In-house bridge design software. Developed by Dr. Ing. A. Aas-Jakobsen A/S, Oslo, Norway (2020)

Drakoulas G, Gortsas TV, Polyzos D (2024) Physics-based reduced order modeling for uncertainty quantification of guided wave propagation using bayesian optimization. *Eng Appl Artif Intell* 133:108531

Eckart C, Young G (1936) The approximation of one matrix by another of lower rank. *Psychometrika* 1(3):211–218

Frazier PI (2018) Bayesian optimization. In: Recent Advances in Optimization and Modeling of Contemporary Problems, 255–278. *Informs, MD, USA*

Frazier PI, Wang J (2016) Bayesian optimization for materials design. *Information science for materials discovery and design*, 45–75

Gardner JR, Kusner MJ, Xu ZE, Weinberger KQ, Cunningham JP (2014) Bayesian optimization with inequality constraints. *ICML 2014:937–945*

Gooijer BM, Havinga J, Geijselaers HJ, Boogaard AH (2021) Evaluation of pod based surrogate models of fields resulting from nonlinear fem simulations. *Advanced Modeling and Simulation in Engineering Sciences* 8(1):25

Guo M, Hesthaven JS (2018) Reduced order modeling for nonlinear structural analysis using gaussian process regression. *Comput Methods Appl Mech Eng* 341:807–826

Helland IS (1988) On the structure of partial least squares regression. *Communications in statistics-Simulation and Computation* 17(2):581–607

Jeong S, Obayashi S (2005) Efficient global optimization (ego) for multi-objective problem and data mining. In: 2005 IEEE Congress on Evolutionary Computation, 3, 2138–2145. IEEE

Jones DR, Schonlau M, Welch WJ (1998) Efficient global optimization of expensive black-box functions. *J Global Optim* 13(4):455–492

Knowles J (2005) Parego: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. *IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION* 10(1)

Liang Y, Lee H, Lim S, Lin W, Lee K, Wu C (2002) Proper orthogonal decomposition and its applications–part i: theory. *J Sound Vib* 252(3):527–544

Lucia DJ, Beran PS, Silva WA (2004) Reduced-order modeling: new approaches for computational physics. *Prog Aerosp Sci* 40(1–2):51–117

Martins JR, Ning A (2021) Engineering Design Optimization. Cambridge University Press, Cambridge, England

Mathern A, Steinholtz OS, Sjöberg A, Önnheim M, Ek K, Rempling R, Gustavsson E, Jirstrand M (2021) Multi-objective constrained bayesian optimization for structural design. *Struct Multidiscip Optim* 63:689–701

McKay MD, Beckman RJ, Conover WJ (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics* 21(2):239–245

Negrin I, Kripka M, Yepes V (2023) Metamodel-assisted design optimization in the field of structural engineering: A literature review. In: *Structures*, 52, 609–631. Elsevier

Nocedal J, Wright SJ (1999) Numerical Optimization. Springer, New York, NY

Norwegian Public Road Administration: *VegLCA*. Online. v5.11B (2022)

NPRA: N400: Bruprosjektering. The Norwegian Public Roads Administration, Norway (2025)

Paul-Dubois-Taine A, Amsallem D (2015) An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models. *Int J Numer Meth Eng* 102(5):1262–1292

Peherstorfer B, Willcox K (2015) Dynamic data-driven reduced-order models. *Comput Methods Appl Mech Eng* 291:21–41

Penadés-Plà V, García-Segura T, Yepes V (2019) Accelerated optimization method for low-embodied energy concrete box-girder bridge design. *Eng Struct* 179:556–565

Rocha PRB, Sousa Almeida JL, Paula Gomes MS, Junior ACN (2023) Reduced-order modeling of the two-dimensional rayleigh-bénard convection flow through a non-intrusive operator inference. Engineering Applications of Artificial Intelligence 126:106923

Røstum H, Gros S, Aas-Jakobsen K (2025) Constrained bayesian optimization for engineering bridge design. *Struct Multidiscip Optim* 68(1):1–14

Saves P, Van EN, Bartoli N, Lefèbvre T, Diouane Y, Morlier J (2024) Gaussian process for bayesian optimization with mixed hierarchical variables: Application to electric-hybrid aircraft eco-design. In: *Workshop on Bayesian Optimization & Related Topics*

Standard Norge: Eurocode 2: Design of Concrete Structures, Part 1-1: General Rules and Rules for Buildings. Standard Norge, Oslo, Norway (2021). EN 1992-1-1:2004+A1+NA

Tfaily A, Diouane Y, Bartoli N, Kokkolaras M (2024) Bayesian optimization with hidden constraints for aircraft design. *Struct Multidiscip Optim* 67(7):123

United Nations Environment Programme, Yale Center for Ecosystems + Architecture: Building Materials and the Climate: Constructing a New Future (2023). https://wedocs.unep.org/20.500.11822/43293

Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Courapeau D, Burovski E, Peterson P, Weckesser W, Bright J, et al (2020) Scipy 1.0: fundamental algorithms for scientific computing in python. *Nature methods* 17(3), 261–272

Williams CK, Rasmussen CE (2006) Gaussian Processes for Machine Learning. MIT press, Cambridge, MA

Yang K, Emmerich M, Deutz A, Bäck T (2019) Efficient computation of expected hypervolume improvement using box decomposition algorithms. *J Global Optim* 75:3–34

Yeses Llorente L, Morlier J, Sridhara S, Suresh K (2024) A hybrid machine learning and evolutionary approach to material selection and design optimization for eco-friendly structures. *Struct Multidiscip Optim* 67(5):69

Zhao H (2021) A reduced order model based on machine learning for numerical analysis: an application to geomechanics. *Eng Appl Artif Intell* 100:104194

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.