Package GPareto

Description

Multi-objective optimization and quantification of uncertainty on Pareto fronts, using Gaussian process models.

Details

Important functions:
GParetoptim
easyGParetoptim
crit_optimizer
plotGPareto
CPF

Note

Part of this work has been conducted within the frame of the ReDice Consortium, gathering industrial (CEA, EDF, IFPEN, IRSN, Renault) and academic (Ecole des Mines de Saint-Etienne, INRIA, and the University of Bern) partners around advanced methods for Computer Experiments. (http://redice.emse.fr/).

The authors would like to thank Yves Deville for his precious advices in R programming and packaging, as well as Olivier Roustant and David Ginsbourger for testing and suggestions of improvements for this package. We would also like to thank Tobias Wagner for providing his Matlab codes for the SMS-EGO strategy.

Author(s)

Mickael Binois, Victor Picheny

References

M. Binois, D. Ginsbourger and O. Roustant (2015), Quantifying Uncertainty on Pareto Fronts with Gaussian process conditional simulations, European Journal of Operational Research, 243(2), 386-394.

O. Roustant, D. Ginsbourger and Yves Deville (2012), DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization, Journal of Statistical Software, 51(1), 1-55, doi:10.18637/jss.v051.i01.

M. T. Emmerich, A. H. Deutz, J. W. Klinkenberg (2011), Hypervolume-based expected improvement: Monotonicity properties and exact computation, Evolutionary Computation (CEC), 2147-2154.

V. Picheny (2015), Multiobjective optimization using Gaussian process emulators via stepwise uncertainty reduction, Statistics and Computing, 25(6), 1265-1280.

T. Wagner, M. Emmerich, A. Deutz, W. Ponweiser (2010), On expected-improvement criteria for model-based multi-objective optimization, Parallel Problem Solving from Nature, 718-727, Springer, Berlin.

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

C. Chevalier (2013), Fast uncertainty reduction strategies relying on Gaussian process models, University of Bern, PhD thesis.

M. Binois, V. Picheny (2019), GPareto: An R Package for Gaussian-Process-Based Multi-Objective Optimization and Analysis, Journal of Statistical Software, 89(8), 1-30, doi:10.18637/jss.v089.i08.

See Also

DiceKriging-package, DiceOptim-package

Examples

## Not run: 
#------------------------------------------------------------
# Example 1 : Surrogate-based multi-objective Optimization with postprocessing
#------------------------------------------------------------
set.seed(25468)

d <- 2 
fname <- P2

plotParetoGrid(P2) # For comparison

# Optimization
budget <- 25 
lower <- rep(0, d) 
upper <- rep(1, d)
     
omEGO <- easyGParetoptim(fn = fname, budget = budget, lower = lower, upper = upper)

# Postprocessing
plotGPareto(omEGO, add= FALSE, UQ_PF = TRUE, UQ_PS = TRUE, UQ_dens = TRUE)


## End(Not run)
#------------------------------------------------------------
# Example 2 : Surrogate-based multi-objective Optimization including a cheap function
#------------------------------------------------------------
set.seed(42)
library(DiceDesign)

d <- 2 

fname <- P1
n.grid <- 19
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
nappr <- 15 
design.grid <- maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))

mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])
model <- list(mf1, mf2)

nsteps <- 1 
lower <- rep(0, d)
upper <- rep(1, d)

# Optimization with fastfun: hypervolume with discrete search

optimcontrol <- list(method = "discrete", candidate.points = test.grid)
omEGO2 <- GParetoptim(model = model, fn = fname, cheapfn = branin, crit = "SMS",
                      nsteps = nsteps, lower = lower, upper = upper,
                      optimcontrol = optimcontrol)
print(omEGO2$par)
print(omEGO2$values)

## Not run:  
plotGPareto(omEGO2)

#------------------------------------------------------------
# Example 3 : Surrogate-based multi-objective Optimization (4 objectives)
#------------------------------------------------------------
set.seed(42)
library(DiceDesign)

d <- 5 

fname <- DTLZ3
nappr <- 25
design.grid <- maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname, nobj = 4))
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])
mf3 <- km(~., design = design.grid, response = response.grid[,3])
mf4 <- km(~., design = design.grid, response = response.grid[,4])

# Optimization
nsteps <- 5 
lower <- rep(0, d) 
upper <- rep(1, d)     
omEGO3 <- GParetoptim(model = list(mf1, mf2, mf3, mf4), fn = fname, crit = "EMI",
                      nsteps = nsteps, lower = lower, upper = upper, nobj = 4)
print(omEGO3$par)
print(omEGO3$values)
plotGPareto(omEGO3)

#------------------------------------------------------------
# Example 4 : quantification of uncertainty on Pareto front
#------------------------------------------------------------
library(DiceDesign)
set.seed(42)

nvar <- 2

# Test function P1
fname <- "P1"

# Initial design
nappr <- 10
design.grid <- maximinESE_LHS(lhsDesign(nappr, nvar, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))

PF <- t(nondominated_points(t(response.grid)))

# kriging models : matern5_2 covariance structure, linear trend, no nugget effect
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])

# Conditional simulations generation with random sampling points 
nsim <- 100 # increase for better results
npointssim <- 1000 # increase for better results
Simu_f1 <- matrix(0, nrow = nsim, ncol = npointssim)
Simu_f2 <- matrix(0, nrow = nsim, ncol = npointssim)
design.sim <- array(0, dim = c(npointssim, nvar, nsim))

for(i in 1:nsim){
  design.sim[,,i] <- matrix(runif(nvar*npointssim), nrow = npointssim, ncol = nvar)
  Simu_f1[i,] <- simulate(mf1, nsim = 1, newdata = design.sim[,,i], cond = TRUE,
                         checkNames = FALSE, nugget.sim = 10^-8)
  Simu_f2[i,] <- simulate(mf2, nsim = 1, newdata = design.sim[,,i], cond = TRUE,
                         checkNames = FALSE, nugget.sim = 10^-8)
}

# Computation of the attainment function and Vorob'ev Expectation
CPF1 <- CPF(Simu_f1, Simu_f2, response.grid, paretoFront = PF)

summary(CPF1)

plot(CPF1)

# Display of the symmetric deviation function
plotSymDevFun(CPF1)

## End(Not run)

Conditional Pareto Front simulations

Description

Compute (on a regular grid) the empirical attainment function from conditional simulations of Gaussian processes corresponding to two objectives. This is used to estimate the Vorob'ev expectation of the attained set and the Vorob'ev deviation.

Usage

CPF(
  fun1sims,
  fun2sims,
  response,
  paretoFront = NULL,
  f1lim = NULL,
  f2lim = NULL,
  refPoint = NULL,
  n.grid = 100,
  compute.VorobExp = TRUE,
  compute.VorobDev = TRUE
)

Arguments

Details

Works with two objectives. The user can provide locations of grid lines for computation of the attainement function with vectors f1lim and f2lim, in the form of regularly spaced points. It is possible to provide only refPoint as a reference for hypervolume computations. When missing, values are determined from the axis-wise extrema of the simulations.

Value

A list which is given the S3 class "CPF".

• x, y: locations of grid lines at which the values of the attainment are computed,

• values: numeric matrix containing the values of the attainment on the grid,

• PF: matrix corresponding to the Pareto front of the observations,

• responses: matrix containing the value of the two objective functions, one objective per column,

• fun1sims, fun2sims: conditional simulations of the first/second output,

• VE: Vorob'ev expectation, computed if compute.VorobExp = TRUE (default),

• beta_star: Vorob'ev threshold, computed if compute.VorobExp = TRUE (default),

• VD: Vorov'ev deviation, computed if compute.VorobDev = TRUE (default),

References

M. Binois, D. Ginsbourger and O. Roustant (2015), Quantifying Uncertainty on Pareto Fronts with Gaussian process conditional simulations, European Journal of Operational Research, 243(2), 386-394.

C. Chevalier (2013), Fast uncertainty reduction strategies relying on Gaussian process models, University of Bern, PhD thesis.

I. Molchanov (2005), Theory of random sets, Springer.

See Also

Methods coef, summary and plot can be used to get the coefficients from a CPF object, to obtain a summary or to display the attainment function (with the Vorob'ev expectation if compute.VorobExp is TRUE).

Examples

library(DiceDesign)
set.seed(42)

nvar <- 2

fname <- "P1" # Test function

# Initial design
nappr <- 10
design.grid <- maximinESE_LHS(lhsDesign(nappr, nvar, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))

# kriging models: matern5_2 covariance structure, linear trend, no nugget effect
mf1 <- km(~., design = design.grid, response = response.grid[,1])
mf2 <- km(~., design = design.grid, response = response.grid[,2])

# Conditional simulations generation with random sampling points 
nsim <- 40
npointssim <- 150 # increase for better results
Simu_f1 <- matrix(0, nrow = nsim, ncol = npointssim)
Simu_f2 <- matrix(0, nrow = nsim, ncol = npointssim)
design.sim <- array(0, dim = c(npointssim, nvar, nsim))

for(i in 1:nsim){
  design.sim[,,i] <- matrix(runif(nvar*npointssim), nrow = npointssim, ncol = nvar)
  Simu_f1[i,] <- simulate(mf1, nsim = 1, newdata = design.sim[,,i], cond = TRUE,
                          checkNames = FALSE, nugget.sim = 10^-8)
  Simu_f2[i,] <- simulate(mf2, nsim = 1, newdata = design.sim[,,i], cond = TRUE, 
                          checkNames = FALSE, nugget.sim = 10^-8)
}

# Attainment and Voreb'ev expectation + deviation estimation
CPF1 <- CPF(Simu_f1, Simu_f2, response.grid)

# Details about the Vorob'ev threshold  and Vorob'ev deviation
summary(CPF1)

# Graphics
plot(CPF1)

Sequential multi-objective Expected Improvement maximization and model re-estimation, with a number of iterations fixed in advance by the user

Description

Executes nsteps iterations of multi-objective EGO methods to objects of class km. At each step, kriging models are re-estimated (including covariance parameters re-estimation) based on the initial design points plus the points visited during all previous iterations; then a new point is obtained by maximizing one of the four multi-objective Expected Improvement criteria available. Handles noiseless and noisy objective functions.

Usage

GParetoptim(
  model,
  fn,
  ...,
  cheapfn = NULL,
  crit = "SMS",
  nsteps,
  lower,
  upper,
  type = "UK",
  cov.reestim = TRUE,
  critcontrol = NULL,
  noise.var = NULL,
  reinterpolation = NULL,
  optimcontrol = list(method = "genoud", trace = 1),
  ncores = 1
)

Arguments

Details

Extension of the function EGO.nsteps for multi-objective optimization.
Available infill criteria with crit are:

• Expected Hypervolume Improvement (EHI) crit_EHI,

• SMS criterion (SMS) crit_SMS,

• Expected Maximin Improvement (EMI) crit_EMI,

• Stepwise Uncertainty Reduction of the excursion volume (SUR) crit_SUR.

Depending on the selected criterion, parameters such as reference point for SMS and EHI or arguments for integration_design_optim with SUR can be given with critcontrol. Also options for checkPredict are available. More precisions are given in the corresponding help pages.

The reinterpolation=TRUE setting can be used to handle noisy objective functions. It works with all criteria and is the recommended option. If reinterpolation=FALSE and noise.var!=NULL, the criteria are used based on a "denoised" Pareto front.

If noise.var="given_by_fn", fn must return a list of two vectors, the first being the objective functions and the second the corresponding noise variances (see examples).

Display of results and various post-processings are available with plotGPareto.

Value

A list with components:

• par: a data frame representing the additional points visited during the algorithm,

• values: a data frame representing the response values at the points given in par,

• nsteps: an integer representing the desired number of iterations (given in argument),

• lastmodel: a list of objects of class km corresponding to the last kriging models fitted.

• observations.denoised: if noise.var!=NULL, a matrix representing the mean values of the km models at observation points. If a problem occurs during either model updates or criterion maximization, the last working model and corresponding values are returned.

References

M. T. Emmerich, A. H. Deutz, J. W. Klinkenberg (2011), Hypervolume-based expected improvement: Monotonicity properties and exact computation, Evolutionary Computation (CEC), 2147-2154.

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

T. Wagner, M. Emmerich, A. Deutz, W. Ponweiser (2010), On expected-improvement criteria for model-based multi-objective optimization. Parallel Problem Solving from Nature, 718-727, Springer, Berlin.

J. D. Svenson (2011), Computer Experiments: Multiobjective Optimization and Sensitivity Analysis, Ohio State university, PhD thesis. V. Picheny and D. Ginsbourger (2013), Noisy kriging-based optimization methods: A unified implementation within the DiceOptim package, Computational Statistics & Data Analysis, 71: 1035-1053.

Examples

set.seed(25468)
library(DiceDesign)

################################################
# NOISELESS PROBLEMS
################################################
d <- 2 
fname <- ZDT3
n.grid <- 21
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
nappr <- 15 
design.grid <- maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))
Front_Pareto <- t(nondominated_points(t(response.grid)))

mf1 <- km(~1, design = design.grid, response = response.grid[, 1], lower=c(.1,.1))
mf2 <- km(~., design = design.grid, response = response.grid[, 2], lower=c(.1,.1))
model <- list(mf1, mf2)

nsteps <- 2
lower <- rep(0, d)
upper <- rep(1, d)

# Optimization 1: EHI with pso
optimcontrol <- list(method = "pso", maxit = 20)
critcontrol <- list(refPoint = c(1, 10))
omEGO1 <- GParetoptim(model = model, fn = fname, crit = "EHI", nsteps = nsteps,
                     lower = lower, upper = upper, critcontrol = critcontrol,
                     optimcontrol = optimcontrol)
print(omEGO1$par)
print(omEGO1$values)

## Not run: 
nsteps <- 10
# Optimization 2: SMS with discrete search
optimcontrol <- list(method = "discrete", candidate.points = test.grid)
critcontrol <- list(refPoint = c(1, 10))
omEGO2 <- GParetoptim(model = model, fn = fname, crit = "SMS", nsteps = nsteps,
                     lower = lower, upper = upper, critcontrol = critcontrol,
                     optimcontrol = optimcontrol)
print(omEGO2$par)
print(omEGO2$values)

# Optimization 3: SUR with genoud
optimcontrol <- list(method = "genoud", pop.size = 20, max.generations = 10)
critcontrol <- list(distrib = "SUR", n.points = 100)
omEGO3 <- GParetoptim(model = model, fn = fname, crit = "SUR", nsteps = nsteps,
                     lower = lower, upper = upper, critcontrol = critcontrol,
                     optimcontrol = optimcontrol)
print(omEGO3$par)
print(omEGO3$values)

# Optimization 4: EMI with pso
optimcontrol <- list(method = "pso", maxit = 20)
critcontrol <- list(nbsamp = 200)
omEGO4 <- GParetoptim(model = model, fn = fname, crit = "EMI", nsteps = nsteps,
                     lower = lower, upper = upper, optimcontrol = optimcontrol)
print(omEGO4$par)
print(omEGO4$values)

# graphics
sol.grid <- apply(expand.grid(seq(0, 1, length.out = 100),
                              seq(0, 1, length.out = 100)), 1, fname)
plot(t(sol.grid), pch = 20, col = rgb(0, 0, 0, 0.05), xlim = c(0, 1),
     ylim = c(-2, 10), xlab = expression(f[1]), ylab = expression(f[2]))
plotGPareto(res = omEGO1, add = TRUE,
            control = list(pch = 20, col = "blue", PF.pch = 17,
                           PF.points.col = "blue", PF.line.col = "blue"))
text(omEGO1$values[,1], omEGO1$values[,2], labels = 1:nsteps, pos = 3, col = "blue")
plotGPareto(res = omEGO2, add = TRUE,
            control = list(pch = 20, col = "green", PF.pch = 17,
                           PF.points.col = "green", PF.line.col = "green"))
text(omEGO2$values[,1], omEGO2$values[,2], labels = 1:nsteps, pos = 3, col = "green")
plotGPareto(res = omEGO3, add = TRUE,
            control = list(pch = 20, col = "red", PF.pch = 17,
                           PF.points.col = "red", PF.line.col = "red"))
text(omEGO3$values[,1], omEGO3$values[,2], labels = 1:nsteps, pos = 3, col = "red") 
plotGPareto(res = omEGO4, add = TRUE,
            control = list(pch = 20, col = "orange", PF.pch = 17,
                           PF.points.col = "orange", PF.line.col = "orange"))
text(omEGO4$values[,1], omEGO4$values[,2], labels = 1:nsteps, pos = 3, col = "orange")
points(response.grid[,1], response.grid[,2], col = "black", pch = 20)
legend("topright", c("EHI", "SMS", "SUR", "EMI"), col = c("blue", "green", "red", "orange"),
 pch = rep(17,4))
 
 
# Post-processing
plotGPareto(res = omEGO1, UQ_PF = TRUE, UQ_PS = TRUE, UQ_dens = TRUE)

################################################
# NOISY PROBLEMS
################################################
set.seed(25468)
library(DiceDesign)
d <- 2 
nsteps <- 3
lower <- rep(0, d)
upper <- rep(1, d)
optimcontrol <- list(method = "pso", maxit = 20)
critcontrol <- list(refPoint = c(1, 10))

n.grid <- 21
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
n.init <- 30
design <- maximinESE_LHS(lhsDesign(n.init, d, seed = 42)$design)$design

fit.models <- function(u) km(~., design = design, response = response[, u],
                             noise.var=design.noise.var[,u])

# Test 1: EHI, constant noise.var
noise.var <- c(0.1, 0.2)
funnoise1 <- function(x) {ZDT3(x) + sqrt(noise.var)*rnorm(n=d)}
response <- t(apply(design, 1, funnoise1))
design.noise.var <- matrix(rep(noise.var, n.init), ncol=d, byrow=TRUE)
model <- lapply(1:d, fit.models)

omEGO1 <- GParetoptim(model = model, fn = funnoise1, crit = "EHI", nsteps = nsteps,
                      lower = lower, upper = upper, critcontrol = critcontrol,
                      reinterpolation=TRUE, noise.var=noise.var, optimcontrol = optimcontrol)
plotGPareto(omEGO1)

# Test 2: EMI, noise.var given by fn
funnoise2 <- function(x) {list(ZDT3(x) + sqrt(0.05 + abs(0.1*x))*rnorm(n=d), 0.05 + abs(0.1*x))}
temp <- funnoise2(design)
response <- temp[[1]]
design.noise.var <- temp[[2]]
model <- lapply(1:d, fit.models)

omEGO2 <- GParetoptim(model = model, fn = funnoise2, crit = "EMI", nsteps = nsteps,
                      lower = lower, upper = upper, critcontrol = critcontrol,
                      reinterpolation=TRUE, noise.var="given_by_fn", optimcontrol = optimcontrol)
plotGPareto(omEGO2)

# Test 3: SMS, functional noise.var
funnoise3 <- function(x) {ZDT3(x) + sqrt(0.025 + abs(0.05*x))*rnorm(n=d)}
noise.var <- function(x) return(0.025 + abs(0.05*x))
response <- t(apply(design, 1, funnoise3))
design.noise.var <- t(apply(design, 1, noise.var))
model <- lapply(1:d, fit.models)

omEGO3 <- GParetoptim(model = model, fn = funnoise3, crit = "SMS", nsteps = nsteps,
                           lower = lower, upper = upper, critcontrol = critcontrol,
                           reinterpolation=TRUE, noise.var=noise.var, optimcontrol = optimcontrol)
plotGPareto(omEGO3)

# Test 4: SUR, fastfun, constant noise.var
noise.var <- 0.1
funnoise4 <- function(x) {ZDT3(x)[1] + sqrt(noise.var)*rnorm(1)}
cheapfn <- function(x) ZDT3(x)[2]
response <- apply(design, 1, funnoise4)
design.noise.var <- rep(noise.var, n.init)
model <- list(km(~., design = design, response = response, noise.var=design.noise.var))

omEGO4 <- GParetoptim(model = model, fn = funnoise4, cheapfn = cheapfn, crit = "SUR", 
                      nsteps = nsteps, lower = lower, upper = upper, critcontrol = critcontrol,
                      reinterpolation=TRUE, noise.var=noise.var, optimcontrol = optimcontrol)
 plotGPareto(omEGO4)
                            
 # Test 5: EMI, fastfun, noise.var given by fn
 funnoise5 <- function(x) {
   if (is.null(dim(x))) x <- matrix(x, nrow=1)
   list(apply(x, 1, ZDT3)[1,] + sqrt(abs(0.05*x[,1]))*rnorm(nrow(x)), abs(0.05*x[,1]))
 }
 
 cheapfn <- function(x) {
   if (is.null(dim(x))) x <- matrix(x, nrow=1)
   apply(x, 1, ZDT3)[2,]
 }
 
 temp <- funnoise5(design)
 response <- temp[[1]]
 design.noise.var <- temp[[2]]
 model <- list(km(~., design = design, response = response, noise.var=design.noise.var))
 
 omEGO5 <- GParetoptim(model = model, fn = funnoise5, cheapfn = cheapfn, crit = "EMI", 
                       nsteps = nsteps, lower = lower, upper = upper, critcontrol = critcontrol,
                       reinterpolation=TRUE, noise.var="given_by_fn", optimcontrol = optimcontrol)
 plotGPareto(omEGO5)               
 
 # Test 6: EHI, fastfun, functional noise.var
 noise.var <- 0.1
 funnoise6 <- function(x) {ZDT3(x)[1] + sqrt(abs(0.1*x[1]))*rnorm(1)}
 noise.var <- function(x) return(abs(0.1*x[1]))
 cheapfn <- function(x) ZDT3(x)[2]
 response <- apply(design, 1, funnoise6)
 design.noise.var <- t(apply(design, 1, noise.var))
 model <- list(km(~., design = design, response = response, noise.var=design.noise.var))
 
 omEGO6 <- GParetoptim(model = model, fn = funnoise6, cheapfn = cheapfn, crit = "EMI", 
                       nsteps = nsteps, lower = lower, upper = upper, critcontrol = critcontrol,
                       reinterpolation=TRUE, noise.var=noise.var, optimcontrol = optimcontrol)
 plotGPareto(omEGO6)    

## End(Not run)

Estimation of Pareto set density

Description

Estimation of the density of Pareto optimal points in the variable space.

Usage

ParetoSetDensity(
  model,
  lower,
  upper,
  CPS = NULL,
  nsim = 50,
  simpoints = 1000,
  ...
)

Arguments

Details

This function estimates the density of Pareto optimal points in the variable space given by the surrogate models. Based on conditional simulations of the objectives at simulation points, Conditional Pareto Set (CPS) simulations are obtained, out of which a density is fitted.

This function relies on the ks-package package for the kernel density estimation.

Value

An object of class kde accounting for the estimated density of Pareto optimal points.

Examples

## Not run:  
#---------------------------------------------------------------------------
# Example of estimation of the density of Pareto optimal points
#---------------------------------------------------------------------------
set.seed(42)
n_var <- 2 
fname <- P1
lower <- rep(0, n_var)
upper <- rep(1, n_var)

res1 <- easyGParetoptim(fn = fname, lower = lower, upper = upper, budget = 15, 
control=list(method = "EHI", inneroptim = "pso", maxit = 20))

estDens <- ParetoSetDensity(res1$model, lower = lower, upper = upper)

# graphics
par(mfrow = c(1,2))
plot(estDens, display = "persp", xlab = "X1", ylab = "X2")
plot(estDens, display = "filled.contour2", main = "Estimated density of Pareto optimal point")
points(res1$model[[1]]@X[,1], res1$model[[2]]@X[,2], col="blue")
points(estDens$x[, 1], estDens$x[, 2], pch = 20, col = rgb(0, 0, 0, 0.15))
par(mfrow = c(1,1))

## End(Not run)

Test functions of x

Description

Multi-objective test functions.

Usage

ZDT1(x)

ZDT2(x)

ZDT3(x)

ZDT4(x)

ZDT6(x)

P1(x)

P2(x)

MOP2(x)

MOP3(x)

DTLZ1(x, nobj = 3)

DTLZ2(x, nobj = 3)

DTLZ3(x, nobj = 3)

DTLZ7(x, nobj = 3)

OKA1(x)

Arguments

Details

These functions are coming from different benchmarks: the ZDT test problems from an article of E. Zitzler et al., P1 from the thesis of J. Parr and P2 from an article of Poloni et al. . MOP2 and MOP3 are from Van Veldhuizen and DTLZ functions are from Deb et al. .

Domains (sometimes rescaled to [0,1]):

• ZDT1-6: [0,1]^d

• P1, P2: [0,1]^2

• MOP2: [0,1]^d

• MOP3: [-3,3], tri-objective, 2 variables

• DTLZ1-3,7: [0,1]^d, m-objective problems, with at least d>m variables.

• OKA1: [0,1]^2, initially [6 sin(pi/12), 6 sin(pi/12) + 2pi cos(pi/12)] x [-2pi sin(pi/12), 6 cos(pi/12)], bi-objective

Value

Matrix of values corresponding to the objective functions, the number of colums is the number of objectives.

References

J. M. Parr (2012), Improvement Criteria for Constraint Handling and Multiobjective Optimization, University of Southampton, PhD thesis.

C. Poloni, A. Giurgevich, L. Onesti, V. Pediroda (2000), Hybridization of a multi-objective genetic algorithm, a neural network and a classical optimizer for a complex design problem in fluid dynamics, Computer Methods in Applied Mechanics and Engineering, 186(2), 403-420.

E. Zitzler, K. Deb, and L. Thiele (2000), Comparison of multiobjective evolutionary algorithms: Empirical results, Evol. Comput., 8(2), 173-195.

K. Deb, L. Thiele, M. Laumanns and E. Zitzler (2002), Scalable Test Problems for Evolutionary Multiobjective Optimization, IEEE Transactions on Evolutionary Computation, 6(2), 182-197.

D. A. Van Veldhuizen, G. B. Lamont (1999), Multiobjective evolutionary algorithm test suites, In Proceedings of the 1999 ACM symposium on Applied computing, 351-357.

T. Okabe, J. Yaochu, M. Olhofer, B. Sendhoff (2004), On test functions for evolutionary multi-objective optimization, International Conference on Parallel Problem Solving from Nature, Springer, Berlin, Heidelberg.

Examples

# ----------------------------------
# 2-objectives test problems
# ---------------------------------- 

plotParetoGrid("ZDT1", n.grid = 21)

plotParetoGrid("ZDT2", n.grid = 21)

plotParetoGrid("ZDT3", n.grid = 21)

plotParetoGrid("ZDT4", n.grid = 21)

plotParetoGrid("ZDT6", n.grid = 21)

plotParetoGrid("P1", n.grid = 21)

plotParetoGrid("P2", n.grid = 21)

plotParetoGrid("MOP2", n.grid = 21)

Prevention of numerical instability for a new observation

Description

Check that the new point is not too close to already known observations to avoid numerical issues. Closeness can be estimated with several distances.

Usage

checkPredict(x, model, threshold = 1e-04, distance = "euclidean", type = "UK")

Arguments

Details

If the distance between x and the closest observations in model is below threshold, x should not be evaluated to avoid numerical instabilities. The distance can simply be the Euclidean distance or the canonical distance associated with the kriging predictive covariance k:

d(x,y) = \sqrt{k(x,x) - 2k(x,y) + k(y,y)}.

The last solution is the ratio between the prediction variance at x and the variance of the process. none can be used, e.g., if points have been selected already.

Value

TRUE if the point should not be tested.

Expected Hypervolume Improvement with m objectives

Description

Multi-objective Expected Hypervolume Improvement with respect to the current Pareto front. With two objectives the analytical formula is used, while Sample Average Approximation (SAA) is used with more objectives. To avoid numerical instabilities, the new point is penalized if it is too close to an existing observation.

Usage

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

Arguments

Details

The computation of the analytical formula with two objectives is adapted from the Matlab source code by Michael Emmerich and Andre Deutz, LIACS, Leiden University, 2010 available here : http://liacs.leidenuniv.nl/~csmoda/code/HV_based_expected_improvement.zip.

Value

The Expected Hypervolume Improvement at x.

References

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

M. T. Emmerich, A. H. Deutz, J. W. Klinkenberg (2011), Hypervolume-based expected improvement: Monotonicity properties and exact computation, Evolutionary Computation (CEC), 2147-2154.

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))
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])

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

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",
               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")
                            }
              )

Expected Maximin Improvement with m objectives

Description

Expected Maximin Improvement with respect to the current Pareto front with Sample Average Approximation. The semi-analytical formula is used in the bi-objective scale if the Pareto front is in [-2,2]^2, for numerical stability reasons. To avoid numerical instabilities, the new point is penalized if it is too close to an existing observation.

Usage

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

Arguments

Details

It is recommanded to scale objectives, e.g. to [0,1]. If the Pareto front does not belong to [-2,2]^2, then SAA is used.

Value

The Expected Maximin Improvement at x.

References

J. D. Svenson & T. J. Santner (2010), Multiobjective Optimization of Expensive Black-Box Functions via Expected Maximin Improvement, Technical Report.

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

See Also

EI from package DiceOptim, crit_EHI, crit_SUR, crit_SMS.

Examples

#---------------------------------------------------------------------------
# Expected Maximin 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 <- 21
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))
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])

EMI_grid <- apply(test.grid, 1, crit_EMI, model = list(mf1, mf2), paretoFront = Front_Pareto,
                     critcontrol = list(nb_samp = 20))

filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
               matrix(EMI_grid, nrow = n.grid), main = "Expected Maximin Improvement",
               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")
                            }
              )

Analytical expression of the SMS-EGO criterion with m>1 objectives

Description

Computes a slightly modified infill Criterion of the SMS-EGO. To avoid numerical instabilities, an additional penalty is added to the new point if it is too close to an existing observation.

Usage

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

Arguments

Value

Value of the criterion.

References

W. Ponweiser, T. Wagner, D. Biermann, M. Vincze (2008), Multiobjective Optimization on a Limited Budget of Evaluations Using Model-Assisted S-Metric Selection, Parallel Problem Solving from Nature, pp. 784-794. Springer, Berlin.

T. Wagner, M. Emmerich, A. Deutz, W. Ponweiser (2010), On expected-improvement criteria for model-based multi-objective optimization. Parallel Problem Solving from Nature, pp. 718-727. Springer, Berlin.

See Also

crit_EHI, crit_SUR, crit_EMI.

Examples

#---------------------------------------------------------------------------
# SMS-EGO 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))
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))
PF <- 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)
critcontrol <- list(refPoint = c(300, 0), currentHV = dominated_hypervolume(t(PF), c(300, 0)))
SMSEGO_grid <- apply(test.grid, 1, crit_SMS, model = model,
                     paretoFront = PF, critcontrol = critcontrol)

filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid),
               matrix(pmax(0, SMSEGO_grid), nrow = n.grid), nlevels = 50,
               main = "SMS-EGO criterion (positive part)", 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")
                            }
              )

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

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")
                            }
              )         

Maximization of multiobjective infill criterion

Description

Given a list of objects of class km and a set of tuning parameters (lower, upper and critcontrol), crit_optimizer performs the maximization of an infill criterion and delivers the next point to be visited in a multi-objective EGO-like procedure.

The latter maximization relies either on a genetic algorithm using derivatives, genoud, particle swarm algorithm pso-package, exhaustive search at pre-specified points or on a user defined method. It is important to remark that the information needed about the objective function reduces here to the vector of response values embedded in the models (no call to the objective functions or simulators (except with cheapfn)).

Usage

crit_optimizer(
  crit = "SMS",
  model,
  lower,
  upper,
  cheapfn = NULL,
  type = "UK",
  paretoFront = NULL,
  critcontrol = NULL,
  optimcontrol = NULL,
  ncores = 1
)

Arguments

Details

Extension of the function max_EI for multi-objective optimization.
Available infill criteria with crit are :

• Expected Hypervolume Improvement (EHI) crit_EHI,

• SMS criterion (SMS) crit_SMS,

• Expected Maximin Improvement (EMI) crit_EMI,

• Stepwise Uncertainty Reduction of the excursion volume (SUR) crit_SUR

Depending on the selected criterion, parameters such as a reference point for SMS and EHI or arguments for integration_design_optim with SUR can be given with critcontrol. Also options for checkPredict are available. More precisions are given in the corresponding help pages.

Value

A list with components:

• par: The best set of parameters found,

• value: The value of expected improvement at par.

References

W.R. Jr. Mebane and J.S. Sekhon (2011), Genetic optimization using derivatives: The rgenoud package for R, Journal of Statistical Software, 42(11), 1-26 doi:10.18637/jss.v042.i11

D.R. Jones, M. Schonlau, and W.J. Welch (1998), Efficient global optimization of expensive black-box functions, Journal of Global Optimization, 13, 455-492.

Examples

## Not run: 
#---------------------------------------------------------------------------
# EHI surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------

set.seed(25468)
library(DiceDesign)

d <- 2
n.obj <- 2 
fname <- "P1" 
n.grid <- 51
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
nappr <- 15 
design.grid <- round(maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design, 1)
response.grid <- t(apply(design.grid, 1, fname))
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)

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

lower <- rep(0, d)
upper <- rep(1, d)     

omEGO <- crit_optimizer(crit = "EHI", model = model,  lower = lower, upper = upper, 
                 optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2),
                 critcontrol = list(refPoint = c(300, 0)))
                 
print(omEGO)
 
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
               matrix(EHI_grid, nrow = n.grid), main = "Expected Hypervolume Improvement", 
               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(omEGO$par, col = "red", pch = 4)
                            }
              )
              
#---------------------------------------------------------------------------
# SMS surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------


SMS_grid <- apply(test.grid, 1, crit_SMS, model = model,
                  critcontrol = list(refPoint = c(300, 0)))

lower <- rep(0, d)
upper <- rep(1, d)     

omEGO2 <- crit_optimizer(crit = "SMS", model = model,  lower = lower, upper = upper, 
                 optimcontrol = list(method="genoud", pop.size = 200, BFGSburnin = 2),
                 critcontrol = list(refPoint = c(300, 0)))
                 
print(omEGO2)
 
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
               matrix(pmax(0,SMS_grid), nrow = n.grid), main = "SMS Criterion (>0)",
               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(omEGO2$par, col = "red", pch = 4)
                            }
              )
#---------------------------------------------------------------------------
# Maximin Improvement surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------


EMI_grid <- apply(test.grid, 1, crit_EMI, model = model,
                  critcontrol = list(nb_samp = 20, type ="EMI"))

lower <- rep(0, d)
upper <- rep(1, d)     

omEGO3 <- crit_optimizer(crit = "EMI", model = model,  lower = lower, upper = upper, 
                 optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
                 
print(omEGO3)
 
filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), nlevels = 50,
               matrix(EMI_grid, nrow = n.grid), main = "Expected Maximin Improvement",
               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(omEGO3$par, col = "red", pch = 4)
                            }
              )
#---------------------------------------------------------------------------
# crit_SUR surface associated with the "P1" problem at a 15 points design
#---------------------------------------------------------------------------
library(KrigInv)

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, n.obj, 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)
EEV_grid <- apply(test.grid, 1, crit_SUR, model=model, paretoFront = paretoFront,
                  critcontrol = critcontrol)
lower <- rep(0, d)
upper <- rep(1, d)
                  
omEGO4 <- crit_optimizer(crit = "SUR", model = model,  lower = lower, upper = upper, 
                 optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
print(omEGO4)

filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid),
               matrix(pmax(0,EEV_grid), n.grid), main = "EEV criterion", nlevels = 50,
               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(omEGO4$par, col = "red", pch = 4)
                            }
              )

# example using user defined optimizer, here L-BFGS-B from base optim
userOptim <- function(par, fn, lower, upper, control, ...){
  return(optim(par = par, fn = fn, method = "L-BFGS-B", lower = lower, upper = upper,
         control = control, ...))
}
omEGO4bis <- crit_optimizer(crit = "SUR", model = model,  lower = lower, upper = upper, 
                 optimcontrol = list(method = "userOptim"))
print(omEGO4bis)


#---------------------------------------------------------------------------
# crit_SMS surface with problem "P1" with 15 design points, using cheapfn
#---------------------------------------------------------------------------

# Optimization with fastfun: SMS with discrete search
# Separation of the problem P1 in two objectives: 
# the first one to be kriged, the second one with fastobj

# Definition of the fastfun
f2 <-   function(x){
  return(P1(x)[2])
}
 
SMS_grid_cheap <- apply(test.grid, 1, crit_SMS, model = list(mf1, fastfun(f2, design.grid)),
                        paretoFront = paretoFront, critcontrol = list(refPoint = c(300, 0)))


optimcontrol <- list(method = "pso")
model2 <- list(mf1)
omEGO5 <- crit_optimizer(crit = "SMS", model = model2,  lower = lower, upper = upper,
                         cheapfn = f2, critcontrol = list(refPoint = c(300, 0)),
                         optimcontrol = list(method = "genoud", pop.size = 200, BFGSburnin = 2))
print(omEGO5)

filled.contour(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid), 
               matrix(pmax(0, SMS_grid_cheap), nrow = n.grid), nlevels = 50,
               main = "SMS criterion with cheap 2nd objective (>0)", 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(omEGO5$par, col = "red", pch = 4)
                            }
              )

## End(Not run)

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

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)
                            }
              )

EGO algorithm for multiobjective optimization

Description

User-friendly wrapper of the function GParetoptim. Generates initial DOEs and kriging models (objects of class km), and executes nsteps iterations of multiobjective EGO methods.

Usage

easyGParetoptim(
  fn,
  ...,
  cheapfn = NULL,
  budget,
  lower,
  upper,
  par = NULL,
  value = NULL,
  noise.var = NULL,
  control = list(method = "SMS", trace = 1, inneroptim = "pso", maxit = 100, seed = 42),
  ncores = 1
)

Arguments

Details

Does not require specific knowledge on kriging models (objects of class km).

The problem considered is of the form: min f(x) = f_1(x), ..., f_p(x). The control argument is a list that can supply any of the following optional components:

• method: choice of multiobjective improvement function: "SMS", "EHI", "EMI" or "SUR" (see crit_SMS, crit_EHI, crit_EMI, crit_SUR),

• trace: if positive, tracing information on the progress of the optimization is produced (1 (default) for general progress, >1 for more details, e.g., warnings from genoud),

• inneroptim: choice of the inner optimization algorithm: "genoud", "pso" or "random" (see genoud and psoptim),

• maxit: maximum number of iterations of the inner loop,

• seed: to fix the random variable generator,

• refPoint: reference point for hypervolume computations (for "SMS" and "EHI" methods),

• 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.

If noise.var="given_by_fn", fn must return a list of two vectors, the first being the objective functions and the second the corresponding noise variances. See examples in GParetoptim.

For additional details or other possible arguments, see GParetoptim.

Display of results and various post-processings are available with plotGPareto.

Value

A list with components:

• par: all the non-dominated points found,

• value: the matrix of objective values at the points given in par,

• history: a list containing all the points visited by the algorithm (X) and their corresponding objectives (y),

• model: a list of objects of class km, corresponding to the last kriging models fitted.

Note that in the case of noisy problems, value and history$y.denoised are denoised values. The original observations are available in the slot history$y.

Author(s)

Victor Picheny (INRA, Castanet-Tolosan, France)

Mickael Binois (Mines Saint-Etienne/Renault, France)

References

M. T. Emmerich, A. H. Deutz, J. W. Klinkenberg (2011), Hypervolume-based expected improvement: Monotonicity properties and exact computation, Evolutionary Computation (CEC), 2147-2154.

V. Picheny (2015), Multiobjective optimization using Gaussian process emulators via stepwise uncertainty reduction, Statistics and Computing, 25(6), 1265-1280.

T. Wagner, M. Emmerich, A. Deutz, W. Ponweiser (2010), On expected-improvement criteria for model-based multi-objective optimization. Parallel Problem Solving from Nature, 718-727, Springer, Berlin.

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

M. Binois, V. Picheny (2019), GPareto: An R Package for Gaussian-Process-Based Multi-Objective Optimization and Analysis, Journal of Statistical Software, 89(8), 1-30, doi:10.18637/jss.v089.i08.

Examples

#---------------------------------------------------------------------------
# 2D objective function, 4 cases
#---------------------------------------------------------------------------
## Not run: 
set.seed(25468)
n_var <- 2 
fname <- ZDT3
lower <- rep(0, n_var)
upper <- rep(1, n_var)

#---------------------------------------------------------------------------
# 1- Expected Hypervolume Improvement optimization, using pso
#---------------------------------------------------------------------------
res <- easyGParetoptim(fn=fname, lower=lower, upper=upper, budget=15, 
                   control=list(method="EHI", inneroptim="pso", maxit=20))
par(mfrow=c(1,2))
plotGPareto(res)
title("Pareto Front")
plot(res$history$X, main="Pareto set", col = "red", pch = 20)
points(res$par, col="blue", pch = 17)

#---------------------------------------------------------------------------
# 2- SMS Improvement optimization using random search, with initial DOE given
#---------------------------------------------------------------------------
library(DiceDesign)
design.init   <- maximinESE_LHS(lhsDesign(10, n_var, seed = 42)$design)$design
response.init <- t(apply(design.init, 1, fname))
res <- easyGParetoptim(fn=fname, par=design.init, value=response.init, lower=lower, upper=upper, 
                       budget=15, control=list(method="SMS", inneroptim="random", maxit=100))
par(mfrow=c(1,2))
plotGPareto(res)
title("Pareto Front")
plot(res$history$X, main="Pareto set", col = "red", pch = 20)
points(res$par, col="blue", pch = 17)

#---------------------------------------------------------------------------
# 3- Stepwise Uncertainty Reduction optimization, with one fast objective function
#---------------------------------------------------------------------------
fname <- camelback
cheapfn <- function(x) {
if (is.null(dim(x))) return(-sum(x))
else return(-rowSums(x))
}
res <- easyGParetoptim(fn=fname, cheapfn=cheapfn, lower=lower, upper=upper, budget=15, 
                   control=list(method="SUR", inneroptim="pso", maxit=20))
par(mfrow=c(1,2))
plotGPareto(res)
title("Pareto Front")
plot(res$history$X, main="Pareto set", col = "red", pch = 20)
points(res$par, col="blue", pch = 17)

#---------------------------------------------------------------------------
# 4- Expected Hypervolume Improvement optimization, using pso, noisy fn
#---------------------------------------------------------------------------
noise.var <- c(0.1, 0.2)
funnoise <- function(x) {ZDT3(x) + sqrt(noise.var)*rnorm(n=2)}
res <- easyGParetoptim(fn=funnoise, lower=lower, upper=upper, budget=30, noise.var=noise.var,
                       control=list(method="EHI", inneroptim="pso", maxit=20))
par(mfrow=c(1,2))
plotGPareto(res)
title("Pareto Front")
plot(res$history$X, main="Pareto set", col = "red", pch = 20)
points(res$par, col="blue", pch = 17)

#---------------------------------------------------------------------------
# 5- Stepwise Uncertainty Reduction optimization, functional noise
#---------------------------------------------------------------------------
funnoise <- function(x) {ZDT3(x) + sqrt(abs(0.1*x))*rnorm(n=2)}
noise.var <- function(x) return(abs(0.1*x))

res <- easyGParetoptim(fn=funnoise, lower=lower, upper=upper, budget=30, noise.var=noise.var,
                     control=list(method="SUR", inneroptim="pso", maxit=20))
par(mfrow=c(1,2))
plotGPareto(res)
title("Pareto Front")
plot(res$history$X, main="Pareto set", col = "red", pch = 20)
points(res$par, col="blue", pch = 17)

## End(Not run)

Fast-to-evaluate function wrapper

Description

Modification of an R function to be used with methods predict and update (similar to a km object). It creates an S4 object which contains the values corresponding to evaluations of other costly observations. It is useful when an objective can be evaluated fast.

Usage

fastfun(fn, design, response = NULL)

Arguments

Value

An object of class fastfun-class.

Examples

########################################################
## Example with a fast to evaluate objective
########################################################
## Not run: 
set.seed(25468)
library(DiceDesign)

d <- 2 

fname <- P1
n.grid <- 21
nappr <- 11 
design.grid <- maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))
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])
model <- list(mf1, mf2)

nsteps <- 5 
lower <- rep(0, d)
upper <- rep(1, d)

# Optimization reference: SMS with discrete search
optimcontrol <- list(method = "pso")
omEGO1 <- GParetoptim(model = model, fn = fname, crit = "SMS", nsteps = nsteps,
                     lower = lower, upper = upper, optimcontrol = optimcontrol)
print(omEGO1$par)
print(omEGO1$values)
plot(response.grid, xlim = c(0,300), ylim = c(-40,0), pch = 17, col = "blue")
points(omEGO1$values, pch = 20, col ="green") 

# Optimization with fastfun: SMS with discrete search
# Separation of the problem P1 in two objectives: 
# the first one to be kriged, the second one with fastobj
f1 <-   function(x){
  if(is.null(dim(x))) x <- matrix(x, nrow = 1) 
  b1 <- 15*x[,1] - 5
  b2 <- 15*x[,2]
  return(  (b2 - 5.1*(b1/(2*pi))^2 + 5/pi*b1 - 6)^2 +10*((1 - 1/(8*pi))*cos(b1) + 1))
}

f2 <-   function(x){
  if(is.null(dim(x))) x <- matrix(x, nrow = 1) 
  b1<-15*x[,1] - 5
  b2<-15*x[,2]
  return(-sqrt((10.5 - b1)*(b1 + 5.5)*(b2 + 0.5))
         - 1/30*(b2 - 5.1*(b1/(2*pi))^2 - 6)^2
         - 1/3*((1 - 1/(8*pi))*cos(b1) + 1)) 
}

optimcontrol <- list(method = "pso")
model2 <- list(mf1)
omEGO2 <- GParetoptim(model = model2, fn = f1, cheapfn = f2, crit = "SMS", nsteps = nsteps,
                     lower = lower, upper = upper, optimcontrol = optimcontrol)
print(omEGO2$par)
print(omEGO2$values)

points(omEGO2$values, col = "red", pch = 15)

## End(Not run)

Class for fast to compute objective.

Description

Class for fast to compute objective.

Usage

## S4 method for signature 'fastfun'
predict(object, newdata, ...)

## S4 method for signature 'fastfun'
update(object, newX, newy, ...)

## S4 method for signature 'fastfun'
simulate(object, nsim, seed, newdata, cond, nugget.sim, checkNames, ...)

Arguments

Methods (by generic)

• predict(fastfun): Predict(by evaluating fun) the result at a new observation.

• update(fastfun): Update the X and y slots with a new design and observation.

• simulate(fastfun): Simulate responses values (for compatibility with methods using [DiceKriging::simulate()])

Slots

d

spatial dimension,

n

observations number,

X

the design of experiments, size n x d,

y

the observations, size n x 1,

fun

the evaluator function.

Objects from the Class

To create a fastfun object, use fastfun. See also this function for more details and examples.

Get design corresponding to an objective target

Description

Find the design that maximizes the probability of dominating a target given by the user.

Usage

getDesign(model, target, lower, upper, optimcontrol = NULL)

Arguments

Value

A list with components:

• par: best design found,

• value: probabilitity that the design dominates the target,

• mean: kriging mean of the objectives at the design,

• sd: prediction standard deviation at the design.

Examples

## Not run: 

#---------------------------------------------------------------------------
# Example of interactive optimization
#---------------------------------------------------------------------------

set.seed(25468)
library(DiceDesign)

d <- 2
n.obj <- 2 
fun <- "P1" 
n.grid <- 51
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
nappr <- 20 
design.grid <- round(maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design, 1)
response.grid <- t(apply(design.grid, 1, fun))
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)
lower <- rep(0, d); upper <- rep(1, d)

sol <- GParetoptim(model, fun, crit = "SUR", nsteps = 5, lower = lower, upper = upper) 

plotGPareto(sol)

target1 <- c(15, -25)
points(x = target1[1], y = target1[2], col = "black", pch = 13)

nDesign <- getDesign(sol$lastmodel, target = target1, lower = rep(0, d), upper = rep(1, d))
points(t(nDesign$mean), col = "green", pch = 20)

target2 <- c(48, -27)
points(x = target2[1], y = target2[2], col = "black", pch = 13)
nDesign2 <- getDesign(sol$lastmodel, target = target2, lower = rep(0, d), upper = rep(1, d))
points(t(nDesign2$mean), col = "darkgreen", pch = 20)

## End(Not run)

Function to build integration points (for the SUR criterion)

Description

Modification of the function integration_design from the package KrigInv-package to be usable for SUR-based optimization. Handles two or three objectives. Available important sampling schemes: none so far.

Usage

integration_design_optim(
  SURcontrol = NULL,
  d = NULL,
  lower,
  upper,
  model = NULL,
  min.prob = 0.001
)

Arguments

Details

The SURcontrol argument is a list with possible entries integration.points, integration.weights, n.points, n.candidates, distrib, init.distrib and init.distrib.spec. It can be used in one of the three following ways:

• A) If nothing is specified, 100 * d points are chosen using the Sobol sequence;

• B) One can directly set the field integration.points (p * d matrix) for prespecified integration points. In this case these integration points and the corresponding vector integration.weights will be used for all the iterations of the algorithm;

• C) If the field integration.points is not set then the integration points are renewed at each iteration. In that case one can control the number of integration points n.points (default: 100*d) and a specific distribution distrib. Possible values for distrib are: "sobol", "MC" and "SUR" (default: "sobol"):

  • C.1) The choice "sobol" corresponds to integration points chosen with the Sobol sequence in dimension d (uniform weight);

  • C.2) The choice "MC" corresponds to points chosen randomly, uniformly on the domain;

  • C.3) The choice "SUR" corresponds to importance sampling distributions (unequal weights).
    When important sampling procedures are chosen, n.points points are chosen using importance sampling among a discrete set of n.candidates points (default: n.points*10) which are distributed according to a distribution init.distrib (default: "sobol"). Possible values for init.distrib are the space filling distributions "sobol" and "MC" or an user defined distribution "spec". The "sobol" and "MC" choices correspond to quasi random and random points in the domain. If the "spec" value is chosen the user must fill in manually the field init.distrib.spec to specify himself a n.candidates * d matrix of points in dimension d.

Value

A list with components:

• integration.points, p x d matrix of p points used for the numerical calculation of integrals,

• integration.weights, a vector of size p corresponding to the weight of each point. If all the points are equally weighted, integration.weights is set to NULL.

References

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

See Also

GParetoptim crit_SUR integration_design

Non-dominated points with respect to a reference

Description

Determines which elements in a set are dominated by reference points

Usage

nonDomSet(points, ref)

Arguments

Examples

## Not run: 
d <- 6
n <- 1000
n2 <- 1000

test <- matrix(runif(d * n), n)
ref <- matrix(runif(d * n), n)
indPF <- which(!is_dominated(t(ref)))

system.time(res <- nonDomSet(test, ref[indPF,,drop = FALSE]))

res2 <- rep(NA, n2)
library(emoa)
t0 <- Sys.time()
for(i in 1:n2){
  res2[i] <- !is_dominated(t(rbind(test[i,, drop = FALSE], ref[indPF,])))[1]
}
print(Sys.time() - t0)

all(res == res2)

## End(Not run)

Plot multi-objective optimization results and post-processing

Description

Display results of multi-objective optimization returned by either GParetoptim or easyGParetoptim, possibly completed with various post-processings of uncertainty quantification.

Usage

plotGPareto(
  res,
  add = FALSE,
  UQ_PF = FALSE,
  UQ_PS = FALSE,
  UQ_dens = FALSE,
  lower = NULL,
  upper = NULL,
  control = list(pch = 20, col = "red", PF.line.col = "cyan", PF.pch = 17, PF.points.col
    = "blue", VE.line.col = "cyan", nsim = 100, npsim = 1500, gridtype = "runif",
    displaytype = "persp", printVD = TRUE, use.rgl = TRUE, bounds = NULL, meshsize3d =
    50, theta = -25, phi = 10, add_denoised_PF = TRUE)
)

Arguments

Details

By default, plotGPareto displays the Pareto front delimiting the non-dominated area with 2 objectives, by a perspective view with 3 objectives and using parallel coordinates with more objectives.

Setting one or several of UQ_PF, UQ_PS and UQ_dens allows to run and display post-processing tools that assess the precision and confidence of the optimization run, either in the objective (UQ_PF) or the variable spaces (UQ_PS, UQ_dens). Note that these options are computationally intensive.

Various parameters can be used for the display of results and/or passed to subsequent function:

• col, pch correspond the color and plotting character for observations,

• PF.line.col, PF.pch, PF.points.col define the color of the line denoting the current Pareto front, the plotting character and color of non-dominated observations, respectively,

• nsim, npsim and gridtype define the number of conditional simulations performed with [DiceKriging::simulate()] along with the number of simulation points (in case UQ_PF and/or UQ_dens are TRUE),

• gridtype to define how simulation points are selected; alternatives are 'runif' (default) for uniformly sampled points, 'LHS' for a Latin Hypercube design using lhsDesign and 'grid2d' for a two dimensional grid,

• f1lim, f2lim can be passed to CPF,

• resolution, option, nintegpoints are to be passed to plot_uncertainty

• displaytype type of display for UQ_dens, see plot.kde,

• printVD logical, if TRUE and UQ_PF is TRUE as well, print the value of the Vorob'ev deviation,

• use.rgl if TRUE, use rgl for 3D plots, else persp is used,

• bounds if use.rgl is TRUE, optional 2*nobj matrix of boundaries, see plotParetoEmp

• meshsize3d mesh size of the perspective view for 3-objective problems,

• theta, phi angles for perspective view of 3-objective problems,

• add_denoised_PF if TRUE, in the noisy case, add the Pareto front from the estimated mean of the observations.

References

M. Binois, D. Ginsbourger and O. Roustant (2015), Quantifying Uncertainty on Pareto Fronts with Gaussian process conditional simulations, European Journal of Operational Research, 243(2), 386-394.

A. Inselberg (2009), Parallel coordinates, Springer.

Examples

## Not run: 
#---------------------------------------------------------------------------
# 2D objective function
#---------------------------------------------------------------------------
set.seed(25468)
n_var <- 2 
fname <- P1
lower <- rep(0, n_var)
upper <- rep(1, n_var)
res <- easyGParetoptim(fn=fname, lower=lower, upper=upper, budget=15, 
control=list(method="EHI", inneroptim="pso", maxit=20))

## Pareto front only
plotGPareto(res)

## With post-processing
plotGPareto(res, UQ_PF = TRUE, UQ_PS = TRUE, UQ_dens = TRUE)

## With noise
noise.var <- c(10, 2)
funnoise <- function(x) {P1(x) + sqrt(noise.var)*rnorm(n=2)}
res2 <- easyGParetoptim(fn=funnoise, lower=lower, upper=upper, budget=15, noise.var=noise.var,
                       control=list(method="EHI", inneroptim="pso", maxit=20))
                       
plotGPareto(res2, control=list(add_denoised_PF=FALSE)) # noisy observations only
plotGPareto(res2)

#---------------------------------------------------------------------------
# 3D objective function
#---------------------------------------------------------------------------
set.seed(1)
n_var <- 3 
fname <- DTLZ1
lower <- rep(0, n_var)
upper <- rep(1, n_var)
res3 <- easyGParetoptim(fn=fname, lower=lower, upper=upper, budget=50, 
control=list(method="EHI", inneroptim="pso", maxit=20))

## Pareto front only
plotGPareto(res3)

## With noise
noise.var <- c(10, 2, 5)
funnoise <- function(x) {fname(x) + sqrt(noise.var)*rnorm(n=3)}
res4 <- easyGParetoptim(fn=funnoise, lower=lower, upper=upper, budget=100, noise.var=noise.var,
                       control=list(method="EHI", inneroptim="pso", maxit=20))
                       
plotGPareto(res4, control=list(add_denoised_PF=FALSE)) # noisy observations only
plotGPareto(res4)   

## End(Not run)

Pareto front visualization

Description

Plot the Pareto front with step functions.

Usage

plotParetoEmp(
  nondominatedPoints,
  add = TRUE,
  max = FALSE,
  bounds = NULL,
  alpha = 0.5,
  ...
)

Arguments

Examples

#------------------------------------------------------------
# Simple example
#------------------------------------------------------------

x <- c(0.2, 0.4, 0.6, 0.8)
y <- c(0.8, 0.7, 0.5, 0.1)

plot(x, y, col = "green", pch = 20) 

plotParetoEmp(cbind(x, y), col = "green")
## Alternative
plotParetoEmp(cbind(x, y), col = "red", add = FALSE)

## With maximization
plotParetoEmp(cbind(x, y), col = "blue", max = TRUE)

## 3D plots
library(rgl)
set.seed(5)
X <- matrix(runif(60), ncol=3)
Xnd <- t(nondominated_points(t(X)))
plot3d(X)
plot3d(Xnd, col="red", size=8, add=TRUE)
plot3d(x=min(Xnd[,1]), y=min(Xnd[,2]), z=min(Xnd[,3]), col="green", size=8, add=TRUE)
X.range <- diff(apply(X,2,range))
bounds <- rbind(apply(X,2,min)-0.1*X.range,apply(X,2,max)+0.1*X.range)
plotParetoEmp(nondominatedPoints = Xnd, add=TRUE, bounds=bounds, alpha=0.5)

Visualisation of Pareto front and set

Description

Plot the Pareto front and set for 2 variables 2 objectives test problems with evaluations on a grid.

Usage

plotParetoGrid(fname = "ZDT1", xlim = c(0, 1), ylim = c(0, 1), n.grid = 100)

Arguments

Examples

#------------------------------------------------------------
# Examples with test functions
#------------------------------------------------------------

plotParetoGrid("ZDT3", n.grid = 21)

plotParetoGrid("P1", n.grid = 21)

plotParetoGrid("MOP2", xlim = c(0, 1), ylim = c(0, 1), n.grid = 21) 

Display the Symmetric Deviation Function

Description

Display the Symmetric Deviation Function for an object of class CPF.

Usage

plotSymDevFun(CPF, n.grid = 100)

Arguments

Details

Display observations in red and the corresponding Pareto front by a step-line. The blue line is the estimation of the location of the Pareto front of the kriging models, named Vorob'ev expectation. In grayscale is the intensity of the deviation (symmetrical difference) from the Vorob'ev expectation for couples of conditional simulations.

References

M. Binois, D. Ginsbourger and O. Roustant (2015), Quantifying Uncertainty on Pareto Fronts with Gaussian process conditional simulations, European Journal of Operational Research, 243(2), 386-394.

Examples

library(DiceDesign)
set.seed(42)

nvar <- 2

# Test function
fname = "P1"

# Initial design
nappr <- 10
design.grid <- maximinESE_LHS(lhsDesign(nappr, nvar, seed = 42)$design)$design
response.grid <- t(apply(design.grid, 1, fname))

ParetoFront <- t(nondominated_points(t(response.grid)))

# kriging models : matern5_2 covariance structure, linear trend, no nugget effect
mf1 <- km(~., design = design.grid, response = response.grid[, 1])
mf2 <- km(~., design = design.grid, response = response.grid[, 2])

# Conditional simulations generation with random sampling points 
nsim <- 10 # increase for better results
npointssim <- 80 # increase for better results
Simu_f1 = matrix(0, nrow = nsim, ncol = npointssim)
Simu_f2 = matrix(0, nrow = nsim, ncol = npointssim)
design.sim = array(0,dim = c(npointssim, nvar, nsim))

for(i in 1:nsim){
  design.sim[,, i] <- matrix(runif(nvar*npointssim), npointssim, nvar)
  Simu_f1[i,] = simulate(mf1, nsim = 1, newdata = design.sim[,, i], cond = TRUE,
                         checkNames = FALSE, nugget.sim = 10^-8)
  Simu_f2[i,] = simulate(mf2, nsim = 1, newdata = design.sim[,, i], cond=TRUE, 
                         checkNames = FALSE, nugget.sim = 10^-8)
}

# Attainment, Voreb'ev expectation and deviation estimation
CPF1 <- CPF(Simu_f1, Simu_f2, response.grid, ParetoFront)

# Symmetric deviation function
plotSymDevFun(CPF1)

Symmetrical difference of RNP sets

Description

Plot the symmetrical difference between two Random Non-Dominated Point (RNP) sets.

Usage

plotSymDifRNP(set1, set2, xlim, ylim, fill = "black", add = "FALSE", ...)

Arguments

Examples

#------------------------------------------------------------
# Simple example
#------------------------------------------------------------
set1 <- rbind(c(0.2, 0.35, 0.5, 0.8),
              c(0.8, 0.6, 0.55, 0.3))

set2 <- rbind(c(0.3, 0.4),
              c(0.7, 0.4))

plotSymDifRNP(set1, set2, xlim = c(0, 1), ylim = c(0, 1), fill = "grey")
points(t(set1), col = "red", pch = 20)
points(t(set2), col = "blue", pch = 20)

Plot uncertainty

Description

Displays the probability of non-domination in the variable space. In dimension larger than two, projections in 2D subspaces are displayed.

Usage

plot_uncertainty(
  model,
  paretoFront = NULL,
  type = "pn",
  lower,
  upper,
  resolution = 51,
  option = "mean",
  nintegpoints = 400
)

Arguments

Details

Function inspired by the function print_uncertainty and print_uncertainty_nd from the package KrigInv-package. Non-dominated observations are represented with green diamonds, dominated ones by yellow triangles.

Examples

## Not run:  
#---------------------------------------------------------------------------
# 2D, bi-objective function
#---------------------------------------------------------------------------
set.seed(25468)
n_var <- 2 
fname <- P1
lower <- rep(0, n_var)
upper <- rep(1, n_var)
res1 <- easyGParetoptim(fn=fname, lower=lower, upper=upper, budget=15, 
control=list(method="EHI", inneroptim="pso", maxit=20))

plot_uncertainty(res1$model, lower = lower, upper = upper)

#---------------------------------------------------------------------------
# 4D, bi-objective function
#---------------------------------------------------------------------------
set.seed(25468)
n_var <- 4
fname <- DTLZ2
lower <- rep(0, n_var)
upper <- rep(1, n_var)
res <- easyGParetoptim(fn=fname, lower=lower, upper=upper, budget = 40, 
control=list(method="EHI", inneroptim="pso", maxit=40))

plot_uncertainty(res$model, lower = lower, upper = upper, resolution = 31)

## End(Not run) 

Predict function for list of km models.

Description

Predict function for list of km models.

Usage

predict_kms(
  model,
  newdata,
  type,
  se.compute = TRUE,
  cov.compute = FALSE,
  light.return = TRUE,
  bias.correct = FALSE,
  checkNames = FALSE,
  ...
)

Arguments

Details

So far only light.return = TRUE handled. For the cov field, a list of cov matrices is returned.