R/trimton.R

In windfarmGA: Genetic Algorithm for Wind Farm Layout Optimization

#' @title Adjust the amount of turbines per windfarm
#' @name trimton
#' @description Legacy repair for binary chromosomes. The GA loop encodes
#'   layouts as \code{n} unique grid IDs, so this function is not called there.
#'   It remains exported for the old 0/1 pipeline
#'   (\code{\link{crossover}} / \code{\link{mutation}}).
#'
#' @export
#'
#' @param mut A binary matrix with the mutated individuals
#' @param nturb A numeric value indicating the amount of required turbines
#' @param allparks A data.frame consisting of all individuals of the current
#'   generation
#' @param nGrids A numeric value indicating the total amount of grid cells
#' @param trimForce If `TRUE`, add or drop turbines using fitness-weighted
#'   probabilities. If `FALSE`, choose cells at random.
#' @param seed Set a seed for comparability. Default is NULL
#'
#' @family Genetic Algorithm Functions
#' @return Returns a binary matrix with the correct amount of turbines per
#'   individual
#' @examples \donttest{
#' ## Create a random rectangular shapefile
#' library(sf)
#' area <- sf::st_as_sf(sf::st_sfc(
#'   sf::st_polygon(list(cbind(
#'     c(0, 0, 2000, 2000, 0),
#'     c(0, 2000, 2000, 0, 0)
#'   ))),
#'   crs = 3035
#' ))
#'
#' ## Create a uniform and unidirectional wind data.frame and plots the
#' ## resulting wind rose
#' ## Uniform wind speed and single wind direction
#' data.in <- as.data.frame(cbind(ws = 12, wd = 0))
#'
#' ## Calculate a Grid and an indexed data.frame with coordinates and grid cell Ids.
#' Grid1 <- grid_area(area = area, size = 200, prop = 1)
#' Grid <- Grid1[[1]]
#' AmountGrids <- nrow(Grid)
#'
#' startsel <- init_population(Grid, 10, 20)
#' wind <- as.data.frame(cbind(ws = 12, wd = 0))
#' wind <- list(wind, probab = 100)
#' fit <- fitness(
#'   population = startsel, reference_height = 100, rotor_height = 100,
#'   surface_roughness = 0.3, area = area, rotor = 20,
#'   wind = wind, terrain = FALSE
#' )
#' allparks <- do.call("rbind", fit)
#' ## selection() returns ID matrices; crossover()/trimton() expect 0/1.
#' sel <- selection(fit, Grid, 2, TRUE, 6, "FIX")
#' ids <- sel[[1]]
#' bins <- matrix(0, nrow(Grid), ncol(ids))
#' for (j in seq_len(ncol(ids))) {
#'   bins[match(ids[, j], Grid[, "ID"]), j] <- 1
#' }
#' selec6best <- list(
#'   data.frame(ID = Grid[, "ID"], bins),
#'   data.frame(ID = 1, t(sel[[2]]))
#' )
#' crossOut <- crossover(selec6best, 2, uplimit = 300, crossPart = "RAN")
#' mut <- mutation(a = crossOut, p = 0.3, NULL)
#' mut1 <- trimton(
#'   mut = mut, nturb = 10, allparks = allparks, nGrids = AmountGrids,
#'   trimForce = FALSE
#' )
#' colSums(mut)
#' colSums(mut1)
#' }
trimton <- function(mut, nturb, allparks, nGrids, trimForce, seed) {
  if (missing(seed)) {
    seed <- NULL
  }
  k <- 0.5
  nGrids1 <- 1:nGrids

  # Probability that a turbine is selected to be eliminated / added.
  indivprop <- subset.matrix(allparks, select = c("Rect_ID", "Parkfitness", "AbschGesamt"))
  indivprop <- aggregate(indivprop[, 2:3], by = list(indivprop[, 1]), FUN = mean)
  colnames(indivprop) <- c("Rect_ID", "Parkfitness", "AbschGesamt")

  lepa <- length(mut[1, ])
  mut1 <- vector("list", lepa)
  for (i in 1:lepa) {
    tmp <- mut[, i]
    e <- tmp == 1
    ## How much turbines are there too many?
    zviel <- sum(e) - nturb
    ## Which grid cell IDs have a turbine
    welche <- which(e)

    propwelche <- cbind(
      RectID = welche,
      Prop = rep(mean(indivprop[, "AbschGesamt"]), length(welche))
    )
    if (trimForce) {
      propexi <- indivprop[indivprop[, "Rect_ID"] %in% welche, ]
      npt <- (1 + ((max(propexi[, "AbschGesamt"]) - propexi[, "AbschGesamt"]) / (1 + max(propexi[, "AbschGesamt"]))))
      npt0 <- (1 + ((max(propexi[, "Parkfitness"]) - propexi[, "Parkfitness"]) / (1 + max(propexi[, "Parkfitness"]))))^k
      NewProb <- 1 / (npt / npt0)
      propwelche[welche %in% indivprop[, "Rect_ID"], "Prop"] <- NewProb
    }

    propwelcheN <- cbind(
      Rect_ID = nGrids1,
      Prop = rep(min(indivprop[, "AbschGesamt"]), nGrids)
    )
    if (trimForce) {
      propexiN <- indivprop[indivprop[, "Rect_ID"] %in% nGrids1, ]
      npt1 <- (1 + ((max(propexiN[, "AbschGesamt"]) - propexiN[, "AbschGesamt"]) /
                      (1 + max(propexiN[, "AbschGesamt"]))))
      npt2 <- (1 + ((max(propexiN[, "Parkfitness"]) - propexiN[, "Parkfitness"]) /
                      (1 + max(propexiN[, "Parkfitness"]))))^k
      NewProb1 <- npt1 / npt2
      propwelcheN[propwelcheN[, "Rect_ID"] %in% indivprop[, "Rect_ID"], "Prop"] <- NewProb1
      if (!all(propwelcheN[, "Rect_ID"] %in% indivprop[, "Rect_ID"] == TRUE)) {
        propwelcheN[!propwelcheN[, "Rect_ID"] %in% indivprop[, "Rect_ID"], "Prop"] <- min(NewProb1)
      }
    }

    propwelcheN <- propwelcheN[!propwelcheN[, "Rect_ID"] %in% welche, , drop = FALSE]
    ## P1 - Deleting Turbines
    prob1 <- propwelche[, "Prop"]
    ## P2 - Adding Turbines
    prob2 <- propwelcheN[, "Prop"]

    if (zviel != 0) {
      if (zviel > 0) {
        if (trimForce) {
          # Delete turbines with Probability
          if (!is.null(seed)) {
            set.seed(as.integer(seed))
          }
          smpra <- base::sample(welche, zviel, replace = FALSE, prob = prob1)
        } else {
          # Delete them randomly
          if (!is.null(seed)) {
            set.seed(as.integer(seed))
          }
          smpra <- base::sample(welche, zviel, replace = FALSE)
        }
        # Delete the 1 entry and make no turbine.
        tmp[smpra] <- 0
        mut1[[i]] <- tmp
      } else {
        if (trimForce) {
          # Add turbines with Probability
          if (!is.null(seed)) {
            set.seed(as.integer(seed))
          }
          smpra <- sample(propwelcheN[, "Rect_ID"], -zviel, replace = FALSE, prob = prob2)
        } else {
          # Add turbines randomly
          if (!is.null(seed)) {
            set.seed(as.integer(seed))
          }
          smpra <- sample(propwelcheN[, "Rect_ID"], -zviel, replace = FALSE)
        }
        # Assign 1 to binary code. So Turbine is created here.
        tmp[smpra] <- 1
        mut1[[i]] <- tmp
      }
    } else {
      mut1[[i]] <- tmp
    }
  }
  mut1 <- do.call("cbind", mut1)
  return(mut1)
}