eas {HyetosMinute}R Documentation

The Evalutionary Annealing-Simplex Method

Description

An enhanced version of the evolutionary annealing-simplex optimization method for the estimation of Bartlett-Lewis model parameters.

Usage

eas(n,m,xmin,xmax,xlow,xup,fn,maxeval=1500,ftol=1.e-07,ratio=0.99,pmut=0.9,
     beta=2,maxclimbs=5,...)

Arguments

n

An integer that specifies the problem dimension (number of control variables).

m

A positive integer, greater than n, that specifies the population size (m >= n + 1).

xmin

A vector, with length equal to n, that defines the interior lower parameter bounds (feasible space of initial population).

xmax

A vector, with length equal to n, that defines the interior upper parameter bounds (feasible space of initial population).

xlow

A vector, with length equal to n, that defines the exterior lower parameter bounds.

xup

A vector, with length equal to n, that defines the exterior upper parameter bounds.

fn

Objective function that is to be optimized. A vector function that takes a real vector as argument and returns the value of the function at that point.

maxeval

A positive integer that specifies the maximum number of function evaluations. Default is maxeval=1500. Suggested value is maxeval>100*n.

ftol

A positive number that specifies the fractional convergence tolerance to be achieved in the function value. Default is ftol=1.e-07.

ratio

A positive number, typically between 0.80-0.99, that specifies the fraction of temperature reduction, when a local minimum is found. Default is ratio=0.99.

pmut

A positive number, between 0.5-0.95, that specifies the probability of accepting an offspring generated via mutation. Default is pmut=0.9. Higher values are suggested for very hard problems, when it is essential to increase randomness.

beta

A positive integer,greater than 1, that specifies the annealing schedule parameter. Default is beta=2.

maxclimbs

A positive integer, typically between 3-5, that specifies the maximum number of uphill steps. Default is maxclimbs=5.

...

Further arguments to be passed in the objective function.

Details

The evalutionary annealing-simplex algorithm (2002), for solving global optimisation problems have been proposed by Efstratiadis and Koutsoyiannis (2002). The evalutionary annealing-simplex algorithm (eas) is a probabilistic heuristic global optimization technique, combining the robustness of simulated annealing in rough response surface, with the efficiency of hill-climbing methods in convex areas.

During one generation, the population evolves as follows: First, a simplex-based pattern is formulated, using random sampling. Next, a candidate individual is selected to die, according to a modified objective function of the form:

g(x) = f(x) + uT

where f is the original objective function, T is the current “temperature” and u is a random number from the uniform distribution. The temperature is gradually reduced, according to an appropriate annealing cooling schedule, automatically adapted during the evolution. Consequently, the probability of replacing individuals with poor performance increases, since the procedure gradually moves from a random walk to a local search.

The recombination operator is based on the well-known downhill simplex transitions (Nelder and Mead, 1965). According to the relatives values of the objective function at the vertices, the simplex is reflected, expanded, contracted or shrinks, where quasi-stochastic sclae factors are employed instead of constant ones. To ensure more flexibility, additional transformations are introduced, namely multiple expansion towards the direction of reflection, when a downhill path (i.e., the gradient of the function) is located, and similar expansions but on the opposite (uphill) direction, in order to escape from the nearest local minimum. If any of the above transitions improves the function value, the new individual is generated through nutation. The related oerator employs a random pertubation scheme outside of the usual range of the population, as determined on the basis of the average and standard deviation values of its coordinates.

R adaption, with some modifications, by Panagiotis Kossieris and Hristos Tyralis, from the original Pascal-Delphi code by Andreas Efstratiadis (2007).

In the current version, the calculation of the reflection point is changed, and the simplex is not reflected toward the direction of the geometrical centroid but towards a weighted centroid, where the weigths are assigned on the basis of the objective fuction value. This allowes for a more accurate estimator of the gradient, which ensures a much more efficient search in multidimensional feasible spaces.

For more details, the reader may refer either to the original work of Efstratiadis (2001) and Efstratiadis and Koutsoyiannis (2002), or to the tutorial style introduction of HyetosR package.

Value

A list with the following components:

bestpar

A vector containing the best set of parameters found.

bestval

The value of fn corresponding to bestpar.

nfeval

Number of function fn evaluations.

niter

Number of iterations taken by algorithm.

ftolpop

Fractional convergence tolerance of population generated at the last iteration.

pop

The population generated at the last iteration.

Author(s)

Kossieris Panagiotis pankoss@hotmail.com, with Hristos Tyralis montchrister@gmail.com and Andreas Efstratiadis A.Efstratiadis@itia.ntua.gr.

References

The methodology of the evolutionary annealing-simplex optimization method and details of its application are described in:

Find the original Pascal-Delphi code of optimisation algorithms on http://itia.ntua.gr/en/docinfo/524/.

Examples

# Historical statistics 
mean1=0.1226;var1 = 0.6323;cov1lag1 = 0.3271;pdr1 = 0.9183
mean6=0.7358;var6 = 10.1490;cov6lag1 = 4.0773;pdr6 = 0.8251
mean12 = 1.4705;var12 = 29.357;cov12lag1 =7.6865;pdr12 = 0.7476
mean24 = 2.9410;var24 = 76.667;cov24lag1 =10.2886;pdr24 = 0.6238

### Fitting the Original Bartlett-Lewis model (with contant parameter n) ###

# Define the objective function that will be used in model fitting:

fopt <- function(x,Weibull) {
  
  l<-x[1];g<-x[2];b<-x[3];n<-x[4];mx<-x[5];sxmx<-x[6]
  
  w1=1;w2=1;w3=1;w4=1
  S1 <- w1*((meanBLRPM(l,g,b,n,mx,h=1)/mean1)-1)^(2)+
    w2*((varBLRPM(l,g,b,n,mx,sxmx,h=1,Weibull)/var1)-1)^(2)+
    w3*((covarBLRPM(l,g,b,n,mx,sxmx,h=1,lag=1,Weibull)/cov1lag1)-1)^(2)+
    w4*((pdrBLRPM(l,g,b,n,h=1)/pdr1)-1)^(2)
  S6 <- w1*((meanBLRPM(l,g,b,n,mx,h=6)/mean6)-1)^(2)+
    w2*((varBLRPM(l,g,b,n,mx,sxmx,h=6,Weibull)/var6)-1)^(2)+
    w3*((covarBLRPM(l,g,b,n,mx,sxmx,h=6,lag=1,Weibull)/cov6lag1)-1)^(2)+
    w4*((pdrBLRPM(l,g,b,n,h=6)/pdr6)-1)^(2)
  S12 <- w1*((meanBLRPM(l,g,b,n,mx,h=12)/mean12)-1)^(2)+
    w2*((varBLRPM(l,g,b,n,mx,sxmx,h=12,Weibull)/var12)-1)^(2)+
    w3*((covarBLRPM(l,g,b,n,mx,sxmx,h=12,lag=1,Weibull)/cov12lag1)-1)^(2)+
    w4*((pdrBLRPM(l,g,b,n,h=12)/pdr12)-1)^(2)
  S24 <- w1*((meanBLRPM(l,g,b,n,mx,h=24)/mean24)-1)^(2)+
    w2*((varBLRPM(l,g,b,n,mx,sxmx,h=24,Weibull)/var24)-1)^(2)+
    w3*((covarBLRPM(l,g,b,n,mx,sxmx,h=24,lag=1,Weibull)/cov24lag1)-1)^(2)+
    w4*((pdrBLRPM(l,g,b,n,h=24)/pdr24)-1)^(2)
  S<-S1+S6+S12+S24
  if(is.infinite(S)) {S<-10^8}
  if(is.na(S)) {S<-10^8}
  return(S)
}

# In the case of EXPONENTIAL distribution for cell intensity, 
# set upper and lower bound of sxmx (last element) equal to 1
xlow<-c(0.001,0.0001,0.001,0.001,0.001,1)
xup<-c(0.1,5,5,5,5,1)

# In the case of GAMMA or WEIBULL distribution for cell intensity:
xlow<-c(0.001,0.0001,0.001,0.001,0.001,1)
xup<-c(0.1,5,5,5,5,5)

# Fitting the model with evalutionary annealing-simplex method
# EXPONENTIAL or GAMMA distribution for cell intensities
modecal<-eas(n=6,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=FALSE)

# WEIBULL distribution for cell intensities
modecal<-eas(n=6,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=TRUE)

l<-modecal$bestpar[1];g<-modecal$bestpar[2];b<-modecal$bestpar[3]
n<-modecal$bestpar[4];mx<-modecal$bestpar[5];sxmx<-modecal$bestpar[6]

# To use the optimised set of parameters in the functions of HyetosMinute 
# please be sure that for parameters mx the length units are millimeters (mm) 
# and for parameters l, g, b, n, and mx the time units are days (d). 
# Due to this, the following unit conversions (hourly to daily) are needed:

l<-l*24;g<-g*24;b<-b*24;n<-n*24
mx<-mx*24;sxmx<-sxmx

# The parameter set used in HyetosMinute functions will be

BLpar=list(lambda=l,phi=g/n,kappa=b/n,alpha=150,v=n,mx=mx,sxmx=sxmx)
# To implement Weibull distribution set: CellIntensityProp=list(Weibull=TRUE,iota=NA)  
# To implement Exponential or Gamma distribution set: CellIntensityProp=list(Weibull=FALSE,iota=NA)  


### Fitting the Random Bartlett-Lewis model (with varying parameter n) ###

# Define the objective function that will be used in model fitting:

fopt <- function(x,Weibull) {
  
  a<-x[1];l<-x[2];v<-x[3];k<-x[4];f<-x[5];mx<-x[6];sxmx<-x[7]
  
  w1=1;w2=1;w3=1;w4=1
  
  S1 <- w1*((meanRPBLRPM(a,l,v,k,f,mx,h=1)/mean1)-1)^(2)+
    w2*((varRPBLRPM(a,l,v,k,f,mx,sxmx,h=1,Weibull)/var1)-1)^(2)+
    w3*((covarRPBLRPM(a,l,v,k,f,mx,sxmx,h=1,Weibull,lag=1)/cov1lag1)-1)^(2)+
    w4*((pdrRPBLRPM(a,l,v,k,f,h=1)/pdr1)-1)^(2)
  S6 <- w1*((meanRPBLRPM(a,l,v,k,f,mx,h=6)/mean6)-1)^(2)+
    w2*((varRPBLRPM(a,l,v,k,f,mx,sxmx,h=6,Weibull)/var6)-1)^(2)+
    w3*((covarRPBLRPM(a,l,v,k,f,mx,sxmx,h=6,Weibull,lag=1)/cov6lag1)-1)^(2)+
    w4*((pdrRPBLRPM(a,l,v,k,f,h=6)/pdr6)-1)^(2)
  S12 <- w1*((meanRPBLRPM(a,l,v,k,f,mx,h=12)/mean12)-1)^(2)+
    w2*((varRPBLRPM(a,l,v,k,f,mx,sxmx,h=12,Weibull)/var12)-1)^(2)+
    w3*((covarRPBLRPM(a,l,v,k,f,mx,sxmx,h=12,Weibull,lag=1)/cov12lag1)-1)^(2)+
    w4*((pdrRPBLRPM(a,l,v,k,f,h=12)/pdr12)-1)^(2)
  S24 <- w1*((meanRPBLRPM(a,l,v,k,f,mx,h=24)/mean24)-1)^(2)+
    w2*((varRPBLRPM(a,l,v,k,f,mx,sxmx,h=24,Weibull)/var24)-1)^(2)+
    w3*((covarRPBLRPM(a,l,v,k,f,mx,sxmx,h=24,Weibull,lag=1)/cov24lag1)-1)^(2)+
    w4*((pdrRPBLRPM(a,l,v,k,f,h=24)/pdr24)-1)^(2)
  S<-S1
  # +S6+S12+S24
  if(is.infinite(S)) {S<-10^8}
  if(is.na(S)) {S<-10^8}
  return(S)
}

# In the case of EXPONENTIAL distribution for cell intensity, 
# set upper and lower bound of sxmx (last element) equal to 1
xlow <- c(1.0001,0.001,0.001,0.001,0.001,0.001,1)
xup <- c(15,0.1,20,20,1,50,1)
# In the case of GAMMA or WEIBULL distribution for cell intensity:
xlow <- c(1.0001,0.001,0.001,0.001,0.001,0.001,1)
xup <- c(15,0.1,20,20,1,50,5)

## Fitting the model with evalutionary annealing-simplex method
# EXPONENTIAL or GAMMA distribution for cell intensities
modecal<-eas(n=7,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=FALSE)

# WEIBULL distribution for cell intensities
modecal<-eas(n=7,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=TRUE)

a<-modecal$bestpar[1];l<-modecal$bestpar[2];v<-modecal$bestpar[3]
k<-modecal$bestpar[4];f<-modecal$bestpar[5];mx<-modecal$bestpar[6]
sxmx<-modecal$bestpar[7]

# To use the optimised set of parameters in the functions of HyetosMinute 
# please be sure that for parameters mx the length units are millimeters (mm) 
# and for parameters l, v, and mx the time units are days (d). 
# Due to this, the following unit conversions (hourly to daily) are needed:

l<-l*24;v<-v/24;mx<-mx*24

# The parameter set used in HyetosMinute functions will be

BLpar=list(lambda=l,phi=f,kappa=k,alpha=a,v=v,mx=mx,sxmx=sxmx)
# To implement Weibull distribution set: CellIntensityProp=list(Weibull=TRUE,iota=NA)  
# To implement Exponential or Gamma distribution set: CellIntensityProp=list(Weibull=FALSE,iota=NA)  


### Random Bartlett-Lewis model with dependent cell intensity-duration ###

## Historical statistics 
mean5min=0.0073 ;var5min=0.00137;cov5minlag1=0.00108; pdr5min=0.89333 
mean1=0.08738 ;var1=0.11229;cov1lag1=0.0662;pdr1=0.82894
mean6=0.52145 ;var6=2.05341 ;cov6lag1=0.81562;pdr6=0.67014
mean12=1.0414 ;var12=5.79845;cov12lag1=1.92132;pdr12=0.55893
mean24=2.07883;var24=15.00291;cov24lag1=4.93788;pdr24=0.42905

## Define the objective function that will be used in model fitting:

fopt <- function(x,Weibull) {
  
  a<-x[1];l<-x[2];v<-x[3];k<-x[4];f<-x[5];i<-x[6];sxmx<-x[7]
  
  w1=1;w2=1;w3=1;w4=1
  S5min <- w1*((meanBLRPRX(l,k,f,i,h=5/60)/mean5min)-1)^(2)+
    w2*((varBLRPRX(a,l,v,k,f,i,sxmx,h=5/60,Weibull)/var5min)-1)^(2)+
    w3*((covarBLRPRX(a,l,v,k,f,i,sxmx,h=5/60,Weibull,lag=1)/cov5minlag1)-1)^(2)+
    w4*((pdrBLRPRX(a,l,v,k,f,h=5/60)/pdr5min)-1)^(2)
  S1h <- w1*((meanBLRPRX(l,k,f,i,h=1)/mean1)-1)^(2)+
    w2*((varBLRPRX(a,l,v,k,f,i,sxmx,h=1,Weibull)/var1)-1)^(2)+
    w3*((covarBLRPRX(a,l,v,k,f,i,sxmx,h=1,Weibull,lag=1)/cov1lag1)-1)^(2)+
    w4*((pdrBLRPRX(a,l,v,k,f,h=1)/pdr1)-1)^(2)
  S6h <- w1*((meanBLRPRX(l,k,f,i,h=6)/mean6)-1)^(2)+
    w2*((varBLRPRX(a,l,v,k,f,i,sxmx,h=6,Weibull)/var6)-1)^(2)+
    w3*((covarBLRPRX(a,l,v,k,f,i,sxmx,h=6,Weibull,lag=1)/cov6lag1)-1)^(2)+
    w4*((pdrBLRPRX(a,l,v,k,f,h=6)/pdr6)-1)^(2)
  S12h <- w1*((meanBLRPRX(l,k,f,i,h=12)/mean12)-1)^(2)+
    w2*((varBLRPRX(a,l,v,k,f,i,sxmx,h=12,Weibull)/var12)-1)^(2)+
    w3*((covarBLRPRX(a,l,v,k,f,i,sxmx,h=12,Weibull,lag=1)/cov12lag1)-1)^(2)+
    w4*((pdrBLRPRX(a,l,v,k,f,h=12)/pdr12)-1)^(2)
  S24h <- w1*((meanBLRPRX(l,k,f,i,h=24)/mean24)-1)^(2)+
    w2*((varBLRPRX(a,l,v,k,f,i,sxmx,h=24,Weibull)/var24)-1)^(2)+
    w3*((covarBLRPRX(a,l,v,k,f,i,sxmx,h=24,Weibull,lag=1)/cov24lag1)-1)^(2)+
    w4*((pdrBLRPRX(a,l,v,k,f,h=24)/pdr24)-1)^(2)
  S<-S5min+S1h+S24h
  if(is.infinite(S)) {S<-10^8}
  if(is.na(S)) {S<-10^8}
  return(S)
}

# In the case of EXPONENTIAL distribution for cell intensity, 
# set upper and lower bound of sxmx (last element) equal to 1
xlow<-c(2.0001,0.0001,0.001,0.001,0.001,0.001,1)
xup<-c(20,20,20,20,20,1,1)
# In the case of GAMMA or WEIBULL distribution for cell intensity:
xlow<-c(2.0001,0.0001,0.001,0.001,0.001,0.001,1)
xup<-c(20,20,20,20,20,1,5)

## Fitting the model with evalutionary annealing-simplex method

# Exponential or Gamma distribution for cell intensities
modecal<-eas(n=7,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=FALSE)
# Weibull distribution for cell intensities
modecal<-eas(n=7,m=30,xmin=xlow,xmax=xup,xlow=xlow,xup=xup,fn=fopt,maxeval=1500,ftol=1.e-07,
             ratio=0.99,pmut=0.9,beta=2,maxclimbs=5,Weibull=TRUE)

a<-modecal$bestpar[1];l<-modecal$bestpar[2];v<-modecal$bestpar[3]
k<-modecal$bestpar[4];f<-modecal$bestpar[5];i<-modecal$bestpar[6]
sxmx<-modecal$bestpar[7]

# To use the optimised set of parameters in the functions of HyetosMinute 
# please be sure that for parameters mx the length units are millimeters (mm) 
# and for parameters l and v the time units are days (d). 
# Due to this, the following unit conversions (hourly to daily) are needed:

l<-l*24;v<-v/24;

# The parameter set used in HyetosMinute functions will be

BLpar=list(lambda=l,phi=f,kappa=k,alpha=a,v=v,mx=NA,sxmx=sxmx)
# To implement Weibull distribution set: CellIntensityProp=list(Weibull=TRUE,iota=i)  
# To implement Exponential or Gamma distribution set: CellIntensityProp=list(Weibull=FALSE,iota=i)  


[Package HyetosMinute version 2.1 Index]