crit_qEHI: Batch Expected Hypervolume Improvement with m objectives

In GPareto: Gaussian Processes for Pareto Front Estimation and Optimization

Batch Expected Hypervolume Improvement with m objectives

Description

Parallel Multi-objective Expected Hypervolume Improvement with respect to the current Pareto front, via Sample Average Approximation (SAA).

Usage

crit_qEHI(
  x,
  model,
  paretoFront = NULL,
  critcontrol = list(nb.samp = 50, seed = 42),
  type = "UK"
)

Arguments

x	a matrix representing the inputs (one per row) for which one wishes to calculate EHI,
model	list of objects of class km, one for each objective functions,
paretoFront	(optional) matrix corresponding to the Pareto front of size [n.pareto x n.obj], or any reference set of observations,
critcontrol	optional list with arguments: nb.samp number of random samples from the posterior distribution (with more than two objectives), default to 50, increasing gives more reliable results at the cost of longer computation time; seed seed used for the random samples (with more than two objectives); refPoint reference point for Hypervolume Expected Improvement; extendper if no reference point refPoint is provided, for each objective it is fixed to the maximum over the Pareto front plus extendper times the range, Default value to 0.2, corresponding to 1.1 for a scaled objective with a Pareto front in [0,1]^n.obj.
type	"SK" or "UK" (by default), depending whether uncertainty related to trend estimation has to be taken into account.

Details

The batch EHI is computed by simulated nb.samp conditional simulation of the objectives on the candidate points, before averaging the corresponding hypervolume improvement of each set of m simulations.

Value

The Expected Hypervolume Improvement at x.

References

J. D. Svenson (2011), Computer Experiments: Multiobjective Optimization and Sensitivity Analysis, Ohio State University, PhD thesis.

See Also

EI from package DiceOptim, crit_EMI, crit_SUR, crit_SMS.

Examples

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

n_var <- 2 
f_name <- "P1" 
n.grid <- 26
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
test.grid <- as.matrix(test.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))
Front_Pareto <- t(nondominated_points(t(response.grid)))
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])

refPoint <- c(300, 0)
EHI_grid <- crit_EHI(x = test.grid, model = list(mf1, mf2),
         critcontrol = list(refPoint = refPoint))

## Create potential batches
q <- 3
ncb <- 500
Xbcands <- array(NA, dim = c(ncb, q, n_var))
for(i in 1:ncb) Xbcands[i,,] <- test.grid[sample(1:nrow(test.grid), q, prob = pmax(0,EHI_grid)),]

qEHI_grid <- apply(Xbcands, 1, crit_qEHI, model = list(mf1, mf2),
  critcontrol = list(refPoint = refPoint))
  
Xq <- Xbcands[which.max(qEHI_grid),,]

## For further optimization (gradient may not be reliable)
# sol <- optim(as.vector(Xq), function(x, ...) crit_qEHI(matrix(x, q), ...),
#   model = list(mf1, mf2), lower = c(0,0), upper = c(1,1), method = "L-BFGS-B",
#   control = list(fnscale = -1), critcontrol = list(refPoint = refPoint, nb.samp = 10000))
# sol <- psoptim(as.vector(Xq), function(x, ...) crit_qEHI(matrix(x, q), ...),
#  model = list(mf1, mf2), lower = c(0,0), upper = c(1,1), 
#  critcontrol = list(refPoint = refPoint, nb.samp = 100),
#  control = list(fnscale = -1))

# Plot EHI surface and selected designs for parallel evaluation
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
               matrix(EHI_grid, n.grid), 
               main = "Expected Hypervolume Improvement surface and best 3-EHI designs",
               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")
                            points(Xq, pch = 20, col = 2)
#                            points(matrix(sol$par, q), col = 4)
                            }
              )

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