The CompModels package for R provides a suite of computer model test functions that can be used for computer model prediction/emulation, uncertainty quantification, and calibration. Moreover, the CompModels package is especially well suited for the sequential optimization of computer models. The package is a mix of real-world physics problems, known mathematical functions, and black-box functions that have been converted into computer models with the goal of Bayesian (i.e., sequential) optimization in mind. Likewise, the package contains computer models that represent either the constrained or unconstrained optimization case, each with varying levels of difficulty. In this paper, we illustrate the use of the package with both real-world examples and black-box functions by solving constrained optimization problems via Bayesian optimization. Ultimately, the package is shown to provide users with a source of computer model test functions that are reproducible, shareable, and that can be used for benchmarking of novel optimization methods.
The CompModels package (Pourmohamad 2020) for R (R Core Team 2020) is a suite of test
functions designed to mimic computer models. Usually deployed when
physical experimentation is not possible, a computer model (or code) is
a mathematical model that simulates a complex phenomena or system under
study via a computer program. For example, weather phenomena, such as
hurricanes or global warming, are not reproducible physical experiments.
Therefore, computer models based on climatology are used to study these
events. At its simplest, a computer model is a mathematical model of the
form
Another typical trait of computer models is that they are often treated
as black-box functions. Here, a black-box computer model is a computer
model where evaluation requires running computer code that reveals
little information about the functional form of the underlying
mathematical function,
The CompModels package can be used to test and develop methods for
computer model emulation (prediction), uncertainty quantification, and
calibration. However, the main focus when developing the package was
placed on building computer models for optimization. Real-world computer
models are often built with the goal of understanding some physical
system of interest, and with that goal usually comes the need to
optimize some output of interest from the computer model. For example,
in hydrology, the minimization of contaminants in rivers and soils is of
interest and so computer models representing pump-and-treat remediation
plans are often used in order to optimize objectives, such as the
associated costs of running pumps for pump-and-treat remediation, while
also ensuring that contaminants do not spread (Pourmohamad and H. K. H. Lee 2019).
Recalling that most computer models are computationally expensive to
run, the need for efficient sequential optimization algorithms (also
known as Bayesian optimization) that do not require many functional
evaluations is high, which is why the focus of the test functions in the
CompModels package is placed on optimization. More specifically, the
CompModels package presents functions to optimize of the following
form
Some of the functions in the CompModels package have known functional
forms, for example, the gram()
and mtp()
functions. However, most
all functions are intended to serve as black-box computer models. All of
the black-box computer model functions within the package are aptly
named bbox
(short for black-box) and followed by a unique integer
value to make the functions discernible. For example, bbox1()
and
bbox2()
are two unique function calls to two different black-box
computer models that can be used for constrained and unconstrained
optimization, respectively. R is an open-source programming language,
and so none of the computer models within the package can ever truly be
a completely black-box function. However, the developers of the
CompModels package have done their best to obscure the analytical
forms of the mathematical functions underlying the computer models. For
example, at the first level of the code, a call to the bbox1()
function tells the user the following:
R> bbox1
function(x1,x2){
if(!is.numeric(x1) | !is.numeric(x2) | length(x1) != 1 | length(x2) !=1){
stop("Input is invalid.")
}else if(x1 < -1.5 | x1 > 2.5 | x2 < -3 | x2 > 3){
stop("Input is outside of the domain.")
}else{
ans <- .C("bbox1c",x1=x1,x2=x2,fx=0,c1x=0,c2x=0)
return(list(obj = ans$fx, con = c(ans$c1x,ans$c2x)))
}
}
The only discernible information that the user can glean from this
output is that the bbox1()
function has an input dimension of fx
, to minimize, and two constraint
functions, c1x
and c2x
, to satisfy. As we see from the .C()
command, the actual source code for the black-box function has been
written using the C programming language. The C programs are publicly
available, but the code within those programs has been heavily
obfuscated to the best of our abilities in order to obscure the source
code such that the computer models remain black-box functions. Moreover,
we believe that a good robust methodology developed for computer models
benefits from being applied to black-box functions and so any attempt to
decipher the black-box computer models is simply a disservice to the
statistician developing the methodology.
When developing the computer models in the package, we kept in mind that
the best computer model examples typically have roots in real
applications. When possible, we tried to develop computer models that
were either based on physics or that appeared in the literature with
real use cases. For example, one computer model, pressure()
, is based
on the real-world engineering problem of minimizing the cost associated
with constructing a pressure vessel (Figure 1). Given
the thickness of the shell (
Likewise, when possible, we sought out real-world problems where
solutions already existed that could be benchmarked to. For example, the
tension spring computer model, tension()
, is designed to minimize the
weight of a tension spring (Figure 2) subject to four
constraints on the shear stress, surge frequency, and deflection. The
three inputs to the computer model are for the wire diameter (
The tension spring problem has been solved many times in the literature, and Table 1 summarizes some of the best solutions.
Optimal Inputs | ||||
Source | Best Solution | |||
(Coello 2000) | 0.051480 | 0.351661 | 11.632201 | 0.012704 |
(He and L. Wang 2007) | 0.051728 | 0.357644 | 11.244543 | 0.012675 |
(Gandomi, X. Yang, A. Alavi, and S. Talatahari 2013) | 0.051690 | 0.356730 | 11.288500 | 0.012670 |
(Mirjalili, S. Mirjalili, and A. Lewis 2014) | 0.051690 | 0.356737 | 11.288850 | 0.012666 |
(Lee and Z. Geem 2005) | 0.051154 | 0.349871 | 12.076432 | 0.012671 |
(Askarzadeh 2016) | 0.051689 | 0.356717 | 11.289012 | 0.012665 |
(Mirjalili, A. Gandomi, Z. Mirjalili, S. Saremi, H. Fairs, and S. Mirjalili 2017) | 0.051207 | 0.345215 | 12.004032 | 0.012676 |
(Li, F. Shuang, P. Zhao, and C. Le 2019) | 0.051618 | 0.355004 | 11.390144 | 0.012665 |
We stress the need for benchmarking in our examples because we believe that benchmarking also helps with allowing for good computer model methodology to be developed. In the computer modeling literature, one tends to see real-world optimization results that stand alone and cannot be compared against or even replicated because practitioners do not have access to the same computer models as others. Being able to benchmark one’s results to others helps discern how well a given optimization method performs and allows for useful internal feedback when developing a method. Thus, a key reason we have developed the CompModels package is so that equitable access to computer models for benchmarking exists. Similarly, a problem with real-world computer models is that they can change over time, and often older versions will be phased out, unsupported, or disappear entirely. For example, the optimization results for the MODFLOW-96 computer model (McDonald and A. Harbaugh 1996) from (Pourmohamad and H. K. H. Lee 2016) was benchmarked to the work in (Lindberg and H. K. H. Lee 2015). However, this computer model is no longer supported by its developers, and so future benchmarking may become infeasible. Thus, the CompModels package also stands as a repository of computer models that should be available to all users for the foreseeable future. Lastly, computer models can often be platform and operating system specific, which ultimately limits the number of potential users of the computer model. Given that R packages, for the most part, tend to be immune to this problem, the CompModels package would be available to as wide of an audience as possible, again providing equitable access to computer models.
The remainder of the paper is organized as follows. Section 2 gives a brief introduction to Bayesian optimization and expected feasible improvement so that the computer models within the CompModels package can be demonstrated. Section 3 illustrates practical applications of package use for optimization, and Section 4 concludes with a discussion.
Tracing its roots as far back as to (Mockus, V. Tiesis, and A. Zilinskas 1978), Bayesian
optimization (BO) is a sequential design strategy for efficiently
optimizing black-box functions in a few steps that does not require
gradient information (Brochu, V. M. Cora, and N. de Freitas 2010). More specifically, BO seeks to
solve the minimization problem
In order to solve the problems in ((2)) and ((3)), an
acquisition function must be chosen for efficiently guiding the search.
Perhaps, one of the most popular acquisition functions for unconstrained
Bayesian optimization is that of expected improvement (EI)
(Jones, M. Schonlau, and W. J. Welch 1998). Originally introduced in the computer
modeling literature, (Jones, M. Schonlau, and W. J. Welch 1998) defined the
improvement statistic at a proposed input
Extending EI to the constrained optimization case, (Schonlau, W. J. Welch, and D. Jones 1998)
defined expected feasible improvement (EFI) as
We illustrate the use and functionality of the computer models in the
CompModels package by solving two constrained optimization problems
using the EFI method outlined in Section 2. We optimize the
tension spring computer model, tension()
, as well as the black-box 1
computer model, bbox1()
. In both cases, we perform Monte Carlo
experiments where we repeat the optimization routine a total of 30 times
to judge the robustness of the solutions. We take advantage of the
function optim.efi()
in the
laGP package (Gramacy 2016) for
running the EFI algorithm. A full list of the available computer models
in the CompModels package is given in the Appendix and is
generalizable to the proceeding examples.
The goal of the tension spring computer model is to minimize the weight
of the tension spring subject to four constraints on the shear stress,
surge frequency, and deflection. Here, the inputs to the tension spring
computer model are the wire diameter (
R> tension(x1 = 1, x2 = 1, x3 = 3)
$obj
[1] 5
$con
[1] 0.9999582 -45.8166667 -0.9995655 0.3333333
All of the computer model functions in the package will return a list
where the first element in the list is the value of the objective
function, and (in the case of constrained optimization) the second
element contains the values of the constraint functions. Here we see
that for a wire diameter of x1 = 1
, mean coil diameter of x2 = 1
,
and x3 = 3
active coils, that the weight of the tension spring is
five. However, the first and last constraint has not been satisfied
since those values of $con
are non-negative. Thus, the input is not a
feasible solution to the problem. The input of optim.efi()
in the laGP package. To be able to use the
optim.efi()
function, we need to first build a wrapper function (which
we call bbox
) for our tension()
function that conforms to the
specifications of the optim.efi()
function.
R> bbox <- function(X){
+ output = tension(X[1], X[2], X[3])
+ return(list(obj = output$obj, c = output$con))
}
Next, we need to create a matrix that encodes the domain of the computer model inputs.
R> B <- matrix(c(.05, .25, 2, 2, 1.3, 15), nrow=3)
We can implement the EFI algorithm by passing our wrapper function and
domain variable as arguments to the optim.efi()
function, and then by
checking the regions where the solution satisfies the constraints.
R> ans <- optim.efi(bbox, B, fhat = TRUE, start = 10, end = 300)
R> constraint <- ifelse(apply(ans$C, 1, max) > 0, "Not Met", "Met")
Here, we see that the optim.efi()
function started with a random input
of 10 data points and sequentially chose 290 more inputs for a total of
300 evaluations. The output of optim.efi()
is a large list storing all
steps of the EFI algorithm. We create the constraint
variable in order
to be able to find where the minimum feasible value exists.
R> min(ans$obj[constraint == "Met"])
[1] 0.0112376
R> ans$X[ans$obj == min(ans$obj[constraint == "Met"])]
[1] 0.05345441 0.45253754 6.69064005
Here, we see that the best feasible value found by the EFI algorithm is
at a weight of 0.0112376, which occurs at an input of
R> S <- 30
R> results <- rep(NA, S)
R> for(i in 1:S){
+ ans <- optim.efi(bbox, B, fhat = TRUE, start = 10, end = 300)
+ constraint <- ifelse(apply(ans$C, 1, max) > 0, "Not Met", "Met")
+ results[i] <- min(ans$obj[constraint == "Met"])
+}
R> summary(results)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.01081 0.01255 0.01302 0.01325 0.01386 0.01859
From the summary of the results, we see that over the 30 Monte Carlo experiments that the EFI algorithm was not able to reliably find as good of a solution over the 300 computer model evaluations. The mean value over the 30 runs was 0.01325, which was much higher than the best solutions presented in Table 1. However, we do see from the summary that the EFI algorithm was able to find at least one more better solution, as compared to the literature, at a spring weight of 0.01081.
Recalling Section 1, the bbox1()
computer model has an
input dimension of fx
, to minimize, and two constraint functions, c1x
and
c2x
, to satisfy. We can, once again, use the optim.efi()
function to
perform the EFI algorithm by creating an appropriate wrapper function
and domain variable.
R> bbox <- function(X){
+ output = bbox1(X[1], X[2])
+ return(list(obj = output$obj, c = output$con))
+}
R> B <- matrix(c(-1.5, -3, 2.5, 3), nrow = 2)
We initialize the optim.efi()
function with an input data set of 10
points and continue to sequentially evaluate the bbox1()
function for
a total of 100 input points.
R> ans <- optim.efi(bbox, B, fhat = TRUE, start = 10, end = 100)
R> constraint <- ifelse(apply(ans$C, 1, max) > 0, "Not Met", "Met")
Checking the EFI algorithm results in the areas where the constraint
functions were satisfied, we obtain a best feasible minimum objective
function value of -4.61008, which occurs at
R> min(ans$obj[constraint == "Met"])
[1] -4.610088
R> xbest <- ans$X[ans$obj == min(ans$obj[constraint == "Met"])]
R> xbest
[1] 0.204649 2.072964
Now, since the bbox1()
function is a black-box computer model, we do
not have any analytical way of checking whether or not our solution to
the optimization problem is a good one. However, the functions in the
CompModels package were not developed with the intent of forcing them
to be computationally expensive if they need not be. Thus, with an input
dimension of bbox1()
function
on a very dense grid to understand what the potential surface of the
objective and constraint functions look like. Doing so does not
guarantee us analytically that our solution is a good one, but we will
be able to tell visually whether or not our solution is a good one.
Plotting the objective and constraint surfaces, we obtain the following
(Figure 3).
R> n <- 200
R> x1 <- seq(-1.5, 2.5, len = n)
R> x2 <- seq(-3, 3, len = n)
R> x <- expand.grid(x1, x2)
R> obj <- rep(NA, nrow(x))
R> con <- matrix(NA, nrow = nrow(x), ncol = 2)
R> for(i in 1:nrow(x)){
+ temp <- bbox1(x[i,1], x[i,2])
+ obj[i] <- temp$obj
+ con[i,] <- temp$con
+}
R> y <- obj
R> y[con[,1] > 0 | con[,2] > 0] <- NA
R> z <- obj
R> z[!(con[,1] > 0 | con[,2] > 0)] <- NA
R> par(ps=15)
R> plot(0, 0, type = "n", xlim = c(-1.5, 2.5), ylim = c(-3, 3),
+ xlab = expression(x[1]), ylab = expression(x[2]), main = "Black-box Function")
R> c1 <- matrix(con[,1], ncol = n)
R> contour(x1, x2, c1, nlevels = 1, levels = 0, drawlabels = FALSE, add = TRUE,
+ lwd = 2)
R> c2 <- matrix(con[,2], ncol = n)
R> contour(x1, x2, c2, nlevels = 1, levels = 0, drawlabels = FALSE, add = TRUE,
+ lwd = 2, lty = 2)
R> contour(x1, x2, matrix(y, ncol = n), nlevels = 10, add = TRUE, col = "forestgreen")
R> contour(x1, x2, matrix(z, ncol = n), nlevels = 20, add = TRUE, col = 2, lty = 2)
R> points(xbest[1], xbest[2], pch = 21, bg = "deepskyblue")
By plotting the objective function surface, and the constraint
functions, we see that the space where the constraints are satisfied are
two disconnected regions where the feasible region with
R> summary(results)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.681 -4.661 -4.645 -4.645 -4.632 -4.602
From the summary of the results, we see that the variation in the results show up in the hundredth decimal point and beyond, which we regard as representing a very robust solution.
The primary goal of the package is to provide users with a source of computer model test functions that are reproducible, shareable, and that can ultimately be used for benchmarking of Bayesian optimization methods. The package will greatly benefit those who do not have access, or connections, to real-world computer models. In time, it is our hope that the package will come to be viewed as a suite of real computer models rather than solely as pseudo ones. Likewise, the CompModels package is not a static package in that we envision it to be a living repository, and so more computer model functions will be expected to be added over time. The success of any R package ultimately comes from the feedback received from its users. We greatly encourage all interested users of the package to please contact the developers in order to provide any insights or examples for new computer models to be added.
Table 2 provides a summary of the current computer models that are available in the CompModels package. The package is a mix of real-world physics problems, known mathematical functions, and black-box functions, as well as a mix of constrained or unconstrained optimization problems.
Function | Input Dimension | Optimization Type | No. of Constraints |
---|---|---|---|
bbox1() |
2 | Constrained | 2 |
bbox2() |
2 | Unconstrained | – |
bbox3() |
2 | Unconstrained | – |
bbox4() |
2 | Constrained | 1 |
bbox5() |
3 | Unconstrained | – |
bbox6() |
1 | Constrained | 2 |
bbox7() |
8 | Constrained | 2 |
gram() |
2 | Constrained | 2 |
mtp() |
2 | Constrained | 2 |
pressure() |
4 | Constrained | 4 |
sprinkler() |
8 | Unconstrained | – |
tension() |
3 | Constrained | 4 |
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Pourmohamad, "CompModels: A Suite of Computer Model Test Functions for Bayesian Optimization", The R Journal, 2021
BibTeX citation
@article{RJ-2021-076, author = {Pourmohamad, Tony}, title = {CompModels: A Suite of Computer Model Test Functions for Bayesian Optimization}, journal = {The R Journal}, year = {2021}, note = {https://rjournal.github.io/}, volume = {13}, issue = {2}, issn = {2073-4859}, pages = {441-449} }