crit_SUR: Analytical expression of the SUR criterion for two or three...

In GPareto: Gaussian Processes for Pareto Front Estimation and Optimization

Analytical expression of the SUR criterion for two or three objectives.

Description

Computes the SUR criterion (Expected Excursion Volume Reduction) at point x for 2 or 3 objectives. To avoid numerical instabilities, the new point is penalized if it is too close to an existing observation.

Usage

crit_SUR(x, model, paretoFront = NULL, critcontrol = NULL, type = "UK")

Arguments

x	a vector representing the input for which one wishes to calculate the criterion,
model	a list of objects of class km (one for each objective),
paretoFront	(optional) matrix corresponding to the Pareto front of size [n.pareto x n.obj], or any reference set of observations,
critcontrol	list with two possible options. A) One can use the four following arguments: integration.points, matrix of integration points of size [n.integ.pts x d]; integration.weights, vector of integration weights of length n.integ.pts; mn.X and sn.X, matrices of kriging means and sd, each of size [n.obj x n.integ.pts]; precalc.data, list of precalculated data (based on kriging models at integration points) for faster computation. B) Alternatively, one can define arguments passed to integration_design_optim: SURcontrol (optional), lower, upper, min.prob (optional). This is slower since arguments of A), used in the function, are then recomputed each time (note that this is not the case when called from GParetoptim and crit_optimizer). Options for the checkPredict function: threshold (1e-4) and distance (covdist) are used to avoid numerical issues occuring when adding points too close to the existing ones.
type	"SK" or "UK" (default), depending whether uncertainty related to trend estimation has to be taken into account.

Value

Value of the criterion.

References

V. Picheny (2014), Multiobjective optimization using Gaussian process emulators via stepwise uncertainty reduction, Statistics and Computing.

See Also

crit_EHI, crit_SMS, crit_EMI.

Examples

#---------------------------------------------------------------------------
# crit_SUR surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
set.seed(25468)
library(DiceDesign)
library(KrigInv)

n_var <- 2 
n.obj <- 2 
f_name <- "P1" 
n.grid <- 14
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
n_appr <- 15 
design.grid <- round(maximinESE_LHS(lhsDesign(n_appr, n_var, seed = 42)$design)$design, 1)
response.grid <- t(apply(design.grid, 1, f_name))
paretoFront <- t(nondominated_points(t(response.grid)))
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])

model <- list(mf1, mf2)


integration.param <- integration_design_optim(lower = c(0, 0), upper = c(1, 1), model = model)
integration.points <- as.matrix(integration.param$integration.points)
integration.weights <- integration.param$integration.weights

precalc.data <- list()
mn.X <- sn.X <- matrix(0, nrow = n.obj, ncol = nrow(integration.points))

for (i in 1:n.obj){
  p.tst.all <- predict(model[[i]], newdata = integration.points, type = "UK", checkNames = FALSE)
  mn.X[i,] <- p.tst.all$mean
  sn.X[i,]   <- p.tst.all$sd
  precalc.data[[i]] <- precomputeUpdateData(model[[i]], integration.points)
}

critcontrol <- list(integration.points = integration.points,
                    integration.weights = integration.weights,
                    mn.X = mn.X, sn.X = sn.X, precalc.data = precalc.data)
## Alternatively: critcontrol <- list(lower = rep(0, n_var), upper = rep(1,n_var))
                
EEV_grid <- apply(test.grid, 1, crit_SUR, model = model, paretoFront = paretoFront,
                  critcontrol = critcontrol)
                  

filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid),
               matrix(pmax(0,EEV_grid), nrow = n.grid), main = "EEV criterion",
               xlab = expression(x[1]), ylab = expression(x[2]), color = terrain.colors,
               plot.axes = {axis(1); axis(2);
                            points(design.grid[,1], design.grid[,2], pch = 21, bg = "white")
                            }
              )         

GPareto documentation built on May 29, 2024, 3:55 a.m.