Discovering Model Structure of Dynamical Systems with Combinatorial Bayesian Optimization
Authors: Lucas Rath
Alexander von Rohr
Andreas Schultze
Sebastian Trimpe
Burkhard Corves
Venue: Accepted by TMLR ()

Abstract: Deciding on a model structure is a fundamental problem in machine learning. In this paper we consider the problem of building a data-based model for dynamical systems from a library of discrete components. In addition to optimizing performance, we consider crash and inequality constraints that arise from additional requirements, such as real-time capability and model complexity. We address this task of model structure selection with a focus on dynamical systems and propose to search over potential model structures efficiently using a constrained combinatorial Bayesian Optimization (BO) algorithm. We propose expressive surrogate models suited for combinatorial domains and an acquisition function that can handle inequality and crash constraints. We provide simulated benchmark problems within the domain of equation discovery of nonlinear dynamical systems. Our method outperforms the state-of-the-art in constrained combinatorial optimization of black-box functions and has a favorable computational overhead compared to other BO methods. As a real-world application example, we apply our method to optimize the configuration of an electric vehicle's digital twin while ensuring its real-time capability for the use in one of the world's largest driving simulators.

Published in Transactions on Machine Learning Research (03/2024)
Discovering Model Structure of Dynamical Systems with Combinatorial Bayesian Optimization
Lucas Rath
lucas. moos.roth@bbmn
Institute for Mechanism Theory, Machine Dynamics and Robotics, RWTH Aachen University
BMW AG
Alexander von Rohr
Instibute for Data Science in Mechanicul Engineering, RWTH Aachen University
Andreas Schultze
andreas.schullze@Bmme.
BMW AG
Scbastian Trimpe
Instibute for Data Science in Mechanical Engineering, RWTH Aachen University
Burkhard Corves
corves Coruesskgmn..uutb-accce.
Institute for Mechanism Theory, Mochine Dynamics and Robolics, RWTH Aachen University
Reviewed on OpenReview: https: // openrevieu net/ net/forum? forum? id= 2400vQmJBKK
Abstract
Deciding model structure is fundamental problem in machine learning. In this paper we consider the problem of building data-based model for dynamical systems from library of discrete components. addition to optimizing performance, we consider crash and inequality constraints that arise from additional requirements, such as real-time capability and model complexity. We address this task of model structure selection with a focus on dynamical systems and propose to search over potential model structures efficiently using constrained combinatorial Bayesian Optimization (BO) algorithm. We propose expressive surrogate models suited for combinatorial domains and acquisition function that can handle inequality and crash constraints. We provide simulated benchmark problems within the domain of equation discovery of nonlinear dynamical systems. Our method outperforms the statc-of-the-art in constrained combinatorial optimization of black-box functions and has a favorable computational overhead compared to other BO methods. a real-world application example, we apply our method optimize the configuration of an electric vehicle's digital twin while ensuring its real-time capability for the in one of the world's largest driving simulators.
1 Introduction
Minimizing the discrepancy between the simulated behavior of digital twins and the mcasured behavior of the observed system is essential in engineering applications and is commonly known system identification (Ljung, 1998) model learning (Nguyen- Tuong & Peters, 2011). In general, system identification comprises two interleaved subtasks, structure identification and paramcter estimation. The former aims to determine the structure of model equations, while the latter focuses on finding the model parameters (Tanevski et al., 2015). While a lot of the machine learning literature focuses on the latter task, we tackle the former.
Structure identification methods explore the space of potential models and select the best-fitting representation of the dynamical system. On the one hand, structure identification can be approached by symbolic-regressio methods (Bongard & Lipson, 2007; Schmidt & Lipson, 2009; Brunton et al., 2016), which explore the space

---

Published in Transactions on Machine Learning Research (03/2024)
Model Library
Measurenennts
Simulation & Performance
Model Struciure
Parameter Estimation
Selection
Evaluation
Figure 1: The process of system identification
of possible arithmetical expressions for building mathematical models and often regularization techniques to introduce bias towards simpler models (Tanevski et al., 2020). These methods can be very general and require only the definition of the mathematical operations and functions that can compose the description of the system's 's behavior. On the other hand, knowledge-driver approaches (Bradley et al., 2001) require experts to encode domain-specific knowledge into model fragments, usually following known relationships, physical laws, constraints, or structural properties of the obscrved system. In both approaches, the idea is to break down the overall model into smaller, more manageable fragments or sub-models. This decomposition allows for modular representation of the dynamics and enables the identification of local relationships, which is important when dealing with complex systems. In addition, both methods provide flexibility in defining the structure and form of the model fragments, which can be combined, composed, or connected in way that reflects the relationships and interactions between the different components of the system.
In Figure 1, we depict an iterative process of system identification. In the structure identification phase, model components are selected from a model library of predefined parameterized templates that align with the system's 's characteristics. In the next step, regression algorithms parameters associated with the selected templates of the model fragments. The final model is then evaluated using criteria such as goodness-offfi fit, model complexity, or predictive performance. This evaluation is used as leedback to the search algorithms to select better model structures for the next iteration.
employed to estimate the coefficients
In practice, structure identification search algorithms enumerate plausible model structures and select the model templates that result in the best model evaluation (Tanevski al., 2020). Due to the discrete nature of the problem, this task often approached as a combinatorial optimization problem. The interplay between components disallows selecting templates independently, and exhaustive search over all potential model structures quickly becomes infeasible due to the combinatorial explosion of the search domain. In addition, evaluating a model structure candidate involves fitting coefficients and simulating dynamical systems, which can be computationally demanding. This motivates the need for sample-efficient Bayesian optimization (BO) methods (Garnett, 2023) that are suitable for expensive-to-evaluate black-box objective functions.
Moreover, the ability to handle inequality constraints and evaluation failures during the search for the best model structure plays an important role in system identification. Inequality constraints can be used to limit the computational budget available for the model, help mitigate overfitting, and prevent models that violate physical principles. In addition, the selection of some structures might lead to invalid composition numerically unstable simulations (Chakrabarty et al., 2021). Such failures or 'crashes" prohibit the assessment of the model performance (Bachoc et al., 2020).
We address the task of model structure selection as a combinatorial BO problem and propose a method specifically designed for models that have additional requirements that can be translated into incquality and crash constraints in the optimization. Our method elficiently searches over potential model structures and is scalable to many design parameters. The proposed method employs expressive surrogate models suited for combinatorial spaces and implements an acquisition function that handles inequality and crash constraints. While our contribution focuses on model discovery for dynamical systems, the proposed method is general and can be used for a wide variety of machine learning tasks.
We empirically evaluate the method on symbolic benchmark problems for equation discovery of nonlincar dynamical systems (Brunton et al., 2016; Mangan et al., 2017). We extend these benchmarks to include inequality constraints in the form of Lasso regularization to combat overfitting, and consider crash constraints
2

---

Published in Transactions on Machine Learning Research (03/2024)
tackling mumerical instabilities and failures. As real-world application example, further optimize knowledge-driven formulation of a multibody driving simulation model.
Our contributions are as follows:
We propose a general formulation for model structure selection as part of a constrained combinatorial optimization problem that considers inequality and crash constraints as a way to handle additional model requirements during the model search.
We present CBDSS, constrained combinatorial Bayesian optimization algorithm that efficiently solve the model structure selection problem and is capable of solving combinatorial optimization problems with up to 10 combinations. CBOSS employs an efficient-to-com co-mpute acquisition function and combines recently proposed surrogate models for discrete domains.
We show empirical evidence that the proposed algorithm outperforms state-of-the-art methods on benchmark problems for oquation discovery and further provide a real-world application example where CBOSS can build digital twin for driving simulator from a large library of modules, while ensuring its real-time capability.
This article continues follows: in Section 2 provide general problem formulation for structure selection. The related work is discussed in Section 3. The proposed combinatorial BO method is described in Section 4. In Section 5, we investigate the performance of our method with various system identification problems and compare to other approaches.
The code for the optimizer and the benchmark problems are publicly http://gghhnb.coo//uccss/////// Molee-Srrnnttrreeeecctooonn
available
2 Problem Statement
We address this structure identification problem constrained combinatorial optimization problem. We require a library of model templates X' and a symbolic-rorressonn knowledge-driven procedure for composing the selected templates. The selection of model templates is to be parameterized with a vector of categorical decision variables denoted as æ E X, defined over the combinatorial domain X' = X2 x X2 x ... x X. with categorical decision variables with respective cardinalities kı....,ku We introduce an objective function f : X' that maps the decision variables a real value indicating the model performance. In practice, this function represents both parameler estimation and model evaluation depicted in Figure and, given noisy time observations, measures the performance of the model structure specified by æ X. Similarly, we define M inequality constraints 9j X' + R with function as h X' + 10. The incquality constraint can be used, for example, to restrict the computational budget, the model complexity or as a regularization, while the binary equality constraint is used to indicate evaluation failures.
MA....,,, and binary equality constraint
We formulate the structure selection problem as the search for the global optimizer x" that fulfills:
x: = arg min aE
f(x)
(1)
s.t.
g, g,(x)) S h(x)) =1
Vj E (1,...,M)
The functions f, 9: and h are all expensive-lo-evaluuee black-box functions and can only be obtained simultaneously. The functions and are noisy and can only be assessed when the experiment is successful, i.c.
(v.c,1)=
(f(x) cgr g(x) + ces (0,
if x is evaluation success
(2)
0
0) if x is evaluation failure
where the noise e, ~ N(0003) and ~ NiO,, diag (2=) are i.i.d and normally distributed. All past observations are collected in dataset D ((x(1) (() c()) 1(0 H where N is the number of experiments.
3

---

Published in Transactions on Machine Learning Research (03/2024)
3 Related Work
In this section, we review some of the existing literature of model structure selection and highlight the potential of employing combinatorial Bayesian optimization promising approach.
Model Structure Selection. Traditional methods for structure selection involve exhaustive search over the space of potential model structures, such as random scarch, A search, and simulated anncaling (Bertsimas & Tsitsiklis, 1993) algorithuns. Others employed evolutionary algorithms (Schmidt & Lipson; 2009; Tanevski et al., 2020) or causal inference (Baumann et al., 2022). However, the main limitation of these methods is their sample efficiency, as the number of potential model structures grows exponentially with the number of discrete variables. Generally, these algorithms are not suitable for expensive-to--eallaae system identification.
A popular non-porametric class of models are kernel methods such as Gaussian processes. Here, the problem of model selection includes the choice of kernels and their composition. These problems have been tackled by greedy search (Duvenaud et al., 2013), a mixture of random walks and BO (Malkomes et al., 2016), and MCMC sampling (Gardner et al., 2017).
Mangan et al. (2017) proposed a sparse symbolic-regression structure selection approach using SINDy (Brunton et al., 2016) as a tool to down-sample parsimonious models from the combinatorially large model Reducing the search space allowed them to brute-force and search for the best model among the remaining ones. SINDy works by recursively performing lincar regression and pruning terms with small model coefficients up to a certain threshold. Although very efficient, this method, as well as other backward- elimination algorithms, over-exploits the model space and does not guarantee to find the optimal one (Guyon Elisscoff, 2003). In addition, SINDy cannot handle constraint problems and is not applicable to knowledge-drieen problems. Another limitation is that the coefficients are pruned based on the regression results obtained by fitting& a prediction model, which is not always the best strategy if the model is to be used as simulation auto-regressive model. BO methods are more general and allow for arbitrary cost functions, whereas SINDy decides on the structure based on the least squares cost.
Combinatorial Bayesian Optimization. The two main challenges in combinatorial Bayesian optimization are the development of (i) better probabilistic regression approaches that can capture the complex interaction between discrete variables and (ii) acquisition function optimizers that efficiently search the combinatorial space using the surrogates. While most methods have been defined over a mixed space of continuous and discrete variables,
focus here on the specific characteristics that are useful for discrete domains. SMAC (Hutter et al., 2011) leveraged tree-based surrogate models and used random walks to obtain local optimum of acquisition function. This method can be applied to mixed input spaces but neglect high-order interaction between variables. BOCS (Baptista & Poloczek, 2018) encoded categorical variables in a combinatorial one-hot binary domain and used polynomial features within sparse Bayesian lincar regression (Carvalho et al., 2010; Makalic & Schmidt,, 2015) combined with Thompson sampling to express the interaction between variables and their effect on the objective function. They optimized the acquisition function using simulated annealing and semi-dcfinite programming. Dadkhahi et al. (2022) extended BOCS and proposed a more compact but still complete and unique encoding that results in fewer monomials. Hase et al. (2018) and Häse et al. (2021) combined Bayesian neural networks and density kernel estimation as a surrogate for BO with categorical variables. Wan et al. (2021) proposed the Casmopolitan algorithm and used a modified Hamming kernel as part of Gaussian process regression, which defines the correlation between two inputs by their Hamming distance the combinatorial graph. To avoid over-expooration due to high-dimenniooal combinatorial spaces they adapted the trust region algorithm from TURBO (Eriksson et al., 2019) to explore only locally near the best location found so far. Further, Oh el al. (2019) proposed COMBO, which used the discrete diffusson kernel built. from the graph Cartesian product of discrete parameters and is able to model high-order interactions between variables. COMBO was improved by Deshwal et al. (2021) (HyBO) by using the closed-form of the discrete diffusion kernel proposed by Imre (2002). In order to tackle high-dimensional combinatorial problems Deshwal et al. (2023) propose an embedding of the high-dimensional discrete search space inlo a low-dimensional continuous one using dictionaries.
Constrained Bayesian optimization is an active area of research for continuous domains (Gardner et al., 2014; Eriksson& & Poloczek, 2021; Ungredda & Branke, 2024; Marco et al., 2021) but rarely investigated for
4

---

Published in Transactions on Machine Learning Research (03/2024)
discrete variables. Daulton et al. (2022) proposed Probabilistic Reparamcterization (PR), the first method for constrained combinatorial BO, which reparameterized the discrete acquisition function optimization problem by introducing discrete probability distributions defined by continuous parameters. This allowed them to optimize the AF using gradient-based methods. We also note that other methods, such as Papalcxopoulos et al. (2022), exist for constrained combinatorial optimization but only consider inexpensive and white-box constraint functions, which is not applicable our problem.
Our proposed method extends to constrained problems and has relatively small computational footprint. Most of the existing methods have high computational demands and therefore scale poorly with the mumber of inputs and data points, since they rely on expensive MCMC methods to infer the posterior distributions. We focus on reducing the high compulational overhead of the state-of-the-ar methods by relying on closed-form models and AF, avoiding expensive approximations.
4 Methods
In the context of structure model selection, Bayesian optimization can be employed to find the optimal model structure that minimizes a chosen evaluation metric. Bayesian optimization is particularly useful when the objective function is expensive-to-evalnate black-box function. Based on past observations, the BO algorithm constructs a probabilistic surrogate model of the objective and the constraints, typically Gaussian Processes (GPs) (Rasmussen et al., 2006). The surrogate model is iteratively updated with new evaluations, allowing for the soquential exploration of the search space. The selection of the next promising evaluation point is guided by an acquisition function, based on the surrogate model's 's predictions and their uncertainty to balance the exploration-exploitttioo trade-off.
Model Structure Selection
Measurements
Rayesian Optimization
Model
Librrr
model structure
Propose next
model stracture candidates x*
Pammette Estimation
Simulation & Performance Evaluation
argmax FRCHEI(X)
Train surrogates
Objective Constraints
model evaluntions (x","")"",1")
Failure
The proposed method for model structure selection using combinatorial Bayesian optimization. Every iteration, the model structure selection algorithu proposes a new model structure æ from the model library. Using measurements, the model's 's parameters are estimated and the composed model performance is evaluated according to objective y f(x) + Ey, inequality constraints c g(x) Eg and crash / h(x) function evaluations.
Figure 2:
We call our method Constrained Combinatorial Bayesian Optimization for model Structure Selection (CBOSS). In the following subsections, we present detailed description of this method. We define (i) probabilistic surrogate models for regression, which will approximate the black-box objective and constraint functions, (ii) probabilistic model for classification, which will approximate the black-box binary function for failures and finally (iii) the acquisition function. An overview of the method is depicted in Figure 2 and in Algorithm 2.
4.1 Probabilistic Surrogate Model for Regression
We consider the setting, where objective function observations are given by y()) f(x(())+ and assume the observation noise 2 ~ NO of) to be an i.i.d. Gaussian random variable with unknown variance of (ef. Equation 2). After making observations, the joint distribution of observations is
y X. oz ~ N(f(X), NTT,).
(3)
5

---

Published in Transactions on Machine Learning Research (03/2024)
where X = (x.)......
....., MiJT XN, f(X)
M() E RN,, y= [v() [y()..... ..... y( R and
IN E NNN is the identity matrix. We assume the objective and constraints to be conditionally independent given X and use the same setup for each j-th dimension of constraint function evaluations oj g1 (x()0) + Ecj Similarly, Cj X ~ N (X), N..N).
(1)
(0)
We now present a surrogate model for the expensive function f a Gaussian process prior over
f fl ~ GP (m, k)
(4)
where the mcan function m : X + and the kernel h: : X x X + are the crucial ingrodients in GP regression since they define the kind of structure that will be captured by the regression model. These components are be defined in the following subsections.
The posterior predictive distribution of function value S(X) at test points X is given analytically (Rasmussen et al., 2006)
f(X') | X X. ~ NN
(5)
μj = m + K(X', XXXXXX +No3)) (x-m(X))
with
Σ₁ 1XKKX,XXX - XX) (NO3)
where m(X): = m(X)) and ('), XJ).
One of the limitations of Gaussian process regression is the light tails of the predictive distribution, which are not robust to outliers (Duvenaud, 2014) and might not capture well enough the discrepancy of the surrogate model (Baptista. & Poloczek, 2018). To allow for better robustness and characterixation of the uncertainty in the prediction while still allowing for a closed-form posterior, we provide formulation of surrogate regressors with hicrarchical priors. A hicrarchical Baycsian model is a model which the prior distribution of some of the model parameters depends on other parameters, which are also assigned prior. We place the following hierarchical GP prior over
oy GP (m, kg oy r (v, om),
(6)
(7)
distribution. To ease the derivations and the hyperparameter tuning, we use more intuitive inverso-gamne paramcterization, whose shape is a = v/2 and scale 6 = ot 2/2 (Taboga, 2017). This parameterization directly reflects the mean and the variance of this distribution since B[1/o2]] = 1/02 and Var[1/o]] = =2/(va Therefore 1/02 can be interpreted as the best guess of the regression precision, while expresses the degree of confidence about the precision. We follow Chen et al. (2023) and place hyperprior distributions over I(Su, and o , Eo2.).
which leads to p(f(X), oz) being inverse-gamma saamma-normal
The predictive prior distribution at X can be caleulated analytically by
f(x) < f x) p(o do T'v, m(), X(X,X)),
(8)
where T(v, 14, Σ) is the multivariate t-distribution, with degrees of freedom v, mean and scale matrix Σ. Furthermore, the predictive posterior distribution at test points X' is also given in closed-form by
f(X') X y, ~ / p(f(x') x', y. o p(o7 X. do = (4. En)
(9)
with
Dy D+N N
f = +
ΣS The (K(',X') -XX, k(X,)))
02 (vo + (y m(X)) (X(X,
(y-m(x) (0- (v+N).

---

Published in Transactions on Machine Learning Research (03/2024)
N
Ev
0g
0₂
fx
N (famme Ment
Eor
A (fumls al Ol)
I (+3)),a)
N (pu
(a) Gaussian Process Regression
(b) T-Process Regression
Figure 3: Graph model of two probabilistic regression models. The array 0g refers to the kernel function k hyperparameters.
which defines the so-called t-Process (TP) regression framework (Shah al., 2013; Traccy & Wolpert, 2018). Interestingly, one can show that as N goes infinity, the TP prior and posterior tends to the GP distributions. Moreover, the equations above show that the parameters v, o, defined in the Gamera distribution of the prediction uncertainty of directly impose how heavy-tailed the predictive distribution will be. The probabilistic graph models for both GP and TP are provided in Figure 3.
4.1.1 Kernel Design
Kernels play a crucial role in GP and TP regression because they define the assumptions and properties of the underlying probabilistic model. The kernel defines correlations between candidates and is required for accurate and meaningful predictions. Notably, common kernels found the literature for continuous spaces rely on natural order of candidates, which is not given in categorical spaces. To address this limitation, many kernels specialized for categorical spaces have been recently proposed in the literature (Imre, 2002; Wan et al., 2021; Deshwal et al., 2021: Oh et al., 2019; Baptista & Poloczek, 2018; Dadkhahi et al., 2022).
In particular, the automatic relevance determination (ARD) discrete diffusion kernel became recently popular and is the discrete analog to the diffusion kernel over continuous spaces, aka the radial basis function (RBF) (Imre, 2002; Deshwal et al., 2021). This kernel defines diffusion over the entire discrete space represented by combinatorial graph, where two nodes are connected if two configurations differ in exactly one variable. It. is given by
d II
R:
8(m,,a())
kain( x²)
e
(10)
1 (k: 1=1
-
where k; is the cardinality of the i-th categorical input variable, 8(2; m is equal zero if x; is equal æy and equal one otherwise, and 3, hyperparameters that control the importance of the i-th discrete variable. Intuitively, this kernel is well suited for globally capturing interactions between configurations. However, one of the main limitations is that it captures similarity based on exact matches and does not consider the degree dissimilarity between categories, which may lead to a loss of information when trying to model relationships or interactions between categorical variables.
Alternatively, Baptista & Poloczek (2018) proposed polynomial kernel that was used in Bayesian linear regression scheme. They observed that. universal surrogale model for binary discrete input. domain X = {0, 1 is given by f(x) Esez as Ilies Ti, where X is the power set of the domain, ag is the coefficient vector, and z E X' is the one-hot encoding of the categorical input æ. In fact, this model describes all the possible configurations and becomes quickly impractical due to the exponential number of monomials.
7

---

Published in Transactions on Machine Learning Research (03/2024)
The idea presented is to truncate monomials up to certain degree
f(z) ao Σ Eaja +
Σ sszi
E
aibkzaaakt +...
(11)
ijiikss
which can be naturally modeled as kernel
d
Kpoly(e, (2, x²) = =a 1 Σ
2 +
Σ 1,21
وHE J -
Σ Tba di kJ>iド>
+...) +.
(12)
where o is the hyperparameter related to the prior variance of the coefficients. Notably, this kernel is more complex than the discrete diffusion kernel, especially when considering higher polynomial degrees, since it can capture higher-order interactions between categorics. However, as the degree increases, the number of parameters becomes prohibitively high. By design, it also suffers from the combinatorial explosion of the discrete domain. Note also that this kernel distinguishes not only the degree of dissimilarity between different. categories but the interaction between them.
In this paper, we propose to combine the strength of multiple kernels (Duvenaud, 2014) in the following way kpolydir( x x') = A (kpooly x²). kann x²)) + (1 - A) (kipoly (x, x') + (aia(a, x'))
(13)
where A ~ (axa 3x) E [0, 1] is a typer-parameter that controls whether the kernels should be added
multiplied together. The idea is to combine the global approximation properties of the discrete diffusion kernel with the high-order interaction and non-stationary properties of polynomial kernels. In the appendix Section A.4, we provide an ablation study of different kernels applied to our benchmark problems.
4.1.2 Hyperprior-parameen estimation
To optimize the hyperparameters, we resort to empirical Bayes, as suggested by Chen et al. (2023), to avoid expensive MCMC samplings. This is possible because both GP and TP regression methods have tractable marginal likelihoods. We achieve this by maximizing the marginal a-posteriori (MMAP) for the prior hyperparameters
0- arg max ply X, 0) II I(0.),
(14)
where the marginal likelihood p(y)X, T(v of + IN and the hyperprior distribution p(O:) is assumed conditionally independent and is specifically designed for each hyperparameter 0,. Hyperpa- rameters without defined hyperpriors are treated as non-informative and are set to the uninformative prior p(0:) & 1. Note that we train all model hyperparameters the model, which for the example of the TP model with kpolydift kernel consists of 0 = {u, on Fa. 8,
4.2 Probabilistic Surrogate Model for Classification
When searching for possible structures of dynamical systems, it is to be expected that simulating the dynamical system might fail due to numerical instabilities. Unfortunately, that lead instabilities a- priori. We model the set of feasible candidates black-box functions and learn the failure regions from data using a probabilistic GP classifier. In this way, the Bayesian optimization algorithm can learn to avoid these regions while scarching for the optimum. We use respectively the following prior and likelihood functions
usually not possible to delimit the regions
GP(m(#)) k(x, x'))
(15)
I h(X) II Bernoulli(t; S(h(x.)))
(16)
where S is the sigmoid function, and I are previous failure observations.
8

---

Published in Transactions on Machine Learning Research (03/2024)
We approximate the intractable Bernoulli posterior distribution p(h(X*) X", .)) with Laplacian approxi- mation, which is simple to implement and can approximate the posterior distribution in closed-form. We use the same kernel as in Equation 13. We implement the alternating optimization procedure from Rasmussen et al. (2006) for finding the posterior mode of the Laplacian approximation Nooton-Raphsoo combined with Strong- Arron-WWolfe linear search) and the model hyperparameters (I-BFGS).
4.3 Acquisition Function
The goal of the acquisition function (AF) is to guide the search for the optimal solution in an optimization problem addressing the tension between exploitation and exploration For the constrained problem, the AF needs to balance the potential utility of a candidate and the probability of feasibility and success. The next candidate model (!!) to be evaluated ileration & is the one that maximizes the acquisition function a for the current model:
xx (t) = arg max a(x'). x'EX
(17)
The AF needs to be optimized in every BО iteration. Generally, this optimization problem is the main computational bottleneck in BO which is why we look for an acquisition function that is cheap to evaluate and given in closed-form. In this work, we combine recent ideas in the field of Bayesian optimization and propose a new acquisítion function called Pailure-Robus Constrained Hierarchical Expected Improvement (FRCHEI), defined follows
a(x')= FRCHEI((') = Psuce (x') Janoe n/N Pieas (aansa")) (x')0su a/N -HEI(x'),
(18)
where
Psuce = p(h((()))) = 1
(19)
m
Pleas Preas(x)= = II ((g)(x')) j HEI(x') = E [(x0))((((()) x) [max (0.0-f((x')
(20)
(21)
The idea of this AF is twofold. First, we consider the hicrarchical expected improvement (HEI) (similar to Chen et al. (2023)). The HEI uses the objective function probabilistic surrogate p(f(='))a' X) to quantify the potential improvement in the objective value over the current best feasible solution found so far y y. Samples with lower predicted mean values and higher uncertainties offer the potential for better solution. HEI modifies the traditional EI by replacing the GP surrogate by TP regression, which improves the robustness against outliers while preserving the closed-form solution of the AF, which is given by
more likely to be explored as they
HEI(z')=00 [ tr (rt:0,,0,1 (T+ D9-1 T(rt:p,.011))
(22)
exploitation
exploration
where ,+ -+ - and Φr is the CDF and T is the PDF of the t-distribution, respectively. The parameters Pf. are the posterior degrees of freedom, mean and standard deviation from the postcrior I-distribution at x'. The derivation of Equation 22 is given in the appendix Section A.1.
Second, avoid unfeasible regions by multiplying HEI with the probability of feasibility Pleas, similar to Gardner et al. (2014) and avoid instability regions by multiplying the HEI with the probability of success Psuce, as proposed in Chakrabarly et al. (2021). Following Hvarfner el al. (2022), we scale these probabilities with Bteas n/N and Bence n/N, where n is the current iteration number, N the maximum number of iterations and Beess Bence are hyperparamcters. The idea is to relax the constraints at the boginning of the optimization when we do not have much knowledge about these black-box functions. As the optimization progresses, we increasingly trust these surrogate models. Note that even though unfeasible models never
9

---

Published in Transactions on Machine Learning Research (03/2024)
Objective
Inequality Constraint
2
1.0
1.0
0.5
0.5
2
0.0
-1
0.0
-2
2
X
PJ Uly HET
X
pig (c)
Prona
pig S 0)
Equality Constraint
Acquisition Function
1.0
1.0 1.0
1.0
0.5
0.5 0.3
0.5
0.0
0.0
0.0
0.0
2
2 -2
X
---h -2 X
Pauce p(h 1
FRCHEI = Pauce Phona HEI
Figure 4: Illustration of the FRCHEI acquisition function for a one-dimensional continuous space example. The true unknown functions are depicted as dashed lines, the observations as crosses and the surrogate predictive posteriors mean and confidence region as purple lines and areas, respectively. Note that no observations could be assessed for the objective and constraint functions in failure regions {x h(z) = 0}. The next candidate point x² is obtained by searching for the maximum of the FRCHEI function. Clearly, the surrogate for the constraint functions and h push down regions of the expected improvement HEI which are likely to be unfeasible or fail. Without these two surrogates, the next candidate would be sampled on an unfcasible and failure-prone region, where the HEI is maximized. Note that the search space for the experiments we consider is discrete and high-dimensional, which can not be well displayed. For illustration purposes, we employed the same method but replaced the combinatorial kernels with a squared exponential.
considered the best experiment, they are still useful to improve the predictive posterior distributions (Gardner et al., 2014), An illustration of this acquisition function for a one-dimensional continuous problem example can be seen in Figure 4.
Note that as By + oo Equation 22, converges to EI
EI(x') = ā, [r ΦN (r+:,,,0,1) +
(23)
which is expected since t-distributions converge to Gaussian distributions as the degrees of freedom parameter V/ + c0. When using EI instead of HEI in Equation 18 we call the acquisition function FRCEI.
4.3.1 Acquisition Function Optimization
Since the input space is fully discrete, optimization of the AF can not be done directly using gradient-based methods. Instead, we use simulated annealing (Bertsimas & Tsitsiklis, 1993), which is standard and performant method for unconstrained optimization over discrete spaces. We use the same simulated annealing approach presented in Dadkhahi et al. (2022), since it simple, computationally efficient, and can optimize over categorical variables. For completeness, we present the algorithm again Algorithm
4.4 Algorithm
We present our final algorithm for black- box constrained combinatorial optimization in Algorithm 2. The algorithm requires an initial dataset D(1) = ((x(*) (+)) c( ((0 ) which can be obtained by evaluating the model structure at randomly sampled initial points z(1) has been shown in recent studies (Wan et al., 2021; Deshwal et al., 2021; Müller al., 2021; Daulton et al., 2022) that local exploitation is beneficial,
10

---

Published in Transactions on Machine Learning Research (03/2024)
especially for high-dimensional discrete spaces. When maximizing the AF, in line 8, we make of this idea and set the initial starting point for the SA optimization method as the best feasible sample evaluated so far. In addition we rerun the AF optimization twice to increase the chance of finding better candidates.
A Algorithm Simulated annealing for categorical variables
function SA(objective function f, categorical domain X, starting point x0), annealing scheduler s(t) for t= = to N do
3:
~~ (mif(d) x - (t-1)
4:
5:
x; Softmax (1-f(x; w, w,z=i)/s())ween)
return (t)
A Algorithm Bayesian Optimization for Model Structure Selection
function CBOSS-FRCHEI (initial dataset D(1))
for t = 1 to N do
3:
Train surrogates:
f < MMAP(D(1)
(Equation 14)
0g (M) (Equation 14)
0: SyplaceModeeMMMPPCH
(Section 4.2)
7:
Optimize acquisition function:
x arg m FRCHEI(>, 0, 0y 01 D(t)) (Algorithm 1 with Equation 18)
y. e. EVALUATTAMOOEESSYRUCUUTRRT
10:
Update data D(++1) ( U ((r,y,c,1)
11:
return lowest feasible evaluation x² = arg min, (y: (,, (,) 0, = 1}
12:
13:
function EAAUUAYEMOORLLTTTUUTTREEEE
14: Estimate model parameters, simulate and evaluate performance U.g.h)
15:
if Failure then
16:
return 10,0,1=00 0}
17:
else
18:
return (g.c,1 =
5
Experiments and Results
In this section, we evaluate the empirical performance of CBOSS on set of constrained equation discovery problems for nonlinear dynamical systems well as knowledge-driven configuration of driving simulator. On the benchmark examples, compare against two simple baselines, random sampling (RS) and simulated annealing (SA). We also compare against the state-of-the-art in constrained combinatorial optimization probabilistic reparameterization (PR) (Daulton et al., 2022). In our experiments, we find:
On the equation discovery problems, CBOSS is either competitive with or outperform the other methods consistently. It always yiclds good equations independent of the random soed. After 500 evaluations, CBOSS is approxinately 10x faster than PR in terms of wall clock time, highlighting the reduced computational overhead due to the closed-form inference and AF,
3. Both constrained BO methods CBOSS and PR perform especially well when finding feasible solutions is more difficult.
CBOSS is able to tune a complex driving simulator with approximately 10 possible configurations towards a driver's preference with the constraint that the simulation must be real-time capable.
While other methods such BOCS, Casmopolitan, COMBO, and HyBO would be interesting for this investigation their very high computational demands make them unsuitable for the the type and size of problems
11

---

Published in Transactions on Machine Learning Research (03/2024)
consider here. In addition, the implementations provided by the authors of these methods do not support constraints of any kind.
5.1 Equation Discovery for Nonlinear Dynamical Systems
In this section, we investigate system identification for the set of low-dimensional nonlinear dynamical systems shown in Figure 5. We use some of the benchmark problems and a similar learning setup from Brunton et al. (2016); Mangan et al. (2017). We investigate a range of dynamical systems: a simple nonlinear damped oscillator, a disease transmission model (SEIR), the chaotic Lorenz Oscillator, and the mean-field model for the cylinder wake in reduced coordinates. More details about these dynamical systems can be found in the appendix section A.3.
1.0
Nonlinear Damped Osc. i - ax 3 Byy
0.8
SEIR
0.6
= μ μs - 3 SI - + a)E SI
0.4
Bux euy
0.0
i= E (
Cylinder Wake
Lorenz Oscillator
D.30
p r w x2
30
i = - x + ay
w wa A yz
20
i p x 1 22
10
-02 0.0
0.0 y
- X 2 X
-10
10
1 = 3 + zy
Figure 5: Equation discovery benchmark problems
2(N7)T Nxd be the matrix of N noisy time measurements of all d states, which are available for learning the system. By differentiating the state measurements, we obtain noisy estimates of the state derivatives over time 2(() (z)) + +e, with e e~ -=(00,02)) We seek to find the structure of the governing nonlinear differential equation fode, which yields the best simulation model.
Let 209 = (2 2/ Rd and Z [x(1)
A common practice in regression problems is to assume that each dimension of fode can be represented by a linear combination of features ((Z) E Nxn, weighted by the coefficient vector = E Rayxd, These nf features compose a large library of possible nonlinear candidate functions that can represent the real system. For the benchmark problems hand, we consider that each ODE dimension can be represented by polynomial terms up to degree p, which is indeed true for all the problems. In vector form equation 24 can be solved for the model coefficients E, using least squares estimator.
(1)
(1)
(2 12)
(1) (1)
)
21
Ed
E
Ea
te (24)
(N)
(z (2
(N) (N)
22
e 2))
Considering large feature library enhances the flexibility and possibly the accuracy of the model up to the point that the error between the predictive and measured gradients Ž --6 (Z)) E approaches zero. However, especially in the presence of noise, complex models have the potential to overfit. the data leading to a deterioration of the simulation performance and even instabilities. To address this issue, we want to consider a large feature library but use structure selection to scarch efficiently and globally for the optimal subset of terms that construct parsimonious models, i.e., the least number of terms that significantly reduce the simulation error (Brunton et al., 2016).
12

---

Published in Transactions on Machine Learning Research (03/2024)
Note that the number of monomials is equal to n, d = di (<t;-), and the total number of possible model structure combinations is given by fm 2 which grows exponentially with the number of monomials, which further grows factorially with the number of dimensions and the polynomial degree. An overview of the search space dimension for the benchmark problems is given in Table 1. Notably, the number of models is the order of quintillions.
Table 1: Equation discovery benchmark problems
Dynamical System
d
Tt
Nonlinear Damped Oscillator
2
5
42
4.4e+11
SEIR
3
3
60
1.2e++18
Cylinder Wake
3
60
1.2e+18
Lorenz Oscillator
3
60
1.2e++88
We define the optimization problem in the same form as Equation 1 as follows
1 N
x' arg min sE
f(x)=
logi0 N
Σ
2- (t)
sim (xx) 2011 1|1 Alog2 E x
(25)
s.t. g(x) Σ [x)(( 55
(26)
h(2²)= = 7em (Zeim) (. -(xx),
(27)
(28)
where {0,1)" are the n, dn binary decision variables that select which of the ny features columns of (Z) will be present in the model for each of the d dimensions. À is hyperparameter, 5 is a constraint threshold, and 2eim (2) is the aulo-regressive simulation trajectory of the system. The coefficients Ea] E Rn of the selected system are estimated by solving the least squares problem for each subset of E,, selected by x. The first right-hand-side term Equation 25 quantifies the mean absolute error between the measurements and the simulated system. The second term is a parsimony- prssimony-bbsed criterion that rewards models with fewer monomials. The inequality in Equation 26 describes the Lasso constraint and is used to restrict the model complexity. If the simulation is unstable, then h(x) 0. The final identified simulation model is given by Equation 28.
E(x)= = [
We run our method CBOSS, as in Algorithm 2, with batch size equal two to reduce the computation time. We investigate CBOSS with the two implementations of our acquisition function: FRCHEI and FRCEI. For RS we pick a configuration uniformly at random, and SA is the same as shown in Algorithm 1, with the difference that we set the objective function infinity if the constraint is violated or the simulation is unstable. PR is used with the recommended settings for discrete binary optimization and batch size equal one. We provide the same randomly sampled 50 evaluations and let the optimizers run for additional 450 model evaluations. We execute each method 10 times and report mean performance and standard error.
Figure 6 shows the best feasible objective function found after a number of evaluations. The proposed methods CBOSS-FRCHEJ and CBOSS-FRCEI outperform in terms of both performance and efficiency. After 500 ovaluations, our methods achieve comparable or better results in mere 3-5 hours, while PR takes around 30-42 hours. The RS and SA methods have a negligible computation time, and simulation takes around 4 seconds. Note that our method runs with a batch size of two, meaning that we choose two candidates in each iteration, which halves the computational overhead. In appendix Section A.4 we show that batch size of two is a good compromise between computation time and performance. For the SEIR and cylinder wake systems, SA never improves from the initial feasible solution. This is because the algorithm changes only one variable at a time, which seems to lead to either constraint violations or crashes. This highlights the need for global scarch methods that can escape from bad initial configurations. Both SA and PR do not perform well for the Lorenz chaotic oscillator and are stuck in local minima, where the states converge fast to a stationary point around the oscillation center (ef. Figure 7a and Figure 7b).
13

---

Published in Transactions on Machine Learning Research (03/2024)
Nonlinear Damped Osc.
SEIR
Cylinder Wake
Lorenz Oscillator
Best feasible objective
-0.6
00.05
-0.5
-1.0
-0.8
-1.0
-0.15
-1.5
100
200
300
400
500
100
300
400
50D
500
Prequency of feasible evaluaaioon
1.0
1.0
1.0
1.0
0.6
0.6
0.6
0.6
0.2
0.2
0.2
0.2
0.0
0.0
100
200
300
400
500
100
200
300
400
500
100
200
300
400
300
100
200
300
400
500
Wall time [in]]
10
10a
10
10
109
100
100
100
10-1
10-1
10
10
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
Number of evaluations
Number of cvaluations
Number of evaluations
Number of cvaluations
RS
SA
PR
CBOSS-FRCEI
CBOSS-FRCHEI
Figure Best feasible objective function evaluation, wall-time per evaluation, and frequency of feasible solutions over the number of iterations for different optimization methods applied to the equation discovery problem.
In Figure 6 we show the frequency of feasible points evaluated by all methods. Our proposed AF (Fquation 18) encourages early exploration and, as the classifier and constraint functions are learned, evaluates more feasible points and finds betier ions. The frequency of infeasible evaluations of the baselines further motivates the benefits of constraint BO. If the search space contains many infeasible points the surrogate model can learn these regions and concentrate its search on feasible regions. This leads to better overall solutions. CBOSS-FRCEI yields comparable but slightly more robust results then CBOSS-FRCHEI.
Figures 7a and 7b show how our method consistently outperforms others in optimizing models across multiple runs and benchmark problems, with low variance crucial for expensive onc-time optimization tasks. Even the pure exploratory policy underperforms. A balanced method like CBOSS strikes the right trade-off between exploration and exploitation, effectively navigating model structures and avoiding local minima.
The appendix Scction A.4 provides ablation studies with further insights and analysis of the proposed method. We investigate the performance of CBOSS when not accounting for failures in the AF. We also consider the effect of using different kernels and show that the proposed combination of popular kernels outperforms others for these benchmark problems. Finally, we inspect the effect of varying the batch size.
5.2 Real-time Multibody Dynamics Model Optimization
In this experiment, aim to optimize digital twin configuration of electric vehicle for driving simulator, which is crucial for developing driving dynamics and advanced driver assistance systems. Driving simulators minimize the need for protolypes and shorten development cycles by enabling early testing of car development. under various realistic conditions before road testing has even started (BMW AG, 2020). Achieving realistic driving dynamics in simulators requires hardware capable of replicating the dynamic range and real-time motion cueing, along with precise digital-twin vehicle models. Companies have invested significantly in cutting-edge simulators to meet. these needs, making the design of accurate simulation models pivotal
14

---

Published in Transactions on Machine Learning Research (03/2024)
RS
SA
PR
CBOSS-FRCEI
CBOSS-FRCHEI
1.0
1.0
10
9C2 SC24
9.5
80
60
1993
1300
1966
1000
1966
3300
1900
320
1082 1900
6.0
60
01
20
-11
-8
-18
-12
-0
10
13
(a) Best models among optimization runs.
RS
SA
PR
CBOSS-FRCEI
CBOSS-FRCHEI
1.0
1.0
3.0
1268
1206
390
6.00
603
200
105
300
20
15
-8
-10
20
13
(b) Worst models among oplimization runs.
Figure 7: Comparison between measurements (orange) depicted without noise and simulated trajectories of the best feasible model after optimization (bluc).
15

---

Published in Transactions on Machine Learning Research (03/2024)
element in virtual vehicle development. Figure 8 shows the high-fidelity simulator at the BMW AG driving simulation site.
Figure 8: High-fidelity driving simulator with 9 degrees of freedom, motion area of nearly 400 square meters and a peak electrical power required up to 6.5 MW. Source BMW AG (2020).
The digital twin a multibody dynamics model, which consist multiple entities (bodies, joints, force and control elements, etc.) (Featherstone, 2014). Every entity has multiple possible implementations, and the goal is to choose suitable one from library while retaining the overall model real-time capability. For different use- of the driving simulator, different implementations twin to be adapted on regular basis. We aim to automate this process using proposed BO method. For example, a damper can be implemented as a simple lookup table or as differential equation. While a complex damper model improves simulation accuracy on a bumpy road, a simple lookup table is sufficient. for longitudinal dynamics maneuvers, In summary, the major challenge in system identification for the driving simulator boils down selecting the appropriate templates for each entity order optimize model simulation accuracy while simultaneously constraining the computational complexity to enable real-time simulations. However, this task is not straightforward. Model designers face overwhelming number of potential configurations arising from the combinatorial explosion when selecting templates. Moreover, adding complexity to one component doesn't 't always enhance model accuracy proportionally increase overall the computational complexity.
necessary. This requires the digital
We investigate automatic template selection for kuowledge-driven structure selection with our proposed com- binatorial BO method CBOSS. We model the vehicle using the Simpack' simulation software and parameterize the model with 46 categorical variables that alter the template of different components. The vehicle for one of the configurations is depicted in Figure 9a. All the template cocfficients were identified in preliminary phase, so there is no need for parameter estimator. The discrete design space X' is defined as follows:
5 binary variables switch between rigid and flexible formulation of bodies.
6 categorical variables switch the number of modes used in the modal-order-rednntioo of various lexible bodies.
2 binary variables switch the tire contact model approach.
4 binary variables switch between rigid (joint) and compliant (force element) motor mounts, suspension rods and stabilizer bushings.
29 binary variables switching the complexity of lookup tables (linear and nonlinear) representing the stiffness of various compliant bushing elements.
The structure identification as a constrained combinatorial optimization problem is defined as
xt arg min
f(x:)
(29)
s.l. g(x) = RTI max (x) 1 0
h(x:) = 1,
where the function f(x) : X + R defines the model performance. In our experiment, the objective function reflects the preference of expert drivers for one model over another. The sil-tiim-i (RTT), is defined as
hhttps://www.c http:://www.dd.coom http://wwww....ooeeeetettttt//as
16

---

Published in Transactions on Machine Learning Research (03/2024)
the ratio of the time required to advance one time step in the simulation divided by the simulation time step size. Real-time models must have RTT below 1 in 99% after the first 2 seconds of 'warming-up". Simulation failures arise from certain input combinations that result in numerical instabilities or in an invalid model, for instance, when it induces kinematic loops that can not be handled by an ODE solver.
We provide the evaluation of 30 initial random samples for the optimizer, with a total evaluation budget of 300 evaluations and run the optimization twice. Figure 9b shows the simulation results for the problem hand. It can be seen that the proposed method finds good real-time capable configurations. Figure 9c depicts the RTI signal for the first 10 seconds of simulation for the best feasible model found during optimization among the two runs. We leave the detailed objective and subjective evaluation of the optimized simulation model for future work.
-120
Best feasible objective
1.4
RTI
RTT max
-140
1,2
160
1.0
0.8
-180
9.6
50
100
150
200
250
300
2
10
Number of evalustions
Simulation time [s]
(a) 3D view of the multibody dy- namical model representing the ve- hicle.
(b) Progress of the optimization runs for the best feasible solution over num- ber of model evaluations.
(c) RTI measurement of the final best vehicle found during optimization.
Figure 9: Optimization results for the multibody dynamics problem with the CEOSS-FRCHEI optimizer.
6 Conclusion and Outlook
In this paper, tackle structure identification in dynamical systems under constraints and potential failures. We introduce constraints to manage computational resources and prevent overfitting. Our approach utilizes combinatorial Bayesian optimization to search for optimal model structures elficiently. We present CBOSS, a combinatorial Bayesian optimization algorithm that effectively manages constraints and failures with low computational overhcad.
We encode the choice of the structure of the dynamical system using binary and categorical decision variables. Our method handles these discrete inputs by employing kernels specific for discrete inputs. We combine recent ideas in kernel design and show that the proposod kernel outperforms state-of-tho-art kernels. Our method handles the black-box constraint and failure functions by learning them with Bayesian regression and classification methods and therefore learns to avoid these regions during search. Finally, we focus on scalability up to a large number of discrete decision variables and make design choices that favour the run time. Our surrogate models and acquisition function are evaluated in closed-form. This allows our mothod to optimize problems with up to 1018 possible combinations.
We provide benchmark problems in the field of symbolic-regressoon that provide evidence that our method outperforms other methods for system identification of variety of nonlinear dynamical systems, such disease models, oscillators and chaotic systems. In addition, provide a complex real-world application example of knowledge-driven system identification where the choice of templates of a multibody dynamical system of an electric vehicle are optimized for accuracy, real-time capabilities and numerical robustness. As the proposed method is not specific to the presented applications, it in principle applicable to other constrained combinatorial problems. Investigating its potential and performance on problems, such as network structure optimization (Du al., 2022), wait-and-judge scenario optimization (Campi & Garatti, 2018), and neural structure optimization with pure categorical variables interesting topics for future research.
17

---

Published in Transactions on Machine Learning Research (03/2024)
Other promising future research directions improved surrogate models that capture higher-order interac- tions of variables leverage the structure of dynamical systems. In addition, there is often a lot of additional prior knowledge about the problem at hand that can be incorporated. For example, expert beliefs (Hvarfner et al., 2022) or qualitative knowlodge about the intcraction terms which can be modeled as graphs.
Acknowledgements
Simulations were performed with computing resources granted by RWTH Aachen University under project rwth1409.
References
Francois Bachoc, Celine Helbert, and Victor Picheny. Gaussian process optimization with failures: classification and convergence proof. Journal of Global Optimization, 2020.
Ricardo Baptista and Matthias Poloczek. Bayesian optimization of combinatorial structures. In International Conference on Machine Learning. PMLR, 2018.
Dominik Baumann, Friedrich Solowjow, Karl Henrik Johansson, and Sebastian Trimpe. Identifying causal structure in dynamical systems. Transactions on Machine Learning Research, 2022.
Dimitris Bertsimas and John Tsitsiklis. Simulated annealing. Statistical science, 1993.
BMW AG. BMW Group sets new standards for driving simulation. http://wwww..pee..bbb.ee
negegeen--oooo---feee--eciiiiii- tt---ee-------riii----ssnntattt- starts-worke waartss-orkllanuaggeeeen November 2020. Accessed: 2023-08-04.
Josh Bongard and Hod Lipson. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 2007.
Elizabeth Bradley, Matthew Easley, and Reinhard Stolle. Reasoning about nonlinear system identification. Artificial Intelligence, 2001.
Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 2016.
Marco C Campi and Simone Garatti. Wait-and-judge scenario optimization. Mothematical Programmingg 2018.
Carlos M Carvalho, Nicholas G Polson, and James G Scott. The horseshoe estimator for sparse signals. Biometrika, 2010.
Ankush Chakrabarty, Scott A Bortoff, and Christopher R Laughman. Simulation failure robust Bayesian optimization for estimating black-box model parameters. In 2021 IEEE International Conference on Systems, Man, and Cybernetics (SMC). IEEE, 2021.
Rick Chartrand. Numerical differontiation of noisy, nonsmooth data. International Scholarly Research Notices, 2011.
Zhehui Chen, Simon Mak, and CF Jeff Wu. A hierarchical expected improvement method for Bayesian optimization. Journal of the American Statistical Association, 2023.
Hamid Dadkhahi, Jesus Rios, Karthikeyan Shanmugam, and Payel Das. Fourier representations for black-box optimization over categorical variables. In Proceedings of the AAAI Conference on Artificial Intelligence, 2022.
18

---

Published in Transactions on Machine Learning Research (03/2024)
Samuel Daulton, Xingchen Wan, David Eriksson, Maximilian Balandat, Michael A Osborne, and Eytan Bakshy. Bayesian optimization over discrete and mixed spaces via probabilistic reparameterization. Aduances in Neural Information Processing Systems, 2022.
Aryan Deshwal, Syrine Belakaria, and Janardhan Rao Doppa. Bayesian optimization over hybrid spaces. In International Conference on Machine Learning. PMLR, 2021.
Aryan Deshwal, Sebastian Ament, Maximilian Balandat, Eytan Bakshy, Janardhan Rao Doppa, and David Eriksson. Bayesian optimization over high-dimensional combinatorial spaces via dictionary-basee embed- dings. In International Conference on Arlificial Intelligence and Statislics. PMLR, 2023.
Wei Du, Gang Li, and Xiaochen He. Network structure optimization for social networks by minimizing the average path length. Computing, 2022.
David Duvenaud. Aulomatic model construction with Gaussian processes. PhD thesis, University of Cambridge, 2014.
David Duvenaud, James Lloyd, Roger Grosse, Joshua Tenenbaum, and Ghahramani Zoubin. Structure discovery in nonparametrie regression through compositional kernel search. In International Conference on Machine Learning. PMLR, 2013.
David Eriksson and Matthias Poloczek. Scalable constrained Bayesian optimization. In International Conference on Artificial Intelligence and Statistics. PMLR, 2021.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. Advances in neural information processing systems, 2019.
Roy Featherstone. Rigid body dynamics algorithms. Springer, 2014.
Jacob Gardner, Chuan Guo, Kilian Weinberger, Roman Garnett, and Roger Grosse. Discovering and exploiting additive structure for Bayosian optimization. In Proceedings of the 20th International Conference Artificial Intelligence and Statislics, 2017.
Jacob R Gardner, Matt Kusner, Zhixiang Eddie Xu, Kilian Q Weinberger, and John P. Cunningham. Bayesian optimization with inequality constraints. In International Conference on Machine Learning. PMLR, 2014.
Roman Garnett. Bayesian optimization. Cambridge University Press, 2023.
David Ginsbourger, Rodolphe Le Riche, and Laurent Carraro. Kriging is well-suited to parallelize optimization. Compulational intelligence in expensive oplimization problems. Springer, 2010.
Isabelle Guyon and André Elisseeff. An introduction to variable and feature selection. Journal of machine learning research, 2003.
Florian Hase, Loic M Roch, Christoph Kreisbeck, and Alán Aspuru-Guzik. Phoenics: Bayesian optimizer for chemistry. ACS central science, 2018.
Florian Häse, Matteo Aldeghi, Riley J Hickman, Loic M Roch, and Alán Aspuru-Guzit Gryffin: An algorithm for Bayesian optimization of categorical variables informed by expert knowledge. Applied Physics Reviews, 2021.
Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown. Sequential model-based optimization for general algorithm configuration. In Learning and Intelligent Optimization. Springer, 2011.
Carl Hvarfner, Danny Stoll, Artur Souza, Luigi Nardi, Marius Lindauer, and Frank Hutter. pii8BO Augmenting acquisition functions with user beliefs for Bayesian optimization. In International Conference on Learning Representations, 2022.
Kondor Risi Imre. Diffusion kernels on graphs and other discrete input spaces. International Conference Machine Learning. PMLR, 2002.
19

---

Published in Transactions on Machine Learning Research (03/2024)
Lennart Ljung. System identification. In Signal analysis and prediction. Springer, 1998.
Enes Makalic and Daniel F Schmidt. A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters, 2015.
Gustavo Malkomes, Chip Schaff, and Roman Garnett. Bayesian optimization for automated model selection. In Proceedings of the Workshop on Aulomatic Machine Learning. PMLR, 2016.
Niall M Mangan, J Nathan Kutz, Steven L Brunton, and Joshua Proctor. Model selection for dynamical systems via sparse regression and information criteria. Proceedings of the Royol Society A: Mathematical, Physical and Engineering Sciences, 2017.
Alonso Marco, Dominik Baumann, Majid Khadiv, Philipp Hennig, Ludovie Righetti, and Sebastian Trimpe. Robot learning with crash constraints. IEEE Robotics and Automation Letters, 2021.
Sarah Müller, Alexander von Rohr, and Sebastian Trimpe. Local policy search with Bayesian optimization. Advances in Neural Information Processing Systems, 2021.
Duy Nguyen-Tuong and Jan Peters. Model lcarning for robot control: survey. Cognitive Processing, 2011.
Changyong Oh, Jakub Tomczak, Efstratios Gavves, and Max Welling, Combinatorial Bayesian optimization using the graph cartesian product. Advances Neural Information Processing Systems, 2019.
Theodore P Papalexopoulos, Christian Tjandraatmadja, Ross Anderson, Juan Pablo Vielma, and David Belanger. Constrained discrete black-box optimization using mixed-integer programming. In International Conference on Machine Learning. PMLR, 2022.
Carl Edward Rasmussen, Christopher KI Williams, et al. Gaussian processes for machine learning. Springer, 2006.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 2009.
Amar Shah, Andrew G Wilson, and Zoubin Ghahramani. Bayesian optimization using student-t processes. NIPS Workshop on Bayesian Optimization, 2013.
Marco Taboga. Lectures on probability theory and mathematical statistics, 2017.
Jovan Tanevski, Ljupčo Todorovski, Yannis Kalaidzidis, and Saão Džeroski. Domain-specific model selection for structural identification of the rab5-rab7 dynamics in endocytosis. BMC Systems Biology, 2015.
Jovan Tanevski, Ljupèo Todorovski, and Sašo Džeroski. Combinatorial search for selecting the structure of models of dynamical systems with equation discovery. Engineering Applications of Artificial Intelligence, 2020.
Brendan D Tracey and David Wolpert. Upgrading from Gaussian processes to student's t processes. In 2018 AIAA Non-Deterministte Approaches Conference, 2018.
Juan Ungredda and Juergen Branke. Bayesian optimisation for constrained problems. ACM Transactions Modeling and Computer Simulation, 2024.
Floris Van Breugel, Yuying Liu, Bingni W Brunton, and J Nathan Kutz. Pynumdiff: A python package for numerical differentiation of noisy time-series data. Journal of Open Source Software, 2022.
Xingchen Wan, Vu Nguyen, Huong Ha, Binxin Ru, Cong Lu, and Michael A. Osbornc. Think global and act local: Bayesian optimisation over high-dimensional categorical and mixed search spaces. In International Conference on Machine Learning. PMLR, 2021.
20

---

Published in Transactions on Machine Learning Research (03/2024)
A Appendix
A.1 Constrained Hierarchical Expected Improvement
In this section, we provide the remaining ingredients to derive the closed-form expression for the acquisiiio function used: the failure-robust constrained hierarchical expected improvement (FRCHET) as defined in Equations 18-22.
Let the hierarchical predictive posterior for the objective and for each i-th constraint function at test point x' be defined as f(x²) f(x') x', y, T(",, of) and (xx)) gx(x") x', y, ~ T(2g Pg. o2) respectively. The definition of the constrained hierarchical improvement function is inspired by Gardner et al. (2014) and is defined as
CHI((') A(x') max y. - f(x')}
(30)
where yt is the lowest feasible objective function observed so far and A(x') Bernoull(()(a')) E {0, 13 is the feasibility indicator function with parameter
Pleas (a') := (ax)
(31)
E (ax'))
(32)
p(ō (x') S ,0...,,mml (x') S 0)
(33)
II ((0)(=)) S 0)
(34)
m
p(u:(ar)) a', X) dg'(x')
(35)
m (0; Pg Mg. o2).
(36)
Note that from Equation 33 to Equation 34 we assumed the simplest case, where are constraints are conditionally independent given x'. In addition, Equation 35 is simple product of univariate cumulative density of the t-distribution कr evaluated at 0, which can be calculated analytically for any x'.
The constrained hierarchical expected improvement is obtained by taking the expectation of Equation 30:
CHEI(=') = E (CHI())))
(37)
E [A(x') E [max(0, y* - f(x')]
(38)
Paeas(x²)) (x') HEI(x'),
(39)
where (xe) is defined as in Equation 31 and HEI can be obtained with the help of the reparameterization trick (Tracey & Wolpert, 2018) with 5 T(p, 0, 1) and 7+ = (y* py)/oyy
HEI(x') (100))))((()) xx) mmax((. (0. y+- [xx')]
(40)
led (y* - ((x')) ((((='): Py, o7) df (a'))
(41)
(y+ - Tos - ),) T(T Dy, 07) dr
(42)
(y+ - m4) 17 (r (rt: 01 1 TT(T: V.. dr
(43)
o, + 1 (rt: Dy, 0,1) 3, ( Dj
11 - )
(1 (1+ r+2 )
(44)
1/
1 + +
=ō1 7 + Φr (7t: 0,1)
T(rt By 0,10
(45)
Dj 11
21

---

Published in Transactions on Machine Learning Research (03/2024)
The other minor modifications to turn CHEI into FRCHEI are defined in Section 4.3.
A.2 Hyperparameters
For each optimizer, we use the same hyperparamettess across all experiments. These hyperparameters have been hand-tuned to reasonable values so structure selection works well for a wide variety of problems. Table 2 shows the hyperparameters used for the optimizers investigated in this work. For the SA algorithm, we use the same temperature scheduler as the one used for AF optimization in CBOSS. The PR method has been used with the default, parameter values suggested by Daulton al. (2022). We set the number of MC samples to 256 as provides a good tradeoff between performance and wall time. We follow their recommendation and disable trust region over the Boolean variables.
Optimizer hyperparameters used for all experiments. The upper-seript (0) denotes the initial value used for hyperparameter optimization. In CBOSS, the same regression hyperparameters to learn the objective and inequality constraint functions.
Table 2:
Alg.
Hyperparameter
Reference
Values
CBOSS
TP regression : o hyperprior
Equation
o [22, 15), 0m 0.15
CBOSS
TP regression: : v hyperprior
Equation
V (22,0.5). /10) 2.5
CBOSS
GP regression: : of hyperprior
Equation 3
o² [(2,, 15). 22 = 0.15
CBOSS
kernel discrete diffusion
Equation 10
(000) for i=1
CBOSS
kernel difffsion
Equation 13
A Beta 1.5), A 0.5
CBOSS
AF exponents
Equation 18
Breas 20, Psuce
CBOSS
AF optimizer: SA scheduler
A Algorithm
s(1) = 0.5 0.95'
CBOSS
AF optimizer: SA max. iterations
Algorithm 1
N = 100
CBOSS
AF optimizer: SA reruns
Algorithm 2
4
SA
SA scheduler
A Algorithm
s(1) = 0.5 0.99*
PR
MC samples
Daulton al. (2022)
256
Figure we show the hyperpriors and initial values used in CBOSS. For all the other hyperparameters, we use uninformative hyperpriors. Given B[[/o]] = 1/a and Var[1/o}] 2/ 2/(uom)) we interpret 1/0 as the estimated regression precision and v as the confidence in this precision. We thus select Gamma hyperpriors for o and with positive support. Before training, normalize the target outputs to standard deviation of one, and set Gamma hyperprior for o with mode at 0.2, indicating 20% of the data variance, while assigning low likelihood for o exceeding half the variance. We assign distinct Gamma hyperprior to v to reflect plausible values for normalized data. For the À variable in Equation 13, which balances kernel multiplication and summation, we choose a Beta distribution with a mode at 0.5 and a broad density range. These hyperpriors represent our beliefs about. the latent variables' true values and guide the optimizer in learning these hyperparameters from the data, detailed in Equation 14.
of zero and
5.0
1.0
3 0.1
0.5
0.0
0.0
0.0
0.0
0.5
1.0
10
20
0
1
2
(a) 2 11(,, 15),
2
20)
-0.15
(b) / 122, 0.5),
(0)_ -25
(c) À (1t5 1.5),
( 0.5
Figure 10: Informative hyperpriors used in CBOSS.
22

---

Published in Transactions on Machine Learning Research (03/2024)
Figure 11 shows the modified probability of feasibility and success, seen that the exponents 3 n/NN relax the constraints at the beginning of the optimization when there is no much knowledge about the inequality and equality black-box functions. As the optimization progresses, we increasingly trust these surrogate models and increase the impact the probability of feasibility in the AF.
in Equation 18. It can be
Figure 12 shows the simulated annealing scheduler s(t) used Algorithm as a function of the number of iterations. In this paper, we use the exponential decay scheduler s(t) To d', where To the initial temperature, d is the temperature decay and is the iteration number. We calculate d such that the temperature decays from To 0.5 to the final temperature of T. 0.005 at the final iteration.
Beeas 20
Bauce
Modifidd probability
1.0
1.0
n/N [%]
0.50
0.5
0.5
0.25
10
0.0
0.0
50
0.00
0.0
0.5
1.0
0.0
0.5
1.0
100
250
500
Probability of feasibility
Probability of sucoess
Iteration
Figure 11: Modified probability of Teasibility and success Pmod pan/N as a function of the probability of feasibility P=Pfeas or success P-P the current iteration n and total number of iterations N.
Figure 12: SA scheduler with ini- tial temperature To=0.5 and de- cay d 0.99.
A.3 Experiment Details
Table provide the experiment details required to reproduce each benchmark problem in Section 5.1. The parameters Oode refer to the differential equations for each benchmark problem, as depicted in Figure 5. The available noisy measurements for parameter estimation and model evaluation consists of single simulation run, starting at an initial state æfto 0), simulated with fixed simulation time step size A, up to the stop time ly. The state observations are corrupted equally for each dimension with zero-mean Gaussian noise with standard deviation ssn- All measurements have been filtered using the Total Variation Regularized Denoising technique (Chartrand, 2011; Van Breugel et al., 2022) with the same regularization parameter Ty = 0.01 and number iterations T, = 10.
Table Benchmark problem parameters
Dynamical System
Dde
2 (tu)
A
t
Nonlinear Damped Osc.
a 0=0. 0=0.1,
22 0]
0.01
35
0.1
SEIR
1-, w - 1,
[0.9995, 4e-4, 1e-4)
0.1
150
0.01
Cylinder Wake
[0.001, 0, 0.1]
0.1
100
0.01
Lorenz Oscillator
[10, 10, 10]
0.01
20
All the benchmark problems have been simulated with a fixed-step RK45 solver, where the 5th stage has been used detect numerical instabilities, which stops the simulation and reports the failure to the main program. We provide selfimplementation of this solver that performs simulations in parallel and is useful to speed up computation when many simulations are executed with the same integration time steps but with different model parameters.
The benchmark experiments ran on Intel Xeon Platinum 8160 Processors "SkyLake" at 2.1 GHz, on 4 isolated physical cores. The optimization and the simulations of the multibody dynamical problem were performed on 2x Intel Xcon Gold 6256 3.6 3.6Hz computer, running a RedHawk Linux RTOS (by Concurrent Real-Time) The first 12 cores were dedicated to the optimizer, while the remaining 12 cores of the second CPU socket were shielded and dedicated to simulations.
23

---

Published in Transactions on Machine Learning Research (03/2024)
A.4 Ablation Studies
In this section, we provide ablation studies that reinforce the design choices made for the proposed algorithm.
Figure 13 shows simulation trajectories for the best feasible and stable models found over the progress of the optimization for the CBOSS-FRCHEI method. It can be seen that the CBOSS-FRCHEI algorithm quickly improves the model simulation performance. The second image from left to right. shows the simulated trajectory of the best model found after 163 model structure evaluations, which is already very close to the measurement trajectory. This illustrates well the efficiency of the proposed method in terms of the number of evaluations needed for system identification.
0.10
0.100
0.100
0.100
0.075
0.375
0.075
0.04
0.050
0.050
0.02
3.025
0.025
100
3.00
0.2
0.2
0.2
0.2
-0.2
-0.2
-0.2
0.2 y
-02
0.0
0.0
Y
x 0.2
0.2
0.2
Figure 13: Comparison between measurements (orange) depicted without noise and simulated trajectories (blue) for the best feasible model found in the course of optimization for the Cylinder Wake problem and CBOSS-FRCHEI. From the left to the right are depicted the models at iteration 8, 163, 200 and 436, respectively.
In Figure 14, we compare CBOSS-PRCHEI to CBOSS-FRCEI and CHEI, two variants of the proposed method. The first is similar to CBOSS-FRCHEI but uses a GP instead of hierarchical TP prior for the regression surrogates. The second does not account for failures in the AF. It can be seen that FRCHEI outperforms CHEI for all benchmarks, which agrees with Chakrabarty et al. (2021) and is indication that learning failure regions improves the optimizer. As expected, CBOSS-FRCHEI has a higher success rate than CBOSS-CHEI for all benchmark problems. Evaluating failures is a waste of resources because no value for the cost function is obtained, and the candidate configuration is never considered the best experiment. In addition, this data point does not contribute to training the surrogates due to the missing target value. For this reason, the optimizer might get stuck into failure regions, where the objective function surrogate considers there is a substantial expected improvement.
It can be further seen in Figure 14 that FRCEI provides similar performance and sometimes slightly outperforms FRCHEI. The idea behind hierarchical priors is to be able to better describe the model discrepancy towards outliers while preserving the closod-form solution of the acquisition function. One side effect of this approach is the inflation of the predictive posterior uncertainty that might lead to over-exploration of the design space, which is critical for Bayesian optimization in high-dimensional discrete problems. Nevertheless, we still offer TP regression formulation in our framework since it it does not increase the implementation complexity and computation time, and may work better for other problems.
Moreover, we investigate the effect of varying the batch size in Figure 15. We choose the batch size equal 2, since it provides a good trade-off between optimization performance and wall time. Throughout this paper, all the other experiments with CBOSS have been conducted with batch size equal 2. It is also important to note that we do not employ more involving batch evaluation strategies, such as Kriging believer strategy (Ginsbourger et al., 2010), which can potentially incrcase even more the performance at the cost of incrcasing the wall time.
Finally, in Figure 16, we compare the performance of the proposed method with different kernels. Evidently, combining the degree 2 polynomial kernel with the discrete diffusion kernel as in Equation 13 provides the best results. Note that the wall time for the polynomial kernel scales almost linearly with the number of samples but is expensive due to the number of multiplications that arise due to the large number of features resulting from the sccond-order polynomial combination of inputs. The discrete diffusion kernel results in a lower wall time at low number of evaluations but increases at a higher rate. This is due to the complexity of calculating pair-wise delta functions between all inputs in the dataset. In terms of computational complexity,
24

---

Published in Transactions on Machine Learning Research (03/2024)
Nonlinear Damped Oscillator
SEIR
Cylinder Wake
Lorenz Oscillator
Best feasible objective
-0.6
-0.05
-0.5
-1.0
-0.8
-0.15
-1.5
100
500
100
50D
Frequency of feasible evaluaaionn
1.0
1.0
1.0
1.0
0.6
0.6
0.6
0.6
0.2
0.2
0.2
0.2
0.0
0.0
100
200
300
400
500
100
200
300
400
500
100
200
300
400
300
100
200
300
400
500
Wall time [min]
10a
10a
10
10
10"
100
109
100
10
10
10
10
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
Number of evaluations
Number of cvaluations
Number of cvaluations
Number of cvaluations
CBOSS-CHEL
CBOSS-PRCEI
CBOSS-PRCHEI
Figure 14:
Optimization results for different variants CBOSS. CBOSS-CHEI is the variant that does not handle failures, and
is the variant that implements GP instead of TP regression.
the final performance should be evaluated individually for the problem at hand and depends the number of samples and the number of input dimensions.
Nonlinear Damped Oscillator
SEIR
Cylinder Wake
Lorenz Oscillator
0.0
0.0
0.00
-0.5
Best feasible objective
0066
-005
0.3
-1.0
-0.8
-1.0
-0.15
-0.9
-1.5
-0.20
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
.000
50D
Wall time [uin]
101
10
101
101
10a
10a
10
10
10
10
100
200
300
400
500
100
200
300
400
500
100
200
300
400
50D
100
200
300
400
500
Number of evaluations
Number of evaluations
Number of evaluations
Number of evaluations
CBOSS-FRCHET bsl
be2
FBOSS-FRCHEL
Figure 15: Optimization results for CBOSS-FRCHEI method for different batch sizes of 1, 2 and 4.
25

---

Published in Transactions on Machine Learning Research (03/2024)
Nonlinear Damped Oscillator
SEIR
Cylinder Wake
Lorenz Oscillator
Best feasible objective
-0.05
0.5
-0.5
00.7
-0.10
-0.8
-0.9
-1.5
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
Wall time [min]
104
10a
104
104
10
100
10
100
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
Number evaluations
Number of evaluations
Number evaluations
Number of evaluations
CBOSS-FRCHEI Kpolyddf
CBOSS-FRCHELKKRE Kdiff
CBOSS-FRCHEL Kpoly
Figure 16: Optimization results for CBOSS-FRCHEI method with different kernel functions: Kpoly is the polynomial kernel, Kdiff the discrete diffusion kernel and Kpolydiff the combination of the two.
26


Reviews:
[{"review_id": "ok4OghXSs0", "reviewer": "Reviewer_1qh9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, the authors' method seems to achieve good performance, and the authors' experimental evaluation is reasonably thorough. I think the modelling choices are mostly fine -- t-Processes are good models if a bit annoying to do hyperparameter learning with sometimes, and the authors' acquisition function looks reasonable if standard. In particular, I haven't seen very many Bayesian optimization papers using t-Processes, thinking of how acquisition functions might change under t-Process models, or comparing them to standard GP based Bayesian optimization. Therefore, the comparison alone is pretty interesting to me.\n\nBeyond my interest in the methodology, the authors' consider a host of fairly interesting real world optimization problems that seem quite challenging. \n\nFor weaknesses, I have a few comments mostly centering around some of the experiments, and the ablation studies in particular are very confusing to me, and indeed seem to contradict some of the authors' messaging about t-Processes in their methods section (see below).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors' propose a combinatorial Bayesian optimization algorithm with both inequality and equality constraints. The authors' motivating application domain is the construction of models for dynamical systems. To achieve this, they assume access to a set of model \"templates\" that can be combined and then parameter fit. The black-box objective optimized using Bayesian optimization is a measure of goodness-of-fit of the resulting model after both model template selection and parameter estimation, under constraints like resource requirements (inequality) and lack of model failure (equality).\n\nThe authors algorithm consists of:\n- A hierarchical Gaussian process model with an inverse Gamma prior over the outputscale (I think -- see my comment below under \"minor\" about the heavy use of \\sigma^2) which leads to multivariate Student's t-distributed predictive distributions. The authors use a linear mixing of a product of the kernels of Deshwal et al., 2021 and Baptista & Poloczek 2018, and an additive kernel of these same two kernels. \n\n- An acquisition function that is a cEI-like combination of a \"hierarchical EI\" acquisition function (eqns. 21 and 22) with the product of feasibility probabilities under the constraint surrogates.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Ablation:\nMaybe I am missing something here, but it seems like the authors' method using a GP surrogate model (FRCEI) has comparable or (possibly statistically significantly!) better performance on three out of four tasks in your ablation study (Figure 11)? The authors even mention this (\"FRCEI provides similar performance and sometimes slightly outperforms FRCHEI\"), but the authors conclude that t processes are nevertheless a \"more robust and flexible probabilistic model.\" Is there actually compelling evidence of this? Certainly not in Figure 11. I'd really like to see the paper do a better job of actually supporting the authors' own modelling choices.\n\nRelated work and baselines:\nThe authors' related works and literature survey is reasonably good for combinatorial Bayesian optimization methods. The only missing relatively recent work that I'm aware of is BODi (Deshwal et al., 2023). BODi has quite good performance in their own results; however, it does not directly handle constraints I think.\n\nI'm not particularly enthusiastic about the authors' claims that other combinatorial bayesopt methods can't be compared to because they don't mention or handle constraints in their own paper. For many of these methods, the key innovation had nothing to do with acquisition, which is often where constraints are handled. For example, BODi uses EI in a straightforward manner -- there's nothing really preventing the use of cEI e.g. as in Gardner et al., 2014 or constrained Thompson sampling as in Eriksson et al., 2019. \n\nWith that said, PR is a reasonably competitive baseline in this area, and the authors' method does seem to perform quite favorably on a few of the tasks they consider.\n\nMinor:\n- The symbol \\sigma^2 is extremely overloaded. On page 6 alone it is used variously to refer to: the likelihood noise for the objective function (\\sigma^2_y), the likelihood noise for the constraint functions (\\sigma^2_c), an outputscale (\\sigma^2), as what I think is a typo again for the likelihood noise (\\sigma^2 -- should be \\sigma^2_y in equation 5?), a scale parameter for an inverse Gamma distribution (\\sigma^2_m). The likelihood noise is then I think dropped in equation (9)? Can the authors take a quick pass through and clean some of that up, or verify correctness here?\n\n- Wan et al., 2021 technically called their method \"Casmopolitan\" (note the a).\n\n- Algorithm 2 is maybe superfluous? An algorithm box that contains \"train the surrogate model, optimize the acquisition function, evaluate objectives and constraints, repeat\" is essentially just a description of Bayesian optimization, which already exists well enough in the text.", "broader_impact_concerns": "I don't see any significant concerns here.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "gVN61kcSJQ", "reviewer": "Reviewer_doWD", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strength**. This paper is quite well-written, and the proposed method is well motivated. The algorithmic decisions are properly contextualized and appear quite intuitive. I like the idea of the T-Process as well as the use of  categorical decision varibles for selecting strctures. The normalization of the P_{succ} and P_{feas} for n/N is also interesting. In the experiment section, I especially appreciate the inclusion of both sets of experiments, which showcase the discrete model search under investigation and its applicability. \n\n**Weakness**. I don't think the proposed methodology is particularly novel. The main surrogate model, T-Process, is from a previous work; and the combination of P_{succ}, P_{feas} and any version of EI is fairly common. \nThe main strenght of this paper is how well they have placed all the components together to solve an interesting problem. \n\nI am also surprise with the high number of binary and categorical variables used here. The authors used simulated annealing (SA) to optimize the acqusition function. I'm wondering if the authors have consider other methodologies for this optimization or could provide any insights into why the SA proposed by Dadkhahi (2022) is a good choice here.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a novel approach to model structure selection, framing it as a black-box constrained combinatorial optimization problem. Crucially, they use a mix of *categorical decision variables* and *continuous variables* to *select model* templates as well as its *free parameters*. Additionally, the proposed method addresses can also handle inequality and failures constraints, which the authors referred by \"crash constraints\" in the paper. \n\nThe proposed algorithm, CBO-FRCHEI, is designed to efficiently solve complex model structure selection problems. The key components of the proposed methodology are:\n* Using a t-process regression instead of a regular Gaussian Process.\n* Developing a custom kernel that combines two previously proposed kernels for discrete spaces: the discrete diffusion kernel and the truncated polynomial kernel.\n* An acquisition function called Failure-Robust Constraint Hierarchical Expected Improvement that uses Expected Improvement using the t-process regression and multiplies it with the probability of success (not experiencing a failure or crash) and the probability of feasibility (satisfying all the constraints).\n\nEmpirical evidence presented in the paper demonstrates that CBO-FRCHEI outperforms existing state-of-the-art methods in benchmark problems related to equation discovery. The paper extends its contributions to practical applications by showcasing the construction of a digital twin for a driving simulator using CBO-FRCHEI.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- I apologize if I overlooked this, but it seems that the proposed acquisition function has hyperparameters βsucc and βfeas. I couldn't find a discussion on how these parameters were determined. Could you kindly provide this information and offer recommendations to assist readers in understanding their impact?\n- Please consider revising the text to ensure that the novelty of the proposed acquisition function is not overly emphasized.\n- For completeness, consider discussing the following related work:\n  1. Gardner, J.R., Guo, C., Weinberger, K.Q., Garnett, R., & Grosse, R. (2017). Discovering and Exploiting Additive Structure for Bayesian Optimization. International Conference on Artificial Intelligence and Statistics (AISTATS 2017).\n  2. Duvenaud, David, Lloyd, James Robert, Grosse, Roger, Tenenbaum, Joshua B., & Ghahramani, Zoubin (2013). International Conference on Machine Learning, 2013.\n  3. Malkomes, G., Schaff, C., & Garnett, R. (2016). Bayesian Optimization for Automated Model Selection. In Advances in Neural Information Processing Systems. (NeurIPS 2016).\n  4. Raissi, M., & Karniadakis, G.E. (2018). Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, Volume 357, 15 March 2018, Pages 125-141.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ibLelYJaOU", "reviewer": "Reviewer_FHKn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n - The paper uses a long list of state of the art methods in machine learning. The usage of these methods is sound and well-motivated, and elevates above similar papers that use rather standard methods.\n - The state of the art methods used for comparison are strong methods (perhaps except SA).\n - The examples for evaluation are strong examples. Repacing SINDy for learning differential equations from data is a major strong point of this paper.\n\nWeaknesses:\n - Written rather wordy. Not critically, but the text might be shortened significantly. In particular, for a journal that strives for short novelty papers, this paper seems much more as a review of the state of the art in several domains and doing a combination of several known techniques. Case in point (this is everywhere else, too) is the last paragraph of page 14, which I could rewrite as \"Figures 7a and 7b show how our method consistently outperforms others in optimizing models across multiple runs and benchmark problems, with low variance – crucial for expensive one-time optimization tasks. Even the pure exploratory policy underperforms. A balanced method like CBO-FRCHEI strikes the right trade-off between exploration and exploitation, effectively navigating model structures and avoiding local minima.\". This is half as long as the original text.\n - Modeling parameters and model structure both in the discrete space X seems rather restrictive to me. First, because a model is only allowed to have discrete parameters. Second, because the model structure does not change the model parameters, which is rather typical in my experience. (Since I see no reasons why the methods could not be combined with classical methods on continuous spaces, a short comment in this direction might easily fix this weakness.)\n - While there are ablation studies, there could be more. Many choices are done in this paper. However, the paper is rather long anyway, and the choices in the paper look rather natural, I see this as a minor point.\n\nQuestions:\n - Why is h a binary equality constraint function? This might make optimization harder. Of course, there are scenarios where this is the case. Perhaps you can comment on which scenarios you have in mind and why/how you are different from the literature. You approach reminds me of some papers I read about safe active learning, perhaps you can comment on differences and connections.\n - EI as acquisition function is no longer state of the art. Can alternative methods be used instead? (This is a rare occasion, where the choise of the method is not state of the art.)\n - How do you construct priors for hyperparameters specifically?\n \nIn general, I think the content of the paper is publishable. The methods are not completely novel, but a non-trivial and strong combination of rather new methods in a non-obvious manner. The language could be improved.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "A version of Bayesian optimization (BO) is suggested, where inequality constraints (this is somewhat easy) and failure of measurements (this is much harder) are taken into account when applying BO on discrete spaces. The BO approach build on state of the art changes to Gaussian process models, combined with changes to the acquisition function. The approaches is tested on two non-trivial applications (differential equations learning and model selection for car dynamics).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Try to rewrite the text to make it more concise.\n\nMinor:\n - Between (7) and (8): Var -> \\operatorname{Var}\n - I do not think that p(theta)=1 is the correct Jeffrey's prior for all kinds of variables.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This manuscript considers the problem of model structure optimization for modeling dynamical systems, a challenging and important problem. The authors propose to apply techniques from Bayesian optimization to this problem, a natural but not at all straightforward idea. The authors outline a well-motivated and feature-rich framework for applying Bayesian optimization to the problem of interest, and validate this approach through a series of well-designed experiments.

The reviewers agree that the work presented in this manuscript is high quality and of interest to the TMLR audience. They are universal in their recommendation of acceptance.

Throughout the review and discussion period, the reviewers provided feedback and suggestions that the authors have faithfully incorporated into the manuscript. I believe these changes have considerably strengthened the work. I commend both authors and reviewers for their engagement throughout this process.

I believe the manuscript is both suitable for and ready for publication.
Decision: Accept (poster)