eas {HyetosMinute} | R Documentation |
An enhanced version of the evolutionary annealing-simplex optimization method for the estimation of Bartlett-Lewis model parameters.
eas(n,m,xmin,xmax,xlow,xup,fn,maxeval=1500,ftol=1.e-07,ratio=0.99,pmut=0.9, beta=2,maxclimbs=5,...)
n |
An integer that specifies the problem dimension (number of control variables). |
m |
A positive integer, greater than |
xmin |
A vector, with length equal to |
xmax |
A vector, with length equal to |
xlow |
A vector, with length equal to |
xup |
A vector, with length equal to |
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 |
ftol |
A positive number that specifies the fractional convergence tolerance to be achieved in the function value. Default is |
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 |
pmut |
A positive number, between 0.5-0.95, that specifies the probability of accepting an offspring generated via mutation. Default is |
beta |
A positive integer,greater than 1, that specifies the annealing schedule parameter. Default is |
maxclimbs |
A positive integer, typically between 3-5, that specifies the maximum number of uphill steps. Default is |
... |
Further arguments to be passed in the objective function. |
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.
A list with the following components:
bestpar |
A vector containing the best set of parameters found. |
bestval |
The value of |
nfeval |
Number of function |
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. |
Kossieris Panagiotis pankoss@hotmail.com, with Hristos Tyralis montchrister@gmail.com and Andreas Efstratiadis A.Efstratiadis@itia.ntua.gr.
The methodology of the evolutionary annealing-simplex optimization method and details of its application are described in:
Efstratiadis, A., and D. Koutsoyiannis, An evolutionary annealing-simplex algorithm for global optimisation of water resource systems, Proceedings of the Fifth International Conference on Hydroinformatics, Cardiff, UK, 1423-1428, International Water Association, 2002. (http://itia.ntua.gr/en/docinfo/524/)
Rozos, E., A. Efstratiadis, I. Nalbantis, and D. Koutsoyiannis, Calibration of a semi-distributed model for conjunctive simulation of surface and groundwater flows, Hydrological Sciences Journal, 49 (5), 819-842, 2004. (http://itia.ntua.gr/en/docinfo/630/)
Find the original Pascal-Delphi code of optimisation algorithms on http://itia.ntua.gr/en/docinfo/524/.
# 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)