R 粒子群演算法 (PSO)

吃飽沒事做把戴老闆的Mathematica的PSO程式改成R

#Chen Yu-Ren PSO in R v1.0 2012/6/5#
#####
function (x)
100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
#####
mypso <-
function (fitfun)
{
    p <- 20
    phi1 <- 2
    phi2 <- 2
    w <- 0.5
    iter <- 300
    lrmax <- 5
    lrmin <- -5
    dims <- 2
    ltemp <- array(runif(dims * p, lrmin, lrmax), dim = c(p,
        dims), dimnames = list(paste("p", 1:p, sep = ""), paste("d",
        1:dims, sep = "")))
    vtemp <- array(runif(dims * p, lrmin/dims, lrmax/dims), dim = c(p,
        dims), dimnames = list(paste("p", 1:p, sep = ""), paste("v",
        1:dims, sep = "")))
    fitv <- apply(ltemp, 1, fitfun)
    pb <- ltemp
    dimnames(pb)[[2]] <- paste("pb", 1:dims, sep = "")
    outtemp <- cbind(ltemp, vtemp, pb, fitv)
    outtemp <- outtemp[order(outtemp[, 3 * dims + 1]), ]
#####
    for (i in 1:iter) {
        vtemp <- w * outtemp[, (dims + 1):(2 * dims)] + phi1 *
            runif(dims * p) * (outtemp[, (2 * dims + 1):(3 *
            dims)] - outtemp[, 1:dims]) + phi2 * runif(dims *
            p) * (array(rep(outtemp[1, 1:dims], each = p), dim = c(p,
            dims)) - outtemp[, 1:dims])
        vtemp <- pmin(pmax(vtemp, vrmin), vrmax)
        ltemp <- outtemp[, 1:dims] + vtemp
        ltemp <- pmin(pmax(ltemp, lrmin), lrmax)
        fitv <- apply(ltemp, 1, fitfun)
        pb <- ifelse(array(rep(fitv < outtemp[, (3 * dims + 1)],
            times = 2), dim = c(p, dims)), ltemp, outtemp[, 1:dims])
        colnames(pb) <- paste("pb", 1:dims, sep = "")
        outtemp <- cbind(ltemp, vtemp, pb, fitv)
        outtemp <- outtemp[order(outtemp[, 3 * dims + 1]), ]
        sdtemp <- sd(fitv)
        if (sdtemp < 1e-05)
            break
    }
#####
    list(par = round(outtemp[1, 1:dims], digits = 7), value = outtemp[1,
        (3 * dims + 1)], Iteration = i,
        `Standard Deviation of Fitness function` = sdtemp)
}
#####
mypso(fit)
optim(runif(2, -5, 5), fit, method = "Nelder-Mead")