BOIS: Bayesian Optimization of Interconnected Systems

Optimization of natural and engineered systems (e.g., chemical processes, biological systems, materials) is often challenging due to incomplete physical knowledge or the high complexity of experiments and available physical models. As a result, it is necessary to device optimization procedures that effectively combine experimental and model data while mitigating complexity, see Conn et al. (2009). These procedures, often referred to as black-box optimization algorithms, treat the system as a black-box, $f(x)$, and use input/output observation data to direct their search for a solution. One of the most popular paradigms to have emerged in this setting is Bayesian optimization (BO), see Mockus (2012). BO uses input/output data to generate a Gaussian process (GP) model that estimates system performance as well as model uncertainty. These estimates are used to construct an acquisition function (AF) that allows for the search to be tuned to emphasize sampling from either high performance (exploitation) or high uncertainty (exploration) regions; this is a key feature that makes BO an especially powerful global optimizer, see Shahriari et al. (2016).

While the black-box assumption substantially facilitates the implementation of BO (no prior knowledge about the system is needed), it ignores the fact that there is often some form of structural knowledge available (e.g., physics or sparse interconnectivity). For example, when dealing with physical systems, several components might be well-modeled and understood, while others might not. For those that are not, the fundamental principles governing their behavior (e.g., conservation laws, equilibrium, value constraints) are, at least qualitatively, understood. Additionally, sparse connectivity, which provides information on how different components affect each other, is also often known. As a result, the system of interest is usually not an unobservable black-box but rather a “grey-box” that is partially observable with a known structure, see Astudillo and Frazier (2021). Several methods exist for performing grey-box BO; studies have shown that they are able to improve algorithm performance by effectively constraining the search space, resulting in lower sampling requirements and better solutions, see Lu et al. (2021); Raissi et al. (2019); Kandasamy et al. (2017). Of these approaches, the use of composite functions has proven to be one of the most widely used methods for incorporating preexisting knowledge into BO, see Astudillo and Frazier (2019).

A composite function expresses the system performance as $f(x,y(x))$ where $x$ are the system inputs, $f$ is a known scalar function, and $y$ is an unknown vector-valued function that describes the behavior of internal system components. This decomposition shifts the modeling task from estimating the performance function directly to estimating the values of $y$ which serve as inputs to $f(x,y(x))$. Additionally, because $f$ is now a known function, derivative information might be available to understand the effects of $x$ and $y$, see Uhrenholt and Jensen (2019). This approach also readily lends itself to the inclusion of constraints as these are often dependent on internal variables which can be captured by $y$ Paulson and Lu (2022). Framing an optimization problem in this manner is therefore a fairly intuitive approach, especially in chemical engineering where the performance metric is usually an economic cost. For example, the cost equations for equipment, material streams, and utilities are often readily available and it is the parameters these equations rely on (compositions, flowrates, duties) that are unknown. Furthermore, traditional unit operations (heat exchangers, distillation columns, compressors) have significantly better mechanistic models available than those that tend to be more niche (bioreactors, non-equilibrium separators, solids-handling). It then makes sense to construct a composite function where the outer function, $f$, returns the price of the system based on the known cost equations while its inputs, $y$, are the mass and energy flows through the system and are estimated via either mechanistic or data-driven models. Constraints can then be incorporated using values estimated for $y$ to ensure that data-driven models obey fundamental physical laws as well as to enforce more traditional requirements such as product specifications, waste generation, utility consumption, and equipment sizing which are often important in process design.

While setting up a composite function optimization problem seems fairly straightforward, implementing it in a BO setting is not a trivial task. As previously stated, one of the main advantages of BO is the inclusion of uncertainty estimates in the surrogate model, which allows for greater exploration of the design space when compared to a deterministic model, see Boukouvala and Floudas (2017). However, when using a composite function, the GP models generated are of $y$ not $f$. Given that $f$ is the performance measure being optimized, it is necessary propagate the predicted uncertainty from $y(x)$ to $f(x,y(x))$ (i.e, the density of $f$ or desired summarizing statistics must be determined). A Gaussian density for $y(x)$ is directly obtained from the GP model; as a result, when $f$ is a linear model, we can make use of the closure of Gaussian random variables under linear operations to generate the density of $f(x,y(x))$ (which is also a Gaussian). When $f$ is nonlinear, however, a closed-form solution is not readily available, and this operation is usually carried out numerically using sampling methods like Monte Carlo (MC), see Balandat et al. (2020).

Given the increased functionality composite functions can provide to BO and the computational intensity of sampling methods, we propose the Bayesian Optimization of Interconnected Systems (BOIS) framework. BOIS provides a novel method for estimating the distribution of $f$ in a significantly more efficient manner than sampling by linearizing the performance function in the neighborhood of a $y(x)$ of interest. This allows us to construct a local Laplace approximation which we can use to generate closed-form expressions for the mean and uncertainty of $f(x,y(x))$. We test the performance of this approach by conducting extensive numerical experiments on a chemical process optimization case study and compare its performance to standard BO as well as MC-driven composite function BO. Our results illustrate that BOIS is able to outperform standard BO while also providing accurate estimates for the distribution of $f$ at a significantly lower computational cost than MC.

The fundamental challenge of composite function BO algorithm is the lack of closed-form expressions for $m_{f}^{\ell}(x)$ and $\sigma_{f}^{\ell}(x)$, which are needed to build the AF. As previously mentioned, obtaining these when $f$ is not a linear mapping of $y$ is generally not possible. However, when $f$ is a once-differentiable mapping, it is possible to conduct a linearization at the current iterate (as is done in standard optimization algorithms such as Newton’s method). If we choose to represent $f$ as

such that $g(x)$ includes the terms that do not depend on $y$ and $h(x,y)$ includes those that do, we can use a first-order Taylor series expansion to linearize $f$ with respect to $y$ around some reference point $y_{0}$:

where

Using this approach, if we then select some point of interest $x^{\ell+1}$ where the mean and covariance of $y(x)$ are given by the GP model as $\hat{y}^{\ell}=m_{y}^{\ell}(x^{\ell+1})$ and $\hat{\Sigma}^{\ell}=\Sigma_{y}^{\ell}(x^{\ell+1})$ respectively and some reference point in the neighborhood of $\hat{y}^{\ell}$, which we denote as $\hat{y}_{0}^{\ell}$, we can obtain the following estimate for $f$ at $x^{\ell+1}$

Note that we make the assumption that $g(x)$ is also a known, easy-to-evaluate function and, therefore, $g(x^{\ell+1})$ is a deterministic variable. Thus, we are now able to derive at a set of closed-form expressions for the mean and uncertainty of the performance function by making use of the closure of normal distributions under linear transformations.

The proposed framework (which we call BOIS and is shown in Figure 3) is built on the expressions derived in (11). The reason we are able to generate these explicit formulations is due to the manner in which linearizing constructs the density of $f$. While the MC approach is agnostic to the nature of the density of $f$, BOIS approximates it using a Gaussian model around the neighborhood of the iterate $\hat{y}_{0}^{\ell}$. In other words, BOIS generates a local Laplace approximation of the performance function by passing the mean and uncertainty estimates of $y(x)$ given by the GP model into (10). This allows us to obtain closed-form expressions for the statistical moments of $f$, such as $m_{f}^{\ell}$ and $\sigma_{f}^{\ell}$, which can be used for constructing AFs. Given that this approximation of the density of $f$ is Gaussian, it is also possible to obtain expressions for the probabilities and quantiles (to construct different types of AFs). This assumption about the shape of $f$ also means that as the distance between $\hat{y}^{\ell}$ and $\hat{y}_{0}^{\ell}$ grows, the Laplace approximation will result in a worse fit, similar to how the linearization itself becomes less accurate. However, it is important to note that the linearization is updated in an adaptive manner (by linearizing around the current iterate). If we compare BOIS to MC-driven approaches, we can see that at any point $x^{\ell+1}$ of interest, BOIS only samples from the GP once to determine $\hat{y}^{\ell}$ and $\hat{\Sigma}^{\ell}$ and evaluates $f$ once to calculate $f(x^{\ell+1},\hat{y}_{0})$; recall that this is done tens to thousands of times in MC. While BOIS also has to compute (9), this is also done only once and has been shown to have a computational cost similar to that of evaluating $f$ when methods like automatic differentiation are used, see Griewank and Walther (2008); Baur and Strassen (1983). As a result, we are able calculate values for $m_{f}^{\ell}$ and $\sigma_{f}^{\ell}$ at a significantly lower computational cost than when using MC.

We presented a framework, which we refer to as BOIS, to enable the use of composite functions $f(x,y(x))$ in BO to allow for the exploitation of structural knowledge (in the form of physics or sparse interconnections). The key contribution of this work is the derivation of a set of closed-form expressions for the moments of the composite, $m_{f}^{\ell}(x)$ and $\sigma_{f}^{\ell}(x)$, based on adaptive linearizations; this overcomes the tractability challenges of sampling based approaches. These expressions are obtained by linearizing $f(x,y(x))$ in order to generate a Laplace approximation around the current iterate. We demonstrate that BOIS outperforms standard BO by making use of the structure conveyed by $f(x,y(x))$. We also show that the statistical moments calculated by BOIS accurately represent the statistical moments of $f$ and that these estimates can be obtained in significantly less time than sampling-based approaches. As part of our future work, we plan to scale up BOIS and deploy it on high-dimensional systems where BO has traditionally not been applied and to obtain alternative types of AFs. Finally, while the GP is the most popular surrogate model choice in BO, the algorithm is not limited to using only GPs, any probabilistic model can be used. Therefore, we would like to explore the use of alternative models such as warped GPs, see Snelson et al. (2004), RNNs, see Thompson et al. (2023), and reference models, see Lu et al. (2021).