Exploiting Separability in Multi-Scale Grey-Box Bayesian Optimization

# Exploiting Separability in Multi-Scale Grey-Box Bayesian Optimization
Source: [https://arxiv.org/html/2608.03045](https://arxiv.org/html/2608.03045)
Tyler A\. Soderstrom2Brian A\. Korgel1,3Michael Baldea1,4Corresponding author: mbaldea@che\.utexas\.edu

###### Abstract

We consider grey\-box optimization problems where the decision variables naturally partition into black\-box variablesxBBx^\{\\mathrm\{BB\}\}\(as arguments to an expensive black\-box function\) and white\-box variablesxWBx^\{\\mathrm\{WB\}\}, governed by a set of explicit, closed\-form equations that also depend on the output of the black\-box function\. We exploit this separability through a bilevel reformulation: an outer Bayesian optimization \(BO\) to optimize the scalar objective as a function ofxBBx^\{\\mathrm\{BB\}\}alone, while an inner problem solves the white\-box subproblem via global optimization\. The Gaussian process surrogate used in BO is therefore definedℝnBB\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}rather thanℝnWB\+nBB\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}\}, and white\-box constraints are satisfied exactly whenever the inner optimizer converges to a feasible point—without penalty functions, chance constraints, or moment approximations\. On a suite of 13 benchmark problems, bilevel BO achieves lower regret, with fewer iterations and wall clock time\. This advantage is robust to initialization set size, exploration parameters, and inner\-solver choice\.

1McKetta Department of Chemical Engineering, The University of Texas at Austin, 200 E\. Dean Keaton St\. Stop C0400, Austin, Texas 78712, USA 2ExxonMobil Technology and Engineering, 22777 Springwoods Village Pkwy, Spring, Texas 77389, USA 3Energy Institute, The University of Texas at Austin, 2304 Whitis Ave\. Stop C2400, Austin, Texas 78712, USA 4Institute for Computational Engineering and Sciences, The University of Texas at Austin, 201 E\. 24th Street, POB 4\.102, Stop C0200, Austin, Texas 78712, USA

*Keywords*Bayesian optimization⋅\\cdotGrey\-box optimization⋅\\cdotBilevel optimization⋅\\cdotSurrogate\-based optimization⋅\\cdotMulti\-scale design⋅\\cdotDimensionality reduction

## 1Introduction

Surrogate\-based optimization is the standard approach for expensive black\-box functions\[[1](https://arxiv.org/html/2608.03045#bib.bib1)\], yet many problems contain known, differentiable substructure that monolithic surrogates waste samples learning\. When a model can be readily separated into an “expensive” black\-box component and a “cheap” white\-box component, building a single surrogate over all variables forces the surrogate model to learn both the unknown and the known components\. The cost of this redundancy grows with dimension: as the size of the combined variable space\[xWB,xBB\]\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]expands, the curse of dimensionality burdens the surrogate for structure it should never have had to discover\.

This situation arises naturally in multi\-scale engineering problems where design and operation are jointly optimized\. For example, a catalyst designer uses density functional theory \(DFT\) simulations to predict kinetic parameters, which are then used in well\-understood reactor mass and energy balances to evaluate reaction yield\. Similarly, a polymer engineer uses molecular simulations to estimate membrane transport properties, then solves the solution\-diffusion equations to size a separation unit\. In each case, the overall objective depends on two groups of decision variables\.*Black\-box variables*xBBx^\{\\mathrm\{BB\}\}are inputs to expensive simulations or experiments\. Conversely,*white\-box variables*xWBx^\{\\mathrm\{WB\}\}are variables present in the closed\-form macroscopic equation\. The black\-box and white\-box elements of the model are*separable*: the expensive computation depends only onxBBx^\{\\mathrm\{BB\}\}, producing outputsy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)that parameterize the white\-box model\. The key structural property is that, for any fixedxBBx^\{\\mathrm\{BB\}\}andyy, the remaining optimization overxWBx^\{\\mathrm\{WB\}\}is a standard NLP solvable by conventional methods assuming that theWB\\mathrm\{WB\}equations are continuous and thatxWBx^\{\\mathrm\{WB\}\}are all real\.

Several lines of work exploit known structure in expensive optimization, ranging from deterministic grey\-box solvers\[[2](https://arxiv.org/html/2608.03045#bib.bib2),[3](https://arxiv.org/html/2608.03045#bib.bib3),[4](https://arxiv.org/html/2608.03045#bib.bib4)\]to composite and structured Bayesian optimization\[[5](https://arxiv.org/html/2608.03045#bib.bib5),[6](https://arxiv.org/html/2608.03045#bib.bib6),[7](https://arxiv.org/html/2608.03045#bib.bib7)\]to data\-driven bilevel methods\[[8](https://arxiv.org/html/2608.03045#bib.bib8),[9](https://arxiv.org/html/2608.03045#bib.bib9)\]\. We review these in Section[3](https://arxiv.org/html/2608.03045#S3)\. The common gap is that no prior method combines a surrogate\-based global search over only the black\-box variables with exact NLP solution of the white\-box subproblem—a configuration that fully exploits variable separability for both dimensionality reduction and exact constraint satisfaction\.

We exploit this separability through a bilevel reformulation: the outer loop uses Bayesian optimization \(BO\) to search overxBBx^\{\\mathrm\{BB\}\}alone, while an inner optimizer solves the white\-box subproblem exactly for each candidate\. The Gaussian process \(GP\) surrogate therefore operates overℝnBB\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}rather thanℝnWB\+nBB\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}\}, and white\-box constraints are satisfied exactly whenever the inner optimizer converges to a feasible point\. The details of this reformulation are presented in Section[4](https://arxiv.org/html/2608.03045#S4)\.

We make three contributions:

1. 1\.A bilevel reformulation of separable grey\-box problems that reduces surrogate dimensionality fromnWB\+nBBn\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}tonBBn\_\{\\mathrm\{BB\}\}by solving the white\-box subproblem via global optimization, with white\-box constraints satisfied exactly whenever the inner optimizer converges to a feasible point—without penalty functions, chance constraints, or moment approximations\.
2. 2\.A benchmark suite of 13 separable grey\-box problems \(2–5 total variables, 0–3 constraints\) considering synthetic test functions and prototype engineering applications—the largest such suite for this problem class\.
3. 3\.Comprehensive empirical evidence from 8,450 independent optimization runs \(8,190 BO \+ 130 NLP \+ 130 BH\): the proposed bilevel strategy achieves 11×\\times–108×10^\{8\}\\timeslower regret than black\-box BO across all 13 problems, with equal or faster wall time and robustness to hyperparameters \(ninitn\_\{\\text\{init\}\},ξ\\xi\) and inner\-solver choice\.

## 2Problem Formulation

We consider optimization problems in which an expensive black\-box model is coupled with known, differentiable equations\. This setting arises in multi\-scale engineering design—molecular simulations feeding into process models—but also in any domain where part of the system is analytically tractable and part is not\. We first formalize the general problem, then identify the structural property our method exploits\.

### 2\.1Problem Setting and Notation

The decision variables comprise two groups:

- •xWB∈ℝnWBx^\{\\mathrm\{WB\}\}\\in\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\}:*white\-box variables*governed by known, differentiable equations \(e\.g\., temperatures, pressures, flow rates in engineering; policy parameters in control; design dimensions in structural optimization\)\. These variables enter the white\-box modelfWBf^\{\\mathrm\{WB\}\}, the objectiveJJ, and the inequality constraintsgg\.
- •xBB∈ℝnBBx^\{\\mathrm\{BB\}\}\\in\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}:*black\-box variables*that parameterize an expensive, non\-differentiable functionfBBf^\{\\mathrm\{BB\}\}\(e\.g\., molecular descriptors requiring DFT, material properties from simulations, hyperparameters of a costly simulator\)\. These variables enterfBBf^\{\\mathrm\{BB\}\}directly and may also appear inJJandgg\.

The black\-box model

y=fBB​\(xBB\),fBB:ℝnBB→ℝny,y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\),\\quad f^\{\\mathrm\{BB\}\}:\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}\\to\\mathbb\{R\}^\{n\_\{y\}\},\(1\)mapsxBBx^\{\\mathrm\{BB\}\}to a vector of intermediate parametersy∈ℝnyy\\in\\mathbb\{R\}^\{n\_\{y\}\}that couple the black\-box and white\-box components\. We assumefBBf^\{\\mathrm\{BB\}\}is expensive to evaluate and that we lack closed\-form expressions or derivative information for it\. Note thatyyrefers exclusively to the outputs of the black\-box function; the white\-box variablesxWBx^\{\\mathrm\{WB\}\}are not “outputs” in this sense but rather decision variables whose optimal values are determined by solving the white\-box subproblem for a givenyy\.

The white\-box model takes the implicit form

fWB​\(xWB,y\)=0,fWB:ℝnWB×ℝny→ℝnWB,f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},y\)=0,\\quad f^\{\\mathrm\{WB\}\}:\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\}\\times\\mathbb\{R\}^\{n\_\{y\}\}\\to\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\},\(2\)encoding known equations \(e\.g\., mass and energy balances\) that relate the white\-box variablesxWBx^\{\\mathrm\{WB\}\}to the black\-box outputsyy\. Althoughyymay appear as a parameter in the white\-box equations, it is fixed for any givenxBBx^\{\\mathrm\{BB\}\}; the white\-box optimization searches only overxWBx^\{\\mathrm\{WB\}\}\. We assumefWBf^\{\\mathrm\{WB\}\}is sufficiently smooth and that the Jacobian∂fWB/∂xWB\\partial f^\{\\mathrm\{WB\}\}/\\partial x^\{\\mathrm\{WB\}\}is available analytically\.

The inequality constraints

g​\(xWB,y,xBB\)≤0,g:ℝnWB×ℝny×ℝnBB→ℝng,g\(x^\{\\mathrm\{WB\}\},y,x^\{\\mathrm\{BB\}\}\)\\leq 0,\\quad g:\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\}\\times\\mathbb\{R\}^\{n\_\{y\}\}\\times\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}\\to\\mathbb\{R\}^\{n\_\{g\}\},\(3\)encodengn\_\{g\}design specifications, safety limits, or feasibility requirements that may depend on both variable groups and the black\-box outputs\.

The complete optimization problem over both variable groups is

minxWB,xBB\\displaystyle\\min\\limits\_\{x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\}\\quadJ​\(xWB,y,xBB\)\\displaystyle J\(x^\{\\mathrm\{WB\}\},y,x^\{\\mathrm\{BB\}\}\)\(4a\)s\.t\.fWB​\(xWB,y\)=0\\displaystyle f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},y\)=0\(4b\)y=fBB​\(xBB\)\\displaystyle y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)\(4c\)g​\(xWB,y,xBB\)≤0\\displaystyle g\(x^\{\\mathrm\{WB\}\},y,x^\{\\mathrm\{BB\}\}\)\\leq 0\(4d\)xWB∈𝒳WB,xBB∈𝒳BB\\displaystyle x^\{\\mathrm\{WB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\},\\quad x^\{\\mathrm\{BB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{BB\}\}\(4e\)
The objectiveJ​\(xWB,y,xBB\)J\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)and inequality constraintsg​\(xWB,y,xBB\)≤0g\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\\leq 0may depend on the white\-box variables, the black\-box outputs, and—in some problems—directly on the black\-box variablesxBBx^\{\\mathrm\{BB\}\}\.

## 3Related Work

Problem \([4](https://arxiv.org/html/2608.03045#S2.E4)\) sits at the intersection of surrogate\-based global optimization, grey\-box Bayesian optimization, bilevel programming, and constrained expensive optimization\. We organize prior work along two axes that follow directly from Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1): \(i\) whether the method distinguishesxBBx^\{\\mathrm\{BB\}\}fromxWBx^\{\\mathrm\{WB\}\}, and \(ii\) whether the white\-box subsystemfWBf^\{\\mathrm\{WB\}\}and constraintsggare exploited exactly or approximated\. Methods that collapse along axis \(i\) pay the curse of dimensionality innWB\+nBBn\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}; methods that collapse along axis \(ii\) trade exactness for the generality of uncertainty propagation\.

### 3\.1Surrogate\-Based Global Optimization

The use of surrogate models to guide expensive optimization traces to the Efficient Global Optimization \(EGO\) algorithm of Jones et al\.\[[10](https://arxiv.org/html/2608.03045#bib.bib10)\], which introduced Expected Improvement \(EI\) as an acquisition function for GP surrogates\. EGO treats the objective as a monolithic black box and builds a single GP over the full decision space𝒳WB×𝒳BB⊂ℝnWB\+nBB\\mathcal\{X\}^\{\\mathrm\{WB\}\}\\times\\mathcal\{X\}^\{\\mathrm\{BB\}\}\\subset\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}\}which is the default in most BO implementations\[[11](https://arxiv.org/html/2608.03045#bib.bib11),[1](https://arxiv.org/html/2608.03045#bib.bib1)\]\.

Within the deterministic surrogate realm, Boukouvala, Hasan, and Floudas\[[2](https://arxiv.org/html/2608.03045#bib.bib2)\]developed a methodology for constrained grey\-box problems that selects surrogates from a diverse set of functional forms and globally optimizes the resulting NLP\. The companion ARGONAUT framework\[[3](https://arxiv.org/html/2608.03045#bib.bib3)\]formalized this into an iterative algorithm with variable selection, bounds tightening, and constrained sampling for problems with up to 100 variables; Kieslich et al\.\[[12](https://arxiv.org/html/2608.03045#bib.bib12)\]extended the approach using Smolyak grids and polynomial surrogates\. These methods share our philosophy of exploiting known equations, but they build surrogates over the joint\(xWB,xBB\)\(x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\)space rather than restricting the surrogate toxBBx^\{\\mathrm\{BB\}\}as Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)permits;fWBf^\{\\mathrm\{WB\}\}is used only to generate training data or to tighten bounds on𝒳WB\\mathcal\{X\}^\{\\mathrm\{WB\}\}, not to eliminatexWBx^\{\\mathrm\{WB\}\}from the surrogate’s domain\. Moreover, the deterministic surrogates they employ lack the principled exploration–exploitation trade\-off of probabilistic models and tend to overexploit early samples under tight evaluation budgets\.

The trust\-region filter algorithms of Eason and Biegler\[[4](https://arxiv.org/html/2608.03045#bib.bib4),[13](https://arxiv.org/html/2608.03045#bib.bib13)\]offer a complementary local approach: they embed surrogates forfBBf^\{\\mathrm\{BB\}\}within a trust region and solve the white\-box NLP inxWBx^\{\\mathrm\{WB\}\}at each iteration—architecturally close to the inner subproblem in Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\), but with only local convergence guarantees inxBBx^\{\\mathrm\{BB\}\}rather than a global exploration mechanism\.

### 3\.2Grey\-Box Bayesian Optimization

Grey\-box BO methods incorporate known structure into the surrogate to reduce its modeling burden\. Astudillo and Frazier\[[5](https://arxiv.org/html/2608.03045#bib.bib5)\]introduced Bayesian Optimization of Composite Functions \(BOCF\), which in our notation modelsJ​\(xWB,fBB​\(xBB\),xBB\)J\(x^\{\\mathrm\{WB\}\},f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\),x^\{\\mathrm\{BB\}\}\)by training a multi\-output GP onfBBf^\{\\mathrm\{BB\}\}and propagating uncertainty throughJJvia sampling; later extensions handle function networks\[[14](https://arxiv.org/html/2608.03045#bib.bib14)\]and partial evaluations\[[15](https://arxiv.org/html/2608.03045#bib.bib15)\]\. Paulson and Lu\[[6](https://arxiv.org/html/2608.03045#bib.bib6)\]extend this to constrained grey\-box problems by combining multivariate GP models with a constrained expected utility function, using sample average approximation and chance constraints to propagate GP uncertainty through known equations\. González and Zavala\[[16](https://arxiv.org/html/2608.03045#bib.bib16)\]propose an interconnected\-systems BO framework with adaptive linearization, Winz et al\.\[[17](https://arxiv.org/html/2608.03045#bib.bib17)\]modify the upper confidence bound to exploit grey\-box derivative information, and Kudva and Paulson\[[7](https://arxiv.org/html/2608.03045#bib.bib7)\]exploit function\-network topology for robust optimization\.

All of these methods place a surrogate on the intermediate outputsyydefined in Equation \([1](https://arxiv.org/html/2608.03045#S2.E1)\) and propagate uncertainty through Equation \([4b](https://arxiv.org/html/2608.03045#S2.E4.2)\), requiring multi\-output GPs and moment approximations or sampling\. When Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\) holds—i\.e\., the white\-box subproblem is solvable to global optimality—using a surrogate model foryyadds modeling complexity without clear benefit\. We instead place a scalar GP on the value function

xBB↦J​\(xWB⁣⋆​\(xBB\),fBB​\(xBB\),xBB\)x^\{\\mathrm\{BB\}\}\\;\\mapsto\\;J\\bigl\(x^\{\\mathrm\{WB\}\\star\}\(x^\{\\mathrm\{BB\}\}\),\\;f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\),\\;x^\{\\mathrm\{BB\}\}\\bigr\)\(5\)inℝnBB\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}, avoiding both multi\-output surrogates and uncertainty propagation \(Table[1](https://arxiv.org/html/2608.03045#S4.T1)\)\.

### 3\.3Bilevel Optimization

Bilevel programming has a rich history in the global optimization literature\. Gümüş and Floudas\[[18](https://arxiv.org/html/2608.03045#bib.bib18)\]developed the first rigorous global optimization approach for nonlinear bilevel problems using theα\\alphaBB framework\. Mitsos, Chachuat, and Barton\[[19](https://arxiv.org/html/2608.03045#bib.bib19)\]extended this to bilevel dynamic optimization, Kleniati and Adjiman\[[20](https://arxiv.org/html/2608.03045#bib.bib20)\]addressed the mixed\-integer case, and Faísca et al\.\[[21](https://arxiv.org/html/2608.03045#bib.bib21)\]proposed parametric global optimization for bilevel programs via multi\-parametric programming\. In the language of Section[2](https://arxiv.org/html/2608.03045#S2), these methods treatxBBx^\{\\mathrm\{BB\}\}as the upper\-level decision andxWBx^\{\\mathrm\{WB\}\}as the lower\-level decision, but assume analytical access to bothfBBf^\{\\mathrm\{BB\}\}andfWBf^\{\\mathrm\{WB\}\}—an assumption that breaks whenfBBf^\{\\mathrm\{BB\}\}is an expensive black\-box simulation\.

Data\-driven approaches address this gap\. The DOMINO framework of Beykal et al\.\[[8](https://arxiv.org/html/2608.03045#bib.bib8)\]solves bilevel mixed\-integer nonlinear problems by sampling the upper\-level objective and solving the lower\-level problem to global optimality at each sample point, using a grey\-box optimization solver \(ARGONAUT/NOMAD\) for the upper\-level search\. DOMINO is the closest non\-BO predecessor to our approach: it solves thexWBx^\{\\mathrm\{WB\}\}subproblem to global optimality at each sampledxBBx^\{\\mathrm\{BB\}\}, matching Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\), but uses deterministic surrogates overxBBx^\{\\mathrm\{BB\}\}\. We instead use a GP with a principled acquisition function that balances exploration and exploitation—an advantage when the evaluation budget is small \(<<250 calls\) and the risk of overexploiting early samples is high\.

Within the BO literature, Kieffer et al\.\[[9](https://arxiv.org/html/2608.03045#bib.bib9)\]first applied BO to bilevel problems, using EI at the outer level and Sequential Least Squares Programming \(SLSQP\) at the inner level\. This nested architecture matches our outer\-over\-xBBx^\{\\mathrm\{BB\}\}/inner\-over\-xWBx^\{\\mathrm\{WB\}\}decomposition, but Kieffer et al\. did not exploit the separability of Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)to reduce surrogate dimensionality fromnWB\+nBBn\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}tonBBn\_\{\\mathrm\{BB\}\}\. Recent work addresses the harder setting where*both*levels are black boxes: Ekmekcioglu et al\.\[[22](https://arxiv.org/html/2608.03045#bib.bib22)\]model both levels as GPs over the joint space; Chew et al\.\[[23](https://arxiv.org/html/2608.03045#bib.bib23)\]construct confidence\-bound trusted sets with regret guarantees; Fu et al\.\[[24](https://arxiv.org/html/2608.03045#bib.bib24)\]provide convergence analysis for BO bilevel problems where the inner loop is unconstrained Stochastic Gradient Descent \(SGD\)\. When the inner problem is a known white\-box NLP, solving it exactly is both cheaper and more accurate than learning a second surrogate\.

Baldea\[[25](https://arxiv.org/html/2608.03045#bib.bib25)\]introduced the multi\-scale bilevel BO framework we extend here, demonstrating it on a single chemical process case study\. We advance that work with a 13\-problem benchmark suite, with systematic hyperparameter sweeps, and detailed comparisons quantifying when and why bilevel decomposition helps\.

### 3\.4Constraint Handling in Expensive Optimization

Standard constrained BO methods treat each component ofg​\(xWB,y,xBB\)g\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)as an independent expensive black box with its own GP\. Gardner et al\.\[[26](https://arxiv.org/html/2608.03045#bib.bib26)\]introduced weighted expected improvement for unknown constraints; Gelbart et al\.\[[27](https://arxiv.org/html/2608.03045#bib.bib27)\]extended this with probability\-of\-feasibility weighting; Picheny et al\.\[[28](https://arxiv.org/html/2608.03045#bib.bib28)\]proposed a slack\-variable augmented Lagrangian approach\. These methods are appropriate when constraint functions are truly unknown, but when components ofggdepend onxWBx^\{\\mathrm\{WB\}\}andyythrough known equations—mass balances, safety envelopes, thermodynamic limits—creating surrogates increases computational cost and introduces approximation error\.

Our inner global optimizer enforces the white\-box components ofggexactly via Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\), avoiding penalty functions, chance constraints, and moment approximations\. The advantage is largest when feasible regions are tight: in such cases\(reflected in the Distillation and Heat\-Exchanger problems in Section[6](https://arxiv.org/html/2608.03045#S6)\), bilevel BO achieves\>\>100,000×\\timeslower regret than penalty\-based black\-box BO\.

## 4Method: Bilevel Bayesian Optimization

Problem \([4](https://arxiv.org/html/2608.03045#S2.E4)\) possesses a structural property that we formalize as an assumption and then exploit algorithmically\.

###### Assumption 1\(Separable grey\-box structure\)\.

Problem \([4](https://arxiv.org/html/2608.03045#S2.E4)\) satisfies:

1. \(i\)The black\-box functionfBBf^\{\\mathrm\{BB\}\}depends only onxBBx^\{\\mathrm\{BB\}\}\. For any fixedxBBx^\{\\mathrm\{BB\}\}andy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\), the resulting optimization overxWBx^\{\\mathrm\{WB\}\}alone is a well\-posed NLP\.
2. \(ii\)For any fixedxBBx^\{\\mathrm\{BB\}\}andy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\), the white\-box subproblem minxWB\\displaystyle\\min\\limits\_\{x^\{\\mathrm\{WB\}\}\}\\;J​\(xWB,y,xBB\)\\displaystyle J\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\(6a\)s\.t\.fWB​\(xWB,y\)=0\\displaystyle f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},y\)=0\(6b\)g​\(xWB,y,xBB\)≤0\\displaystyle g\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\\leq 0\(6c\)xWB∈𝒳WB\\displaystyle x^\{\\mathrm\{WB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\}\(6d\)is solvable to global optimality \(of the subproblem\) by a suitable global optimizer\.

Under Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1), the optimization can be decomposed as follows: for any fixedxBBx^\{\\mathrm\{BB\}\}, the black\-box evaluationy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)and the inner optimization overxWBx^\{\\mathrm\{WB\}\}are self\-contained\. This decoupling is the structural feature we exploit: rather than optimizing all variables jointly, we can search over the low\-dimensional black\-box space while solving the white\-box subproblem exactly at each step\. Note thatfBBf^\{\\mathrm\{BB\}\}can always be augmented with identity mappings \(e\.g\.,yny\+1=x1BBy\_\{n\_\{y\}\+1\}=x^\{\\mathrm\{BB\}\}\_\{1\}\) to absorb any directxBBx^\{\\mathrm\{BB\}\}dependence inJJorgg, so the distinction between “depends onxBBx^\{\\mathrm\{BB\}\}throughyy” and “depends onxBBx^\{\\mathrm\{BB\}\}directly” is a modeling choice, not a structural limitation\.

### 4\.1Bilevel Reformulation

Based on Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1), we reformulate \([4](https://arxiv.org/html/2608.03045#S2.E4)\) as a bi\-level optimization problem that separates what we know from what we must learn:

minxBB\\displaystyle\\min\\limits\_\{x^\{\\mathrm\{BB\}\}\}\\quadJ​\(xWB⁣⋆,y,xBB\)\\displaystyle J\(x^\{\\mathrm\{WB\}\\star\},y,x^\{\\mathrm\{BB\}\}\)\(7a\)s\.t\.y=fBB​\(xBB\)\\displaystyle y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)\(7b\)xWB⁣⋆=arg⁡minxWB⁡J​\(xWB,y,xBB\)\\displaystyle x^\{\\mathrm\{WB\}\\star\}=\\arg\\min\\limits\_\{x^\{\\mathrm\{WB\}\}\}\\;J\(x^\{\\mathrm\{WB\}\},y,x^\{\\mathrm\{BB\}\}\)\(7c\)s\.t\.fWB​\(xWB,y\)=0\\displaystyle\\quad\\quad\\text\{s\.t\.\}\\quad f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},y\)=0\(7d\)g​\(xWB,y,xBB\)≤0\\displaystyle\\quad\\quad\\quad\\quad\\;g\(x^\{\\mathrm\{WB\}\},y,x^\{\\mathrm\{BB\}\}\)\\leq 0\(7e\)xWB∈𝒳WB\\displaystyle\\quad\\quad\\quad\\quad\\;x^\{\\mathrm\{WB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\}\(7f\)xBB∈𝒳BB\\displaystyle x^\{\\mathrm\{BB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{BB\}\}\(7g\)
Theinner problem\([7c](https://arxiv.org/html/2608.03045#S4.E7.3)\)–\([7f](https://arxiv.org/html/2608.03045#S4.E7.6)\) is a standard nonlinear program \(NLP\)\. Given fixedxBBx^\{\\mathrm\{BB\}\}andy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\), it optimizes only the white\-box variablesxWBx^\{\\mathrm\{WB\}\}using the white box model;xBBx^\{\\mathrm\{BB\}\}andyyenter as parameters\. Because the white\-box equations are closed\-form and inexpensive to evaluate, this subproblem can be solved to global optimality using a global optimizer such as Basin\-Hopping \(BH\) with an SLSQP local minimizer\.

Theouter problem\([7a](https://arxiv.org/html/2608.03045#S4.E7.1)\)–\([7b](https://arxiv.org/html/2608.03045#S4.E7.2)\) searches over black\-box decisionsxBBx^\{\\mathrm\{BB\}\}\. For each candidatexBBx^\{\\mathrm\{BB\}\}, we evaluate the black box to obtainyy, solve the inner problem to global optimality to getxWB⁣⋆x^\{\\mathrm\{WB\}\\star\}, and record the resulting objective valueJ​\(xWB⁣⋆,y,xBB\)J\(x^\{\\mathrm\{WB\}\\star\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\. This outer loop is where we apply Bayesian optimization\. In practice, the number of outer\-loop iterations—each requiring one evaluation offBBf^\{\\mathrm\{BB\}\}—is the primary driver of total campaign time, since each black\-box evaluation may involve hours of computation or physical experimentation\. The inner NLP solve adds negligible overhead by comparison \(see Remark 4 and Section[7](https://arxiv.org/html/2608.03045#S7)\)\.

### 4\.2Algorithm, Surrogate Model, and Acquisition Function

Algorithm 1Bilevel Bayesian Optimization1:Black\-box bounds

𝒳BB\\mathcal\{X\}^\{\\mathrm\{BB\}\}, white\-box bounds

𝒳WB\\mathcal\{X\}^\{\\mathrm\{WB\}\}, models

fBBf^\{\\mathrm\{BB\}\},

fWBf^\{\\mathrm\{WB\}\},

JJ,

gg; initial samples

ninitn\_\{\\text\{init\}\}, iterations

NN, exploration parameter

ξ\\xi
2:Best

xBB⁣⋆x^\{\\mathrm\{BB\}\\star\},

xWB⁣⋆x^\{\\mathrm\{WB\}\\star\},

J⋆J^\{\\star\}
3:// Initialization

4:Sample

\{xiBB\}i=1ninit\\\{x^\{\\mathrm\{BB\}\}\_\{i\}\\\}\_\{i=1\}^\{n\_\{\\text\{init\}\}\}via Latin hypercube in

𝒳BB\\mathcal\{X\}^\{\\mathrm\{BB\}\}
5:for

i=1,…,niniti=1,\\ldots,n\_\{\\text\{init\}\}do

6:Evaluate

yi=fBB​\(xiBB\)y\_\{i\}=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\_\{i\}\)
7:Solve inner problem \(BH\):

xiWB⁣⋆=arg⁡minxWB⁡J​\(xWB,yi,xiBB\)x^\{\\mathrm\{WB\}\\star\}\_\{i\}=\\arg\\min\_\{x^\{\\mathrm\{WB\}\}\}J\(x^\{\\mathrm\{WB\}\},y\_\{i\},x^\{\\mathrm\{BB\}\}\_\{i\}\)s\.t\.

fWB​\(xWB,yi\)=0f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},y\_\{i\}\)=0,

g​\(xWB,yi,xiBB\)≤0g\(x^\{\\mathrm\{WB\}\},y\_\{i\},x^\{\\mathrm\{BB\}\}\_\{i\}\)\\leq 0
8:Record

Ji=J​\(xiWB⁣⋆\)J\_\{i\}=J\(x^\{\\mathrm\{WB\}\\star\}\_\{i\}\)\(or penalty if infeasible\)

9:endfor

10:// BO loop

11:for

t=ninit\+1,…,ninit\+Nt=n\_\{\\text\{init\}\}\+1,\\ldots,n\_\{\\text\{init\}\}\+Ndo

12:Fit GP surrogate

𝒢​𝒫\\mathcal\{GP\}to

\{\(xiBB,Ji\)\}i=1t−1\\\{\(x^\{\\mathrm\{BB\}\}\_\{i\},J\_\{i\}\)\\\}\_\{i=1\}^\{t\-1\}
13:Select

xtBB=arg⁡maxxBB∈𝒳BB⁡αEI​\(xBB;𝒢​𝒫,ξ\)x^\{\\mathrm\{BB\}\}\_\{t\}=\\arg\\max\_\{x^\{\\mathrm\{BB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{BB\}\}\}\\alpha\_\{\\text\{EI\}\}\(x^\{\\mathrm\{BB\}\};\\mathcal\{GP\},\\xi\)
14:Evaluate

yt=fBB​\(xtBB\)y\_\{t\}=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\_\{t\}\)
15:Solve inner problem \(BH\)

→xtWB⁣⋆\\to x^\{\\mathrm\{WB\}\\star\}\_\{t\}, record

JtJ\_\{t\}
16:Update dataset

17:endfor

18:returnbest

\(xBB⁣⋆,xWB⁣⋆,J⋆\)\(x^\{\\mathrm\{BB\}\\star\},x^\{\\mathrm\{WB\}\\star\},J^\{\\star\}\)

Given some dataset𝒟=\{xiB​B,Ji∗\}​i=1​…​n\\mathcal\{D\}=\\\{x\_\{i\}^\{BB\},J\_\{i\}^\{\*\}\\\}\\,i=1\\dots n, create a surrogate modelJ^=𝒢​𝒫​\(xB​B\)\\hat\{J\}=\\mathcal\{GP\}\(x^\{BB\}\)Estimatexn\+1B​Bx\_\{n\+1\}^\{BB\}Outer Problem \(BO\)Black\-box Functionyn\+1=fB​B​\(xn\+1B​B\)y\_\{n\+1\}=f^\{BB\}\(x\_\{n\+1\}^\{BB\}\)Solve givenxn\+1B​Bx\_\{n\+1\}^\{BB\},yn\+1y\_\{n\+1\}:xn\+1W​B⁣∗=arg⁡minxW​B⁡J​\(xW​B,yn\+1,xn\+1B​B\)s\.t\.fW​B​\(xW​B,yn\+1\)=0g​\(xW​B,yn\+1,xn\+1B​B\)≤0\\begin\{aligned\} x\_\{n\+1\}^\{WB\*\}=\\arg\\min\_\{x^\{WB\}\}J\(x^\{WB\},y\_\{n\+1\},x\_\{n\+1\}^\{BB\}\)\\\\ \\text\{s\.t\.\}\\quad f^\{WB\}\(x^\{WB\},y\_\{n\+1\}\)&=0\\\\ g\(x^\{WB\},y\_\{n\+1\},x\_\{n\+1\}^\{BB\}\)&\\leq 0\\end\{aligned\}Inner Problem \(BH\)xn\+1B​Bx\_\{n\+1\}^\{BB\}yn\+1y\_\{n\+1\}Jn\+1∗​\(xn\+1W​B⁣∗\)J\_\{n\+1\}^\{\*\}\(x\_\{n\+1\}^\{WB\*\}\)Figure 1:Flow of the multi\-scale Bayesian optimization framework\. The outer loop uses BO to select black\-box variablesxBBx^\{\\mathrm\{BB\}\}, which are evaluated through the black\-box model to produce parametersyy\. The inner loop solves the white\-box optimization overxWBx^\{\\mathrm\{WB\}\}via Basin\-Hopping givenxBBx^\{\\mathrm\{BB\}\}andyy, returning the optimal objectiveJ​\(xWB⁣⋆,y,xBB\)J\(x^\{\\mathrm\{WB\}\\star\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)to update the BO surrogate\.We describe the proposed approach in Algorithm[1](https://arxiv.org/html/2608.03045#alg1)and Figure[1](https://arxiv.org/html/2608.03045#S4.F1)and provide a discussion as follows: The GP surrogate in Algorithm[1](https://arxiv.org/html/2608.03045#alg1)models the scalar mapping \([5](https://arxiv.org/html/2608.03045#S3.E5)\)—a function ofxBBx^\{\\mathrm\{BB\}\}alone\. In our implementation, we use a Gaussian process with an automatic relevance determination \(ARD\) radial basis function \(RBF\) kernel, where a separate length scale per dimension allows the GP to adapt to anisotropic landscapes\. The outer loop selects the next evaluation by maximizing Expected Improvement \(EI\) with exploration parameterξ≥0\\xi\\geq 0\. Kernel and acquisition function details for our implementation are in Appendices[C](https://arxiv.org/html/2608.03045#A3)and[B](https://arxiv.org/html/2608.03045#A2); here we emphasize that both are standard—the gains reported in Section[6](https://arxiv.org/html/2608.03045#S6)come from problem structure, not from novel surrogate or acquisition choices, although we expect that adopting them would yield further benefits\.

### 4\.3Key Properties

The bilevel reformulation \([7](https://arxiv.org/html/2608.03045#S4.E7)\) has several properties worth highlighting\. The first is formalized as a proposition\.

###### Proposition 1\.

Under Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1), the bilevel reformulation \([7](https://arxiv.org/html/2608.03045#S4.E7)\) preserves the global optimum of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) and reduces the GP surrogate domain fromℝnWB\+nBB\\mathbb\{R\}^\{n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}\}toℝnBB\\mathbb\{R\}^\{n\_\{\\mathrm\{BB\}\}\}\.

The proof is in Appendix[G](https://arxiv.org/html/2608.03045#A7)\. The key consequences are threefold:

*Exact constraint satisfaction\.*Constraintsg​\(xWB,y,xBB\)≤0g\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\\leq 0are handled by the inner global optimizer, not by penalty functions or chance constraints\. The GP never needs to learn the constraint boundary\. When the inner optimizer converges to a globally feasible point \(Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\)\), feasibility is guaranteed\. In our implementation, Basin\-Hopping \(BH\)\[[29](https://arxiv.org/html/2608.03045#bib.bib29),[30](https://arxiv.org/html/2608.03045#bib.bib30),[31](https://arxiv.org/html/2608.03045#bib.bib31),[32](https://arxiv.org/html/2608.03045#bib.bib32),[33](https://arxiv.org/html/2608.03045#bib.bib33)\]with Sequential Least Squares Programming \(SLSQP\)\[[34](https://arxiv.org/html/2608.03045#bib.bib34),[35](https://arxiv.org/html/2608.03045#bib.bib35),[36](https://arxiv.org/html/2608.03045#bib.bib36)\]local minimization provides high\-confidence global solutions\. We also compare a simpler solution using SLSQP with multiple restarts as a more time\-efficient method to solve the reduced\-dimension inner problem\. Further guarantees of global optimization could be obtained with deterministic global solvers \(e\.g\., BARON\), as discussed in Section[7](https://arxiv.org/html/2608.03045#S7)\. Details on these solvers and implementation are available in Appendix[E](https://arxiv.org/html/2608.03045#A5)\.

*Scalar surrogate\.*The method requires only one GP for the scalar optimal objectiveJ​\(xWB⁣⋆,y,xBB\)J\(x^\{\\mathrm\{WB\}\\star\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\. Approaches that surrogate the black\-box outputsy∈ℝnyy\\in\\mathbb\{R\}^\{n\_\{y\}\}instead require a multi\-output GP, with the associated challenges of cross\-covariance specification andO​\(ny3\)O\(n\_\{y\}^\{3\}\)scaling\.

Additional implementation details—infeasibility handling and solver modularity—are discussed in Appendix[E](https://arxiv.org/html/2608.03045#A5)\.

### 4\.4Relationship to Existing Grey\-Box Methods

Table 1:Comparison of grey\-box BO approaches\. Bilevel BO \(this work\) is the only method that both separates variables and solves the white\-box subproblem exactly\.Table[1](https://arxiv.org/html/2608.03045#S4.T1)summarizes the distinctions between the proposed method and a subset of the methods reviewed in Section[3](https://arxiv.org/html/2608.03045#S3)that are most directly related and relevant\. Both bilevel BO and COBALT achieve the same input dimensionality reduction tonBBn\_\{\\mathrm\{BB\}\}\. The distinction is in the surrogate output \(scalar vs\. multi\-output GP\) and constraint handling \(exact global optimization vs\. moment approximation\), not the input dimension\. COBALT and BOCF surrogate intermediate black\-box outputsyyand propagate uncertainty through the white\-box equations, requiring multi\-output GPs and moment approximations\. Our method surrogates the scalar optimal objectiveJ​\(xWB⁣⋆\)J\(x^\{\\mathrm\{WB\}\\star\}\)directly, at the cost of an inner global optimization solve per outer iteration\. When the white\-box model is moderately sized and the black\-box evaluation is the primary computational expense, this overhead is negligible\. The approaches are complementary: COBALT handlesyy\-dependent constraints natively and enables risk\-averse decisions via uncertainty propagation; BOCF handles non\-separable composite functions\. Our method trades their generality for simplicity and exactness when Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)holds\.

## 5Benchmark Problem Suite

We propose a suite of 13 problems spanning synthetic test functions and engineering domains, summarized in Table[2](https://arxiv.org/html/2608.03045#S5.T2)\. Full mathematical formulations are in Appendix[A](https://arxiv.org/html/2608.03045#A1)\. The suite is designed around five principles: \(1\) every problem has a verified global optimum \(Appendix[F](https://arxiv.org/html/2608.03045#A6)\); \(2\) constraints are physically motivated; \(3\) black\-box functions encode realistic structure\-property relationships \(Sabatier volcano curves, Langmuir isotherms, Robeson upper bounds, Hansen solubility\); \(4\) problems span a range of difficulty \(unconstrained to 3 constraints, 2D to 5D, min and max\); and \(5\) all problems are cheap to evaluate, enabling the 8,450\-run statistical comparison in Section[6](https://arxiv.org/html/2608.03045#S6)\.

To our knowledge, no existing benchmark suite targets a separable grey\-box structure that reflects Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1): an outer black\-box mappingy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)coupled with an inner white\-box optimization overxWBx^\{\\mathrm\{WB\}\}, subject to inequality constraints\. Existing problem suites for composite BO\[[5](https://arxiv.org/html/2608.03045#bib.bib5)\]do not separate variables; those used in grey\-box BO\[[6](https://arxiv.org/html/2608.03045#bib.bib6)\]contain only a handful of problems; and standard global optimization testbeds such as the IEEE Congress on Evolutionary Computation \(CEC\) benchmarks\[[37](https://arxiv.org/html/2608.03045#bib.bib37)\]and the Comparing Continuous Optimizers \(COCO\) platform with its Black\-Box Optimization Benchmarking \(BBOB\) suite\[[38](https://arxiv.org/html/2608.03045#bib.bib38)\]lack bilevel structure entirely\.

Table 2:Summary of the 13\-problem benchmark suite\.nWBn\_\{\\mathrm\{WB\}\}: white\-box variables,nBBn\_\{\\mathrm\{BB\}\}: black\-box variables,nyn\_\{y\}: black\-box outputs,ngn\_\{g\}: inequality constraints\.ProblemnWBn\_\{\\mathrm\{WB\}\}nBBn\_\{\\mathrm\{BB\}\}nyn\_\{y\}ngn\_\{g\}TypeDomainSmall\-Feasible\-Region 11111minSyntheticSmall\-Feasible\-Region 21111minSyntheticRastrigin2110minMultimodalToy\-Hydrology1112minHydrologyRosen\-Suzuki2223minConstrainedCSTR3232maxCatalysisHeat\-Exchanger3232minHeat Exchanger NetworkPSA3232minAdsorptionBatch\-Reactor3232minKineticsDistillation3233minSeparationEvaporator3232minEvaporationMembrane3232minMembrane sep\.Williams\-Otto3222maxProcess opt\.

### 5\.1Representative Problems

We highlight two problems; full formulations for all 13 appear in Appendix[A](https://arxiv.org/html/2608.03045#A1)\. Adapted from Gardner et al\.\[[26](https://arxiv.org/html/2608.03045#bib.bib26)\]and Ariafar et al\.\[[39](https://arxiv.org/html/2608.03045#bib.bib39)\], these 2D constrained minimization problems feature small, disconnected feasible regions\. The two\-dimensional structure permits direct visualization of the search behavior \(Section[6\.2](https://arxiv.org/html/2608.03045#S6.SS2)\)\.

## 6Computational Experiments

We evaluate bilevel BO against full\-space baselines on the 13\-problem benchmark suite\. Section[6\.1](https://arxiv.org/html/2608.03045#S6.SS1)specifies the solvers, metrics, and hyperparameter sweep\. Section[6\.2](https://arxiv.org/html/2608.03045#S6.SS2)then uses the two Small\-Feasible\-Region problems to visualize how bilevel and black\-box BO explore the search space differently\. Section[6\.3](https://arxiv.org/html/2608.03045#S6.SS3)reports final regret across all 13 benchmarks, convergence speed, and robustness to hyperparameter and dimension\.

### 6\.1Experimental Design

This section specifies the four solvers compared \(Section[6\.1\.1](https://arxiv.org/html/2608.03045#S6.SS1.SSS1)\) and the four metrics reported \(Section[6\.1\.2](https://arxiv.org/html/2608.03045#S6.SS1.SSS2)\)\. Per\-problem hyperparameter settings and full sweep results are deferred to Appendix[J](https://arxiv.org/html/2608.03045#A10)\.

#### 6\.1\.1Solvers Compared

We compare four solvers that represent different ways of handling the grey\-box structure in problem \([4](https://arxiv.org/html/2608.03045#S2.E4)\)\.

The first solver,black\-box BH, applies SciPy’s Basin\-Hopping \(BH\) solver over the full variable space\[xWB,xBB\]\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]with an SLSQP local minimizer\. Basin\-Hopping performs a global search by alternating random perturbations \(uniform steps clipped to variable bounds\) with local SLSQP minimizations, accepting or rejecting each local minimum via a Metropolis criterion \(T=100T=100\)\. This solver provides a global\-search baseline: Basin\-Hopping is gradient\-free at the global level and explores the full feasible region\. However, it evaluatesfBBf^\{\\mathrm\{BB\}\}at every local minimization, so the total number of black\-box evaluations grows with the number of Basin\-Hopping iterations\.

The second solver,black\-box BO \(BB\-BO\), applies standard Bayesian optimization over the full\(nWB\+nBB\)\(n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}\)\-dimensional space\. A single GP surrogate models the mapping\[xWB,xBB\]↦J\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]\\mapsto J\. Constraints are handled via a penalty function\. This solver ignores the separable structure entirely and serves as the BO baseline\.

The third solver,bilevel BO with multi\-start SLSQP \(Bi\-BO \(SLSQP\)\), is the proposed method with a gradient\-based inner solver\. It places the GP surrogate over only thenBBn\_\{\\mathrm\{BB\}\}\-dimensional black\-box space and solves the constrained white\-box subproblem using 50 Latin hypercube restarts of SLSQP, returning the best feasible solution found, as described in Algorithm[1](https://arxiv.org/html/2608.03045#alg1)\.

The fourth solver,bilevel BO with Basin\-Hopping \(Bi\-BO \(BH\)\), is identical to the third except that it uses Basin\-Hopping \(100 iterations, SLSQP local minimizer, Metropolis temperatureT=100T=100\) as the inner\-loop global solver\. Both bilevel variants share the same outer\-loop GP surrogate and EI acquisition function; they differ only in the inner\-loop solver used for the white\-box subproblem\.

All four solvers use Expected Improvement \(EI\) as the acquisition function for the BO variants\. All three BO solvers share the same ARD RBF kernel and GP configuration \(Appendix[C](https://arxiv.org/html/2608.03045#A3)\)\. We sweep over initial sample sizesninit∈\{5,25,50\}n\_\{\\text\{init\}\}\\in\\\{5,25,50\\\}and exploration parametersξ∈\{0\.001,…,1\.0\}\\xi\\in\\\{0\.001,\\ldots,1\.0\\\}, yielding 8,450 independent runs across all 13 problems \(8,190 BO \+ 130 NLP \+ 130 BH; Appendix[J](https://arxiv.org/html/2608.03045#A10)\)\.

#### 6\.1\.2Metrics

We report four metrics\.*Simple regret*rt=\|Jbest​\(t\)−J⋆\|r\_\{t\}=\|J\_\{\\text\{best\}\}\(t\)\-J^\{\\star\}\|measures the gap between the best objective found afterttevaluations and the verified global optimumJ⋆J^\{\\star\}\.*Convergence speed*is the number of iterations needed to reduce regret to 1% of its initial value\.*Wall time*is the total elapsed time forninit\+200n\_\{\\text\{init\}\}\+200evaluations\. Finally, we count the total number of*black\-box evaluations*\(calls tofBBf^\{\\mathrm\{BB\}\}\) to assess sample efficiency\.

### 6\.2Mechanism Visualization

We begin with the two Small\-Feasible\-Region problems because their two\-dimensional structure permits direct visualization of how bilevel BO searches differently from black\-box BO\. Both problems have one white\-box variable \(xWBx^\{\\mathrm\{WB\}\}\) and one black\-box variable \(xBBx^\{\\mathrm\{BB\}\}\), with small, disconnected feasible regions defined by nonlinear constraints\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x1.png)Figure 2:Small\-Feasible\-Region 1: search behavior of black\-box BO \(left\) and bilevel BO \(right\) after 50 initial samples \+ 200 BO iterations\. The colored background shows the objective valueJ​\(xWB,y\)J\(x^\{\\mathrm\{WB\}\},y\)where feasible \(yellow = low, purple = high\)\. White regions are*infeasible*: no combination of\(xWB,xBB\)\(x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\)in those areas satisfies the constraints\. The hatched islands are the only feasible regions, and the⋆\\starmarks the global optimum within them\. Black\-box BO scatters queries across the full\(xWB,xBB\)\(x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\)space, wasting many evaluations in infeasible white regions\. Bilevel BO operates along thexBBx^\{\\mathrm\{BB\}\}axis only—each query selects anxBBx^\{\\mathrm\{BB\}\}value, and the inner BH solver finds the best feasiblexWBx^\{\\mathrm\{WB\}\}, concentrating evaluations in or near the feasible islands\.Figure[2](https://arxiv.org/html/2608.03045#S6.F2)shows the key difference\. In SFR\-1, the black\-box function is the identity \(y=xBBy=x^\{\\mathrm\{BB\}\}\), so the bilevel decomposition reduces to: given a candidateyy\-value from the outer loop, find thexWBx^\{\\mathrm\{WB\}\}that minimizessin⁡\(xWB\)\+y\\sin\(x^\{\\mathrm\{WB\}\}\)\+ywithin the feasible set\. The inner BH solver handles constraint satisfaction exactly, so every bilevel BO query lands in a feasible region\. Black\-box BO, by contrast, must jointly learn the objective function landscape*and*the constraint boundary in two dimensions, wasting many samples in infeasible regions\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x2.png)Figure 3:Small\-Feasible\-Region 2: same layout as Figure[2](https://arxiv.org/html/2608.03045#S6.F2)\. Here the black\-box function is nonlinear \(y=xBB−5\.52​exp⁡\(xBB/2\)y=\\frac\{x^\{\\mathrm\{BB\}\}\-5\.5\}\{2\}\\exp\(x^\{\\mathrm\{BB\}\}/2\)\), distorting the relationship betweenxBBx^\{\\mathrm\{BB\}\}and the feasible set\. Bilevel BO still operates near the optimum because the GP surrogate learns the scalar mappingxBB↦J​\(xWB⁣⋆\)x^\{\\mathrm\{BB\}\}\\mapsto J\(x^\{\\mathrm\{WB\}\\star\}\)in one dimension\.SFR\-2 \(Figure[3](https://arxiv.org/html/2608.03045#S6.F3)\) introduces a nonlinear black\-box function, making the mapping fromxBBx^\{\\mathrm\{BB\}\}to the feasible objective landscape more complex\. The same mechanism applies: the outer loop proposesxBBx^\{\\mathrm\{BB\}\}, the black box function is used to computeyy, and the inner BH solver finds the best feasiblexWBx^\{\\mathrm\{WB\}\}\. The GP surrogate now models a one\-dimensional functionxBB↦J​\(xWB⁣⋆\)x^\{\\mathrm\{BB\}\}\\mapsto J\(x^\{\\mathrm\{WB\}\\star\}\)rather than the two\-dimensional objective\-plus\-constraint landscape that black\-box BO must learn\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x3.png)Figure 4:Regret convergence for the two Small\-Feasible\-Region problems\. Each curve usesninit=50n\_\{\\text\{init\}\}=50initial Latin hypercube samples and the EI exploration parameterξ\\xithat yielded the lowest final regret\. The dashed line marks the end of initialization\. Solid lines show the mean over 10 independent repetitions; shaded bands show±1\\pm 1standard deviation\. Bilevel BO \(red\) converges 1–2 orders of magnitude below black\-box BO \(blue\) on both problems\.The convergence curves \(Figure[4](https://arxiv.org/html/2608.03045#S6.F4)\) confirm the visual impression\. On SFR\-1, Bi\-BO \(SLSQP\) reduces regret by an order of magnitude relative to black\-box BO—an 11×\\timesimprovement\. On SFR\-2, Bi\-BO \(SLSQP\) reaches the global optimum \(regret<10−5<10^\{\-5\}\) while black\-box BO stalls near10−1\.610^\{\-1\.6\}—a 7,872×\\timesimprovement\. In both cases, the mechanism is the same: bilevel BO needs only to learn a one\-dimensional function, while black\-box BO must simultaneously model a two\-dimensional objective and discover the small feasible islands\.

### 6\.3Main Results

The pattern observed on the SFR problems extends to all 13 benchmarks\. Table[3](https://arxiv.org/html/2608.03045#S6.T3)reports final regret after 200 BO iterations; Figure[5](https://arxiv.org/html/2608.03045#S6.F5)shows convergence curves for the remaining 11 problems\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x4.png)Figure 5:Regret convergence for the 11 benchmark problems not shown in Figure[4](https://arxiv.org/html/2608.03045#S6.F4)\. Bi\-BO \(SLSQP\) \(red\) converges to values that are 1–8 orders of magnitude below black\-box BO \(blue\) on every problem\. Each curve usesninit=50n\_\{\\text\{init\}\}=50initial Latin hypercube samples and the EI exploration parameterξ\\xithat yielded the lowest final regret for that problem \(see Appendix[J](https://arxiv.org/html/2608.03045#A10)for per\-problem values\)\. Solid lines show the mean over 10 independent repetitions; shaded bands show±1\\pm 1standard deviation\.Table 3:Mean final regret and mean wall time \(seconds\) after 200 iterations \(ninit=50n\_\{\\text\{init\}\}=50, bestξ\\xiper method, 10 repetitions; see Appendix[J](https://arxiv.org/html/2608.03045#A10)for per\-problemξ\\xivalues\)\. Bi\-BO \(SLSQP\) achieves lower regret on all 13 problems; Bi\-BO \(BH\) on 12 of 13 \(†BH inner solver is worse than BB\-BO on SFR\-1\)\.Both bilevel variants achieve lower regret than black\-box BO\. Bi\-BO \(SLSQP\) wins on all 13 problems \(11×\\times–108×10^\{8\}\\timeslower; geometric mean 3,192×\\times\) at wall time comparable to BB\-BO\. Bi\-BO \(BH\) wins on 12 of 13 \(geometric mean 4,840×\\times\), but Basin\-Hopping’s random perturbations overshoot the narrow feasible region of SFR\-1 and leave it worse than BB\-BO there\. Gains scale with constraint tightness: Heat\-Exchanger \(108×10^\{8\}\\times\) and Distillation \(106×10^\{6\}\\times\) expose the penalty method’s failure in five dimensions, while the smallest SLSQP gains \(SFR\-1: 11×\\times, Evaporator: 26×\\times\) occur where the white\-box subproblem itself has a narrow feasible region or multiple local optima\.

Wall\-time differences concentrate on problems with large infeasible regions inxBBx^\{\\mathrm\{BB\}\}\. When a sampledxBBx^\{\\mathrm\{BB\}\}admits no feasible process variables, Basin\-Hopping keeps perturbing in search of a feasible point rather than returning quickly, inflating runtime to 826 s on SFR\-1 and 3,130 s on Williams\-Otto \(4–10×\\timesthe SLSQP inner solver\)\. SLSQP exits on infeasibility and lets the outer BO loop advance\.

The full\-space BH baseline fails on Rastrigin \(multimodal joint space\) and performs poorly on Distillation; Bi\-BO matches or approaches its quality on the remaining problems using∼\\sim250 black\-box evaluations versus the substantially higher counts BH requires\. All 13 SLSQP pairwise comparisons are statistically significant \(Wilcoxon signed\-rank, Holm–Bonferroni adjustedp<0\.05p<0\.05; sign test 13/13,p=0\.000244p=0\.000244\); 12 of 13 hold for BH \(sign test,p=0\.0034p=0\.0034\)\. A Friedman test across the four methods confirms differences \(χ2=17\.03\\chi^\{2\}=17\.03,p=0\.0007p=0\.0007\), with Nemenyi post\-hoc placing BB\-BO below both Bi\-BO variants and NLP \(critical difference=1\.30=1\.30\)\.

When using BH as the inner solver for Bi\-BO, the performance on the small feasible region problem \(SFR\-1\) is worse than the black\-box BO baseline\. Similarly, the solution times for the small feasible region problem and the Williams\-Otto problem are significantly higher for the Bi\-BO \(BH\) variant compared to both the Bi\-BO \(SLSQP\) variant and the black\-box BO baseline\. Both of these problems are characterized by narrow feasible regions, and regions which are completely infeasible for certainxBBx^\{\\mathrm\{BB\}\}values\. The random perturbations in the BH algorithm cause it to overshoot these narrow feasible regions, or to iterate until satisfied that no feasible solution exists\. This highlights a key tradeoff between the two inner solvers: BH provides stronger global optimization guarantees but can be inefficient in problems with narrow or disconnected feasible regions, while SLSQP is more efficient but may struggle with multimodality\. The choice of inner solver should therefore be informed by problem characteristics\.

Figure[6](https://arxiv.org/html/2608.03045#S6.F6)visualizes the relationship between final regret for black\-box BO and Bi\-BO \(SLSQP\) across all 273 \(problem,ninitn\_\{\\text\{init\}\},ξ\\xi\) configurations\. Every point lies below the diagonal, confirming that Bi\-BO outperforms BB\-BO regardless of hyperparameter choice\. The wide range of regret ratios \(11×\\times–108×10^\{8\}\\times\) is visible in the vertical spread of points, with the largest gains occurring on problems with tight constraints and high dimensionality\. The greater vertical spread in Bi\-BO final regret \(compared to the relatively narrow horizontal spread in BB\-BO\) reflects the sensitivity of the inner solver to problem landscape: on problems with multimodal or non\-convex inner subproblems \(e\.g\., Evaporator, Rosen\-Suzuki\), the multi\-start SLSQP inner solver occasionally converges to local optima, producing variable final regret across hyperparameter configurations\. On problems with well\-behaved inner landscapes \(e\.g\., Heat\-Exchanger, Membrane\), Bi\-BO regret is consistently near zero\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x5.png)Figure 6:Scatter plot of final regret: black\-box BO vs\. bilevel BO \(SLSQP variant\)\. Each point is one \(problem,ninitn\_\{\\text\{init\}\},ξ\\xi\) configuration across 273 total configurations\. All points lie below the diagonal, confirming that Bi\-BO \(SLSQP\) outperforms BB\-BO regardless of hyperparameter choice\.

## 7Discussion

Bi\-BO achieves 11×\\times–108×10^\{8\}\\timeslower regret than black\-box BO on every problem tested \(SLSQP variant; geometric mean 3,192×\\times\), with equal or faster wall time on 12 of 13 problems\. The source of this advantage is structural, not algorithmic: we use the same GP kernel, the same acquisition function, and the same optimizer as the black\-box baseline\. The gains come entirely from exploiting the separable problem structure—reducing the surrogate dimension and satisfying the white\-box constraints exactly\. Three additional analyses—feasibility rates, the relationship between regret ratio and problem characteristics, and sample efficiency crossover—provide mechanistic insight into when and why the bilevel decomposition helps\.

Constraint satisfaction and wasted evaluations\.A key mechanism underlying the regret gap is that black\-box BO wastes evaluations on infeasible queries\. Figure[7](https://arxiv.org/html/2608.03045#S7.F7)reports the fraction of all 250 evaluations that land in feasible regions, reconstructed from stored search histories\. Black\-box BO wastes 6–72% of its budget on infeasible points depending on the problem, with Small\-Feasible\-Region \(72%\), Batch\-Reactor \(71%\), and Evaporator \(58%\) being the worst cases\. Bilevel BO achieves 100% feasibility on 8 of 12 constrained problems because the inner BH solver handles constraint satisfaction exactly\. On the remaining four problems \(Batch\-Reactor, Evaporator, Small\-Feasible\-Region, PSA\), bilevel BO feasibility ranges from 59–91%\. This occurs for two reasons: \(1\) the inner BH may return solutions where the total constraint violationmaxi⁡gi​\(xWB,y\)\\max\_\{i\}g\_\{i\}\(x^\{\\mathrm\{WB\}\},y\)is positive but small \(on the order of10−410^\{\-4\}–10−210^\{\-2\}\), placing the solution near but outside the feasible boundary; or \(2\) certainy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)values yield an inner problem with no feasible solution, because noxWB∈𝒳WBx^\{\\mathrm\{WB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\}satisfiesg​\(xWB,y\)≤0g\(x^\{\\mathrm\{WB\}\},y\)\\leq 0\. This is not a limitation of the bilevel architecture itself but of the inner solver; deterministic global solvers \(e\.g\., BARON\) would restore the 100% guarantee\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x6.png)Figure 7:Fraction of evaluations landing in feasible regions across 12 constrained problems \(ninit=50n\_\{\\text\{init\}\}=50, bestξ\\xi, mean over 10 repetitions\)\. Black\-box BO \(blue\) wastes 6–72% of its budget on infeasible queries; bilevel BO \(red and orange\) achieves≥\\geq59% feasibility on all problems and 100% on 8 of 12\.Factors influencing the regret ratio\.Factors beyond raw problem dimensionality such as constraint geometry, inner\-problem multimodality, and black\-box landscape complexity determine the magnitude of improvement from BB\-BO to Bi\-BO\. For example, Rastrigin \(unconstrained, 3D total\) achieves 698×\\timesimprovement with the SLSQP inner solver \(Table[3](https://arxiv.org/html/2608.03045#S6.T3)\) purely from reducing the GP surrogate from 3D to 1D\. Meanwhile, Evaporator and Rosen\-Suzuki have the smallest SLSQP ratios \(26×\\timesand 201×\\times, respectively\), because their non\-convex inner landscapes cause the multi\-start SLSQP to find local optima\. Heat\-Exchanger has a loose constraint set \(93% feasible\) but the largest gain \(108×10^\{8\}\\times\), driven by dimensionality reduction from 5D to 2D combined with a well\-behaved inner landscape\. The analysis confirms that when the inner optimizer converges reliably, dimensionality reduction alone produces large gains; exact constraint satisfaction provides additional benefit but is neither necessary nor sufficient\.

Sample efficiency\.Beyond final regret, practitioners are interested in how quickly a method reaches useful solutions\. Figure[8](https://arxiv.org/html/2608.03045#S7.F8)reports the number of evaluations at which bilevel BO first matches the final regret that black\-box BO achieves after all 250 evaluations\. On 5 of 13 problems \(Heat\-Exchanger, Membrane, Williams\-Otto, Distillation, CSTR\), bilevel BO matches or surpasses black\-box BO’s final performance within the initial sampling phase alone—effectively from the first BO iteration\. Even in the hardest case \(Evaporator\), crossover occurs at 47 evaluations \(19% of the total budget\)\. This has direct practical implications: when each black\-box evaluation costs hours of compute time, bilevel BO can deliver black\-box BO\-quality solutions using 80–100% fewer expensive evaluations\.

![Refer to caption](https://arxiv.org/html/2608.03045v1/x7.png)Figure 8:Sample efficiency crossover: number of evaluations for bilevel BO to match black\-box BO’s final regret \(at 250 evaluations\)\. Percentages show the fraction of the total budget required\. On 10 of 13 problems, bilevel BO matches black\-box BO’s final performance within 10% of the evaluation budget\.Connection to decomposition in global optimization\.As discussed in Section[3](https://arxiv.org/html/2608.03045#S3), the bilevel reformulation belongs to a broader family of decomposition strategies—block coordinate descent\[[40](https://arxiv.org/html/2608.03045#bib.bib40)\], trust\-region filter methods\[[4](https://arxiv.org/html/2608.03045#bib.bib4),[13](https://arxiv.org/html/2608.03045#bib.bib13)\], and data\-driven bilevel solvers like DOMINO\[[8](https://arxiv.org/html/2608.03045#bib.bib8)\]\. Our results provide empirical support for a principle from bilevel optimization theory\[[41](https://arxiv.org/html/2608.03045#bib.bib41)\]: solving the inner problem exactly improves outer\-level convergence\. In the surrogate\-based setting, the inner global optimizer eliminates noise from outer\-level observations, giving the GP a cleaner signal\. This suggests that bilevel decomposition may be beneficial whenever part of the optimization landscape admits an efficient exact solver, even beyond the grey\-box problems studied here\.

Inner\-solver choice: SLSQP vs\. Basin\-Hopping\.The choice of inner\-loop solver affects both solution quality and computational cost\. Across all 13 problems, multi\-start SLSQP and Basin\-Hopping achieve comparable final regret \(family\-level sign test: SLSQP wins 6, BH wins 7,p=1\.0p=1\.0; geometric mean SLSQP/BH ratio: 1\.52×\\times\)\. However, the two solvers exhibit complementary strengths\. BH excels on multimodal inner landscapes: on Rastrigin, Bi\-BO \(BH\) achieves7\.2×107×7\.2\\times 10^\{7\}\\timesimprovement over BB\-BO, while Bi\-BO \(SLSQP\) achieves 698×\\times—a gap of five orders of magnitude attributable entirely to SLSQP’s inability to escape local optima in the multimodal Rastrigin inner problem\. Conversely, SLSQP excels on problems with narrow feasible regions: on SFR\-1, Bi\-BO \(SLSQP\) achieves 11×\\timesimprovement while Bi\-BO \(BH\) is 0\.1×\\times\(*worse*than BB\-BO\), because Basin\-Hopping’s random perturbations overshoot the small feasible set\. SLSQP is also faster on all 13 problems \(geometric mean 1\.74×\\times\), making it the recommended default\. BH is preferred when the inner landscape is known or suspected to be multimodal\.

Complementarity with existing grey\-box BO\.As detailed in Section[3](https://arxiv.org/html/2608.03045#S3), COBALT\[[6](https://arxiv.org/html/2608.03045#bib.bib6)\]and BOCF\[[5](https://arxiv.org/html/2608.03045#bib.bib5)\]address settings our method does not: non\-separable composite functions and risk\-averse decisions via uncertainty propagation\. Our method trades their generality for simplicity and exactness when Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)holds\. A direct comparison on the same benchmark suite is planned for future work\.

Baseline constraint handling\.The black\-box BO baseline uses a fixed penalty \(10610^\{6\}\) for constraint violations\. More sophisticated constrained BO methods—feasibility\-weighted expected improvement\[[26](https://arxiv.org/html/2608.03045#bib.bib26)\], augmented Lagrangian approaches\[[42](https://arxiv.org/html/2608.03045#bib.bib42)\]—would likely narrow the gap on tightly constrained problems such as Heat\-Exchanger and Distillation, where the penalty method is least effective\. However, the dimensionality reduction advantage of bilevel BO would persist regardless of constraint handling, as the surrogate dimension reduction fromnWB\+nBBn\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}tonBBn\_\{\\mathrm\{BB\}\}is independent of the acquisition function\. We intentionally use vanilla methods for both BO and NLP to isolate the effect of the bilevel structure rather than specific algorithmic choices\.

Limitations and open problems\.We identify six limitations, each pointing to an open problem\. \(1\) The inner optimizer must converge to the global optimum \(Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\)\)\. For non\-convex inner problems, the inner solver may converge to local optima—as observed in Section[6](https://arxiv.org/html/2608.03045#S6), this explains the modest 11×\\timesimprovement on SFR\-1 \(SLSQP variant\) and the sub\-100% feasibility on Batch\-Reactor and Evaporator \(Figure[7](https://arxiv.org/html/2608.03045#S7.F7)\)\. Replacing the stochastic inner solver with a deterministic global solver \(e\.g\., BARON\) would provide optimality certificates but at higher computational cost\. \(2\) The inner solver can perform*worse*than the black\-box baseline on problems with very narrow feasible regions: on SFR\-1, the BH inner solver’s random perturbations generate predominantly infeasible starting points, yielding a 0\.1×\\timesratio\. Multi\-start SLSQP avoids this failure mode on SFR\-1 but lacks BH’s global exploration for multimodal problems\. \(3\) Expensive inner optimizations can dominate wall time\. With the BH inner solver, Williams\-Otto requires 3,130s vs\. 273s for BB\-BO \(11\.5×\\timesoverhead\); even the faster SLSQP variant incurs a 1\.19×\\timesoverhead on this problem\. Warm\-starting the inner solver from previous solutions or using reduced\-order inner models could mitigate this\. \(4\) Pure black\-box constraints—those depending directly onfBBf^\{\\mathrm\{BB\}\}outputs without passing through the white\-box model—are not handled; extending the outer loop with constrained BO methods\[[26](https://arxiv.org/html/2608.03045#bib.bib26),[27](https://arxiv.org/html/2608.03045#bib.bib27)\]would address this\. \(5\) All black\-box functions in our suite are closed\-form surrogates of what would be expensive computations; validation on a real DFT or molecular dynamics black box would strengthen the practical case\. \(6\) Formal convergence guarantees for the outer BO loop—which follow from standard results since it is BO over a compact domain with a GP surrogate—are not analyzed here; characterizing how the bilevel structure affects dimension\-dependent regret bounds is future work\.

Broader implications\.Expensive black\-box optimization problems frequently contain exploitable structure—known physics, analytical submodels, or differentiable components—that monolithic surrogates ignore\. The bilevel decomposition is one instance of a general principle:*solve what you can solve exactly, and approximate only what you must*\. As surrogate\-based methods scale to higher\-dimensional search spaces, structural decomposition becomes not merely helpful but necessary to overcome the curse of dimensionality\. The 13\-problem benchmark suite provides a testbed for developing and comparing methods that exploit this structure\.

## 8Conclusion

When grey\-box optimization problems exhibit variable separability—the expensive black\-box depends only on one variable group while the remaining variables enter known, differentiable equations—solving the known subproblem exactly and approximating only the black\-box component yields 11×\\times–108×10^\{8\}\\timeslower regret than monolithic surrogate\-based optimization\. The advantage is robust to hyperparameters and inner\-solver choice, and requires no algorithmic novelty beyond the decomposition itself\. A suite of 13 benchmark problems with verified global optima provides a standardized testbed for future methods that exploit this structure\.

The opening observation of this paper—that monolithic surrogates waste samples learning what is already known—is nearly self\-evident once the separable structure is recognized\. Yet no prior method systematically exploits this structure for dimensionality reduction with exact constraint satisfaction \(in principle; subject to inner\-solver global optimality in practice\)\. The empirical evidence across 8,450 independent optimization runs confirms that the payoff is large and consistent\. As expensive black\-box optimization scales to higher dimensions, identifying and exploiting known substructure will be essential; the bilevel decomposition studied here is one principled way to do so\.

## Acknowledgements

Computational resources were provided by the Texas Advanced Computing Center \(TACC\) at The University of Texas at Austin\.

## Declarations

- •Funding\.This work was supported by the ExxonMobil, whose financial support is acknowledged with gratitude\.
- •Conflict of interest\.The authors declare that they have no conflict of interest\.
- •Data availability\.All data generated during this study are reproducible from the code described below\. No external datasets were used\.
- •
- •Author contributions\.Joshua E\. Hammond: conceptualization, methodology, software, investigation, writing—original draft\. Tyler A\. Soderstrom: methodology, investigation, writing—review & editing\. Brian A\. Korgel: supervision, writing—review & editing\. Michael Baldea: conceptualization, supervision, funding acquisition, writing—review & editing\.
- •Use of generative AI\.Generative AI tools \(Claude Opus 4\.6, Anthropic\) were used to assist with code development, data analysis, and manuscript drafting\. All AI\-generated content was reviewed, verified, and edited by the authors, who take full responsibility for the accuracy and integrity of the work\. No AI tool was used to generate scientific claims, interpret results, or design experiments without human oversight\.

## Appendix ATest Problem Definitions

Full mathematical formulations for all 13 benchmark problems\. Each problem specifies the white\-box and black\-box models, constraints, variable bounds, and verified global optimum\.

### A\.1Small Feasible Region Problem 1

Gardner et al\.\[[26](https://arxiv.org/html/2608.03045#bib.bib26)\]and Ariafar et al\.\[[39](https://arxiv.org/html/2608.03045#bib.bib39)\]propose a constrained optimization test problem with a small feasible region that can be formulated as follows:

minx\\displaystyle\\min\\limits\_\{x\}sin⁡x1\+y1,\\displaystyle\\sin\{x\_\{1\}\}\+y\_\{1\},\(8\)s\.t\.g1​\(x\):=sin⁡\(x1\)​sin⁡\(y1\)\+0\.95≤0,\\displaystyle~~g\_\{1\}\(x\):=\\sin\(x\_\{1\}\)\\sin\(y\_\{1\}\)\+0\.95\\leq 0,y1=fBB​\(x\):=x2,\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\(x\):=x\_\{2\},0≤xi≤6,∀i∈\{1,2\}\\displaystyle~~0\\leq x\_\{i\}\\leq 6,~~~~\\forall i\\in\\\{1,2\\\}
The global minimum isJ⋆=0\.2532J^\{\\star\}=0\.2532and is located atx1WB⁣⋆=3​π2x^\{\\mathrm\{WB\}\\star\}\_\{1\}=\\frac\{3\\pi\}\{2\},xBB⁣⋆=x2=1\.2532x^\{\\mathrm\{BB\}\\star\}=x\_\{2\}=1\.2532\.

### A\.2Modified Small Feasible Region Problem 2

Gardner et al\.\[[26](https://arxiv.org/html/2608.03045#bib.bib26)\]and Ariafar et al\.\[[39](https://arxiv.org/html/2608.03045#bib.bib39)\]propose a constrained optimization test problem with a small feasible region\. We add an additional black\-box functiony1=fBB​\(x2\)y\_\{1\}=f^\{\\mathrm\{BB\}\}\(x\_\{2\}\)that forms a grey\-box problem as follows:

minx\\displaystyle\\min\\limits\_\{x\}cos⁡\(2​x1\)​cos⁡\(y1\)\+sin⁡\(x1\),\\displaystyle\\cos\(2x\_\{1\}\)\\cos\(y\_\{1\}\)\+\\sin\(x\_\{1\}\),\(9\)s\.t\.g1​\(x\):=cos⁡\(x1\)​cos⁡\(y1\)−sin⁡\(x1\)​sin⁡\(y1\)−0\.5≤0,\\displaystyle~~g\_\{1\}\(x\):=\\cos\(x\_\{1\}\)\\cos\(y\_\{1\}\)\-\\sin\(x\_\{1\}\)\\sin\(y\_\{1\}\)\-0\.5\\leq 0,y1=fBB​\(x\):=\(x2−5\.5\)2​exp⁡\(x22\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\(x\):=\\frac\{\(x\_\{2\}\-5\.5\)\}\{2\}\\exp\\left\(\\frac\{x\_\{2\}\}\{2\}\\right\),0≤xi≤6,∀i∈\{1,2\}\\displaystyle~~0\\leq x\_\{i\}\\leq 6,~~~~\\forall i\\in\\\{1,2\\\}
The global minimumJ⋆=−2\.0J^\{\\star\}=\-2\.0is located atxBB⁣⋆=5\.5,xWB⁣⋆=3​π/2x^\{\\mathrm\{BB\}\\star\}=5\.5,x^\{\\mathrm\{WB\}\\star\}=3\\pi/2\.

### A\.3Rastrigin

We use a modified formulation of the Rastrigin function\[[43](https://arxiv.org/html/2608.03045#bib.bib43)\]that can be formulated as a multi\-scale grey\-box problem as follows:

minx,y\\displaystyle\\min\\limits\_\{x,y\}30\+x12−10​cos⁡\(2​π​x1\)\+x22−10​cos⁡\(2​π​x2\)\+y1,\\displaystyle~~30\+x\_\{1\}^\{2\}\-10\\cos\(2\\pi x\_\{1\}\)\+x\_\{2\}^\{2\}\-10\\cos\(2\\pi x\_\{2\}\)\+y\_\{1\},\(10\)s\.t\.y1=fBB​\(x\):=x32−10​cos⁡\(2​π​x3\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\(x\):=x\_\{3\}^\{2\}\-10\\cos\(2\\pi x\_\{3\}\),−5\.12≤xi≤5\.12,∀i∈\{1,2,3\},\\displaystyle~~\-5\.12\\leq x\_\{i\}\\leq 5\.12,~~~~~\\forall i\\in\\\{1,2,3\\\},x=\[xWB,xBB\],xWB=\[x1,x2\],xBB=\[x3\]\.\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],~~~~~x^\{\\mathrm\{WB\}\}=\[x\_\{1\},x\_\{2\}\],~~~~~x^\{\\mathrm\{BB\}\}=\[x\_\{3\}\]\.The global minimum is equal to0withxWB⁣⋆=\[0,0\],xBB⁣⋆=\[0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[0,0\],x^\{\\mathrm\{BB\}\\star\}=\[0\]^\{\\top\}\.

### A\.4Toy\-Hydrology

We use a modified formulation of the Toy Hydrology problem\[[42](https://arxiv.org/html/2608.03045#bib.bib42)\]that can be formulated as a multi\-scale grey\-box problem as follows:

minx,y\\displaystyle\\min\\limits\_\{x,y\}x1\+x2,\\displaystyle~~x\_\{1\}\+x\_\{2\},\(11\)s\.t\.g1​\(x,y\):=1\.5−x1−2​x2−0\.5​sin⁡\(−4​π​x2\+y1\)≤0,\\displaystyle~~g\_\{1\}\(x,y\):=1\.5\-x\_\{1\}\-2x\_\{2\}\-0\.5\\sin\(\-4\\pi x\_\{2\}\+y\_\{1\}\)\\leq 0,g2​\(x,y\):=x12\+x22−1\.5≤0,\\displaystyle~~g\_\{2\}\(x,y\):=x\_\{1\}^\{2\}\+x\_\{2\}^\{2\}\-1\.5\\leq 0,y1=fBB​\(x1\):=2​π​x12,\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\(x\_\{1\}\):=2\\pi x\_\{1\}^\{2\},0≤xi≤1,∀i∈\{1,2\},\\displaystyle~~0\\leq x\_\{i\}\\leq 1,~~~~~\\forall i\\in\\\{1,2\\\},x=\[xWB,xBB\],xBB=\[x1\],xWB=\[x2\]\.\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],~~~~~x^\{\\mathrm\{BB\}\}=\[x\_\{1\}\],~~~~~x^\{\\mathrm\{WB\}\}=\[x\_\{2\}\]\.The global minimum is0\.59980\.5998withxWB⁣⋆=0\.4047,xBB⁣⋆=0\.1951x^\{\\mathrm\{WB\}\\star\}=0\.4047,x^\{\\mathrm\{BB\}\\star\}=0\.1951\.

### A\.5Rosen\-Suzuki

We use a modified formulation of the Rosen\-Suzuki function\[[44](https://arxiv.org/html/2608.03045#bib.bib44)\]that can be formulated as a multi\-scale grey\-box problem as follows:

minx,y\\displaystyle\\min\\limits\_\{x,y\}x12\+x22\+x42−5​x1−5​x2\+y1,\\displaystyle~~x\_\{1\}^\{2\}\+x\_\{2\}^\{2\}\+x\_\{4\}^\{2\}\-5x\_\{1\}\-5x\_\{2\}\+y\_\{1\},\(12\)s\.t\.g1​\(x,y\):=−\(8−x12−x22−x32−x42−x1\+x2−x3\+x4\)≤0,\\displaystyle~~g\_\{1\}\(x,y\):=\-\(8\-x\_\{1\}^\{2\}\-x\_\{2\}^\{2\}\-x\_\{3\}^\{2\}\-x\_\{4\}^\{2\}\-x\_\{1\}\+x\_\{2\}\-x\_\{3\}\+x\_\{4\}\)\\leq 0,g2​\(x,y\):=−\(10−x12−2​x22−y2\+x1\+x4\)≤0,\\displaystyle~~g\_\{2\}\(x,y\):=\-\(10\-x\_\{1\}^\{2\}\-2x\_\{2\}^\{2\}\-y\_\{2\}\+x\_\{1\}\+x\_\{4\}\)\\leq 0,g3​\(x,y\):=−\(5−2​x12−x22−x32−2​x1\+x2\+x4\)≤0,\\displaystyle~~g\_\{3\}\(x,y\):=\-\(5\-2x\_\{1\}^\{2\}\-x\_\{2\}^\{2\}\-x\_\{3\}^\{2\}\-2x\_\{1\}\+x\_\{2\}\+x\_\{4\}\)\\leq 0,y1=f1BB​\(z\):=2​x32−21​x3\+7​x4,\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(z\):=2x\_\{3\}^\{2\}\-21x\_\{3\}\+7x\_\{4\},y2=f2BB​\(z\):=x32\+2​x42,\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(z\):=x\_\{3\}^\{2\}\+2x\_\{4\}^\{2\},−2≤xi≤2,∀i∈\{1,…,4\}\\displaystyle~~\-2\\leq x\_\{i\}\\leq 2,~~~~~\\forall i\\in\\\{1,\\ldots,4\\\}x=\[xWB,xBB\],xWB=\[x1,x2\],xBB=\[x3,x4\]\.\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],~~~~~x^\{\\mathrm\{WB\}\}=\[x\_\{1\},x\_\{2\}\],~~~~~x^\{\\mathrm\{BB\}\}=\[x\_\{3\},x\_\{4\}\]\.The global minimum is−44\-44withxWB⁣⋆=\[0,1\]⊤,xBB⁣⋆=\[2,−1\]⊤x^\{\\mathrm\{WB\}\\star\}=\[0,1\]^\{\\top\},x^\{\\mathrm\{BB\}\\star\}=\[2,\-1\]^\{\\top\}\.

### A\.6Catalytic Reactor Design

We use a catalytic reactor design problem that combines a white\-box CSTR process model with black\-box catalyst properties akin to those derived from DFT calculations\. The objective is to*maximize*the yield of intermediate productBBin consecutive reactionsA→B→CA\\to B\\to C\. The problem can be formulated as a multi\-scale grey\-box problem as follows:

maxT,P,τ,Eb,Δ​Ea\\displaystyle\\max\\limits\_\{T,P,\\tau,E\_\{b\},\\Delta E\_\{a\}\}YB\\displaystyle~~Y\_\{B\}\(13\)s\.t\.𝐟WB=\[YB−k1​τ\(1\+k1​τ\)​\(1\+k2​τ\)k1​\(T,P,Eb,Δ​Ea\)−A1​v​\(Eb\)​exp⁡\(−Ea,1R​T\)​\(PPref\)β1​1\(1\+Kads​P\)2k2​\(T,P,Eb,Δ​Ea\)−A2​exp⁡\(−Ea,2R​T\)​\(PPref\)β2​1\(1\+Kads​P\)2Kads​\(Eb\)−K0​exp⁡\(κ​\(Eb−Eb∗\)\)\],\\displaystyle~~\{\\bf f\}^\{\\mathrm\{WB\}\}=\\left\[\\begin\{array\}\[\]\{l\}Y\_\{B\}\-\\frac\{k\_\{1\}\\tau\}\{\\left\(1\+k\_\{1\}\\tau\\right\)\\left\(1\+k\_\{2\}\\tau\\right\)\}\\\\ k\_\{1\}\\left\(T,P,E\_\{b\},\\Delta E\_\{a\}\\right\)\-A\_\{1\}v\\left\(E\_\{b\}\\right\)\\exp\\left\(\-\\frac\{E\_\{a,1\}\}\{RT\}\\right\)\\left\(\\frac\{P\}\{P\_\{\\mathrm\{ref\}\}\}\\right\)^\{\\beta\_\{1\}\}\\frac\{1\}\{\\left\(1\+K\_\{\\mathrm\{ads\}\}P\\right\)^\{2\}\}\\\\ k\_\{2\}\\left\(T,P,E\_\{b\},\\Delta E\_\{a\}\\right\)\-A\_\{2\}\\exp\\left\(\-\\frac\{E\_\{a,2\}\}\{RT\}\\right\)\\left\(\\frac\{P\}\{P\_\{\\mathrm\{ref\}\}\}\\right\)^\{\\beta\_\{2\}\}\\frac\{1\}\{\\left\(1\+K\_\{\\mathrm\{ads\}\}P\\right\)^\{2\}\}\\\\ K\_\{\\mathrm\{ads\}\}\\left\(E\_\{b\}\\right\)\-K\_\{0\}\\exp\\left\(\\kappa\\left\(E\_\{b\}\-E\_\{b\}^\{\*\}\\right\)\\right\)\\end\{array\}\\right\],\(18\)𝐟BB=\[v​\(Eb\)−\[exp⁡\(−12​\(Eb−Eb∗σE\)2\)\]Ea,1−\[Ea,1\(0\)−α1​Δ​Ea\]Ea,2−\[Ea,2\(0\)\+α2​Δ​Ea\]\],\\displaystyle~~\{\\bf f\}^\{\\mathrm\{BB\}\}=\\left\[\\begin\{array\}\[\]\{l\}v\\left\(E\_\{b\}\\right\)\-\\left\[\\exp\\left\(\-\\frac\{1\}\{2\}\\left\(\\frac\{E\_\{b\}\-E\_\{b\}^\{\*\}\}\{\\sigma\_\{E\}\}\\right\)^\{2\}\\right\)\\right\]\\\\ E\_\{a,1\}\-\\left\[E\_\{a,1\}^\{\(0\)\}\-\\alpha\_\{1\}\\Delta E\_\{a\}\\right\]\\\\ E\_\{a,2\}\-\\left\[E\_\{a,2\}^\{\(0\)\}\+\\alpha\_\{2\}\\Delta E\_\{a\}\\right\]\\end\{array\}\\right\],\(22\)g1​\(x,y\):=YC−0\.25≤0,\\displaystyle~~g\_\{1\}\(x,y\):=Y\_\{C\}\-0\.25\\leq 0,g2​\(x,y\):=PPmax−T−TminTmax−Tmin−0\.3≤0,\\displaystyle~~g\_\{2\}\(x,y\):=\\frac\{P\}\{P\_\{\\max\}\}\-\\frac\{T\-T\_\{\\min\}\}\{T\_\{\\max\}\-T\_\{\\min\}\}\-0\.3\\leq 0,
and the conversion and byproduct yield are:

XA\\displaystyle X\_\{A\}=k1​τ1\+k1​τ,\\displaystyle=\\frac\{k\_\{1\}\\tau\}\{1\+k\_\{1\}\\tau\},YC\\displaystyle Y\_\{C\}=XA−YB\.\\displaystyle=X\_\{A\}\-Y\_\{B\}\.
The black\-box functionfBBf^\{\\mathrm\{BB\}\}maps catalyst propertiesΔ​Ea\\Delta E\_\{a\},EbE\_\{b\}to kinetic parameters via a volcano\-type activity model:

Table 4:Parameters for the catalytic reactor design problem\.The global maximum yield given the parameters in Table[4](https://arxiv.org/html/2608.03045#A1.T4)is92\.32%92\.32\\%withxWB⁣⋆=\[647\.3,20\.0,5\.0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[647\.3,20\.0,5\.0\]^\{\\top\},xBB⁣⋆=\[−1\.047,0\.200\]⊤x^\{\\mathrm\{BB\}\\star\}=\[\-1\.047,0\.200\]^\{\\top\}\.

### A\.7Heat Exchanger with Novel Working Fluid

We consider the design of a heat exchanger system where the working fluid’s thermodynamic properties \(heat capacitycpc\_\{p\}, viscosityμ\\mu, and thermal conductivitykk\) depend on two molecular descriptors \(m1m\_\{1\},m2m\_\{2\}\) in a black\-box function\. White\-box equations govern the heat exchanger design and operation \(heat transfer areaAA, mass flow ratem˙\\dot\{m\}, and log\-mean temperature differenceΔ​Tl​m\\Delta T\_\{lm\}\) with the goal of minimizing the total cost \(capital \+ operating\) while satisfying minimum heat duty\[[45](https://arxiv.org/html/2608.03045#bib.bib45)\]\. Capital cost scales with area according to the six\-tenths rule\[[46](https://arxiv.org/html/2608.03045#bib.bib46)\]while operating cost is a function of pumping power which depends on fluid viscosity and mass flow rate according to the relation provided by\[[47](https://arxiv.org/html/2608.03045#bib.bib47)\]\. We use a simplified form for the overall heat transfer coefficientUUbased the resistance model from Ravagnani and Caballero\[[48](https://arxiv.org/html/2608.03045#bib.bib48)\]\.

The problem can be formulated as:

minA,m˙,Δ​Tl​m,m1,m2\\displaystyle\\min\\limits\_\{A,\\dot\{m\},\\Delta T\_\{lm\},m\_\{1\},m\_\{2\}\}Ccapital\+Coperating\\displaystyle~~C\_\{\\text\{capital\}\}\+C\_\{\\text\{operating\}\}\(23\)s\.t\.Ccapital=1000⋅x1WB0\.6,\\displaystyle~~C\_\{\\text\{capital\}\}=1000\\cdot\{x^\{\\mathrm\{WB\}\}\_\{1\}\}^\{0\.6\},Coperating=500⋅x2WB3⋅y2x1WB,\\displaystyle~~C\_\{\\text\{operating\}\}=500\\cdot\\frac\{\{x^\{\\mathrm\{WB\}\}\_\{2\}\}^\{3\}\\cdot y\_\{2\}\}\{x^\{\\mathrm\{WB\}\}\_\{1\}\},U=y30\.01\+0\.001/y1,\\displaystyle~~U=\\frac\{y\_\{3\}\}\{0\.01\+0\.001/y\_\{1\}\},g1​\(x,y\):=100−U⋅x1WB⋅x3WB≤0,\\displaystyle~~g\_\{1\}\(x,y\):=100\-U\\cdot x^\{\\mathrm\{WB\}\}\_\{1\}\\cdot x^\{\\mathrm\{WB\}\}\_\{3\}\\leq 0,g2​\(x,y\):=4000−4​x2WBπ⋅0\.05⋅y2≤0,\\displaystyle~~g\_\{2\}\(x,y\):=4000\-\\frac\{4x^\{\\mathrm\{WB\}\}\_\{2\}\}\{\\pi\\cdot 0\.05\\cdot y\_\{2\}\}\\leq 0,y1=f1BB​\(x1BB,x2BB\):=2\.0\+0\.5​sin⁡\(π​x1BB50\)\+0\.3​x2BB,\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(x^\{\\mathrm\{BB\}\}\_\{1\},x^\{\\mathrm\{BB\}\}\_\{2\}\):=2\.0\+0\.5\\sin\\left\(\\frac\{\\pi x^\{\\mathrm\{BB\}\}\_\{1\}\}\{50\}\\right\)\+0\.3x^\{\\mathrm\{BB\}\}\_\{2\},y2=f2BB​\(x1BB,x2BB\):=0\.001⋅exp⁡\(x1BB−10050\)⋅\(1\+0\.5​x2BB\),\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(x^\{\\mathrm\{BB\}\}\_\{1\},x^\{\\mathrm\{BB\}\}\_\{2\}\):=0\.001\\cdot\\exp\\left\(\\frac\{x^\{\\mathrm\{BB\}\}\_\{1\}\-100\}\{50\}\\right\)\\cdot\(1\+0\.5x^\{\\mathrm\{BB\}\}\_\{2\}\),y3=f3BB​\(x1BB,x2BB\):=0\.1\+0\.05​\(1−\(x1BB−12575\)2\)\+0\.02​x2BB,\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(x^\{\\mathrm\{BB\}\}\_\{1\},x^\{\\mathrm\{BB\}\}\_\{2\}\):=0\.1\+0\.05\\left\(1\-\\left\(\\frac\{x^\{\\mathrm\{BB\}\}\_\{1\}\-125\}\{75\}\\right\)^\{2\}\\right\)\+0\.02x^\{\\mathrm\{BB\}\}\_\{2\},A∈\[1,50\],m˙∈\[0\.1,5\],Δ​Tl​m∈\[5,50\],\\displaystyle~~A\\in\[1,50\],\\quad\\dot\{m\}\\in\[0\.1,5\],\\quad\\Delta T\_\{lm\}\\in\[5,50\],m1∈\[50,200\],m2∈\[0,1\],\\displaystyle~~m\_\{1\}\\in\[50,200\],\\quad m\_\{2\}\\in\[0,1\],x=\[xWB,xBB\],xWB=\[A,m˙,Δ​Tl​m\],xBB=\[m1,m2\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[A,\\dot\{m\},\\Delta T\_\{lm\}\],\\quad x^\{\\mathrm\{BB\}\}=\[m\_\{1\},m\_\{2\}\],
Table 5:Variables, physical meaning, bounds, and units \(Heat exchanger\)VariablePhysical meaningLowerUpperUnitsx1WBx^\{\\mathrm\{WB\}\}\_\{1\}Heat transfer area \(AA\)150m2x2WBx^\{\\mathrm\{WB\}\}\_\{2\}Mass flow rate \(m˙\\dot\{m\}\)0\.15kg/sx3WBx^\{\\mathrm\{WB\}\}\_\{3\}Log\-mean temperature difference \(Δ​Tl​m\\Delta T\_\{lm\}\)550Kx1BBx^\{\\mathrm\{BB\}\}\_\{1\}Molecular descriptor \(m1m\_\{1\}\)50200\(unitless\)x2BBx^\{\\mathrm\{BB\}\}\_\{2\}Molecular descriptor \(m2m\_\{2\}\)01\(unitless\)y1y\_\{1\}Heat capacity \(cpc\_\{p\}\)N/AN/AkJ/kg⋅\\cdotKy2y\_\{2\}Viscosity \(μ\\mu\)N/AN/APa⋅\\cdotsy3y\_\{3\}Thermal conductivity \(kk\)N/AN/AW/m⋅\\cdotKwherey1y\_\{1\}is the heat capacity \[kJ/kg⋅\\cdotK\],y2y\_\{2\}is the viscosity \[Pa⋅\\cdots\], andy3y\_\{3\}is the thermal conductivity \[W/m⋅\\cdotK\]\. The constraintg1g\_\{1\}ensures the heat dutyQ=U⋅A⋅Δ​Tl​m≥100Q=U\\cdot A\\cdot\\Delta T\_\{lm\}\\geq 100kW is satisfied, andg2g\_\{2\}enforces turbulent flow with Reynolds numberRe≥4000\\text\{Re\}\\geq 4000\. The global minimum cost given the parameters in Table[5](https://arxiv.org/html/2608.03045#A1.T5)is1000\.01000\.0withxWB⁣⋆=\[1\.0,0\.1,19\.3\]⊤x^\{\\mathrm\{WB\}\\star\}=\[1\.0,0\.1,19\.3\]^\{\\top\},xBB⁣⋆=\[50\.0,0\.0\]⊤x^\{\\mathrm\{BB\}\\star\}=\[50\.0,0\.0\]^\{\\top\}\.

### A\.8Pressure Swing Adsorption with Novel Adsorbent

We consider the design of a pressure swing adsorption \(PSA\) cycle for gas separation, a problem that naturally exhibits grey\-box structure\[[3](https://arxiv.org/html/2608.03045#bib.bib3),[49](https://arxiv.org/html/2608.03045#bib.bib49),[50](https://arxiv.org/html/2608.03045#bib.bib50)\]\. The adsorbent properties—maximum loadingqmaxq\_\{\\max\}, adsorption equilibrium constant for the target componentKAK\_\{A\}, and equilibrium constant for the competing speciesKBK\_\{B\}—are computed from molecular simulations\[[51](https://arxiv.org/html/2608.03045#bib.bib51),[52](https://arxiv.org/html/2608.03045#bib.bib52)\]that serve as black\-box functions of adsorbent structural parameters \(pore sizeσ\\sigmaand interaction energyϵ\\epsilon\)\[[53](https://arxiv.org/html/2608.03045#bib.bib53)\]\.

The white\-box model captures the PSA cycle performance through Langmuir isotherm\-based working capacity calculations\[[54](https://arxiv.org/html/2608.03045#bib.bib54),[55](https://arxiv.org/html/2608.03045#bib.bib55)\]\. The working capacityΔ​q\\Delta qrepresents the difference in loading between high\-pressure adsorption and low\-pressure desorption steps\[[56](https://arxiv.org/html/2608.03045#bib.bib56),[50](https://arxiv.org/html/2608.03045#bib.bib50)\]\. Selectivityα=KA/KB\\alpha=K\_\{A\}/K\_\{B\}determines separation performance\[[57](https://arxiv.org/html/2608.03045#bib.bib57)\], while productivity is defined as working capacity divided by adsorption time\[[58](https://arxiv.org/html/2608.03045#bib.bib58)\]\. The objective minimizes the trade\-off between maximizing productivity and minimizing compression costs, where isothermal compression work scales withln⁡\(PH/PL\)\\ln\(P\_\{H\}/P\_\{L\}\)\[[59](https://arxiv.org/html/2608.03045#bib.bib59)\]\.

The problem can be formulated as:

minPH,PL,tads,σ,ϵ\\displaystyle\\min\\limits\_\{P\_\{H\},P\_\{L\},t\_\{\\text\{ads\}\},\\sigma,\\epsilon\}−Prod\+0\.1⋅ln⁡\(PHPL\)\\displaystyle~~\-\\text\{Prod\}\+0\.1\\cdot\\ln\\left\(\\frac\{P\_\{H\}\}\{P\_\{L\}\}\\right\)\(24\)s\.t\.α=y2y3,\\displaystyle~~\\alpha=\\frac\{y\_\{2\}\}\{y\_\{3\}\},Δ​q=y1⋅\(y2⋅PH1\+y2⋅PH−y2⋅PL1\+y2⋅PL\),\\displaystyle~~\\Delta q=y\_\{1\}\\cdot\\left\(\\frac\{y\_\{2\}\\cdot P\_\{H\}\}\{1\+y\_\{2\}\\cdot P\_\{H\}\}\-\\frac\{y\_\{2\}\\cdot P\_\{L\}\}\{1\+y\_\{2\}\\cdot P\_\{L\}\}\\right\),Prod=Δ​qtads,\\displaystyle~~\\text\{Prod\}=\\frac\{\\Delta q\}\{t\_\{\\text\{ads\}\}\},g1​\(x,y\):=3−α≤0,\\displaystyle~~g\_\{1\}\(x,y\):=3\-\\alpha\\leq 0,g2​\(x,y\):=0\.95−α⋅PH1\+α⋅PH≤0,\\displaystyle~~g\_\{2\}\(x,y\):=0\.95\-\\frac\{\\alpha\\cdot P\_\{H\}\}\{1\+\\alpha\\cdot P\_\{H\}\}\\leq 0,y1=f1BB​\(σ,ϵ\):=5\.0⋅exp⁡\(−\(σ−6\)24\)⋅\(1\+0\.1​ϵ\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(\\sigma,\\epsilon\):=5\.0\\cdot\\exp\\left\(\-\\frac\{\(\\sigma\-6\)^\{2\}\}\{4\}\\right\)\\cdot\\left\(1\+0\.1\\epsilon\\right\),y2=f2BB​\(σ,ϵ\):=0\.5⋅exp⁡\(ϵ−105\)⋅\(1−0\.05​\(σ−5\)2\),\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(\\sigma,\\epsilon\):=0\.5\\cdot\\exp\\left\(\\frac\{\\epsilon\-10\}\{5\}\\right\)\\cdot\\left\(1\-0\.05\(\\sigma\-5\)^\{2\}\\right\),y3=f3BB​\(σ,ϵ\):=0\.1⋅exp⁡\(ϵ−158\)⋅\(1\+0\.03​\(σ−7\)2\),\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(\\sigma,\\epsilon\):=0\.1\\cdot\\exp\\left\(\\frac\{\\epsilon\-15\}\{8\}\\right\)\\cdot\\left\(1\+0\.03\(\\sigma\-7\)^\{2\}\\right\),PH∈\[2,10\],PL∈\[0\.1,1\],tads∈\[10,120\],\\displaystyle~~P\_\{H\}\\in\[2,10\],\\quad P\_\{L\}\\in\[0\.1,1\],\\quad t\_\{\\text\{ads\}\}\\in\[10,120\],σ∈\[3,10\],ϵ∈\[5,25\],\\displaystyle~~\\sigma\\in\[3,10\],\\quad\\epsilon\\in\[5,25\],x=\[xWB,xBB\],xWB=\[PH,PL,tads\],xBB=\[σ,ϵ\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[P\_\{H\},P\_\{L\},t\_\{\\text\{ads\}\}\],\\quad x^\{\\mathrm\{BB\}\}=\[\\sigma,\\epsilon\],The black\-box functionsfBBf^\{\\mathrm\{BB\}\}model the adsorbent properties as functions of the Lennard\-Jones parameters: pore sizeσ\\sigma\[Å\] and interaction energyϵ\\epsilon\[kJ/mol\]\[[53](https://arxiv.org/html/2608.03045#bib.bib53),[52](https://arxiv.org/html/2608.03045#bib.bib52)\]\. These relationships capture how molecular\-scale adsorbent design affects macroscopic separation performance\[[51](https://arxiv.org/html/2608.03045#bib.bib51)\]\.

Table 6:Variables, physical meaning, bounds, and units \(Pressure Swing Adsorption\)The constraintg1g\_\{1\}ensures minimum selectivityα≥3\\alpha\\geq 3for effective separation\[[57](https://arxiv.org/html/2608.03045#bib.bib57),[56](https://arxiv.org/html/2608.03045#bib.bib56)\], andg2g\_\{2\}enforces minimum product purity based on the Langmuir isotherm equilibrium\[[54](https://arxiv.org/html/2608.03045#bib.bib54),[50](https://arxiv.org/html/2608.03045#bib.bib50)\]\. The global minimum given the parameters in Table[6](https://arxiv.org/html/2608.03045#A1.T6)is−0\.6433\-0\.6433withxWB⁣⋆=\[3\.97,0\.1,10\.0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[3\.97,0\.1,10\.0\]^\{\\top\},xBB⁣⋆=\[6\.04,19\.55\]⊤x^\{\\mathrm\{BB\}\\star\}=\[6\.04,19\.55\]^\{\\top\}\.

### A\.9Batch Reactor with Catalyst Optimization

We consider the optimization of a batch reactor for consecutive reactionsA→B→CA\\to B\\to C, where the objective is to maximize the yield of the intermediate productBB\[[60](https://arxiv.org/html/2608.03045#bib.bib60)\]\. This classic selectivity problem\[[61](https://arxiv.org/html/2608.03045#bib.bib61)\]gains a multi\-scale dimension when catalyst properties must be co\-optimized with operating conditions\.

The white\-box model consists of standard Arrhenius kinetics\[[62](https://arxiv.org/html/2608.03045#bib.bib62)\]governing the batch reactor mass balances, which admit analytical solutions for first\-order consecutive reactions\. TemperatureTT, batch timetft\_\{f\}, and initial concentrationCA,0C\_\{A,0\}are the white\-box decision variables\.

The black\-box model implements a Sabatier volcano relationship\[[63](https://arxiv.org/html/2608.03045#bib.bib63),[64](https://arxiv.org/html/2608.03045#bib.bib64),[65](https://arxiv.org/html/2608.03045#bib.bib65)\]that captures how catalyst activity depends on binding energyEbE\_\{b\}\. The volcano principle states that optimal catalysts bind reactants neither too strongly \(limiting desorption\) nor too weakly \(limiting adsorption\), with peak activity at an intermediate binding energy\. The catalyst particle sizeddprovides an additional degree of freedom that affects the pre\-exponential factor through available surface area\.

The problem can be formulated as:

minT,tf,CA,0,Eb,d\\displaystyle\\min\\limits\_\{T,t\_\{f\},C\_\{A,0\},E\_\{b\},d\}−CB\+0\.001⋅T⋅tf\\displaystyle~~\-C\_\{B\}\+0\.001\\cdot T\\cdot t\_\{f\}\(25\)s\.t\.k1′=y1⋅exp⁡\(−2000T\),\\displaystyle~~k\_\{1\}^\{\\prime\}=y\_\{1\}\\cdot\\exp\\left\(\\frac\{\-2000\}\{T\}\\right\),k2′=y2⋅exp⁡\(−2400T\),\\displaystyle~~k\_\{2\}^\{\\prime\}=y\_\{2\}\\cdot\\exp\\left\(\\frac\{\-2400\}\{T\}\\right\),CA=CA,0⋅exp⁡\(−k1′⋅tf\),\\displaystyle~~C\_\{A\}=C\_\{A,0\}\\cdot\\exp\(\-k\_\{1\}^\{\\prime\}\\cdot t\_\{f\}\),CB=CA,0⋅k1′k2′−k1′⋅\(exp⁡\(−k1′⋅tf\)−exp⁡\(−k2′⋅tf\)\),\\displaystyle~~C\_\{B\}=C\_\{A,0\}\\cdot\\frac\{k\_\{1\}^\{\\prime\}\}\{k\_\{2\}^\{\\prime\}\-k\_\{1\}^\{\\prime\}\}\\cdot\\left\(\\exp\(\-k\_\{1\}^\{\\prime\}\\cdot t\_\{f\}\)\-\\exp\(\-k\_\{2\}^\{\\prime\}\\cdot t\_\{f\}\)\\right\),g1​\(x,y\):=0\.8−\(1−CACA,0\)≤0,\\displaystyle~~g\_\{1\}\(x,y\):=0\.8\-\\left\(1\-\\frac\{C\_\{A\}\}\{C\_\{A,0\}\}\\right\)\\leq 0,g2​\(x,y\):=0\.7−CBCA,0−CA≤0,\\displaystyle~~g\_\{2\}\(x,y\):=0\.7\-\\frac\{C\_\{B\}\}\{C\_\{A,0\}\-C\_\{A\}\}\\leq 0,y1=f1BB​\(Eb,d\):=10⋅d⋅exp⁡\(−\(Eb\+1\)20\.25\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(E\_\{b\},d\):=10\\cdot d\\cdot\\exp\\left\(\-\\frac\{\(E\_\{b\}\+1\)^\{2\}\}\{0\.25\}\\right\),y2=f2BB​\(Eb,d\):=0\.5⋅d⋅exp⁡\(−\(Eb\+0\.5\)20\.5\),\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(E\_\{b\},d\):=0\.5\\cdot d\\cdot\\exp\\left\(\-\\frac\{\(E\_\{b\}\+0\.5\)^\{2\}\}\{0\.5\}\\right\),y3=f3BB​\(Eb,d\):=exp⁡\(−2⋅Eb−1\),\(not used in white\-box equations\)\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(E\_\{b\},d\):=\\exp\\left\(\-2\\cdot E\_\{b\}\-1\\right\),\\quad\\text\{\(not used in white\-box equations\)\}T∈\[300,500\],tf∈\[0\.5,10\],CA,0∈\[0\.5,5\],\\displaystyle~~T\\in\[300,500\],\\quad t\_\{f\}\\in\[0\.5,10\],\\quad C\_\{A,0\}\\in\[0\.5,5\],Eb∈\[−2,0\],d∈\[0\.5,2\],\\displaystyle~~E\_\{b\}\\in\[\-2,0\],\\quad d\\in\[0\.5,2\],x=\[xWB,xBB\],xWB=\[T,tf,CA,0\],xBB=\[Eb,d\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[T,t\_\{f\},C\_\{A,0\}\],\\quad x^\{\\mathrm\{BB\}\}=\[E\_\{b\},d\],The Gaussian functions inf1BBf^\{\\mathrm\{BB\}\}\_\{1\}andf2BBf^\{\\mathrm\{BB\}\}\_\{2\}encode the volcano\-shaped activity curves:k1k\_\{1\}peaks atEb=−1E\_\{b\}=\-1eV whilek2k\_\{2\}peaks atEb=−0\.5E\_\{b\}=\-0\.5eV, creating a selectivity landscape where catalyst design can favor the desired reaction over the consecutive side reaction\. The outputy3y\_\{3\}does not appear in any white\-box equation or constraint; it is included intentionally to test robustness to irrelevant black\-box outputs, since in practice molecular simulations often return quantities not all of which enter the process model\.

Table 7:Variables, physical meaning, bounds, and units \(Batch Reactor\)The constraintg1g\_\{1\}ensures minimum conversionXA≥0\.8X\_\{A\}\\geq 0\.8, andg2g\_\{2\}enforces selectivitySB=CB/\(CA,0−CA\)≥0\.7S\_\{B\}=C\_\{B\}/\(C\_\{A,0\}\-C\_\{A\}\)\\geq 0\.7, requiring that at least 70% of converted reactant forms the desired intermediate\. The global minimum given the parameters in Table[7](https://arxiv.org/html/2608.03045#A1.T7)is−1\.7488\-1\.7488withxWB⁣⋆=\[500\.0,4\.394,5\.0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[500\.0,4\.394,5\.0\]^\{\\top\},xBB⁣⋆=\[−1\.01,2\.0\]⊤x^\{\\mathrm\{BB\}\\star\}=\[\-1\.01,2\.0\]^\{\\top\}\.

### A\.10Extractive Distillation with Novel Solvent

We consider the integrated design of an extractive distillation column and entrainer selection, a canonical example of simultaneous solvent and process optimization\[[66](https://arxiv.org/html/2608.03045#bib.bib66),[67](https://arxiv.org/html/2608.03045#bib.bib67)\]\. The white\-box model employs classical shortcut methods: the Fenske equation\[[68](https://arxiv.org/html/2608.03045#bib.bib68)\]determines minimum stagesNminN\_\{\\min\}from the relative volatility and separation specifications, the Underwood equations\[[69](https://arxiv.org/html/2608.03045#bib.bib69),[70](https://arxiv.org/html/2608.03045#bib.bib70)\]yield minimum refluxRminR\_\{\\min\}, and operating constraints follow from the Gilliland correlation\[[71](https://arxiv.org/html/2608.03045#bib.bib71),[72](https://arxiv.org/html/2608.03045#bib.bib72)\]\.

The black\-box model represents COSMO\-RS or other molecular thermodynamic predictions that map solvent molecular descriptors—here parameterized by Hansen solubility parameters for hydrogen bondingδH\\delta\_\{H\}and polar interactionsδP\\delta\_\{P\}—to the modified relative volatilityα12mod\\alpha\_\{12\}^\{\\text\{mod\}\}in the presence of the entrainer\[[73](https://arxiv.org/html/2608.03045#bib.bib73)\]\. The entrainer also affects operating costs through its densityρ\\rho\(determining column sizing\) and viscosityμ\\mu\(affecting tray hydraulics and heat transfer\)\.

The problem can be formulated as:

minR,N,S/F,δH,δP\\displaystyle\\min\\limits\_\{R,N,S/F,\\delta\_\{H\},\\delta\_\{P\}\}1000⋅N\+500⋅R\+200⋅\(S/F\)⋅y3\\displaystyle~~1000\\cdot N\+500\\cdot R\+200\\cdot\(S/F\)\\cdot y\_\{3\}\(26\)s\.t\.Rmin=1y1−1⋅\(xDxF−y1⋅1−xD1−xF\),\\displaystyle~~R\_\{\\min\}=\\frac\{1\}\{y\_\{1\}\-1\}\\cdot\\left\(\\frac\{x\_\{D\}\}\{x\_\{F\}\}\-y\_\{1\}\\cdot\\frac\{1\-x\_\{D\}\}\{1\-x\_\{F\}\}\\right\),Nmin=ln⁡\(xD​\(1−xB\)xB​\(1−xD\)\)ln⁡\(y1\),\\displaystyle~~N\_\{\\min\}=\\frac\{\\ln\\left\(\\frac\{x\_\{D\}\(1\-x\_\{B\}\)\}\{x\_\{B\}\(1\-x\_\{D\}\)\}\\right\)\}\{\\ln\(y\_\{1\}\)\},g1​\(x,y\):=1\.2⋅Rmin−R≤0,\\displaystyle~~g\_\{1\}\(x,y\):=1\.2\\cdot R\_\{\\min\}\-R\\leq 0,g2​\(x,y\):=1\.5⋅Nmin−N≤0,\\displaystyle~~g\_\{2\}\(x,y\):=1\.5\\cdot N\_\{\\min\}\-N\\leq 0,g3​\(x,y\):=\(S/F\)⋅y2−5000≤0,\\displaystyle~~g\_\{3\}\(x,y\):=\(S/F\)\\cdot y\_\{2\}\-5000\\leq 0,y1=f1BB​\(δH,δP\):=2\.5\+1\.0⋅sin⁡\(π​\(δH−10\)20\)⋅cos⁡\(π​\(δP−10\)15\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(\\delta\_\{H\},\\delta\_\{P\}\):=2\.5\+1\.0\\cdot\\sin\\left\(\\frac\{\\pi\(\\delta\_\{H\}\-10\)\}\{20\}\\right\)\\cdot\\cos\\left\(\\frac\{\\pi\(\\delta\_\{P\}\-10\)\}\{15\}\\right\),y2=f2BB​\(δH,δP\):=800\+50⋅δH−20⋅δP,\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(\\delta\_\{H\},\\delta\_\{P\}\):=800\+50\\cdot\\delta\_\{H\}\-20\\cdot\\delta\_\{P\},y3=f3BB​\(δH,δP\):=0\.5\+0\.1⋅δH\+0\.05⋅δP2/100,\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(\\delta\_\{H\},\\delta\_\{P\}\):=0\.5\+0\.1\\cdot\\delta\_\{H\}\+0\.05\\cdot\\delta\_\{P\}^\{2\}/100,R∈\[1,10\],N∈\[5,50\],S/F∈\[0\.5,5\],\\displaystyle~~R\\in\[1,10\],\\quad N\\in\[5,50\],\\quad S/F\\in\[0\.5,5\],δH∈\[5,25\],δP∈\[5,20\],\\displaystyle~~\\delta\_\{H\}\\in\[5,25\],\\quad\\delta\_\{P\}\\in\[5,20\],x=\[xWB,xBB\],xWB=\[R,N,S/F\],xBB=\[δH,δP\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[R,N,S/F\],\\quad x^\{\\mathrm\{BB\}\}=\[\\delta\_\{H\},\\delta\_\{P\}\],wherexD=0\.99x\_\{D\}=0\.99\(distillate purity\),xF=0\.5x\_\{F\}=0\.5\(feed composition\), andxB=0\.01x\_\{B\}=0\.01\(bottoms impurity\)\.

Table 8:Variables, physical meaning, bounds, and units \(Extractive Distillation\)The constraintg1g\_\{1\}ensures the reflux ratio exceeds the minimum by at least 20% for operational flexibility,g2g\_\{2\}requires sufficient stages above the Fenske minimum, andg3g\_\{3\}limits solvent inventory based on density to control capital costs\. The global minimum cost for the parameters in Table[8](https://arxiv.org/html/2608.03045#A1.T8)is11,75811\{,\}758withxWB⁣⋆=\[1\.0,11\.0,0\.5\]⊤x^\{\\mathrm\{WB\}\\star\}=\[1\.0,11\.0,0\.5\]^\{\\top\},xBB⁣⋆=\[19\.8,9\.99\]⊤x^\{\\mathrm\{BB\}\\star\}=\[19\.8,9\.99\]^\{\\top\}\.

### A\.11Multi\-Effect Evaporator with Phase\-Change Material

We consider a prototype test problem for the design of a multi\-effect evaporator \(MEE\) system integrated with a novel phase\-change material \(PCM\) as the heat transfer medium\. Multi\-effect evaporator synthesis follows established energy balance formulations\[[74](https://arxiv.org/html/2608.03045#bib.bib74)\], where steam economy improves with additional effects but with diminishing returns due to tighter temperature driving forces\. The efficiency degradation across effects is captured by the factor0\.9n−10\.9^\{n\-1\}, consistent with typical industrial values of 0\.7–0\.9 pounds of evaporation per pound of steam per effect\[[75](https://arxiv.org/html/2608.03045#bib.bib75)\]\.

The white\-box model includes energy balances relating the heat transfer fluid \(HTF\) flow rate and properties to evaporation capacity\[[76](https://arxiv.org/html/2608.03045#bib.bib76)\]\. The latent heat of vaporization for water \(2260 kJ/kg\) determines the evaporation rate from available thermal energy\. The number of effectsnn, HTF mass flow ratem˙HTF\\dot\{m\}\_\{\\text\{HTF\}\}, and temperature drop per effectΔ​Teffect\\Delta T\_\{\\text\{effect\}\}are the white\-box decision variables\.

The black\-box model represents molecular simulation or group\-contribution predictions of PCM thermophysical properties\[[77](https://arxiv.org/html/2608.03045#bib.bib77),[78](https://arxiv.org/html/2608.03045#bib.bib78)\]as functions of melting temperatureTmT\_\{m\}and a molecular structure parameterλ\\lambdathat controls the length of alkyl chains \(affecting both latent heat and heat capacity\)\. The key outputs are latent heat of fusionΔ​Hfus\\Delta H\_\{\\text\{fus\}\}\(typical range 150–300 kJ/kg for organic PCMs\), liquid heat capacitycpliqc\_\{p\}^\{\\text\{liq\}\}\(typically 2–3 kJ/\(kg⋅\\cdotK\)\), and the actual melting temperatureTmeltT\_\{\\text\{melt\}\}\. The functional forms infBBf^\{\\mathrm\{BB\}\}are intended to capture qualitative structure\-property relationships rather than precise physical models\.

The problem can be formulated as:

minn,m˙HTF,Δ​Teffect,Tm,λ\\displaystyle\\min\\limits\_\{n,\\dot\{m\}\_\{\\text\{HTF\}\},\\Delta T\_\{\\text\{effect\}\},T\_\{m\},\\lambda\}−m˙evap\+0\.1⋅n2\+0\.01⋅m˙HTF2\\displaystyle~~\-\\dot\{m\}\_\{\\text\{evap\}\}\+0\.1\\cdot n^\{2\}\+0\.01\\cdot\\dot\{m\}\_\{\\text\{HTF\}\}^\{2\}\(27\)s\.t\.QHTF=m˙HTF⋅\(y1\+y2⋅10\),\\displaystyle~~Q\_\{\\text\{HTF\}\}=\\dot\{m\}\_\{\\text\{HTF\}\}\\cdot\(y\_\{1\}\+y\_\{2\}\\cdot 10\),m˙evap=QHTF⋅n⋅0\.9n−12260,\\displaystyle~~\\dot\{m\}\_\{\\text\{evap\}\}=\\frac\{Q\_\{\\text\{HTF\}\}\\cdot n\\cdot 0\.9^\{n\-1\}\}\{2260\},g1​\(x,y\):=100\+n⋅Δ​Teffect−y3≤0,\\displaystyle~~g\_\{1\}\(x,y\):=100\+n\\cdot\\Delta T\_\{\\text\{effect\}\}\-y\_\{3\}\\leq 0,g2​\(x,y\):=5−m˙evap≤0,\\displaystyle~~g\_\{2\}\(x,y\):=5\-\\dot\{m\}\_\{\\text\{evap\}\}\\leq 0,y1=f1BB​\(Tm,λ\):=150⋅λ⋅\(1\+0\.2​sin⁡\(π​Tm100\)\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(T\_\{m\},\\lambda\):=150\\cdot\\lambda\\cdot\\left\(1\+0\.2\\sin\\left\(\\frac\{\\pi T\_\{m\}\}\{100\}\\right\)\\right\),y2=f2BB​\(Tm,λ\):=2\.0\+0\.5⋅λ−0\.002⋅Tm,\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(T\_\{m\},\\lambda\):=2\.0\+0\.5\\cdot\\lambda\-0\.002\\cdot T\_\{m\},y3=f3BB​\(Tm,λ\):=Tm\+10⋅sin⁡\(π​λ2\),\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(T\_\{m\},\\lambda\):=T\_\{m\}\+10\\cdot\\sin\\left\(\\frac\{\\pi\\lambda\}\{2\}\\right\),n∈\[2,6\]​\(treated as continuous; integer rounding would apply in practice\),\\displaystyle~~n\\in\[2,6\]\\text\{ \(treated as continuous; integer rounding would apply in practice\)\},m˙HTF∈\[1,20\],Δ​Teffect∈\[5,20\],\\displaystyle~~\\dot\{m\}\_\{\\text\{HTF\}\}\\in\[1,20\],\\quad\\Delta T\_\{\\text\{effect\}\}\\in\[5,20\],Tm∈\[50,150\],λ∈\[0\.5,2\],\\displaystyle~~T\_\{m\}\\in\[50,150\],\\quad\\lambda\\in\[0\.5,2\],x=\[xWB,xBB\],xWB=\[n,m˙HTF,Δ​Teffect\],xBB=\[Tm,λ\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[n,\\dot\{m\}\_\{\\text\{HTF\}\},\\Delta T\_\{\\text\{effect\}\}\],\\quad x^\{\\mathrm\{BB\}\}=\[T\_\{m\},\\lambda\],
Table 9:Variables, physical meaning, bounds, and units \(Multi\-Effect Evaporator\)SymbolDescriptionRange / ValueUnitsWhite\-box decision variablesnnNumber of evaporator effects\[2, 6\]—m˙HTF\\dot\{m\}\_\{\\text\{HTF\}\}Heat transfer fluid mass flow rate\[1, 20\]kg/sΔ​Teffect\\Delta T\_\{\\text\{effect\}\}Temperature drop per effect\[5, 20\]∘CBlack\-box decision variablesTmT\_\{m\}Target melting temperature\[50, 150\]∘Cλ\\lambdaMolecular structure parameter\[0\.5, 2\]—Black\-box outputsy1=Δ​Hfusy\_\{1\}=\\Delta H\_\{\\text\{fus\}\}Latent heat of fusion—kJ/kgy2=cpliqy\_\{2\}=c\_\{p\}^\{\\text\{liq\}\}Liquid phase heat capacity—kJ/\(kg⋅\\cdotK\)y3=Tmelty\_\{3\}=T\_\{\\text\{melt\}\}Actual melting temperature—∘CFixed parametersΔ​Hvap,water\\Delta H\_\{\\text\{vap,water\}\}Latent heat of water vaporization2260kJ/kgηeffect\\eta\_\{\\text\{effect\}\}Per\-effect efficiency factor0\.9—The constraintg1g\_\{1\}ensures temperature feasibility: the PCM melting temperature must exceed the cumulative temperature drop across all effects plus a 100∘C baseline\. The constraintg2g\_\{2\}enforces a minimum evaporation rate of 5 kg/s to meet production requirements\. The global minimum for parameters in Table[9](https://arxiv.org/html/2608.03045#A1.T9)is−1\.9587\-1\.9587withxWB⁣⋆=\[4\.09,19\.07,5\.0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[4\.09,19\.07,5\.0\]^\{\\top\},xBB⁣⋆=\[120\.47,2\.0\]⊤x^\{\\mathrm\{BB\}\\star\}=\[120\.47,2\.0\]^\{\\top\}\.

### A\.12Membrane Separation with Designed Polymer

We consider the integrated design of a membrane gas separation process and the polymer membrane material\[[79](https://arxiv.org/html/2608.03045#bib.bib79)\]\. This problem exemplifies the permeability\-selectivity trade\-off characterized by the Robeson upper bound\[[80](https://arxiv.org/html/2608.03045#bib.bib80),[81](https://arxiv.org/html/2608.03045#bib.bib81)\], which establishes that polymeric membranes exhibit a fundamental inverse relationship between gas permeability and selectivity explained by free volume theory\[[82](https://arxiv.org/html/2608.03045#bib.bib82)\]\.

The white\-box model follows the solution\-diffusion transport mechanism\[[83](https://arxiv.org/html/2608.03045#bib.bib83)\], where the permeate fluxJAJ\_\{A\}is proportional to permeability, pressure driving force, and inversely proportional to membrane thicknessδ\\delta\. The membrane areaAmA\_\{m\}, thicknessδ\\delta, and transmembrane pressure differenceΔ​p\\Delta pare the white\-box decision variables, with the objective of maximizing product recovery while minimizing capital \(area\) and operating \(compression\) costs\.

The black\-box model represents structure\-property relationships from molecular simulations or machine learning predictions\[[84](https://arxiv.org/html/2608.03045#bib.bib84)\]that map polymer descriptors—fractional free volume \(FFV\) and chain spacingdspacingd\_\{\\text\{spacing\}\}—to transport properties\. The Robeson trade\-off is encoded in the black\-box functions: increasing FFV enhances permeability but reduces selectivity, while chain spacing affects both the upper bound position and mechanical integrity\. The mechanical strength constraint reflects the requirement that the membrane must withstand the transmembrane pressure difference during module operation\[[85](https://arxiv.org/html/2608.03045#bib.bib85)\]\.

The problem can be formulated as:

minAm,δ,Δ​p,FFV,dspacing\\displaystyle\\min\\limits\_\{A\_\{m\},\\delta,\\Delta p,\\text\{FFV\},d\_\{\\text\{spacing\}\}\}−FA⋅y2\+0\.1⋅Am\+10⋅Δ​p\\displaystyle~~\-F\_\{A\}\\cdot y\_\{2\}\+0\.1\\cdot A\_\{m\}\+10\\cdot\\Delta p\(28\)s\.t\.JA=y1⋅Δ​pδ⋅10−10,\\displaystyle~~J\_\{A\}=\\frac\{y\_\{1\}\\cdot\\Delta p\}\{\\delta\}\\cdot 10^\{\-10\},FA=JA⋅Am,\\displaystyle~~F\_\{A\}=J\_\{A\}\\cdot A\_\{m\},g1​\(x,y\):=10−6−JA≤0,\\displaystyle~~g\_\{1\}\(x,y\):=10^\{\-6\}\-J\_\{A\}\\leq 0,g2​\(x,y\):=Δ​p⋅Amy3−100≤0,\\displaystyle~~g\_\{2\}\(x,y\):=\\frac\{\\Delta p\\cdot A\_\{m\}\}\{y\_\{3\}\}\-100\\leq 0,y1=f1BB​\(FFV,dspacing\):=1000⋅exp⁡\(FFV−0\.150\.05\)⋅\(1\+0\.1⋅dspacing\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(\\text\{FFV\},d\_\{\\text\{spacing\}\}\):=1000\\cdot\\exp\\left\(\\frac\{\\text\{FFV\}\-0\.15\}\{0\.05\}\\right\)\\cdot\\left\(1\+0\.1\\cdot d\_\{\\text\{spacing\}\}\\right\),y2=f2BB​\(FFV,dspacing\):=50⋅exp⁡\(−FFV−0\.150\.1\)⋅exp⁡\(−\(dspacing−5\)24\),\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(\\text\{FFV\},d\_\{\\text\{spacing\}\}\):=50\\cdot\\exp\\left\(\-\\frac\{\\text\{FFV\}\-0\.15\}\{0\.1\}\\right\)\\cdot\\exp\\left\(\-\\frac\{\(d\_\{\\text\{spacing\}\}\-5\)^\{2\}\}\{4\}\\right\),y3=f3BB​\(FFV,dspacing\):=100⋅\(0\.4−FFV\)⋅\(10−dspacing\),\\displaystyle~~y\_\{3\}=f^\{\\mathrm\{BB\}\}\_\{3\}\(\\text\{FFV\},d\_\{\\text\{spacing\}\}\):=100\\cdot\(0\.4\-\\text\{FFV\}\)\\cdot\(10\-d\_\{\\text\{spacing\}\}\),Am∈\[10,1000\],δ∈\[0\.1,10\],Δ​p∈\[1,20\],\\displaystyle~~A\_\{m\}\\in\[10,1000\],\\quad\\delta\\in\[0\.1,10\],\\quad\\Delta p\\in\[1,20\],FFV∈\[0\.1,0\.3\],dspacing∈\[3,8\],\\displaystyle~~\\text\{FFV\}\\in\[0\.1,0\.3\],\\quad d\_\{\\text\{spacing\}\}\\in\[3,8\],x=\[xWB,xBB\],xWB=\[Am,δ,Δ​p\],xBB=\[FFV,dspacing\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[A\_\{m\},\\delta,\\Delta p\],\\quad x^\{\\mathrm\{BB\}\}=\[\\text\{FFV\},d\_\{\\text\{spacing\}\}\],Refer to the parameter definitions in Table[10](https://arxiv.org/html/2608.03045#A1.T10)\.

Table 10:Variables, physical meaning, bounds, and units \(Membrane Separation\)The white\-box flux equationJA=y1⋅Δ​p/δJ\_\{A\}=y\_\{1\}\\cdot\\Delta p/\\deltafollows directly from the solution\-diffusion model\[[83](https://arxiv.org/html/2608.03045#bib.bib83)\], where permeate flux is proportional to permeability and pressure driving force, and inversely proportional to membrane thickness\. The factor10−1010^\{\-10\}converts from Barrer \(the conventional permeability unit\) to SI flux units\.

The black\-box functions are constructed to capture established structure\-property relationships from polymer physics\. The permeability functiony1y\_\{1\}follows the exponential free volume dependence derived by Cohen and Turnbull\[[86](https://arxiv.org/html/2608.03045#bib.bib86)\]and extended by Fujita\[[87](https://arxiv.org/html/2608.03045#bib.bib87)\], who showed that diffusivity \(and hence permeability\) scales asP∝exp⁡\(B⋅FFV\)P\\propto\\exp\(B\\cdot\\text\{FFV\}\)whereBBis a polymer\-penetrant parameter\. This exponential form arises because molecular transport occurs through transient gaps in the polymer matrix, with the probability of sufficiently large gaps increasing exponentially with free volume fraction\[[88](https://arxiv.org/html/2608.03045#bib.bib88)\]\. The chain spacing term provides a secondary contribution reflecting the effect of interchain distance on diffusion pathways\[[89](https://arxiv.org/html/2608.03045#bib.bib89)\]\.

The selectivity functiony2y\_\{2\}encodes the Robeson upper bound trade\-off\[[80](https://arxiv.org/html/2608.03045#bib.bib80),[81](https://arxiv.org/html/2608.03045#bib.bib81)\]: as FFV increases, permeability rises but selectivity decreases because larger free volume elements become less size\-discriminating\. Freeman\[[82](https://arxiv.org/html/2608.03045#bib.bib82)\]provided a theoretical basis showing that the slope of the permeability\-selectivity trade\-off depends only on penetrant size ratios\. The Gaussian dependence ondspacingd\_\{\\text\{spacing\}\}reflects an optimal chain spacing for size\-sieving selectivity\[[85](https://arxiv.org/html/2608.03045#bib.bib85)\]\.

The mechanical strength proxyy3y\_\{3\}captures the physical principle that increased free volume and chain spacing reduce polymer density and chain entanglement, thereby weakening mechanical integrity\[[90](https://arxiv.org/html/2608.03045#bib.bib90)\]\. While the linear functional form is a simplification, it correctly represents the inverse relationship between transport\-enhancing microstructure and mechanical robustness that constrains practical membrane design\[[91](https://arxiv.org/html/2608.03045#bib.bib91)\]\.

The constraintg1g\_\{1\}ensures minimum fluxJA≥10−6J\_\{A\}\\geq 10^\{\-6\}mol/\(m⋅2\{\}^\{2\}\\cdots\) for practical separation rates, andg2g\_\{2\}enforces mechanical stability by limiting the pressure\-area product relative to membrane strength\.

The global optimum for parameters in Table[10](https://arxiv.org/html/2608.03045#A1.T10)occurs atxWB⁣∗=\[Am,δ,Δ​p\]=\[10,0\.1,1\]x^\{\\mathrm\{WB\}\*\}=\[A\_\{m\},\\delta,\\Delta p\]=\[10,0\.1,1\]andxBB⁣∗=\[FFV,dspacing\]=\[0\.3,5\.13\]x^\{\\mathrm\{BB\}\*\}=\[\\text\{FFV\},d\_\{\\text\{spacing\}\}\]=\[0\.3,5\.13\], yieldingJ∗=10\.997J^\{\*\}=10\.997\. The optimizer drives all white\-box variables to their lower bounds—minimizing membrane area \(capital cost\) and pressure difference \(operating cost\) while maximizing flux through minimum thickness—and selects maximum FFV for permeability while tuning chain spacing to optimize the selectivity\-permeability trade\-off near the Gaussian peak atdspacing≈5d\_\{\\text\{spacing\}\}\\approx 5Å\.

### A\.13Williams\-Otto Process

We consider the Williams\-Otto process\[[92](https://arxiv.org/html/2608.03045#bib.bib92),[93](https://arxiv.org/html/2608.03045#bib.bib93)\], a benchmark chemical process consisting of a continuously stirred tank reactor \(CSTR\) with three reactions:

A\+B\\displaystyle\\text\{A\}\+\\text\{B\}→C,\\displaystyle\\to\\text\{C\},C\+B\\displaystyle\\text\{C\}\+\\text\{B\}→P\+E,\\displaystyle\\to\\text\{P\}\+\\text\{E\},P\+C\\displaystyle\\text\{P\}\+\\text\{C\}→G,\\displaystyle\\to\\text\{G\},where P is the desired product and G is waste\. We formulate this as a multi\-scale grey\-box problem where catalyst properties \(black\-box from DFT calculations\) affect kinetic selectivity, while the reactor mass balances at steady state are known white\-box equations\. The problem can be formulated as:

maxT,Ff​B,μ,ε,σc\\displaystyle\\max\\limits\_\{T,F\_\{fB\},\\mu,\\varepsilon,\\sigma\_\{c\}\}Fp​P\\displaystyle~~F\_\{pP\}\(29\)s\.t\.ki​\(T\)=aiρ​exp⁡\(−biT\),i∈\{1,2,3\},\\displaystyle~~k\_\{i\}\(T\)=\\frac\{a\_\{i\}\}\{\\rho\}\\exp\\left\(\-\\frac\{b\_\{i\}\}\{T\}\\right\),\\quad i\\in\\\{1,2,3\\\},r1=y1⋅k1​\(T\)⋅XA​XB,r2=y2⋅k2​\(T\)⋅XB​XC,r3=k3​\(T\)⋅XC​XP,\\displaystyle~~r\_\{1\}=y\_\{1\}\\cdot k\_\{1\}\(T\)\\cdot X\_\{A\}X\_\{B\},\\quad r\_\{2\}=y\_\{2\}\\cdot k\_\{2\}\(T\)\\cdot X\_\{B\}X\_\{C\},\\quad r\_\{3\}=k\_\{3\}\(T\)\\cdot X\_\{C\}X\_\{P\},0=Ff​A−μ​XA−r1⋅m,\\displaystyle~~0=F\_\{fA\}\-\\mu X\_\{A\}\-r\_\{1\}\\cdot m,0=Ff​B−μ​XB−r1⋅m−r2⋅m,\\displaystyle~~0=F\_\{fB\}\-\\mu X\_\{B\}\-r\_\{1\}\\cdot m\-r\_\{2\}\\cdot m,0=−μ​XC\+2​r1⋅m−2​r2⋅m−r3⋅m,\\displaystyle~~0=\-\\mu X\_\{C\}\+2r\_\{1\}\\cdot m\-2r\_\{2\}\\cdot m\-r\_\{3\}\\cdot m,0=−μ​XE\+2​r2⋅m,\\displaystyle~~0=\-\\mu X\_\{E\}\+2r\_\{2\}\\cdot m,0=−μ​XP\+r2⋅m−0\.5​r3⋅m,\\displaystyle~~0=\-\\mu X\_\{P\}\+r\_\{2\}\\cdot m\-0\.5r\_\{3\}\\cdot m,0=−μ​XG\+1\.5​r3⋅m,\\displaystyle~~0=\-\\mu X\_\{G\}\+1\.5r\_\{3\}\\cdot m,Fp​P=μ​XP−0\.1​μ​XE,Fw​G=μ​XG,\\displaystyle~~F\_\{pP\}=\\mu X\_\{P\}\-0\.1\\mu X\_\{E\},\\quad F\_\{wG\}=\\mu X\_\{G\},g1​\(x,y\):=Fw​G−1\.0≤0,\\displaystyle~~g\_\{1\}\(x,y\):=F\_\{wG\}\-1\.0\\leq 0,g2​\(x,y\):=0\.5−XA−XB−XC−XE−XP−XG≤0,\\displaystyle~~g\_\{2\}\(x,y\):=0\.5\-X\_\{A\}\-X\_\{B\}\-X\_\{C\}\-X\_\{E\}\-X\_\{P\}\-X\_\{G\}\\leq 0,y1=f1BB​\(ε,σc\):=exp⁡\(−\(ε−15\)250\)⋅\(1\+0\.1​\(σc−5\)\),\\displaystyle~~y\_\{1\}=f^\{\\mathrm\{BB\}\}\_\{1\}\(\\varepsilon,\\sigma\_\{c\}\):=\\exp\\left\(\-\\frac\{\(\\varepsilon\-15\)^\{2\}\}\{50\}\\right\)\\cdot\\left\(1\+0\.1\(\\sigma\_\{c\}\-5\)\\right\),y2=f2BB​\(ε,σc\):=1\.5⋅exp⁡\(−\(ε−12\)240\)⋅exp⁡\(−\(σc−6\)28\),\\displaystyle~~y\_\{2\}=f^\{\\mathrm\{BB\}\}\_\{2\}\(\\varepsilon,\\sigma\_\{c\}\):=1\.5\\cdot\\exp\\left\(\-\\frac\{\(\\varepsilon\-12\)^\{2\}\}\{40\}\\right\)\\cdot\\exp\\left\(\-\\frac\{\(\\sigma\_\{c\}\-6\)^\{2\}\}\{8\}\\right\),T∈\[400,700\]∘​R,Ff​B∈\[5,50\]​klb/h,μ∈\[50,200\]​klb/h,\\displaystyle~~T\\in\[400,700\]~^\{\\circ\}\\text\{R\},\\quad F\_\{fB\}\\in\[5,50\]~\\text\{klb/h\},\\quad\\mu\\in\[50,200\]~\\text\{klb/h\},ε∈\[5,25\]​kJ/mol,σc∈\[3,10\]​Å,\\displaystyle~~\\varepsilon\\in\[5,25\]~\\text\{kJ/mol\},\\quad\\sigma\_\{c\}\\in\[3,10\]~\\text\{\\AA \},x=\[xWB,xBB\],xWB=\[T,Ff​B,μ\],xBB=\[ε,σc\],\\displaystyle~~x=\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\],\\quad x^\{\\mathrm\{WB\}\}=\[T,F\_\{fB\},\\mu\],\\quad x^\{\\mathrm\{BB\}\}=\[\\varepsilon,\\sigma\_\{c\}\],whereXi=mi/mX\_\{i\}=m\_\{i\}/mare mass fractions andmmis total reactor mass\. Refer to the parameters in Table[11](https://arxiv.org/html/2608.03045#A1.T11)\.

Table 11:Variables, physical meaning, bounds, and units \(Williams\-Otto Process\)The black\-box outputsy1y\_\{1\}andy2y\_\{2\}represent catalyst activity factors for reactions 1 and 2, modeled as volcano\-type functions of catalyst surface energyε\\varepsilonand pore sizeσc\\sigma\_\{c\}\. The constraintg1g\_\{1\}limits waste production \(Fw​G≤1\.0F\_\{wG\}\\leq 1\.0klb/h\) andg2g\_\{2\}ensures mass balance closure \(∑iXi≥0\.5\\sum\_\{i\}X\_\{i\}\\geq 0\.5\)\. The global maximum net product flow is5\.4705\.470klb/h withxWB⁣⋆=\[581\.9,50\.0,50\.0\]⊤x^\{\\mathrm\{WB\}\\star\}=\[581\.9,50\.0,50\.0\]^\{\\top\},xBB⁣⋆=\[12\.6,6\.11\]⊤x^\{\\mathrm\{BB\}\\star\}=\[12\.6,6\.11\]^\{\\top\}\.

## Appendix BAcquisition Function Details

All BO experiments use Expected Improvement \(EI\) as the acquisition function\. We also implemented four alternatives for completeness; the EI results are reported in the main text because EI is the most widely used baseline and enables direct comparison with prior work\.

### B\.1Expected Improvement

For minimization, Expected Improvement is defined as

αEI​\(x\)=\(y^−μ​\(x\)−ξ\)​Φ​\(Z\)\+σ​\(x\)​ϕ​\(Z\),Z=y^−μ​\(x\)−ξσ​\(x\),\\alpha\_\{\\text\{EI\}\}\(x\)=\(\\hat\{y\}\-\\mu\(x\)\-\\xi\)\\,\\Phi\(Z\)\+\\sigma\(x\)\\,\\phi\(Z\),\\quad Z=\\frac\{\\hat\{y\}\-\\mu\(x\)\-\\xi\}\{\\sigma\(x\)\},\(30\)whereμ​\(x\)\\mu\(x\)andσ​\(x\)\\sigma\(x\)are the GP posterior mean and standard deviation,y^\\hat\{y\}is the best observed objective value,Φ\\Phiandϕ\\phiare the standard normal CDF and PDF, andξ≥0\\xi\\geq 0is an exploration parameter\. Whenσ​\(x\)=0\\sigma\(x\)=0, we setαEI​\(x\)=0\\alpha\_\{\\text\{EI\}\}\(x\)=0\. For maximization, the improvement direction is reversed:Z=\(μ​\(x\)−y^−ξ\)/σ​\(x\)Z=\(\\mu\(x\)\-\\hat\{y\}\-\\xi\)/\\sigma\(x\)\.

The exploration parameterξ\\xicontrols the trade\-off between exploitation \(ξ=0\\xi=0, selecting points where the GP mean is low\) and exploration \(ξ\>0\\xi\>0, favoring regions of high uncertainty even if the mean is not promising\)\. We sweepξ∈\{0\.001,0\.01,0\.05,0\.1,0\.2,0\.5,1\.0\}\\xi\\in\\\{0\.001,0\.01,0\.05,0\.1,0\.2,0\.5,1\.0\\\}in all experiments\.

### B\.2Probability of Improvement

Probability of Improvement measures the probability that a candidate point improves over the current best:

αPI​\(x\)=Φ​\(Z\),Z=y^−μ​\(x\)−ξσ​\(x\)\.\\alpha\_\{\\text\{PI\}\}\(x\)=\\Phi\(Z\),\\quad Z=\\frac\{\\hat\{y\}\-\\mu\(x\)\-\\xi\}\{\\sigma\(x\)\}\.\(31\)PI is a simpler alternative to EI that ignores the magnitude of improvement\. It tends to exploit more aggressively than EI for the sameξ\\xi\.

### B\.3Lower Confidence Bound

The Lower Confidence Bound \(LCB\) acquisition function for minimization is

αLCB​\(x\)=−\(μ​\(x\)−κ​σ​\(x\)\),\\alpha\_\{\\text\{LCB\}\}\(x\)=\-\\bigl\(\\mu\(x\)\-\\kappa\\,\\sigma\(x\)\\bigr\),\(32\)whereκ\>0\\kappa\>0controls exploration\. We negate LCB so that all acquisition functions follow the convention that higher values are better \(i\.e\., the next point is selected byarg⁡maxx⁡α​\(x\)\\arg\\max\_\{x\}\\alpha\(x\)\)\. For maximization, the Upper Confidence BoundαUCB​\(x\)=μ​\(x\)\+κ​σ​\(x\)\\alpha\_\{\\text\{UCB\}\}\(x\)=\\mu\(x\)\+\\kappa\\,\\sigma\(x\)is used\. We setκ=2\.0\\kappa=2\.0by default\.

### B\.4Modified Watson\-Barnes 2

The modified Watson\-Barnes 2 \(mWB2\) acquisition function was proposed in the COBALT framework\[[6](https://arxiv.org/html/2608.03045#bib.bib6)\]:

αmWB2​\(x\)=sn⋅αEI​\(x\)−μ​\(x\),sn=−3​\(1−nN\),\\alpha\_\{\\text\{mWB2\}\}\(x\)=s\_\{n\}\\cdot\\alpha\_\{\\text\{EI\}\}\(x\)\-\\mu\(x\),\\quad s\_\{n\}=\-3\\left\(1\-\\frac\{n\}\{N\}\\right\),\(33\)wherennis the current iteration andNNis the total budget\. The dynamic scaling factorsns\_\{n\}starts large \(emphasizing EI\-driven exploration\) and decreases toward zero \(shifting weight to the GP mean for exploitation\)\. For maximization, the sign ofμ​\(x\)\\mu\(x\)is reversed\.

### B\.5Thompson Sampling

Thompson Sampling draws a sample functionf~\\tilde\{f\}from the GP posterior and selects the point that optimizes the sample:

xnext=arg⁡minx∈𝒳⁡f~​\(x\),f~∼𝒢​𝒫​\(μ,k\)\.x\_\{\\text\{next\}\}=\\arg\\min\\limits\_\{x\\in\\mathcal\{X\}\}\\tilde\{f\}\(x\),\\quad\\tilde\{f\}\\sim\\mathcal\{GP\}\(\\mu,k\)\.\(34\)In practice, we evaluatef~\\tilde\{f\}at a discrete set of candidate points by drawing from the multivariate normal𝒩​\(μ​\(Xgrid\),K​\(Xgrid,Xgrid\)\)\\mathcal\{N\}\(\\mu\(X\_\{\\text\{grid\}\}\),K\(X\_\{\\text\{grid\}\},X\_\{\\text\{grid\}\}\)\)\. Exploration arises naturally from the randomness of the posterior draw rather than from an explicit exploration parameter\.

### B\.6Acquisition Optimization

For all acquisition functions, we maximizeα​\(x\)\\alpha\(x\)over a random grid of 1,000 candidate points drawn uniformly from the search domain \(either𝒳BB\\mathcal\{X\}^\{\\mathrm\{BB\}\}for bilevel BO or𝒳WB×𝒳BB\\mathcal\{X\}^\{\\mathrm\{WB\}\}\\times\\mathcal\{X\}^\{\\mathrm\{BB\}\}for black\-box BO\)\. This grid search approach avoids the gradient\-based optimization of the acquisition function that can be problematic in low dimensions with multi\-modal acquisition landscapes\. The candidate with the highest acquisition value is selected as the next evaluation point\.

## Appendix CGaussian Process Configuration

Both BO solvers usescikit\-learn’sGaussianProcessRegressorwith identical configuration\. The kernel is a product of a constant kernel and an ARD RBF kernel, plus additive white noise:

k​\(x,x′\)=σf2​∏d=1Dexp⁡\(−\(xd−xd′\)22​ℓd2\)\+σn2​δ​\(x,x′\),k\(x,x^\{\\prime\}\)=\\sigma\_\{f\}^\{2\}\\prod\_\{d=1\}^\{D\}\\exp\\\!\\left\(\-\\frac\{\(x\_\{d\}\-x\_\{d\}^\{\\prime\}\)^\{2\}\}\{2\\ell\_\{d\}^\{2\}\}\\right\)\+\\sigma\_\{n\}^\{2\}\\,\\delta\(x,x^\{\\prime\}\),\(35\)whereDDis the input dimension \(nBBn\_\{\\mathrm\{BB\}\}for bilevel BO,nWB\+nBBn\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}for black\-box BO\),σf2\\sigma\_\{f\}^\{2\}is the signal variance,ℓd\\ell\_\{d\}is the length scale for dimensiondd\(automatic relevance determination\), andσn2\\sigma\_\{n\}^\{2\}is the noise variance\.

Table[12](https://arxiv.org/html/2608.03045#A3.T12)summarizes the hyperparameter bounds and optimization settings\.

Table 12:GP hyperparameter configuration\.The bilevel BO surrogate fits a GP overnBB=1n\_\{\\mathrm\{BB\}\}=1–2 input dimensions, while the black\-box BO surrogate fits overnWB\+nBB=2n\_\{\\mathrm\{WB\}\}\+n\_\{\\mathrm\{BB\}\}=2–5 dimensions\. Since GP fitting scales asO​\(N3\)O\(N^\{3\}\)in the number of observationsNN\(which is the same for both methods\) and the per\-prediction cost isO​\(N2​D\)O\(N^\{2\}D\), the lower dimensionality of bilevel BO reduces prediction cost but does not change the dominantO​\(N3\)O\(N^\{3\}\)training cost\. In practice, withN≤250N\\leq 250andD≤5D\\leq 5, GP fitting is fast for both methods \(typically<<0\.1 seconds per iteration\)\.

## Appendix DBH Baseline Details

The black\-box BH baseline solves the full problem \([4](https://arxiv.org/html/2608.03045#S2.E4)\) directly using SciPy’sbasinhoppingover the full variable space\[xWB,xBB\]\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]with an SLSQP local minimizer, a Metropolis temperature ofT=100T=100, full\-range stepsize \(uniform perturbation clipped to variable bounds\), and 1,000 Basin\-Hopping iterations\.

### D\.1Constraint Handling

Constraints are handled by the SLSQP local minimizer within Basin\-Hopping, which enforcesg​\(xWB,y,xBB\)≤0g\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\\leq 0as inequality constraints at each local minimization\. After optimization, feasibility is verified: solutions satisfying all constraints to within10−610^\{\-6\}are considered feasible\. If no feasible solution is found, the solver returns the solution with minimum constraint violation\.

In contrast, both BO methods track feasibility explicitly\. For black\-box BO, constraints are handled via a penalty of10610^\{6\}added to the objective when anygi​\(x\)\>0g\_\{i\}\(x\)\>0\. For bilevel BO, the inner BH solver checks constraint satisfaction and returns a penalty value when the inner problem is infeasible for a givenyy\.

### D\.2Sample Efficiency Comparison

Table[13](https://arxiv.org/html/2608.03045#A4.T13)compares the number of black\-box evaluations across solvers\. The BH baseline evaluatesfBBf^\{\\mathrm\{BB\}\}at every local minimization within Basin\-Hopping, resulting in substantially morefBBf^\{\\mathrm\{BB\}\}evaluations than either BO method\.

Table 13:Black\-box evaluations \(fBBf^\{\\mathrm\{BB\}\}calls\) per solver\.All three BO methods callfBBf^\{\\mathrm\{BB\}\}exactly once per iteration, so the total sample budget isninit\+Nn\_\{\\text\{init\}\}\+NwhereN=200N=200is the number of BO iterations\. Withninit=5n\_\{\\text\{init\}\}=5\(the smallest setting\), this gives 205 evaluations\. Both bilevel BO variants additionally solve an inner optimization problem at each iteration, but this uses only the white\-box modelfWBf^\{\\mathrm\{WB\}\}\(which is cheap\) and does not require additional calls tofBBf^\{\\mathrm\{BB\}\}\.

### D\.3Failure Modes

The BH baseline fails on Rastrigin, where multimodality in the full combined space challenges even global solvers: the BH baseline achieves a final regret of 3\.38, worse than all BO methods\. It also performs poorly on Distillation, where the tight constraint set and complex objective landscape limit convergence\. On all other problems, the BH baseline approaches the global optimum but requires substantially morefBBf^\{\\mathrm\{BB\}\}evaluations than the BO methods\.

## Appendix EInner\-Loop and Full\-Space Solver Benchmarks

Basin\-Hopping \(BH\) was selected as the unified global solver for both the inner\-loop white\-box subproblem and the full\-space NLP baseline\. To justify this choice, we benchmarked five solvers across all 13 problems in both contexts: multi\-start SLSQP with 50 Latin hypercube restarts \(SLSQP\-50\), differential evolution \(DE\)\[[94](https://arxiv.org/html/2608.03045#bib.bib94)\], simplicial homology global optimization \(SHGO\)\[[95](https://arxiv.org/html/2608.03045#bib.bib95)\], dual annealing \(DA\)\[[96](https://arxiv.org/html/2608.03045#bib.bib96)\], and Basin\-Hopping with 100 iterations \(BH\)\. We additionally ran a high\-iteration Basin\-Hopping variant with 1,000 iterations \(BH\-1000\) to assess convergence with extended computation\.

### E\.1Inner\-Loop Solver Comparison

Table[14](https://arxiv.org/html/2608.03045#A5.T14)reports regret \(distance from verified global optimum\) for each solver on the inner white\-box subproblem, wherexBBx^\{\\mathrm\{BB\}\}is fixed and onlyxWBx^\{\\mathrm\{WB\}\}is optimized\. All five solvers have access to analytical gradients ofJJandggwith respect toxWBx^\{\\mathrm\{WB\}\}\.

Table 14:Inner\-loop solver regret\. A dash \(—\) indicates the solver returned an infeasible solution\. All feasible regret values are near machine precision \(<10−7<10^\{\-7\}\) except for DA, which fails to converge on six problems\.SLSQP\-50 achieves the lowest regret on 7 of 13 problems, leveraging gradient information most directly\. However, it is inapplicable to the derivative\-free full\-space problem \(Section[E\.2](https://arxiv.org/html/2608.03045#A5.SS2)\)\. SHGO and DA suffer from feasibility failures: SHGO returns infeasible solutions on 2 problems \(Small\-Feasible\-Region\-1 and Distillation\), while DA fails on 6 problems \(Toy\-Hydrology, Rosen\-Suzuki, Batch\-Reactor, Distillation, Evaporator, and Williams\-Otto\)\. DE is reliable \(12/13 feasible\) but misses feasibility on one problem\. BH and BH\-1000 both achieve feasibility on all 13 problems\. The regret differences among feasible solvers are negligible in practice—on the order of10−1410^\{\-14\}to10−810^\{\-8\}—well below any meaningful engineering threshold\.

Table[15](https://arxiv.org/html/2608.03045#A5.T15)reports the corresponding wall times\. SLSQP\-50 and DA are the fastest \(median<0\.1<0\.1s\), while SHGO is the slowest \(median 3\.1 s\), with several problems exceeding 20 s\. BH is fast at 100 iterations \(median 0\.20 s\) and remains practical at 1,000 iterations \(median 0\.45 s\), except for Williams\-Otto \(30\.9 s\) and Small\-Feasible\-Region\-1 \(27\.0 s\)\.

Table 15:Inner\-loop solver wall time \(seconds\)\.
### E\.2Full\-Space Solver Comparison

Table[16](https://arxiv.org/html/2608.03045#A5.T16)reports regret when each solver is applied to the full joint problem over\[xWB,xBB\]\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]\. SLSQP\-50 is inapplicable in this context because the black\-box functionfBBf^\{\\mathrm\{BB\}\}does not provide analytical gradients\.

Table 16:Full\-space solver regret\. A dash \(—\) indicates the solver returned an infeasible solution or is inapplicable\. BH \(100 iterations\) fails catastrophically on Rastrigin due to multimodality in the combined space, but BH\-1000 recovers\.The full\-space results exhibit sharper solver differentiation\. SHGO fails on 7 of 13 problems, with 5 hitting the 60\-second timeout on problems with higher dimensionality or tight constraints \(Rastrigin, CSTR, Heat\-Exchanger, PSA, Distillation, Membrane\)\. DA fails on 9 of 13 problems and returns large\-regret solutions even when feasible \(e\.g\.,3\.9×1063\.9\\times 10^\{6\}on PSA\)\. DE is fully reliable \(13/13 feasible\) with uniformly low regret, but is the slowest among the feasible solvers\.

BH achieves feasibility on all 13 problems with near\-zero regret on 12 of them\. The exception is Rastrigin, where BH at 100 iterations converges to a poor local optimum with a regret of9\.9×1099\.9\\times 10^\{9\}\. This failure is resolved by increasing to 1,000 iterations \(BH\-1000\), which achieves zero regret\. This is the only problem where BH requires more than 100 iterations to find the global optimum in the full space\.

### E\.3Solver Selection Rationale

Table[17](https://arxiv.org/html/2608.03045#A5.T17)summarizes the aggregate performance\. BH is the only solver that achieves feasibility on all 13 problems in both contexts while also being applicable to derivative\-free problems\. Although SLSQP\-50 wins on 7 of 13 inner\-loop problems, it cannot be used for the full\-space baseline\. DE is the next most reliable solver \(13/13 feasible in both contexts\) but is slower \(median 0\.97 s vs\. 0\.39 s for BH in the full space\)\. SHGO and DA are unsuitable as general\-purpose solvers due to frequent feasibility failures\.

Table 17:Aggregate solver statistics across all 13 benchmark problems\. Feasibility counts and median regret are computed over feasible solutions only\. BH provides the best combination of universal feasibility, near\-exact optimality, and broad applicability\.The regret gap between BH and the best solver on any given problem is negligible—typically10−1410^\{\-14\}to10−810^\{\-8\}—confirming that the choice of BH does not compromise solution quality\. Because BH wraps an SLSQP local minimizer inside a global exploration loop, it inherits gradient\-based precision on the inner problem while providing the global search capability needed for the full\-space baseline\. This combination of universal feasibility, near\-exact optimality, and applicability to both derivative\-available and derivative\-free contexts makes BH the most defensible single\-solver choice for the benchmark suite\.

## Appendix FGlobal Optimum Verification

Each problem’s global optimumJ⋆J^\{\\star\}was verified through a multi\-stage process designed to provide high confidence that the reported optima are correct\.

### F\.1Verification Procedure

We used two complementary approaches, each repeated across 10 independent random seeds\. First, we ran Basin\-Hopping \(BH\) over the full variable space\[xWB,xBB\]\[x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\]with an SLSQP local minimizer, a Metropolis temperature ofT=100T=100, full\-range stepsize, and 2,000 BH iterations\. Constraints were enforced by the SLSQP local minimizer\. Second, we ran multi\-start SLSQP with 500 Latin hypercube starting points per seed\.

For each problem, we compared the best feasible solutions across all 20 runs \(10 BH \+ 10 multi\-start NLP\)\. All 13 problems achieved agreement to within2×10−42\\times 10^\{\-4\}relative error, confirming the reported optima\. For constrained problems, we verified that the optimal solution satisfies all inequality constraintsgi​\(xWB⁣⋆,y⋆\)≤0g\_\{i\}\(x^\{\\mathrm\{WB\}\\star\},y^\{\\star\}\)\\leq 0to within a tolerance of10−610^\{\-6\}\.

### F\.2Verified Optima

Table[18](https://arxiv.org/html/2608.03045#A6.T18)reports the verified global optimum for each problem, along with the optimal decision variables\.

Table 18:Verified global optima for all 13 benchmark problems\. Optima verified via Basin\-Hopping \(SLSQP local minimizer,T=100T=100, 2000 iterations, 10 seeds\) and multi\-start SLSQP \(500 Latin hypercube starts, 10 seeds\)\.

## Appendix GProof of Proposition[1](https://arxiv.org/html/2608.03045#Thmtheorem1)

###### Proof\.

We show thatinf\(bilevel\)=inf\(original\)\\inf\(\\text\{bilevel\}\)=\\inf\(\\text\{original\}\)by constructing mappings in both directions\.

*\(i\) Any feasible point of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) maps to a feasible outer\-level point with equal or worse objective\.*Let\(xWB,xBB\)\(x^\{\\mathrm\{WB\}\},x^\{\\mathrm\{BB\}\}\)be feasible for \([4](https://arxiv.org/html/2608.03045#S2.E4)\) withy=fBB​\(xBB\)y=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\}\)\. By Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\), the inner problem \([7c](https://arxiv.org/html/2608.03045#S4.E7.3)\)–\([7f](https://arxiv.org/html/2608.03045#S4.E7.6)\) at thisxBBx^\{\\mathrm\{BB\}\}andyyhas a global solutionxWB⁣⋆x^\{\\mathrm\{WB\}\\star\}satisfyingJ​\(xWB⁣⋆,y,xBB\)≤J​\(xWB,y,xBB\)J\(x^\{\\mathrm\{WB\}\\star\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\\leq J\(x^\{\\mathrm\{WB\}\},\\allowbreak y,\\allowbreak x^\{\\mathrm\{BB\}\}\)\. ThusxBBx^\{\\mathrm\{BB\}\}is outer\-feasible with objective value no worse than the original point\.

*\(ii\) Any outer\-level feasible point maps back to a feasible point of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) with equal objective\.*Letx~BB\\tilde\{x\}^\{\\mathrm\{BB\}\}be outer\-feasible withy~=fBB​\(x~BB\)\\tilde\{y\}=f^\{\\mathrm\{BB\}\}\(\\tilde\{x\}^\{\\mathrm\{BB\}\}\)and inner solutionx~WB⁣⋆\\tilde\{x\}^\{\\mathrm\{WB\}\\star\}\. By construction,x~WB⁣⋆\\tilde\{x\}^\{\\mathrm\{WB\}\\star\}satisfiesfWB​\(x~WB⁣⋆,y~\)=0f^\{\\mathrm\{WB\}\}\(\\tilde\{x\}^\{\\mathrm\{WB\}\\star\},\\allowbreak\\tilde\{y\}\)=0,g​\(x~WB⁣⋆,y~,x~BB\)≤0g\(\\tilde\{x\}^\{\\mathrm\{WB\}\\star\},\\allowbreak\\tilde\{y\},\\allowbreak\\tilde\{x\}^\{\\mathrm\{BB\}\}\)\\leq 0, andx~WB⁣⋆∈𝒳WB\\tilde\{x\}^\{\\mathrm\{WB\}\\star\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\}\. Together withx~BB∈𝒳BB\\tilde\{x\}^\{\\mathrm\{BB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{BB\}\}, the pair\(x~WB⁣⋆,x~BB\)\(\\tilde\{x\}^\{\\mathrm\{WB\}\\star\},\\tilde\{x\}^\{\\mathrm\{BB\}\}\)is feasible for \([4](https://arxiv.org/html/2608.03045#S2.E4)\) with objectiveJ​\(x~WB⁣⋆,y~,x~BB\)J\(\\tilde\{x\}^\{\\mathrm\{WB\}\\star\},\\tilde\{y\},\\tilde\{x\}^\{\\mathrm\{BB\}\}\)\.

*\(iii\) Infeasibility\.*If the inner problem is infeasible for somey~∈range​\(fBB\)\\tilde\{y\}\\in\\text\{range\}\(f^\{\\mathrm\{BB\}\}\)—i\.e\., noxWB∈𝒳WBx^\{\\mathrm\{WB\}\}\\in\\mathcal\{X\}^\{\\mathrm\{WB\}\}satisfiesfWB​\(xWB,y~\)=0f^\{\\mathrm\{WB\}\}\(x^\{\\mathrm\{WB\}\},\\allowbreak\\tilde\{y\}\)=0andg​\(xWB,y~,x~BB\)≤0g\(x^\{\\mathrm\{WB\}\},\\allowbreak\\tilde\{y\},\\allowbreak\\tilde\{x\}^\{\\mathrm\{BB\}\}\)\\leq 0—then no feasible point of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) exists at thisx~BB\\tilde\{x\}^\{\\mathrm\{BB\}\}either, so excluding it from the outer search does not lose optimality\.

From \(i\) and \(ii\), the feasible sets project onto each other with preserved objectives, soinf\(bilevel\)=inf\(original\)\\inf\(\\text\{bilevel\}\)=\\inf\(\\text\{original\}\)\. To see how a global optimum of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) appears in \([7](https://arxiv.org/html/2608.03045#S4.E7)\): let\(xWB⁣⋆⋆,xBB⁣⋆⋆\)\(x^\{\\mathrm\{WB\}\\star\\star\},x^\{\\mathrm\{BB\}\\star\\star\}\)be a global minimizer of \([4](https://arxiv.org/html/2608.03045#S2.E4)\) withy⋆⋆=fBB​\(xBB⁣⋆⋆\)y^\{\\star\\star\}=f^\{\\mathrm\{BB\}\}\(x^\{\\mathrm\{BB\}\\star\\star\}\)\. By step \(i\),xBB⁣⋆⋆x^\{\\mathrm\{BB\}\\star\\star\}is outer\-feasible in \([7](https://arxiv.org/html/2608.03045#S4.E7)\), and the inner problem atxBB⁣⋆⋆x^\{\\mathrm\{BB\}\\star\\star\}returns somex^WB\\hat\{x\}^\{\\mathrm\{WB\}\}withJ​\(x^WB,y⋆⋆,xBB⁣⋆⋆\)≤J​\(xWB⁣⋆⋆,y⋆⋆,xBB⁣⋆⋆\)J\(\\hat\{x\}^\{\\mathrm\{WB\}\},\\allowbreak y^\{\\star\\star\},\\allowbreak x^\{\\mathrm\{BB\}\\star\\star\}\)\\leq J\(x^\{\\mathrm\{WB\}\\star\\star\},\\allowbreak y^\{\\star\\star\},\\allowbreak x^\{\\mathrm\{BB\}\\star\\star\}\)\. Since\(xWB⁣⋆⋆,xBB⁣⋆⋆\)\(x^\{\\mathrm\{WB\}\\star\\star\},x^\{\\mathrm\{BB\}\\star\\star\}\)is already globally optimal for \([4](https://arxiv.org/html/2608.03045#S2.E4)\), the inequality must hold with equality:x^WB=xWB⁣⋆⋆\\hat\{x\}^\{\\mathrm\{WB\}\}=x^\{\\mathrm\{WB\}\\star\\star\}\(up to alternative optima with the same objective value\)\. Thus the bilevel formulation recovers both the optimalxBB⁣⋆⋆x^\{\\mathrm\{BB\}\\star\\star\}at the outer level and the optimalxWB⁣⋆⋆x^\{\\mathrm\{WB\}\\star\\star\}at the inner level\.

Assumption[1](https://arxiv.org/html/2608.03045#Thmassumption1)\(ii\) is used in step \(i\) to ensure the inner global optimum is attained, and in step \(ii\) to ensure the inner solution is globally optimal \(not merely a local minimum\)\. The surrogate in Algorithm[1](https://arxiv.org/html/2608.03045#alg1)models the mapping \([5](https://arxiv.org/html/2608.03045#S3.E5)\), which is a function ofnBBn\_\{\\mathrm\{BB\}\}variables\. ∎

## Appendix HImplementation Details

*Infeasibility handling\.*If the inner problem is infeasible for someyy—because the black\-box outputs lead to a white\-box problem with no feasible solution—we return a large penalty value\. The acquisition function naturally steers away from these regions without requiring explicit constraint modeling at the outer level\.

*Modularity\.*The inner optimizer can be swapped \(Basin\-Hopping, SLSQP, IPOPT, BARON\) without changing the outer BO framework\. Similarly, any acquisition function \(EI, PI, LCB\) can be used in the outer loop\.

## Appendix IConvergence Speed Details

Table 19:Median iterations to reach 1% of initial regret for Bi\-BO \(SLSQP\) \(ninit=50n\_\{\\text\{init\}\}=50, bestξ\\xi, 10 repetitions\)\. “\>\>200” means the target was not reached within 200 iterations\. Speedup values marked with≥\\geqare lower bounds \(BB\-BO did not converge\)\.Table 20:Post\-initialization BO iterations required for Bi\-BO \(SLSQP\) to reach BB\-BO’s median final regret, and vice versa \(ninit=50n\_\{\\text\{init\}\}=50, 200\-iteration BO budget, 10 repetitions, bestξ\\xiper problem per method\)\. An iteration count of0means the target is reached at the end of the initial sampling phase\. “\>\>200” indicates the target is never reached within the budget in the majority of repetitions\. Bi\-BO’s NLP\-assisted initial samples already dominate BB\-BO’s full\-budget final on 11 of 13 problems\.
## Appendix JHyperparameter Sensitivity Details

We sweep over two key BO hyperparameters: the number of initial samplesninit∈\{5,25,50\}n\_\{\\text\{init\}\}\\in\\\{5,25,50\\\}drawn from a Latin hypercube design, and the EI exploration parameterξ∈\{0\.001,0\.01,0\.05,0\.1,0\.2,0\.5,1\.0\}\\xi\\in\\\{0\.001,0\.01,0\.05,0\.1,0\.2,0\.5,1\.0\\\}\. Each configuration is repeated 10 times with independent seeds 42–51 to quantify variability\. Every BO run uses 200 iterations after initialization\. In total, the experiment comprises 8,450 independent optimization runs:13×3×7×10×3=8,19013\\times 3\\times 7\\times 10\\times 3=8\{,\}190BO runs \(one black\-box and two bilevel variants\) plus13×10=13013\\times 10=130NLP runs and13×10=13013\\times 10=130full\-space BH runs \(one each per problem\-seed combination, since NLP and BH results do not depend onninitn\_\{\\text\{init\}\}orξ\\xi\)\.

Table 21:Bi\-BO \(SLSQP\) mean final regret acrossninitn\_\{\\text\{init\}\}values \(bestξ\\xiper problem perninitn\_\{\\text\{init\}\}, 10 repetitions\)\.Table 22:Best EI exploration parameterξ\\xiper problem for Bi\-BO \(SLSQP\), selected by lowest mean final regret across allninitn\_\{\\text\{init\}\}values\. These are theξ\\xivalues used to report results in Table[3](https://arxiv.org/html/2608.03045#S6.T3)\.![Refer to caption](https://arxiv.org/html/2608.03045v1/x8.png)Figure 9:Sensitivity to the EI exploration parameterξ\\xi\. Left: bilevel BO; right: black\-box BO\. Rows are problems, columns areξ\\xivalues, aggregated over allninitn\_\{\\text\{init\}\}\. Cell values show log10\(median final regret\); bold indicates bestξ\\xiper problem\. No singleξ\\xiis universally optimal, but bilevel BO outperforms BB\-BO regardless ofξ\\xi\.
## References

- \[1\]Peter I\. Frazier\.A tutorial on Bayesian optimization\.arXiv preprint arXiv:1807\.02811, 2018\.
- \[2\]Fani Boukouvala, M\. M\. Faruque Hasan, and Christodoulos A\. Floudas\.Global optimization of general constrained grey\-box models: New method and its application to constrained PDEs for pressure swing adsorption\.Journal of Global Optimization, 67\(1\):3–42, 2017\.
- \[3\]Fani Boukouvala and Christodoulos A\. Floudas\.ARGONAUT: AlgoRithms for global optimization of constrained grey\-box computational problems\.Optimization Letters, 11:895–913, 2017\.
- \[4\]John P\. Eason and Lorenz T\. Biegler\.A trust region filter method for glass box/black box optimization\.AIChE Journal, 62\(9\):3124–3136, 2016\.
- \[5\]Raul Astudillo and Peter I\. Frazier\.Bayesian optimization of composite functions\.InProceedings of the 36th International Conference on Machine Learning \(ICML\), pages 354–363\. PMLR, 2019\.
- \[6\]Joel A\. Paulson and Congwen Lu\.COBALT: COnstrained Bayesian optimizAtion of computationaLly expensive grey\-box models exploiting derivative information\.Computers & Chemical Engineering, 160:107700, 2022\.
- \[7\]Akshay Kudva and Joel A Paulson\.Bonsai: Structure\-exploiting robust bayesian optimization for networked black\-box systems under uncertainty\.Computers & Chemical Engineering, page 109393, 2025\.
- \[8\]Burcu Beykal, Styliani Avraamidou, Ioannis P\. E\. Pistikopoulos, Melis Onel, and Efstratios N\. Pistikopoulos\.DOMINO: Data\-driven optimization of bi\-level mixed\-integer nonlinear problems\.Journal of Global Optimization, 78:1–36, 2020\.
- \[9\]Emmanuel Kieffer, Grégoire Danoy, Pascal Bouvry, and Anass Nagih\.Bayesian optimization approach of general bi\-level problems\.InProceedings of the Genetic and Evolutionary Computation Conference Companion \(GECCO\), pages 1614–1621\. ACM, 2017\.
- \[10\]Donald R\. Jones, Matthias Schonlau, and William J\. Welch\.Efficient global optimization of expensive black\-box functions\.Journal of Global Optimization, 13\(4\):455–492, 1998\.
- \[11\]Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P\. Adams, and Nando de Freitas\.Taking the human out of the loop: A review of Bayesian optimization\.Proceedings of the IEEE, 104\(1\):148–175, 2016\.
- \[12\]Chris A\. Kieslich, Fani Boukouvala, and Christodoulos A\. Floudas\.Optimization of black\-box problems using Smolyak grids and polynomial approximations\.Journal of Global Optimization, 71\(4\):845–869, 2018\.
- \[13\]John P\. Eason and Lorenz T\. Biegler\.Advanced trust region optimization strategies for glass box/black box models\.AIChE Journal, 64\(11\):3934–3943, 2018\.
- \[14\]Raul Astudillo and Peter I\. Frazier\.Bayesian optimization of function networks\.InAdvances in Neural Information Processing Systems, volume 34, pages 14463–14475, 2021\.
- \[15\]Poompol Buathong, Jiayue Wan, Raul Astudillo, Sam Daulton, Maximilian Balandat, and Peter I\. Frazier\.Bayesian optimization of function networks with partial evaluations\.InProceedings of the 41st International Conference on Machine Learning \(ICML\), pages 4752–4784, 2024\.
- \[16\]Leonardo D\. González and Victor M\. Zavala\.BOIS: Bayesian Optimization of Interconnected Systems\.IFAC\-PapersOnLine, 58\(14\):446–451, 2024\.
- \[17\]Joschka Winz, Florian Fromme, and Sebastian Engell\.Bayesian optimization of gray\-box process models using a modified upper confidence bound acquisition function\.Computers & Chemical Engineering, 194:108976, 2025\.
- \[18\]Zeynep H\. Gümüş and Christodoulos A\. Floudas\.Global optimization of nonlinear bilevel programming problems\.Journal of Global Optimization, 20\(1\):1–31, 2001\.
- \[19\]Alexander Mitsos, Benoît Chachuat, and Paul I\. Barton\.Towards global bilevel dynamic optimization\.Journal of Global Optimization, 45:63–93, 2009\.
- \[20\]Polyxeni\-Margarita Kleniati and Claire S\. Adjiman\.Branch\-and\-sandwich: a partially relaxed branch\-and\-bound algorithm for bi\-level problems\. Part I: General algorithm\.Journal of Global Optimization, 60\(3\):425–458, 2014\.
- \[21\]Nuno P\. Faísca, Vivek Dua, Berç Rustem, Pedro M\. Saraiva, and Efstratios N\. Pistikopoulos\.Parametric global optimisation for bilevel programming\.Journal of Global Optimization, 38\(4\):609–623, 2007\.
- \[22\]Omer Ekmekcioglu, Nursen Aydin, and Juergen Branke\.Bayesian optimization of bilevel problems\.arXiv preprint arXiv:2412\.18518, 2024\.
- \[23\]Ruth W\. T\. Chew, Quoc Phong Nguyen, and Bryan Kian Hsiang Low\.BILBO: BIlevel Bayesian Optimization\.InProceedings of the 42nd International Conference on Machine Learning \(ICML\), 2025\.arXiv:2502\.02121\.
- \[24\]Shi Fu, Fengxiang He, Xinmei Tian, and Dacheng Tao\.Convergence of Bayesian bilevel optimization\.InThe Twelfth International Conference on Learning Representations \(ICLR\), 2024\.
- \[25\]Michael Baldea\.A multiscale Bayesian optimization framework for process and material codesign\.AIChE Journal, 2026\.
- \[26\]Jacob Gardner, Matt Kusner, Zhixiang Xu, Kilian Weinberger, and John Cunningham\.Bayesian optimization with inequality constraints\.In Eric P\. Xing and Tony Jebara, editors,Proceedings of the 31st International Conference on Machine Learning, volume 32 ofProceedings of Machine Learning Research, pages 937–945, Bejing, China, 22–24 Jun 2014\. PMLR\.
- \[27\]Michael A\. Gelbart, Jasper Snoek, and Ryan P\. Adams\.Bayesian optimization with unknown constraints\.InProceedings of the 30th Conference on Uncertainty in Artificial Intelligence \(UAI\), pages 250–259, 2014\.
- \[28\]Victor Picheny, Robert B\. Gramacy, Stefan Wild, and Sébastien Le Digabel\.Bayesian optimization under mixed constraints with a slack\-variable augmented Lagrangian\.InAdvances in Neural Information Processing Systems, volume 29, pages 1435–1443, 2016\.
- \[29\]DJ Wales\.Energy landscapes \(cambridge molecular science\), 2003\.
- \[30\]David J Wales and Jonathan PK Doye\.Global optimization by basin\-hopping and the lowest energy structures of lennard\-jones clusters containing up to 110 atoms\.The Journal of Physical Chemistry A, 101\(28\):5111–5116, 1997\.
- \[31\]Zhenqin Li and Harold A Scheraga\.Monte carlo\-minimization approach to the multiple\-minima problem in protein folding\.Proceedings of the National Academy of Sciences, 84\(19\):6611–6615, 1987\.
- \[32\]David J Wales and Harold A Scheraga\.Global optimization of clusters, crystals, and biomolecules\.Science, 285\(5432\):1368–1372, 1999\.
- \[33\]Brian Olson, Irina Hashmi, Kevin Molloy, and Amarda Shehu\.Basin hopping as a general and versatile optimization framework for the characterization of biological macromolecules\.Advances in Artificial Intelligence, 2012\(1\):674832, 2012\.
- \[34\]Jorge Nocedal and Stephen J Wright\.Numerical optimization\.Springer, 2006\.
- \[35\]Dieter Kraft\.A software package for sequential quadratic programming\.Forschungsbericht\- Deutsche Forschungs\- und Versuchsanstalt fur Luft\- und Raumfahrt, 1988\.
- \[36\]Charles L Lawson and Richard J Hanson\.Solving least squares problems\.SIAM, 1995\.
- \[37\]P\. N\. Suganthan, N\. Hansen, J\. J\. Liang, K\. Deb, Y\.\-P\. Chen, A\. Auger, and S\. Tiwari\.Problem definitions and evaluation criteria for the CEC 2005 special session on real\-parameter optimization\.Technical report, Nanyang Technological University, 2005\.
- \[38\]Nikolaus Hansen, Anne Auger, Raymond Ros, Olaf Mersmann, Tea Tušar, and Dimo Brockhoff\.COCO: A platform for comparing continuous optimizers in a black\-box setting\.Optimization Methods and Software, 36\(1\):114–144, 2021\.
- \[39\]Setareh Ariafar, Jaume Coll\-Font, Dana Brooks, and Jennifer Dy\.Admmbo: Bayesian optimization with unknown constraints using admm\.Journal of Machine Learning Research, 20\(123\):1–26, 2019\.
- \[40\]Stephen J\. Wright\.Coordinate descent algorithms\.Mathematical Programming, 151\(1\):3–34, 2015\.
- \[41\]Benoît Colson, Patrice Marcotte, and Gilles Savard\.An overview of bilevel optimization\.Annals of Operations Research, 153\(1\):235–256, 2007\.
- \[42\]R\. B\. Gramacy, G\. A\. Gray, S\. Le Digabel, H\. K\. H\. Lee, P\. Ranjan, G\. Wells, and S\. M\. Wild\.Modeling an augmented lagrangian for blackbox constrained optimization\.Technometrics, 58\(1\):1–11, 2016\.
- \[43\]L\. A\. Rastrigin\.Systems of extremal control\.Nauka, 1974\.
- \[44\]H\. Rosenbrock\.An automatic method for finding the greatest or least value of a function\.The Computer Journal, 3\(3\):175–184, 1960\.
- \[45\]Timothy F\. Yee and Ignacio E\. Grossmann\.Simultaneous optimization models for heat integration—II\. Heat exchanger network synthesis\.Computers & Chemical Engineering, 14\(10\):1165–1184, 1990\.
- \[46\]Max S\. Peters and Klaus D\. Timmerhaus\.Plant Design and Economics for Chemical Engineers\.McGraw\-Hill, 4th edition, 1991\.
- \[47\]Fabiana T\. Mizutani, Fernando L\. P\. Pessoa, Eduardo M\. Queiroz, Steinar Hauan, and Ignacio E\. Grossmann\.Mathematical programming model for heat\-exchanger network synthesis including detailed heat\-exchanger designs\.Industrial & Engineering Chemistry Research, 42\(17\):4009–4018, 2003\.
- \[48\]Mauro A\. S\. S\. Ravagnani and José A\. Caballero\.A MINLP model for the rigorous design of shell and tube heat exchangers using the TEMA standards\.Chemical Engineering Research and Design, 85\(10\):1423–1435, 2007\.
- \[49\]Anshul Agarwal, Lorenz T\. Biegler, and Stephen E\. Zitney\.Superstructure\-based optimal synthesis of pressure swing adsorption cycles for precombustion CO2 capture\.Industrial & Engineering Chemistry Research, 49\(11\):5066–5079, 2010\.
- \[50\]Douglas M\. Ruthven, Shamsuzzaman Farooq, and Kent S\. Knaebel\.Pressure Swing Adsorption\.VCH Publishers, New York, 1994\.
- \[51\]Tyler D\. Burns, Kasturi Nagesh Pai, Sai Gokul Subraveti, Sean P\. Collins, Mykhaylo Krykunov, Arvind Rajendran, and Tom K\. Woo\.Process\-level modelling and optimization to evaluate metal\-organic frameworks for post\-combustion capture\.Molecular Systems Design & Engineering, 5:1456–1469, 2020\.
- \[52\]Berend Smit and Theo L\. M\. Maesen\.Molecular simulations of zeolites: Adsorption, diffusion, and shape selectivity\.Chemical Reviews, 108\(10\):4125–4184, 2008\.
- \[53\]David Dubbeldam and Randall Q\. Snurr\.Design, parameterization, and implementation of atomic force fields for adsorption in nanoporous materials\.Advanced Theory and Simulations, 2\(11\):1900135, 2019\.
- \[54\]Irving Langmuir\.The adsorption of gases on plane surfaces of glass, mica and platinum\.Journal of the American Chemical Society, 40\(9\):1361–1403, 1918\.
- \[55\]Yiwen Wang, Jie Chen, Xiang Li, and Wei Zhang\.Adsorption process optimization and adsorbent evaluation based on Langmuir isotherm model\.Langmuir, 39\(46\):16404–16414, 2023\.
- \[56\]Ralph T\. Yang\.Gas Separation by Adsorption Processes\.Butterworths, Boston, 1987\.
- \[57\]M\. M\. Faruque Hasan, Eric L\. First, and Christodoulos A\. Floudas\.Cost\-effective CO2 capture based on in silico screening of zeolites and process optimization\.Physical Chemistry Chemical Physics, 15\(40\):17601–17618, 2013\.
- \[58\]Kaitlin T\. Leperi, Randall Q\. Snurr, and Fengqi You\.Surrogate models based on artificial neural networks to simulate and optimize pressure swing adsorption cycles for CO2 capture\.Industrial & Engineering Chemistry Research, 58\(39\):18241–18252, 2019\.
- \[59\]J\. M\. Smith, Hendrick C\. Van Ness, Michael M\. Abbott, and Mark T\. Swihart\.Introduction to Chemical Engineering Thermodynamics\.McGraw\-Hill Education, 8th edition, 2018\.
- \[60\]Satish J\. Parulekar\.Yield optimization for multiple reactions\.Chemical Engineering Science, 43\(8\):2131–2139, 1988\.
- \[61\]Rein Luus, Jens Dittrich, and Frerich J\. Keil\.Towards practical optimal control of batch reactors\.Chemical Engineering Science, 54\(17\):4137–4145, 1999\.
- \[62\]H\. Scott Fogler\.Elements of Chemical Reaction Engineering\.Prentice Hall, 5th edition, 2016\.
- \[63\]Ali Hussain Motagamwala and James A\. Dumesic\.Microkinetic modeling: A tool for rational catalyst design\.Chemical Reviews, 121\(2\):1049–1076, 2021\.
- \[64\]Hideshi Ooka, Jun Huang, and Kai S\. Exner\.The Sabatier principle in electrocatalysis: Basics, limitations, and extensions\.Frontiers in Energy Research, 9:654460, 2021\.
- \[65\]Matthew D\. Wodrich, Benjamin Sawatlon, Michael Busch, and Clémence Corminboeuf\.Microkinetic molecular volcano plots for enhanced catalyst selectivity and activity predictions\.ACS Catalysis, 14:3523–3532, 2024\.
- \[66\]Stefan Kossack, Korbinian Kraemer, Rafiqul Gani, and Wolfgang Marquardt\.A systematic synthesis framework for extractive distillation processes\.Chemical Engineering Research and Design, 86\(7\):781–792, 2008\.
- \[67\]Jan Scheffczyk, Lorenz Fleitmann, André Schwarz, Maximilian Hoppe, Kai Leonhard, and André Bardow\.COSMO\-CAMPD: A framework for integrated design of molecules and processes based on COSMO\-RS\.Molecular Systems Design & Engineering, 3\(4\):645–657, 2018\.
- \[68\]Merrell R\. Fenske\.Fractionation of straight\-run Pennsylvania gasoline\.Industrial & Engineering Chemistry, 24\(5\):482–485, 1932\.
- \[69\]Arthur J\. V\. Underwood\.Fractional distillation of multicomponent mixtures\.Chemical Engineering Progress, 44\(8\):603–614, 1948\.
- \[70\]Ramkumar S\. Kamath, Lorenz T\. Biegler, and Ignacio E\. Grossmann\.An equation\-oriented approach for handling thermodynamics based on group methods in process optimization\.Computers & Chemical Engineering, 34\(12\):2085–2096, 2010\.
- \[71\]Edwin R\. Gilliland\.Multicomponent rectification: Estimation of the number of theoretical plates as a function of the reflux ratio\.Industrial & Engineering Chemistry, 32\(9\):1220–1223, 1940\.
- \[72\]Ignacio E\. Grossmann, Pio A\. Aguirre, and Matías Barttfeld\.Optimal synthesis of complex distillation columns using rigorous models\.Computers & Chemical Engineering, 29\(6\):1203–1215, 2005\.
- \[73\]Stefan Bruggemann and Wolfgang Marquardt\.Shortcut methods for nonideal multicomponent distillation: 3\. Extractive distillation columns\.AIChE Journal, 50\(6\):1129–1149, 2004\.
- \[74\]James B\. Hillenbrand and Arthur W\. Westerberg\.The synthesis of multiple\-effect evaporator systems using minimum utility insights\.Computers & Chemical Engineering, 12\(6\):611–624, 1988\.
- \[75\]Swenson Technology\.Multiple effect evaporators\.[https://swensontechnology\.com/multiple\-effect\-evaporators/](https://swensontechnology.com/multiple-effect-evaporators/), 2024\.Accessed: 2024\.
- \[76\]Elvis Ahmetović, Zdravko Kravanja, and Ignacio E\. Grossmann\.Simultaneous optimisation and heat integration of evaporation systems including mechanical vapour recompression and background process\.Energy, 158:1160–1191, 2018\.
- \[77\]Atul Sharma, V\. V\. Tyagi, C\. R\. Chen, and D\. Buddhi\.Review on thermal energy storage with phase change materials and applications\.Renewable and Sustainable Energy Reviews, 13\(2\):318–345, 2009\.
- \[78\]Belén Zalba, José M\. Marín, Luisa F\. Cabeza, and Harald Mehling\.Review on thermal energy storage with phase change: Materials, heat transfer analysis and applications\.Applied Thermal Engineering, 23\(3\):251–283, 2003\.
- \[79\]Rongxin Qi and Michael A\. Henson\.Membrane system design for multicomponent gas mixtures via MINLP optimization\.Computers & Chemical Engineering, 24\(12\):2719–2737, 2000\.
- \[80\]Lloyd M\. Robeson\.Correlation of separation factor versus permeability for polymeric membranes\.Journal of Membrane Science, 62\(2\):165–185, 1991\.
- \[81\]Lloyd M\. Robeson\.The upper bound revisited\.Journal of Membrane Science, 320\(1–2\):390–400, 2008\.
- \[82\]Benny D\. Freeman\.Basis of permeability/selectivity tradeoff relations in polymeric gas separation membranes\.Macromolecules, 32\(2\):375–380, 1999\.
- \[83\]Johannes G\. Wijmans and Richard W\. Baker\.The solution\-diffusion model: A review\.Journal of Membrane Science, 107\(1–2\):1–21, 1995\.
- \[84\]Jason W\. Barnett, Connor R\. Bilchak, Yiwen Wang, Brian C\. Benicewicz, Laura A\. Murdock, Taner Berber, and Sanat K\. Kumar\.Designing exceptional gas\-separation polymer membranes using machine learning\.Science Advances, 6\(20\):eaaz4301, 2020\.
- \[85\]Ho Bum Park, Jovan Kamcev, Lloyd M\. Robeson, Menachem Elimelech, and Benny D\. Freeman\.Maximizing the right stuff: The trade\-off between membrane permeability and selectivity\.Science, 356\(6343\):eaab0530, 2017\.
- \[86\]Morrel H\. Cohen and David Turnbull\.Molecular transport in liquids and glasses\.The Journal of Chemical Physics, 31\(5\):1164–1169, 1959\.
- \[87\]Hiroshi Fujita\.Diffusion in polymer\-diluent systems\.Fortschritte der Hochpolymeren\-Forschung, 3:1–47, 1961\.
- \[88\]Scott Matteucci, Yuri Yampolskii, Benny D\. Freeman, and Ingo Pinnau\.Transport of gases and vapors in glassy and rubbery polymers\.Materials Science of Membranes for Gas and Vapor Separation, pages 1–47, 2006\.
- \[89\]W\. M\. Lee\.Selection of barrier materials from molecular structure\.Polymer Engineering & Science, 20\(1\):65–69, 1980\.
- \[90\]A\. Bondi\.van der Waals volumes and radii\.The Journal of Physical Chemistry, 68\(3\):441–451, 1964\.
- \[91\]Natsayi Chiwaye, Thokozani Majozi, and Michael O\. Daramola\.Superstructure\-based optimization of membrane gas separation processes: A review\.Industrial & Engineering Chemistry Research, 64:13905–13919, 2025\.
- \[92\]T\. J\. Williams and R\. E\. Otto\.A generalized chemical processing model for the investigation of computer control\.Transactions of the American Institute of Electrical Engineers, Part I: Communication and Electronics, 79\(5\):458–473, 1960\.
- \[93\]Jochen Schmid, Katrin Teichert, Moncef Chioua, Thorsten Schindler, and Michael Bortz\.Simulation and optimal control of the Williams\-Otto process using Pyomo\.Chemie Ingenieur Technik, 92\(11\):1728–1740, 2020\.
- \[94\]Rainer Storn and Kenneth Price\.Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces\.Journal of Global Optimization, 11\(4\):341–359, 1997\.
- \[95\]Stefan C\. Endres, Carl Sandrock, and Walter W\. Focke\.A simplicial homology algorithm for Lipschitz optimisation\.Journal of Global Optimization, 72\(2\):181–217, 2018\.
- \[96\]Y\. Xiang, D\. Y\. Sun, W\. Fan, and X\. G\. Gong\.Generalized simulated annealing algorithm and its application to the Thomson model\.Physics Letters A, 233\(3\):216–220, 1997\.