From 238ffc40ef407e0c4f94b3a778cfead4b8a5e724 Mon Sep 17 00:00:00 2001
From: Mauricio Zambrano-Bigiarini <hzambran@users.noreply.github.com>
Date: Fri, 21 Sep 2012 16:18:14 +0000
Subject: [PATCH] 'hydroPSO': new 'normalise' argument for 'control'
---
NEWS | 14 +++++++++--
R/PSO_v2012.R | 67 ++++++++++++++++++++++++++++++++-----------------
man/hydroPSO.Rd | 23 ++++++++++-------
3 files changed, 70 insertions(+), 34 deletions(-)
diff --git a/NEWS b/NEWS
index 69ca1be..158d6a0 100755
--- a/NEWS
+++ b/NEWS
@@ -1,8 +1,14 @@
NEWS/ChangeLog for hydroPSO
--------------------------
0.1-59 (under-development & not working yet)
- o 'hydroPSO' : -) towards SPSO 2011 capable...
+ o 'hydroPSO' : -) towards SPSO 2011 capable, which includes:
+ * new default value for swarm size (npart=40),
+ * new equation for updating the velocity ( V[t+1] = w*V + Gr -x + r*alea_sphere(Gr, ||Gr-x||) ),
+ * new initialization of the velocity ( Vini = U(lower-Xini, upper-Xini) )
+ * optional normalisation of parameter values
+ -) new 'normalise' parameter for the 'control' variable, in order to improve the performance when the search space is not
+ a hypercube
-) argument 'method' now allows the following 2 new values: 'spso2007' and 'spso2011', which set the values of the PSO
engine to the ones of the SPSO 2007 and SPSO 2011, respectively.
-) argument 'method'. The old (default) value 'pso' was replaced by the new (default) value 'spso2011'
@@ -16,8 +22,12 @@ NEWS/ChangeLog for hydroPSO
method == 'spso2007' => npart=10+2*[sqrt(n)]
method != 'spso2007' => npart=40
-) in the documentation, default values are now mentioned for each argument of 'control'.
+ -) the 'PSO_logfile.txt' now contains information about: 'best.update', 'random.update', 'Xini.type', 'Vini.type'.
+ In addition, information about 'IW.type' and 'IW.exp' is now written only when 'length(IW.w) > 1'
- o Results obtained when running 'hydroPSO' >= 0.1-59 with default settings WILL BE DIFFERENT from those obtained with 'hydroPSO' <= 0.1-58, due to the following changes in default values:
+ o Running 'hydroPSO' >= 0.1-59 with default settings will produce DIFFERENT RESULTS from those obtained with 'hydroPSO' <= 0.1-58, due to the
+ following changes in default values:
+ ver <= 0.1-58 -> ver >= 0.1-59
-) npart : 10+2*[sqrt(n)] -> 40
-) Vini.type : 'lhs2007' -> 'random2011'
-) boundary.wall: 'reflecting' -> 'absorbing'
diff --git a/R/PSO_v2012.R b/R/PSO_v2012.R
index 938409d..757ac78 100755
--- a/R/PSO_v2012.R
+++ b/R/PSO_v2012.R
@@ -537,7 +537,7 @@ position.update.and.boundary.treatment <- function(x, v, x.MinMax, boundary.wall
if ( length(byd.min.pos) > 0) {
if ( boundary.wall == "absorbing") {
x.new[byd.min.pos] <- x.min[byd.min.pos]
- v.new[byd.min.pos] <- 0*v[byd.min.pos]
+ v.new[byd.min.pos] <- -0.0*v[byd.min.pos]
} else if ( boundary.wall == "reflecting") {
x.new[byd.min.pos] <- 2*x.min[byd.min.pos] - x.new[byd.min.pos]
v.new[byd.min.pos] <- v[byd.min.pos]
@@ -555,7 +555,7 @@ position.update.and.boundary.treatment <- function(x, v, x.MinMax, boundary.wall
if ( length(byd.max.pos) > 0 ) {
if ( boundary.wall == "absorbing") {
x.new[byd.max.pos] <- x.max[byd.max.pos]
- v.new[byd.max.pos] <- 0*v[byd.max.pos]
+ v.new[byd.max.pos] <- -0.0*v[byd.max.pos]
} else if ( boundary.wall == "reflecting") {
x.new[byd.max.pos] <- 2*x.max[byd.max.pos] - x.new[byd.max.pos]
v.new[byd.max.pos] <- v[byd.max.pos]
@@ -1559,6 +1559,7 @@ hydroPSO <- function(
maxfn=Inf,
c1= 0.5+log(2),
c2= 0.5+log(2),
+ use.IW= TRUE, IW.w=1/(2*log(2)), IW.type=c("linear", "non-linear", "runif", "aiwf", "GLratio"), IW.exp= 1,
use.CF= FALSE,
lambda= 1,
@@ -1572,8 +1573,9 @@ hydroPSO <- function(
topology=c("random", "gbest", "lbest", "vonNeumann"), K=3,
iter.ini=0, # only used when 'topology=lbest'
ngbest=4, # only used when 'method=ipso'
+
+ normalise=FALSE,
- use.IW = TRUE, IW.w=1/(2*log(2)), IW.type=c("linear", "non-linear", "runif", "aiwf", "GLratio"), IW.exp= 1,
use.TVc1= FALSE, TVc1.rng= c(1.28, 1.05), TVc1.type=c("linear", "non-linear", "GLratio"), TVc1.exp= 1.5,
use.TVc2= FALSE, TVc2.rng= c(1.05, 1.28), TVc2.type=c("linear", "non-linear"), TVc2.exp= 1.5,
use.TVlambda=FALSE, TVlambda.rng= c(1, 0.25), TVlambda.type=c("linear", "non-linear"), TVlambda.exp= 1,
@@ -1623,6 +1625,9 @@ hydroPSO <- function(
maxfn <- con[["maxfn"]]
c1 <- con[["c1"]]
c2 <- con[["c2"]]
+ use.IW <- as.logical(con[["use.IW"]])
+ IW.w <- con[["IW.w"]]
+ IW.exp <- con[["IW.exp"]]
use.CF <- con[["use.CF"]]
lambda <- con[["lambda"]]
abstol <- con[["abstol"]]
@@ -1630,10 +1635,8 @@ hydroPSO <- function(
random.update <- as.logical(con[["random.update"]])
K <- con[["K"]]
iter.ini <- con[["iter.ini"]]
- ngbest <- con[["ngbest"]]
- use.IW <- as.logical(con[["use.IW"]])
- IW.w <- con[["IW.w"]]
- IW.exp <- con[["IW.exp"]]
+ ngbest <- con[["ngbest"]]
+ normalise <- as.logical(con[["normalise"]])
use.TVc1 <- as.logical(con[["use.TVc1"]])
TVc1.rng <- con[["TVc1.rng"]]
TVc1.exp <- con[["TVc1.exp"]]
@@ -1744,6 +1747,14 @@ hydroPSO <- function(
} # ELSE end
n <- nrow(X.Boundaries)
+
+ if (normalise) {
+ X.Boundaries[,1] <- rep(0, n)
+ X.Boundaries[,2] <- rep(1, n)
+
+ lower.mat <- matrix( rep(X.Boundaries[,1], npart), nrow=npart, byrow=TRUE)
+ upper.mat <- matrix( rep(X.Boundaries[,2], npart), nrow=npart, byrow=TRUE)
+ } # IF end
if (is.null(rownames(X.Boundaries))) {
param.IDs <- paste("Param", 1:n, sep="")
@@ -1890,7 +1901,7 @@ hydroPSO <- function(
} else if ( (class(par)=="matrix") | (class(par)=="data.frame") ) {
tmp <- ncol(par)
if ( tmp != n )
- stop( paste("Invalid argument: 'ncol(par) != n' (",tmp, "!=", n, ")", sep="") )
+ stop( "Invalid argument: 'ncol(par) != n' (",tmp, "!=", n, ")" )
tmp <- nrow(par)
X[1:tmp,] <- par
} # ELSE end
@@ -2034,6 +2045,18 @@ hydroPSO <- function(
writeLines(c("c2 :", c2), PSOparam.TextFile, sep=" ")
writeLines("", PSOparam.TextFile)
} # ELSE end
+ writeLines(c("use.IW :", use.IW), PSOparam.TextFile, sep=" ")
+ writeLines("", PSOparam.TextFile)
+ if (use.IW) {
+ writeLines(c("IW.w :", IW.w), PSOparam.TextFile, sep=" ")
+ writeLines("", PSOparam.TextFile)
+ if ( length(IW.w) > 1 ) {
+ writeLines(c("IW.type :", IW.type), PSOparam.TextFile, sep=" ")
+ writeLines("", PSOparam.TextFile)
+ writeLines(c("IW.exp :", IW.exp), PSOparam.TextFile, sep=" ")
+ writeLines("", PSOparam.TextFile)
+ } # IF end
+ } # IF end
if (use.TVlambda) {
writeLines(c("use.TVlambda :", use.TVlambda), PSOparam.TextFile, sep=" ")
writeLines("", PSOparam.TextFile)
@@ -2047,16 +2070,6 @@ hydroPSO <- function(
writeLines(c("lambda :", lambda), PSOparam.TextFile, sep=" ")
writeLines("", PSOparam.TextFile)
} # ELSE end
- if (use.IW) {
- writeLines(c("use.IW :", use.IW), PSOparam.TextFile, sep=" ")
- writeLines("", PSOparam.TextFile)
- writeLines(c("IW.type :", IW.type), PSOparam.TextFile, sep=" ")
- writeLines("", PSOparam.TextFile)
- writeLines(c("IW.w :", IW.w), PSOparam.TextFile, sep=" ")
- writeLines("", PSOparam.TextFile)
- writeLines(c("IW.exp :", IW.exp), PSOparam.TextFile, sep=" ")
- writeLines("", PSOparam.TextFile)
- } # IF end
writeLines(c("maxfn :", maxfn), PSOparam.TextFile, sep=" ")
writeLines("", PSOparam.TextFile)
writeLines(c("abstol :", abstol), PSOparam.TextFile, sep=" ")
@@ -2284,9 +2297,12 @@ hydroPSO <- function(
} # IF end
##########################################################################
+
+ if (normalise) X <- X * (upper.mat - lower.mat) + lower.mat
+
# 3.a) Evaluate the particles fitness
if ( fn.name != "hydromod" ) {
-
+
# Evaluating an R Function
GoF <- apply(X, fn, MARGIN=1, ...)
@@ -2322,7 +2338,6 @@ hydroPSO <- function(
} # ELSE end
-
if ( best.update == "sync" ) {
tmp <- sync.update.pgbests(x=X,
xt.fitness= Xt.fitness[iter, ],
@@ -2574,7 +2589,7 @@ hydroPSO <- function(
gbest.fit.bak <- ngbest.fit
} # IF end
- if (verbose) message("[Re-grouping particles in the swarm (iter: ", iter, ") ...]")
+ if (verbose) message("[ Re-grouping particles in the swarm (iter: ", iter, ") ... ]")
tmp <- RegroupingSwarm(x=X,
gbest= X.best.part[gbest.pos, ],
@@ -2665,13 +2680,17 @@ hydroPSO <- function(
# File 'BestParamPerIter.txt' #
GoF <- gbest.fit
if(is.finite(GoF)) {
+
+ ifelse(normalise, temp <- X.best.part[gbest.pos, ] * (upper - lower) + lower,
+ temp <- X.best.part[gbest.pos, ] )
+
writeLines( as.character( c(iter,
formatC(GoF, format="E", digits=digits, flag=" "),
- formatC(X.best.part[gbest.pos, ], format="E", digits=digits, flag=" ")
+ formatC(temp, format="E", digits=digits, flag=" ")
) ), BestParamPerIter.TextFile, sep=" ")
} else writeLines( as.character( c(iter,
"NA",
- formatC(X.best.part[gbest.pos, ], format="E", digits=digits, flag=" ")
+ formatC(temp, format="E", digits=digits, flag=" ")
) ), BestParamPerIter.TextFile, sep=" ")
writeLines("", BestParamPerIter.TextFile)
flush(BestParamPerIter.TextFile)
@@ -2700,6 +2719,8 @@ hydroPSO <- function(
} # WHILE end
########################## END Main Iteration Loop #########################
+
+ if (normalise) X.best.part <- X.best.part * (upper.mat - lower.mat) + lower.mat
if (write2disk) {
close(OFout.Text.file)
diff --git a/man/hydroPSO.Rd b/man/hydroPSO.Rd
index e05f60e..9df5adc 100755
--- a/man/hydroPSO.Rd
+++ b/man/hydroPSO.Rd
@@ -107,6 +107,15 @@ By default \code{c1= 0.5 + log(2)}
\item{c2}{
numeric, social acceleration coefficient. Encourages the exploration of the current global best and reflects how much the particle is influenced by the best-known optimum of the swarm \cr
By default \code{c2= 0.5 + log(2)}
+}
+ \item{use.IW}{
+logical, indicates if an inertia weight (\env{w}) will be used to avoid swarm explosion, i.e. particles flying around their best position without converging into it (see Shi and Eberhart, 1998) \cr
+By default \code{use.IW=TRUE}
+}
+ \item{IW.w}{
+OPTIONAL. Used only when \code{use.IW= TRUE \& IW.type!='GLratio'} \cr
+numeric, value of the inertia weight(s) (\env{w} or \env{[w.ini, w.fin]}). It can be a single number which is used for all iterations, or it can be a vector of length 2 with the initial and final values (in that order) that \env{w} will take along the iterations \cr
+By default \code{IW.w=1/(2*log(2))}
}
\item{use.CF}{
logical, indicates if the Clerc's Constriction Factor (see Clerc, 1999; Eberhart and Shi, 2000; Clerc and Kennedy, 2002) is used to avoid swarm explosion \cr
@@ -180,12 +189,13 @@ OPTIONAL. Only used when \code{method=='ipso'} \cr
numeric, number of particles considered in the computation of the global best \cr
By default \code{ngbest=4} (see Zhao, 2006)
}
- \item{use.IW}{
-logical, indicates if an inertia weight (\env{w}) will be used to avoid swarm explosion, i.e. particles flying around their best position without converging into it (see Shi and Eberhart, 1998) \cr
-By default \code{use.IW=TRUE}
+ \item{randomise}{
+logical, indicates whether the parameter values have to be normalised to the [0,1] interval during the optimisation or not\cr
+This option appears in the C and Matlab version of SPSO-2011 (See \url{http://www.particleswarm.info/standard_pso_2011_c.zip}) and there it is recommended to use this option when the search space is not an hypercube. If the search space is an hypercube, it is better not normalise (there is a small difference between the position without any normalisation and the de-normalised one)\cr
+By default \code{randomise=FALSE}
}
\item{IW.type}{
-OPTIONAL. Used only when \code{use.IW= TRUE} \cr
+OPTIONAL. Used only when \code{use.IW= TRUE} \& \code{length(IW.w)>1}\cr
character, defines how the inertia weight \env{w} will vary along iterations. Valid values are: \cr
-)\kbd{linear}: \env{w} varies linearly between the initial and final values specified in \code{IW.w} (see Shi and Eberhart, 1998; Zheng et al., 2003). This is the DEFAULT option \cr
-)\kbd{non-linear}: \env{w} varies non-linearly between the initial and final values specified in \code{IW.w} with exponential factor \kbd{IW.exp} (see Chatterjee and Siarry, 2006) \cr
@@ -193,11 +203,6 @@ character, defines how the inertia weight \env{w} will vary along iterations. Va
-)\kbd{aiwf}: adaptive inertia weight factor, where the inertia weight is varied adaptively depending on the goodness-of-fit values of the particles (see Liu et al., 2005) \cr
-)\kbd{GLratio}: \env{w} varies according to the ratio between the global best and the average of the particle's local best (see Arumugam and Rao, 2008) \cr
By default \code{IW.type='linear'}
-}
- \item{IW.w}{
-OPTIONAL. Used only when \code{use.IW= TRUE \& IW.type!='GLratio'} \cr
-numeric, value of the inertia weight(s) (\env{w} or \env{[w.ini, w.fin]}). It can be a single number which is used for all iterations, or it can be a vector of length 2 with the initial and final values (in that order) that \env{w} will take along the iterations \cr
-By default \code{IW.w=1/(2*log(2))}
}
\item{IW.exp}{
OPTIONAL. Used only when \code{use.IW= TRUE} AND \code{IW.type= 'non-linear'} \cr
--
GitLab