R/globpso.R

In globpso: Swarm Intelligence Optimization

#' Particle Swarm Optimization Algorithms for Minimization Problems
#'
#' The general-purpose implementation of particle swarm Optimization algorithms for minimizing an user-defined objective function.
#'
#' @param objFunc The R or Rcpp compiled objective function. See the example.
#' @param lower The vector of finite lower bounds of the search domain.
#' Its length should be equal to the dimension of the domain space.
#' @param upper The vector of finite upper bounds of the search domain.
#' Its length should be equal to the dimension of the domain space.
#' @param init The vector of initial swarm.
#' Its length should be equal to the dimension of the domain space.
#' When there are more than one initial vectors, specify \code{init} as
#' a matrix. Each row vector represents one initial point.
#' The default for \code{init} is \code{NULL}.
#' @param fixed The vector of real values and NA values that controls PSO to search only 
#' for the NA-valued components. 
#' @param PSO_INFO The list of PSO parameters generated by \code{getPSOInfo()}.
#' @param seed The random seed that controls initial swarm of PSO. The default is \code{NULL}.
#' @param verbose The logical value controls if PSO would reports the updating progress. The default is \code{TRUE}.
#' @param ... Further arguments to be passed to \code{objFunc}.
#' @return An List.
#' \describe{
#' \item{par}{ the global best particle.}
#' \item{val}{ the objective function value of the global best particle.}
#' \item{history}{ a vector of objective function values of the global best particle in PSO search history.}
#' \item{cputime}{ the computational time in seconds.}
#' }
# @details
# TBD
#' @examples
#' # three-dimensional function
#' objf <- function(x, loc) {
#'   val <- 0
#'   for (i in 1:length(x)) {
#'     val <- val + (x[i] - loc)^2
#'   }
#'   return(val)
#' }
#' 
#' upp_bound <- rep(5, 3)
#' low_bound <- rep(-5, 3)
#' loc_shift <- 1
#'
#' alg_setting <- getPSOInfo(nSwarm = 32, maxIter = 100, psoType = "basic")
#' res <- globpso(objFunc = objf, lower = low_bound, upper = upp_bound, 
#'                PSO_INFO = alg_setting, loc = loc_shift)
#' res$par
#' res$val
#' 
#' # C++ function example
#' \donttest{
#' library(Rcpp)
#' library(RcppArmadillo)
#' objf_c <- cppFunction('double objf_c(SEXP x, SEXP loc) {
#'     double val = 0;
#'     double loc_c = (double)Rcpp::as<double>(loc);
#'     arma::rowvec x_c = (arma::rowvec)Rcpp::as<arma::rowvec>(x);
#'     for (arma::uword i = 0; i < x_c.n_elem; i++) {
#'       val += (x_c(i) - loc_c)*(x_c(i) - loc_c);
#'     }
#'     return val;
#'   }', depends = "RcppArmadillo")
#'   
#' alg_setting <- getPSOInfo(nSwarm = 32, maxIter = 100, psoType = "quantum")
#' res_c <- globpso(objFunc = objf_c, lower = low_bound, upper = upp_bound, 
#'                  PSO_INFO = alg_setting, loc = 1)
#' res_c$par
#' res_c$val
#' }
#' @references 
#' \enumerate{
#'   \item Bonyadi, M. R., & Michalewicz, Z. (2014). A locally convergent rotationally invariant particle swarm optimization algorithm. Swarm intelligence, 8(3), 159-198.
#'   \item Cheng, R., & Jin, Y. (2014). A competitive swarm optimizer for large scale optimization. IEEE transactions on cybernetics, 45(2), 191-204.
#   \item Eberhart, R. & Kennedy, J. (1995). A new optimizer using particle swarm theory. In The 6th International Symposium on Micro Machine and Human Science, pages 39-43. IEEE.
#'   \item Shi, Y., & Eberhart, R. (1998, May). A modified particle swarm optimizer. In Evolutionary Computation Proceedings, 1998. IEEE World Congress on Computational Intelligence., The 1998 IEEE International Conference on (pp. 69-73). IEEE.
#'   \item Stehlík, M., Chen, P. Y., Wong, W. K., and Kiseľák, J. (2024). A double exponential particle swarm optimization with non-uniform variates as stochastic tuning and guaranteed convergence to a global optimum with sample applications to finding optimal exact designs in biostatistics. Applied Soft Computing, 163, 111913.
#'   \item Sun, J., Feng, B., and Xu, W. (2004a). Particle swarm optimization with particles having quantum behavior. In Evolutionary Computation, 2004. CEC2004. Congress on, volume 1, pages 325-331. IEEE.
#' }
#' @name globpso
#' @rdname globpso
#' @importFrom methods hasArg
#' @importFrom Rcpp evalCpp cppFunction sourceCpp
#' @export
globpso <- function(objFunc, lower, upper, init = NULL, fixed = NULL, 
	PSO_INFO = NULL, seed = NULL, verbose = TRUE, ...) {

  stopifnot(all(is.finite(lower)), all(is.finite(upper)), 
            length(lower) == length(upper), all(upper >= lower)
            )
  hasInitSwarm <- 0
  if (!is.null(init)) {
    if (is.vector(init)) {
      stopifnot(length(lower) == length(init), all(upper >= init), 
                all(is.finite(init)), all(init >= lower))
      init <- matrix(init, 1, length(init))
      hasInitSwarm <- 1
    } else if (is.matrix(init)) {
      for (i in 1:nrow(init)) {
        stopifnot(length(lower) == length(init[i,]), all(upper >= init[i,]), 
                  all(is.finite(init[i,])), all(init[i,] >= lower))
      }
      hasInitSwarm <- 1
    }
  }
  if (is.null(fixed)) { 
    fixed <- rep(NA, length(lower))
  } else {
    stopifnot(length(fixed) == length(lower))
  }
  
	if (is.null(PSO_INFO)) {
	  nSwarm <- 32
	  if (hasInitSwarm) { nSwarm = max(32, nrow(init)) }
		PSO_INFO <- getPSOInfo(nSwarm = nSwarm, maxIter = 100)
		if (verbose) message(paste0("Use the default settings for PSO. See '?getPSOInfo'."))
	}
	stopifnot(all(names(PSO_INFO) == names(getPSOInfo())))
	if (hasInitSwarm) { stopifnot(PSO_INFO$nSwarm >= nrow(init)) }
	
	environment <- new.env()

	# Fill the rest of PSO options in PSO_INFO
	PSO_INFO$varUpper <- upper
	PSO_INFO$varLower <- lower
	PSO_INFO$dSwarm <- length(upper)
	PSO_INFO$hasInitSwarm <- hasInitSwarm
	PSO_INFO$initSwarm <- init
	PSO_INFO$fixedDims <- fixed
	# Start
	set.seed(seed)
	cputime <- system.time(
		psoOut <- cppPSO(objFunc, PSO_INFO, environment, FALSE, verbose)
	)[3]
	if (verbose) message(paste0("CPU time: ", round(cputime, 2), " seconds."))

	return(
		list(par = psoOut$GBest, val = psoOut$fGBest, 
		     pbest = psoOut$PBest,
		     fpbest = psoOut$fPBest,
		     history = psoOut$fGBestHist, cputime = cputime)
	)
}


#' Generation function of PSO parameter settings
#' 
#' Create a list with PSO parameters for Minimization.
#' 
#' @param nSwarm A integer number of swarm size in PSO algorithm.
#' @param maxIter A integer number of maximal PSO iterations.
# @param checkConv A logical value which controls whether PSO checks the stopping criterion during updating procedure.
# Specify \code{TRUE} for PSO to compute the stopping criterion \eqn{|f'-f|<\varepsilon}
# where \eqn{f'} and \eqn{f} are the objective function values in the previous and current iterations, respectively.
# The default is \code{FALSE}.
#' @param freeRun A number between \eqn{[0,1]} that controls the percentage of PSO iterations which are free from examining the
#' stopping criterion, \eqn{|f'-f|<\varepsilon}
#' where \eqn{f'} and \eqn{f} are the objective function values in the previous and current iterations, respectively.
#' The default is 1.0 implying the PSO will completely ignore the stopping criterion.
#' Otherwise, the PSO checks the stopping criterion after free iterations.
#' @param tol A small value for the tolerance, \eqn{\varepsilon}, in the stopping criterion.
#' For \code{freeRun} smaller than 1.0, the default is \code{1e-6}. Otherwise, this value would not affect the algorithm.
#' @param psoType string. The type of PSO. This package current supports the following types:
#' \describe{
#' \item{"basic"}{ Linearly Decreasing Weight PSO (Eberhart & Kennedy, 1995)}
# \item{1}{ GCPSO (van den Bergh, F. &	Engelbrecht, A. P., 2002)}
#' \item{"quantum"}{ Quantum PSO (Sun et al., 2004)}
#' \item{"lcri"}{ LcRiPSO (Bonyadi & Michalewicz, 2014)}
#' \item{"comp"}{ Competitive Swarm Optimization (Cheng & Jin, 2014)}
#' \item{"dexp"}{ DExPSO (Stehlík et al., 2024)}
#' }
#' @param c1 The value of cognitive parameter in PSO updating procedure. The default is 2.05.
#' @param c2 The value of social parameter in PSO updating procedure. The default is 2.05.
#' @param w0 The value of starting inertia weight in PSO updating procedure. The default is 1.2.
#' @param w1 The value of ending inertia weight in PSO updating procedure. The default is 0.2.
#' @param w_var A number between \eqn{[0,1]} that controls the percentage of iterations that PSO linearly decreases the inertia weight
#' from \code{w0} to \code{w1}. The default is 0.8.
#' @param vk The value of velocity clamping parameter. The default is 4.
#' @param Q_cen_type The type of the center position in QPSO updating procedure (\code{0}: local attractor, default; \code{1}: mean best).
#' @param Q_a0 The value of starting contraction-expansion (CE) coefficient in QPSO updating procedure. The default is 1.7.
#' @param Q_a1 The value of ending contraction-expansion (CE) coefficient in QPSO updating procedure. The default is 0.7.
#' @param Q_a_var A number between \eqn{[0,1]} that controls the percentage of iterations that QPSO linearly decreases the CE coefficient
#' from \code{Q_a0} to \code{Q_a1}. The default is 0.8.
#' @param LcRi_L The value of random number generator based on normal density social parameter in LcRiPSO updating procedure. The default is 0.01. (for \code{psoType = c("lcri")} only)
#' @param CSO_phi The value of social parameter in CSO updating procedure. The default is 0.1. (for \code{psoType = c("comp", "cdexp")} only)
#' @param TE_b The value of random number generator based on double-exponential density. The default is 2.0. (for \code{psoType = c("dexp", "qdexp", "cdexp")} only).
#' @return A list of PSO parameter settings.
#' @examples
#' PSO_INFO <- getPSOInfo(nSwarm = 32, maxIter = 100)
#' @name getPSOInfo
#' @rdname getPSOInfo
#' @export
getPSOInfo <- function(nSwarm = 32, maxIter = 100, psoType = "basic",
  #dSwarm = NULL, varUpper = NULL, varLower = NULL, checkConv = 0,
  freeRun = 1.0, tol = 1e-6, c1 = 2.05, c2 = 2.05,
  w0 = 1.2, w1 = 0.2, w_var = 0.8, vk = 4, #, chi = NULL,
  #GC_S_ROOF = 5, GC_F_ROOF = 15, GC_RHO = 1,
  Q_cen_type = 1, Q_a0 = 1.7, Q_a1 = 0.7, Q_a_var = 0.8, LcRi_L = 0.01,
  CSO_phi = 0.1, TE_b = 2
  ) {

  stopifnot(length(nSwarm) == 1)
  if (!(psoType %in% c("basic", "quantum", "lcri", "comp", "dexp", "qdexp", "cdexp"))) {
		stop("Currently the function supports: \npsoType = 'basic' for basic PSO algorithm.\npsoType = 'quantum' for quantum PSO algorithm. \npsoType = 'lcri' for LcRi PSO algorithm. \npsoType = 'comp' for Competitive Swarm Optimizer algorithm. \npsoType = 'dexp' for DExPSO algorithm.   
		")
	} else {
	  if (psoType == "basic") {
	    typePSO <- 0
	  } else if (psoType == "quantum") {
	    typePSO <- 2
	  } else if (psoType == "lcri") {
	    typePSO <- 4
	  } else if (psoType == "comp") {
	    typePSO <- 5
	  } else if (psoType == "dexp") {
	    typePSO <- 2024
	  } else if (psoType == "qdexp") {
	    typePSO <- 20241
	  } else if (psoType == "cdexp") {
	    typePSO <- 20242
	  }
	}

  nLoop <- length(nSwarm)
  #if (length(dSwarm))     dSwarm     <- numeric(nLoop)
  #if (length(varUpper))   varUpper   <- matrix(0, nLoop)
  #if (length(varLower))   varLower   <- matrix(0, nLoop)
  #if (length(maxIter))    maxIter    <- rep(100   , nLoop)
  #if (length(checkConv) < nLoop)  checkConv  <- rep(checkConv, nLoop)
  if (length(typePSO) < nLoop)    typePSO    <- rep(typePSO, nLoop)
  if (length(freeRun) < nLoop)    freeRun    <- rep(freeRun, nLoop)
  if (length(tol) < nLoop)        tol        <- rep(tol, nLoop)
  if (length(c1) < nLoop)         c1         <- rep(c1, nLoop)
  if (length(c2) < nLoop)         c2         <- rep(c2, nLoop)
  if (length(w0) < nLoop)         w0         <- rep(w0, nLoop)
  if (length(w1) < nLoop)         w1         <- rep(w1, nLoop)
  if (length(w_var) < nLoop)      w_var      <- rep(w_var, nLoop)
  #if (length(chi) < nLoop)        chi        <- rep(0.729 , nLoop)
  if (length(vk) < nLoop)         vk         <- rep(vk, nLoop)
  #if (length(GC_S_ROOF) < nLoop)  GC_S_ROOF  <- rep(5     , nLoop)
  #if (length(GC_F_ROOF) < nLoop)  GC_F_ROOF  <- rep(15    , nLoop)
  #if (length(GC_RHO) < nLoop)     GC_RHO     <- rep(1     , nLoop)
  if (length(Q_cen_type) < nLoop) Q_cen_type <- rep(Q_cen_type, nLoop)
  if (length(Q_a0) < nLoop)       Q_a0       <- rep(Q_a0, nLoop)
  if (length(Q_a1) < nLoop)       Q_a1       <- rep(Q_a1, nLoop)
  if (length(Q_a_var) < nLoop)    Q_a_var    <- rep(Q_a_var, nLoop)
  if (length(LcRi_L) < nLoop)     LcRi_L     <- rep(0.01  , nLoop)
  if (length(CSO_phi) < nLoop)    CSO_phi    <- rep(0.1   , nLoop)
  if (length(TE_b) < nLoop)       TE_b       <- rep(2.    , nLoop)

  list(
    nSwarm = nSwarm, dSwarm = "autogen", varUpper = "autogen", varLower = "autogen", initSwarm = "autogen", 
    hasInitSwarm = "autogen", fixedDims = "autogen", maxIter = maxIter, 
  	typePSO = typePSO, #checkConv = checkConv,
    freeRun = freeRun, tol = tol, c1 = c1, c2 = c2, w0 = w0, w1 = w1, w_var = w_var, #chi = chi,
    vk = vk,  #GC_S_ROOF = GC_S_ROOF, GC_F_ROOF = GC_F_ROOF,
    #GC_RHO = GC_RHO, 
    Q_cen_type = Q_cen_type, Q_a0 = Q_a0, Q_a1 = Q_a1, Q_a_var = Q_a_var,
    LcRi_L = LcRi_L, CSO_phi = CSO_phi,  TE_b = TE_b
  )
}